Severity of Lung Function Impairment Drives Transcriptional Phenotypes of COPD and Relates to Immune and Metabolic Processes

Introduction

Chronic obstructive pulmonary disease (COPD) is a major global health problem with a high illness burden and rising incidence,¹ and is expected to remain a challenge for health-care systems well into the 21st century.² Previously, the disease management of COPD was mainly based on the severity of airflow obstruction, as assessed by forced expiratory volume in 1 second (FEV₁).^3,4 However, research in the field has shown that the one-dimensional analysis of FEV₁ fails to adequately address the multiple dimensions of COPD.⁴ Heterogeneity of COPD is represented in the broad spectrum of respiratory and systemic symptoms, histopathological and radiographic findings, pulmonary physiology and likely different underlying molecular mechanisms (endotypes). There is much recent research interest surrounding the presence of disease phenotypes of COPD that may have different underlying mechanisms and the need for alternative management and treatment approaches.

The use of analytical approaches, such as cluster analysis, has advanced the study of disease phenotypes.⁵ Several possible phenotypes of COPD have been identified in various studies that have used cluster analysis.^5,6 Some of these COPD phenotypes include non-exacerbator, exacerbation prone, emphysematous COPD, non-smoking COPD and the comorbidity or systemic phenotype.^5,7,8 As Barnes et al⁸ noted, “it is likely that molecular phenotyping may lead to the identification of distinct phenotypes linked to specific molecular processes, which can be recognized by specific biomarkers and may be treated more precisely with new or repurposed therapies”.

Transcriptomics can provide useful information relating to disease heterogeneity, mechanisms of pathogenesis, treatment responses and classifications of COPD, along with an opportunity for novel biomarker discovery.^9,10 Over the past few years, a number of studies have examined transcriptomics in COPD cohorts using samples obtained from lung tissue,^11–16 blood,^16–22 and spontaneous²³ and induced sputum.^{16,20–22,24–27} However, these studies have not investigated unbiased clustering of whole transcriptomic profiles, how these clusters relate to clinical and inflammatory phenotypes, and what underlying mechanisms and features are driving these transcriptional clusters. Induced sputum has proven to be very valuable for molecular profiling of COPD, likely due to its less invasive nature compared to obtaining lung tissue.^22,24–27 Gene expression patterns derived from induced sputum have been associated with disease severity,²⁷ COPD inflammatory phenotypes²⁵ and treatment response.²⁰

Studies that utilize cluster analysis combining transcriptomic datasets with COPD-related clinical characteristics, comorbidities and biomarker are important in helping better understand mechanisms underlying the disease as well as strengthening the robustness of any identified COPD phenotypes.²⁸ This study used unsupervised hierarchical clustering of induced sputum gene expression profiles of 72 stable COPD patients from Newcastle area (Australia) to identify distinct and clinically relevant transcriptional COPD phenotypes, and the driving factors behind these cluster phenotypes.

Materials and Methods

Study Design and Population

A cross-sectional analytical study was conducted involving 72 participants with stable COPD and 15 healthy controls. The data for this sub study were obtained from our previously published studies.^29–31 Participants with remaining stored sputum samples available for RNA extraction, where the RNA was of appropriate quality and yield, were included in the present study. All participants provided written informed consent, and ethics approval was obtained from the Human Ethics Research Committees of the Hunter New England Local Health District (06/12/13/3.08, 05/12/07/3.11 and 08/08/20/3.10) and the University of Newcastle (H-2008-0272). This study complies with the ethical principles of the Declaration of Helsinki.

COPD diagnosis was confirmed by incompletely reversible airflow limitation (post-bronchodilator forced expiratory volume in 1 second [FEV1] <80% predicted and FEV1 to forced vital capacity [FVC] ratio of <0.7). Stable COPD was defined as no increase in bronchodilator use, no use of oral corticosteroids or antibiotics, no unscheduled doctor’s visit, or no hospitalization due to COPD in the past 4 weeks. Healthy control participants were recruited by advertisement and were eligible if they had no previous COPD diagnosis and normal lung function (FEV1 >80% predicted, FEV1/FVC >0.7).

Clinical Assessment

