This article is mentioned in:

Introduction

Lung cancer is among the most common types of cancer, accounting for ~13% of all cancer cases (1). The generally poor prognosis of lung cancer renders it a leading cause of cancer-associated mortality worldwide (2). In 2010, 1.5 million mortalities due to lung cancer were reported, representing 19% of all cancer-associated mortality (3). The incidence of lung cancer has doubled in China over the past decade due to issues including the aging population, smoking and the reduced air quality (4). Lung cancer is initiated by the activation of oncogenes or the inactivation of tumor suppressor genes (5). Despite advances in diagnosis and treatment, the prognosis of lung cancer remains relatively poor. The identification of reliable biomarkers and novel genes involved in lung cancer carcinogenesis is important for improving the ability to predict the prognosis and to guide the therapy of lung cancer.

Forkhead box D3 (FOXD3) is a member of the FOX transcription factor family, which is characterized by a distinct forkhead domain (6). FOXD3 acts as a transcriptional repressor or activator (7). The abnormal expression of FOXD3 has been reported to participate in tumor onset and progression in non-small cell lung cancer tumor cells (8). Other studies have indicated tumor suppressive activities for FOXD3, including the inhibition of cell growth and invasion in various types of cancer, including gastric cancer and melanoma (9,10). A number of genes associated with tumorigenesis have been reported to be targets of FOXD3. One study demonstrated that FOXD3 regulated RND3 expression and migration properties in melanoma cells (11). Another reported that FOXD3 exhibited tumor suppressive activity that affected the growth, aggressiveness and angiogenesis of neuroblastoma through the transcriptional regulation of NDRG1 (12). However, the role of FOXD3 in lung cancer remains uncharacterized.

In this study, differentially expressed genes (DEGs) and alternative splicing genes (ASGs) were identified in FOXD3-knockout samples compared with normal samples. Functional and pathway enrichment analyses of the DEGs and ASGs were performed. A protein-protein interaction (PPI) network was constructed based on the overlaps between the DEGs and ASGs. An improved understanding of FOXD3 in regulating the process of lung cancer was obtained, which may allow the development of novel strategies for the diagnosis and therapy of lung cancer.

Materials and methods

Datasets

The gene expression profile GSE64513 was downloaded from the Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/) database. The data set contained the RNAseq data from 6 samples, including 3 FOXD3-knockout A549 lung cancer cell samples and 3 normal A549 cell samples.

Screening of DEGs and ASGs

The data was first analyzed using FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc), a java-based high-throughput data quality control software. Reads with base quality scores <20 were discarded, and reads longer than 30 bp were selected for further investigation. The remaining reads were mapped to the GRCh37/hg19 genome based on the Tophat2 program (13). The number of reads mapped to the exons of each gene was counted with the HTSeq-Count tool (14) and regarded as the expression profile of each gene. Differently expressed genes (DEGs) in FOXD3 knockout lung cancer samples compared with normal samples were identified using the edge R package (15) with the following thresholds: False discovery rate <0.01 and |log (fold change)| >1. The hierarchical clustering of DEGs was performed using the heatmap.2 function of the gplots package in Various R Programming tool (version 2.12) (16). The alternative splicing genes (ASGs) in the FOXD3 knockout samples were identified using the replicate multivariate analysis of transcript splicing (rMATS) program, a computer program designed to detect differential alternative splicing from replicate RNA-Seq data (17).

Functional and pathway enrichment analysis

