Identification of novel candidate genes involved in the progression of emphysema by bioinformatic methods

Introduction

COPD, characterized by persistent respiratory symptoms and airflow limitation, is the third leading cause of mortality worldwide.¹ Airflow limitation is mainly due to small airway obstruction and emphysema, which have distinct physiopathologic mechanisms.^2,3 Most patients with COPD have pathological alterations of both emphysema and small airway obstruction, while some have only one or no obvious change.⁴ Therefore, the two pathological phenotypes are regarded as potential subtypes of COPD.⁵

Contrary to the feature of small airway remodeling, emphysema is due to decreased deposition and excessive destruction of extracellular matrix, leading to loss of alveolar septum and attachment.^6,7 However, many studies show that the pathogenesis and progression mechanism of emphysema are complex and heterogeneous, which need to be further elucidated.^6,8,9

As a noninvasive tool to measure morphological indices, quantitative computed tomography (QCT) is an effective approach to determine the severity of COPD and distinguishing the above subtypes.¹⁰ Its assessment of emphysema has been demonstrated to be reliable, correlating well with indices of lung function, microscopic manifestations of emphysema, and clinical status of COPD patients. In addition, its assessment of small airway obstruction is also well associated with FEV₁%.^11,12

After searching the Gene Expression Omnibus (GEO) database, which is one of the largest gene expression databases in the world, we found the original gene expression profile GSE76925 with records of QCT index.¹³ To investigate the inherent molecular mechanisms in emphysema subtype of COPD, by using several bioinformatics methods,^14–16 we constructed the interacting network of differentially expressed genes (DEGs) in this profile and known COPD driver genes to identify novel candidate genes promoting progression of emphysema.

Materials and methods

Acquisition of microarray data