Adults (n = 72) with stable physician-diagnosed COPD and healthy controls (n = 15) attended the research centre, data collection included demographic information, lung function, smoking and medical history, medication use, comorbidities (Charlson Comorbidity index (CCI)),³² dyspnea (modified Medical Research Council [mMRC] scale),³³ self-reported prior year exacerbation and hospitalization history, and health-related quality of life (St George’s Respiratory Questionnaire [SGRQ]).³⁴ Airflow limitation was assessed using spirometry (Medgraphics, CPFS/D™ USB Spirometer, BreezeSuite v7.1, MGC Diagnostics, Saint Paul, MN, USA) to measure pre- and post-bronchodilator FEV1, FVC, and FEV1/FVC ratio according to the standards of the American Thoracic Society.³⁵ Severity of COPD was graded according to the Global Initiative for Chronic Obstructive Lung Disease (GOLD) classification system of airflow limitation based on post-bronchodilator FEV₁% predicted (ie, GOLD grade 1, FEV₁ ≥80%; GOLD grade 2, 50%≤ FEV₁ <80%; GOLD grade 3, 30%≤ FEV₁ <50%; GOLD grade 4, FEV₁ <30%).¹ Participants were evaluated using GOLD quadrants according to the refined GOLD “ABCD” assessment tool (using mMRC for symptom assessment).¹ A 6-minute walk test was performed³⁶ and the BODE index (Body mass index [BMI], airflow Obstruction, Dyspnoea and Exercise capacity) calculated.³⁷ Peripheral venous blood was collected and serum high-sensitivity C-reactive protein (hs-CRP) was measured using enzyme-linked immunosorbent assay.

Sputum Induction and Analysis

Sputum collection and analysis were performed as previously described^38–40 and detailed methods are provided in the Methods Section of Supplementary Information.

Sputum Transcriptomics

Transcriptomic data were generated as previously described.⁴¹ Briefly, sputum RNA was extracted using the RNeasy Mini Kit (Qiagen, Hilden, Germany) and reverse transcribed into cRNA and biotin-UTP labelled using the Illumina TotalPrep RNA Amplification Kit (Ambion, Thermo Fisher Scientific, Scoresby, VIC, Australia) and hybridized to the Illumina Sentrix HumanRef-8 Version 2 Expression BeadChips (Illumina, San Diego, CA, USA). Samples and gene profiling results were included in the analysis if the sample was of suitable purity (OD 260/280 1.7–2.1) and was successfully amplified (sufficient cRNA generated) and hybridized (95th percentile of fluorescence score >500) and the data passed quality controls in GeneSpring GX (correlation coefficients and principle component analysis plots; Agilent Technologies, Santa Clara, Calif).

Statistical Analysis

Clinical and cell count data were analysed using Stata 15 (StataCorp, College Station, TX). Results are reported as mean ± standard deviation (SD) for normally distributed data and as median (quartile 1, quartile 3) for nonparametric data. Student’s t-test was used for 2 group comparisons of normally distributed data, and Wilcoxon rank sum test was used for non-parametric data. Comparisons between multiple groups were assessed using Kruskal Wallis for non-parametric data. Comparison of categorical data was done using Fisher’s exact test. P-values <0.05 were considered significant.

Sputum transcriptomic data analysis was performed as previously described,⁴¹ and is described in more detail in the Methods Section of Supplementary Information. The data reported in this paper is deposited in the National Center for Biotechnology Information’s Gene Expression Omnibus (GSE212331). To investigate the natural groupings created by the similarities and differences of the sputum from participants, gene expression profiles were subjected to unsupervised agglomerative hierarchical clustering by using the Euclidean algorithm with Ward’s linkage. Differentially expressed genes that were found to be linked with the only significant factor driving our cluster groupings (ie, GOLD stage) were further validated by investigating their differences in the ECLIPSE data set GSE22148²⁷ which consists of 148 ex-smokers with GOLD stage 2–4 COPD (GOLD stage 2, n = 74 and GOLD stage 3 and 4, n = 74). A brief summary of the ECLIPSE cohort²⁷ is provided in the Methods Section of Supplementary Information.

Results

Clinical Characteristics

