Genetic variants of genes in the NER pathway associated with risk of breast cancer: A large-scale analysis of 14 published GWAS datasets in the DRIVE study.

The overall analysis included 28758 BC cases and 24349 controls from 14 studies ( Supplementary Table 2 ). All subjects were women. The median age of controls was 60 years. The age distribution was statistically different between cases and controls ( P <0.0001), with the control group being younger than the case group (≤ 60 years: 52.47% versus 49.73%). The proportion of estrogen receptor (ER) positive patients was 84%, and the proportion of patients with invasive tumors was 92% after deleting the missing data. The significant PCs among the first 20 in each study ( Supplementary Table 4 ) were included in the analyses of associations between SNPs and BC risk. Therefore, age and PCs were adjusted for as possible confounders in the multivariate logistic regression analysis. The workflow of the present study is shown in Figure 1 . In total, we used 7345 SNPs that passed QC in the analysis, including 442 genotyped SNPs and 6903 imputed SNPs. Multivariate logistic regression and meta-analysis results showed that there were 666 SNPs significantly associated with BC risk ( P <0.05), of which 101 SNPs remained significant after FDR correction < 0.05. The associations between SNPs of genes involved in the NER pathway and BC risk in the DRIVE study are shown in Figure 2 . These 101 top SNPs were mapped to BIVM-ERCC5 , GTF2H4, COPS2 , ELL and COPS4 ( Supplementary Table 5 ). COPS4 rs75870305 was deleted due to mapping or clustering errors on the NCBI website. Seven tagSNPs remained for additional analysis after removal of SNPs in high pairwise LD ( Supplementary Figure 1 ). To search for functional SNPs, we used stringent criteria with RegulomeDB scores ≤ 3 and with eQTL evidence in breast tissues or blood cells ( Supplementary Table 6 and Supplementary Table 7 ). As a result, BIVM-ERCC5 rs12870541 and ELL rs34664859 were left out, because of no functional annotation. Finally, the stepwise analysis kept four independent, potentially functional SNPs for further analysis ( Table 1 ). Supplementary Figure 2 presents the forest plots of the meta-analysis of the four independent SNPs. The results showed that SNP rs1323697 G>C and rs141308737 C>deletion were associated with a significantly increased risk of BC (OR = 1.06, 95% CI = 1.03–1.10, P = 2.66 × 10 −4 ; OR = 1.06, 95% CI = 1.03–1.09, P = 3.60 × 10 −4 , respectively), while two other SNPs were associated with a significantly decreased BC risk (rs1264308 C>T: OR = 0.93, 95% CI = 0.89–0.97, P = 2.21 × 10 –4 ; and rs1469412 T>C: OR = 0.93, 95% CI = 0.90–0.96, P = 3.09 × 10 −6 ). There was no heterogeneity observed for the effect estimates of these four SNPs from the 14 GWASs. Dropping any one of the studies in the DRIVE study, did not change the pooled ORs and their 95% CIs ( Supplementary Table 8 ). The association results from different genetic models for each SNP, including additive and dominant models, showed that all of the SNPs were significantly associated with BC risk in all of the genetic models ( Table 2 ). Although the SNP rs4808801 (in the chromosome region 19p13.11 where ELL is located) has been previously reported by a GWAS 31 , the ELL rs1469412 that we identified in the present study has a moderate LD with rs4808801 (r 2 =0.471). The effect size (beta) values of each SNPs are very similar, from 0.056 to 0.074, so we did not consider the weight of each SNP in our analysis of the combined risk genotypes. Using a dominant model, we combined risk genotypes of rs1323697 GC+CC, rs141308737 C-/--, rs1264308 CC and rs1469412 TT into a genetic score as the number of unfavourable genotypes (NUGs). The trend test indicated a significant association between an increased NUGs and an increased risk of BC ( P <0.0001, Table 3 ). Stratified analyses were performed to assess subgroups defined by age, ER status and invasiveness. We found that the risk associated with NUGs was more evident in the younger group (OR = 1.14, 95% CI = 1.08–1.20, P < 0.0001, Supplementary Table 9 ), but no heterogeneity or interaction were observed between these strata ( P = 0.227 and 0.274, respectively, Supplementary Table 9 ). Subgroups analysis (ER status and histological type) also showed similar results by age among patients with ER + and invasive tumors ( Supplementary Table 9 ). Additionally, we found no significant differences between ER + and ER − patients ( P = 0.990) or between invasive and in situ carcinomas ( P = 0.945). The in silico eQTL analysis among 373 European descendants with both SNP genotype and mRNA expression data showed that BIVM-ERCC5 rs1323697 C allele demonstrated a significant association with increased mRNA expression levels of BIVM in both additive ( P = 0.022) and dominant models ( P = 0.025) ( Figure 3a and 4b ). The ELL rs1469412 C allele also demonstrated a significant association with increased mRNA expression levels of ELL in all genetic models ( Figure 3f , 4g and 4h : P = 2.67E-22, 1.14E-17 and 3.01E-11, respectively). However, no significant associations between the other two SNPs and corresponding mRNA expression levels were found ( Figure 3d and 4e ). In addition, GTF2H4 rs114596632, the same SNP with rs1264308, has been reported significantly associated with a decreased mRNA expression levels in 270 lymphoblastoid cell lines from HapMap 32 . To further examine the correlation between the significant SNPs and mRNA expression levels, we searched GTEx as well and found that BIVM-ERCC5 rs1323697, GTF2H4 rs1264308 and ELL rs1469412 were correlated with their specific mRNA expression levels in the whole blood cells ( P =0.003, 0.032 and <0.0001, respectively), but COPS2 rs141308737 was unrelated to its gene expression levels. In addition, COPS2 rs141308737 ( P =0.026) and BIVM-ERCC5 rs1323697 ( P =0.001) had a positive correlation with their gene specific mRNA expression levels in breast tissues ( Supplementary Table 10 ).

