Alternatively Expressed Transcripts Analysis of Non-Small Cell Lung Cancer Cells under Different Hypoxic Microenvironment.

Globally, non-small cell lung cancer (NSCLC) is the most fatal form of malignancy. Numerous studies have shown that people living at high altitudes are at a higher risk for cancer. Hypoxia is one of the most important features in high altitude area. Compared with normal cells, cancer cells are more adapted to hypoxia atmosphere. However, at high altitudes, hypoxic conditions are also accompanied by other altered environmental conditions. To identify the single influence of hypoxia, we performed second-generation sequencing to identify gene expression changes triggered by the different oxygen concentrations. We identified 782 genes in A549 cells and 1122 genes in H520 cells that showed altered expression by the combined analysis in 5% oxygen concentration group and 1% oxygen concentration group control group. We further analyzed these targets and found 113 genes altered in both cell lines. Interestingly, we found KxD1 was the only one in both top 10 lists. Further analysis revealed KxD1 to be significantly elevated in NSCLC patients and negatively correlated with prognosis in stage I and II NSCLC patients. Moreover, this correlation reversed in stage III patients. Additionally, compared with patients who only received clean margin operation or chemotherapy, patients who received radiotherapy also showed opposite result. Thus, KxD1 may be a promising target for the treatment of NSCLC in high-altitude areas.

1. Introduction

Lung cancer is one of the most deleterious forms of cancers worldwide, contributing to approximately 25.3% of all cancer-related deaths [1]. In 2020, the number of new cases of lung and bronchial cancer in the United States had risen to 228820, and the estimated death toll is 135720 [2]. Non-small cell lung cancer (NSCLC) accounts for approximately 85% of all lung cancer cases, of which lung adenocarcinoma (LUAD) and lung squamous cell carcinoma (LUSC) account for the largest proportion of NSCLC [3, 4]. At present, although epidemiological studies have confirmed that environmental pollution and smoking are closely related to the occurrence and development of NSCLC, further investigations are needed to identify additional factors [5].

The concept of tumor hypoxia was first presented in 1955, while studying tumor specimens obtained from lung cancer patients [6]. Over the next 60 years, scientists have gradually confirmed that hypoxia is a widespread characteristic in tumors and is closely related to tumor differentiation, angiogenesis, energy metabolism, invasion, and drug resistance [7-11]. Compared with normal cells, cancer cells are more adapted to grow in anoxic environments [12]. A 100 m increase in altitude leads to an approximately 1.2 mmHg decrease in the partial pressure of oxygen. Therefore, the oxygen content of plateaus is approximately 20%–40% lower than that at sea level [13, 14], and numerous studies have found that the incidence and mortality rates of cancer are higher at high altitudes than those at sea level [15, 16].

At present, the effect of hypoxia on gene expression in lung cancer is still lacking; second generation sequencing has become a common experimental technology, which can provide us with detailed, accurate, and specific information pertaining to gene expression [17, 18]. Thus, in this study, we conducted high-throughput transcriptome sequencing of NSCLC cell lines, A549 and H520, under hypoxic conditions of different oxygen concentrations to identify consequent gene expression changes, screen new coding genes and transcription factors, and perform a detailed survival analysis of NSLSC cells. We provide evidence for potential targets that can be used in future clinical treatment of lung cancer patients living at high altitudes.

2. Materials and Methods

2.1. Reagents

Dulbecco's Modified Eagle's medium (DMEM), fetal bovine serum (FBS), and penicillin-streptomycin dual antibiotic solution were purchased from Gibco (USA). The Transwell system and matrix adhesive were purchased from Corning (USA). TB Green Premix Ex Taq II and PrimeScript RT Master Mix were purchased from Takara (Japan). TRIzol reagent was purchased from Invitrogen (USA). Trypsin/EDTA solution was purchased from Thermo Fisher (USA).

2.2. Cell Culture