Details of the study participants are provided in Table 1. Participants with COPD (n = 72) had a median age of 70 (64, 75.5) years and mean post-bronchodilator predicted FEV₁ of 54 (15.8) %. There were 34 (47.2%) females. Fifty-seven (79.2%) participants were ex-smokers, with median (Q1-Q3) pack-years of 33 (15–63). In terms of GOLD grades, 4 (5.6%) patients were in GOLD I, 39 (54.2%) in GOLD II, 26 (36.1%) in GOLD III, and 3 (4.2%) in GOLD IV. Almost half (45.7%) were “prone to exacerbation”, having had two or more exacerbations in the past 12 months. Most (n = 65, 93%) of the participants were prescribed maintenance inhaled corticosteroids (ICS) or ICS and long-acting β2 agonist (LABA) combination therapy (ICS/LABA) with a median (Q1-Q3) daily dose of 500 (400–500) mg beclomethasone equivalents/day. Healthy controls without COPD (n = 15) were age 25 to 65 and almost half were females (Table 1). Participants with COPD were older, had greater history of smoking, higher body mass index (BMI), more comorbidity burden (higher Charlson Comorbidity Index [CCI] score) and lower lung function. Both percentage and proportion of sputum neutrophils and eosinophils as well as total cell counts were significantly higher in COPD patients compared to healthy controls. On the contrary, percentage sputum macrophage was significantly lower in COPD participants. Details of differences between COPD and healthy controls gene expression profiles are provided in the Results Section Of Supplementary Material (Tables S1 and S2).

Table 1 Demographics and Clinical Characteristics of the Study Population

Unsupervised Hierarchical Clustering of Sputum Gene Expression Profiles in COPD

Unsupervised hierarchical clustering of gene expression profiles showed 2 clear clusters that were quite even with 35 (48.6%) participants in the first cluster and 37 participants in the second cluster (51.4%) (Figure 1). The two clusters further divide into two distinct sub-clusters each. There were 14 (40%) participants in sub-cluster 1.1 and 21 (60%) in sub-cluster 1.2. Sub-clusters 2.1 and 2.2 had 20 (54%) and 17 (46%) participants, respectively.

Figure 1 Phenotypes of COPD identified from unsupervised hierarchical clustering of gene expression profiles (18,511 entities). This approach groups expression profiles that are similar which are joined together to form this tree structure. The root of the tree consists of a single cluster containing all samples, and the leaves correspond to the individual sample profiles, which separate into 2 main clusters and 4 sub-clusters.

Clinical Characteristics of Transcriptional Phenotypes of COPD

The clinical characteristics of the transcriptional COPD phenotypes are described in Table 2. Overall, the proportion of participants in each cluster was comparable and both groups were of similar age and sex. Compared to cluster 1, cluster 2 exhibited worse lung function, more severe disease (higher proportion of GOLD stage 3 and 4), more severe breathlessness (higher mMRC score), increased proportion of sputum neutrophils, and decreased proportion of macrophages. There were no relevant differences in health status, pack-years exposure, CCI, anxiety and depression score, ICS use and other measured inflammatory markers between the two clusters. Although not statistically significant, participants in Cluster 2 had slightly higher BODE score, lower BMI, and fewer proportion of people with comorbid cardiovascular diseases compared to those in cluster 1. There were more people with comorbid eye disease in Cluster 2 compared to Cluster 1, even though this did not reach statistical significance.

Table 2 Clinical Characteristics and Sputum Cell Counts of the Transcriptional COPD Phenotypes

Differences in clinical characteristics and sputum cell counts between the sub-clusters are shown in Table S3 in Supplementary Information. The parameters that were significantly different between the four sub-clusters were lung function (FEV₁) and proportion of people with cardiovascular disease (CVD) comorbidity only (Figure 2).

Figure 2 Differences in (A) lung function (forced expiratory volume in 1 second [FEV1]) and (B) proportion of people with comorbid cardiovascular disease between the four Sub-clusters 1.1, 1.2, 2.1 and 2.2. Notes: *Kwallis2 p = 0.002 vs Sub-cluster 1.2, **Fisher's exact test p < 0.013 vs Sub-clusters 1.2, 2.1 and 2.2.

Differential Gene Expression Associated with Key Clinical and Inflammatory Differences Driving Transcriptional Clusters