To determine whether genetic variants in the NER pathway genes contribute to BC susceptibility, we performed association analyses between 7345 SNPs in 138 genes and BC risk with a large sample size of 28758 cases and 24349 controls. As a result, we identified four novel susceptibility variants, BIVM-ERCC5 rs1323697 at 13q33.1, GTF2H4 rs1264308 at 6p21.33, COPS2 rs141308737 at 15q21.2, and ELL rs1469412 at 19p13.11. In addition, the eQTL analysis results revealed that BIVM-ERCC5 rs1323697 C allele was associated with an increased mRNA expression levels, as was the ELL rs1469412 C allele. These results indicate that these two SNPs might influence mRNA expression levels and thus the functions of the genes, a possible mechanism underlying the observed associations. These findings suggest that variants in the NER pathway genes play an important role in the development of BC possibly by influencing mRNA expression. The NER pathway is a mechanism that recognizes and repairs bulky DNA damage caused by chemical compounds, environmental carcinogens, and exposure to UV-light. The repair of damaged DNA involves at least 30 polypeptides within two different sub-pathways of NER known as transcription-coupled repair (TC-NER) and global genome repair (GG-NER) 33 . The TCR and GGR processes are different in terms of damage recognition: RNA polymerase II (RNAP II) is needed in TC-NER, while XPC-hHR23B complexes together with XPE complex are needed in GG-NER. In general, genes of GG-NER have been associated with cancer predisposition 15 . However, the present study indicated some genes of TC-NER also might be involved in BC susceptibility, such as ELL, COPS2 and GTF2H4 , but their exact mechanisms involved in the BC etiology need to be further studied. BIVM-ERCC5 rs1323697 is located on 13q33.1, which has not been reported by any of the GWASs included in the present analysis. Based on the NCBI website ( https://www.ncbi.nlm.nih.gov/gene/100533467 ), this locus represents naturally occurring read-through transcription between the neighbouring BIVM (basic, immunoglobulin-like variable motif containing) and ERCC5 (excision repair cross-complementing rodent repair deficiency, complementation group 5) genes on chromosome 13. The read-through transcript encodes a fusion protein that shares sequence identity with the products of each individual gene ( Supplementary Figure 3 ). Because the present study mainly indicated that rs1323697 was correlated with BIVM gene expression levels, the discussion will focus on the function of BIVM only. Previous studies have shown that BIVM possesses virtually no sequence similar to any currently described protein, making the prediction of a function challenging 34 . It is highly likely that BIVM is essential for some aspect of basic cellular functioning and is expressed in a nearubiquitous manner 34 . The presence of a CpG island at the 5’ end of BIVM and its wide tissue distribution suggest that it may function as a housekeeping gene 34 , 35 . While others think it is likely that the immunoglobulin-like motif in BIVM may have functions similar to an immunoglobulin, but this remains to be experimentally confirmed 36 . Furthermore, we found that SNP rs1323697 is located at the LUN-1 motif, as shown by the position weight matrix (PWM) based Sequence Logo ( Supplementary Figure 4 and Supplementary Table 7 ). GTF2H4 , known as a general transcription factor IIH subunit 4, encodes a subunit of transcription factor IIH ( TFIIH ), a helicase that is responsible for unwinding DNA structure, allowing repair of the damaged DNA, and it is involved in both NER process and transcription control interacting with variable factors important in carcinogenesis 37 . The TFIIH complex has both ATPase and helicase activities and opens DNA at sites of DNA distorting damage, and the TFIIH4 subunit may regulate the ATPase activity of the TFIIH subunit (XPB, a protein coded by ERCC3) 38 . Previous studies have found that some GTF2H4 SNPs were significantly associated with lung cancer risk and survival, multiple sclerosis risk and cervical cancer 32 , 39 – 41 , but there is no report on the associations between GTF2H4 SNPs and BC risk to date. There is an interesting finding that BC and lung cancer risk was associated with the same SNP, GTF2H4 rs1264308 32 . GWAS catalog results indicated that some of the adjacent genes shared the same location 6p21.33, including ABCF1 , PPP1P18 and LOC105375013 , have a high LD with rs1264308 ( Supplementary Table 11 ). As an intron SNP, GTF2H4 rs1264308 may have an effect on the disease by changing motif FOXJ3 ( Supplementary Figure 4 ) or by mechanisms of interacting with other genes above-mentioned. However, none of the other genes have a known function in NER. In the present study, COPS2 rs141308737 has no functional clues from mRNA expression levels. However, a study showed that over-expression of COPS2 was linked to chromosome instability 42 . Functional prediction software shows that rs141308737 is located at the ER motif ( Supplementary Figure 4 ) and can bind to the CJUN protein ( Supplementary Table 7 ). It has been reported that endogenous c-Jun plays a key role in ErbB2-induced migration and invasion of mammary epithelial cells and mediates the expansion of a self-renewing population of mammary tumor stem cells via the production of CCL5 and SCF to enhance BC tumor invasiveness 43 . As for ELL rs1469412, although some SNPs in this region have been reported by GWASs, it is necessary to include this SNP, because it was only in moderate-to-low LD with other reported SNPs ( Supplementary Table 11 ) and lack of functional analysis in the previously published study. ELL is known as an elongation factor for RNA polymerase II, which is an important gene in the TC-NER sub-pathway. One study reported that ELL encoded an elongation factor that could increase the catalytic rate of RNA polymerase II transcription by suppressing transient pausing by polymerase at multiple sites along the DNA 44 . Another study showed that ELL was a key regulator of transcriptional elongation, suggesting that, as an E3 ubiquitin ligase for c-Myc and a potential tumor suppressor, ELL may function as a partner of steroid receptors, hypoxiainducible factor 1-alpha (HIF-1α), E2F1 and the TFIIH complex, modulating their binding partner’s activity 45 . The present study showed that the ELL rs1469412 C allele was associated with an increase in mRNA expression levels, exerting a protective effect on BC risk. However, further studies are needed to investigate biological mechanisms underlying the observed associations between ELL rs1469412 and BC risk. It should also be mentioned that the present study has some limitations. Firstly, due to the limited access to phenotypes of the published GWAS datasets with many PCs included in the analysis, we could not adjust for some known risk factors, such as smoking, menstrual, reproductive and lactational history 46 , and the findings need to be verified in other BC studies with more detailed information about the known risk factors. Secondly, we did not have access to the target tissues collected by the participating GWAS studies, and we only did in silico analysis using published data for the functional prediction of the identified SNPs. Therefore, the biological mechanisms by which the four SNPs of the genes may influence BC risk remain unclear. Thirdly, the study populations were of non-Hispanic whites, and thus the findings may not generalizable to other ethnic groups, and thus additional studies in other ethnic groups are warranted. In conclusion, this large-scale meta-analysis of 14 published GWASs among 53107 subjects of European descent identified four novel BC susceptibility loci in the NER pathway genes (i.e., BIVM-ERCC5 rs1323697, GTF2H4 rs1264308, COPS2 rs141308737 and ELL rs1469412) and also provided some evidence for their functional relevance. Further studies on the exact biological mechanisms and functional analysis of these SNPs in the BC etiology are needed.

This meta-analysis included a sub-dataset of SNPs in the NER pathway genes from each of 14 previously published BC GWASs for a total of 28758 BC cases and 24349 controls of European ancestry from the DRIVE study (phs001265.v1.p1), which is different from the DRIVE-Genome-Wide Association meta-analysis (phs001263.v1.p1) previously used by others 16 ( Supplementary Table 1 ). The DRIVE study (phs001265.v1.p1), which included 17 GWASs, was one of the five projects funded in 2010 as part of the NCI’s Genetic Associations and Mechanisms in Oncology (GAME-ON) initiative. For this meta-analysis, we excluded three studies including the “Women of African Ancestry Breast Cancer Study (WAABCS)”, which is a study of African ancestry, and “The Sister Study (SISTER)” and “The Two Sister Study (2 SISTER)”, which had a different study design that used cases’ sisters as the controls. These 14 GWAS studies consist of Breast Oncology Galicia Network (BREOGAN); Copenhagen General Population Study (CGPS); Cancer Prevention Study-II Nutrition Cohort (CPSII); European Prospective Investigation Into Cancer and Nutrition (EPIC); Melbourne Collaborative Cohort Study (MCCS); Multiethnic Cohort (MEC); Nashville Breast Health Study (NBHS); Nurses’ Health Study (NHS); Nurses’ Health Study 2 (NHS2); NCI Polish Breast Cancer Study (PBCS); The Prostate, Lung, Colorectal and Ovarian Cancer Screening Trial (PLCO); Study of Epidemiology and Risk factors in Cancer Heredity (SEARCH); Swedish Mammography Cohort (SMC); and Women’s Health Initiative (WHI). The details of case and control recruitment and their characteristics are summarized in Supplementary Table 2 . For all of the GWAS datasets, Illumina Infinium OncoArray-500k BeadChip genotyping platforms were used, and only two main de-identified variables (sex and age at interview) were available to us. For the cases, other three de-identified variables (age at diagnosis, estrogen receptor status, and histology type) were available. Each of the 14 studies was reviewed and approved by the corresponding Institutional Review Board and thus exempted by Duke Institutional Review Board. Candidate genes in the NER pathway were selected according to the online datasets “Molecular Signatures Database v6.1 (MsigDB)” ( http://software.broadinstitute.org/gsea/msigdb/search.jsp ) and “PathCards” ( http://pathcards.genecards.org/ ) using the key words “nucleotide excision repair”. In total, we selected 138 candidate genes from eight NER-related pathways after excluding duplicate genes, pseudo genes and withdrawn genes (LOC652672 and LOC652857) in the National Center for Biotechnology Information (NCBI). The detailed genes selection results are listed in Supplementary Table 3 . To avoid poor quality markers to be included in the imputation, we performed stringent quality control before imputation by including the following criteria: the minor allelic frequency (MAF) ≥1%, genotyping rate ≥95%, missing rate ≤90%, Hardy-Weinberg equilibrium (HWE) ≥1×10 −6 , All SNPs were flipped to forward strand and aligned with the reference genome data. Ambiguous SNPs with A-T or G-C alleles that are hard to determine the strand orientation by allele frequency were removed. According to the multi-population reference panels from the 1000 Genomes Project Phase 3, SNPs within the aforementioned 138 candidate genes and their ± 500 kb flanking regions were also extracted, and imputation for each study was performed using IMPUTE2 software 25 . Imputed SNPs within 2-kb up- and down-stream of each gene’s region were extracted for further analysis. After imputation, SNPs that met the following quality control criteria were included in further analysis: imputation SNPs with information score ≥0.80 in IMPUTE2; a minor allele frequency (MAF) ≥5%; and a P value for the Hardy–Weinberg Equilibrium test ≥10 −6 . Due to differences between the 14 studies, 8433 to 9016 common SNPs remained in each study for further analysis. The final analysis included 7345 SNPs that were common to all 14 studies. To investigate the functions of candidate SNPs, we searched for functional annotation of the SNPs in three online functional prediction website: RegulomeDB ( http://regulomedb.org/ ), HaploReg ( http://archive.broadinstitute.org/mammals/haploreg/haploreg.php ) and SNPinfo ( https://snpinfo.niehs.nih.gov/snpinfo/snpfunc.html ). In addition, we performed the expression quantitative trait loci (eQTL) analysis by using data from multiple sources: lymphoblastoid cellline data of 373 subjects from the European Variation in Health and Disease Study (GEUVADIS) and the 1000 Genomes Project (phase I integrated release 3, March 2012) 26 . Furthermore, we used Genotype-Tissue Expression project (GTEx) results to obtain the corresponding mRNA expression in whole blood and breast tissues ( https://gtexportal.org/home/) 27 . For each study and the combined dataset, principal components (PCs) were calculated using the Genome-wide Complex Trait Analysis (GCTA) on the LD-pruned subset of the whole-genome-typed dataset 28 . The top 20 PCs were assessed for their associations with BC risk using univariate logistic regression analysis. Those PCs with significant associations in each study were included as covariates in further analyses of associations between SNPs and BC risk. For each SNP, we estimated odds ratios (ORs) and 95% confidence intervals (CIs) by unconditional logistic regression of case/control groups with adjustment for age and PCs. We performed the meta-analysis by using the inverse variance method to combine the results of the 14 studies. We defined heterogeneity as a Cochran’s Q test P ≤ 0.10 or I 2 >50.0%. We used fixed-effects models, if no heterogeneity existed among the 14 studies, and random effects models were used, when heterogeneity existed. To assess the robustness of the results, we performed a sensitivity analysis by omitting each study one by one 29 . The false discovery rate (FDR) with a critical cut-off value of 0.05 using the linear step-up method of Benjamini and Hochberg was mainly used to correct for multiple comparisons to reduce the probability of false-positive findings 30 . To observe the combined effect of significant SNPs, we used the number of unfavorable genotypes (NUGs) of the significant SNPs as a genetic score to assess classification performance of the model. According to the frequency of each group and the effect values, we also dichotomized all the individuals into a low-risk group (0–2 NUGs) and a high-risk group (3–4 NUGs). In the eQTL analysis, we calculated the correlations between SNPs and specific mRNA expression levels by using a general linear regression model. Statistical analyses were performed using PLINK (version 1.9), SAS (version 9.3; SAS Institute, Cary, NC, USA) and R (version 3.0.2). The Manhattan plots and linkage disequilibrium (LD) plots were generated by Haploview v4.2, and regional association plots were constructed by LocusZoom ( http://locuszoom.sph.umich.edu/locuszoom/ ).

Breast cancer (BC) is the most frequently diagnosed cancer and the leading cause of cancer deaths among women worldwide, with an estimated 1.7 million cases and 521 900 deaths in 2012, accounting for 25% of all cancer cases and 15% of all cancer deaths among women 1 . Despite the declining mortality rate due to early screenings and advanced medical therapies, the incidence rate of BC has remained steady over the past two decades in the US ( https://seer.cancer.gov/statfacts/html/breast.html ). Therefore, it is necessary to identify additional genetic factors that can be used for defining susceptible individuals at risk for BC. Although the mechanisms of breast carcinogenesis are still not fully understood, a variety of risk factors have already been identified 2 – 5 . Some studies have shown that mammalian cells can convert estrogen into related compounds that not only generate free radicals capable of damaging DNA but also bind to DNA, causing the loss of a nucleotide base, a process known as depurination. The resulting mutations can convert a normal cell into a cancerous one 6 – 8 . Another putative risk factor is smoking. Although there are no consistent results about the association between smoking and BC risk, there are carcinogens in tobacco smoke such as polycyclic aromatic hydrocarbons (PAH), aromatic amines, and nitrosoamines, and these carcinogens might cause DNA damage and adduct formation in mammary epithelial cells 9 , 10 . In addition, many epidemiologic studies reported a positive association between BC risk and alcohol consumption. Animal models of BC, although not entirely consistent, do provide the support for an enhancing action of ethanol on mammary carcinogenesis 11 . Overall, evidence from human studies, animal studies and cell culture experiments support some biologically plausible mechanisms, such as an increase in circulating estrogens and androgens, enhancement of mammary gland susceptibility to carcinogenesis, increased mammary carcinogen-induced DNA damage, and a greater potential for invasiveness of BC cells 11 , 12 . These mechanisms are all likely involved in DNA damage leading to the initiation of mutations and carcinogenesis, and thus the DNA repair system plays a critical role in protecting against mutations, maintaining genomic integrity and preventing carcinogenesis of the breasts 13 . One of the DNA repair pathways is nucleotide excision repair (NER), a highly versatile and sophisticated DNA damage removal mechanism that counteracts the deleterious effects of a multitude of DNA lesions, including major types of damage induced by environmental mutagens and carcinogens. The most relevant lesions to be repaired by NER are cyclobutane pyrimidine dimers (CPDs) and 6–4 photoproducts (6–4PPs) produced by the shortwave UV component of sunlight. In addition, numerous bulky chemical adducts are eliminated by this repair process as well 14 , 15 . Given the importance of NER in the repair of UV-induced DNA damage, it seems that the NER pathway may not be relevant to BC risk, because there is no evidence that UV light may cause BC; however, it is likely that tobacco smoke may cause DNA damage in breast tissues. A recent large study with pathway-level analysis using hierarchical modelling across five cancers, including 11 DRIVE (the Discovery, Biology, and Risk of Inherited Variants in Breast Cancer) GWAS datasets of 33832 BC study subjects of European descent, did not find any specific risk-associated SNPs in genes involved in NER, but the limited study power did not allow the investigators to find an association with the overall genetic variation of the NER pathway 16 . Other prior studies also investigated associations between SNPs in DNA repair pathway genes and BC risk, but these studies had relatively small sample sizes without a focus on the NER pathway, although they have some notable findings, such as XRCC3 and ERCC4 . 17 – 24 . Therefore, we hypothesize that genetic variants in the NER pathway genes are associated with BC risk. To assess the role of functional SNPs of the NER pathway genes in the BC etiology, we performed a much larger meta-analysis of 14 previously published DRIVE GWAS datasets with 53107 study subjects of European descent. In contrast to the previously published studies, the present analysis had a much larger sample size to focus on functional SNPs in the NER pathway genes. Hence, using a hypothesis-driven pathway-based approach with a much increased study power, we expected to identify some susceptibility loci in the NER pathway genes that have biologically relevant functions and thus play a role in the BC etiology.