Distinct skeletal muscle molecular responses to pulmonary rehabilitation in chronic obstructive pulmonary disease: a cluster analysis.

BACKGROUND: Pulmonary rehabilitation (PR) is a cornerstone in the management of chronic obstructive pulmonary disease (COPD), targeting skeletal muscle to improve functional performance. However, there is substantial inter-individual variability in the effect of PR on functional performance, which cannot be fully accounted for by generic phenotypic factors. We performed an unbiased integrative analysis of the skeletal muscle molecular responses to PR in COPD patients and comprehensively characterized their baseline pulmonary and physical function, body composition, blood profile, comorbidities, and medication use.
METHODS: Musculus vastus lateralis biopsies were obtained from 51 COPD patients (age 64 ± 1 years, sex 73% men, FEV1 , 34 (26-41) %pred.) before and after 4 weeks high-intensity supervised in-patient PR. Muscle molecular markers were grouped by network-constrained clustering, and their relative changes in expression values-assessed by qPCR and western blot-were reduced to process scores by principal component analysis. Patients were subsequently clustered based on these process scores. Pre-PR and post-PR functional performance was assessed by incremental cycle ergometry and 6 min walking test (6MWT).
RESULTS: Eight molecular processes were discerned by network-constrained hierarchical clustering of the skeletal muscle molecular rehabilitation responses. Based on the resulting process scores, four clusters of patients were identified by hierarchical cluster analysis. Two major patient clusters differed in PR-induced autophagy (P < 0.001), myogenesis (P = 0.014), glucocorticoid signalling (P < 0.001), and oxidative metabolism regulation (P < 0.001), with Cluster 1 (C1; n = 29) overall displaying a more pronounced change in marker expression than Cluster 2 (C2; n = 16). General baseline characteristics did not differ between clusters. Following PR, both 6 min walking distance (+26.5 ± 8.3 m, P = 0.003) and peak load on the cycle ergometer test (+9.7 ± 1.9 W, P < 0.001) were improved. However, the functional improvement was more pronounced in C1, as a higher percentage of patients exceeded the minimal clinically important difference in peak workload (61 vs. 21%, P = 0.022) and both peak workload and 6 min walking test (52 vs. 8%, P = 0.008) upon PR.
CONCLUSIONS: We identified patient groups with distinct skeletal muscle molecular responses to rehabilitation, associated with differences in functional improvements upon PR.

Introduction

Chronic obstructive pulmonary disease (COPD) is increasingly recognized as a complex chronic disease with extra‐pulmonary manifestations and comorbidities, including skeletal muscle dysfunction.1 Skeletal muscle dysfunction in COPD arises from structural and metabolic alterations such as a loss of muscle mass, a shift in fibre type distribution, decreased oxidative capacity, and mitochondrial dysfunction.2 Importantly, besides affecting physical functioning and health‐related quality of life,3 skeletal muscle dysfunction is a predictor of mortality independent of lung function.4, 5

In clinically stable COPD patients, exercise training is currently the most potent intervention to improve skeletal muscle function. As such, exercise training forms a cornerstone in pulmonary rehabilitation (PR) in patients with COPD. However, despite the overall positive effects of exercise training on muscle mass and function at the group level, there is substantial inter‐individual variability in the effect of PR on skeletal muscle function.2, 6, 7, 8, 9, 10 The sources of this variability currently remain unexplained but may be better understood when we improve our comprehension of the underlying skeletal muscle molecular response pattern. In the older population, which is prone to the loss of muscle mass and strength (i.e. sarcopenia), the exercise response can be blunted due to impaired skeletal muscle regeneration.11 This potentially results from an impaired induction of skeletal muscle remodelling‐related processes, as reflected by a blunted induction of protein degradation,12 anabolic resistance,13 and a delayed or blunted satellite cell response.11 Similarly, the exercise response in COPD patients may be limited by impaired satellite cell mediated skeletal muscle regeneration.14 Furthermore, we recently showed an elevated basal activity of skeletal muscle remodelling‐related processes in COPD patients.15 This may hinder the normal induction of these processes during skeletal muscle recovery after disuse or upon exercise training.12