To further investigate the mechanisms underlying the transcriptional clusters, a supervised analysis was performed to analyze differences in gene expression between GOLD stage groupings, GOLD quadrant groupings, BODE index quartiles, presence of breathlessness, history of exacerbation, eosinophilic airway inflammation, neutrophilic airway inflammation, and the presence of comorbidities and systemic inflammation. Out of these comparisons, it was only GOLD stage that had differentially expressed genes identified based on our strict statistical criteria. Given this, worsening lung function appears to be the main significant driver responsible for the shifts in sputum transcriptional profiles of COPD that lead to our cluster groupings.

Clinical Features, Inflammatory Characteristics, and Differential Gene Expression Between GOLD Stages

We performed comparisons of clinical features and inflammatory characteristics between GOLD stage 3 and 4 against GOLD stage 1 and 2 (Table 3). GOLD stages were grouped into two groups for the comparison as there were small number of participants in GOLD stage 1 (n = 4) and GOLD stage 4 (n = 3). As shown in Table 3, smoking history, mMRC and BODE scores were significantly different clinical parameters between the two GOLD stage groupings. GOLD stage 3 and 4 participants had significantly higher percentage neutrophil counts, but lower absolute number and proportion of sputum macrophages when compared to their GOLD stage 1 and 2 counterparts.

Table 3 Clinical Features and Inflammatory Characteristics of Study Participants Categorised by 2 GOLD Stage Groupings

There were 4507 out of 23,337 entities significantly differentially expressed between the two GOLD stage groups, with 758 entities that were ≥2-fold change (452 genes after removing duplicates and predicted, 189 genes upregulated and 263 genes downregulated, Table S4 in Supplementary information). A heat map cluster showing the patterns of the genes significantly differentially expressed between the two GOLD stage groupings can be seen in Figure 3, with healthy control levels of expression added as a reference. Figure 3 shows the patterns in the gene expression profiles across the 3 groups, with the GOLD stages 3 and 4 being vastly different to the healthy controls and GOLD stage 1 and 2 groupings.

Figure 3 Heat map cluster showing the pattern of genes differentially expressed in the two GOLD stage groupings and healthy controls with red resembling high expression and blue showing lower expression.

To validate our differentially expressed genes, we further investigated these changes in the ECLIPSE data set (GSE22148).²⁷ We confirmed that out of the 452 genes with ≥2-fold change in our data set, 281 genes (62.2%) were also found to be significantly differentially expressed (adjusted p < 0.05) between these GOLD stages in the ECLIPSE data set (Table S4 in Supplementary information). These 281 validated genes, related to 40 GO biological processes (Table 4) and 5 interaction networks (Figure 4). Biological processes represented by the differentially expressed genes predominantly related to (i) immune and metabolic processes, (ii) metabolic processes in oxidation/reduction and mitochondrial function, (iii) antigen processing and presentation, regulation of immune responses, (iv) lipoprotein and cholesterol metabolism and (v) regulation of complement activation and innate immune responses.

Table 4 Gene Ontology Categories Significantly Overrepresented in the 281 Overlapping Genes

Figure 4 Five protein–protein interaction networks (A–E) among the 281 overlapping differentially expressed genes between our data and the ECLIPSE data set. Network A contains 10 genes relating to lipoprotein and cholesterol metabolism. Network B containing 16 genes relates to metabolic processes in oxidation/reduction and mitochondrial function. Network C contains 10 genes relating to antigen processing and presentation, regulation of immune responses (including T cells, neutrophils). Network D, with 5 genes, shows a network relating to regulation of complement activation and innate immune responses. Network E is the largest containing 25 genes relating to immune and metabolic responses.

Discussion

