Roundabout signaling pathway involved in the pathogenesis of COPD by integrative bioinformatics analysis.

Purpose: To explore the potential mechanism underpinning the development of chronic obstructive pulmonary disease (COPD) and to investigate the role of the Roundabout signaling pathway in COPD.
Methods: Three microarray datasets (GSE1650, GSE38974 and GSE76925) including 139 cases of severe COPD and 52 cases of normal smokers without carcinoma, were integrated to screen differentially expressed genes (DEGs) using bioinformatics methods. Gene ontology (GO) annotations and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway of the DEGs were performed by a DAVID online tool. Finally, a cigarette smoke (CS)- induced emphysema mice model was established, the lung mRNA expression levels of genes associated with Slit guidance ligand 2 (SLIT2) -Roundabout (ROBO) signaling pathway were detected by reverse transcription quantitative polymerase chain reaction (RT-qPCR), and the protein level of SLIT2 was examined by immunohistochemistry staining.
Results: A total of 315 DEGs were identified in three databases. GO and KEGG pathway analyses suggested that the inflammatory response, extracellular matrix disassembly, immune response, the apoptotic signaling pathway, ubiquitination and the Roundabout signaling pathway all together were involved in the development of COPD. The genes SLIT2 and ROBO2 were decreased in patients with COPD and these decreases were significantly negatively correlated with the disease stages of COPD. Consistently, the mRNA expression levels of SLIT2, ROBO1 and ROBO2, and the protein level of SLIT2 were revealed to be lower in the lungs of CS-induced emphysema mice compared with the air-exposed control mice. In addition, the SLIT2 protein level was negatively associated with alveolar mean linear intercept.
Conclusion: Integrated bioinformatics analysis may provide novel insights into the complicated pathogenesis of COPD, and to the best of our knowledge, this study is the first to provide evidence to suggest that the Roundabout signaling pathway may be involved in the pathogenesis of COPD.

Introduction

As a typical chronic respiratory disease, chronic obstructive pulmonary disease (COPD), with a mean prevalence rate of 13.1% worldwide,1 is characterized by progressive and irreversible airflow restriction and recurrent respiratory symptoms. As many as 3.2 million patients with COPD succumbed to mortality in 2015-, and COPD was the eighth leading cause of global disability.2 It is widely acknowledged that cigarette smoke (CS) is the predominant risk factor of COPD;3 however, only a few smokers finally develop emphysema, the exact molecular and cellular pathogenesis of this complicated process has yet to be fully elucidated. Additionally, effective prevention and treatment methods available for the disease are still insufficient. Therefore, there is an urgent need to further investigate the molecular mechanisms of COPD to identify novel targets for the treatment of the disease.

In recent decades, mRNA microarrays have been widely applied in cases of multiple tumors, as well as COPD.4–6 Furthermore, an ever-increasing number of chip data recordings are available on public databases, and integrating these microarray datasets should enable us to further explore genes associated with diseases, and thereby to predict the molecular targets for precision therapy.7 Recently, a few investigators have conducted mRNA gene expression profiles on COPD,4,5 and identified hundreds of differentially expressed genes (DEGs) that are likely to be involved in the pathogenesis of COPD. However, the DEGs were found to be inconsistent or even contradictory based on the comparison of the different studies, and the reasons underlying this phenomenon are complex. The lung tissues included in the dataset were mostly from COPD patients that accompanied with lung cancer, which may have introduced bias into the analysis. Besides, different individual samples, different methods of analysis, different disease stages, and different chip platforms also contributed to the difficulties in interpreting the data. However, integrating microarray data from several datasets may circumvent the limitations associated with single datasets, and enable the identification of reliable and efficacious molecular markers.

The Roundabout signaling pathway, consisting predominantly of slit guidance ligand (SLIT) and roundabout (ROBO), was initially shown to be involved in neural development by inhibiting the migration of axons.8 Subsequently, another study demonstrated that SLIT was also able to inhibit the migration of immune cells, including T lymphocytes, neutrophils, macrophages, and dendritic cells,9 and therefore SLIT was regarded as a potential anti-inflammatory molecule.9 T lymphocytes, neutrophils, macrophages and dendritic cells are all involved in the development of COPD,10–12 however, the role of the Roundabout signaling pathway in COPD has yet to be elucidated.

Patients with severe COPD usually suffer from obvious respiratory symptoms that result in high rates of mortality, which places a heavy burden on society. On this basis, the present study focused on the differences in the genes and pathways comparing between cases of severe COPD and smokers to explore the novel mechanism underlying COPD. In the present study, the DEGs from three NCBI Gene Expression Omnibus (GEO) databases were screened for COPD and integrated to explore the potential pathogenesis of COPD according to bioinformatics analysis. By analyzing the biological processes of the differences between COPD and smokers, the present study has attempted to unravel the mechanism to explain how only a few smokers develop COPD, and also to investigate the role of the Roundabout signaling pathway in COPD and cigarette smoke (CS)- induced emphysema mice model.

Materials and methods

mRNA microarray chip datasets

Three original microarray datasets [GSE1650,13 GSE389744 and GSE7692514] were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/), which consisted of 139 patients with COPD and 52 healthy smokers.

