Identification of Key Genes and Pathways in Gefitinib-Resistant Lung Adenocarcinoma using Bioinformatics Analysis.

Gefitinib resistance is a serious threat in the treatment of patients with non-small cell lung cancer (NSCLC). Elucidating the underlying mechanisms and developing effective therapies to overcome gefitinib resistance is urgently needed. The differentially expressed genes (DEGs) were screened from the gene expression profile GSE122005 between gefitinib-sensitive and resistant samples. GO and KEGG analyses were performed with DAVID. The protein-protein interaction (PPI) network was established to visualize DEGs and screen hub genes. The functional roles of CCL20 in lung adenocarcinoma (LUAD) were examined using gene set enrichment analysis (GSEA). Functional analysis revealed that the DEGs were mainly concentrated in inflammatory, cell chemotaxis, and PI3K signal regulation. Ten hub genes were identified based on the PPI network. The survival analysis of the hub genes showed that CCL20 had a significant effect on the prognosis of LUAD patients. GSEA analysis showed that CCL20 high expression group was mainly enriched in cytokine-related signaling pathways. In conclusion, our analysis suggests that changes in inflammation and cytokine-related signaling pathways are closely related to gefitinib resistance in patients with lung cancer. The CCL20 gene may promote the formation of gefitinib resistance, which may serve as a new biomarker for predicting gefitinib resistance in patients with lung cancer.

Introduction

Based on recent statistics, lung cancer is not only the most commonly diagnosed cancer (12.9% of total) but also the most common cause of cancer death (19.4% of total) around the world. In China alone, about 580 000 people die of lung cancer each year. Non-small-cell lung cancer (NSCLC) is the most frequent type of lung cancer, which accounts for nearly 80% of lung cancer cases, causing millions of deaths per year worldwide. Many studies have shown that a significant proportion of NSCLC is associated with genetic mutations or abnormal expression, including epidermal growth factor receptor (EGFR), KRAS proto-oncogene (KRAS), and ALK receptor tyrosine kinase (ALK) genes.

In recent years, due to the improved comprehension of NSCLC molecular mechanism, EGFR tyrosine kinase inhibitors, such as Gefitinib, are applying to clinical treatment, which acquiring significant benefits in overall survival and progression-free survival through inhibiting the tyrosine kinase activity of EGFR by reversibly competing with adenosine triphosphate (ATP) at the ATP-binding site within the EGFR protein.

Although gefitinib can bring significant improvement to NSCLC patients, the vast majority of patients develop resistance after receiving treatment for about a year.[6,7] Studies have shown that a secondary mutation in EGFR (T790M) and amplification of MET were closely related to gefitinib resistance.[7,8] These 2 mechanisms account for about 50% and 20% of the total number of gefitinib resistance cases, respectively. However, the mechanisms of the other cases of gefitinib resistance remain unclear. In addition, few genetic markers that predict gefitinib resistance in lung cancer have been reported. Therefore, it is significant to study the potential molecular mechanisms of gefitinib resistance and to find reliable molecular markers which can predict tumor response to gefitinib and clinical prognosis for NSCLC patients.

Gene expression profiling, also known as microarray analysis has been shown to play an important role in predic-ting tumor prognosis and revealing mechanisms of drug resistance.[10 -12] As most of the gene expression profile data is deposited in public databases, integration and re-analysis of these data may contribute to discovering novel clues for further study. In recent years, several studies have used gene expression profiling to analyze the molecular mechanism of gefitinib resistance.[13,14] The molecular mechanisms found in these studies and biomarkers that predict gefitinib resistance have been shown to improve survival rate and prognosis in patients with NSCLC. Therefore, further research on the molecular mechanism of gefitinib resistance and the exploitation of potential biomarkers can be of great benefit for the treatment of NSCLC.

In this study, we used bioinformatics methods to analyze gene expression profile data in Gene Expression Omnibus (GEO) database to identify differentially expressed genes (DEGs) between gefitinib-sensitive and gefitinib-resistant human lung adenocarcinoma cells. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were conducted to analyze the signaling pathways represented by the differential genes. Furthermore, a protein-protein interaction (PPI) network was constructed to identify hub genes. Bioinformatics analysis of crucial genes and pathways in the current study revealed novel molecular mechanisms of gefitinib resistance, which may assist clinicians to be more targeted in clinical medication.

Materials and Methods

Microarray data

Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo) is a public functional genomic database that stores microarray and high-throughput sequencing data. The gene expression dataset GSE122005, based on GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array and submitted by Wu et al, was downloaded from the GEO database. The GSE122005 dataset has 6 samples, including 3 gefitinib-sensitive lung cancer cell line samples and 3 acquired gefitinib-resistant lung cancer cell line samples.

Identification of DEGs

