Different Biological Pathways Are Up-regulated in the Elderly With Asthma: Sputum Transcriptomic Analysis.

BACKGROUND: Elderly asthma (EA) is increasing, but the pathogenesis is unclear. This study aimed to identify EA-related biological pathways by analyzing genome-wide gene expression profiles in sputum cells.
METHODS: A total of 3,156 gene probes with significantly differential expressions between EA and healthy elderly controls were used for a hierarchical clustering of genes to identify gene clusters. Gene set enrichment analysis provided biological information, with replication from Gene Expression Omnibus expression profiles.
RESULTS: Fifty-five EA patients and 10 elderly control subjects were enrolled. Two distinct gene clusters were found. Cluster 1 (n = 35) showed a lower eosinophil proportion in sputum and less severe airway obstruction compared to cluster 2 (n = 20). The replication data set also identified 2 gene clusters (clusters 1' and 2'). Among 5 gene sets significantly enriched in cluster 1 and 3 gene sets significantly enriched in cluster 2, we confirmed that 2 were significantly enriched in the replication data set (OXIDATIVE_PHOSPHORYLATION gene set in cluster 1 and EPITHELIAL MESENCHYMAL TRANSITION gene set in cluster 2').
CONCLUSIONS: The findings of 2 distinct gene clusters in EA and different biological pathways in each gene cluster suggest 2 different pathogenesis mechanisms underlying EA.

INTRODUCTION

Global society is aging.1 Accordingly, elderly asthma (EA) involving asthma in people aged 65 years and over is expected to increase. The estimated prevalence of EA in developed countries is 6%-10%.2 Previous reports suggest that EA is phenotypically different from non-elderly asthma (NEA).345 However, the pathogenesis underlying EA has not been clearly elucidated.

Induced sputum is a reliable non-invasive assessment of bronchial inflammation in asthma.6 The evaluation of gene expression profiles of sputum cells has been successfully applied to understand asthma pathogenesis.78 The present study searched biological pathways related to EA using genome-wide gene expression profiles of sputum cells from EA. We selected genes showing significantly differential expressions between EA and healthy elderly controls and identified 2 distinct clusters by hierarchical clustering of these genes. We performed gene set enrichment analysis (GSEA) to gain further biological insights of each cluster. Finally, we confirmed that our findings were replicated in a dependent gene expression profile of sputum cells obtained from Gene Expression Omnibus (GEO), a publicly available database of gene expression profiles. To the best of our knowledge, this is the first study to identify peculiar biological pathways contributing to the EA pathogenesis based on sputum gene expression profiles.

MATERIALS AND METHODS

Discovery data set

Participants aged 65 years or older were recruited from the Seoul National University Hospital (Seoul, Korea). Asthma was diagnosed according to the Global Initiative for Asthma guideline on the basis of current (past 12 months) episodic respiratory symptoms and demonstrated evidence of airway hyperresponsiveness to methacholine or positive bronchodilator (BD) response.9 Healthy elderly controls without asthma (n = 10) were recruited by advertisement. None of them had respiratory symptoms, sputum eosinophilia, abnormality on chest radiography and obstructive pulmonary function test results. Three of them were smokers (1 current smoker and 2 ex-smokers). Exclusion criteria included recent (past month) respiratory tract infection, change in maintenance therapy and asthma exacerbation. However, to reflect real situations, a smoking status was not considered an exclusion criteria. Sputum induction and processing were performed as previously described.10

Gene expression arrays

RNA was extracted from induced sputum samples by using the RNeasy Mini Kit (Qiagen, Hilden, Germany). Gene expression levels were measured using the GeneChip Human Gene 2.0 ST (Affymetrix, Santa Clara, CA, USA). We removed probes with bad chromosome annotation and probes in X or Y chromosome. We then did variance stabilizing transformation and quantile normalization, respectively, to reduce the effects of technical noises and to make the distribution of expression level for each array closer to a normal distribution.

Statistical analysis

A total of 3,156 gene probes showing significantly different expressions between EA and healthy elderly controls (P < 0.05) were used to further analysis. To search for meaningful information patterns and dependencies in gene expression data, we performed hierarchical clustering using the pvclust package in R version 3.4.3 (www.r-project.org; R Foundation for Statistical Computing, Vienna, Austria). This package provides an approximately unbiased P value generated by multi-scale bootstrap resampling. The P value indicates how strong the cluster is supported by the data.11 An approximately unbiased P value greater than 95% was used to define a cluster. We next performed GSEA using the GSEA software (version 3.0) provided by the Broad Institute (Boston, MA, USA).12 We used the hallmark gene sets (H collection) from the Molecular Signatures Database (MSigDB, version 6.0) and defined a significantly enriched gene set when a false discovery rate threshold was less than 0.05.

Ethics approval

This study was approved by the Seoul National University Hospital Review Board (1608-101-786), and informed consent was obtained from all study participants.

Replication data set

A dependent gene expression profile of sputum cells (GSE41863) obtained from GEO (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE41863) was used to replicate our results. To identify markers associated with various asthma subtypes, sputum samples were collected from asthmatics and healthy controls and were subjected to expression profiling using Affymetrix HG-U133Plus2.0. From this profile, we selected 20 subjects aged 65 years or older (15 asthmatics and 5 healthy controls) and identified 3,264 gene probes showing significantly (P < 0.05) different expressions. Hierarchical clustering and GSEA were performed using these gene probes.

RESULTS

Fifty-five patients with EA were enrolled. Based on the differential gene expression patterns of sputum cells, 2 distinct clusters were identified (Fig. 1A and Supplementary Fig. S1). Cluster 1 consisted of 35 patients with EA. Cluster 1 featured a significantly lower proportion of eosinophils in the sputum and less severe airway obstruction as measured by the post-BD ratio of the forced expiratory volume in 1 second and forced vital capacity (FEV1/FVC) compared to cluster 2. Detailed characteristics of the 2 clusters are provided in Table 1.

Fig. 1

Two gene clusters identified in the discovery and replication dataset. (A) Discovery dataset. Three outliers (Pt5, Pt21, and Pt28) were excluded from analysis. (B) Replication dataset.

Pt, Patient.

Table 1

Baseline characteristics (discovery dataset)

Characteristics		Cluster 1 (n = 35)	Cluster 2 (n = 20)	P value
Male, No. (%)		16 (45.7)	9 (45)	0.959
Age (yr)		72.6 (5.4)	75.1 (5.4)	0.103
BMI (kg/m2)		24.9 (3.7)	26.4 (4.1)	0.205
Smoking status				0.841
	Current smoker, No. (%)	4 (11.4)	2 (10)
	Ex-smoker, No. (%)	12 (34.3)	8 (40)
	Never-smoker, No. (%)	19 (54.3)	10 (50)
Age of symptom onset (yr)		67.6 (6.8)	69.1 (6.9)	0.430
Age of symptom onset < 65 yr, No. (%)		11 (31.4)	6 (30)	0.912
Symptom duration (yr)		4.9 (6.0)	5.9 (4.7)	0.524
Chronic rhinosinusitis, No. (%)		18 (51.7)	8 (40)	0.514
Atopy, No. (%)		4 (44.4)	2 (22.2%)	0.620
Sputum eosinophil (%)		2.8 (4.3)	5.9 (8.2)	0.033
Sputum neutrophil (%)		35.9 (14.3)	37.9 (14.8)	0.580
Serum uric acid (mg/L)		5.2 (1.2)	6.1 (2.3)	0.083
Peripheral eosinophil count		209.6 (192.3)	197.2 (193.7)	0.819
Serum total IgE		91.0 (143.1)	130.0 (226.7)	0.616
Post-BD FVC (L)		2.75 (0.80)	2.74 (0.93)	0.974
Post-BD FVC (% pred.)		96.4 (16.0)	104.8 (15.8)	0.068
Post-BD FEV1 (L)		1.80 (0.44)	1.51 (0.37)	0.011
Post-BD FEV1 (% pred.)		91.0 (18.7)	82.0 (15.9)	0.065
Post-BD FEV1/FVC (%)		66.3 (11.8)	56.3 (13.0)	0.007

Data were presented as mean ± standard deviation except those indicated otherwise.

BD, bronchodilator; BMI, body mass index; FEV1, forced expiratory volume in 1 second; FVC, forced vital capacity.

GSEA results (Table 2A and Fig. 2) revealed the significant enrichment of 5 gene sets in cluster 1 (OXIDATIVE_PHOSPHORYLATION [OXPHOS], UNFOLDED_PROTEIN_RESPONSE [UPR], MYC_TARGETS_V1, DNA_REPAIR, and ADIPOGENESIS) and 3 gene sets in cluster 2 (EPITHELIAL_MESENCHYMAL_TRANSITION [EMT], MYOGENESIS, and KRAS_SIGNALING_DN). In the replication analysis, hierarchical clustering of gene expression patterns also revealed 2 clusters (cluster 1' with 8 EA patients and cluster 2' with 7 EA patients, Fig. 1B). Clinical characteristics of each cluster are shown in Supplementary Table S1. GSEA identified 8 gene sets (OXPHOS, ESTROGEN_RESPONSE_EARLY, E2F_TARGETS, MYC_TARGETS_V1, FATTY_ACID_METABOLISM, MYC_TARGETS_V2, ESTROGEN_RESPONSE_LATE, and DNA_REPAIR) were significantly enriched in cluster 1' and 3 gene sets (TNFA_SIGNALING_VIA_NFKB, INFLAMMATORY_RESPONSE, and EMT) in cluster 2' (Table 2B and Supplementary Fig. S2). The OXPHOS gene set was significantly enriched in clusters 1 and 1' and the EMT gene set in both clusters 2 and 2'. To assess the potential functional relevance of the replicated gene sets, the leading edge genes of gene sets enriched in both discovery and replication datasets (Supplementary Table S2) were summarized into a single metagene metric using a principal component (PC) analysis. Then, the correlations between PC1, and serum uric acid levels and post-BD FEV1/FVC values were measured in each cluster. Circulating uric acid is a major antioxidant that might help protect against oxidative stress.13 The post-BD FEV1/FVC value is an indirect indicator of airway remodeling.14 Levels of serum uric acid were lower in cluster 1 with a borderline significance (P = 0.083) compared to cluster 2, whereas post-BD FEV1/FVC values were significantly higher in cluster 1 (P = 0.008) (Supplementary Fig. S3). PC1 of leading edge genes in the OXPHOS gene set showed a negative correlation with serum uric acid levels (P = 0.075) only in cluster 1 (Fig. 3). Meanwhile, PC1 of leading edge genes in the EMT gene set showed a significantly negative correlation with post-BD FEV1/FVC values (P = 0.005) only in cluster 2 (Fig. 3). For individual gene, we found 4 genes (MRPS11, HSPA9, NDUFB4, and ACAT1) in the leading edge genes of the OXPHOS gene set and 2 genes (SNTB1 and FUCA1) in leading edge genes of the EMT gene set belonged to genes which showed more than 1.5 log 2-fold expression difference between clusters 1 and 2 with P values less than 0.01 (Fig. 4).

Table 2

Gene sets enriched significantly in each cluster

Gene set name			No. of overlapped genes	Normalized enrichment score	FDR P value
Discovery dataset
	Cluster 1
		OXPHOS	28	2.99	<0.001
		UPR	16	2.55	<0.001
		MYC_TARGETS_V1	22	1.83	0.027
		DNA_REPAIR	15	1.81	0.025
		ADIPOGENESIS	31	1.80	0.021
	Cluster 2
		EMT	31	−3.29	<0.001
		MYOGENESIS	27	−2.67	<0.001
		KRAS_SIGNALING_DN	29	−1.81	0.24
Replication dataset
	Cluster 1'
		OXPHOS	85	1.96	0.005
		ESTROGEN_RESPONSE_EARLY	24	1.89	0.005
		E2F_TARGETS	37	1.75	0.019
		MYC_TARGETS_V1	67	1.65	0.048
		FATTY_ACID_METABOLISM	38	1.64	0.041
		MYC_TARGETS_V2	28	1.64	0.034
		ESTROGEN_RESPONSE_LATE	26	1.64	0.031
		DNA_REPAIR	39	1.61	0.038
	Cluster 2'
		TNFA_SIGNALING_VIA_NFKB	27	−2.59	<0.001
		INFLAMMATORY_RESPONSE	36	−2.05	0.002
		EMT	28	−1.74	0.019

Data are shown as number (%).

FDR, false discovery rate; OXPHOS, OXIDATIVE_PHOSPHORYLATION; UPR, UNFOLDED_PROTEIN_RESPONSE; EMT, EPITHELIAL_MESENCHYMAL_TRANSITION.

Fig. 2

Gene sets enriched in each cluster identified in the discovery dataset with FDR P values less than 0.001. (A) Cluster 1, (B) Cluster 2.

FDR, false discovery rate; OXPHOS, OXIDATIVE_PHOSPHORYLATION; UPR, UNFOLDED_PROTEIN_RESPONSE, EMT, EPITHELIAL_MESENCHYMAL_TRANSITION.

Fig. 3

Association between clinical variables and PC1 of the leading edge genes from gene sets enriched in both discovery and replication datasets. (A) Cluster 1, (B) Cluster 2.

Figures without P values denote statistically insignificant associations.

PC, principal component; OXPHOS, OXIDATIVE_PHOSPHORYLATION; EMT, EPITHELIAL_MESENCHYMAL_TRANSITION.

Fig. 4

Volcano plot displaying differential expressed genes. The 4 genes (MRPS11, HSPA9, NUDF4, and ACTA1) belong to the leading edge genes of the OXPHOS gene set and 2 genes (SNTB1 and FUCA1) belong to the leading edge genes of the EMT gene set.

OXPHOS, OXIDATIVE_PHOSPHORYLATION; EMT, EPITHELIAL_MESENCHYMAL_TRANSITION.

DISCUSSION

Using gene expression profiles on sputum cells, we could identify 2 molecular clusters. Cluster 1 showed a significantly lower proportion of eosinophils in sputum and less severe airway obstruction compared to cluster 2. GSEA revealed that different biologic pathways were enriched in each cluster, implying that peculiar and discriminative pathogenesis mechanisms might exist. In replication analysis, we also identified 2 distinct molecular clusters. In addition, similar to the discovery dataset, the OXPHOS gene set was significantly enriched in one cluster and the EMT gene set was significantly enriched in the other cluster.

Approximately 70% of EA enrolled in the present study displayed symptom onset after 65 years old and the mean symptom duration was 5 years. These findings suggest that EA is not merely a prolongation of NEA but rather that aging itself may pose susceptibility to the development of asthma, as discussed elsewhere.15 In addition to anatomical changes, repetitive and long-standing exposures to environmental noxious stimuli cause the pro-inflammatory state in elderly subjects which results in aging-associated diseases in combination with oxidative stress.1617 The OXPHOS system embedded in mitochondria is the final biochemical pathway to produce adenosine triphosphate. A defect in the coupling between oxidation and phosphorylation causes various pathologic conditions including airway disorders.181920 At least 70 of the OXPHOS subunits are encoded by the nuclear genes whose expressions were measured in the present study and some defects of these genes were related with many diseases.21 The OXPHOS system is an important source of reactive oxygen species (ROS) within cells. ROS production contributes to mitochondrial damage in a range of pathologies and is also important in redox from the organelle to the rest of the cell.2223 Oxidative stress has also been linked to endoplasmic reticulum (ER) stress and to the activation of the UPR, leading to the activation of different inflammatory responses and dysregulation of the innate immune functions in the airways.2425 These prior observations match well with our findings that both the OXPHOS and UPR gene sets were significantly enriched in cluster 1. Given that previous suggestions that the OXPHOS and UPR systems have specific roles in the pathogenesis of asthma,2526 age-related changes in the OXPHOS and UPR systems may influence the development of asthma in the elderly. A decrease in serum uric acid levels in cluster 1 was also an interesting finding. As anti-oxidants may scavenge ROS, an imbalance between ROS and anti-oxidant capacity must be one of the important factors determining the effect of oxidative stress. Currently, it is not clear that patients with EA in cluster 1 have an intrinsic defect in anti-oxidant capacity. However, combined with changes in the OXPHOS and UPR systems, reduced anti-oxidant capacity may contribute to the pathogenesis of asthma in the elderly.

Cluster 2 as characterized by the enrichment of EMT gene sets. EMT describes a situation that epithelial and mesenchymal cells can, under certain conditions, switch their phenotypes27 and emerging evidence suggests that EMT is an important mechanism contributing to airway remodeling in asthma,28 which can explain lower post-BD FEV1 and post-BD FEV1/FVC ratios in cluster 2 compared to those in cluster 1. Senescent cells secrete a variety of proteins collectively known as the senescence-associated secretory phenotype (SASP), which can induce cellular plasticity and tissue change in a paracrine manner.29 Recently, it has been reported that SASP might play a causal role in the pathogenesis of chronic obstructive pulmonary disease.3031 The role of SASP in the pathogenesis of asthma in the elderly has been totally unknown. However, cellular senescence is believed to contribute to cancer progression via EMT.3233 Cellular senescence observed in the airway may induce EMT and increase the severity of asthma in the elderly.

Interestingly, we did not find any specific pathway directly related with eosinophil biology, although a significantly higher sputum eosinophil proportion was observed in cluster 2. However, given that a previous report showed that eosinophil promoted EMT of bronchial epithelial cells,34 airway eosinophilia might be complicated in the pathogenesis underlying cluster 2.

We aimed to identify biological pathways explaining a whole feature of EA, and thus we did not consider asthma-COPD overlap (ACO) as a subtype of EA in this study. If we define ACO as 1) post-BD FEV1/FVC < 0.7 and 2) classic symptoms of chronic bronchitis or suggesting signs of emphysema with chest radiography or pulmonary function test,35 the proportion of ACO is higher in cluster 2 (28.6% vs 55.5%). Although the difference showed a borderline significance (P = 0.083), the present study raised new insight in understanding ACO in the elderly. As is known, ACO is a common clinical problem in the elderly.36 Therefore, large scale studies to exam possible associations between biological pathways enriched in cluster 2 and ACO are necessary.