A549 and H520 cells were purchased from the Chinese Academy of Sciences Cell Bank (Shanghai, China). H520 and A549 cells were cultured in DMEM supplemented with 10% FBS, 100 μg/mL streptomycin, and 100 U/mL penicillin. The cells were seeded in six-well culture plates in the medium for 24 h. Each well was designated randomly into three groups: normal (21%) oxygen group (N), 5% oxygen group (L), and 1% oxygen (S). Mixed gas (CO2 + 5% Bal N2) was infused into the closed cell incubator (MIC-101, Hangzhou, China) for 24 h to regulate oxygen concentration.

2.3. Invasion Assay

Cell invasion ability was determined by a matrix gel invasion assay using a Transwell system. The filter surface of the upper chamber (8 μm aperture) was coated with a 1 mg/mL of matrix, while the lower compartment contained complete medium. After incubation for 24 h, the upper noninvasive cells were removed using a cotton swab. Migrated cells were fixed with 4% paraformaldehyde and stained with crystal violet. Cell infiltration was observed and imaged under a microscope (CKX41, Olympus, Japan), and five random field counts on the submembrane surface were used to quantify cell invasion.

2.4. Wound Healing Assay

A549 and H520 cells in the logarithmic phase of growth were digested with trypsin and cultured in six-well plates with DMEM. The cells were cultured for 12 h after they attached to the walls of the plate. A 1 ml plastic pipette tip was used to scratch vertical lines in the middle of the well. After washing with phosphate-buffered saline at 37°C, fresh 2% FBS prepared in DMEM was added. The migration distance from the edge to the center of the scratch for each group was observed at 0 h, 12 h, and 24 h under an inverted relative ratio microscope. The scratch area was measured using the ImageJ software (National Institutes of Health).

2.5. Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR)

Total RNA from the culture cells was extracted using TRIzol reagent, and its concentration was measured using a Nanodrop (Thermo Fisher Scientific, USA). Total RNA was reverse-transcribed using a cDNA synthesis kit. qRT-PCR was performed using the SYBR Green PCR Master Mix in the SteponePlus RT-PCR system (Thermo Fisher Scientific). The 2 ΔΔ Ct method was used to analyze relative gene expression. The primer sequences used are as follows: induced VEGFR Forward (Fw): CAGAATCATCACGAAGTGGTG, Reverse (Rev): GAAGATGTCCACCAGG GTC; FN1 Fw: AAATGGCCAGATGATGAGC, Rev: TAACACGTTGCCTCATGA G; GLUT1 Fw: GAAGGTCACCATCCTGGAG, Rev:AATAGAAGACAGCGTTG ATGC; and GAPDH Fw: TCAAGATCATCAGCAATGCC, Rev: CGATACCAAA GTTGTCATGGA.

2.6. Preparation of Transcriptome Library and Sequencing

To determine the integrity of the RNA and the extent of contamination by DNA, the samples were subjected to agarose gel electrophoresis. Then, a nanophotometer was used to estimate the purity of the samples, wherein OD 260/280 and OD 260/230 ratios were determined, and a Agilent 2100 BioAnalyzer was used to determine the RNA integrity, after which mRNAs with polyA tail were enriched by Oligo(dT) magnetic beads. After end repair, a tail was added, the sequencing connector was connected, and cDNA of about 250–300 bp was screened using AMPure XP beads. PCR amplification was performed, and the AMPure XP beads were used to purify the PCR products and obtain a library. qRT-PCR was used to accurately quantify the library effective concentration (library effective concentration was higher than 2 nM) and ensure library quality. The constructed libraries were sequenced with Illumina X.

2.7. Bioinformatics Analysis

All the analyses were performed on high-quality clean data. Hisat2v2.0.5 was used to construct the reference genome index, and Stringtie (1.3.3b) (Mihaela Pertea et al. 2015) was used to predict new genes. Featurecounts (1.5.0-P3) calculated the readings mapped to each gene and performed differential expression analysis (two biological replicates per group) between the groups using Deseq2 software (1.16.1). The ClusterProfiler (3.4.4) software was used to perform gene ontology (GO) enrichment analysis of differentially expressed genes and statistical enrichment of differentially expressed genes in the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway.