Several studies assessed molecular markers of remodelling‐related processes upon exercise training in COPD patients,10, 16, 17, 18 but none report a clear differential rehabilitation response between predefined phenotypic categories, for example, based on body composition. This may be partially due to the limited sample size in most of these studies or the individual assessment of a limited selection of markers, which often display a large variation. Nevertheless, age, body composition, and comorbidities, such as metabolic and heart diseases, have been weakly associated with the functional response to rehabilitation.16, 19, 20 Other studies, however, report inconclusive or even contradicting results,6, 7, 8, 9, 17, 21, 22 strongly indicating that these phenotypic factors do not fully account for the variability in rehabilitation responses.

In the current study, we take a clustering approach to—for the first time—unbiasedly address inter‐individual heterogeneity in the skeletal muscle molecular response to PR of patients with advanced COPD.

Materials and methods

Study design and participants

We selected COPD patients from a prospective observational study conducted in the in‐patient PR unit at the University Clinic Golnik. The study design has previously been published23 (ClinicalTrials.gov identifier: NCT02550808). The study was approved by the Slovenian National Medical Ethics Committee (Ljubljana, Slovenia) and was carried out in accordance with the ethical standards laid down in the 1964 Declaration of Helsinki and its later amendments, and all participants provided written informed consent. Patient selection for the current analysis was based on completion of the 4 week rehabilitation programme and availability of both a pre‐rehabilitation and post‐rehabilitation muscle biopsy specimens, yielding 51 patients.

Pulmonary and physical function

Spirometry was used to obtain FEV1, FVC, and their ratio (FEV1/FVC).24 Physical function was assessed with the 6 min walking test (6MWT)25 and by measurement of peak load (W) and peak oxygen consumption (VO2, mL/kg/min) by an incremental load cycling cardiopulmonary exercise test.26 In addition, the improvements in physical function were categorized by the minimal clinically important difference (MCID; 6MWT: 25 m,27 peak load: 10 W28).

Body composition

Whole body dual‐energy X‐ray absorptiometry was used to assess total and appendicular (i.e. arms and legs) fat mass and fat‐free mass at baseline, as previously described.23 Fat mass index (FMI), fat‐free mass index (FFMI), and appendicular skeletal muscle mass index were calculated as the ratio of mass to height2 (kg/m2). Presence of muscle mass depletion was defined by the cut‐offs for appendicular skeletal muscle mass index (<7.23 kg/m2 for men; <5.76 kg/m2 for women).29 Bioelectrical impedance (BIA; QuadScan 4000; Bodystat, Douglas, United Kingdom) was performed pre‐PR and post‐PR as previously described23 and was used to assess the change in FMI and FFMI upon PR.

Skeletal muscle molecular analyses

Pre‐PR and post‐PR biopsies were obtained from the vastus lateralis muscle of the dominant leg by needle biopsy, one day (≥20 h) after the last exercise test. Muscle tissue was snap‐frozen in liquid nitrogen and stored at −80°C until processed for molecular analyses; mRNA levels and protein abundance of mediators and regulators of muscle mass and metabolism were determined by RT–qPCR and western blotting, respectively. Details are provided in Supporting Information.

Data analysis

Clustering

Details on the clustering analyses are provided in Supporting Information. Briefly, molecular rehabilitation responses, calculated as the percentage change from pre‐rehabilitation values, were grouped by network‐constrained hierarchical clustering of the markers. The resulting processes were named based on enrichment of gene ontology (GO) biological processes. Data dimensionality was reduced by computing a single score for each process using principle component analysis. Subsequently, patients were hierarchically clustered based on their individual process scores. Analyses were performed in Python using the scikit‐learn library, and results were visualized using Graph‐tool.

