Association Analysis of Driver Gene–Related Genetic Variants Identified Novel Lung Cancer Susceptibility Loci with 20,871 Lung Cancer Cases and 15,971 Controls

Lung cancer has been one of the most commonly diagnosed malignancies and the leading cause of cancer-related death worldwide (1). The development of lung cancer is a multistep process that involves both genetic and environmental factors (2–4). Genome-wide association studies (GWAS) have been proven to be a powerful approach to dissect genetic architectures of complex diseases. To date, GWASs have identified 51 lung cancer susceptibility loci in various populations (5, 6). However, the information provided by GWAS remains inadequate. The heritability of lung cancer was estimated to be 20.6% in European populations (7), while only a small proportion of lung cancer heritability could be explained by risk loci that were identified in previous lung cancer GWASs (8). Therefore, more risk loci for lung cancer are needed to be identified.

Several waves of technology have facilitated the identification of lung cancer driver genes (lung CDG), which are improving our understanding of oncogenic process for lung cancer. On the basis of The Cancer Genome Atlas (TCGA) research on lung cancer, the most commonly mutated oncogenes in lung adenocarcinoma included KRAS, EGFR, BRAF, PIK3CA, and MET; mutations in tumor suppressors such as TP53, STK11, KEAP1, NF1, RB1, and CDKN2A were also frequently detected in lung adenocarcinoma (9–11). Although TP53, RB1, ARID1A, CDKN2A, PIK3CA, and NF1 were significantly mutated in both lung adenocarcinoma and lung squamous cell carcinoma (lung SqCC), significantly mutated genes like NOTCH1 and HRAS were only identified in lung squamous cell carcinoma (10–12). In addition to somatic mutations, somatic copy number alterations (SCNA) and rearrangements also play important roles in lung cancer development. Amplification of TERT and EGFR, as well as fusions involving ALK and ROS1, were commonly identified in lung adenocarcinoma. Deletions of CDKN2A have been identified in both lung adenocarcinoma and squamous cell carcinoma (9, 10, 12).

Emerging evidence has shown that a substantial proportion of CDGs are also cancer predisposition genes (13). The TCGA PanCanAtlas Germline Working Group identified 44 genes that showed coclustering or colocalization of pathogenic germline variants with recurrent somatic mutations, implying shared oncogenic processes in germline and somatic genomes (14). In addition, susceptibility variants could regulate the functions of nearby cancer driver genes. For example, rs2736100, a risk variant of lung cancer, is located in the first intron of driver gene TERT, and was associated with increased expression of TERT in lung tumors (15). However, the associations between common genetic variants in lung CDGs and lung cancer risk have rarely been explored. Therefore, we integrated lung CDGs, genetics of gene expression, and functional annotation databases with large-scale lung cancer GWAS datasets to systematically investigate the associations between lung CDG–related genetic variants and lung cancer risk.

This study utilized data from two existing GWASs of European descent: the OncoArray dataset (16) and Division of Cancer Epidemiology and Genetics (DCEG) Lung Cancer Study (17). The OncoArray dataset was derived from the Transdisciplinary Research of Cancer in Lung of the International Lung Cancer Consortium (TRICL-ILCCO) and the Lung Cancer Cohort Consortium (LC3). Quality control and imputation processes were described previously (16), resulting in 18,444 cases and 14,027 controls remained. The DCEG Lung Cancer GWAS data were obtained from dbGap phs000336.v1.p1 (17). Detailed quality control and imputation processes have been described previously (18). We further excluded individuals in the DCEG Lung Cancer Study that overlapped with or were related to individuals from the OncoArray dataset based on identity by descent (IBD) analysis (IBD > 0.45). As a result, a total of 2,427 cases and 1,944 controls from the DCEG Lung Cancer Study remained. All participants signed informed consents and study protocols were approved by the ethical review boards of each institution.

Genes were annotated as lung CDGs if they fulfilled any of the following criteria: (i) lung cancer–related genes in the COSMIC Cancer Gene Census (v78; ref. 19); (ii) mutational drivers, SCNA drivers, and fusion drivers detected by the IntOGen pipeline in lung tumors (20); and (iii) significantly mutated genes (SMG) and candidate CDGs with significant SCNAs that were identified in lung adenocarcinoma and/or lung squamous cell carcinoma by the TCGA projects (10).