2.8. Validation of Differentially Expressed Genes

The identified differentially expressed genes that were associated with NSCLC in a hypoxic microenvironment were validated using datasets expressed in The Cancer Genome Atlas (TCGA) and related databases. TNMplot (https://www.tnmplot.com/) was used to analyze differential gene and transcript expression in tumor and normal tissues. The Human Protein Atlas [19] (https://www.proteinatlas.org/) provided the localization and qualitative data for the KXD1 protein.

2.9. Survival Analysis

The Oncomine database (https://www.oncomine.org/) [11] and the Kaplan–Meier Plotter (https://kmplot.com/analysis/) [20] was used to determine the effect of differentially expressed genes (DEGs) under different concentrations of oxygen on the overall survival of NSCLC patients. The survival curve was downloaded from the following data links: TCGA CAARRAY GSE31210, GSE14814, GSE19188, GSE29013, GSE30219, GSE3141, GSE31908, GSE37745, GSE43580, GSE4573, GSE50081, GSE31908, and GSE8894. The influence of these genes on tumor prognosis was analyzed, and differences were considered statistically significant at P < 0.05.

2.10. RNA-Sequencing (RNA-Seq) Datasets

The original RNA-seq data have been deposited to NCBI sequence read archive (PRJNA656891). The code required for the replication of differential expression and differential splicing analysis will be published when this manuscript is accepted.

2.11. Statistical Analysis

The SPSS 22.0 software was used for statistical analysis and GraphPad Prism 7.0 was used for graphical representation of the data. The two-tailed Student's t-test was used for two-group comparisons, while the one-way analysis of variance was used for multigroup comparisons. Differences were considered statistically significant at P < 0.05.

3. Results

3.1. Hypoxia Significantly Increased the Migration and Invasion of NSCLC Cells In Vitro

Wound healing assay revealed that the migration rate of A549 and H520 cells increased gradually with the aggravation of hypoxia in group N (21% O2) (Figure 1(a)), while invasion assay revealed that cells of group S (5% O2) and L (1% O2) showed a significantly higher degree of invasion than that observed in group N (Figure 1(b)).

Figure 1

Hypoxia promotes metastasis and invasion of NSCLC. (a) Cellular wound healing assays of tumor cells (100x), N 21% O2, L 5% O2, S 1% O2; the results of wound healing assays. (b) Cell invasion was measured by Transwell assays (100x). The numbers of invading cells was counted from five fields of view in each group. (c) The mRNA level of VEGFR, FN1, and GLUT1 was detected by qRT-PCR.ns. P > 0.05, P < 0.05, P < 0.01, and  P < 0.001.

qRT-PCR was used to estimate the levels of mRNA expression of VEGFA, FN1, and GLUT1 (Figure 1(c)), which are all genes regulating the migration and invasion ability of lung cancer cells, and we found that their expression levels were significantly higher under hypoxic conditions than under control conditions. These results showed that hypoxia significantly increased the migration and invasion of both LUAD and LUSC cells in vitro.

3.2. Hypoxia Promotes a More Abundant Set of Differentially Expressed Transcripts (DETs) than DEGs

We analyzed the gene and transcript expression levels of A549 and H520 cells at different oxygen concentrations and found that the correlation between the transcriptional expression levels of samples from each group of the two lung cancer cell lines was greater than 0.9, indicating a close correlation between the groups (Figure 2(a)).The transcription discrepancy between the RNA-seq data of both A549 and H520 is more obvious with the severity of hypoxia, wherein the transcription expression is higher than the gene expression (Figure 2(b)). Combining the analysis of DEGs and DETs of the three experimental groups of both cell lines, the corresponding heatmaps validating the aforementioned results are shown in Figure 2(c). A comparison between LUSC (H520) cells grown under 21% oxygen and 5% oxygen revealed 5,244 upregulated transcripts and 2,801 downregulated transcripts, while a comparison between LUSC cells grown under 21% oxygen and 1% oxygen revealed 7,284 upregulated transcripts and 3160 downregulated transcripts. Comparing LUAD (A549) cells grown under 21% and 5% oxygen revealed that 7,249 transcripts were upregulated and 3,096 transcripts were downregulated, while comparing LUAD cells grown under 21% and 1% oxygen revealed that 8,957 transcripts were upregulated and 4,082 transcripts were downregulated (Figure 2(d)). These results suggest that hypoxia mainly exerts its effect at the posttranscriptional level rather than the transcriptional level.

Figure 2

The results of NSCLC RNA-seq. (a) Heat map of correlation between samples. (b) Analysis and statistics of gene and Transcription expression difference. The horizontal axis is the name of difference comparison, and the vertical axis is the number of up-down difference genes. Among them, the green is downregulated counts and the red is upregulated counts. (c) Heatmaps of differential expressed genes (DEGs) and differential expressed transcripts (DETs) in NSCLC cells under different oxygen concentrations. On the left is the tree diagram of gene clustering. (d) Transcripts volcano maps of H520 and A549 cells compared with 21% O2VS 5% O2 and 21% O2VS 1% O2.

3.3. Functional GO Enrichment and KEGG Analysis of Hypoxia-Induced Transcriptome

We found 782 coexpressed DETs in A549, 1122 coexpressed DETs in H520, and 117 coexpressed transcripts in both A549 and H520 cells, under different oxygen concentrations (Figure 3(a)). Using the league Sichuan Biological Cloud Platform [17] (https://www.omicstudio.cn/index), we performed an enrichment analysis for KEGG pathways and GO terms of the DETs. Enriched GO terms upon functional analysis of the DETs were namely “cells”, “cell structure”, “intracellular processes”, “cellular processes”, and “binding” (Figure 3(b)). Enriched KEGG pathways included focal adhesion, HIF-1, AMPK, and mTOR in A549 cells, while the PI3K-Akt signaling pathway, MAPK signaling pathway, and regulation of the actin cytoskeleton were enriched in H520 cells (Figure 3(c)).

Figure 3

Functional enrichment and KEGG analysis of DETS. (a) Venn diagram of differential transcripts. (b) Scatter plot of top 20 enriched GO terms of molecular function (MF), biological process (BP), and cellular component (CC) separately. (c) Scatter plot of top 30 enriched KEGG pathways.

3.4. KXD1 Was the Top-Most DET in Both Cell Lines

DETs common to both LUAD and LUSC cells were selected, amounting to a total of 113 transcripts (Figure 4(a)), of which KxD1 was the top most gene that was chosen from our analysis of hypoxia-induced differential expression of the protein network (Figure 4(c)). In the protein-protein interaction network, only RPL34 was found to interact with KXD1, indicating that functional studies of KxD1 are still lacking (Figure 4(b)).

Figure 4

KXD1 is the most valuable protein according to DETs analysis. (a)Venn diagram of the most significant difference transcripts between H520 and A549. (b) PPI networks among the target genes of 113. (c) Venn diagram of top 10 DETs in A549 and H520 induced by hypoxia (N/L, 21% O2/5% O2; N/S, 21% O2/1% O2).

3.5. KXD1 Was Overexpressed in Both LUSC and LUAD Cells

According to the immunohistochemistry results downloaded from The Human Protein Atlas, the expression of KXD1 was significantly upregulated in both LUSC and LUAD cells compared to that in control cells (Figure 5(a)). A comparison of 486 LUAD patients and 524 normal patients showed that KXD1 expression was elevated in adenocarcinoma, and similarly, a comparison of 476 LUAD patients and 501 normal patients also showed that KXD1 was elevated in LUAD patients (Figure 5(b)). We conclude that KXD1 is closely associated with NSCLC and has the potential to be a primary screening and prognostic marker. We also speculate that elevated KXD1 caused by hypoxia is key to the increasing NSLSC migration and invasion.

Figure 5

KXD1 was overexpressed both in LUSC and LUAD. (a) IHC stained tissues representing KXD1 expression is shown in NSCLC or normal lung tissues, respectively. (b) Expression level of KXD1 in normal lung tissue and NSCLC tissue.

3.6. Survival Analysis of KxD1

According to the Kaplan–Meier analysis, we found that overexpression of KxD1 was significantly correlated with poorer OS in NSLSC patients (Figure 6(a)). We also found that this trend is consistent in stage I and stage II NSLSC; however, this trend reverses for stage III NSLSC patients (the IV stage survival analysis was excluded because only 4 patients were included) (Figure 6(b)). Furthermore, we also found that, compared with patients who received surgery or chemotherapy, the OS was negatively correlated with KxD1 expression levels, but patients who received radiotherapy showed opposite trend (Figure 6(c)).

Figure 6

Survival analysis of KxD1. (a) Survival analysis of KxD1 in the NSCLC population. (b) Survival analysis of KxD1 in different stages of the NSCLC population. (c) Survival analysis of KxD1 in NSCLC patients treated with different treatments.

4. Discussion

Physiologically appropriate functioning of cellular process depends on adequate oxygen and energy supply [21]. Hypoxia has been reported to promote metastasis and invasion of NSCLC [22-24]. In this study, we found that hypoxia significantly increased the migration and invasion ability of NSCLC cells. Moreover, these changes were more obvious in the 1% oxygen concentration group than in the 5% oxygen concentration group in both LUSC and LUAD cells. This may be because, under hypoxic conditions, the transcriptional instability of tumor cells may cause the activation of some cancer survival-related factors, resulting in the enhancement of tumor migration and invasion and promotion of cancer development [25, 26].

With the increase of altitude, it is often accompanied by the changes of inhalable particles, sunlight exposure, air pressure, and other factors [27, 28]. This experiment aims to explore the single effect of hypoxia on the incidence of lung cancer in high altitude areas; therefore, we used the in vitro model of hypoxia to analyze the changes of gene expression.

In this study, we found a total of 113 transcripts that were expressed differently in LUSC and LUAD cells with increasing hypoxic conditions. Among them, KxD1 was the top DET among both cell lines, and it may be the key molecule induced by hypoxia in NSCLC. KxD1 is a BLOS1-interacting protein that is correlated with the biogenesis of lysosome-related organelles [29]. In KxD1 knockout mice, BLOS1 expression has been reported to decrease [29]. At present, the pathophysiological function of KxD1 remains unknown; however, an increase BLOS1 degradation can lead to elevated endoplasmic reticulum stress and accumulation of ubiquitinated protein aggregates, significantly influencing cellular stress [30]. Therefore, we speculated that KxD1 may also be a cellular stress modulator.

The correlation between KxD1 and NSCLC has not been previously reported. In this study, we found that KxD1 was markedly elevated in NSCLC and negatively correlated with the prognosis of patients with NSCLC. However, we also found that although in stage I and stage II the prognosis trend was consistent with the overall population, this trend was reversed in stage III. Furthermore, compared with patients who received surgery or chemotherapy, patients who received radiotherapy also showed the opposite trend. Our results indicate that KxD1 is not only positively correlated with cancer invasion and metastasis, but also negatively correlated to nonoperational anticancer therapy resistance. Hence, elevated KxD1 levels may be an indicator of radiotherapy sensitivity.

5. Conclusion

We have demonstrated the influence of different hypoxic concentrations on gene expression in NSCLC. In this study, we also found a promising molecule, KxD1, which may be beneficial for the treatment of NSCLC in high altitude areas. However, the exact function and the correlation with radiotherapy of KxD1 require further research.