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

PURPOSE: By reanalyzing the gene expression profile GSE76925 in the Gene Expression Omnibus database using bioinformatic methods, we attempted to identify novel candidate genes promoting the development of emphysema in patients with COPD. PATIENTS AND METHODS: According to the Quantitative CT data in GSE76925, patients were divided into mild emphysema group (%LAA-950<20%, n=12) and severe emphysema group (%LAA-950>50%, n=11). Differentially expressed genes (DEGs) were identified using Agilent GeneSpring GX v11.5 (corrected P-value <0.05 and |Fold Change|>1.3). Known driver genes of COPD were acquired by mining literatures and retrieving databases. Direct protein-protein interaction network (PPi) of DEGs and known driver genes was constructed by STRING.org to screen the DEGs directly interacting with driver genes. In addition, we used STRING.org to obtain the first-layer proteins interacting with DEGs' products and constructed the indirect PPi of these interaction proteins. By merging the indirect PPi with driver genes' PPi using Cytoscape v3.6.1, we attempted to discover potential pathways promoting emphysema's development.
RESULTS: All the patients had COPD with severe airflow limitation (age=62±8, FEV1%=28±12). A total of 57 DEGs (including 12 pseudogenes) and 135 known driving genes were identified. Direct PPi suggested that GPR65, GNB4, P2RY13, NPSR1, BCR, BAG4, and IMPDH2 were potential pathogenic genes. GPR65 could regulate the response of immune cells to the acidic microenvironment, and NPSR1's expression on eosinophils was associated with asthma's severity and IgE level. Indirect merging PPi demonstrated that the interacting network of TP53, IL8, CCR2, HSPA1A, ELANE, PIK3CA was associated with the development of emphysema. IL8, ELANE, and PIK3CA were molecules involved in the pathological mechanisms of emphysema, which also in return proved the role of TP53 in emphysema.
CONCLUSION: Candidate genes such as GPR65, NPSR1, and TP53 may be involved in the progression of emphysema.

Introduction

