Genetic analyses of differences between solid and nonsolid predominant lung adenocarcinomas.

BACKGROUND: Solid predominant lung adenocarcinomas (LUAD) have distinct histopathological and clinical characteristics compared with nonsolid subtypes. A comprehensive comparison of altered genes found in solid and nonsolid subtypes has not previously been performed. In this study, we analyzed differences in gene expression, genetic mutations, and DNA methylation to better understand the risk factors for these two subtypes of LUAD.
METHODS: Differentially expressed genes (DEGs) and differentially mutated genes (DMGs) were analyzed from RNA-seq data downloaded from The Cancer Genome Atlas (TCGA) and Broad Institute database. To understand the functional significance of molecular changes, we examined the DEGs and DMGs with Gene Ontology and Kyoto Encyclopedia of Genes and Genomes pathway analysis.
RESULTS: A total of 184 patients in the TCGA cohort and 140 patients in the Broad Institute cohort were included in this study. We identified 75 DEGs, of which 15 were upregulated and 56 downregulated in the solid group relative to the nonsolid group. The DEGs were mainly involved in the regulation of water and fluid transport. We discovered 38 significantly differentially expressed genes that overlapped in the two groups. The DMGs were mainly enriched for pathways involved in cell-cell adhesion, cell adhesion, biological adhesion, and hemophilic cell adhesion. We additionally discovered nine significantly methylated genes between solid and nonsolid LUAD.
CONCLUSIONS: Our study identified distinct DEGs, DMGs, and methylation genes for solid and nonsolid LUAD subtypes. These findings improve our understanding of the different carcinogenesis mechanisms in LUAD and will help to develop new therapeutic strategies.

Introduction

Lung cancer remains the leading cause of cancer‐related death worldwide.1 Adenocarcinoma is the most common histological type, and presents extensive intratumor heterogeneity. More than 80–90% of lung adenocarcinomas (LUADs) demonstrate mixed morphologic patterns.2 According to the novel LUAD classification system proposed by the International Association for the Study of Lung Cancer (IASLC), the American Thoracic Society (ATS), and the European Respiratory Society (ERS) in 2011, after performing comprehensive histological subtyping with semiquantitative assessment of each subtype in 5% increments, invasive LUAD is classified into one of the following subtypes: lepidic, acinar, papillary, micropapillary, and solid.3 The IASLC/ATS/ERS histological classification system is a reliable and powerful tool to predict prognosis.4, 5, 6 Survival outcomes are significantly different among the five subtypes. The lepidic subtype predicts an excellent prognosis, the acinar and papillary subtypes exhibit intermediate clinical behavior, and the solid and micropapillary subtypes are associated with poor prognosis. In particular, the presence of a solid pattern in the nonsolid subtype is correlated with reduced disease‐free survival (DFS).7, 8

The solid predominant subtype exhibits distinct molecular and clinical characteristics compared to the nonsolid subtypes. The solid subtype presents a unique molecular profile, in that patients with the solid pattern have a less frequent rate of EGFR mutations but are more likely to harbor KRAS mutations.9, 10 Clinically, the solid predominant pattern is associated with poor prognostic factors, including a higher rate of lymph node metastasis,11 tumor spread through air spaces (STAS),12, 13 early recurrence, and a high incidence of extrathoracic and multiple‐site recurrence.14 However, explanations as to why the solid predominant subtype is associated with aggressive biological behavior are limited to driver‐mutation genes. Therefore, comprehensive investigations into the differences in the molecular characteristics between the solid and nonsolid LUAD subtypes are imperative, which will lead to a deeper understanding of the pathogenic mechanisms of solid subtype LUAD.

In this study, we explored differences in gene expression, mutated genes, DNA methylation, Gene Ontology (GO) biological annotations, and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways between solid and nonsolid LUAD by analyzing clinical samples derived from The Cancer Genome Atlas (TCGA) and the Broad Institute (BI) database. We then analyzed survival curves for carriers with low and high expression of the most distinct differentially expressed genes (DEGs) in a selected TCGA cohort. Our study aimed to provide a comprehensive perspective into the underlying molecular mechanisms, prognostic predictive biomarkers, and therapeutic targeted genes for solid predominant LUAD.

Methods

Patient cohort