Statistics

After checking for normality, rehabilitation responses were tested by paired‐sample t‐test. Between‐group differences were tested by independent sample t‐test (two groups), by one‐way analysis of variance with Bonferroni post‐hoc comparisons (>2 groups), or by their non‐parametric or categorical equivalents. Correlations between molecular processes were tested using Pearson's correlation coefficient, and between‐group differences in correlations were tested by a Fisher's r‐to‐z transformation. Analyses were performed using spss Statistics (version 22.0, IBM Corp., Armonk, NY, USA). A P value <0.05 was considered statistically significant.

Results

General characteristics

Fifty‐one patients were included in the analyses, with subject characteristics displayed in Table 1. The patient group consisted mostly out of men, current smokers, and 76% of patients completed the PR without an exacerbation. Based on the Global Initiative for Chronic Obstructive Lung Disease classification,30 most patients had severe (stage III; 53%) to very severe COPD (stage IV; 33%), and 60% of patients had muscle mass depletion.

Table 1

Baseline subject characteristics

	Whole group	Cluster 1	Cluster 2	P valuea
	n = 51	n = 29	n = 16
Demographics
Age, years	64 ± 1	64 ± 2	65 ± 2	ns
Sex, m/f (% male)	37/14 (73)	21/8 (72)	11/5 (69)	ns
Smoking status				<0.1
Current, n (%)	44 (88)	24 (83)	16 (100)
Former, n (%)	6 (12)	5 (17)	0 (0)
Completion of PR				ns
Without exacerbation, n (%)	39 (76)	23 (79)	11 (69)
With exacerbation, n (%)	12 (24)	6 (21)	5 (31)
Pulmonary function
FEV1, % predictedb	34 [26–41]	34 [26–41]	38 [29–44]	ns
FVC, % predicted	79 ± 3	80 ± 3	77 ± 5	ns
FEV1/FVC, %b	33 [26–42]	31 [26–40]	33 [28–46]	ns
Physical function
6MWT, m	344.68 ± 15.44	352.59 ± 22.11	348.00 ± 20.94	ns
Peak load, W	62.96 ± 3.30	62.81 ± 4.81	65.20 ± 5.51	ns
Peak VO2, mL/kg/min	11.75 ± 0.58	12.31 ± 0.85	11.35 ± 0.77	ns
Body composition
BMI, kg/m2	24.8 ± 0.6	24.2 ± 0.7	25.5 ± 1.3	ns
FMI, kg/m2	7.8 ± 0.3	7.4 ± 0.4	8.3 ± 0.7	ns
FFMI, kg/m2	16.9 ± 0.4	16.7 ± 0.5	16.9 ± 0.9	ns
ASMI, kg/m2	6.6 ± 0.2	6.5 ± 0.2	6.9 ± 0.4	ns
Muscle mass depletion, n (%)	30 (60)	20 (69)	7 (47)	ns
Blood profile
Triglycerides, mmol/L	1.28 ± 0.08	1.16 ± 0.10	1.46 ± 0.14	<0.1
HDL cholesterol, mmol/L	1.77 ± 0.09	1.80 ± 0.10	1.73 ± 0.16	ns
LDL cholesterol, mmol/L	3.01 ± 0.15	2.77 ± 0.19	3.25 ± 0.26	ns
LDL/HDL, ratio	1.94 ± 0.14	1.65 ± 0.14	2.21 ± 0.27	<0.1
HOMA‐IRb	3.40 [1.80–6.70]	3.70 [1.95–5.75]	2.90 [1.85–8.00]	ns
Comorbidities
Long‐term oxygen therapy, n (%)	20 (39)	11 (38)	5 (31)	ns
Heart failure, n (%)	8 (16)	5 (17)	2 (13)	ns
Ischemic heart disease, n (%)	4 (8)	3 (10)	0 (0)	ns
Atrial fibrillation, n (%)	1 (2)	1 (3)	0 (0)	ns
Arterial hypertension, n (%)	21 (41)	10 (34)	8 (50)	ns
Cholesterolaemia, n (%)	14 (27)	6 (21)	5 (31)	ns
Type 2 diabetes mellitus, n (%)	7 (14)	5 (17)	1 (6)	ns