Searching on the GEO database with the keywords “COPD mRNA”, an emphasis was placed on those datasets which comprised patients with severe COPD and normal smokers without carcinoma. Global initiative for chronic obstructive lung disease (GOLD) stages 3 and 4 were defined as severe COPD. The platform for GSE1650 was GPL96, [HG-U133A] Affymetrix Human Genome U133A array, which contained 18 lung samples from cases with severe COPD, and 3 lung specimens from normal smokers. The platform for GSE38974 was GPL4133, Agilent-014850 Whole Human Genome microarray 4x44K G4112F (feature number version), which comprised 10 lung samples from cases of severe COPD and 9 lung specimens from normal smokers. Finally, the platform for GSE76925 was GPL10558, Illumina HumanHT-12 V4.0 expression bead-chip, which included 111 lung samples from cases of severe COPD and 40 lung specimens from normal smokers. Detailed information of the three datasets is shown in Table 1.

Table 1

Details for GEO COPD data

Reference	Sample	GEO	Platform	Smoker	Severe COPD
Spira A et al (2004)13	Lung tissue	GSE1640	GPL96	3	18
Ezzie ME et al (2012)4	Lung tissue	GSE38974	GPL4133	9	10
Morrow JD et al (2017)14	Lung tissue	GSE76925	GPL10558	40	111

Abbreviations: GEO, gene expression omnibus; COPD, chronic obstructive pulmonary disease.

Screening for DEGs, and subsequent integration of DEGs in three datasets