The GEO2R (https://www.ncbi.nlm.nih.gov/geo/geo2r/) was used to screen DEGs between HCC827 and HCC827/gef cells. GEO2R, which is based on R language limma package, is a web-based tool designed for identification the DEGs between different treatment groups in a GEO series. The probes lacking relevant gene symbols or the genes with more than 1 probe ID were deleted or averaged, respectively. The cutoff criteria: |log2 fold change (FC)| > 2 and adjusted P value < .01 were considered statistically significant. Additionally, ImageGP software (http://www.ehbio.com/ImageGP/index.php/Home/Index/index.html) was applied to present the volcano plot and heat map to visualize DEGs between 2 groups.

KEGG and GO enrichment analyses of DEGs

The Database for Annotation Visualization and Integrated Discovery ([DAVID] https://david.ncifcrf.gov/) is considered as a typical valuable approach to perform the Gene ontology analysis (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG).[16,17] The GO Ontologies offer a systematic terminology for the explanation of characteristics of genes and gene products in 3 key domains that are shared by all organisms, namely molecular function, biological process, and cellular component. KEGG is a database for the methodical study of gene functionalities regarding the systems of genes and molecules. It is a data source intended for comprehending top-level functions of the biological system from molecular-level details, particularly massive molecular datasets produced through genome sequencing as well as other high-throughput measures. To evaluate the function of DEGs, GO, and KEGG analyses were performed through DAVID software. P < .05 was considered as the cut-off value with a statistical difference. R language ggplot2 package is used to visualize related analyses.

Protein-protein interaction (PPI) network construction and module analysis

PPI network of the DEGs was investigated using Search Tool for the Retrieval of Interacting Genes ([STRING], http://www.string-db.org/) and mapped with Cytoscape V_3.7.1 software (San Diego, CA). The STRING database is intended for gathering, rating as well as combining all openly accessible resources of protein-protein interaction data, furthermore, matching these by computational predictions. Additionally, we conducted module analysis through Molecular Complex Detection (MCODE) using the default parameters (K-core = 2, Node score cutoff = 0.2, Max depth = 100, and Degree cutoff = 2). MCODE is a Cytoscape app that detects closely related regions within a specific network grounded on topology.

Hub genes selection and analysis

The hub genes in the network were chosen using CytoHubba, based on connection degree. CytoHubba is a Cytoscape application that predicts and investigates significant hubs and subnetworks in a given system by a few topological calculations. Gene Expression Profiling Interactive Analysis (GEPIA) is a web tool for investigating the RNA sequencing expression information of tumors and relevant normal samples. Many functions could be achieved with GEPIA, for instance, similar gene detection, patient survival analysis, differential expression analysis, correlation analysis and so on. In the current study, the expression level, overall survival and disease-free survival analyses of hub genes were performed utilizing GEPIA.

Gene set enrichment analysis

Gene Set Enrichment Analysis (GSEA, http://software.broadinstitute.org/gsea/index.jsp) is a computational strategy that decides if a previously characterized set of genes indicates statistically significant, concordant differences between 2 biological states (eg, phenotypes). In this study, 6 LUAD samples from GSE122005 were partitioned into 2 groups (high vs low) based on the expression level of CCL20. GSEA was performed between the 2 groups to examine the potential function of CCL20. FDR < .05 and normal P-value < .01 were set as the cut-off criteria.

Results

There were 6 samples in this study, 3 samples were Gefitinib-sensitive lung cancer cell lines, and the other 3 were acquired gefitinib-resistant lung cancer cell lines. We used the GEO2R analysis tool of the GEO database to find differentially expressed genes between the control group and the experimental group, using adjusted P value < .01 and |logFC| ⩾2 as cut-off criteria. A total of 403 differentially expressed genes were obtained, among which 181 were down-regulated while 222 were up-regulated (Figure 1a). The expression level of the top 50 DEGs with fold change 2 was displayed in Figure 1b.

Figure 1.

Analysis of GSE122005. (a) Volcano plot of genes detected in samples. Red means up-regulated DEGs; green means down-regulated DEGs; grey means no difference. (b) Heat map of top 50 DEGs.

Abbreviation: DEGs, differentially expressed genes.

GO function and KEGG pathway enrichment analysis

In order to obtain a more comprehensive and in-depth knowledge of those chosen DEGs, GO function and KEGG pathway enrichment analysis were employed via DAVID. After importing all the DEGs to the DAVID software, we discovered up-regulated DEGs and down-regulated DEGs by GO analysis (Figure 2a). To be more specific, these DEGs were mainly enriched in biological processes (BP), involving signal transduction, positive regulation of GTPase activity, cell adhesion, immune response, inflammatory response, negative regulation of cell proliferation, response to drug, intracellular signal transduction, positive regulation of gene expression, and cellular response to lipopolysaccharide. As for cell component (CC), DEGs were mainly implicated with plasma membrane, extracellular exosome, extracellular space, extracellular region, integral component of plasma membrane, endoplasmic reticulum, cell surface, perinuclear region of cytoplasm, apical plasma membrane and extracellular matrix. In addition, GO molecular function (MF) analysis uncovered that DEGs were principally enriched in calcium ion binding, receptor binding, oxidoreductase activity, signal transducer activity, protein tyrosine kinase activity, phospholipid binding, protease binding, chemokine activity, calcium-dependent phospholipid binding as well as phosphatidylinositol-4,5-bisphosphate 3-kinase activity (Table 1).

Figure 2.

Functional analysis of integrated DEGs. (a) Gene ontology analysis and significant enriched GO terms related with gefitinib resistance. (b) Significantly enriched KEGG pathway terms related with gefitinib resistance.

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

Table 1.

Gene ontology analysis of differentially expressed genes associated with gefitinib resistance in lung adenocarcinoma.

Expression	Category	Term	Description	Count	%	P value
Upregulated	BP	GO:0006955	Immune response	15	0.06	6.34E−05
	BP	GO:0071222	Cellular response to lipopolysaccharide	7	0.03	8.19E−04
	BP	GO:0001558	Regulation of cell growth	6	0.02	.001084
	BP	GO:0006954	Inflammatory response	12	0.05	.001205
	BP	GO:0008285	Negative regulation of cell proliferation	12	0.05	.001707
	BP	GO:0060326	Cell chemotaxis	5	0.02	.003658
	BP	GO:0070098	Chemokine-mediated signaling pathway	5	0.02	.005019
	BP	GO:0006935	Chemotaxis	6	0.02	.006756
	BP	GO:0042060	Wound healing	5	0.02	.007641
	BP	GO:0007165	Signal transduction	21	0.08	.008856
	BP	GO:0071230	Cellular response to amino acid stimulus	4	0.02	.010653
	BP	GO:0008045	Motor neuron axon guidance	3	0.01	.01437
	BP	GO:0070207	Protein homotrimerization	3	0.01	.019048
	BP	GO:0090023	Positive regulation of neutrophil chemotaxis	3	0.01	.019048
	BP	GO:0051216	Cartilage development	4	0.02	.019634
	BP	GO:0032496	Response to lipopolysaccharide	6	0.02	.022005
	BP	GO:0071356	Cellular response to tumor necrosis factor	5	0.02	.022337
	BP	GO:0071347	Cellular response to interleukin-1	4	0.02	.031749
	BP	GO:0044344	Cellular response to fibroblast growth factor stimulus	3	0.01	.034107
	BP	GO:0071395	Cellular response to jasmonic acid stimulus	2	0.01	.03827
	BP	GO:0051414	Response to cortisol	2	0.01	.03827
	BP	GO:0010628	Positive regulation of gene expression	7	0.03	.042811
	BP	GO:0001501	Skeletal system development	5	0.02	.044583
	BP	GO:0009409	Response to cold	3	0.01	.047575
	BP	GO:0048469	Cell maturation	3	0.01	.047575
	BP	GO:0042420	Dopamine catabolic process	2	0.01	.047607
	BP	GO:0043547	Positive regulation of GTPase activity	11	0.04	.049048
	CC	GO:0005615	Extracellular space	34	0.13	2.25E−07
	CC	GO:0005576	Extracellular region	31	0.12	1.70E−04
	CC	GO:0005887	Integral component of plasma membrane	23	0.09	.011618
	CC	GO:0031012	Extracellular matrix	8	0.03	.020731
	CC	GO:0005886	Plasma membrane	51	0.19	.020821
	CC	GO:0005783	Endoplasmic reticulum	15	0.06	.022148
	CC	GO:0005622	Intracellular	21	0.08	.022291
	CC	GO:0048786	Presynaptic active zone	3	0.01	.029464
	MF	GO:0008009	Chemokine activity	5	0.02	.001209
	MF	GO:0005544	Calcium-dependent phospholipid binding	4	0.02	0.017917
	MF	GO:0005509	Calcium ion binding	14	0.05	.019895
	MF	GO:0047086	Ketosteroid monooxygenase activity	2	0.01	.028342
	MF	GO:0047718	Indanol dehydrogenase activity	2	0.01	.028342
	MF	GO:0047115	Trans-1,2-dihydrobenzene-1,2-diol dehydrogenase activity	2	0.01	.03761
	MF	GO:0018636	Phenanthrene 9,10-monooxygenase activity	2	0.01	.03761
	MF	GO:0016491	Oxidoreductase activity	6	0.02	.042989
Downregulated	BP	GO:0008544	Epidermis development	7	0.02	1.21E−04
	BP	GO:0008203	Cholesterol metabolic process	5	0.02	.003354
	BP	GO:0016540	Protein autoprocessing	3	0.01	.005947
	BP	GO:0043410	Positive regulation of MAPK cascade	5	0.02	.00626
	BP	GO:0045599	Negative regulation of fat cell differentiation	4	0.01	.006447
	BP	GO:0007507	Heart development	7	0.02	.006465
	BP	GO:0042552	Myelination	4	0.01	.008305
	BP	GO:0032355	Response to estradiol	5	0.02	.009384
	BP	GO:0042542	Response to hydrogen peroxide	4	0.01	.011025
	BP	GO:0045216	Cell-cell junction organization	3	0.01	.011325
	BP	GO:0050680	Negative regulation of epithelial cell proliferation	4	0.01	.014204
	BP	GO:0016266	O-glycan processing	4	0.01	.017083
	BP	GO:0042632	Cholesterol homeostasis	4	0.01	.020264
	BP	GO:0061098	Positive regulation of protein tyrosine kinase activity	3	0.01	.022949
	BP	GO:0014065	Phosphatidylinositol 3-kinase signaling	3	0.01	.02464
	BP	GO:0030855	Epithelial cell differentiation	4	0.01	.025606
	BP	GO:0035330	Regulation of hippo signaling	2	0.01	.026912
	BP	GO:0060672	Epithelial cell morphogenesis involved in placental branching	2	0.01	.026912
	BP	GO:0042127	Regulation of cell proliferation	6	0.02	.02692
	BP	GO:0071549	Cellular response to dexamethasone stimulus	3	0.01	.028169
	BP	GO:0001843	Neural tube closure	4	0.01	.0327
	BP	GO:0030198	Extracellular matrix organization	6	0.02	.032913
	BP	GO:0014066	Regulation of phosphatidylinositol 3-kinase signaling	4	0.01	.033789
	BP	GO:0070830	Bicellular tight junction assembly	3	0.01	.033814
	BP	GO:0043588	Skin development	3	0.01	.037801
	BP	GO:0006919	Activation of cysteine-type endopeptidase activity involved in apoptotic process	4	0.01	.039512
	BP	GO:2000647	Negative regulation of stem cell proliferation	2	0.01	.044453
	BP	GO:0038128	ERBB2 signaling pathway	3	0.01	.046278
	BP	GO:0006953	Acute-phase response	3	0.01	.048496
	CC	GO:0070062	Extracellular exosome	46	0.14	4.64E−05
	CC	GO:0016323	Basolateral plasma membrane	9	0.03	2.19E−04
	CC	GO:0016324	Apical plasma membrane	11	0.03	2.97E−04
	CC	GO:0031410	Cytoplasmic vesicle	9	0.03	.001282
	CC	GO:0032587	Ruffle membrane	5	0.02	.006282
	CC	GO:0005615	Extracellular space	22	0.07	.008953
	CC	GO:0009986	Cell surface	12	0.04	.009641
	CC	GO:0048471	Perinuclear region of cytoplasm	13	0.04	.010077
	MF	GO:0030296	Protein tyrosine kinase activator activity	3	0.01	.002927
	MF	GO:0004716	Receptor signaling protein tyrosine kinase activity	3	0.01	.003637
	MF	GO:0005543	Phospholipid binding	5	0.02	.00829
	MF	GO:0004872	Receptor activity	7	0.02	.015578
	MF	GO:0046934	Phosphatidylinositol-4,5-bisphosphate 3-kinase activity	4	0.01	.019676
	MF	GO:0005515	Protein binding	94	0.29	.033865
	MF	GO:0004713	Protein tyrosine kinase activity	5	0.02	.034974
	MF	GO:0004364	Glutathione transferase activity	3	0.01	.04138
	MF	GO:0038132	Neuregulin binding	2	0.01	.045365

Abbreviation: GO, gene ontology.

Table 2 displayed the most significantly enriched KEGG pathway of the upregulated and downregulated DEGs. These downregulated DEGs were enriched in ErbB signaling pathway, amoebiasis, renal cell carcinoma, complement and coagulation cascades, and biosynthesis of amino acids, while the upregulated DEGs were enriched in cytokine-cytokine receptor interaction, TNF signaling pathway, alcoholism, NOD-like receptor signaling pathway, amoebiasis, transcriptional misregulation in cancer, chemokine signaling pathway, serotonergic synapse, and adipocytokine signaling pathway. Figure 2b gives a KEGG pathway enrichment plot.

Table 2.

KEGG pathway analysis of differentially expressed genes associated with gefitinib resistance in lung adenocarcinoma.

Expression	Term	Count	%	P value	Genes
Upregulated	hsa04060:Cytokine-cytokine receptor interaction	11	5.392157	3.08E−04	CXCL1, CSF2, TNFRSF1B, TNFSF10, CCL2, CCL20, CRLF2, CXCL2, TNFSF15, CXCL8, IL7R
	hsa04668:TNF signaling pathway	7	3.431373	.0011132	CXCL1, CSF2, TNFRSF1B, CCL2, CCL20, CXCL2, AKT3
	hsa05034:Alcoholism	8	3.921569	.0032513	HIST1H2AC, HIST2H2BE, MAOA, MAOB, GNB4, H2AFJ, PKIA, HIST1H4H
	hsa04621:NOD-like receptor signaling pathway	5	2.45098	.0032727	CXCL1, CCL2, CXCL2, CXCL8, CARD6
	hsa05146:Amoebiasis	6	2.941176	.0061389	CSF2, ARG2, COL3A1, SERPINB2, CXCL8, COL4A6
	hsa05202:Transcriptional misregulation in cancer	7	3.431373	.0100728	PLAT, CSF2, CDKN2C, CXCL8, ETV1, RUNX2, ETV5
	hsa04062:Chemokine signaling pathway	7	3.431373	.0164744	CXCL1, CCL2, CCL20, CXCL2, CXCL8, GNB4, AKT3
	hsa04726:Serotonergic synapse	5	2.45098	.0339256	MAOA, MAOB, GNB4, ALOX5, ITPR1
	hsa04920:Adipocytokine signaling pathway	4	1.960784	.0418645	TNFRSF1B, AKT3, CPT1A, ACSL5
Downregulated	hsa04012:ErbB signaling pathway	5	2.857143	.0127849	PAK6, PAK3, ERBB3, ERBB2, GAB1
	hsa05146:Amoebiasis	5	2.857143	.0246201	IL1R2, LAMA3, PLCB4, ITGB2, COL1A1
	hsa05211:Renal cell carcinoma	4	2.285714	.0314588	PAK6, PAK3, GAB1, EGLN3
	hsa04610:Complement and coagulation cascades	4	2.285714	.0352341	F12, FGB, F3, SERPINA1
	hsa01230:Biosynthesis of amino acids	4	2.285714	.0392337	ASS1, ALDOC, PHGDH, CBS

PPI network construction

The DEG expression products in gefitinib-resistant lung cancer were constructed using the STRING database (http://string-db.org) to construct PPI networks, with a total of 403 DEGs, including 222 upregulated genes and 181 downregulated genes. The network contained 355 nodes and 209 edges. The average node degree was 1.18, and average local clustering coefficient was 0.258. PPI enrichment P-value was 4.3e-07 (Figure 3).

Figure 3.

Protein–protein interaction networks.

The circles represent genes, the structure within the circle is the 3D protein structure, and the connection between the 2 genes indicates protein interaction. The different colors of the lines represent different kinds of protein interactions.

Module screening from the PPI network

In order to detect the most significant modules in this PPI network, we employed the MCODE plug-in. The top 3 modules were selected (Figure 4).

Figure 4.

Top 3 modules from the protein-protein interaction network: (a) module 1, (b) module 2, and (c) module 3.

Hub gene selection and analysis

Based on the information of the STRING protein query from public databases, we constructed the PPI network of the top 10 hub genes according to the degree of connectivity (Figure 5). The top 10 hub genes with a higher degree of connectivity are as follows: GNB4, LPAR1, CXCL1, PMCH, NMU, CXCL8, CXCL2, GPER1, CCL20, CXCL16. The names, abbreviations and functions for these hub genes are shown in Table 3.

Figure 5.

Protein-protein interaction network of the top 10 hub genes.

Table 3.

Functional roles of 10 hub genes with degree ⩾10.

No.	Gene symbol	Full name	Function
1	GNB4	G protein subunit beta 4	Heterotrimeric guanine nucleotide-binding proteins (G proteins), which integrate signals between receptors and effector proteins, are composed of an alpha, a beta, and a gamma subunit. These subunits are encoded by families of related genes. This gene encodes a beta subunit. Beta subunits are important regulators of alpha subunits, as well as of certain signal transduction receptors and effectors.
2	LPAR1	Lysophosphatidic acid receptor 1	The integral membrane protein encoded by this gene is a lysophosphatidic acid (LPA) receptor from a group known as EDG receptors. These receptors are members of the G protein-coupled receptor superfamily. Utilized by LPA for cell signaling, EDG receptors mediate diverse biologic functions, including proliferation, platelet aggregation, smooth muscle contraction, inhibition of neuroblastoma cell differentiation, chemotaxis, and tumor cell invasion. Two transcript variants encoding the same protein have been identified for this gene
3	CXCL1	C-X-C motif chemokine ligand 1	This antimicrobial gene encodes a member of the CXC subfamily of chemokines. The encoded protein is a secreted growth factor that signals through the G-protein coupled receptor, CXC receptor 2. This protein plays a role in inflammation and as a chemoattractant for neutrophils. Aberrant expression of this protein is associated with the growth and progression of certain tumors. A naturally occurring processed form of this protein has increased chemotactic activity. Alternate splicing results in coding and non-coding variants of this gene. A pseudogene of this gene is found on chromosome 4.
4	PMCH	Pro-melanin concentrating hormone	This gene encodes a preproprotein that is proteolytically processed to generate multiple protein products. These products include melanin-concentrating hormone (MCH), neuropeptide-glutamic acid-isoleucine (NEI), and neuropeptide-glycine-glutamic acid (NGE). Melanin-concentrating hormone is a 19-amino acid neuropeptide that stimulates hunger and may additionally regulate energy homeostasis, reproductive function, and sleep. Pseudogenes of this gene have been identified on chromosome 5.
5	NMU	Neuromedin U	This gene encodes a member of the neuromedin family of neuropeptides. The encoded protein is a precursor that is proteolytically processed to generate a biologically active neuropeptide that plays a role in pain, stress, immune-mediated inflammatory diseases and feeding regulation. Increased expression of this gene was observed in renal, pancreatic and lung cancers. Alternative splicing results in multiple transcript variants encoding different isoforms. Some of these isoforms may undergo similar processing to generate the mature peptide.
6	CXCL8	C-X-C motif chemokine ligand 8	The protein encoded by this gene is a member of the CXC chemokine family and is a major mediator of the inflammatory response. The encoded protein is secreted primarily by neutrophils, where it serves as a chemotactic factor by guiding the neutrophils to the site of infection. This chemokine is also a potent angiogenic factor. This gene is believed to play a role in the pathogenesis of bronchiolitis, a common respiratory tract disease caused by viral infection. This gene and other members of the CXC chemokine gene family form a gene cluster in a region of chromosome 4q.
7	CXCL2	C-X-C motif chemokine ligand 2	This antimicrobial gene is part of a chemokine superfamily that encodes secreted proteins involved in immunoregulatory and inflammatory processes. The superfamily is divided into four subfamilies based on the arrangement of the N-terminal cysteine residues of the mature peptide. This chemokine, a member of the CXC subfamily, is expressed at sites of inflammation and may suppress hematopoietic progenitor cell proliferation.
8	GPER1	G protein-coupled estrogen receptor 1	This gene encodes a multi-pass membrane protein that localizes to the endoplasmic reticulum and a member of the G-protein coupled receptor 1 family. This receptor binds estrogen and activates multiple downstream signaling pathways, leading to stimulation of adenylate cyclase and an increase in cyclic AMP levels, while also promoting intracellular calcium mobilization and synthesis of phosphatidylinositol 3,4,5-trisphosphate in the nucleus. This protein therefore plays a role in the rapid nongenomic signaling events widely observed following stimulation of cells and tissues with estrogen. This receptor has been shown to play a role in diverse biological processes, including bone and nervous system development, metabolism, cognition, male fertility and uterine function.
9	CCL20	C-C motif chemokine ligand 20	This antimicrobial gene belongs to the subfamily of small cytokine CC genes. Cytokines are a family of secreted proteins involved in immunoregulatory and inflammatory processes. The CC cytokines are proteins characterized by two adjacent cysteines. The protein encoded by this gene displays chemotactic activity for lymphocytes and can repress proliferation of myeloid progenitors. Two transcript variants encoding different isoforms have been found for this gene.
10	CXCL16	C-X-C motif chemokine ligand 16	Acts as a scavenger receptor on macrophages, which specifically binds to OxLDL (oxidized low density lipoprotein), suggesting that it may be involved in pathophysiology such as atherogenesis (by similarity). Induces a strong chemotactic response. Induces calcium mobilization. Binds to CXCR6/Bonzo.

We obtained the prognostic information of the top 10 hub genes in http://gepia.cancer-pku.cn/. It was demonstrated that expression of CCL20 (HR high = 1.4, logrank P = .021) was associated with worse overall survival (OS) and disease-free survival (DFS) for LUAD patients (Figure 6). However, the other hub genes had no obvious difference on the survival curve of LUAD patients.

Figure 6.

Prognostic value of CCL20 gene in LUAD patients.

Then, we verified the expression level of 10 hub genes in lung tissues between LUAD and healthy people using GEPIA, and Figure 7 showed that compared to the normal group, the expression level of NMU and CCL20 significantly elevated in LUAD patients, while CXCL2 and GPER1 significantly down-regulated. Intriguingly, only CCL20 had a significant influence on the prognosis of LUAD patients.

Figure 7.

Expression level of (a) CCL20, (b) CXCL2, (c) GPER1, and (d) NMU in LUAD and normal tissues.

Subsequently, the expression profile of CCL20 in human tissue was displayed using GEPIA. We found that CCL20 mRNA in esophagus, lung, liver, cholangio, pancreas, stomach and colon tumors displayed higher levels as compared with the matched normal tissues (Figure 8).

Figure 8.

Expression profile of CCL20 in human tissues (tumor vs normal tissues).

Gene set enrichment analysis (GSEA)

To acquire further insight into the function of hub gene CCL20, GSEA was conducted to map into GO analysis and KEGG pathways database. According to the median of CCL20 gene expression value, we divided the samples into high expression group and low expression group.

Using GO as gene sets database, GSEA analysis suggests high expression of CCL20 is enriched in chemokine mediated signaling pathway, DNA packaging complex, CCR chemokine receptor binding, proximal-distal pattern formation, nuclear nucleosome, chemokine activity, taste receptor activity, monocyte chemotaxis, negative regulation of leukocyte migration, bitter taste receptor activity, water transport, chromatin silencing, protein DNA complex, chromatin silencing at rDNA, regulation of lymphocyte migration, cell-cell recognition (Figure 9). While, low expression of CCL20 is enriched in gas transport, cornified envelope, neuropeptide receptor activity, keratin filament, neuropeptide signaling pathway. Using KEGG as gene sets database, GSEA analysis suggests high expression of CCL20 is enriched in systemic lupus erythematosus. While in low expression of CCL20, no pathways are significantly enriched.

Figure 9.

Gene set enrichment analysis (GSEA). Listed pictures are 16 functional gene sets enriched in CCL20 highly expressed groups.

Discussion

Currently, gefitinib is the first-line drug for the treatment of locally advanced or metastatic NSCLC patients that have previously received chemotherapy. Therefore, gefitinib resistance is a global problem that needs to be solved urgently. While many previous related studies have focused on a single genetic mutation site, this study is based on a holistic perspective. We analyzed in depth the differences in gene expression between gefitinib-sensitive and gefitinib-resistant NSCLC samples and found 403 DEGs that were significantly related to gefitinib-resistance. After that, both KEGG and GO analyses were carried out to illustrate the underlying biological functions and pathways associated with the gefitinib resistance. GO enrichment in our study indicated that upregulated DEGs were mainly enriched in inflammatory response and cell chemotaxis.

Interstitial lung disease (ILD) is a rare, but fatal, serious adverse reaction caused by gefitinib and it could be resolved by applying high doses of glucocorticoids. The period from beginning gefitinib treatment to the outbreak of ILD was short. It is conjectured that a Th1 type of lung tissue inflammation could be associated with gefitinib-induced ILD. Inoue et al reported that pretreatment with gefitinib in the mouse enhanced LPS-induced lung inflammation through inhibiting the expression of SP-A on the alveolar type II cells, therefore, raises the sensitivity to pathogens, which might explain a part of the underlying mechanisms of gefitinib-induced ILD.

Most patients accepting EGFR inhibitors treatment usually experience skin toxicity, such as papulopustular rash. Kanazawa et al reported that lymphocytes were activated to release inflammatory cytokine in the skin after gefitinib treatment, which evoked and regulated tumor-specific T-cell immunity in the NSCLC patients. Inflammatory cytokines such as PGE2 in the tumor microenvironment may promote EGFR tyrosine kinase inhibitors (TKI) resistance in NSCLC by reducing E-cadherin expression. In order to improve the effectiveness of EGFR-TKI, it is essential to discover the particular inflammatory signals regulating the process of epithelial to mesenchymal transition (EMT).

In addition, chronic inflammation is also closely related to the occurrence and development of tumors, which has become a consensus. Smoking is the major etiologic risk factor for lung carcinoma, and growing proof revealed that cigarette smoke exposure facilitates extensive inflammatory and mutagenic reactions in the lungs which initiate a pro-cancer immune response. Breathing harmful particulates from tobacco smoke results in recurrent damage to lung epithelial cells, which causes continuous recruitment of host inflammatory cells to pulmonary airways and alveolus. The lungs of tobacco users are featured with higher amounts and activity of activated lymphocytes, macrophages, dendritic cells, and granulocytes, consequently generating an atmosphere for the multiplication of malignantly transformed cells.

In non-smokers, inflammation is also closely related to the development of lung cancer. Patients with a history of tuberculosis have a 50% increased risk of developing lung cancer, and in male patients, lung tumors tend to occur on the same side of previous tuberculosis infections. Large-scale research based on population performed in China suggested a positive correlation between previous pulmonary tuberculosis and lung cancer, which was independent of smoking status. Accumulating evidence suggests that human papillomavirus (HPV) in lung tissue contributes significantly to lung carcinogenesis. Studies in Taiwan have shown that the detection rate of HPV 6 in lung tissues of lung cancer patients (28.4%) is significantly higher than healthy people (1.7%), indicating that HPV 6 virus infection has a correlation with lung cancer. HPV 33 was more prevalent in patients with NSCLC in Korea than HPV 16 and HPV 18.

Cell chemotaxis is a directional response of cells to chemical stimuli in the external environment. Goswami et al reported that compared to stationary tumor cells, EGF-responsive chemotactic tumor cells are more resistant to chemotherapeutic drugs. By applying human cDNA array technology, Duan et al found that the progression of paclitaxel resistance is along with various modifications in gene expression involving consistent changes in particular chemokine and cytokine expression. Moreover, chemotaxis of tumor cells in the neighboring microenvironment plays a key role in tumor propagation during development and transition. KEGG pathway analysis further confirmed that up-regulated DEGs are associated with the chemokine signaling pathway.

Gene ontology study further suggested that the downregulated DEGs were mainly associated with the modification of phosphatidylinositol kinase involving phosphatidylinositol-4,5-bisphosphate 3-kinase activity and regulation of phosphatidylinositol 3-kinase (PI3K) signaling. The PIK3CA gene encodes the catalytic subunit of PI3K. It mutates or amplifies in a variety of tumors, including lung cancer. Mutations in the PI3KA gene are often co-existing with mutations in the KRAS and EGFR genes. Meanwhile, shorter overall survival in patients with single PIK3CA mutation has been identified compared with wildtype PIK3CA and PIK3CA-EGFR/KRAS co-mutation. Taken together, Gene ontology analysis revealed several critical biological processes, cellular components and molecular functions which could be associated with the tumor progression and gefitinib resistance in NSCLC.

KEGG pathway enrichment uncovered that up-regulated DEGs were enriched in cytokine-cytokine receptor interaction, TNF signaling pathway, NOD-like receptor signaling pathway, and so forth.

Cytokines are molecules secreted by cells, including chemokines, growth factors, extracellular proteases, and the like. Cytokines are closely related to the occurrence, development and metastasis of tumors. Recent studies have shown that cytokines play an important role in drug resistance and tumor progression. Vascular endothelial growth factor (VEGF), also known as vascular permeability factor (VPF), is the most potent and specific angiogenic factor inducing tumor angiogenesis. VEGF is secreted by different tumor cells, regulates the permeability of blood vessels, and promotes the migration and proliferation of vascular endothelial cells. Different methods of intervention or inhibition of VEGF can inhibit tumor growth.[41,42]

Tumor necrosis factor (TNF) is the first cytokine used in tumor biotherapy, mainly produced by activated macrophages, NK cells and T lymphocytes. TNF kills certain tumor cells in vivo and in vitro or inhibits proliferation. TNF is a pro-inflammatory cytokine and plays an important role in a variety of cellular processes, such as cell necrosis, apoptosis, and differentiation. Most TNF-induced cellular responses are mediated by 2 receptors, TNF-R1 and TNF-R2.

NOD-like receptor (NLR) belongs to the intracellular pattern recognition receptor and consists of 3 parts: the N-terminal effector domain binds to downstream effector molecules, the middle oligomeric domain mediates its own oligomerization, and the C-terminal leucine-rich repeat sequence (LRR) is able to recognize ligands. So far, studies found that the human NLR family contains 23 members. NLR plays an important role in the process of immune recognition and pro-inflammatory factors production, and can lead to the occurrence of certain diseases. NLRP12, belongs to NLR family, can reduce the expression of inflammatory cytokines. By knocking out NLRP12 in mice, Zaki et al found that this receptor plays a critical role in the development of colon tumors and colitis. In fact, mice lacking NLRP12 are highly susceptible to colon infection and colon cancer. In addition, NLRP6 has also been found to inhibit inflammation in colon and colitis-associated tumor by regulating tissue repair. Therefore, NLRP6 and NLRP12 are promising drug targets through which drugs can be developed to prevent colon inflammation and colon tumorigenesis.

In addition, the down-regulated DEGs are enriched in the KEGG pathway such as ErbB signaling pathway, complement and coagulation cascades. ErbB is a typical receptor tyrosine kinase with 4 members in the family. ErbB2 is associated with the development of a variety of human cancers. In particular, ErbB2 was found to have consistent amplification in 20% of sporadic breast tumors, and this was the first gene found to be widely mutated in breast cancer. Drug development targeting ErbB receptors has been underway for many years, and drugs have been marketed and are effective. Among them, gefitinib and erlotinib are receptor tyrosine kinase inhibitors targeting ErbB, which are used to treat EGFR-mutated NSCLC. In addition, studies have shown that the ErbB receptor family of cisplatin-resistant ovarian cancer cells shows a change in their ligand response, and the gene expression of the ErbB receptor family has changed, partly explaining the occurrence mechanisms of cisplatin resistance.

Complement is a serum protein found in human and vertebrate serum and tissue fluids that mediates immune and inflammatory responses. Studies have shown that tumor survival is always accompanied by activation of the complement generation pathway. Complement recognizes and destroys tumor cells, but tumor cells have multiple mechanisms to inhibit the toxic effects of complement. For example, tumor cells secrete soluble complement inhibitors, and tumor cell surface expresses complement regulatory proteins, which inhibit complement activity. Additionally, Kubisch et al found that tumor resistance in cyclophosphamide treatment is closely related to coagulation-related gene expression.

PPI refers to the process by which 2 or more protein molecules form a protein complex through non-covalent bonds. PPI analysis can provide us with a direct physical connection between proteins. A holistic study of the PPI network will also obtain knowledge about the approaches in the area of drug discovery. From module screening of PPI network, we uncovered that the interactions among the proteins in gefitinib-resistant lung cancer cells mainly concentrated on pathways implicated with inflammation and chemotaxis, which is consistent with GO enrichment. Chronic infections and inflammation caused by bacteria or viruses are crucial factors in tumorigenesis and progression. The pro-inflammatory properties of chemokines that persist in specific sites can result in the advancement of chronic inflammation. Studies have shown that the nuclear factor of kappa light polypeptide gene enhancer in B cells (NF-κB) complexes play a significant role in chronic inflammation, and activation of the NF-κB pathway leads to chemoresistance of tumor cells.

Meanwhile, several hub genes had been identified base upon the degree via PPI network building and module screening. An in-depth study of these genes may provide new insights into gefitinib resistance. Half of these hub genes are chemokine genes, including CCL20, CXCL1, CXCL2, CXCL8, CXCL16.

CCL20 (C-C motif chemokine ligand 20) is a subfamily of the small cytokine CC gene. Chemokines secreted by tumor cells can recruit lymphocytes into the tumor microenvironment, allowing the tumor to form an immune tolerance state. CCL20 gene consists of 3 introns and 4 exons. CCL20-encoded proteins are capable of chemotaxis of lymphocytes and inhibit the proliferation of myeloid progenitors. The mRNA and protein of CCL20 are abundantly expressed in breast carcinoma samples, while it is not expressed in normal breast epithelial cells. Consistent with our findings, CCL20 has been proved to serve as a drug-resistance promoting factor in several studies recently. A study of triple-negative breast cancer (TNBC) showed that paclitaxel-induced CCL20 protein mediates chemoresistance by up-regulating ABCB1. In addition, activation of the NF-κB pathway increases the expression of CCL20. Thus, blocking CCL20 or its downstream pathway may be a new approach to address paclitaxel resistance in breast cancer patients. Liu et al reported that CCL20 was up-regulated significantly in cisplatin-stimulated classically activated macrophages (CAMs) thus promoting ovarian cell cancer migration. Meanwhile, CCL20’s receptor CCR6 was also activated which leads to a poor prognosis in clinical. Additionally, blockade of CCL20 and CCR6 effectively reduced the cell migration induced by cisplatin-stimulated CAMs. In a recent study, researchers found the expression of CCL20 is closely related to terminal stages of tumor and decreased survival rate. CCR6—the specific receptor of CCL20 was found to express abundantly in microvascular endothelial cells, which contributes to angiogenesis of tumors. On the flip side, CCR6-deficient mice exhibited significantly reduced tumor growth as well as tumor-associated vascularisation. In patients with NSCLC, the CCL20 gene and protein are overexpressed, and autocrine of CCL20 can promote the migration and proliferation of lung cancer cells. However, the relationship between CCL20 and gefitinib resistance has not been reported. We found that CCL20 expression was significantly up-regulated in gefitinib-resistant lung cancer cells and survival analysis showed that LUAD patients with high expression of CCL20 tend to have shorter survival time. It suggests that up-regulation of CCL20 is associated with gefitinib resistance, and CCL20 may serve as a biomarker for predicting gefitinib resistance. Furthermore, the development of inhibitors for CCL20 might become an alternative to resolve resistance to gefitinib and other chemotherapeutic agents and further improve cancer therapy.

GSEA results also showed that chemokine-related signaling pathways were significantly enriched in samples with high CCL20 expression. However, further experiments are needed to verify the role of CCL20 in causing gefitinib resistance in NSCLC patients.

GPER1 encodes a G protein-coupled receptor that binds to estrogen and activates multiple downstream signaling pathways. Estrogen binding to the GPER receptor results in the disintegration of heterotrimeric G protein aggregates. High expression of GPER in the breast leads to acquired resistance to tamoxifen, which is an agonist of GPER. In TNBC, GPER is usually highly expressed and is closely related to tumor recurrence. Therefore, the high expression of this gene can be used as a marker for the poor prognosis of breast tumor patients. In addition, studies have shown that the use of gefitinib reduces the expression of GPCRs in TNBC by inhibiting EGFR, and further inhibits the growth of tumor cells. However, no studies have reported the role of GPER1 in gefitinib resistance in patients with LUAD. We found that GPER1 expression was elevated in gefitinib-resistant lung cancer cells and was not associated with the prognosis of LUAD.

CXCL16 encodes a chemokine that binds to CXCR6 and activates CD8 effector T cells. In tumor cells receiving radiation therapy, CXCL16 expression is up-regulated and binds to the CXCR6 receptor on T cells, causing activated T cells to aggregate to the site of the irradiated tumor. We found that CXCL16 expression was down-regulated in gefitinib-resistant lung cancer samples. This is consistent with the findings of Xu et al, who found that CXCL16 was down-regulated 67 times in erlotinib-resistant lung cancer cells relative to erlotinib-sensitive lung cancer cells.

CXCL2 has also been shown to be associated with tumor resistance. The researchers found that exogenously added CXCL2 significantly offset the inhibition of lung cancer cell migration induced by anlotinib and attenuated lung cancer cell apoptosis induced by anlotinib. Interestingly, gefitinib reduced the expression of CXCL2 in gastric epithelial cells infected by Helicobacter pylori, and we found that CXCL2 expression was up-regulated in gefitinib-resistant lung cancer cells, indicating that up-regulation of CXCL2 can predict resistance to gefitinib after treatment.

CXCL1 is a paralog of CXCL2 and is capable of recruiting neutrophils in inflammation. The receptor for CXCL1 is CXCR2, both of which are expressed in human ovarian tumors. The expression of CXCR2 promotes the expression of proinflammatory chemokines such as CXCL1, CXCL2, CXCL3, CXCL6. Gefitinib can inhibit the proliferation of ovarian cancer cells induced by CXCL1. Like CXCL2, gefitinib also reduced CXCL1 expression in gastric epithelial cells infected with Helicobacter pylori. In our study, we found that CXCL1 expression was up-regulated in gefitinib-resistant cells.

CXCL8 is the main mediator of the inflammatory response and is thought to promote the growth of a variety of tumor cells. Normal cells generally do not express CXCL8 and its receptors, but when infected by viruses or bacteria, CXCL8 and its receptors are induced to express. Initially, CXCL8 was considered to be an inflammatory chemokine, which was later shown to promote angiogenesis. As an angiogenic factor and mitogenic factor, CXCL8 promotes the progression of many tumors. CXCL8 binds to its receptors CXCR1 and CXCR2 to exercise biological functions, and its receptor belongs to the 7 transmembrane G-protein-coupled receptor (GPCR) family. Studies have shown that gefitinib can inhibit tumor cells by blocking the expression of CXCL8. Our study shows that CXCL8 expression is up-regulated in resistant cells compared to gefitinib-sensitive lung cancer cells. Therefore, blocking the expression of CXCL8 and its receptor in cancer cells may alleviate gefitinib resistance.

Conclusions

In summary, we comprehensively analyzed the gene expression profiles of gefitinib-sensitive and resistant lung cancer samples, which will help us better understand the mechanism of occurrence and development of gefitinib resistance in NSCLC patients. The hub genes we analyzed provide several novel biomarkers for the prognosis of lung cancer patients treated with gefitinib. In the meantime, they also provide potential therapeutic targets for managing gefitinib resistance.