6MWT, 6 min walk test; ASMI, appendicular skeletal muscle mass index; BMI, body mass index; FEV1, forced expiratory volume in 1 s; FFMI, fat‐free mass index; FMI, fat mass index; FVC, forced vital capacity; HDL, high‐density lipoprotein; HOMA‐IR, homeostatic model of insulin resistance; IQR, interquartile range; LDL, low‐density lipoprotein; PR, pulmonary rehabilitation. Due to the insufficient sample size, individual values for Cluster 3 (n = 4) and Cluster 4 (n = 2) are not depicted in this table. Data expressed as mean ± SEM unless indicated otherwise.

Significance of between group (Cluster 1 vs. Cluster 2) comparisons assessed by independent sample t‐test, Mann–Whitney U‐test, or χ2 test.

Median [IQR].

Functional and skeletal muscle molecular rehabilitation responses in chronic obstructive pulmonary disease patients

Following PR, both 6 min walking distance (+26.5 ± 8.3 m, P < 0.01) and peak load on the cycle ergometer test (+9.7 ± 1.9 W, P < 0.001) were improved. Moreover, 25 patients (51%) exceeded the MCID on the 6MWT, and 17 patients (40%) exceeded the MCID for peak load upon PR.

Using data from all 51 patients, eight molecular processes (i.e. P1–P8) were discerned by network‐constrained hierarchical clustering of the skeletal muscle molecular rehabilitation responses (Figure 1A, Supporting Information, ) and were named based on gene ontology term overrepresentation as indicated in Supporting Information, . In general, P1 (autophagy) related mRNA expression was decreased, whereas autophagy‐related protein phosphorylation and expression were increased after PR (Figure 1B). Furthermore, P2 (AKT/mTOR signalling) related mRNA expression, and protein phosphorylation and expression were predominantly increased upon rehabilitation (Figure 1B). Similarly, P3 (myogenesis) related mRNA and protein expression and P4 [oxidative phosphorylation capacity (OXCAP)] and P5 (mitophagy, FUNDC1/DNM1L) related protein expression were increased upon rehabilitation, whereas P6 (mitophagy, PRKN) related protein expression was decreased (Figure 1B). Furthermore, P7 (glucocorticoid signalling) related inhibitory phosphorylation seemed to increase upon rehabilitation, and in line with this, glucocorticoid signalling‐related downstream mRNA expression was decreased (Figure 1B). Despite the significant decrease in mRNA expression of transcriptional regulators of mitochondrial biogenesis (TFAM and PPRC1) and the increase in a regulator of mitochondrial quality (PINK1), P8 (oxidative metabolism regulation) related markers displayed no consistent overall pattern of change upon rehabilitation (Figure 1B).

Figure 1

Molecular network and network‐constrained hierarchical clustering of molecular rehabilitation responses. (A) Network‐constrained clusters revealing eight distinct processes (P), as indicated with different colours. (B) Molecular rehabilitation responses. The literature‐based molecular network is indicated as lines between molecular markers. Circles represent mRNA markers; pentagons represent protein markers. Numbers correspond to individual markers as depicted in Supporting Information, .

Identification of patient clusters with differential skeletal muscle molecular rehabilitation responses

Both the functional and skeletal muscle molecular rehabilitation responses displayed a substantial inter‐individual variation. To gain insight in the determinants of this variation, an unbiased approach was used to classify patients based on their molecular rehabilitation responses. To this end, one process score per identified molecular process was computed by principle component analysis, with component loadings indicated in Supporting Information, . Four clusters (C1–4) of patients were identified by hierarchical cluster analysis based on process scores, using a cut‐off based on the local optimum in silhouette score (Figure 2A and 2B). C3 and C4 had an insufficient sample size (n = 4 and n = 2, respectively) for the statistical detection of relevant group differences and were therefore omitted in further analyses.