To investigate functional variants in lung CDGs, we included single-nucleotide polymorphisms (SNP) if they satisfied either of the following criteria: (i) SNPs that were associated with expressions of lung CDGs (expression-related SNP, or eSNPs) in normal lung tissues based on the Genotype-Tissue Expression Project (GTEx, v6p release; P < 0.05; ref. 21) or (ii) nonsynonymous variants of lung CDGs identified using Variant Effect Predictor (22). The selected eSNPs and nonsynonymous variants were extracted from the two GWAS datasets. SNP with imputation INFO < 0.8, minor allele frequency (MAF) in controls <0.005, Hardy–Weinberg equilibrium (HWE) test P in controls < 1 × 10⁻⁷, or HWE test P in cases < 1 × 10⁻¹² was excluded from the analysis.

We performed logistic regression to generate odds ratios and confidence intervals (CI) for each SNP. The OncoArray dataset was used in the discovery stage with age, gender, and the first three principal components (PC) adjusted (16). Variations with association P < 0.001 were further tested in the DCEG Lung Cancer Study (the validation stage), and we adjusted age, gender, and the first PC in logistic regression model (23). SNPTEST v2.5 was used for the association analysis, taking dosage format of imputed genotypes. For variations with P < 0.05 in the validation stage, meta-analysis that combined effect estimates from the two datasets was performed using GWAMA v2.0.2 (24). The index of heterogeneity (I²) and P value based on Cochran Q test were calculated to assess the heterogeneity between studies. Fixed-effect model was used for absent of heterogeneity between studies (P_{heterogeneity} > 0.05); otherwise random-effect model was adopted. Variations with the same direction of effect in both GWAS datasets and P < 1 × 10⁻⁴ in the meta-analysis were considered as suggestive risk SNPs (Supplementary Fig. S1).

In addition to the overall lung cancer, we also investigated the associations of lung CDG–related SNPs with risk of lung adenocarcinoma and lung squamous cell carcinoma. As the DCEG Lung Cancer Study lacked information of histologic types, we performed association analysis using logistic regression model in the OncoArray dataset. To control the false discovery rate (FDR), we used Benjamini–Hochberg step-down method to calculate FDR for each variation. Variations with FDR < 0.01 were considered as suggestive risk SNPs.

We mapped suggestive risk SNPs to lung CDGs based on the GTEx v6p release, and performed functional prediction for significant nonsynonymous variants using SIFT (25) and PolyPhen2 (26), which were implemented in ANNOVAR (27). For lung CDGs with multiple risk SNPs, conditional and joint association analysis were performed to identify independent signals using genome-wide complex trait analysis (GCTA; ref. 28). During the model selection process, the testing SNP was not selected if its regression R² on the selected SNPs was greater than 0.1. The threshold P value of 0.0001 was adopted to identify significant independent hits. SNPs that were significant after the multiple testing corrections and that were not in linkage disequilibrium (LD, r² < 0.1) with and were located at least 500 kilobases apart from known risk variants were considered as novel susceptibility SNPs.

Expression data on 56,238 genes for 320 normal lung tissues were downloaded from the GTEx website (21). Genome-wide expression correlation analysis was performed using a linear regression model to identify genes coexpressed with significant lung CDGs. Significant coexpressed genes that satisfied the Bonferroni correction were used for pathway enrichment analysis. We downloaded pathways in the Kyoto Encyclopedia of Genes and Genomes (KEGG) database from MSigDB (29–31), and performed pathway enrichment analysis using “PHYPER” function as implemented in R software (version 3.4.1), which computes a P value for each pathway based on hypergeometric distribution.

TCGA datasets of lung adenocarcinoma and lung squamous cell carcinoma were used to model the association between independent SNPs and somatic alterations in lung CDGs (9, 12). Germline genotype data generated using Affymetrix Genome-Wide Human SNP Array 6.0 were applied for and approved in Feb, 2015. Standard quality control and genotype imputation process have been described previously (32).