COPD, characterized by persistent respiratory symptoms and airflow limitation, is the third leading cause of mortality worldwide.1 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.4 Therefore, the two pathological phenotypes are regarded as potential subtypes of COPD.5

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.10 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 FEV1%.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.13 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.13 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 Oncology17 by retrieving the Database for Annotation, Visualization and Integrated Discovery (DAVID)18 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,3 peer-reviewed literatures,19–23 Online Mendelian Inheritance (OMIM) database24 (https://www.ncbi.nlm.nih.gov/omim/, July 4, 2017), and Genetic Association Database25 (GAD, http://geneticassociationdbnih.gov/). The GOLD guideline illustrated some mainstream mecha nisms 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.16 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.526 (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, neighbor-hood, gene fusion, and co-occurrence. The interaction network was further visualized by Cytoscape27 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.14 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 Genomes28 by the Functional Annotation tool in the DAVID database. Then, the merge tool of Cytoscape software v3.6.127 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.29 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 database30 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.3 Despite relatively small sample size, a significant difference of many characteristics between two groups was observed, like the ratio of FEV1/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

Characteristics	Severe emphysema group (n=11)	Mild emphysema group (n=12)	P-value
Age (years)	61.6±5.9	63.1±10.4	.0.05
Male/female	8/3	6/6	.0.05
BMI (kg/m2)	23.0±3.3	27.9±5.3	0.0102
Smoking history (pack-years)	66.5±21.2	58.6±29.2	0.0013
%LAA-950	52.7±2.0	7.3±5.3	<0.0001
FEV1 (%predicted)	23.0±8.4	32.7±13.3	0.051
FEV1/FVC (%)	25.5±4.8	42.3±14.2	0.0012

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

Functional category	Gene symbol	Dysregulation	P-value	Fold change
Transcriptional regulation	KANK1PHF1PHF6TADA2ATRIM34ZFHX3ZNF322ZNF451	DownDownUpUpUpDownUpUp	5.58E–054.45E–052.08E–058.32E–053.85E–064.31E–052.88E–068.66E–05	−1.82−1.632.012.421.54−1.541.582.54
Membrane receptor andsignal pathway	BAG4BCRFYB1GNB4GPR65NPSR1NPHP4P2RY13RNF213ZFP106ZC3HAV1	UpDownUpUpUpDownDownUpUpDownDown	0.0000616.05E–055.24E–067.66E–057.58E–057.72E–058.48E–053.03E–056.32E–052.71E–058.41E–05	1.71−1.592.432.492.54−1.73−1.934.22−1.5−1.43
Metabolism	DPM3ELOVL3ETNK2IMPDH2	DownDownDownDown	7.09E–054.72E–052.63E–057.49E–05	−1.44−1.47−1.74−1.5
Cilium	IFT140TMEM80	DownDown	1.65E–055.48E–05	−1.93−1.75
Protein modification	USP33NUP58PARP16	UpUpDown	2.29E–053.89E–056.93E–05	2.642.96−1.73
Others	DNAJB14ZBTB8OSEHBP1ATP1B2FAM149ATLN1SF3A1FAM168BCYB5D2KCNJ4ZCCHC3MRPS24swi5SERPINI1SVEP1VPS28OGFOD3	UpUpDownDownDownDownDownDownDownDownDownDownDownDownDownDownDown	4.21E–056.47E–059.96E–077.85E–067.86E–062.68E–052.46E–054.84E–053.72E–058.54E–057.39E–054.46E–053.21E–057.31E–055.94E–056.99E–057.25E–05	2.161.65−1.42−3.09−3.32−1.47−1.49−1.71−1.33−1.77−1.57−1.36−1.34−1.58−1.67−1.52−1.57

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

Synthesis and degradation of ECM (n=37)		Oxidative stress (n=12)	Abnormal inflammation (n=36)		Others (n=50)
ELN	COL1A1	GSTP1	TNF	CCL5	LTA4H	CLASP1	VEGFA
FBLN4	COL1A2	GSTM1	TNFRSF1A	IL17F	MUC5AC	ADRB2	GABPA
FBLN5	FBN2	HMOX1	TNFRSF1B	IL1RN	MUC5B	GC	MTHFR
FBN1	FBN3	NOS2	IL-17A	IFNG	SLC6A4	DEFB1	SPAR
ATP7A	COL3A1	NOS3	IL-18	TSLP	EGF	CYP21A2	HSPA1B
TGFB1	COL8A1	SOD2	IL1b	CCR2	EGFR	CFTR	CAT
TGFBR3	COL4A1	MMAC1	HDAC2	IL8RB	FGF10	APOE	OGG1
LTBP4	FN4	SOD3	IL12	IL13	CHRNA3	AGTR1	PDE4D
SERPINE2	FN1	PIK3CA	IL21	IL11	CHRNA5	ADRB3	TCEAL1
ELANE	DCN	PIK3R1	IL22	CCR5	IREB2	TF	BCL2
MMP1	BGN	NFE2L2	IL-23	CCR6	FAM13A	SFTPB	DBP
TIMP1	TGFBR1	EPHX1	IL27	CXCL8	FTO	SERPINE1	HCK
MMP2	SMAD3		IL32	CXCR1	BICD1	SERPINA1
MMP3	SMAD7		IL-4	CXCR2	HHIP	NAT2
MMP8	VCAN		IL-6	CXCR3	ACE	LTA
MMP9	TNC		IL-10	TLR9	KCNIP4	HSPA1L
MMP10	SPP1		CCL11	IL8RA	CRHR1	HSPA1A
MMP14	TIMP2		CLL2	CCL1	CYP1A2	HRAS
MMP12					TP53	SCGB1A1

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.31 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,14 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,35 where eosinophil showed decreased apoptosis and increased viability in a GPR65-depedent manner.36

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

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 β-cells42 and differentiation of marrow stem cells into osteoblasts.43 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,46 especially the part of senescence and senescence-associated secretary phenotype (SASP).47 Senescence is an irreversible cell state, at which a cell is deprived of its replicative capacity with cell cycle arrest.48 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.49 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).49 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.6 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.50

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.6 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.51 Some studies confirmed that TP53 was overexpressed in the emphysematous lung tissue.52 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.53 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.54

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 FEV1 increased with growth of age.3 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.58 However, we did not reproduce the DEGs of other studies for emphysema. On one hand, it was due to different grouping methods;59 on the other hand, their samples mainly came from patients with moderate COPD (FEV1% 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,62 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.

Table S1

GO annotation of DEGs between severe and mild emphysema groups

Functional category	Gene Symbol	GO annotation
Transcriptionalregulation	KANK1PHF1PHF6TADA2ATRIM34ZFHX3ZNF322ZNF451	Positive regulation of Wnt signaling pathway, negative regulation of actin filament polymerization and so onInvolved in regulation of histone H3-K27 methylation and cellular response to DNA damage stimulusNegative regulation of transcription from RNA polymerase II promoter, an oncogeneA transcriptional activator adaptor; acetylating and destabilizing nucleosomesInvolved in interferon signaling pathwaysTranscription factor activity, RNA polymerase II distal enhancer sequence-specific bindingRegulate transcriptional activation in MAPK signaling pathwaysNegative regulation of transcription initiation from RNA polymerase II promoter, histone H3-K9 acetylation and TGF-β signaling pathway
Membrane receptor and signal pathway	BAG4BCRFYB1GNB4GPR65NPSR1NPHP4P2RY13RNF213ZFP106ZC3HAV1	Negative regulation of apoptotic process, response to TNFαGTPase activator activity and Rho guanyl-nucleotide exchange factor activityInvolved in TCR signaling pathways, the expression of IL-2 and process of NLS-bearing protein importing into nucleusA subunit of heterotrimeric guanine nucleotide-binding proteins involved in cellular response to glucagon stimulusInvolved in G-protein coupled receptor signaling pathway, actin cytoskeleton reorganization, and apoptotic processNeuropeptide and vasopressin receptor activity, increased expression in lung for asthmaInvolved in actin cytoskeleton organization, hippo signaling, and negative regulation of canonical Wnt signaling pathwayG-protein coupled purinergic nucleotide receptor and negative regulation of adenylate cyclase activity signaling pathwayATPase activity, ubiquitin-protein transferase activity, and negative regulation of noncanonical Wnt signaling pathwayInsulin receptor signaling pathwayDefense response to virus
Metabolism	DPM3ELOVL3ETNK2IMPDH2	GPI anchor biosynthetic process and protein mannosylationFatty acid elongase activity providing precursors for synthesis of sphingolipids and ceramidesA member of choline/ethanolamine kinase family that catalyses phosphatidylethanolamine biosynthetic processPurine ribonucleoside monophosphate biosynthetic process, neutrophil degranulation, and oxidation-reduced process
Cilium	IFT140TMEM80	Intraciliary transport involved in cilium assemblyIntegral component of membrane
Protein modification	USP33NUP58PARP16	Protein deubiquitination and involved in slit-dependent cell migration and beta-2 adrenergic receptor signalingA component of the nuclear pore complex playing a role of nucleocytoplasmic transporter activityNAD+ ADP-ribosyltransferase activity and protein serine/threonine kinase activator activity
Others	DNAJB14ZBTB8OSEHBP1ATP1B2FAM149ATLN1SF3A1FAM168BCYB5D2KCNJ4ZCCHC3MRPS24swi5SERPINI1SVEP1VPS28OGFOD3	Hsp70 protein binding and chaperone cofactor-dependent protein refoldingtRNA splicing via endonucleolytic cleavage and ligationEndocytosis and its mutation associated with prostate cancerATP hydrolysis coupled transmembrane transport and cell adhesionAssociated with acute mountain sicknessIntegrin-mediated signaling pathway, cell–cell and cell–substrate junction assembly such as actinA component of the mature U2 snRNP playing a role of pre-mRNA splicingMyelin-associated neurite-outgrowth inhibitorPositive regulation of neuron differentiationA member of the inward rectifier potassium channel familyRNA bindingA structural constituent of ribosome related to mitochondrial translationDNA repair protein swi5 homologSerine-type endopeptidase inhibitor activity and association with central and peripheral nervous system developmentA ligand for integrin α9β1 and involved in cell adhesionEndosomal transport, macroautophagy, negative regulation of protein ubiquitination and viral buddingOxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen

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

Translational regulation	Signaling pathway	Matrix and cell adhesion	Virus-associated genes	Others
RPL12	TP53	PIK3CA	CHMP6	KCNJ2	RBBP7	ZNF91	GPR37L1
RPL13A	PIK3CA	MAPK1	RAE1	KCNJ4	AEBP2	AP4E1	NPHP4
RPL18	MAPK1	FYN	ZC3HAV1	FOS	CYB5D2	NEURL	PIP5K1C
RPL18A	GRB2	CRKL	DDX1	GART	PARP14	TEX10	SEL1L
RPL8	SHC1	GRB2	DDX58	IFT122	RIPK4	AP4M1	TMEM67
RPS15	PDGFRB	SHC1	IRF3	MKS1	ZCCHC3	DPAGT1	ATIC
RPS3	ITGB1	PXN	IRF7	P2RY13	AFP	NEURL1B	EHD2
RPS9	ITGB3	PDGFRB	RPL12	WDR19	IMPDH1	APBB1IP	NPLOC4
MRPS10	ITGB5	ITGB1	IRF9	ACACA	PARP16	PI15	PISD
MRPS12	VWF	ZYX	RPL13A	CIDEA	RNF19B	TLX1	SELV
PAIP1	GNAI1	ITGB2	RPL18	GCC2	ZFHX3	ARC	TMEM80
MRPS14	GNB1	TLN1	RPL18A	IFT140	AGBL3	NMRAL1	EHD3
MRPS15	GNB2	ITGB3	RPL8	MOCS1	IMPDH2	TLX3	GRAP2
MRPS22	GNB3	ITGB5	RPS15	WDR35	RNF213	ARHGEF6	KIF21A
MRPS23	GNG2	VASP	RPS3	ACACB	TADA2A	ITPA	NPSR1
MRPS24	GNAO1	VCL	RPS9	COL14A1	ZFP106	NNMT	PLAT
MRPS31	GNB4	VWF	TMEM48	GM2A	AGXT2L1	PIGM	ATP1A1
MRPS33	GNG10		ATG7	IFT172	DGAT2	TMEM171	ELANE
MRPS5	GNG13		TNF	MOCS2	IQCB1	ARHGEF7	KNG1
RPL10L	GNG3	Substance transport	KPNB1	P2RY4	PCYT2	DR1	NUDT11
SEPSECS	GNG4	RAE1	NUP107	ACE2	RPE65	NOL10	PLAU
CHCHD1	GNG5	TMEM48	NUP188	CPB1	TADA3	PIGV	SERP1NI1
	GNG7	NUP107	TRIM34	GMPS	ZNF322	TMEM218	TNFRSF1A
	JAK2	NUP188	NUP205	IFT52	AKAP9	DTX3L	ATP1A2
Protein modification and folding	GNGT2	NUP205	NUP35	P2RY8	DMRT1	JARID2	ELAVL4
RAE1	GNA11	NUP35	PML	CPS1	IRF2	NOL12	GRM6
TMEM48	GNAQ	NUP62	TSG101	IFT57	RPGRIP1L	PIK3C2B	PLEKHG7
NUP107	ADRBK1	NUP93	NUP54	SUPT3H	TAF1	SASH1	SF3A1
NUP188	IL8	NUPL1	HERC5	YEATS2	ZNF385D	TMEM222	ATP1A4
NUP205	CCR2	SUMO1	NUP62	ADRA1	ALG1	EED	ELOVL3
NUP35	ABL1	SUMO2	NUP93	GNA14	DNAJB14	GOLT1A	H2AFV
TSG101	P2RY14	SUMO3	NUPL1	IFT80	MSRB1	KAT2A	LACE
NUP62	GNA15	SUMO4	OAS1	PARL	PFAS	NOL3	PLG
NUP93	SYK	RANGAP1	PPIA	SUZ12	TAF9	SCD	SF3A2
NUPL1	MYC	PAIP1	OAS2	YIPF6	ALG3	TMEM247	ATP1B1
VPS37A	PHLPP1	STAT3	OAS3	ADRB2	DOCK8	ATG3	ENAH
VPS4A	PHLPP2	KPNB1	OASL	CTPS1	PHF19	EHBP1	H2AFZ
SUMO1	GPR65	NUP54	VPS28	IFT81	TCEB2	GPCPD1	PLRG1
SUMO2	LAT	UBE2I	VPS37A	RBBP4	ZNF768	KAT2B	SF3A3
SUMO3	TSHB	NUTF2	VPS37B	SVEP1	DOLK	NPHP1	TRMT10C
SUMO4	LPAR1		VPS37C	ZBTB8OS	NBEAL1	PIP5K1A	ATP1B3
RANGAP1	LPAR2		VPS4A	IFT88	PHF6	SDCCAG8	ENG
PARP1	OXGR1			PARP10	TCTN1	EHD1	HERC1
ZNF451	FPR1			VPS36	MIB2	CEP290	SIRPA
ALG5	FPR2		DNA repair	CEPT1	OS9	OGFOD3	USP33
DOLPP1	ADCY3		SUMO1	HSPA6	SLTM	SKAP2	CAD
DPM1	CRKL		PARP1	CERS2	FCGR1B	UXT	FAM98A
NFATC2IP	CXCR2		RPS3	FXYD1	HSPA2	CCP110	HSPA14
DPM2	GNAI3		RAD51	STAT5A	MICAL1	FAM98C	MBIP
PIAS3	ARRB1		RAD51B	FXYD2	SNF8	HSPA1B	OBSCN
DPM3	ARRB2		RAD51C	FXYD6	VHL	MIB1	PVR
TP53	LYN		XRCC2	HUWE1	CDKL5	OPHN1	SKAP1
FAM125A	PXN		RAD51D	FXYD7	HSPA4	SLC20A2	USP48
FAM125B	STAT3		XRCC3	HSPA5	MKI67IP	FASN	CCDC101
GNAI1	STAT5B		SWI5	SPATA5	SPATA2	HSPA1L	FAM98B
GNAI3			PARP2	SH3GL1	LATS1	LLPH	B9D2
GNB1			ATM	C1orf177	SF3B1	POTEI	EZH2
GNB2			SFPQ	FAM155B	AWAT1	SF3B4	HERC6
GNB3			SFR1	HSD17B12	ETNK2	UBA7	HSPA13
GNG2			CDC5L	PRPF19	HERC4	HIST1H1A	MAX
GNAO1				USP20	LCP2	POTEJ	PTPRF
PPIA				C22orf28	POTEE	BCR	FAM149A
BAG4				FAM168B	SF3B2	HMG20A	HPRT1
SIL1				HSPA12B	B4GALNT1	POU1F1	LRRK2
MESDC2				PRPF6	EVL	UBR2	POTEF
HSPA8				USP22	LIAS	C14orf166	SF3B3
HSPA9				C2orf49	HSPA1A	FAM71E1	TTC21B

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

Degree	Name	GO annotation
86	PIK3CA	Protein serine/threonine kinase activity and signaling pathway
74	TP53	Transcription factor activity
72	MAPK1	Protein serine/threonine kinase activity and signaling pathway
68	HSPA8	Chaperone and protein folding
64	ACACA	Acetyl-CoA carboxylase activity
63	CAD	Aspartate carbamoyltransferase activity
62	ACACB	Acetyl-CoA carboxylase activity
61	POTEF	Retina homeostasis
60	LRRK2	Protein serine/threonine kinase activity
58	IL8	Neutrophil chemotaxis
54	POTEE	Retina homeostasis
53	POTEI	Retina homeostasis
53	POTEJ	Retina homeostasis
53	PHLPP1	Protein dephosphorylation and signaling pathway
53	RIPK4	Protein serine/threonine kinase activity
53	PHLPP2	Protein dephosphorylation
52	GART	Purine nucleobase biosynthetic process
50	MYC	Transcription factor activity
48	GMPS	Purine nucleobase biosynthetic process
47	OAS2	Purine nucleobase biosynthetic process

Abbreviations: GO, Gene Oncology; PPi, protein–protein interaction.

Table S4

Transcript factors of common key genes screened by indirect PPi

TP53		IL-8		PIK3CA	ELANE	CCR2		CXCR2
AHR	MEIS1	ALX3	MYB	AP2GAMMA	SPIB	ALX3	TEF1	AHRHIF	PMX1
AML1	MTF1	AP1	NF1C	BCL6	PARP	AML1	THAP1	AML1	POU4F3
AP1	MYB	BCL6	NFE2	CBF	NR2C2	AP1	TORC2	AP2REP	PRRX2
AP2GAMMA	MYOGENIN	CDX1	NKX61	CDX		AP2GAMMA	USF	BARHL1	RORA1
BEN	MZF1	CDX2	NMYC	CDX2		BRCA	USF2	BARHL2	RPC155
BRCA	NEUROD	CDXA	NURR1	CDXA	HSPA1L	CEBP	ZBTB2	BEN	RXRA
CACD	NFAT1	CETS1	PARP	CHCH	KID3	CEBPB	ZEB1	BRCA	SATB1
CDX1	NFAT3	CFOS	PEA3	CPBP	EGR1	CEBPD	ZNF333	BRN1	SMAD2
CDX2	NFAT4	CMYB	PLZFB	EGR3	GABPA	CMYB		CDX1	SOX10
CEBP	NKX25	CPBP	PMX1	ETF		CPBP		CEBPA	SOX17
CEBPA	NKX2B	CREB	PRRX2	ETS1		CPEB1		CEBPB	SOX18
CETS1	NMYC	ELF1	RAX	FOS	HSPA1A	DLX3		CHCH	SOX4
CHCH	NR1B2	ELK1	RORBETA	FOSL1	CHCH	E2A		CP2	SOX5
CJUN	OSX	EMX1	SATB1	FOXM1	CPBP	E2F1		CPBP	SPI1
CMAF	PARP	ETS2	SF1	FOXO1A	E2F3	EGR3		DMBX1	SPIB
CMYB	PAX4	ETV7	SOX17	GABPA	EGR1	EMX1		E2A	SREBP2
CPBP	PEA3	EVX1	SOX4	GATA3	KLF7	EVX1		E2F1	SRY
CPEB1	PITX2	EVX2	SPI1	GKLF	KLF8	EVX2		E2F3	STAT3
DR4	PLAGL2	FLI1	SPIB	HFH8	LKLF	FOXK1		ELF1	TAL1
E2F3	PRRX2	FOSL2	SREBP1	HNF3A	MOVOB	FOXP3		ETS	TCF1
EBF1	PU1	FOXC1	TAF1	HOXA13	PLAGL2	GATA1		ETV7	TCF11
EGR1	PUR1	FOXD2	TATA	ISL2	SP2	GKLF		FOXC1	TCFE2A
EHF	RELA	FOXD3	TBP	KID3	SP3	GR		FOXO1A	TFAP2C
ELF1	RORBETA	FOXG1	TEF1	KLF	SP6	HMGIY		GATA2	TORC2
ELF5	SALL2	FOXI1	TTF1	LKLF	VMYB	HOXA10		GKLF	USF2
ELK1	SATB1	FOXK1	ZNF333	MAZR	ZBTB2	HOXA13		GMEB2	VDRRXRALPHA
ETS2	SMAD	FOXL1		MITF		HOXA2		GR	ZBTB2
ETV7	SMAD2	FOXM1		NFAT1		HOXD3		HOXA13	ZFP770
FLI1	SMAD5	FOX01		NFAT4		HSF2		HSF4	ZNF333
FOXA1	SOX10	FOXO1A		NFE4		IK		IK	ZNF536
FOXJ3	SOX11	FOXO3		NR2E1		IK2		ING4
FOXL1	SOX17	FOXO4		PLZFB		KID3		IPF1
FOXM1	SOX18	FOXO6		PMX1		KLF		IRF7
FOXO1	SOX30	FOXP3		RXRA		LHX2		ISL2
FOXO1A	SOX4	FRA1		SALL2		LKLF		IKD3
FOXO3A	SOX9	FXR		SATB1		MAX		LBP1
FOXP3	SPI1	GATA1		SMAD5		MEOX2		LMX1
FRA1	SPIB	GATA3		SOX4		MYB		LPOLYA
GABPA	SREBP1	GATA4		SOX5		MYCMAX		MEF2C
GATA1	SRY	GATA5		SPIB		NF1A		MOVOB
GATA2	STAT	GSX1		SREBP2		NFATC2		MTF1
GATA3	STAT3	GSX2		SRY		NMYC		MYC
GATA4	TCF4	HMGIY		TCF		2-Oct		MYOD
GKLF	TEL2	HMX3		TFE		PARP		MYOGENIN
HDAC1	THAP1	HOXA1		WT1		PAX5		MZF1
HMGIY	TORC2	HOXA13		ZAC		PEBP2B		NANOG
HNF3A	USF	HOXA2		ZFP641		PMX1		NEUROD
HNF3G	USF2	HOXB13		ZNF333		PR		NKX32
HOXA13	VMYB	HOXB5		ZNF641		PRRX2		NMYC
HSF1	WT1	HOXC13				RELA		NR1B2
HSF4	YY1	HOXD13				SALL2		NR2C2
IK	ZBTB44	JUNB				SATB1		1-Oct
ING4	ZEB1	KID3				SMAD2		OG2
IRF1	ZFP532	LBX2				SOX10		OTX
IRF7	ZFP536	LMX1A				SOX17		P300
KID3	ZIC1	LRH1				SOX18		PARP
KLF	ZIC3	MEF2D				SPIB		PAX5
KLF17	ZNF333	MEOX2				SREBP2		PEA3
LHX2	ZNF367	MIXL1				TCF1		PEBP2B
LKLF	ZNF515	MSX2				TCF11		PIT1

Note: Key genes are highlighted in bold.

Abbreviation: PPi, protein–protein interaction.