Figure 2

Hierarchical clustering of patients. (A) Silhouette scores for n clusters (mean ± SEM). Clustering based on raw molecular rehabilitation responses (marker values; i.e. no data reduction) is indicated in grey; clustering based on eight process scores is indicated in black. (B) Patients individual silhouette coefficient values per cluster. (c) Dendrogram and clustered heatmap of individual marker rehabilitation responses (Z‐scores) and process scores.

C1 (n = 29) and C2 (n = 16) differed in rehabilitation‐induced modulation of autophagy, myogenesis, glucocorticoid signalling, and oxidative metabolism regulation and tended to differ by regulation of OXCAP (Figure 2C, Supporting Information, ). Specifically, the difference in rehabilitation‐induced change in autophagy markers between C1 and C2 was mainly reflected by a differential change in SQSTM1, MAP 1LC3B, BECN1, BNIP3, ULK1, GABARAPL1, and OPTN mRNA expression (Supporting Information, ). Although both clusters displayed an overall increase in autophagy‐related protein expression, a differential change in MAP 1LC3BI and a tendency towards a differential change in MAP 1LC3BII/I were observed (Supporting Information, ). Conversely, the differential regulation of myogenesis was only reflected by a significant difference in the rehabilitation‐induced change in MYOD1 and CDH15 mRNA expression (Supporting Information, ). In contrast to myogenesis, the differential regulation of glucocorticoid signalling consisted of a differential rehabilitation‐induced change in all related markers, although FOXO3 mRNA did not reach significance (Supporting Information, ). The difference in the rehabilitation‐induced modulation of oxidative metabolism regulation was reflected by a differential change in PINK1, OPA1, FIS1, COX4I1, MYH2, MFN1, TFAM, PPARGC1A, PPRC1, PPARGC1B, PRKN, MYH1, and SOD2 mRNA expression (Supporting Information, ). Interestingly, although the tendency towards a differential regulation in OXCAP was not apparent at the marker level, a significant rehabilitation‐induced increase in NDUFB8, SDHB, MT‐CO1, ATP5F1A, and TOMM20 protein expression was observed for C1, while these markers remained unaltered in C2 (Supporting Information, ).

Association between muscle molecular rehabilitation responses

To assess whether a differential response at the process and marker level resulted in a differential coordination between processes, the correlations between scores of the eight identified processes were computed and compared between the two clusters (Figure 3).

Figure 3

Associations between molecular rehabilitation responses. Each square displays a scatter plot and regression line of the processes indicated on the x and y axis, for both Cluster 1 and Cluster 2. Differences between correlations were tested by a Fisher's r‐to‐z transformation and indicated with a red outline (P < 0.05) or orange outline (P < 0.1). Pearson correlation coefficients per cluster are depicted for differential correlations and when P < 0.1. # P < 0.1, * P < 0.05, ** P < 0.01, *** P < 0.001.

Differential correlations for C1 and C2 were observed between autophagy and oxidative metabolism regulation, between AKT/mTOR signalling and OXCAP, and between mitophagy (FUNDC1/DNM1L) and mitophagy (PRKN) (all P < 0.05) (Figure 3). In addition, AKT/mTOR signalling tended to be differentially correlated with myogenesis (P = 0.06), mitophagy (PRKN) (P = 0.07), and glucocorticoid signalling (P = 0.05) (Figure 3). Moreover, oxidative metabolism regulation tended to be differentially correlated with OXCAP (P = 0.07) and glucocorticoid signalling (P = 0.09) (Figure 3).

Cluster characteristics and functional rehabilitation responses