We downloaded Mutation Annotation Format derived from whole-exome sequencing, as well as somatic copy numbers calculated using GISTIC2 from the Broad Institute Genome Data Analysis Center (GDAC) Firehose portal (stamp analyses_2016_01_28) (33). For each patient, a lung CDG was considered mutated if one or more somatic mutations mapped to this gene. We also assessed truncation mutations (frame shift insertion/deletion, nonsense, nonstop, and splice site mutations; ref. 34) in pathways that are recurrently altered in lung cancer, including cell cycle, spliceosome, Notch signaling pathway, transcriptional misregulation in cancer, Ras signaling pathway, and PI3K–Akt signaling pathway (11). A pathway was considered as mutated if one or more truncation mutations were observed in this pathway. We used logistic regression models to evaluate the association between independent SNPs and mutational status of lung CDGs or pathways. In the analysis of SCNAs, somatic copy number of lung CDG was used as outcome, and we used linear regression to model the association between independent SNPs and SCNAs. Age, gender, smoking status, clinical stage, and the first 10 PCs were adjusted as covariates. The association analysis between independent SNPs and somatic alterations were performed in lung adenocarcinoma and lung squamous cell carcinoma, separately. Benjamini–Hochberg step-down method was used to calculate FDR for each SNP-lung CDG (or SNP pathway) pair to control the FDR. Association analysis was conducted using the R software (version 3.4.1).

The OncoArray dataset included 18,444 cases and 14,027 controls. The mean (± SE) age of the subjects was 63.79 ± 10.44 for cases and 61.77 ± 10.29 for controls. For the DCEG Lung Cancer Study, a total of 2,427 cases and 1,944 controls were included. Among participants across both studies with known histologic types, there were 6,819 lung adenocarcinomas and 4,490 lung squamous cell carcinomas. Detailed characteristics and clinical features of participants in each data set were shown in Supplementary Table S1.

A total of 348 protein-coding lung CDGs were included from published data (Supplementary Table S2). We identified 139,666 eSNPs and 2,041 nonsynonymous variants of lung CDGs. Among SNPs that passed the quality control process, a total of 234 SNPs were identified in the OncoArray dataset (P < 0.001) and validated in the DCEG Lung Cancer Study (P < 0.05), which were mapped to five lung CDGs (Supplementary Table S3). After conditional analysis, seven independent signals were identified. Among these loci, rs78062588, which was mapped to TPM3 in chromosome 1q21.3, was a new lung cancer susceptibility locus [OR = 0.87, 95% CI: 0.81–0.92, P = 1.55 × 10⁻⁵ in the OncoArray dataset; OR = 0.82, 95% CI: 0.68–0.98, P = 3.11 × 10⁻² in the DCEG Lung Cancer Study; and OR = 0.86, 95% CI: 0.81–0.91, P = 1.65 × 10⁻⁶ in the meta-analysis; Tables 1 and 2; Supplementary Table S3]. In addition, rs71658797 in FUBP1 (1p31.1), rs1655931 and rs2517586 in HLA-A (6p22.1), rs2887532 in KDM5A (12p13.33), rs7359276 and rs7161774 in IREB2 (15q25.1) had been reported by previous GWASs as lung cancer susceptibility loci (Tables 1 and 2; Supplementary Table S3; refs. 5, 6).

Stratified analyses in lung adenocarcinoma and lung squamous cell carcinoma found another nine susceptibility genes, including seven genes that were identified in lung adenocarcinoma and two genes that were identified only in lung squamous cell carcinoma (Fig. 1A and B; Supplementary Table S4). Independent variants derived from conditional analysis are shown in Supplementary Table S5. Of these loci, rs2700389 in KALRN (3q21.1), rs79518818 in MGA (15q15.1), and rs62054832 in EFTUD2 (17q21.31) were first identified as risk loci for lung adenocarcinoma, while rs148797791 in IRF6 (1q32.2) was found as a novel risk locus for lung squamous cell carcinoma. SNPs rs7823498 in NRG1 (8p12), rs10757256 and rs1011970 in CDKN2A (9p21.3), rs79040073 in COPS2 (15q21.1), rs2281925 in ARFGAP1 (20q13.33), and rs17879961 in CHEK2 (22q12.1) had been reported by previous GWASs as lung cancer susceptibility loci (5, 6).

Among 234 significant SNPs in overall lung cancer, three were nonsynonymous variants. Two additional nonsynonymous variants (rs1136688 in HLA-A and rs17879961 in CHEK2) were identified in lung squamous cell carcinoma (Supplementary Table S6). We predicted functional consequence of nonsynonymous variants using SIFT and Polyphen-2 (25, 26). Notably, risk variant rs707910 in HLA-A (NM_001242758, c.G203A) was predicted as deleterious by SIFT and possibly damaging by Polyphen-2. SNP rs17879961 in CHEK2 (NM_007194, c.T470C) was predicted as tolerated by SIFT and possibly damaging by Polyphen-2.