The .txt files of the DEGs for the three microarray datasets were obtained using GEO2R online analysis (https://www.ncbi.nlm.nih.gov/geo/geo2r/), and the DEGs in the three datasets were ranked by logFC respectively. The significant DEGs in the three datasets were defined as such when the corrected p-value was <0.05, and the absolute folder change (FC) was >1.5. The raw data were corrected for background noise normalized by the Robust Microarray Analysis (RMA) method. Subsequently, DEGs for the three microarray datasets were integrated using the Robust Rank Aggregation (RRA) method (http://www.icesi.edu.co/CRAN/web/packages/RobustRankAggreg/)15 in R version 3.5.3 (http://cran.rproject.org/). The significant integrated DEGs were defined when the corrected p-value was <0.05.

Protein-protein interaction (PPI) network analysis

PPI networks of integrated DEGs were constructed using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) online database (http://stringdb.org/), and the Cytoscape software 3.5.116 and CytoHubba Apps17 were applied to reveal the molecular interactions and to identify the top 20 hub genes associated with COPD.

Gene ontology (GO) and kyoto encyclopedia of genes and genomes (KEGG) pathway analysis for DEGs

GO annotations and KEGG pathway analysis of the integrated DEGs were performed by the database for annotation, visualization, and integrated discovery (DAVID) online tool. The results were downloaded into a .txt file, and subsequently, GOPlot18 and ggplot2 packages in R were used for the visual analysis. The genes associated with the SLIT2- ROBO pathway were further analyzed in COPD patients in the dataset GSE38974.

CS-induced emphysema mice model

Male C57BL/6 mice (6–8 weeks old, 20–25 g body weight) were purchased from the Guangxi Medical University Laboratory Animal Centre (Nanning, China) and were exposed to air (control) or cigarette smoke for 24 weeks (n=6 per group), as previously described.11 Following 24 weeks’ exposure to air or CS, the mice were sacrificed after anesthesia, and the lungs were subsequently removed from the chest cavity, the left lungs were used for RNA analysis, and the right lungs were used for pathological and immunohistochemistry examination. All animal experiments protocols were approved by the Animal Research of Guangxi Medical University Ethics Committee.

Quantifications of mean linear intercept (Lm) and immunohistochemistry staining

The right lungs were fixed in 10% formalin for 24 h before embedding in paraffin; midsagittal sections were stained with hematoxylin and eosin (HE) for histological analysis, we determined enlargement of alveolar spaces by quantifying the Lm as previously described.11 Lung tissue sections were deparaffinized with xylene, hydrated with alcohol, repaired with sodium citrate in high temperature and then incubated with antibody against mouse Slit2 (ab134166; Abcam, Cambridge, UK) overnight at 4 °C. Subsequently, secondary Goat Anti-Rabbit IgG H&L (HRP) antibody (ab205718, Abcam, Cambridge, UK) was stained at room temperature for 20 mins and finally reacted with diaminobenzidine (Solarbio, DA1010, China) for coloration. Fluorescence pathological microscope (Olympus */BX53+DP80, Japan) was used to take photos for all the sections. The mean density of SLIT2 staining was quantified using Image-Pro Plus software 6.0 (Media Cybernetics), and the average score of 2 different viewers was taken for analysis.

RNA extraction and reverse transcription- quantitative polymerase chain reaction (RT-qPCR)

The lung mRNA of air or CS-exposed mice (n=6 per group) were extracted using Eastep® Super Total RNA Extraction kit (Promega Corporation), and then cDNA was prepared using GoScript™ Reverse Transcription mix (Promega Corporation). Quantitative PCR was performed using TB Green® Premix Ex Taq ™ II (Takara Biotechnology Co., Ltd). All experiments were performed by the kit instructions. Data are presented relative to β-actin for target genes based on calculations of 2−ΔΔCt. The mRNA expression results are reported as the quantity of mRNA relative to control. Primers were synthesized by Sangon Biotech Company (Shanghai Co., Ltd.). Primers used are shown in Table 2.

Table 2

Details for primers of mice

Gene name	Primers (5ʹ—3ʹ)
SLIT2	FORWORD: AGCCTTGTCACACTTAGCGATTGGREVERSE: ACAGCGTGCAATTCCAGGTTCC
ROBO2	FORWORD: TTGGAGCAAGTTCACGGGAGREVERSE: TAAGCCGCTCTGTTAGTCGG
ROBO1	FORWARD: CCTTTTTCACTGACGCTGATTTREVERSE: GATCACAACAGCACTCATGTTG
CDC42	FORWARD: CAGACTACGACCGCTAAGTTATREVERSE: CAGCAGTCTCTGGAGTAATAGG
RAC2	FORWARD: CGAGAAGCTGAAGGAGAAGAAGREVERSE: GAGTGCAGAACATTCCAAGTAC
B-ACTIN	FORWARD: GTGCTATGTTGCTCTAGACTTCGREVERSE: ATGCCACAGGATTCCATACC

Abbreviations: SLIT2, slit guidance ligand 2; ROBO, roundabout guidance receptor; RAC2, Rac family small GTPase 2; CDC42, cell division cycle 42; B-ACTIN, actin beta.

Statistical analysis

The GraphPad Prism 5, and R, version 3.5.3 software were used for statistical analysis. Data are presented as mean ± SD and compared between groups using unpaired Student’s t-tests. Correlation analysis is calculated using Pearson or Spearman rank correlation test. p<0.05 was considered to indicate a statistically significant value.

Results

Identification of DEGs in three datasets

After screening with the threshold of a corrected p-value <0.05 and absolute folder change >1.5, 52 DEGs (37 upregulated and 15 downregulated) were identified in the GSE1650 dataset, 4957 DEGs (1989 upregulated and 2968 downregulated) were identified in GSE38974, and 1649 DEGs (367 upregulated and 1282 downregulated) were identified in GSE76925. Figure 1 shows the heatmaps of the top 200 DEGs in the three datasets.

Figure 1

Hierarchical clustering heatmap of DEGs screened on the basis of |fold change| >1.5 and a corrected P-value <0.05.

Notes: (A) GSE1650 data, (B) GSE38974 data, (C) GSE76925 data. Red indicates that the expression of genes is relatively upregulated, green indicates that the expression of genes is relatively downregulated, and black indicates no significant changes in gene expression. Orange indicates samples from severe COPD patients, and blue indicates samples from normal smokers.

Abbreviations: DEGs, differentially expressed genes; COPD, chronic obstructive pulmonary disease.

Identification of integrated DEGs in severe COPD using the RRA method

The DEGs for the three microarray datasets were sorted according to their logFC values and subsequently were further integrated using the RRA method (adjusted p-value <0.05). At the end of this analysis, 315 DEGs, including 119 upregulated genes and 195 downregulated genes, were identified. Table 3 presents the integrated DEGs in detail. The top 20 up- and down-regulated hub genes were drawn into a heatmap using R–heatmap software, as shown in Figure 2.

Table 3

Integrated DEGs in COPD by RRA

DEGs	Gene names
Upregulated	APOL6,ADAMTS4,FOXE1,PTX3,TEP1,CA3,CACNA1E,HNRNPA0,METTL7A,NKD2,GNL3L,FGG,ARNTL2,IL6,SPP1,ARHGAP26,TUBB3,LMNB1,HTR2B,SYNPO,CHIT1,MEDAG,CCL19,ZNF652,MACC1,PHLDA2,CXCL8,FLJ35934,MMP9,SLCO4A1,HRASLS5,MT1M,S100A9,DPYS,PSORS1C2,POU2AF1,SFN,SMG1P5,IL1R2,GFPT2,TNFRSF17,CHMP4B,FBRSL1,CXCL13,CSF3,PLA2G2A,TMEM201,TM4SF19,TAF5L,SLC39A7,TRIM35,IRX3,ANKDD1A,COL10A1,KIAA2012,DLST,PRKAR2A,CYCS,NEUROG1,COMP,HSD17B1,MMP1,CH25H,ITPKC,LOC649294,EFS,SLC7A5,PRRG4,HPR,MT1JP,CCDC71,HMOX1,HS3ST2,SP5,TNFRSF1A,MT1L,SNAI1,IGHM,SOX8,CAND1,MYCN,PVR,SLC20A2,ACBD4,FCRLA,RDH10,ZSCAN10,FDCSP,MT1G,ITPRIPL2,DAPL1,LMF1,RNF125,RND1,KRT16,MT1H,LIF,TSPAN10,SPPL2B,NRSN2,SLC43A2,CLDN1,PDAP1,CD19,CLECL1,F2RL3,ATM,C5AR1,KRT7,CRHR2,THBS1,CBLN4,HIF1A,UBQLN1,BHLHA15,CLMP,SHROOM1,LAD1,SERPIND1
Downregulated	CD86,LRRC2,KLF5,CABP7,PFKFB3,SUMO1,PPP1R3C,GDF10,ID4,ESM1,PPP6R1,DST,AGR3,PLN,AGBL1,SLC12A2,IRX5,MAF,TDRD10,ADH1A,SLIT2,PTCH1,FRY,HIAT1,ITLN2,NOV,ADI1,CYP26B1,PCDH18,DNASE1L3,DNTTIP2,OGN,RAPH1,P2RY14,B3GNT5,SOSTDC1,NR2F2,INMT,ADH1C,KIAA1033,DDR2,HP1BP3,PSMC6,LOC100128164,OGT,TTRAP,DACT3,SEMA6D,SAMSN1,ACTG2,L1TD1,PPIB,ROBO2,ANKRA2,SFTPA1,PAPSS1,IL33,HMCN1,PRDM6,MYH11,HSPC324,PLEKHH2,SLC51B,CNN1,LRIG1,GGH,HOXA4,MYZAP,VPS18,LIMS2,AQP4,HMGB1,TMEM100,FAM162B,SFRS3,ZNF177,TBC1D31,KRT4,SOX4,CD164,ATP1A2,COLGALT2,C7orf11,IQCD,SDC1,TANK,SFRS11,ACTA2,EFCC1,LINC01089,SMAD5,PI16,PABPC4L,RUNX1T1,C1QTNF7,HERPUD2,ASPN,PRRT2,MYOC,ITM2A,NOTCH3,NEK3,FHL1,SPAST,EDNRB,SH3GL3,AGER,TAF1D,LTBP2,RAB23,MAOB,POSTN,CYP3A4,SMAD7,SBSPON,CD302,SBF2,CX3CR1,ARL6IP1,ANKMY2,SLCO1A2,LSM5,ZBED8,ANGPT1,ROCK1,IFIT1,BCHE,RERG,PAX4,TSNAX,MEGF6,TCEAL4,MED21,GRB14,PPP1R15B,SLITRK6,ZNF575,CARD16,FAM117B,FCER1A,FABP5,TUBE1,SLC1A1,PPP3CA,SEMA5A,ENPP2,AZIN1,SCN4B,NPIPB3,ADGRG6,U2AF1,TNXB,SGCE,RABL2A,FRAS1,PNPLA8,HSPB3,MS4A2,ZNF23,ECM2,DKK2,AKAP12,GPR20,CDH6,ALG13,SLC35A1,PPP3R1,MST1,FENDRR,PDZD2,PDCD10,HOXA2,CCDC91,LUC7L3,CROCCP3,MRPL47,LDB2,MAT2B,CIRBP,EIF4E,NOTCH2NL,PPM1K,FOLH1,ST6GALNAC5,SPOPL,FIBIN,NPY1R,SPAG1,GNG13,LOC100507311,KRT18P17,CDH13,MYRIP,PCDH20,TSN

Abbreviations: DEGs, differentially expressed genes; COPD, chronic obstructive pulmonary disease; RRA, robust rank aggregation.

Figure 2

LogFC heatmap of top 20 up- and down-regulated DEGs of three microarrays screening by RRA method.

Notes: The abscissa is the GEO ID, and the ordinate is the gene name. Red represents logFC >0, green represents logFC <0, and the values in the box represent the logFC values.

Abbreviations: FC, fold change; GEO, Gene Expression Omnibus; DEGs, differentially expressed genes.

PPI network

The PPI network for integrated DEGs is shown in Figure 3. The 20 hub genes that extensively interacted with other genes were IL6, CXCL8, MMP9, SNAI1, CYCS, HIF1A, SPP1, THBS1, MMP1, POSTN, CXCL13, CD86, CD19, SDC1, HMOX1, ATM, CX3CR1, GNG13, ANGPT1 and TNFRSF1A (Figure 4).

Figure 3

PPI network of DEGs.

Notes: Circles represent genes, lines represent the interaction of proteins between genes, and the results within the circle represent the structure of proteins.

Abbreviations: PPI, protein–protein interaction; DEGs, differentially expressed genes.

Figure 4

The 20 top hub genes in PPI network.

Notes: The diamond icon represents the downregulated hub genes, the circle icon represents the upregulated hub genes, the size of the icon represents the strength of nodes that interact with other genes, lines represent the interaction of proteins between genes.

Abbreviation: PPI, protein–protein interaction.

GO analysis of integrated DEGs

Three functional groups in the GO analysis of integrated DEGs were sought after. Namely, cell composition (CC), biological processes (BP), and molecular function (MF). The significant GO enrichment results are shown in Tables 4 and 5, and Figures 5 and 6. With the BP functional group, the up-regulated DEGs were concentrated in the categories of negative regulation of growth, inflammatory response, extracellular matrix (ECM) disassembly, immune response, intrinsic apoptotic signaling pathway in response to DNA damage and positive regulation of the extracellular-signal- regulated (ERK) 1 and ERK2 cascade. By contrast, the downregulated DEGs were concentrated in the categories of negative regulation of protein ubiquitination, Roundabout signaling pathway, cell adhesion, cell fate determination, and negative regulation of neuron differentiation. With the MF functional group, the upregulated DEGs were concentrated in the categories of structural molecule activity, C-C motif chemokine receptor 10 (CCR10) chemokine receptor binding, and chemokine activity. In addition, the downregulated DEGs were mainly associated with the categories of heparin- binding, protein kinase A binding, collagen binding, calcium ion binding, protein serine/threonine phosphatase activity, and integrin binding. Finally, with the CC functional group, the upregulated DEGs were concentrated in the categories of an extracellular region, extracellular exosome, and extracellular space, the downregulated DEGs were mainly associated with the plasma membrane, proteinaceous ECM and ECM. Considering all these results together, the majority of the DEGs were found to be enriched in the categories of the inflammatory response, immune response, apoptotic signaling pathway, ECM disassembly, ubiquitination, and the Roundabout signaling pathway.

Table 4

GO analysis of upregulated genes associated with COPD

Term	Description	Count	p-Value
Upregulated
GO:0005576	Extracellular region	27	1.80E-06
GO:0045926	Negative regulation of growth	5	4.64E-06
GO:0005615	Extracellular space	23	1.12E-05
GO:0006954	Inflammatory response	12	1.90E-05
GO:0071294	Cellular response to zinc ion	4	1.96E-04
GO:0006955	Immune response	11	2.41E-04
GO:0008630	Intrinsic apoptotic signaling pathway in response to DNA damage	4	0.002906525
GO:0070374	Positive regulation of ERK1 and ERK2 cascade	6	0.004286116
GO:0071276	Cellular response to cadmium ion	3	0.004682164
GO:0061436	Establishment of skin barrier	3	0.005246663
GO:0045766	Positive regulation of angiogenesis	5	0.005299858
GO:0090023	Positive regulation of neutrophil chemotaxis	3	0.007797363
GO:0071347	Cellular response to interleukin-1	4	0.009225227
GO:0022617	Extracellular matrix disassembly	4	0.011102727
GO:0051897	Positive regulation of protein kinase B signaling	4	0.014534642
GO:0045765	Regulation of angiogenesis	3	0.015149905
GO:0005198	Structural molecule activity	6	0.015754473
GO:0009306	Protein secretion	3	0.016096591
GO:0031735	CCR10 chemokine receptor binding	2	0.017667432
GO:0048469	Cell maturation	3	0.020127137

Abbreviations: GO, gene ontology; COPD, chronic obstructive pulmonary disease.

Table 5

GO analysis of downregulated genes associated with COPD

Term	Description	Count	p-Value
Downregulated
GO:0008201	Heparin binding	9	1.86E-04
GO:0005578	Proteinaceous extracellular matrix	11	2.98E-04
GO:0001657	Ureteric bud development	5	5.47E-04
GO:0005509	Calcium ion binding	18	7.13E-04
GO:0031012	Extracellular matrix	10	0.002482215
GO:0005518	Collagen binding	5	0.002865969
GO:0007162	Negative regulation of cell adhesion	4	0.005955399
GO:0086004	Regulation of cardiac muscle cell contraction	3	0.006148552
GO:0035385	Roundabout signaling pathway	3	0.006148552
GO:0007155	Cell adhesion	12	0.006554914
GO:0031397	Negative regulation of protein ubiquitination	4	0.007936104
GO:0001709	Cell fate determination	3	0.013704383
GO:0005886	Plasma membrane	54	0.014038038
GO:0051018	Protein kinase A binding	3	0.014727798
GO:0045665	Negative regulation of neuron differentiation	4	0.018512605
GO:0004722	Protein serine/threonine phosphatase activity	4	0.019415924
GO:0005178	Integrin binding	5	0.020021107
GO:0009611	Response to wounding	4	0.025216936
GO:0002040	Sprouting angiogenesis	3	0.025673813
GO:0008217	Regulation of blood pressure	4	0.027341577

Abbreviations: GO, gene ontolog; COPD, chronic obstructive pulmonary disease.

Figure 5

GO enrichment analysis of DEGs in severe COPD.

Notes: (A) GO analysis divided DEGs into three functional groups: molecular function, biological processes and cell composition. (B) GO enrichment significance items of DEGs in different functional groups.

Abbreviations: DEGs, differentially expressed genes; GO, gene ontology; COPD, chronic obstructive pulmonary disease.

Figure 6

Distribution of DEGs in severe COPD for different GO-enriched functions.

Abbreviations: DEGs, differentially expressed genes; GO, gene ontology; COPD, chronic obstructive pulmonary disease.

KEGG pathway analysis of integrated DEGs

As shown in Figure 7 and Table 6, the KEGG pathway analysis of the DEGs revealed that they were concentrated in the categories of mineral absorption, cytokine-cytokine receptor interaction, axon guidance, transcriptional misregulation in cancer, retinol metabolism and the hypoxia-inducible factor-1 (HIF-1) signaling pathway.

Figure 7

KEGG pathway enrichment analysis for DEGs.

Notes: The x-axis shows the -log10 (p-value) of each term and y-axis shows the KEGG pathway terms.

Abbreviations: DEGs, differentially expressed genes; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Table 6

KEGG pathway analysis of DEGs associated with COPD

ID	Pathway	Count	p-Value	Genes
hsa05144	Malaria	6	0.002186	CSF3, IL6, SDC1, COMP, CXCL8, THBS1
hsa04978	Mineral absorption	5	0.009218	MT1M, HMOX1, ATP1A2, MT1H, MT1G
hsa04360	Axon guidance	8	0.010029	SEMA5A, RND1, ROCK1, SEMA6D, PPP3R1, ROBO2, PPP3CA, SLIT2
hsa05202	Transcriptional misregulation in cancer	9	0.013515	MAF, IL1R2, IL6, CD86, MMP9, RUNX1T1, CXCL8, ATM, MYCN
hsa00830	Retinol metabolism	5	0.032375	CYP3A4, RDH10, CYP26B1, ADH1C, ADH1A
hsa04066	HIF-1 signaling pathway	6	0.034808	IL6, EIF4E, HIF1A, PFKFB3, HMOX1, ANGPT1
hsa04060	Cytokine-cytokine receptor interaction	10	0.040111	LIF, CSF3, IL1R2, TNFRSF1A, IL6, CXCL13, CX3CR1, CXCL8, CCL19, TNFRSF17
hsa05219	Bladder cancer	4	0.041745	MMP9, CXCL8, THBS1, MMP1

Abbreviations: DEGs, differentially expressed genes; KEGG, Kyoto Encyclopedia of Genes and Genomes; COPD, chronic obstructive pulmonary disease.

The SLIT2-ROBO pathways in COPD

In the GO enrichment analysis, it was revealed that the genes SLIT2 and ROBO2 were enriched in the Roundabout signaling pathway, which has not previously been associated with the pathogenesis of COPD, and furthermore, both genes were downregulated in severe COPD; therefore, the present study sort to investigate the role of the SLIT2-ROBO pathway in all disease stages of the COPD patients. The dataset GSE38974 included not only smokers and patients with severe COPD but also patients with mild and moderate COPD; therefore, the expression matrix of GSE38974, including a total of 31 samples, was downloaded. The mRNA expression levels of SLIT2 and ROBO2 were found to be lower in all stages of COPD compared with smokers (Figure 8A and B), and both were negatively correlated with the GOLD stages in COPD (Figure 8D and E). However, no statistical differences in the expression levels of ROBO1 were identified between COPD and smokers (Figure 8C and F). The expression levels of SLIT2 and ROBO2 were independent of age and smoking status (Figure 8G–J). In addition, the expression of SLIT2 was positively correlated with ROBO2, rather than with ROBO1 (Figure 8K and L). The mRNA expression levels of cell division cycle 42 (Cdc42) and Rac family small GTPase 2 (Rac2) were increased in COPD, and both were significantly positively correlated with the GOLD stages in COPD (Figure 9A–D). The expression of SLIT2 was negatively correlated with Cdc42 and Rac2, whereas the expression of ROBO2 was negatively correlated with Rac2, rather than with Cdc42 (Figure 9E–H).

Figure 8

The SLIT2 and ROBO expression in GSE38974.

Notes: The mRNA expression of SLIT2 (A) ROBO2 (B) and ROBO1 (C) in GSE38974 microarray dataset. (D) Scatter plot of correlation analysis for SLIT2and GOLD stage. (E) Scatter plot of correlation analysis for ROBO2 and GOLD stage. (F) Scatter plot of correlation analysis for ROBO1 and GOLD stage. The GOLD 0 represents smokers without COPD, GOLD 1–4 represents COPD patients with GOLD stages 1–4. (G) Scatter plot of correlation analysis for SLIT2 and age. (H) Scatter plot of correlation analysis for ROBO2 and age. (I) Scatter plot of correlation analysis for SLIT2 and the pack year of smoking. (J) Scatter plot of correlation analysis for ROBO2 and the pack year of smoking. (K) Scatter plot of correlation analysis for SLIT2 and ROBO2. (L) Scatter plot of correlation analysis for SLIT2 and ROBO1.

Abbreviations: GOLD, Global initiative for chronic obstructive lung disease; SLIT2, slit guidance ligand 2; ROBO, roundabout guidance receptor; COPD, chronic obstructive pulmonary disease.

Figure 9

The Cdc42 and Rac2 expression in GSE38974.

Notes: The mRNA expression of Cdc42 (A) and Rac2 (B) in GSE38974 microarray dataset. (C) Scatter plot of correlation analysis for Cdc42 and GOLD stage. (D) Scatter plot of correlation analysis for Rac2 and GOLD stage. The GOLD 0 represents smokers without COPD, GOLD 1–4 represents COPD patients with GOLD stages 1–4. (E) Scatter plot of correlation analysis for SLIT2 and Cdc42. (F) Scatter plot of correlation analysis for ROBO2 and Cdc42 (G) Scatter plot of correlation analysis for SLIT2 and Rac2. (H) Scatter plot of correlation analysis for ROBO2 and Rac2.

Abbreviations: GOLD, Global initiative for chronic obstructive lung disease; SLIT2, slit guidance ligand 2; ROBO, roundabout guidance receptor; COPD, chronic obstructive pulmonary disease; Rac2, Rac family small GTPase 2; Cdc42, cell division cycle 42.

The SLIT2-ROBO signal pathway is downregulated in CS-induced emphysema mice model

As it proved to be difficult for us to obtain the lungs from COPD patients, CS-induced emphysema mice model was alternatively established to investigate whether the SLIT2–ROBO signal pathway is involved in the pathogenesis of emphysema. In accordance with our previous studies,11 the Lm of CS-exposed mice was significantly higher in than that of control mice (Figure 10A and B), which indicated that the emphysema model had been success established. As shown in Figure11A and B, in accordance with the microarray results, compared with the control mice exposed to only air, the mRNA expression levels of SLIT2 and ROBO2 both were lower in the lungs of CS-exposed mice (0.21 folder change, p=0.000 and 0.03 folder change, p=0.001 respectively). Unexpectedly, the mRNA expression levels of ROBO1 was also lower in CS-induced emphysema mice (0.35 folder change, p=0.004, Figure 11C), however, there were no statistical differences in the mRNA expression levels of Cdc42 and Rac2(Figure 11D and E). Furthermore, immunohistochemistry staining of SLIT2 was recruited to delineate the role of SLIT2 in emphysema; the SLIT2 protein expression was consistently decreased in CS-induced emphysema mice and the SLIT2 was mainly expressed in bronchial epithelial cells (Figure10C and D). In addition, the mean density level of SLIT2 was negatively correlated with Lm (r= −0.645, p=0.024, Figure 10E).

Figure 10

The mean linear intercept (Lm) and immunohistochemistry staining of SLIT2 in lung of AIR-exposed and CS-exposed mice.

Notes: Representative photomicrographs of hematoxylin and eosin stained lung tissue of 24 weeks AIR-exposed and CS-exposed mice, magnification ×20 (A) The Lm of alveolar in 24 weeks AIR-exposed and CS-exposed mice (B) The photomicrographs for immunohistochemistry staining of SLIT2 in lung of AIR-exposed and CS-exposed mice, magnification ×40 (C) The mean density value of SLIT2 staining in the lung of AIR-exposed and CS-exposed mice (D) Scatter plot of correlation analysis for mean density of SLIT2 staining and Lm in lung of mice (E) AIR-24W represents mice exposed to air for 24 weeks, and CS-24w represents mice exposed to CS for 24 weeks. **p<0.01, ***p<0.001.

Abbreviations: Lm, the mean linear intercept; CS, cigarette smoke; SLIT2, slit guidance ligand 2.

Figure 11

The SLIT2-ROBO signal pathway is downregulated in CS-induced emphysema mice model.

Notes: The relative mRNA expression of SLIT2 (A) ROBO2 (B) ROBO1 (C) Rac2 (D) and Cdc42 (E) in lungs of mice exposed to air or CS for 24 weeks. AIR-24W represents mice exposed to air for 24 weeks, and CS-24w represents mice exposed to CS for 24 weeks. ** p<0.01, ***p<0.001, ns: not significant.

Abbreviations: SLIT2, slit guidance ligand 2; ROBO, roundabout guidance receptor; CS, cigarette smoke; Rac2, Rac family small GTPase 2; Cdc42, cell division cycle 42.

Discussion

In the present study, 314 DEGs were identified, among which 119 genes were upregulated, and 195 genes were downregulated (Table 2). Subsequently, a PPI network of the DEGs was constructed, and the top 20 hub genes were identified (Figure 4). The top 20 hub genes are mainly associated with the downstream inflammatory response and alveolar destruction that has been widely accepted, however, these hub genes, in themselves, are insufficient to explain the mechanism by which the immune reaction is initiated in patients with COPD. Chemokine (C-X-C motif) ligand 8 (CXCL8) secreted by lung macrophages is able to attract an accumulation of polymorphonuclear neutrophils to the lungs, thereby stimulating the airway epithelium to contract and increase its permeability to inflammatory cells.19 Interleukin-6 (IL-6), CXCL13 and Chemokine (C-C motif) ligand 19 (CCL19), B-cell recruitment and immunoglobulin transcription genes, along with CD19, revealed that the B cell signaling pathway has an important role in the development of COPD. CXCL13 is increased in the pulmonary B-cell regions of lymphoid follicles of patients with COPD and CS-exposed mice, and therapeutic targeting of CXCL13 was shown to absolutely reverse lymphoid follicle formation induced by CS in mice.20 CS is able to induce pulmonary antinuclear antibodies and anti-tissue antibodies, resulting in autoimmune processes and alveolar damage.21 Tumor necrosis factor receptor superfamily member 1A (TNFRSF1A) is associated with symptom burden in COPD,22 which is connected with the apoptotic signals in airway epithelium cells.23 CX3CR1 is essential for the generation of TNF-α and IL-6 by mononuclear phagocytes exposed to CS, and CS-induced emphysema was attenuated in the CX3CR1 knockout mice.24 Increasing expression levels of matrix metalloproteinase (MMP)-1 and MMP-9 were shown to lead to alveolar destruction, and this was demonstrated to be positively correlated with COPD severity.25 HIF1A, C5AR1, heme oxygenase 1 (HMOX1), CXCL8 and thrombospondin 1 (THBS1) are associated with the oxidative stress reaction and are enriched in the positive regulation of angiogenesis (data not shown). In our previous studies, the expression level of CD86 in dendritic cells was elevated both in CS-induced emphysema mice and COPD patients;10 however, in the present study, compared with healthy smokers, the expression of CD86 was reduced in patients with severe COPD, although the reason for this is unclear. Possibly, the patients with severe COPD acquired an impaired immune function following repeated bouts of respiratory infection. As the disease progresses, the ability to dispose of damage signals, including bacteria and viruses, may also be impaired and exhausted. Secreted phosphoprotein 1 (SPP1), which encodes osteopontin, is a pleiotropic molecule that has been previously shown to be linked with Th17 responses in COPD.26 Accumulating evidence has supported the proposal that osteopontin was significantly higher in lung tissue and bronchoalveolar lavage fluid of patients with COPD, and furthermore, this was significantly correlated with the levels of neutrophils in sputum and the extent of emphysema.27 A recently published study identified that osteopontin could competitively bind to several antibacterial proteins in the airways, increasing the susceptibility of COPD patients to infections.28

Subsequently, the GO annotations of 314 DEGs were performed on David database. The enriched functions of these genes were, in part, consistent with the pathogenesis of COPD demonstrated in previous studies.11,29,30 For example, an imbalanced activity of proteases and antiproteases identified in lung parenchyma provide one of the important etiological processes in emphysema,29 and the present study showed that the genes enriched in ECM disassembly were upregulated (Table 4), whereas the genes encoding ECM proteins were downregulated (Table 5). In addition, genes encoding collagen- binding proteins were also downregulated (Table 5). In recent years, an increasing number of studies have suggested that the initiation and imbalance of acquired immunity contribute to the amplification and persistence of chronic inflammation in COPD,11,12 and it has been proposed that COPD may also be an autoimmune disease. Similarly, the results of the present study have also demonstrated that genes associated with the immune response and inflammation were significantly differentially expressed in COPD (Table 4). The apoptosis of both endothelial and epithelial cells has been shown to be involved in the development of emphysema,31 and the apoptosis-associated genes were upregulated in the present study (Table 4). Genes associated with negative regulation of protein ubiquitination were decreased in our analysis. The ubiquitin-proteasome system regulates the degradation and function of proteins. CS can induce cytotoxicity via the ubiquitin-proteasome system,30 in addition, airway epithelial cells and macrophages exposed to CS extract may result in more histone deacetylase 2 (HDAC2) degradation via the ubiquitin-proteasome system,32 which may promote the transcription of inflammatory genes. The ubiquitination of contractile proteins was also shown to be involved in the formation of muscle atrophy in patients with COPD.33 In conjunction with the results reported above, this corroborates the credibility of our method of screening for DEGs.

In the present study, the genes Osteoglycin (OGN), ROBO2 and SLIT2 involved in the Roundabout signaling pathway were downregulated in patients with severe COPD (Table 5), a finding which was unanticipated, and which had not been identified in terms of the etiology of COPD previously. However, the role of ROBO and SLIT2 are poorly understood in terms of COPD. In the present study, it was observed that the expression levels of SLIT2 and ROBO2 were decreased in all stages of COPD, and this was significantly negatively correlated with disease progression (Figure 8A–D). As the SLIT and ROBO genes are associated with tissue morphogenesis,34 the present study investigated whether their expression levels were associated with age, and the findings revealed that the expression changes were independent of age and smoking status (Figure 8G–H). In accordance with the microarray results, the protein level of SLIT2, the mRNA expression levels of SLIT2 and ROBO2, were all decreased in the CS-induced emphysema mice compared to the control mice, and the SLIT2 protein level was negatively associated with the degree of emphysema (Figures 10 and 11). In addition, the mRNA expression level of ROBO1 was also decreased in CS-induced emphysema mice. Therefore, it was possible to confirm that SLIT2-ROBO signaling pathway is involved in the pathogenesis of emphysema induced by CS in the mice model.

The ROBO proteins comprise a group of transmembrane receptors belonging to the immunoglobulin superfamily; SLIT proteins, acting as ligands of ROBO, are secreted glycoproteins.35 Four ROBO proteins (ROBO1, ROBO2, ROBO3, ROBO4) have been identified in humans.36 SLIT proteins are endogenously cleaved into a variety of N- and C-terminal fragments, which can inhibit chemoattraction of many immune cells, including T lymphocytes, dendritic cells, macrophages, and neutrophils.9,37,38 Our previous study identified that CS might promote the formation of neutrophil extracellular traps (NETs), which can activate dendritic cells and initiate acquired immunity.12 Therefore, it was possible to hypothesize that the downregulation of the SLIT2- ROBO signaling pathway may result in the attraction of neutrophils and T lymphocytes, which ultimately destroy lung tissue, and forming emphysema; however, this hypothesizes still requires further investigation to verify.

The SLIT2-ROBO pathway is able to inhibit the migration of neutrophils by suppressing the two small GTPases, Cdc42 and Rac2.9,39 Given this knowledge, the role of Cdc42 and Rac2 in patients with COPD was further explored. Consistently, the expression levels of Cdc42 and Rac2 were significantly increased in COPD, and these increases were positively correlated with disease progression (Figure 9A–D). Furthermore, the expression of SLIT2 was negatively correlated with Cdc42 and Rac2 (Figure 9F and H). However, there were no statistical differences in the mRNA expressions of Cdc42 and Rac2 in the lung of the CS-induced emphysema mice compared to control mice. The reason for this is complicated, as the majority cells in lung tissue are able to express Cdc42 and Rac2, it would be better to isolate the neutrophils and then detect the expressions of Cdc42 and Rac2.

The Roundabout signaling pathway was initially found to be involved in neural development via inhibiting the migration of axons.8 In the KEGG analysis of the DEGs, the downregulated genes were also shown to be enriched in the axon guidance (Table 6 and Figure 7). Branchfield et al40 demonstrated that ROBO inactivation in the lung led to an inability of pulmonary neuroendocrine cells to aggregate into sensory organoids, and induced an increased in neuropeptide production which could attract immune cells, and, as a consequence, the ECM was remodeled, and alveoli numbers were reduced. In the CS-induced emphysema mice, the ROBO1 and ROBO2 mRNA expressions both were decreased; however, whether the pulmonary neuroendocrine cells aggregation is also decreased owing to this phenomenon remains incompletely understood, and further experiments are required to further clarify this matter.

Conclusion

The present study has demonstrated that integrated bioinformatics analysis may provide novel insights into elucidating the complicated pathogenetic processes of COPD. The B cell signaling pathway, osteopontin differentiation, inflammatory response, ECM disassembly, immune response, apoptotic signaling pathway, ubiquitination and the Roundabout signaling pathway all together play an important role in the development of COPD. The SLIT2-ROBO signaling pathway was shown to be downregulated in both patients with COPD and CS-induced emphysema mice, and these changes were significantly negatively correlated with disease progression. To our best known, this study is the first to suggest that the Roundabout signaling pathway may be involved in the pathogenesis of COPD.