C1 and C2 did not differ significantly in demographics, pulmonary function, or body composition (Table 1) nor comorbidities or use of medication (Table 1, Supporting Information, ). Likewise, C1 and C2 did not differ in baseline values of the functional parameters (Table 1) and average rehabilitation‐induced change in FFMI or FMI assessed by BIA (Supporting Information, ). Furthermore, these clusters did not differ significantly in the average rehabilitation‐induced change in 6MWT or peak load (Figure 4A and 4B). However, a pronounced and significant increase in 6MWT was observed in C1, whereas the change in 6MWT did not reach significance in C2 (Figure 4A). Similarly, a pronounced and significant rehabilitation‐induced increase in peak load was observed in C1 but not in C2 (Figure 4B).

Figure 4

Functional rehabilitation responses per cluster. Individual rehabilitation‐induced changes. (A) Distance (metres) walked in 6 min walk test (6MWT), n = 28/15. (B) Peak load (W) on a cycle ergometer test, n = 23/14. Data expressed as mean ± SD. (C) Percentage of patients with a change in 6MWT, peak load, or both, exceeding the minimal clinically important difference (MCID; 6MWT: 25 m, peak load: 10 W). * P < 0.05, ** P < 0.01, *** P < 0.001, indicating significance of within‐group rehabilitation responses or significance of differences between indicated groups.

Categorization of rehabilitation responses based on the MCID revealed that the distribution of change in peak load differed between C1 and C2, with a larger amount of patients exceeding the MCID in C1 (61 vs. 21%, P < 0.05) (Figure 4C). In accordance, C1 contained more patients that exceeded the MCID for both the 6MWT and peak load than C2 (52 vs. 8%, P < 0.01) (Figure 4C).

Discussion

To our knowledge, this study is the first to perform unbiased clustering of COPD patients based on their skeletal muscle molecular response to PR. Using this approach, we show that clustering based on the molecular response to PR leads to the identification of patient groups who differ by their functional response to PR.

Variability in the response to pulmonary rehabilitation

Study participants followed a supervised, high‐intensity in‐patient PR programme, which targeted skeletal muscle using a comprehensive exercise training programme.31 In line with previous studies,32 PR improved physical performance assessed by 6MWT and peak workload on a cycle ergometer test, but the effect was highly variable. We measured a panel of molecular markers that are important for exercise‐induced skeletal muscle mass and metabolic plasticity33, 34, 35, 36 and, similarly, observed high variability in the rehabilitation‐induced changes in skeletal muscle molecular marker expression levels. Therefore, molecular markers were clustered into processes and assessed within the context of co‐clustered markers. Subsequently, an integrative clustering approach was used for unbiased identification of patient groups with a differential skeletal muscle molecular response to PR.

Patient clusters with differential skeletal muscle molecular rehabilitation responses

We identified two major patient clusters, which differed in both the direction and the magnitude of the molecular response to PR. Overall, C1 displayed the most pronounced changes in markers of autophagy, myogenesis, glucocorticoid signalling, oxidative metabolism regulation, and OXCAP in response to PR.

Specifically, C1 displayed a change in autophagy‐related marker expression reflective of a decrease in autophagy and a corresponding decrease in glucocorticoid signalling‐related mRNA expression upon PR, as is in line with the PR‐response in previous studies.16 Furthermore, corresponding to previous observations,10, 16, 18 AKT/mTOR signalling was induced upon PR, resulting in an increased capacity and stimulatory signalling of protein synthesis. In keeping with an induction of anabolic signalling upon PR in C1, proliferation marker expression was mostly increased. However, in contrast to previous studies,10, 16 mRNA expression of the myogenic differentiation markers MYOD1 and CDH15 were reduced, while mRNA expression of the fusion marker MYMK was induced upon PR in C1. Furthermore, a decrease in expression of oxidative metabolism regulation markers was observed in C1. Nevertheless, this does not negate the previously reported importance of oxidative metabolism regulation in the coordination of exercise training‐induced metabolic alterations.37 Indeed, as anticipated upon exercise training,38 OXCAP protein expression was increased upon PR in C1. Moreover, the PR‐induced modulation of mitophagy‐related markers in C1 is in agreement with an important role for mitochondrial dynamics and mitophagy in the regulation of mitochondrial quality upon exercise training, as reviewed by Drake et al.37 Together, these data indicate a relative decrease in catabolic signalling, a relative increase in anabolic signalling, and an increase in mediators of oxidative metabolism upon PR in C1, which strikingly reflects a normal response to exercise training.39 Importantly, the ‘normal’ molecular response observed in C1 was accompanied by a strong improvement in physical functioning in this cluster, as determined by an increase in peak load on the cycle ergometer test and an improvement in walking distance in the 6MWT.