This study demonstrates the use of unsupervised hierarchical clustering of induced sputum gene expression profiles to identify 2 main clusters and 4 sub-clusters in the studied population. The 2 distinct transcriptional COPD phenotypes relate to both clinical characteristics (lung function and breathlessness) and the type of airway inflammation present (sputum neutrophils). Cluster 1 was characterized by less severe lung function and breathlessness. Cluster 2, on the other hand, was neutrophil-enriched group with worsened airflow obstruction, more severe disease and breathlessness, and lower proportion of macrophages. With regard to the sub-clusters, interestingly, sub-cluster 1.1 had significantly higher proportion of people with comorbid cardiovascular disease compared to the remaining 3 sub-clusters. Through supervised analysis we determined that the predominant factor driving gene expression differences in our transcriptional clusters was lung function (GOLD stage). There were 452 genes (adjusted p < 0.05 and ≥2 fold) altered in GOLD stage 3 and 4 versus 1 and 2. Out of the 452 genes, 281 (62%) were also significantly differentially expressed between the two GOLD groupings in the ECLIPSE data set. Of importance, the validated differentially expressed genes connected in 5 networks and related to the following themes (i) lipoprotein and cholesterol mechanism (including APOE, APOC1 and SCD), (ii) metabolic processes in oxidation/reduction and mitochondrial function (NDUFs), (iii) antigen processing and presentation (including HLA complexes (MHC class II)), (iv) regulation of complement activation and innate immune responses (C1Qs) and (v) immune and metabolic processes (including PSMA5, PSMC2 and PSMDB6).

The past decade in the COPD arena has seen an enormous interest in personalized, precision medicine. To this end, the use of techniques such as cluster analysis to identify groups of COPD patients with similar clinical or physiological characteristics has been the focus of studies.^6,28 For instance, in a study that utilized hierarchical cluster analysis of clinical, functional and imaging data of stable COPD patients, Burgel et al⁴² identified 3 distinct phenotypes of COPD with varying COPD severity and risk of mortality. Garcia-Aymerich et al⁴³ also identified three clusters of COPD, namely moderate, severe and systemic COPD. Studies that have utilized cluster analysis integrating genetic⁴⁴ or transcriptomic data^15,45–47 in COPD cohorts, like ours, have great benefit in that they provide the opportunity to analytically and jointly assess COPD-related clinical characteristics, comorbidities, and biomarker data to strengthen the robustness of the COPD phenotypes as well as to better understand the underlying biological mechanisms of the condition.²⁸ Our study has linked transcriptomic profiles underlying 2 main COPD clusters, and with thorough investigations we identified the main driving factor to these differences in airway gene expression patterns being the degree of lung function impairment. A better understanding of phenotypes of COPD on a deep cellular and molecular level will lead to the development of more targeted strategies for personalized COPD treatment and management. This study has increased the knowledge base on potential mechanistic themes further underlying lung function impairment that warrant further investigation as targets for potential intervention.

In our study, there were several mechanistic themes uncovered through network analysis associated with potential importance to lung function impairment in COPD. These themes were largely related to immune (ie, antigen presentation, complement and innate immune responses) and metabolic processes (ie, lipoprotein and cholesterol metabolism, oxidation/reduction and mitochondrial function). Of interest, all genes in the networks were downregulated in GOLD stages 3 and 4 compared with GOLD stages 1 and 2. This implies that as COPD worsens, there is a shift in gene expression that is downregulation of some important immune and metabolic functions within the cells of the airway lumen.

Antigen presentation serves to ensure adaptive immune responses are initiated to invading microorganisms.⁴⁸ Dysfunction of this vital immune response may lead to increased susceptibility to viral infection and bacterial colonization, which would lead to increased inflammation. HLA-complexes corresponding to MHC Class II were one of the top downregulated class of genes with worsening of GOLD stage. MHC Class II receptors are involved in the presentation of antigens to CD4(+) T-lymphocytes, and thereby are critical for the initiation of the antigen-specific immune response.⁴⁹ Previously we have reported the association of low-dose azithromycin treatment with the downregulation of several human HLA molecules in the airways of neutrophilic COPD patients, resulting in suppression of inflammatory response most likely by limiting the immune system’s capacity to identify and respond to antigen.²²

Our analysis also shows that genes in the complement pathway including C1Q are highly altered. C1Q is a pattern recognition protein that is involved in the recognition and clearance of apoptotic cells. C1Q also directs macrophage polarization and limits inflammasome activity during the uptake of apoptotic cells.⁵⁰ Defect in these key macrophage functions contribute to dysregulated resolution of the downstream inflammatory responses and impaired tissue repair.⁵¹ Recently, Akata et al⁵¹ have shown reduced gene expression of the complement components C1QA and C1QB in non-polarized macrophage cells isolated from BAL liquids of stable COPD patients. An earlier study by Yuan et al⁵² has shown that cigarette smoke downregulates C1Q produced by antigen-presenting cells isolated from emphysematous humans, and that it relates with the severity of lung obstruction, which is in line with our findings.