The GEO database (http://www.ncbi.nlm.nih.gov/geo, May 18, 2017) was retrieved to obtain gene expression profiles of lung tissues of COPD patients. The dataset GSE76925, the only one with QCT indices, was downloaded.¹³ Tests for these surgically resected lung tissue samples in GSE76925 dataset were performed using the GPL10558 platform, Illumina HumanHT-12 V4.0 expression beadchip.

Group division and statistical analysis

After screening samples’ phenotype information (from GSM2040796 to GSM2040942), samples without records of percents of low attenuation areas <−950 Hounsfield unit on inspiratory CT (%LAA-950) were ruled out. Based on the value of %LAA-950, we divided the remaining samples into two groups, severe emphysema group (%LAA-950>50%, n=11) and mild emphysema group (%LAA-950<20%, n=12).

All the continuous variables were expressed as mean ± standard deviation, and t-tests were applied to make comparison between the two groups. The categorical variables were described by constituent ratio and analyzed by Pearson chi-squared test. All statistical analyses were performed using GraphPad Prism 7 (GraphPad Software Inc, La Jolla, CA, USA). A two-side P<0.05 was considered to be statistically significant.

Identification of DEGs between severe and mild emphysema groups

To explore the underlying genes, we filtered DEGs between severe and mild emphysema groups, using GeneSpring GX software v11.5 (Agilent technologies, Santa Clara, CA, USA) at the cutoff value of corrected P-value <0.05 and |Fold Change|>1.3. We annotated them with Gene Oncology by manually retrieving Gene database (http://www.ncbi.nlm.nih.gov/gene, July 16, 2017) and roughly classified them according to the section of biological process in Gene Oncology¹⁷ by retrieving the Database for Annotation, Visualization and Integrated Discovery (DAVID)¹⁸ v6.8 (https://david.ncifcrf.gov/, October 9, 2018).

Retrieval of COPD driver genes

There has been a variety of known COPD-related genes in Global Initiative for Chronic Obstructive Pulmonary Disease (GOLD) guideline,³ peer-reviewed literatures,^19–23 Online Mendelian Inheritance (OMIM) database²⁴ (https://www.ncbi.nlm.nih.gov/omim/, July 4, 2017), and Genetic Association Database²⁵ (GAD, http://geneticassociationdb.nih.gov/). The GOLD guideline illustrated some mainstream mechanisms of COPD and emphysema, such as protease–antiprotease imbalance, which guided us to further search for specific genes in some canonical reviews. In addition, OMIM and GAD are open access databases, providing a comprehensive and authoritative compendium of genetic alterations associated with disease phenotypes. Based on the keywords of COPD or emphysema, we retrieved the above literatures and databases and identified driver genes of COPD.

Direct protein–protein interaction network of DEGs and known driver genes

The topological and functional analysis of protein interaction network is helpful in the identification of key genes and functional modules that participate in disease onset and progression.¹⁶ In network pharmacology, merging the interaction networks of drug predicted targets and driver genes of disease is an effective and original method to identify the concrete genes or pathways by which drug affects the disease.^15,16 Enlightened by this analytical method, we tried to analyze the interacted relationship between DEGs and accepted mechanism of COPD in order to identify more credible DEGs participating in emphysema development.

STRING v10.5²⁶ (https://www.string-db.org/, July 20, 2017), a web database recording physical and functional protein–protein interaction (PPi) information, was used to predict the interacted relationship between driver genes and DEGs. A variety of active interaction sources in STRING were included into our search strategy, such as text mining, experiment record, database record, coexpression, neighborhood, gene fusion, and co-occurrence. The interaction network was further visualized by Cytoscape²⁷ v3.6.1 which is an open access software aimed at annotating and visualizing biological pathways and molecular interaction networks.

Indirect PPi of DEGs and known driver genes

Weighed protein–protein interaction network analysis has been regarded as a novel approach to highlight key functional genes of complex disorders like frontotemporal dementia.¹⁴ It indicates that analyzing disease-spectrum genes, also known as first-layer interacting proteins of key genes, is a greatly potential approach to validate previous findings and explore novel disease-related mechanisms.

Thus, we retrieved STRING database v10.5 to obtain the first-layer proteins associated with DEGs products and constructed an indirect PPi of these proteins. The first-layer interacting proteins were roughly classified according to the clustering annotation of Gene Oncology and Kyoto Encyclopedia of Genes and Genomes²⁸ by the Functional Annotation tool in the DAVID database. Then, the merge tool of Cytoscape software v3.6.1²⁷ was applied to merge the indirect PPi with driver genes’ PPi to discover the interconnected and intersected functional modules and target the core genes.

In addition, highly connected nodes with a great number of edges in the network are likely to be significantly functional in the disease context and defined as hub genes.²⁹ The number of each gene node’s edges in the indirect PPi network was ranked to identify hub genes with functional significance in emphysema by Cytoscape software.

Identification of candidate transcription factors

TRANSFAC^® Professional database³⁰ is an authoritative and paid database, recording comprehensive information of transcription factor (TF), their regulated genes and binding sites prediction profiles. We performed the TF prediction of core genes by using Gene Radar tool on the GCBI website (Genminix Informatics Ltd., Shanghai, China). Based on all transcripts of each gene (Ensembl database GRCh38 version), the Gene Radar tool could acquire comprehensive TF prediction results from the TRANSFAC Professional database. In addition, Gene Radar tool could screen out high-recommended TFs by integrating the scores from the TRANSFAC database, the existence of single-nucleotide polymorphism (SNP) loci and methylation modification in TF binding sites. Therefore, we identified the candidate TFs of core genes with high recommendation grade.

Results

Baseline characteristic between severe and mild emphysema groups

As Table 1 shows, all patients were former smokers and presented with severe to very severe airflow limitation according to the GOLD guideline.³ Despite relatively small sample size, a significant difference of many characteristics between two groups was observed, like the ratio of FEV₁/FVC and body mass index, proving the credibility of %LAA-950-dependent grouping method.

Table 1 Comparisons of baseline characteristics suggested the credibility of %LAA-950-dependent grouping method Abbreviations: BMI, body mass index; %LAA-950, percents of low attenuation areas <−950 Hounsfield unit on inspiratory CT.

DEGs between two groups and the list of COPD driver genes

We identified 57 DEGs including 15 upregulated genes, 30 downregulated genes, and 12 pseudogenes (unlisted) in severe emphysema group, compared with the mild emphysema group (shown in Table 2). The Gene Oncology annotations of 45 genes were shown in Table S1.

Table 2 Forty-five DEGs were identified between severe and mild emphysema groups Note: DEGs were roughly classified according to the BP and MF terms of Gene Oncology by using the Functional Annotation tool in the DAVID database. Abbreviations: BP, biological process; DAVID, the Database for Annotation, Visualization and Integrated Discovery; DEGs, differentially expressed genes; MF, molecular function.

According to involved pathways, 135 retrieved COPD driver genes were separately placed in extracellular matrix-associated column, oxidative stress column, inflammation column, and others column (shown in Table 3).

Table 3 Known driver genes of COPD as grouped into four categories Abbreviation: ECM, extracellular matrix.

Candidate genes directly interacted with driver genes

A total of 180 genes (45 DEGs +135 driver genes) were recruited to construct the network, 7 were withdrawn for failed identification of gene symbol and 147 were found to have interaction with others. Eight of the 45 DEGs were found to have interaction relationship with driver genes: G-protein coupling receptor 65 (GPR65), Neuropeptide S receptor 1 (NPSR1), purinergic receptor P2RY13, RhoGEF and GTPase activating protein (BCR), G protein subunit β4 (GNB4), BCL2-associated athanogene 4 (BAG4), inosine monophosphate dehydrogenase 2 (IMDPH2), and Hsp40 member 14 (DNAJB14; shown in Figure 1).

Figure 1 Candidate genes were screened by direct PPi of DEGs and COPD driver genes. Notes: The left panel shows the 147 genes-constructed interaction network. Eight driver genes-associated DEGs are highlighted at the center of circle and the red lines identify the interaction relationship of the eight highlighted DEGs and the corresponding driver genes. The right panel amplifies the mutual relationship of eight DEGs and their interacted COPD driver genes. Abbreviations: DEGs, differentially expressed genes; PPi, protein–protein interaction.

A relatively separate interacting set was composed of GPR65, NPSR1, P2RY13, and GNB4. In addition, BCR and BAG4, IMDPH2 and DNAJB14 had separately bilateral relationship. When the cutoff value of the combined interaction score was set at 0.9, we found that GNB4 and P2RY13 mostly interacted with chemokines and chemokine receptors, such as CXCR1, CCR2, CXCR2, IL8, CXCR3, CCR5, CCL5, and CCR6. BAG4 interacted with TNFα and its receptor as well as heat shock protein (HSP) family. In addition, PIK3CA and PIK3R1 may play an important role by interacting with GPR65, GNB4, BCR, NPSR1, and BAG4.

Common key genes and their TFs filtered by merging of indirect PPi and driver PPi

A total of 422 first-layer interacting proteins were attained by retrieving STRING database v10.5 (shown in Table S2). Among them, 375 proteins were recruited to construct the indirect PPi and the remaining proteins were withdrawn due to failed identification or isolation from interaction network. According to the number of each node’s edges in the topological network, 375 proteins in indirect PPi were ranked and the top 20 are shown in Table S3. PIK3CA, TP53, and MAPK1, the top three genes in the rank of topological nodes of indirect PPi, had separately 86, 74, and 72 interacting nodes, which showed their potentially predominant and interconnected roles in the mechanism of emphysema progression.

The merged network illustrated in Figure 2 shows a total of 10 genes that constituted the intersection network of the two networks: TP53, IL8, CCR2, CXCR2, PIK3CA, ELANE, HSPA1A, HSPA1B, HSPA1L, and ADRB2.

Figure 2 Common key genes are screened by merging of indirect PPi and driver PPi. Notes: The left panel represents the 125 driver genes-constructed interaction network, and the right panel shows the 375 genes-constructed network of first-layer proteins of DEGs. The top three genes in the rank of topological network node stand out at the center of right circle and the red lines identify their interaction relationship with other first-layer proteins. The upper panel shows the merge network of the lower two PPis, representing the common genes and pathways involved in the two networks. Abbreviations: DEGs, differentially expressed genes; PPi, protein–protein interaction.

TFs that could bind to promoter region of the above eight genes were retrieved and shown in Table S4. Because ADRB2 was independent from the network and none of the TFs with high recommendation score was retrieved for HSPA1B, they were omitted for presentation. What’s more, SPIB, CPBP, SATB1, ZNF333, HOXA13, KID3, SOX4, and FOXO1A were potentially meaningful TFs, which could regulate no less than half genes of the above nine genes.

Discussion

We identified eight novel candidate genes (GPR65, GNB4, P2RY13, NPSR1, BCR, BAG4, IMPDH2, and TP53) promoting the progression of emphysema by means of network analysis of DEGs and COPD driver genes.

This is the first study that QCT index was applied to classify emphysema for analyzing DEGs, and known COPD driver genes were retrieved to construct interacting networks with DEGs. Our method of direct and indirect network analysis has some merit. For analysis of DEGs, it is a difficult problem to interpret the biological role of the identified single gene in the pathogenic mechanism. Performing external and experimental validation for all DEGs is cumbersome and inefficient. Incorporating driver genes into direct network analysis with DEGs is helpful in quickly highlighting causative DEGs and excluding random DEGs caused by covariates, making the role of identified DEGs more credible. In addition, protein function is regulated not only at transcriptional level but also at posttranscription level which DNA microarray could not detect. A previous study of breast cancer demonstrated that known driver genes, with their expression profiles not changed, were still capable of interconnecting many transcriptionally dysregulated genes in the protein interacting network.³¹ Therefore, we innovatively used the method of first-layer protein interacting network to further explore the indirect effects of DEGs and seek potential ignored genes. Furthermore, for the polygenic complex disease, a single gene is incapable of comprehensively illustrating the molecular mechanism of phenotypes. By merging the indirect PPi and driver PPi, we could efficiently extract the candidate protein networks involved in a specific disease phenotype. This method,¹⁴ of which the efficacy has been confirmed in a study of frontotemporal dementia, could be applied in exploring other complex disorders or extended to other phenotypes of COPD, such as airway remodeling. By comparing difference of the critical protein interactome in different phenotypes, we could unveil the different molecular mechanisms promoting complex pathological processes, which was crucial to promote biomarker and drug discovery.

As a proton-sensing receptor, GPR65 could regulate the immune response of T cells and macrophages and induce the production of MMP3 in the acidic microenvironment.^32–34 Asthma, another chronic airway disease with obstructive airflow limitation, was demonstrated to have local acidic microenvironment,³⁵ where eosinophil showed decreased apoptosis and increased viability in a GPR65-depedent manner.³⁶

The SNP of NPSR1 was associated with the decline of FEV₁ after adjusting for covariates in normal aging population.³⁷ Moreover, DNA methylation status of NPSR1 in adult severe asthma population and childhood allergic asthma population was distinct from that of control population.³⁸ NPSR1’s expression on peripheral blood eosinophils was positively correlated with asthma’s severity and serum IgE level.³⁹

Asthma and COPD have many common traits in terms of risk factors, inflammatory responses, clinical features, and therapeutic methods.^3,40 Furthermore, the role of eosinophils in the pathogenesis and treatment of COPD is gradually recognized.^3,41 Therefore, we speculated that the above genes related to asthma were highly likely to be involved in the pathogenesis of COPD and emphysema.

As an extracellular ADP receptor, P2RY13 participated in purinergic signaling pathway, resulting in the apoptosis of pancreatic β-cells⁴² and differentiation of marrow stem cells into osteoblasts.⁴³ Since the roles of extracellular adenosine ATP and its receptor P2RX in COPD have been confirmed,^44,45 ADP, the intermediate in purinergic metabolic pathways, may also have pathogenic effects on COPD.

In addition, common key genes identified by the indirect method matched well with the two-hit hypothesis of COPD,⁴⁶ especially the part of senescence and senescence-associated secretary phenotype (SASP).⁴⁷ Senescence is an irreversible cell state, at which a cell is deprived of its replicative capacity with cell cycle arrest.⁴⁸ The p53 (encoded by TP53)/p21 pathway participated in all types of senescence mechanisms, arresting cell cycle at the G1/S and G2/M check points.⁴⁹ SASP refers to the alteration of aging cell’s secretome toward more production of proinflammatory cytokines, including IL-8 and monocyte chemotactic protein 1 (MCP-1).⁴⁹ IL-8 and its receptor CXCR2, with neutrophil chemotactic ability, are just one of the most important chemokine-receptor pairs in COPD pathogenesis, as well as MCP-1 encoded by CCR2.⁶ Moreover, phosphoinositide 3 kinase (PI3K), the product of PIK3CA, was also known as a pro-senescent kinase by inactivating HDAC-2 which is an antiaging molecule, because knockdown of HDAC-2 could induce cellular senescence by enhancing p53-dependent transcriptional responses.⁵⁰

Based on these evidence, we hypothesized that TP53 might play a central role in promoting progression of emphysema. Firstly, beside IL-8, CXCR2 and CCR2, elastase, the products of ELANE, and PI3K are also well-recognized COPD driver genes playing an important role in protease–antiprotease imbalance and chronic inflammation of COPD.⁶ The involvement of these genes supports our results and in return proves the role of TP53.

Secondly, TP53 can induce cell cycle arrest, apoptosis, senescence, DNA repair, or metabolic alterations, in response to oxidative stress and DNA damage.⁵¹ Some studies confirmed that TP53 was overexpressed in the emphysematous lung tissue.⁵² A Genome-Wide Association Study for 365 patients with emphysema proved the association of TP53’s SNP with apoptotic signaling and smoking-related emphysematous changes in smoker’s lungs.⁵³ Furthermore, a RNA-sequencing study of COPD patients’ lung tissues identified the enrichment of p53/hypoxia pathway and the phenomenon of much frequent molecule’s alternative splicing in this pathway.⁵⁴

Thirdly, the role of TP53 in senescence might reveal its effects in COPD. Many evidences have shown the association between senescence and pathogenesis of COPD. Cellular experiments proved that alveolar epithelial and endothelial cell as well as fibroblast underwent accelerated senescence in emphysematous lung.^55,56 Epidemiological surveys indicated that the incidence of COPD and the decline of FEV₁ increased with growth of age.³ Moreover, airway and parenchyma of the patients with COPD and healthy senior citizens had similar structural changes.^50,57

There are many studies searching for key genes associated with emphysema. In a research on seeking differently expressed miRNAs of emphysema, the miR-638 was identified as an effector molecule and it could regulate accelerated senescence, which was partially consistent with our hypothesis.⁵⁸ However, we did not reproduce the DEGs of other studies for emphysema. On one hand, it was due to different grouping methods;⁵⁹ on the other hand, their samples mainly came from patients with moderate COPD (FEV₁% was about 60%),^60,61 so their results mainly explained the early mechanisms of emphysema progression.

Our study has some limitations. Firstly, the sample size is relatively small,⁶² which is due to limited numbers of accessible datasets in GEO database. Secondly, the selection of COPD driver genes is potentially biased and incomplete so that some meaningful DEGs may be ignored. Thirdly, PPi prediction has false positives and false negatives. The web tool STRING v10.5 defines PPi by the standard of text mining, experiment record, database record, coexpression, neighborhood, gene fusion, and co-occurrence, which may have a bit of controversy. In addition, interactions proved by experiments in vitro may also have differences compared with those in vivo.

Despite these limitations, our study has put forward some novel candidate genes, and following experiments or larger databases are needed to testify the role of the above candidate genes in the mechanism of emphysema progression.

Conclusion

We have identified several novel candidate genes promoting emphysema, like GPR65, NPSR1, and TP53, which may be helpful in filling in the gap of knowledge in the field of COPD.

Acknowledgments

The authors would like to acknowledge Dr Xiao Shi who critically reviewed the article for English writing. This study was supported by the National Key Research and Development Program of China (grant Nos 2017YFC1309303 and 2017YFC1309300) and the National Natural Science Foundation of China (grant Nos 81670030 and 81470231).

Disclosure

The authors report no conflicts of interest in this work.

---

References

1.		Global Health Estimates 2016: Deaths by Cause, Age, Sex, by Country and by Region, 2000–2016. Geneva: World Health Organization; 2018.
2.		Hogg JC, Timens W. The pathology of chronic obstructive pulmonary disease. Annu Rev Pathol. 2009;4:435–459.
3.		Global Initiative for Chronic Obstructive Pulmonary Disease: Global Strategy for the Diagnosis, Management and Prevention of COPD; revised; 2017. Available from: https://goldcopd.org/gold-2017-global-strategy-diagnosis-management-prevention-copd/. Accessed April 4, 2017.
4.		Subramanian DR, Gupta S, Burggraf D, et al. Emphysema- and airway-dominant COPD phenotypes defined by standardised quantitative computed tomography. Eur Respir J. 2016;48(1):92–103.
5.		Ziegler-Heitbrock L, Frankenberger M, Heimbeck I, et al. The EvA study: aims and strategy. Eur Respir J. 2012;40(4):823–829.
6.		Chung KF, Adcock IM. Multifaceted mechanisms in COPD: inflammation, immunity, and tissue repair and destruction. Eur Respir J. 2008;31(6):1334–1356.
7.		Vermylen JH, Kalhan R. Revealing the complexity of chronic obstructive pulmonary disease. Transl Res. 2013;162(4):203–207.
8.		Gharib SA, Manicone AM, Parks WC. Matrix metalloproteinases in emphysema. Matrix Biol. 2018;73:34–51.
9.		Tuder RM, Petrache I. Pathogenesis of chronic obstructive pulmonary disease. J Clin Invest. 2012;122(8):2749–2755.
10.		Lynch DA, Al-Qaisi MA. Quantitative computed tomography in chronic obstructive pulmonary disease. J Thorac Imaging. 2013;28(5):284–290.
11.		Nambu A, Zach J, Schroeder J, et al. Quantitative computed tomography measurements to evaluate airway disease in chronic obstructive pulmonary disease: relationship to physiological measurements, clinical index and visual assessment of airway disease. Eur J Radiol. 2016;85(11):2144–2151.
12.		Castaldi PJ, San José Estépar R, Mendoza CS, et al. Distinct quantitative computed tomography emphysema patterns are associated with physiology and function in smokers. Am J Respir Crit Care Med. 2013;188(9):1083–1090.
13.		Morrow JD, Zhou X, Lao T, et al. Functional interactors of three genome-wide association study genes are differentially expressed in severe chronic obstructive pulmonary disease lung tissue. Sci Rep. 2017;7:44232.
14.		Ferrari R, Lovering RC, Hardy J, Lewis PA, Manzoni C. Weighted protein interaction network analysis of frontotemporal dementia. J Proteome Res. 2017;16(2):999–1013.
15.		Li J, Zhao P, Li Y, Tian Y, Wang Y. Systems pharmacology-based dissection of mechanisms of Chinese medicinal formula Bufei Yishen as an effective treatment for chronic obstructive pulmonary disease. Sci Rep. 2015;5:15290.
16.		Zanzoni A, Soler-López M, Aloy P. A network medicine approach to human disease. FEBS Lett. 2009;583(11):1759–1765.
17.		Huntley RP, Sawford T, Mutowo-Meullenet P, et al. The GOA database: gene Ontology annotation updates for 2015. Nucleic Acids Res. 2015;43(Database issue):D1057–D1063.
18.		Huang da W, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4(1):44–57.
19.		Marciniak SJ, Lomas DA. Genetic susceptibility. Clin Chest Med. 2014;35(1):29–38.
20.		Barnes PJ. The cytokine network in chronic obstructive pulmonary disease. Am J Respir Cell Mol Biol. 2009;41(6):631–638.
21.		Mcguinness AJ, Sapey E. Oxidative stress in COPD: sources, markers, and potential mechanisms. J Clin Med. 2017;6(2):E21.
22.		Brusselle GG, Joos GF, Bracke KR. New insights into the immunology of chronic obstructive pulmonary disease. Lancet. 2011;378(9795):1015–1026.
23.		Burgstaller G, Oehrle B, Gerckens M, et al. The instructive extracellular matrix of the lung: basic composition and alterations in chronic lung disease. Eur Respir J. 2017;50(1):1601805.
24.		Mckusick VA. Mendelian inheritance in man and its online version, OMIM. Am J Hum Genet. 2007;80(4):588–604.
25.		Becker KG, Barnes KC, Bright TJ, Wang SA. The genetic association database. Nat Genet. 2004;36(5):431–432.
26.		Szklarczyk D, Morris JH, Cook H, et al. The STRING database in 2017: quality-controlled protein-protein association networks, made broadly accessible. Nucleic Acids Res. 2017;45(D1):D362–D368.
27.		Shannon P, Markiel A, Ozier O, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–2504.
28.		Kanehisa M, Furumichi M, Tanabe M, Sato Y, Morishima K. KEGG: new perspectives on genomes, pathways, diseases and drugs. Nucleic Acids Res. 2017;45(D1):D353–D361.
29.		Puig-Butille JA, Gimenez-Xavier P, Visconti A, et al. Genomic expression differences between cutaneous cells from red hair color individuals and black hair color individuals based on bioinformatic analysis. Oncotarget. 2017;8(7):11589–11599.
30.		Matys V, Fricke E, Geffers R, et al. TRANSFAC: transcriptional regulation, from patterns to profiles. Nucleic Acids Res. 2003;31(1):374–378.
31.		Chuang HY, Lee E, Liu YT, Lee D, Ideker T. Network-based classification of breast cancer metastasis. Mol Syst Biol. 2007;3:140.
32.		Onozawa Y, Fujita Y, Kuwabara H, et al. Activation of T cell death-associated gene 8 regulates the cytokine production of T cells and macrophages in vitro. Eur J Pharmacol. 2012;683(1–3):325–331.
33.		Radu CG, Nijagal A, McLaughlin J, Wang L, Witte ON. Differential proton sensitivity of related G protein-coupled receptors T cell death-associated gene 8 and G2A expressed in immune cells. Proc Natl Acad Sci U S A. 2005;102(5):1632–1637.
34.		Xu H, Chen X, Huang J, et al. Identification of GPR65, a novel regulator of matrix metalloproteinases using high through-put screening. Biochem Biophys Res Commun. 2013;436(1):96–103.
35.		Hunt JF, Fang K, Malik R, et al. Endogenous airway acidification. Implications for asthma pathophysiology. Am J Respir Crit Care Med. 2000;161(3 Pt 1):694–699.
36.		Kottyan LC, Collier AR, Cao KH, et al. Eosinophil viability is increased by acidic pH in a cAMP- and GPR65-dependent manner. Blood. 2009;114(13):2774–2782.
37.		Poon AH, Houseman EA, Ryan L, et al. Variants of asthma and chronic obstructive pulmonary disease genes and lung function decline in aging. J Gerontol A Biol Sci Med Sci. 2014;69(7):907–913.
38.		Reinius LE, Gref A, Sääf A, et al. DNA methylation in the Neuropeptide S Receptor 1 (NPSR1) promoter in relation to asthma and environmental factors. PLoS One. 2013;8(1):e53877.
39.		Ilmarinen P, James A, Moilanen E, et al. Enhanced expression of neuropeptide S (NPS) receptor in eosinophils from severe asthmatics and subjects with total IgE above 100 IU/ml. Peptides. 2014;51:100–109.
40.		Silva GE, Sherrill DL, Guerra S, Barbee RA. Asthma as a risk factor for COPD in a longitudinal study. Chest. 2004;126(1):59–65.
41.		Vedel-Krogh S, Nielsen SF, Lange P, Vestbo J, Nordestgaard BG. Blood eosinophils and exacerbations in chronic obstructive pulmonary disease. The Copenhagen General Population Study. Am J Respir Crit Care Med. 2016;193(9):965–974.
42.		Tan C, Salehi A, Svensson S, Olde B, Erlinge D. ADP receptor P2Y(13) induce apoptosis in pancreatic beta-cells. Cell Mol Life Sci. 2010;67(3):445–453.
43.		Biver G, Wang N, Gartland A, et al. Role of the P2Y13 receptor in the differentiation of bone marrow stromal cells into osteoblasts and adipocytes. Stem Cells. 2013;31(12):2747–2758.
44.		Polosa R. Adenosine-receptor subtypes: their relevance to adenosine-mediated responses in asthma and chronic obstructive pulmonary disease. Eur Respir J. 2002;20(2):488–496.
45.		Pelleg A, Schulman ES, Barnes PJ. Extracellular adenosine 5′-triphosphate in obstructive airway diseases. Chest. 2016;150(4):908–915.
46.		Aoshiba K, Tsuji T, Yamaguchi K, Itoh M, Nakamura H. The danger signal plus DNA damage two-hit hypothesis for chronic inflammation in COPD. Eur Respir J. 2013;42(6):1689–1695.
47.		Kumar M, Seeger W, Voswinckel R. Senescence-associated secretory phenotype and its possible role in chronic obstructive pulmonary disease. Am J Respir Cell Mol Biol. 2014;51(3):323–333.
48.		Campisi J, D’Adda di Fagagna F. Cellular senescence: when bad things happen to good cells. Nat Rev Mol Cell Biol. 2007;8(9):729–740.
49.		Abbadie C, Pluquet O, Pourtier A. Epithelial cell senescence: an adaptive response to pre-carcinogenic stresses? Cell Mol Life Sci. 2017;74(24):4471–4509.
50.		Ito K, Barnes PJ. COPD as a disease of accelerated lung aging. Chest. 2009;135(1):173–180.
51.		Vousden KH, Prives C. Blinded by the light: the growing complexity of p53. Cell. 2009;137(3):413–431.
52.		Morissette MC, Vachon-Beaudoin G, Parent J, Chakir J, Milot J. Increased p53 level, Bax/Bcl-x(L) ratio, and TRAIL receptor expression in human emphysema. Am J Respir Crit Care Med. 2008;178(3):240–247.
53.		Mizuno S, Ishizaki T, Kadowaki M, et al. p53 Signaling pathway polymorphisms associated with emphysematous changes in patients with COPD. Chest. 2017;152(1):58–69.
54.		Kusko RL, Brothers JF, Tedrow J, et al. Integrated genomics reveals convergent transcriptomic networks underlying chronic obstructive pulmonary disease and idiopathic pulmonary fibrosis. Am J Respir Crit Care Med. 2016;194(8):948–960.
55.		Tsuji T, Aoshiba K, Nagai A. Alveolar cell senescence in patients with pulmonary emphysema. Am J Respir Crit Care Med. 2006;174(8):886–893.
56.		Dagouassat M, Gagliolo JM, Chrusciel S, et al. The cyclooxygenase-2-prostaglandin E2 pathway maintains senescence of chronic obstructive pulmonary disease fibroblasts. Am J Respir Crit Care Med. 2013;187(7):703–714.
57.		Mercado N, Ito K, Barnes PJ. Accelerated ageing of the lung in COPD: new concepts. Thorax. 2015;70(5):482–489.
58.		Christenson SA, Brandsma CA, Campbell JD, et al. miR-638 regulates gene expression networks associated with emphysematous lung destruction. Genome Med. 2013;5(12):114.
59.		Spira A, Beane J, Pinto-Plata V, et al. Gene expression profiling of human lung tissue from smokers with severe emphysema. Am J Respir Cell Mol Biol. 2004;31(6):601–610.
60.		Faner R, Cruz T, Casserras T, et al. Network analysis of lung transcriptomics reveals a distinct B-cell signature in emphysema. Am J Respir Crit Care Med. 2016;193(11):1242–1253.
61.		Francis SM, Larsen JE, Pavey SJ, et al. Expression profiling identifies genes involved in emphysema severity. Respir Res. 2009;10:81.
62.		Hart SN, Therneau TM, Zhang Y, Poland GA, Kocher JP. Calculating sample size estimates for RNA sequencing data. J Comput Biol. 2013;20(12):970–978.

Supplementary materials

Table S1 GO annotation of DEGs between severe and mild emphysema groups Note: Some COPD-associated genes are highlighted in bold. Abbreviations: DEGs, differentially expressed genes; GO annotation, Gene Oncology annotation.

Table S2 The list of first-layer interacting proteins associated with DEGs between two groups Note: The first-layer interacting proteins were roughly classified according to the BP terms of Gene Oncology and KEGG by the Functional Annotation tool in the DAVID database. Abbreviations: BP, biological process; DAVID, the Database for Annotation, Visualization and Integrated Discovery; DEGs, differentially expressed genes; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Table S3 Top 20 topological network nodes in first-layer’s PPi Abbreviations: GO, Gene Oncology; PPi, protein–protein interaction.

Table S4 Transcript factors of common key genes screened by indirect PPi Note: Key genes are highlighted in bold. Abbreviation: PPi, protein–protein interaction.