To explore biological processes underlying significant lung CDGs, we performed genome-wide coexpression and KEGG pathway enrichment analysis. We identified essential pathways in lung carcinogenesis such as apoptosis, MAPK signaling pathway, spliceosome, cell cycle, and nucleotide excision repair (Supplementary Table S7; ref. 11).

We investigated the associations between independent SNPs and somatic alterations in lung CDGs. The protective rs78062588[C] allele (TPM3 in 1q21.3) was associated with increased expression of TPM3 in normal lung tissues (OR = 1.14, P = 0.04) and elevated somatic copy number of TPM3 in TCGA lung adenocarcinomas (OR = 1.16, P = 0.02; Supplementary Fig. S2). However, the analysis of somatic mutations in lung CDGs did not identify any association with P < 0.05. As the mutational frequencies of lung CDGs are relatively low, we further analyzed the associations between independent risk SNPs and truncation mutations at the pathway level. Among patients with lung squamous cell carcinoma, we found that rs1611182 (HLA-A in 6p22.1), a risk SNP for lung adenocarcinomas, was associated with decreased frequency of truncation mutations in the transcriptional misregulation in cancer pathway (OR = 0.66, 95% CI: 0.50–0.85, P = 1.76 × 10⁻³, FDR < 0.25; Table 3; Supplementary Table S8; Supplementary Fig. S3).

This study comprehensively incorporated lung cancer GWASs, lung CDGs, genetics of gene expression, somatic alterations in lung tumors, and functional annotation databases to investigate the associations of CDG-related genetic variants with lung cancer risk. We identified five lung CDGs in overall lung cancer. Subgroup analysis by histologic types further identified seven and two genes in lung adenocarcinoma and lung squamous cell carcinoma, respectively. Genes coexpressed with the identified lung CDGs were involved in essential pathways including cell cycle, MAPK signaling, and nucleotide excision repair pathways. Incorporation of somatic alterations identified lung cancer risk variants that were associated with somatic alterations in lung CDGs or recurrently mutated pathways.

TPM3 is included in the COSMIC Cancer Gene Census. Translocation of TPM3 could form oncogenic fusion proteins, such as TPM3-ROS1 observed in advanced lung adenocarcinoma (35). Previously conducted functional assessment in NIH3T3 cells showed that TPM3-ALK fusion protein can interact with endogenous tropomyosin, which may induce changes in cell morphology and cytoskeleton organization and further bestowed higher metastatic capacities (36). Our results found that the protective allele of rs78062588 was associated with increased TPM3 expression as well as increased somatic copy number alterations of TPM3 in lung adenocarcinomas. However, reaching a better understanding of the functional impact of TPM3 on lung cancer development warrants further investigation.

CDKN2A in 9p21.3 encodes several alternatively spliced transcripts, among which are p16 and ARF. p16 is a tumor suppressor that functions as an inhibitor of CDK4 and CDK6 (37). Another tumor suppressor protein, ARF, functions as a stabilizer of the tumor suppressor protein p53. Both p16 and ARF have functionality in cell-cycle G₁ control. CDKN2A is recognized as an important tumor suppressor gene. Deletion of CDKN2A was frequently identified in lung tumors (10). In addition, CDKN2A has been identified as susceptibility gene for lung adenocarcinoma (16). We validated this locus and identified a second signal within CDKN2A. Consistently, the risk alleles of independent SNPs were associated with decreased expression of CDKN2A in normal lung tissues.

The transcription factor interferon regulatory factor 6 (IRF6) was identified as significantly mutated gene in TCGA lung squamous cell carcinomas (10). IRF6 has essential role in epidermal development. It is induced in differentiation through a Notch-dependent mechanism. Downregulation of IRF6 in epithelial squamous cell carcinomas promotes ras-induced tumor formation and reintroduction of IRF6 strongly inhibits cell growth (38, 39). The tumor suppressor role of IRF6 has also been demonstrated in vulvar squamous cell carcinoma (40). In addition, elevated IRF6 expression in nasopharyngeal carcinomas suppressed cell proliferation and growth (41). We identified IRF6 as a susceptibility gene for lung squamous cell carcinoma. Consistent with the tumor suppressor role of IRF6, the risk allele of rs148797791 was associated with decreased expression of IRF6 in normal lung tissues. These results indicate that germline variant might contribute to lung cancer risk by downregulation of IRF6.