There are a few limitations to generalize our findings. One is the small number of participants. However, our results were replicated in an independent data set. This was a cross-sectional study, and thus we cannot know whether that the pathways identified were the primary cause or the secondary effect of asthma. Secondly, we used the uric acid level as an indirect biomarker to reflect oxidative stress status as it was easily available from routine laboratory examinations. Additional biomarkers need to be measure to confirm our observations. Finally, studies to elucidate mechanisms need to be followed.

Among leading genes identified in this study, NUDFB4 and FUCA1 are of special interest. NUDFB4 is a subunit of NADH dehydrogenase and its expression increases in smoking associated severe asthmatics.3738 Similar to the present study, it has been reported that a considerable portion of elderly asthmatics were former or current smokers.45 FUCA1 codes α-L-fucosidase that can block wound repair of primary airway epithelium.39 As well known, a delay in airway epithelial damage is associated with airway remodeling.40

In this study, we identified 2 distinct molecular clusters using gene expression profiles in sputum cells from patients with EA. The OXOPHOS and UPR gene sets were significantly enriched in one cluster with lower serum levels of uric acid. The EMT gene sets were significantly enriched in one cluster with airway remodeling defined by lower post-BD FEV/FVC ratios. This is the first study to show biologic pathways possibly related to EA pathogenesis. We postulate that oxidative stress and cellular senescence-associated with aging may be important in the development or progress of EA.