The Database for Annotation, Visualization and Integrated Discovery (DAVID; https://david.ncifcrf.gov/) is a web-based tool for genomic functional annotations (18). To further explore the biological functions of the DEGs and ASGs, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed using DAVID, with the threshold of P<0.05.

Construction of a PPI network

The overlapping DEGs and ASGs were analyzed using the Search Tool for the Retrieval of Interacting Genes (STRING; http://string-db.org/) (19,20). A PPI network to illustrate the identified interactions was constructed and visualized using Cytoscape 3.4 (21).

Results

Identification of DEGs and ASGs

The total number of reads, the number of mapped reads and the mapping rate of each sample is provided in Table I. A total of 1,853 DEGs were identified, of which 382 were upregulated and 1,471 were downregulated. The top 20 DEGs are listed in Table II. Fig. 1 demonstrates the hierarchical clustering results for each sample graphically (Fig. 1A), the fold-change trend of the expression of the identified DEGs (Fig. 1B) and the hierarchical cluster analysis of the samples based on the DEGs (Fig. 1C). A total of 2,249 genes with alternative splicing were identified in FOXD3-knockout lung cancer samples compared with normal A549 cell samples, including 545 with an alternative 3′ splice site, 412 with an alternative 5′ splice site, 1,629 with mutually exclusive exons and 67 with retained introns.

Figure 1. DEGs between 3 expression profiles from A549 cells with FOXD3 knockout and 3 from normal A549 cells. (A) Hierarchical clustering analysis of the 6 samples by MDS; (B) the fold change trend of the DEGs; and (C) the hierarchical clustering analysis of the samples based on the expression of the DEGs. DEGs, differentially expressed genes; FOXD3, forkhead box D3; MDS, multidimensional scaling; DEG, differentially expressed gene; FC, fold change.

Table I. Total reads, the number of mapped reads, and the mapping rates of each sample.

Table I.

Total reads, the number of mapped reads, and the mapping rates of each sample.

Sample	Total reads	Mapped reads	Mapping rate, %
SRR1734826	10,339,232	9,252,784	89.5
SRR1734827	10,472,212	9,298,751	88.8
SRR1734828	10,868,010	9,651,798	88.8
SRR1734829	11,224,483	9,666,662	86.1
SRR1734830	10,548,877	9,241,415	87.6
SRR1734831	11,578,464	10,104,885	87.3

Table II. Top 20 differentially expressed genes of the forkhead Box D3-knockout lung cancer A549 cell samples compared with normal A549 cells.

Table II.

Top 20 differentially expressed genes of the forkhead Box D3-knockout lung cancer A549 cell samples compared with normal A549 cells.

Gene symbol	False discovery rate	P-value	Log fold change
SAA1	1.02×10−426	8.46×10−313	−7.70213
C3	3.82×10−478	8.15×10−562	−2.97374
GAS6	9.26×10−283	2.01×10−286	−2.52365
CFB	4.97×10−271	1.44×10−274	−3.25002
LCN2	7.38×10−254	2.67×10−257	−3.10027
TGM2	8.36×10−244	3.63×10−247	−2.45302
SAT1	3.53×10−238	1.79×10−241	−2.46616
PDZK1IP1	3.63×10−230	2.11×10−233	−4.94302
PLAU	7.60×10−226	4.96×10−229	−2.07647
SAA2	6.42×10−194	4.65×10−197	−7.47593
TNIP1	2.23×10−189	1.78×10−192	−2.05325
SPP1	2.23×10−189	1.94×10−192	2.29076
S100A8	1.27×10−168	1.20×10−171	−3.79768
TMEM132A	1.42×10−166	1.44×10−169	−1.93297
ASNS	6.48×10−156	7.04×10−159	−2.09409
SERPINE1	2.92×10−149	3.38×10−152	−2.55607
PHLDB2	4.19×10−147	5.17×10−150	−1.54635
ICAM1	7.56×10−144	9.87×10−147	−2.03788
CXCL8	1.08×10−138	1.48×10−141	−5.25318

Enriched GO terms and KEGG pathways of DEGs and ASGs

The DEGs were enriched in 338 GO terms and 21 KEGG pathways. The ASGs were enriched in 470 GO terms and 22 KEGG pathways. The top 10 GO terms for the ASGs and DEGs are listed in Fig. 2A and B, respectively. Table III lists the enriched KEGG pathways for the ASGs and DEGs. The DEGs were predominately enriched in ‘graft-vs.-host disease’, ‘hematopoietic cell lineage’, ‘ECM-receptor interaction’ and ‘NOD-like receptor signaling pathway’. The ASGs were predominately enriched in ‘ubiquitin mediated proteolysis’, ‘chronic myeloid leukemia’, ‘aminoacyl-tRNA biosynthesis’ and ‘mTOR signaling pathway’.

Figure 2. The top 10 GO terms in which (A) the DEGs, and (B) the alternative splicing genes were enriched. GO, Gene Ontology; DEGs, differentially expressed genes.

Table III. Continued.

Table III.

Continued.

A, Enriched KEGG pathways for DEGs
Pathway name	Genes, n	P-value
Graft-vs.-host disease	12	0.0006
Hematopoietic cell lineage	18	0.0019
NOD-like receptor signaling pathway	14	0.0038
ECM-receptor interaction	17	0.0038
Cell adhesion molecules	23	0.0044
Allograft rejection	10	0.0045
Cytokine-cytokine receptor interaction	38	0.0055
Glycine, serine and threonine metabolism	9	0.0059
MAPK signaling pathway	38	0.0075
p53 signaling pathway	14	0.0085
Natural killer cell mediated cytotoxicity	22	0.0100
Toll-like receptor signaling pathway	18	0.0106
Viral myocarditis	14	0.0122
Nitrogen metabolism	7	0.0157
Arginine and proline metabolism	11	0.0217
Complement and coagulation cascades	13	0.0230
Pathways in cancer	42	0.0266
Axon guidance	20	0.0274
Pathogenic Escherichia coli infection	11	0.0344
B cell receptor signaling pathway	13	0.0413
Small cell lung cancer	14	0.0438
B, Enriched KEGG pathways for ASGs
Pathway name	Gene, n	P-value
Ubiquitin mediated proteolysis	33	0.0002
Chronic myeloid leukemia	20	0.0014
Aminoacyl-tRNA biosynthesis	13	0.0029
mTOR signaling pathway	15	0.0032
Renal cell carcinoma	18	0.0040
Pancreatic cancer	18	0.0054
Neurotrophin signaling pathway	26	0.0079
Pathways in cancer	55	0.0133
Ribosome	19	0.0173
Pyrimidine metabolism	20	0.0208
Acute myeloid leukemia	14	0.0215
Small cell lung cancer	18	0.0250
Wnt signaling pathway	28	0.0276
Cell cycle	24	0.0297
VEGF signaling pathway	16	0.0373
Lysine degradation	11	0.0377
Insulin signaling pathway	25	0.0387
Glioma	14	0.0403
Prostate cancer	18	0.0414
Lysosome	22	0.0465
Endometrial cancer	12	0.0484
N-glycan biosynthesis	11	0.0495

[i] KEGG, Kyoto Encyclopedia of Genes and Genomes; DEG, differentially expressed genes; ASG, alternative splicing gene; NOD, nucleotide oligomerization; ECM, extracellular matrix; MAPK, mitogen-activated protein kinase; mTOR, mechanistic target of rapamycin; tRNA, transfer RNA; VEGF, vascular endothelial proliferation factor.

PPI network

A total of 199 overlaps between the DEGs and the ASGs were identified, and the PPI network constructed from the 199 overlapping genes contained 97 nodes and 115 pairs (Fig. 3). Table IV lists the top 20 pairs with highest combined scores, and Table V lists the top 10 nodes according to connectivity degree.

Figure 3. PPI network demonstrating the interactions between overlapping DEGs and alternative splicing genes. PPI, protein-protein interaction; DEGs, differentially expressed genes.

Table IV. Top 20 pairs of the protein-protein interaction network as determined by the highest combined score.

Table IV.

Top 20 pairs of the protein-protein interaction network as determined by the highest combined score.

Gene 1	Gene 2	Combined score
SMC4	SMC2	0.999
NOP56	DKC1	0.997
TFRC	OPTN	0.994
SRM	ODC1	0.989
CDA	TYMP	0.987
HNRNPM	U2AF2	0.979
ZWILCH	CENPE	0.977
CEACAM6	CEACAM1	0.970
TGM2	FN1	0.970
NOS3	AKT3	0.964
DTYMK	NME1	0.954
CENPE	KIF18A	0.946
NOP56	NOC2L	0.941
ATG7	ULK1	0.937
ARHGAP11A	RAC2	0.936
ZWILCH	KIF18A	0.930
OASL	IRF1	0.925
HDAC9	PHF21A	0.925
HDAC9	NOS3	0.923
DKC1	NOC2L	0.911

Table V. The top 10 nodes of the protein-protein interaction network as determined by the highest connectivity degree.

Table V.

The top 10 nodes of the protein-protein interaction network as determined by the highest connectivity degree.

Gene	Degree
AURKA	11
NOS3	8
NOC2L	8
CENPE	7
AKT3	7
NSUN2	6
SOD2	5
SMC4	5
RAC2	5
NOP56	5

Discussion

Lung cancer is a serious threat to human health and survival (22). Despite progress in diagnosis and treatment, the 5-year survival rate of patients with lung cancer is only 9–20% (23). FOXD3 has been suggested to be a tumor suppressor in various types of cancer (8–10). However, the underlying mechanism of FOXD3 activity in lung cancer remains unclear. In the present study, DEGs and ASGs between FOXD3-knockout and normal lung cancer A549 cells were identified, and functional enrichment analysis was performed to identify the associated biological processes involved in lung cancer. Finally, a PPI network of the most significant genes was constructed. These results may contribute to the understanding of the role of FOXD3 in lung cancer.

The most enriched GO terms for the DEGs were ‘response to wounding’, ‘extracellular region’, ‘plasma membrane’ and ‘immune response’. The ASGs were mainly enriched in ‘cytosol’, ‘intracellular organelle lumen’, ‘organelle lumen’ and ‘membrane-enclosed lumen’ (Fig. 2). The wound response involves clotting and coagulation, tissue remodeling, cellular migration and proliferation, and angiogenesis (24). The majority of these processes also serve important roles in the progression of cancer. One study reported that the upregulation of factors associated with the ‘wound response’ term was highly prognostic of breast cancer survival, and revealed a strong association between the pathogenic conditions identified by this signature and those identified using serum-treated fibroblasts (25). In lung cancer, the upregulation of genes associated with the ‘wound response’ term has been demonstrated as predictive of poor overall survival time and increased risk of metastasis (26).

The cell membrane is a biological membrane that separates the interior of cells from the outside environment (27). Plasma membrane fluidity depends on the composition of the lipids and proteins in the membrane, and has been demonstrated to be significantly associated with the malignant potential of cancer cells (28), with alterations in the plasma membrane fluidity of cancer cells associated with their capacity to form metastases (29). In lung cancer, studies reported that patients with high plasma membrane fluidity had poorer prognoses than those with less fluid membranes, and the fluidity variable may be used as an independent additional prognostic factor (28,30,31).

Cytosol is the fluid within cells, a component essential to the process of cytokinesis, a critical stage in cell proliferation (32,33). Another major function of cytosol is to transport metabolites; most tumor cells demonstrate different metabolic pathways to normal cells (34). One study indicated that metabolism contributed to the tumor proliferation, migration, and metastasis of lung cancer (35).

Other enriched GO terms, e.g., ‘organelle lumen’, have also been associated with tumorigenesis. Jingye et al (36) reported that a disordered pH in the organelle lumen is a common characteristic of cancer cells. Despite a number of studies reporting the FOXD3-mediated inhibition of the growth, invasion and migration of tumor cells in various types of cancer, including lung cancer (37–39), limited data is available regarding the association between FOXD3 and these GO terms. As discussed, the identified GO terms have been associated with the growth, invasion and migration of tumor cells, thus it is speculated that FOXD3 may affect the progression of lung cancer indirectly by regulating these biological processes.

From the identified KEGG pathways, the mechanistic target of rapamycin (mTOR) signaling pathway has also been associated with the growth and proliferation of tumor cells, and the deregulation of multiple elements of the mTOR pathway has been reported in numerous types of cancer (40). The NOD-like receptor signaling pathway is involved in the formation of inflammasomes, and numerous types of cancer are associated with inflamed tissue (41). However, the associations between FOXD3 and the identified KEGG pathways require further exploration.

A total of 199 overlaps between the DEGs and the ASGs were identified, from which the PPI network was constructed (Fig. 3). The top 5 nodes of the PPI network, with the highest degree, were aurora kinase A (AURKA), nitric oxide synthase 3 (NOS3), NOC2-like nucleolar associated transcriptional repressor (NOC2L), centromere protein E (CENPE) and AKT3. The majority of these genes have been previously associated with tumorigenesis. AURKA and NOS3 serve important roles in the development of various types of cancer, including lung cancer; AURKA is a cell cycle-regulated kinase involved in spindle formation and chromosome segregation (42). Various types of cancer exhibit the overexpression of AURKA, which is associated with chromosomal instability, centrosomal amplification/aneuploidy, therapeutic resistance, cell-cycle progression and anti-apoptosis. As an oncogene, AURKA is an important therapeutic target in lung cancer, and cell proliferation, apoptosis and cell cycle progression are associated with the expression of AURKA (43). NOS3 encodes an enzyme that regulates the production of nitric oxide and contributes to uncontrollable cell growth in a number of cancer types (44). Various studies have demonstrated associations between NOS3 and cancer processes. For example, Arıkan et al (45) reported that the NOS3 Glu298Asp polymorphism may be associated with the risk and progression of colorectal cancer. Lee et al (46) reported that genetic polymorphisms in NOS3 modified individual susceptibility to invasive breast cancer with lymph node involvement in Korean women. Furthermore, the expression of NOS3 has been reported to contribute to the tumor angiogenesis and lymph metastasis of human non-small cell lung cancer (47).

The expression of other genes, including CENPE, NOC2L and AKT3 has also been associated with tumorigenesis (48–50). CENPE was identified as a novel therapeutic candidate in neuroblastoma (50), and the selective activation of the AKT3 protein promoted cell survival and tumor development in non-familial melanomas in one study (48). To the best of our knowledge, there is no experimental evidence of the direct association between FOXD3 and these genes. However, the biological functions associated with these genes in the context of cancer correspond with the regulating mechanism of FOXD3 in lung cancer. FOXD3 acts as a tumor suppressor by regulating the expression of the target genes, thus inhibiting the growth, invasion and migration of tumor cells (51). Few specific targets for FOXD3 in lung cancer have been reported, whereas AURKA and NOS3 serve critical roles in the growth, invasion and migration of tumor cells in lung cancer. Therefore, we speculate that AURKA and NOS3 may be the targets of FOXD3 that execute its effect in lung cancer. Confirmation of these conclusions and further exploration of the specific mechanism of FOXD3 regulation in lung cancer are required.

In conclusion, FOXD3 serves an important role in regulating the growth, migration and proliferation of lung cancer cells. Genes such as AURKA and NOS3 may be targets of FOXD3, mediating its effect in lung cancer. The present study contributes to the existing understanding of the molecular mechanism of lung cancer and may provide data to contribute towards novel strategies for improving the diagnosis and therapy of lung cancer.

Acknowledgements

The present study was supported by the Municipal Science and Technology Commission of Tianjin (grant nos. 15ZLZLZF00440 and 16ZLZXZF00120) and the Health Bureau Science and Technology Foundation of Tianjin (grant no. 2014KZ102).

References