Messenger RNA expression profiles and DNA methylation data (combining level 3 data from Illumina GA and HTSeq platforms), as well as clinical data of lung adenocarcinoma patients were downloaded from TCGA (https://portal.gdc.cancer.gov/). The histologic subtypes of cases were obtained from the supplementary data of previously published studies (http://www.nature.com/nature/journal/v511/n7511/full/nature13385.html#supplementaryinformation). DNA variant data was downloaded from TCGA (https://tcga-data.nci.nih.gov/tcga/findArchives.htm) an comprised the TCGA cohort and data from cbioPortal (http://www.cbioportal.org/study.do?cancer_study_id=luad_broad) comprised the BI cohort.

Lung adenocarcinoma was classified according to the 2011 IASLC/ATS/ERS classification system. Invasive adenocarcinomas were classified into lepidic, acinar, papillary, micropapillary, and solid subtypes based on the predominant histological pattern present in the tumor. Patients were divided into solid and nonsolid (lepidic, acinar, papillary, micropapillary) groups. Invasive adenocarcinoma variant subtypes and cases for which RNAseq and gene mutation data were not available were excluded. Finally, a total of 184 (57 solid and 127 nonsolid) patients were included in TCGA cohort and 140 (46 solid and 94 nonsolid) in the BI cohort.

RNA‐seq data preprocessing

Human gene annotations were downloaded from GENCODE (v25; http://www.gencodegenes.org). Expression profiles were measured as fragments per kilobase per million (FPKM) values using the FPKM function in the “DESeq2” package (http://www.bioconductor.org/packages/release/bioc%20/html/DESeq2.html) and were then log2 transformed. Genes were considered robustly expressed if their raw read counts were larger than 0 in at least half of all samples. If a gene contained multiple transcripts, the expression level of the longest transcript was chosen. Only protein coding genes were used in the downstream analysis.

Identification of differentially expressed genes (DEGs)

Differential expression analysis across all samples between the solid and nonsolid cohorts was performed using the Wald significance test defined by the nbinom Wald test function, as implemented in DESeq2 package in R, which is based on raw read counts for each gene and is a robust method for analyzing RNA‐seq data. The statistical threshold for significance was a false discovery rate (FDR) < 0.05 and fold change > 2.

Survival analysis

In order to verify whether the upregulated genes in the solid group were associated with poor survival, the patients were divided into two groups per gene according to the expression value: a high expression group (> median expression level across all samples) and a low expression cohort (≤ median expression level across all samples). Kaplan–Meier curves of the high and low expression cohorts were created and compared using the log‐rank test. The Cox proportional hazard model was utilized for multivariate analysis to identify factors associated with survival. Adjusted hazard ratios (HRs) and 95% confidence intervals (CIs) are reported. The significance level for all tests was two‐sided at P < 0.05. The procedure was performed using the “survival” package (https://cran.r-project.org/web/packages/survival/index.html) in R.

Identification of differentially mutated genes (DMGs)

Only 104 cases in the TCGA cohort had complete mutation data, including 34 cases with the solid subtype and 70 with the nonsolid subtype. Somatic mutation data on 104 cases in the TCGA cohort and 140 in the BI cohort were integrated to analyze the differentially mutated genes (DMGs) between the solid and nonsolid groups. Genes carrying silent mutations were discarded. The Fisher's exact test was used to compare the difference in the mutation ratio for each gene and a P value of ≤ 0.05 was considered statistically significant.

Identification of differentially methylated genes

There were 24 nonsolid and 22 solid cases in the TCGA cohort containing DNA methylation data. To calculate the methylation level of expressed genes, we screened the C‐phosphate‐G (CpG) sites in the gene promoter regions. For a gene containing multiple CpG sites, the methylation level of a gene was estimated with the average methylation level of all CpG sites. For each CpG site, methylation must have occurred in at least 80% of the samples. For missing values, we used the knnImputation function of DMwR package in R to replace the missing values. A differentially methylated gene was defined as the absolute difference in the mean methylation level of the gene in solid and nonsolid subtypes ≥ 0.1 and P value < 0.05 in a Wilcoxon test.

Functional enrichment analyses

Gene Ontology and KEGG analysis was applied to identify the main functions of the DEGs and DMGs. Functional annotation was implemented using the Database for Annotation, Visualization and Integrated Discovery (https://david.ncifcrf.gov/). The background was set to all protein‐coding genes in humans. The statistical significance threshold for all of the inclusive GO terms and KEGG pathways was P < 0.05.

Results

The TCGA cohort included 10 (5.4%) lepidic, 68 (40.0%) acinar, 26 (14.1%) papillary, 23 (12.5%) micropapillary, and 57 (31.0%) solid subtype patients. The numbers of patients in stages I, II, III, and IV were 94 (51.1%), 40 (21.7%), 36 (19.6%), and 11 (6.0%), respectively. TNM stage data was not available for 3 (1.6%) patients. The BI cohort included 13 (9.3%) lepidic, 46 (32.9%) acinar, 18 (12.9%) papillary, 17 (12.0%) micropapillary and 46 (32.9%) solid subtype patients. The numbers of patients in stages I, II, III, and IV were 70 (50.0%), 27(19.3%), 21(15.0%), and 9 (6.4%), respectively. TNM stage data was not available for 13 (9.3%) patients.

Identification of DEGs

The DEGs between the solid and nonsolid groups were selected with strict criteria of fold change > 2 and an FDR < 0.05. As a result, a total of 75 genes were identified: 19 were significantly upregulated and 56 were significantly downregulated. A volcano plot and a heatmap were constructed to show all of the genes that were differentially expressed between solid and nonsolid LUAD (Figs 1, 2). The top 10 elevated DEGs in the solid group were: NTS, S100A7, COL2A1, SBSN, CALML3, GABRA2, KRT6A, IGFBP1, CYP1A1, and SLC6A15. The top 10 downregulated DEGs in the solid group were: PCSK2, PGC, REG1A, LGALS4, TMEM229A, SLC14A2, GKN2, SCGB3A1, and CAPN6.

Figure 1

Volcano plots showing all genes in the solid and nonsolid subtypes according to P value and fold changes: red dots represent genes that were not differentially expressed, while blue dots represent differentially expressed genes. Statistical significance was set at a false discovery rate < 0.05 and fold change > 2. () P < 0.05 and () P ≥ 0.05.

Figure 2

Heatmap of differentially expressed genes between solid and nonsolid lung adenocarcinoma subtypes.

To further analyze the correlation between elevated DEGs and prognosis, we performed survival analysis using TCGA data. Patients who harbored a high CTCFL expression level had a strong correlation with poor OS (HR 1.83, 95% CI 1.17–2.88; P = 0.008). Patients with high expression of HES2 (HR 1.55, 95% CI 0.99–2.41; P = 0.054), SBSN (HR 1.45, 95% CI 0.94–2.25; P = 0.095), and S100A7 (HR 1.42, 95% CI 0.92–2.21; P = 0.117) had significantly shorter OS than patients with low expression (Table 1).

Table 1

Correlations between differentially expressed gene expression levels and overall survival outcome

Gene	Log2FC	P	HR	95% CI	P
CTCFL	2.44	< 0.001	1.83	1.17–2.88	0.008
HES2	1.95	< 0.001	1.55	0.99–2.41	0.054
SBSN	3.34	< 0.001	1.45	0.94–2.25	0.095
S100A7	4.75	< 0.001	1.42	0.92–2.21	0.117

CI, confidence interval; FC, fold change; HR, hazard ratio.

Identification of DMGs

We next compared the DMGs between the solid and nonsolid groups. The results showed a significantly increased mutational load in the solid compared to the nonsolid subtypes. In the TCGA cohort, we identified 356 significantly differentially mutated genes between solid and nonsolid groups. Genes with P < 0.01 in the Fisher's test are shown in Figure S1. The top five DMGs with significantly increased mutation frequency in the solid group were: COL14A1, PRKCB, TP53, GRP158, and PTPRT.

In the BI cohort, we identified 639 significantly differentially mutated genes between the solid and nonsolid groups. Genes with P < 0.01 in the Fisher's test are shown in Figure S2. The top five DMGs with significantly increased mutation frequency in the solid group were MYH8, NF1, PCDH17, SYNE1, and TTN. Figure 3 shows the 38 significantly DMGs in both cohorts, including CSMD3, MUC16, FATS, TP53, and ATG2B.

Figure 3

Genetic mutation profiles in solid and nonsolid lung adenocarcinoma detected in both The Cancer Genome Atlas and Broad Institute cohorts. The results showed a higher frequency of genetic mutation in the solid than nonsolid subtypes. P < 0.05.

We identified nine differentially methylated genes between the solid and nonsolid samples, including LOR, DRD4, GAL3ST3, IGFBP6, SOX1, CLDN9, SMTNL2, FCAR, and DUSP22 (Fig 4).

Figure 4

Heatmap of differentially methylated genes between solid and nonsolid lung adenocarcinoma subtypes.

Enriched Gene Ontology and pathways

Gene Ontology and KEGG analyses were used to understand the altered biological functions and pathways of the identified DEGs and DMGs. For DEGS, there were 15 enriched functional categories. The most significantly enriched GO and KEGG terms are listed in Figure 5. The DEGs were mainly enriched for the regulation of water and fluid transport.

Figure 5

Gene Ontology and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analyses performed using the Database for Annotation, Visualization and Integrated Discovery on differentially expressed genes between solid and nonsolid lung adenocarcinoma. GnRH, gonadotropin‐releasing hormone.

For the DMGs, we included the 38 DMEs identified in both cohorts. There were nine enriched functional categories, and the most significantly enriched GO and KEGG terms are listed in Figure 6. The DMGs were mainly enriched in cell–cell adhesion, cell adhesion, biological adhesion, chemical homeostasis, cation homeostasis, and hemophilic cell adhesion.

Figure 6

Gene Ontology and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analyses performed using the Database for Annotation, Visualization and Integrated Discovery on differentially mutated genes between solid and nonsolid lung adenocarcinoma.

Discussion

Recent studies have indicated that the solid subtype of LUAD has distinct clinical characteristics compared to nonsolid subtypes, and tremendous genomic diversity among LUAD subtypes has been shown.4, 15, 16, 17 In order to further discover critical molecular and cellular mechanisms driving the solid pattern, LUAD initiation, maintenance, recurrence, and metastasis, we analyzed the differences in DEGs, DMGs, and differently methylated genes between solid and nonsolid LUAD samples using TCGA and BI data. To further elucidate the functions of these DEGs and DMGs, we performed GO and KEGG pathway analyses. However, further independent validation with experimental data is still required.

It is well known that a solid predominant pattern in LUAD is strongly associated with poor prognosis, especially in early‐stage LUAD patients. Our previous study indicated that for stage IB LUAD patients, the five‐year overall survival was only 60% in the solid group, while it was 80% in the acinar and papillary groups.5 The risk of postsurgical recurrence peaked significantly earlier for patients with stage I solid subtypes (within 12 months) than nonsolid subtypes (within approximately 24 months), and the majority of postoperative recurrences were detected in an asymptomatic condition during the regular follow‐up period.14, 18 Advanced analysis demonstrated that even a minor solid pattern mixed in nonsolid subtypes predicts a significantly poorer prognosis compared to a solid‐absent pattern, and a solid minor pattern is also associated with a higher incidence of lymph node metastasis,11, 19 revealing that the solid pattern has more aggressive behavior.

Based on these clinical phenomena, we hypothesized that the DEGs and DMGs between the solid and nonsolid subtypes may play a crucial role in the pathogenesis, development, tumor recurrence, and metastasis of solid subtype LUAD. Under strict selection conditions, we identified 75 DEGs between the solid and nonsolid groups. Among the top 10 upregulated genes in the solid group, we discovered that high expression of CTCFL, HES2, S100A7, and SBSN were associated with decreased OS. S100A7 is a calcium‐binding protein that belongs to the S100 family.20 Proteins in the S100 family are involved in the regulation of a number of cellular processes, such as cell cycle progression and differentiation.21 Previous studies have revealed that elevated S100A7 is associated with tumor aggressiveness and metastasis.22, 23, 24 Hu et al. reported that S100A7 expression is significantly increased in lung cancer, and is associated with poor prognosis.25 Nasser et al. studied the expression level of S100A7 in a breast cancer model and discovered that S100A7 enhanced breast cancer cell growth and metastasis and may act via activating proinflammatory and metastatic pathways.26 Although little research has focused on SBSN, Alam et al. reported that SBSN was involved in the angiopoiesis of tumor endothelial cells.27 Additional experiments should evaluate the therapeutic value of S100A7 and SBSN.

Several previous studies have revealed that solid subtype LUAD is more likely to harbor KRAS and TP53 mutations but less likely to harbor EGFR mutations compared to the nonsolid subtypes.9 However, there is still a lack of effective therapeutic agents targeting KRAS because the blockage of key KRAS effectors may lead to the activation of compensatory or parallel pathways.28 In order to discover more druggable mutated genes, we identified 356 significantly differentially mutated genes in the TCGA cohort and 639 in the BI cohort. Only 38 genes overlapped in the two groups. MUC16 encodes CA‐125, which was frequently mutated in both cohorts. Li et al. reported that MUC16 is also frequently mutated in gastric cancer, and MUC16 mutations are associated with a higher tumor mutation load (TML), and better survival outcomes and immune response in patients with gastric cancer.29 Recent studies have demonstrated that solid subtype LUAD is associated with an increased tumor mutation burden and PD‐L1 expression level compared to nonsolid subtypes. We hypothesize that MUC16 mutation may also be a useful marker to reflect TML in lung adenocarcinoma, and these findings may shed light on immunotherapies beneficial to LUAD patients with the solid subtype.

Some of the DMGs that were only detected in one of our study cohorts require further research. PRKCB is a family of serine‐specific and threonine‐specific protein kinases that can be activated by calcium. The protein encoded by PRKCB may be involved in apoptosis induction, endothelial cell proliferation, and several other cellular signaling pathways.30 However, the role of this gene in the pathogenesis of lung cancer remains unclear. The protein encoded by PTPRT is a family of the protein tyrosine phosphatase (PTP) family.31 Proteins in this family participate in cell growth, differentiation, the mitotic cycle, and oncogenic transformation. A previous report indicated that PTPRT mutations have been identified in 10% of lung cancers; our results revealed that the mutation rate of PTPRT in the solid group was 18.2%, while it was only 2.5% in the nonsolid group.32 PTPRT largely contains missense mutations in tumors, while only a small portion of PTPRT mutations are nonsense, insertion, and deletion mutations, and most of the tumor‐derived extracellular domain mutations of PTPRT impair cell–cell adhesion.33 We hypothesized that PTPRT mutation may play an important role in the process of tumor metastasis of solid LUAD, while further analyses may valid PTPRT as a therapeutic target for LUAD.

Interestingly, we found that only one gene was both differentially expressed and differentially mutated between the solid and nonsolid groups. CASR encodes a plasma membrane G protein‐coupled receptor that senses small changes in circulating calcium concentration. CASR is deregulated in solid LUAD compared to the nonsolid subtypes. Previous reports have indicated that CASR dysfunction is associated with increased aggressiveness and unfavorable outcomes in breast cancer34 however, the role of CARS in the carcinogenesis of LUAD remains unknown.

To further explore the underlying biological processes that DEGs and DMGs are associated with, we performed functional enrichment analysis. The enriched GO and KEGG terms from the DEGs were primarily associated with intracellular and extracellular material transport, including water, fluid, urea, and amide transport. This may reflect differences in metabolism between these two groups. A distinct feature of solid subtype LUAD is early, multisite, and extrathoracic metastasis. GO and KEGG enrichment analysis indicated that the DMGs were mainly associated with cell–cell adhesion, cell adhesion, biological adhesion, chemical homeostasis, cation homeostasis, and hemophilic cell adhesion. These results reveal that cell adhesion malfunction may play an important role in the process of solid LUAD metastasis.

There are several limitations to our study. First, because the results were derived from TCGA and the BI database, the DEGs and DMGs might be different between Asians and Caucasians. Second, the median follow‐up time in TCGA cases was inadequate, which may affect the accuracy of survival outcomes. Furthermore, the prognostic predictive value of selected genes has not been validated using our tumor specimens.

The present study investigated the differences in DEGs, DMGs, and differentially methylated genes between LUAD solid and nonsolid subtypes, and explored their different biological characteristics using GO and KEGG analysis. Further experiments are required to validate our findings and further functional investigations of targeted genes are needed to explore the molecular mechanisms underlying the aggressive behavior of solid subtype LUAD.

Disclosure

No authors report any conflict of interest.

Figure S1. Genetic mutation profiles in solid and nonsolid lung adenocarcinoma detected in The Cancer Genome Atlas. The results showed a higher frequency of genetic mutation in the solid than in the nonsolid subtypes. P < 0.01.

Click here for additional data file.

Figure S2. Genetic mutation profiles in solid and nonsolid lung adenocarcinoma detected in the Broad Institute cohort. The results showed a higher frequency of genetic mutation in the solid than in the nonsolid subtypes. P < 0.01.

Click here for additional data file.