Decreased HLF Expression Predicts Poor Survival in Lung Adenocarcinoma.

BACKGROUND Lung adenocarcinoma (LUAD) is a type of non-small cell carcinoma. Its pathogenesis is being explored and there is no cure for the disease. MATERIAL AND METHODS The Gene Expression Omnibus (GEO) was searched to obtain data on expression of messenger RNA. GEO2R, an interactive web tool, was used to calculate the differentially expressed genes (DEGs) in LUAD. All the DEGs from different datasets were imported into VENNY 2.1 (https://bioinfogp.cnb.csic.es/tools/venny/index.html) to identify the intersection of the DEGs. An online analysis tool, the Database for Annotation, Visualization, and Integrated Discovery (DAVID), was used to help understand the biological meaning of DEG enrichment in LUAD. Cytoscape 3.7.2 was used to perform centrality analysis and visualize hub genes and related networks. Furthermore, the prognostic value of the hub genes was evaluated with the Kaplan-Meier plotter survival analysis tool. RESULTS The GEO database was used to obtain RNA sequencing information for LUAD and normal tissue from the GSE118370, GSE136043, and GSE140797 datasets. A total of 376 DEGs were identified from GSE118370, 248 were identified from GSE136403, and 718 DEGs were identified from GSE140797. The 10 genes with the highest degrees of expression - the hub genes - were CAV1, TEK, SLIT2, RHOJ, DGSX, HLF, MEIS1, PTPRD, FOXF1, and ADRB2. In addition, Kaplan-Meier survival evaluation showed that CAV1, TEK, SLIT2, HLF, MEIS1, PTPRD, FOXF1, and ADRB2 were associated with favorable outcomes for LUAD. CONCLUSIONS CAV1, TEK, SLIT2, HLF, MEIS1, PTPRD, FOXF1, and ADRB2 are hub genes in the DEG interaction network for LUAD and are involved in the development of and prognosis for the disease. The mechanisms underlying these genes should be the subject of further studies.

Background

Non-small cell lung carcinoma (NSCLC) includes large cell carcinoma, squamous cell carcinoma (SCC), and adenocarcinoma [1]. Unlike squamous cell lung cancer, lung adenocarcinoma (LUAD) is more common in non-smokers and women [2]. The incidence of LUAD is lower than for SCC and undifferentiated cancer and patients diagnosed with it are younger [2]. In general, there are no obvious clinical symptoms in the early stage of the disease and it often is diagnosed on a chest X-ray. LUAD tends to spread through the bloodstream in early stages, whereas metastasis through the lymphatic system occurs much later [3].

The pathogenesis of lung cancer is complex and involves multiple causes and genes [4]. Many studies are focusing on the pathological mechanism of LUAD [1,4]. However, the cause of lung cancer is still unclear, and there is no cure for it [2]. Exploring the process of LUAD at the gene and protein level can provide detailed and useful information for a comprehensive understanding of the disease.

Based on big data, bioinformatics methods are an increasingly popular way of exploring the role of differentially expressed genes (DEGs) in disease [5]. Hub genes can be identified by analyzing the DEG interaction network. Bioinformatics has been used to identify some diseases, especially in oncology [6]. However, few studies have combined assessment of DEG interaction networks and survival analysis for LUAD. Therefore, the aim of the present study was to use bioinformatics methods to retrieve the DEGs for LUAD, analyze the interaction network of the DEGs, identify hub genes in LUAD, and perform a Kaplan-Meier survival analysis.

Material and Methods

Search of Gene Chip Data

The Gene Expression Omnibus (GEO) is a gene expression database created and maintained by the National Center for Biotechnology Information of the National Library of Medicine in the United States. We searched this database for gene chip data that compared tissue from LUAD with normal tissue. The keywords used were “lung adenocarcinoma.” Matrix data on messenger RNA (mRNA) expression were downloaded for analysis.

Identification of Differentially Expressed Genes

GEO2R is online software used to analyze data from GEO. We imported the data into the tool and analyzed it with the default settings. The DEG information was provided by GEP2R, which sorted the DEGs with LogFC. P<0.05 and LogFC >2 were considered statistically significant. DEGs from different datasets were imported into Venn. Venn was used to identify DEGs that were common to the different datasets ().

Gene Ontology and Kyoto Encyclopedia of Genes and Genomes Pathway Enrichment Analysis of Genes

We used the Database for Annotation, Visualization and Integrated Discovery (DAVID) 6.8 to explore the biological roles of involved genes [7,8]. Gene ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were performed with DAVID. DEGs were imported into DAVID with the setting “gene list” and “gene symbol,” and the associated GO and KEGG terms were obtained. The top 3 terms for biological process, cellular component, and molecular function, and the KEGG pathway with the highest P value were selected.

Protein–Protein Interaction Network Analysis

The Protein–Protein Interaction (PPI) network can be used to show the functional relationship between common DEGs. Upregulated and downregulated genes were imported into the Search Tool for the Retrieval of Interacting Genes (STRING) to establish the network. The file “string_interactions.tsv” was downloaded. Then the PPI network “string_interactions.tsv” was imported into Cytoscape software (version 3.7.2) and visualized [9]. Using the software Centiscape 2.2 in Cytoscape, the DEGs were ranked by degree centrality. That information for the top 10 DEGs in the network was visualized with ggplot2 [10]. The enrichment results for the top 10 DEGs and the GO and KEGG terms were visualized with GOplot [11].

Kaplan-Meier Survival Analysis of Hub Genes

All the hub genes were imported into the Kaplan-Meier plotter web tool () [12]. The Kaplan-Meier plotter can assess the effect of 54 000 genes on survival in 21 cancer types, helping us discover and validate survival biomarkers [13]. Kaplan-Meier survival analysis results for each hub gene were obtained with default settings: “lung adenocarcinoma,” “auto select best cutoff,” and “all the follow-up threshold.” P<0.01 was defined as a statistically significant difference.

Cell Culture and Transfection

The LUAD cell line A549 was purchased from the Shanghai Cell Bank. Dulbecco’s modified Eagle’s medium and 10% fetal bovine serum were used to culture all the cells at 37°C in 5% CO2. For transfection, the negative control small interfering RNA (siRNA) and HLF siRNA were designed by Suzhou GenePharma Biotechnology Co., Ltd. A549 cells were transfected with 100 nM NC siRNA or HLA siRNA using Lipofectamine® 3000.

Quantitative Reverse Transcription Polymerase Chain Reaction

Invitrogen TRIzol was used to extract the total RNA from the cells and tissues. Quantitative real-time reverse transcription polymerase chain reaction (qRT-PCR) was performed to determine the HLF mRNA level. The primers were designed by Shanghai Bio Tech, as follows:

HLF: forward, 5′-CGGAGGTGTTCTATGAGCTGG-3′ and

reverse, 5′-AGCTTGTGTGTTCGCAGGAA-3′; and

GAPDH: forward, 5′-ACTGCGAATGGCTCATTAAATCA-3′ and

reverse, 5′-AGCTCTAGAATTACCACAGTTATCCAAGT-3′.

The 2−ΔΔCq approach was used to quantify the HLF mRNA level.

Results

Identification of DEGs

Gene expression profiles from the GSE118370, GSE136043, and GSE140797 datasets were obtained from the GEO database. GSE118370 is a LUAD dataset; GSE136043 is an mRNA microarray from lung cancer not otherwise specified, and GSE140797 is 7 pairs of LUAD tissue and normal tissue. There were 6 LUAD and 6 normal specimens in GSE118370, 5 LUAD and 5 normal in GSE136043, and 7 LUAD and 7 normal specimens in GSE140797 (Table 1). Next, 245 upregulated and 686 downregulated DEGs were identified in GSE118370; 617 upregulated and 687 downregulated DEGs were identified in GSE136043; and 1073 upregulated and 1281 downregulated DEGs were identified in GSE140797. Finally, 35 common upregulated and 142 common downregulated DEGs were obtained from the 3 datasets (Figure 1).

Table 1

Information about the 3 datasets obtained from GEO.

ID	LUAD	Normal	Total
GSE118370	6	6	12
GSE136043	5	5	10
GSE140797	7	7	14

GEO – Gene Expression Omnibus; LUAD – lung adenocarcinoma.

Figure 1

Common differentially expressed genes (DEGs) in GSE118370, GSE136043, and GSE140797. (A) There were 35 upregulated DEGs common to the 3 datasets. (B) There were 142 downregulated DEGs common to the 3 datasets.

Gene Ontology and Kyoto Encyclopedia of Genes and Genomes Enrichment Analysis of DEGs

For all the common DEGs, GO and KEGG analysis results showed that during biological processes, common DEGs were significantly enriched in cell adhesion, morphogenesis of a branching structure, and angiogenesis. In molecular function, common DEGs were significantly enriched in protein binding, lipid binding, transcriptional activator activity, and RNA polymerase II core promoter proximal region sequence-specific binding. In the cellular component, common DEGs were significantly enriched as an integral component of the plasma membrane, the membrane raft, and the plasma membrane. In the KEGG pathway, common DEGs were significantly enriched in extracellular matrix-receptor interaction, protein digestion and absorption, and axon guidance (Table 2).

Table 2

Functional and pathway enrichment analysis of the DEGs.

Term	Name	Count	P value	Genes
Biological processes
GO: 0007155	Cell adhesion	17	9.7E-6	TNXB, PODXL, FERMT1, SEMA5A, KIAA1462, WISP2, CD36, LAMA3, COL6A6, ITGA8, CD24, PDZD2, MFAP4, EMP2, AOC3, RS1, MYH10
GO: 0001763	Morphogenesis of a branching structure	4	4.8E-5	TCF21, COL13A1, FOXF1, NOTCH4
GO: 0001525	Angiogenesis	10	3.2E-4	TAL1, CAV1, EMCN, HLF, TEK, TBX4, ROBO4, RASIP1, TMEM100, MEIS1
Molecular functions
GO: 0005515	Protein binding	94	1.1E-2	RHOJ, STEAP3, NDST1, ELF3, TNNC1, EIF5, FERMT2, SPINK1, PTPN21, ASPA, DGSX, HHIP, TMEM100, SLC1A1, GOLM1, AHNAK, COL10A1, GHR, LIMK1, TOX3, SLIT2, NEBL, SH2D3C, EYA4, SDC1, ADRB2, CD36, PTRF, SGCG, PLEKHH2, CFL2, CLIC5, SIX1, RYR2, JAM2, MFAP4, EMP2, AOC3, CAV2, TSHZ3, CAV1, HMGB3, ALDH18A1, ADAMTSL3, ABLIM3, PEAK1, BEX1, MME, AFAP1L1, MEIS1, TIMP3, STARD13, TAL1, GPM6A, HOXA5, DISP1, MYZAP, TEK, P2RY1, ETV1, MSI2, VMP1, CD24, FGD5, MYOC, ICA1, PTPRD, ERG, ADARB1, TNXB, HLF, TBX2, COL13A1, PODXL, SMAD6, MAL, DACH1, CDC25C, AGER, SHANK3, CCDC68, LIN7A, TMPRSS4, DOCK4, CDKN1C, UACA, NOTCH4, DEPTOR, SYNM, SH3D19, GRK5, ARAP3, AGR2, MYH10
GO: 0008289	Lipid binding	6	1.2E-2	APOL3, CD36, STARD9, STARD8, MAL, STARD13
GO: 0001077	Transcriptional activator activity, RNA polymerase II core promoter proximal region sequence-specific binding	7	1.9E-2	TCF21, ELF3, HLF, HOXA5, SIX1, ETV1, MEIS1
Cellular component
GO: 0005887	Integral component of plasma membrane	29	7.7E-5	STEAP3, CAV2, CAV1, MME, CD52, SLC19A3, PTGIR, DGSX, TEK, P2RY1, SLC39A8, HHIP, SLC1A1, GOLM1, GHR, PTPRD, PODXL, MAL, AGER, CLEC1A, ADRB2, SDC1, CD36, SSTR1, P2RY14, NOTCH4, ST14, ABCC3, JAM2
GO: 0045121	Membrane raft	10	1.2E-4	CAV2, CAV1, CD36, PTRF, PODXL, TEK, MAL, CD24, EMP2, AHNAK
GO: 0005886	Plasma membrane	59	1.9E-4	RHOJ, EIF5, DUOX1, PTGIR, DGSX, KCNQ3, TMEM100, SCN7A, DPP6, SLC1A1, AHNAK, GHR, SLIT2, ADRB2, SDC1, CD36, SGCG, STXBP6, PLEKHH2, SSTR1, CLIC6, ST14, RYR2, JAM2, EMP2, AOC3, CAV2, EMCN, CAV1, CLDN18, PEAK1, MME, SLC19A3, GNG11, AFAP1L1, SEMA5A, STK32A, GPM6A, P2RY1, TEK, SLC39A8, VMP1, PTPRD, COL13A1, PODXL, AGER, SHANK3, GJB2, DOCK4, LIN7A, P2RY14, ITGA8, NOTCH4, ABCC3, NLN, SH3D19, GRK5, ARAP3, MYH10
KEGG pathway
hsa04512	ECM-receptor interaction	6	2.5E-3	SDC1, LAMA3, CD36, TNXB, COL6A6, ITGA8
Hsa04974	Protein digestion and absorption	5	1.5E-2	COL6A6, COL13A1, MME, SLC1A1, COL10A1
hsa04360	Axon guidance	5	4.9E-2	SEMA5A, LIMK1, CFL2, ABLIM3, SLIT2

ECM – extracellular membrane; KEGG – Kyoto Encyclopedia of Genes and Genomes. When more than 3 enriched terms were identified in each category, the top 3 terms were selected according to P value.

PPI Network Construction and Module Analysis

The network of protein interactions was obtained using STRING and then the top 10 hub genes with high degrees of expression were identified using Cytoscape (Figure 2). Those hub genes were CAV1, TEK, SLIT2, RHOJ, DGSX, FOXF1, HLF, MEIS1, PTPRD, and ADRB2, all of which were downregulated in LUAD. The degree centrality information for the hub genes is shown in Figure 3 and Table 3. Results of GO and KEGG enrichment analysis of the hub genes are shown in Figure 4. GO and KEGG analysis is one of the most important parts of bioinformatics. It can help us systematically understand the ways in which DEGs, including hub genes, affect biological functions. In the present study, the results of GO and KEGG analysis revealed the function of DEGs in LUAD, as shown in Table 2.

Figure 2

Protein-protein interaction network of the differentially expressed genes (DEGs) constructed with Cytoscape 3.7.2. DEGs shown in purple are upregulated and those in green are downregulated. The 10 hub genes ranked by Degree Centrality, which are located in the middle, are CAV1, TEK, SLIT2, RHOJ, DGSX, FOXF1, HLF, MEIS1, PTPRD, and ADRB2.

Figure 3

Information about the degree, betweenness, and closeness of hub genes. The size of the bubbles represents betweenness. The depth of color indicates closeness. The length of the column represents degree. As depicted in the figure, CAV1 has the greatest degree of expression, betweenness, and closeness.

Table 3

Information for the top 10 genes with high degrees of expression.

Gene symbol	Gene name	Degree of expression
CAV1	Caveolin 1	42
TEK	TEK receptor tyrosine kinase	34
SLIT2	Slit guidance ligand 2	31
RHOJ	Ras homolog family member J	30
DGSX	Glypican 3	27
FOXF1	Forkhead box F1	25
HLF	HIF-1-alpha-like factor	25
MEIS1	Meis homeobox 1	25
PTPRD	Protein tyrosine phosphatase receptor Type D	25
ADRB2	Adrenoceptor beta 2	24

Figure 4

Enrichment results for the hub genes based on assessment of the top Gene Ontology and Kyoto Encyclopedia of Genes and Genomes terms. GO: 0001763, morphogenesis of a branching structure; GO: 0001525, angiogenesis; GO: 0005515, protein binding; GO: 0001077, transcriptional activator activity, RNA polymerase II core promoter proximal region sequence-specific binding; GO: 0005887, integral component of plasma membrane; GO: 45121, membrane raft; GO: 0005886, plasma membrane; hsa: 04360, axon guidance.

Kaplan-Meier Survival Analysis of hub Genes

The relationship between hub genes and prognostic conditions was evaluated with the Kaplan-Meier plotter online analysis tool. CAV1, TEK, SLIT2, FOXF1, HLF, MEIS1, PTPRD, and ADRB2 were found to be associated with favorable overall survival (OS) in patients with LUAD (Figure 5). The probes for hub genes that we used are shown in Table 4. The schema for this study is shown in Figure 6.

Figure 5

Results of Kaplan-Meier overall analyses for hub genes identified in the study. CAV1, TEK, SLIT2, FOXF1, HLF, MEIS1, PTPRD, and ADRB2 were found to be associated with favorable overall survival in patients with LUAD.

Table 4

Hub genes in the Kaplan-Meier plotter database and corresponding probes.

Gene symbol	Probe ID
CAV1	212097_at
TEK	206702_at
SLIT2	228850_s_at
RHOJ	235131_at
DGSX	209220_at
HLF	204754_at
MEIS1	204069_at
PTPRD	214043_at
FOXF1	205935_at
ADRB2	206170_at

Figure 6

Schema of the study.

Inhibition of HLF Promotes A549 Cell Proliferation

To determine the expression of HLF with the lowest log rank P value, qRT-PCR analysis was used to compare the HLF expression between LUAD and adjacent normal tissue. A significantly lower level of HLF mRNA was found in the LUAD (Supplementary Figure 1A). Next, to study the function of HLF in LUAD, HLF siRNA was used to knock down its expression in A549 cells. The qRT-PCR result indicated an acceptable knockdown efficiency (Supplementary Figure 1B). In addition, to explore the impact of HLF on the proliferation of A549 cells, a qRT-PCR assay was performed to determine the expression of Cyclin D1 and Cyclin D3. Compared with the siRNA-NC group, the expression of Cyclin D1/3 was significantly increased (Supplementary Figure 1C).

Discussion

LUAD is a malignant tumor that originates from bronchial mucosal glandular epithelium and which accounts for about 33% of all lung carcinomas [14]. Lung cancer may not cause any symptoms, especially in the early stages of the disease. Patients with LUAD may develop a cough and chest symptoms. In the late stage, the tumor can block off an airway and interfere with breathing. Current treatments for LUAD include targeted therapy, radiotherapy, immunotherapy, chemotherapy, and surgery, but there is no criterion standard for the disease. The mechanism of LUAD is still unclear and there is controversy about what causes the disease. DEG analysis based on bioinformatic methods can help us explore the pathogenesis of LUAD and find novel targeted treatments for it.

CAV1 is an isoform in the caveolin family that participates in some cellular processes, such as cellular signaling and immune responses [15]. In fact, CAV1 has been proven to be involved in some cancers, but there is controversy about the role it plays [16]. It is upregulated in some cancers, such as those of the breast, pancreas, and prostate, but downregulated in ovarian and colorectal cancers [17-20]. In the present study, CAV1 was found to be downregulated in LUAD, and those data combined with the Kaplan-Meier survival analysis results indicate that it can inhibit LUAD. This result is consistent with some previous studies, which have shown that CAV1 can positively regulate lung cancer by promoting cell proliferation and decreasing cell apoptosis [15].

TEK reportedly may have prognostic potential in head and neck SCCs treated with radiation or chemoradiation [21]. TEK was also proven to be associated with gliomas, other brain tumors, and breast cancer [22,23]. Some preliminary data show that TEK suppresses LUAD cell phenotypes by interacting with miR-19a-3p [24]. Taking those findings together with the results from our study, it can be concluded that TEK plays a positive role in the development of LUAD.

As reported in some studies, SLIT2 may suppress lung cancer progression [25,26]. The underlying mechanism is unclear. We found that expression of SLIT2 is downregulated in LUAD and survival analysis shows that it is associated with favorable prognosis. Therefore, SLIT2 may have a tumor suppressor role in LUAD.

The FOXF1 transcription factor regulates E-cadherin expression downstream of the P53 family members, modulating cancer cell migration and invasion [27]. In NSCLC, FOXF1 may inhibit cancer growth by inducing tumor suppressor and G1-phase cell-cycle arrest [28]. In LUAD in our study, downregulation of FOXF1 and positive survival results were observed. Therefore, we concluded that FOXF1 is a prognostic factor for LUAD [29].

HLF is an alias for EPAS1. Low EPAS1 expression in cancer, including LUAD, has been reported [30]. In NSCLC, EPAS1 can promote peritoneal carcinomatosis by enhancing the mesothelial-mesenchymal transition [31]. In the present study, downregulation of EPAS1 was further confirmed and EPAS1 was associated with favorable prognosis for LUAD, indicating that it may be a promising biomarker and therapeutic target.

MEIS1 may play a role in limiting the proliferation of NSCLC cells [32]. In fact, MEIS1 may mediate suppression of metastasis in many types of cancer, including LUAD [33]. Low expression of MEIS1 in LUAD was consistently observed in the present study, and its favorable effect on LUAD was verified by survival analysis.

PTPRD was been proven to be a potential candidate tumor suppressor gene not only in lung cancer but also in some other cancers [34-36]. Deletion of PTPRD and CDKN2A cooperate to accelerate tumorigenesis [37]. Indeed, in LUAD, we found that expression of PTPRD is downregulated and high expression of it led to better survival results.

The literature regarding the relationship between ADRB2 and LUAD is limited. ADRB2 reportedly plays an essential role in gastric cancer, promoting progression and metastasis [38]. However, a case-control study did not support a major independent role for ADRB2 polymorphisms in LUAD risk [39]. In our study, ADRB2 was downregulated in LUAD, and high expression of it was associated with longer survival time.

In the present study, CAV1, TEK, SLIT2, FOXF1, HLF, MEIS1, PTPRD, and ADRB2 were found to be associated with favorable OS in patients with LUAD. Moreover, inhibition of HLF was found to promote A549 cell proliferation, indicating that these genes are involved in the pathological development of LUAD and that they can be used as predictors of prognosis in patients.

There are some limitations of the present study. First, there was no in vitro or in vivo validation. Second, only 3 datasets were used. More data could provide more convincing evidence about our results. Nevertheless, we believe that this bioinformatics-based study provides some useful information for further research on LUAD.

Conclusions

In this bioinformatics-based study, CAV1, TEK, SLIT2, HLF, MEIS1, PTPRD, FOXF1, and ADRB2 were identified as hub genes in the DEG interaction network for LUAD. They were involved in pathological development of LUAD and associated with favorable OS. These genes can be used as predictors of prognosis in patients with LUAD. Furthermore, these findings can facilitate further investigation into the mechanism of LUAD.

Inhibition of HLF promotes A549 cell proliferation. (A) qRT-PCR results of the HLF expression in LUAD tissue and the adjacent normal tissue. (B) qRT-PCR results of the HLF expression in A549 cells after siRNA-NC and siRNA-HLF treatments. (C) The expression of proliferation-related genes Cyclin D1 and Cyclin D3 in the different groups was assessed by qRT-PCR analysis. * p<0.05, ** p<0.01, *** p<0.001.