Macrophages constitute a heterogeneous cell population, and pro-inflammatory M1 and anti-inflammatory M2 and M2-like cells represent the extremities of a pattern of macrophage polarization.⁵³ Studies in the literature have shown that both of these macrophage functions are altered in COPD.⁵³ In relation to this, our analysis demonstrated that genes previously shown to be involved in macrophage polarization, either in COPD or other inflammatory diseases, were downregulated, along with a lower proportion of sputum macrophages. These genes included peroxisome proliferator activated receptor gamma (PPARG)⁵⁴ and the CD74 molecule. The PPARs are known to be involved in the regulation of inflammation and lipid metabolism, important in asthma and COPD.⁵⁵ A recent transcriptomic study of asthma-COPD overlap (ACO) has reported a cluster of ACO that had features of mitochondria and peroxisome dysfunction.⁵⁶ Certainly, the relationship between mitochondrial function, peroxisomes, oxidative stress and inflammation in COPD warrants further investigation in its role in lung function decline.

Another interesting finding within our 4 identified subclusters was the increased rate of cardiovascular comorbidities in sub-cluster 1.1, with 93% of participants within this cluster having cardiovascular conditions. Garcia-Aymerich et al⁴³ also identified a systemic COPD cluster that exhibited high rates of cardiovascular comorbidities, in keeping with our sub-cluster 1.1. COPD and cardiovascular disease frequently occur together, and the coexistence of these conditions may impact disease risk and prognosis.⁵⁷ On differential gene expression analysis between GOLD stages, we also note that genes related to lipoprotein and cholesterol metabolism are higher in the GOLD stage 1 and 2 group compared with the more severe GOLD stage 3 and 4. One of the differentially expressed genes in this network was apolipoprotein E (APOE). APOE plays a critical role in lipid metabolism and has been linked to cardiovascular diseases.⁵⁸ Of note, in our string-network (Figure 4A), APOE is connected with ABCG1. ABCG1 is a member of the superfamily of ATP-binding cassette (ABC) transporter, which is expressed in a number of tissues including lung and involved in movement of cholesterol and phospholipids.⁵⁹ Results obtained from mouse models have shown that ABCG1 plays a critical role in maintaining cellular lipid homeostasis and controlling pulmonary inflammation. As such, understanding the relationship between cholesterol and inflammation in the lung, and the role that ABC transporters like ABCG1 play in this may lead to the identification of new treatment targets.⁶⁰

Our study does have limitations, including a small sample size which will need to be further explored in a larger COPD cohort. However, our lung function differential gene expression findings were validated in a second cohort of COPD and led to identification of key immune and metabolic themes for future investigation. The smaller sample size limited the analysis for the comparison of sub-clusters may have been underpowered. RNA was extracted from mixed cell populations in sputum with bulk RNA-transcriptomic profiles generated, and thus groupings may reflect cell populations present, and signals from rare cell populations will potentially be missed. The majority of our population was taking inhaled corticosteroids, which may impact cell activity and transcription, and may limit the extension of the findings to patients with COPD that are not taking these medications.

Conclusion

In summary, this study has found that gene expression profiling provides a means to investigate the molecular mechanisms and classifications of COPD phenotypes. We identified 2 main clusters and 4 sub-clusters of COPD that relate to important differences in clinical characteristics, inflammatory cells, gene expression, and metabolic and immune function pathways. We found one cluster that was neutrophil-enriched with more severe disease and breathlessness. In our studied population, we have also found a sub-cluster of COPD that had high rates of cardiovascular comorbidity. The most significant characteristics that was driving the transcriptional phenotypes of our study population was the degree of lung function impairment, classified by GOLD stage. Differences between GOLD stages were largely related to genes involved in lipoprotein and cholesterol mechanism, metabolic processes in oxidation/reduction and mitochondrial function, antigen processing and presentation, regulation of complement activation and innate immune responses and immune and metabolic processes. The clinical relevance of our work lies in the fact that it has the potential to assist in clinical investigations of the underlying aetiology of COPD and in guiding effective treatment strategies given that it incorporates not only clinical, inflammatory cell, and comorbidity data but also molecular aspects of COPD.