In C2, both the functional and molecular responses to PR were less pronounced. This cluster significantly differed from C1 in the rehabilitation‐induced regulation of autophagy, myogenesis, glucocorticoid signalling, and oxidative metabolism regulation. Specifically, expression of autophagy‐related mRNA and glucocorticoid signalling markers remained unaltered, or even increased upon PR in C2, which is in line with previous observations,10, 18 and may be due to disease‐related factors such as hypoxia or systemic inflammation.40, 41 However, an early transient induction of catabolic processes is also evident upon exercise and may play an important role in muscle regeneration.42 Furthermore, the differential regulation of myogenesis is reflected by a differential change in the expression of specific markers. Interestingly, the reduction in the myogenesis markers MYOD1 and CDH15, and induction in MYMK mRNA expression upon PR in C1, was previously shown to be reflective of late‐stage muscle regeneration.43, 44 Conversely, the apparent increase in MYOD1 and CDH15 mRNA expression, together with only a tendency towards an increase in MYMK upon PR in C2, is reflective of an earlier phase of muscle regeneration.43, 44 Despite differences in the regulation of autophagy, myogenesis, and glucocorticoid signalling between C1 and C2, these processes correlate similarly in both clusters, which shows that there is no cluster‐specific dissociation between these processes. This may indicate that, with respect to muscle mass plasticity, both clusters display a normal physiological response to PR but may reflect distinct phases of PR‐induced muscle remodelling.

Interestingly, we found a tendency towards a differential association between oxidative metabolism regulation and OXCAP in Clusters 1 and 2 and a differential correlation between protein turnover and mitochondrial turnover regulatory processes. This supports previous studies showing a correlation between tissue remodelling and bioenergetics pathways that were altered in COPD36, 43, 45 and adds that this alteration may be present only in a subgroup of patients. Indeed, as expected upon exercise training,39 C1 displayed an increase in OXCAP markers upon PR, reflecting a restoration of the oxidative metabolic machinery. This is corroborated by the decreased expression of oxidative metabolism regulation markers, as a normalization of metabolism‐related mRNA has previously been shown in late‐stage muscle regeneration.43, 46 Conversely, in C2, oxidative metabolism regulation and OXCAP were statistically unaltered, although oxidative metabolism regulation seemed to increase upon PR. Furthermore, C2 displayed a strong induction in PINK1 mRNA expression, which plays an important role in the regulation of mitochondrial quality.47 Together, this suggests that the oxidative metabolic machinery is incompletely restored in C2, which may indicate that this process is impaired—as previously suggested in COPD48—or yet ongoing.

Potential drivers of cluster differences

Compellingly, and in line with previous studies,6, 7, 9, 22, 49 baseline variables did not seem predictive of the muscle molecular or functional response, as patient clusters did not differ in baseline values of functional parameters, demographics, pulmonary function, body composition, presence of comorbidities, or use of medication. Furthermore, both the current physiological and molecular data do not provide clear indications of a ‘defective’ response in either cluster. Rather, we propose that molecular differences between clusters are reflective of distinct phases of muscle remodelling, which may be due to differences in the speed by which the remodelling process occurs and/or the timing of initiation of muscle remodelling‐related processes upon PR. Interestingly, a differential response to PR may reflect the underlying cause of exercise limitation,50 including reduced muscle capillarization and neuromuscular function as reported upon ageing,11, 12 which delays or slows PR‐induced muscle mass and metabolic remodelling processes. This may be influenced by lifestyle factors, which is corroborated by the tendency towards a difference in smoking status and blood lipid profile between the clusters. Alternatively, the molecular changes reflecting an early stage of muscle remodelling may originate from an impaired progression through regulating processes, preventing their late‐stage normalization, as may be the case for myogenesis.51 However, the drivers of differential clustering remain elusive and will be an important target for future investigations.