Genes coexpressed with the identified lung CDGs were enriched in essential pathways such as apoptosis, MAPK signaling pathway, spliceosome, cell cycle, and nucleotide excision repair. A comprehensive molecular profiling of lung adenocarcinoma demonstrated recurrent somatic alterations in cell cycle and MAPK signaling pathway (9, 42). In addition, deregulated RNA Splicing is involved in lung adenocarcinoma, and cell-cycle pathway is involved in both lung adenocarcinoma and lung squamous cell carcinoma (11, 42).

We comprehensively collected 348 lung CDGs from three databases, and tested associations between functional SNPs of lung CDGs and risk of lung cancer in large-scale lung cancer GWASs of Europeans. We identified five novel susceptibility loci of lung cancer, and validated nine loci that had been reported by previous lung cancer GWASs. These results showed that genetic variants in lung CDGs contribute to lung cancer susceptibility. Our findings might help to unravel biological functions of lung cancer susceptibility loci.

No potential conflicts of interest were disclosed.

The content is solely the responsibility of the authors and does not necessarily represent the official views of the NIH. Where authors are identified as personnel of the International Agency for Research on Cancer/World Health Organization, the authors alone are responsible for the views expressed in this article and they do not necessarily represent the decisions, policy, or views of the International Agency for Research on Cancer/World Health Organization.

Conception and design: Y. Wang, M. Zhu, O. Melander, S. Zienolddiny, H. Shen

Development of methodology: Y. Wang

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): D. Albanes, S. Lam, A. Tardon, C. Chen, G.E. Goodman, S.E. Bojesen, M.T. Landi, M. Johansson, A. Risch, H.-E. Wichmann, H. Bickeboller, D.C. Christiani, G. Rennert, S.M. Arnold, P. Brennan, J.K. Field, S. Shete, L. Le Marchand, O. Melander, H. Brunnstrom, G. Liu, R.J. Hung, A.S. Andrew, L.A. Kiemeney, S. Zienolddiny, K. Grankvist, M. Johansson, N.E. Caporaso, P.J. Woll, P. Lazarus, M.B. Schabath, M.C. Aldrich, C.I. Amos, H. Shen

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): Y. Wang, I.P. Gorlov, M. Zhu, J. Dai, G.E. Goodman, D.C. Christiani, G. Rennert, S.M. Arnold, R.J. Hung, N.E. Caporaso, H. Shen

Writing, review, and/or revision of the manuscript: Y. Wang, O.Y. Gorlova, J. Dai, D. Albanes, S. Lam, A. Tardon, C. Chen, G.E. Goodman, S.E. Bojesen, M.T. Landi, M. Johansson, A. Risch, H.-E. Wichmann, H. Bickeboller, D.C. Christiani, G. Rennert, S.M. Arnold, P. Brennan, J.K. Field, S. Shete, L. Le Marchand, O. Melander, H. Brunnstrom, G. Liu, R.J. Hung, L.A. Kiemeney, K. Grankvist, M. Johansson, N.E. Caporaso, P.J. Woll, M.B. Schabath, M.C. Aldrich, V.L. Stevens, H. Ma, G. Jin, Z. Hu, C.I. Amos, H. Shen

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): Y. Wang, G.E. Goodman, M. Johansson, D.C. Christiani, J.K. Field, G. Liu, S. Zienolddiny, M. Johansson, M.B. Schabath, C.I. Amos, H. Shen

Study supervision: R.J. Hung, L.A. Kiemeney, G. Jin, H. Shen

The authors thank the study participants and research staff for their contributions and commitment to this study. H. Shen was awarded grants from National Natural Science of China (81820108028, 81521004). C.I. Amos was awarded grants from the NIH (U19 CA148127, U19 CA203654) and the Cancer Prevention Research Institute of Texas (RR170048). H. Ma was awarded a grant from National Natural Science of China (81922061). D.C. Christiani was awarded a grant from the NIH (U01 CA209414). C. Chen was awarded grants from the NIH (U01-CA063673, UM1-CA167462, and U01-CA167462). M.B. Schabath was awarded the Moffitt Cancer Center Support Grant (P30 CA076292) and SPORE in Lung Cancer (P50 CA119997).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.