Strengths and limitations

We acknowledge that future longitudinal studies are needed to discern a potential defect in muscle mass and metabolic plasticity from a differential progression of these processes. Nevertheless, the use of a relatively short intervention period was likely critical in the identification of clusters with molecular differences reflective of distinct phases of the muscle remodelling process.

As we aimed to address the variability in the skeletal muscle rehabilitation response between COPD patients, no healthy control group was included. Furthermore, although no PR‐induced change in FFMI was detected by BIA, we did not measure changes in muscle mass using a sensitive measure such as CT or MRI. Therefore, we cannot confirm our assumption that C1 displays a more ‘healthy’ response to PR than C2. However, it is unlikely that the delivered training load differed between clusters, as the training programme was conducted in a controlled, supervised in‐patient setting. Moreover, the training modality was based on baseline characteristics, which did not differ between clusters.

The current study population was characterized by severe COPD and a high prevalence of muscle mass depletion. Therefore, distribution of patients among clusters cannot be directly extrapolated to the general COPD population. In line, it should be noted that aside from the two extensively discussed major clusters, two other patient clusters have been identified. We chose to omit these two very small clusters in further analysis, rather than recursively merging clusters based on their hierarchical structure, because they showed a distinct molecular response to PR. However, these clusters may thus be clinically meaningful and call for future investigation.

To support identification of robust, clinically relevant subgroups, patients were clustered based on process scores. The physiological relevance of the identified processes is substantiated by the separation of two mitophagy processes that have been reported in literature52 and was achieved by inclusion of prior knowledge of the molecular network into the marker clustering.53 To this end, we relied on well‐established knowledge of their interactions and previous unbiased network analyses of ‘omics’ data.36, 45 In contrast to these studies, we have performed a targeted analysis of a selection of molecular markers, which can be confined in future studies as several markers show a low contribution to the process scores. Nevertheless, component loadings did not imply single biomarkers for the identified processes. Moreover, assessment of single marker expression proved important for the interpretation of observed differences in processes, certainly in the context of temporal responses, and therefore remains necessary for obtaining insight in the mechanisms underlying muscle dysfunction in COPD.

Conclusions

Based on the cluster‐based analysis, no specific cluster with ‘non‐responders’ was identified, reinforcing that exercise training should form a cornerstone of PR to improve skeletal muscle function in advanced COPD. However, this study identifies two major clusters of COPD patients with a differential muscle molecular rehabilitation response, which may be reflective of distinct phases of muscle remodelling, and corresponds to a differential gain in physical functioning.

Identification of the drivers of these differential temporal responses will provide an opportunity for personalized interventions optimizing the exercise training response. We propose that, next to detailed analyses of baseline lifestyle factors such as smoking, diet and physical activity, and lung function impairments and blood gasses, future studies should investigate muscle intrinsic differences that determine the kinetics of the muscle molecular response to PR. This could open up new avenues for pharmacological and nutritional interventions to optimize the skeletal muscle functional improvements upon PR.

Conflict of interest

The authors declare that they have no conflict of interest.

Table S1. Compositions of network‐constrained clusters of molecular rehabilitation responses.

Table S2. Component loadings of markers on each process

Table S3. Rehabilitation responses per cluster, process scores

Table S4. Rehabilitation responses per cluster, molecular markers

Table S5. Medication use

Click here for additional data file.