Association of genetic variants of FBXO32 and FOXO6 in the FOXO pathway with breast cancer risk.

The flowchart for the analysis is shown in Figure 1 . The results of the top 20 PCs of the datasets are shown in Table S3 . Because of the differences in genotyping platforms used by the 14 studies, there were a range between 6,163 and 6,429 SNPs in each individual study for further analysis. In total, the final meta-analysis of the 14 studies included 5,214 SNPs that passed the quality control, including 365 genotyped SNPs and 4,849 imputed SNPs. The distribution of information score for those selected SNPs in each study is shown in Figure S1 . The meta-analysis showed that 338 SNPs were significantly associated with BC risk in an additive genetic model ( P < 0.05), and after multiple testing corrections by BFDP, five SNPs were still statistically noteworthy with BFDP < 0.8 ( Figure 2 ). To further identify SNPs as independent predictors of BC risk, stepwise logistic regression analyses were performed to evaluate the independent effects of the five significant SNPs on BC risk with adjustment for age, significant PCs and another four previously published risk-associated SNPs in the same DRIVE study 20 . As a result, three independent SNPs (i.e., FBXO32 rs10093411 A>G, FOXO6 rs61229336 C>T and FBXO32 rs62521280 C>T) remained statistically significant in association with BC risk ( P = 0.0008, 0.0011 and 0.0017, respectively), which were then used for further analyses ( Table 1 ). The types of the three SNPs and their associations with BC risk are presented in Table S4 , of which SNP rs10093411 was genotyped, and the other two were imputed. There was no heterogeneity existed among the 14 GWASs for the effects of these independent SNPs. The forest plots of the three SNPs by the meta-analysis are summarized in Figure S2 . The results showed that two SNPs in the same gene were associated with a significantly increased risk of BC ( FBXO32 rs10093411 A>G: OR = 1.06, 95% CI = 1.03–1.09, and P = 1.55 × 10 –4 ; FBXO32 rs62521280 C>T: OR = 1.07, 95% CI = 1.03–1.11, and P = 8.35 × 10 –4 ), while the other SNP was associated with a significantly decreased BC risk ( FOXO6 rs61229336 C>T: OR = 0.96, 95% CI = 0.94–0.99, and P = 5.53 × 10 –3 ). Regional association plots of these three independent SNPs in the 200 kb up- and down-stream regions are presented in Figure S3 . As shown in Table 2 , the effects of the FBXO32 rs10093411 G, FOXO6 rs61229336 T and FBXO32 rs62521280 T alleles on BC risk were statistically significant (the trend test in multivariate analysis: P = 0.0002, 0.001 and 0.0005, respectively). We subsequently combined risk genotypes of FBXO32 rs10093411 AG+GG, FOXO6 rs61229336 CC and FBXO32 rs62521280 CT+TT into a genetic risk score as the NRG to assess the joint effect of the genotypes of these three independent SNPs on BC risk. All the participants were divided into four groups of zero, one, two and three risk genotypes. The trend test indicated that the increased NRG was significantly associated with an increased BC risk ( P < 0.0001, Table 2 ). According to the effect values and the frequency of each group, we further dichotomized all the participants into two groups: low-risk (0–1 NRG) and high-risk (2–3 NRG), and we found that the risk associated with NRG was more evident in the high-risk group (in multivariate analysis: OR = 1.11, 95% CI = 1.07–1.15, P < 0.0001, Table 2 ). Further stratified analyses by subgroups of age, ER status and invasiveness showed that risk associated with NRG was more evident in the subgroup of age ≤60 (in multivariate analysis: OR = 1.13, 95% CI = 1.07–1.19, P 0.05) ( Table 3 ). Functional prediction by RegulomeDB showed that FBXO32 rs10093411 A>G, FOXO6 rs61229336 C>T and FBXO32 rs62521280 C>T had a RegulomeDB score of 3a, 7 and 5, respectively and that these SNPs may be located at transcription factor binding sites or DNase I regulating sites. We also searched for SNPs in high LD (r 2 ≥ 0.8) with these three independent SNPs, and their functional prediction was made as well by using HaploReg. The results suggest that FBXO32 rs10093411 A>G is located in an intron and may change the motifs of MZF1; FOXO6 rs61229336 C>T is also located in an intron with the selected eQTL for three hits and may change the motifs of CACD, PU.1 and Zbtb3, which may be markers of promoter histone in THYM and enhancer histone in 9 tissues; whereas FBXO32 rs62521280 C>T is located in the potential enhancer region with histone methylation in 5 tissues ( Table S4 ). We further assessed potential functions of these three independent SNPs by using data from the ENCODE Project, which suggests that FOXO6 rs61229336 C>T is located in DNase I hypersensitive sites and has considerable levels of H3K4Me1 acetylation ( Figure S4 ). We further explored correlations between genotypes of the three independent SNPs and their corresponding mRNA expression levels in the publically available RNA-seq data of lymphoblastoid cell lines generated from 373 European descendants in the 1000 Genomes Project. Unfortunately, no significant results were found ( Figure S5 ). We further performed the correlation analysis by using data from the GTEx Project, and the results showed that rs10093411 A>G was not significantly associated with FBXO32 mRNA expression levels in breast tissue ( P = 0.100, Figure 3a ) and the whole blood ( P = 0.260, Figure 3b ); however, the rs62521280 C>T was significantly associated with decreased levels of mRNA expression of FBXO32 in breast tissues ( P = 0.030, Figure 3c ) but not for the whole blood ( P = 0.540, Figure 3d ). Although rs61229336 C>T was not significantly associated with levels of mRNA expression of FOXO6 in breast tissue ( P = 0.180, Figure 3e ), it was significantly associated with decreased levels of mRNA expression of FOXO6 in the whole blood ( P = 8.70 × 10 −4 , Figure 3f ).

We performed a case-control meta-analysis of the participants from 14 out of the 17 previously published BC GWASs from the DRIVE study (phs001265.v1.p1), which is different from the DRIVE-Genome-Wide Association meta-analysis (phs001263.v1.p1) previously used by other researchers, and the differences between the two datasets have been described in detail elsewhere 20 . In brief, the major differences are as follows: 1) different sources of studies and 2) different data types: The DRIVE study we used has genotyping data with detailed SNP information, but the DIRVE study previously used by other researchers had the summary data only good for meta-analyses. The DRIVE study used in the present study was one of five projects funded in 2010 by the National Cancer Institute-supported Genetic Associations and Mechanisms in Oncology (GAME-ON). Among the 17 studies, “Women of African Ancestry Breast Cancer Study (WAABCS)” was an African ancestry study and “The Sister Study (SISTER)” and “the Two Sister Study (2 SISTER)” had different research designs from others and used cases’ sisters as controls; therefore, our meta-analysis excluded these three studies. As a result, a total of 28,758 BC cases and 24,349 controls in 14 GWAS studies of European ancestry participants were included in the final analysis, whose characteristics are summarized in Table S1 . These 14 DRIVE GWASs included: 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). Illumina Infinium OncoArray-500k BeadChip genotyping platforms were used for all of the GWAS datasets, and only sex and age at interview for all the participants were available to us, while three other variables including age at diagnosis, estrogen receptor status, and tumor histology type were available to the cases; for age variables, we adopted the age at interview for the controls and the age at diagnosis for the cases. Each of the original studies approved by the Institutional Review Boards of the Participating institutions received written informed consent from the participants. We selected candidate genes in the FOXO pathway according to the databases of KEGG, BIOCARTA, REACTOME, Canonical pathways and Gene Ontology (GO) in the “Molecular Signatures Database v7.0 (MsigDB)” ( http://software.broadinstitute.org/gsea/msigdb/search.jsp ) used by the keyword “FOXO”. In total, we identified 55 candidate genes after excluding 16 duplicate genes ( Table S2 ). We performed quality control before imputation to avoid poor quality markers to be included with the following stringent criteria: (1) the minor allelic frequency (MAF) ≥ 1%, (2) missing rate ≤ 10%, (3) genotyping success rate ≥ 95%, and (4) Hardy-Weinberg equilibrium (HWE) P ≥ 1×10 −6 . SNPs located in the aforementioned 55 candidate genes and their ± 500 kb flanking regions were extracted by using IMPUTE2 software with the reference panel from the 1000 Genomes Project data (phase 3) 21 . After quality control, imputed SNPs within 2-kb up- and down-stream of genes in the FOXO pathway were extracted for further analysis. SNPs for the final meta-analysis were selected with the following quality control criteria: (1) a genotyping call rate ≥ 95%; (2) Imputed SNPs with an information score ≥ 0.80; (3) MAF ≥ 5%; and (4) HWE P ≥ 10 –6 . Principal components (PCs) were calculated for each of the GWASs, and their combined dataset was evaluated by using the Genome-wide Complex Trait Analysis 22 . The associations between the top 20 PCs and BC risk were evaluated by using univariate logistic regression analysis. Significant PCs together with age as covariates in the final model were adjusted for in further SNP association analysis. We calculated odds ratios (ORs) and 95% confidence intervals (CIs) for each SNP by unconditional logistic regression with adjustment for covariates (age and significant PCs). The four previously published SNPs from the same DRIVE study 20 were also adjusted for in a logistic regression model to identify additional significant SNPs. A meta-analysis was further performed with the inverse variance method by combining the results of a log-additive model of the 14 studies. If the Cochran’s Q test P -value ≤ 0.1 or I 2 > 50%, a random-effects model was used; otherwise a fixed-effects model was employed. The results were first corrected by false discovery rate for multiple testing correction. Because many SNPs under investigation were a high linkage disequilibrium (LD) as a result of imputation, the Bayesian false-discovery probability (BFDP) approach was also used for multiple test correction in substitution for the false discovery rate, with a cut-off value of 0.8 as recommended 23 . We used a prior probability of 0.01 to detect an upper bound of 3.0 for an association with variant genotypes or minor alleles of the SNPs. The number of risk genotypes (NRGs) of the independent SNPs was counted as a genetic score and subsequently used to evaluate combined effects of the SNPs. Additionally, Manhattan and LD plots were constructed by Haploview v4.2 24 , and regional association plots for independent SNPs were produced by LocusZoom 25 . Other statistical analyses were performed with SAS software Version 9.4 (SAS Institute, Cary, NC), R (version 3.5.0) and PLINK (version 1.90), if not specified otherwise. Finally, we predicted potential functions of the identified independent SNPs in online functional prediction website: RegulomeDB ( http://www.regulomedb.org/ ) and HaploReg 26 ( http://archive.broadinstitute.org/mammals/haploreg/haploreg.php ). In addition, we performed an expression quantitative trait locus (eQTL) analysis to assess the associations between the SNPs and mRNA expression levels of their corresponding genes by using the genotyping and expression data from the lymphoblastic cell lines of 373 European descendants available in the 1000 Genomes Project 27 and the Genotype-Tissue Expression (GTEx) project v8.p2 database ( https://gtexportal.org/home/ ) 28 , which included genomic data from both whole blood and breast tissues. The samples used for the eQTL analysis are publically available and not from the DRIVE study.

Some published studies have used the data from different GWAS studies to identify the associations of different pathway genes and BC risk 1 , 4 , 29 , 30 , but few have investigated the association between FOXO pathway genes and BC risk. One study found that AKT1 rs3803304 may be a predictive biomarker for BC risk in the high-altitude Ecuadorian mestizo population 18 . In another Egyptian population, SIRT1 rs3758391 and rs12778366 polymorphisms were found to be associated with BC risk and prognosis 19 . In the present study, we have investigated whether SNPs of FOXO pathway-related genes are associated with BC risk by using available genotyping data from the previously published 14 GWASs in the DRIVE study, and the study population was different from the above-mentioned two studies. We identified three independent SNPs in two genes (i.e., FBXO32 rs10093411 and rs62521280 at 8q24.13, and FOXO6 rs61229336 at 1p34.2) to be associated with BC risk. There was a significant effect of the combined genotypes of these three SNPs on BC risk. In stratified analyses, all the results with P inter > 0.05 suggested that there was no interaction between age and NRG in all subgroups. Further eQTL analyses showed that the FBXO32 rs62521280 T allele was associated with decreased mRNA expression levels in breast tissues, while the FOXO6 rs61229336 T allele was found to be associated with decreased mRNA expression levels in the whole blood cells. These eQTL results provide some support for biological plausibility of the observed associations. FBXO32 , located on chromosome 8q24.13, encodes an F-box only protein 32 (FBXO32), an E3 ubiquitin ligase that is essential for hallmark phenotypic changes and gene expression underlying epithelial-mesenchymal transition as well as involved in the process of tumorigenesis 31 – 34 . Several studies have indicated that FBXO32 may be a functional tumor suppressor in some cancer types via promoter methylation 35 – 38 . Even with the known importance of FBXO32 in cancers, few reported studies have investigated the roles of FBXO32 in BC risk. One study found that an increased expression of FBXO32 facilitated apoptosis of multiple breast cancer cell lines 39 , and another study showed that FBXO32 deficiency in breast cancer cells led to the accumulation of Krüppel-like factor 4 and promoted tumorigenesis 40 . A more recent GWAS study identified a locus in RNTSKP155 that shares the same location of 8q24.13 as FBXO32 41 , but the previously reported rs58847541 SNP in RNTSKP155 was not in LD with rs10093411 and rs62521280 ( Figure S6a ). Interestingly, we found in the present study that the rs62521280 T allele might down-regulate the expression of FBXO32 , suggesting that FBXO32 may play a protective role in BC risk, but additional experimental studies are needed to explain how FBXO32 rs62521280 C>T influences BC risk. FOXO6 , located on chromosome 1p34.2, encodes a forkhead box protein O6 (FOXO6) 42 , but little is known about the roles of FOXO6 in tumorigenesis, especially in BC. One study showed that FOXO6 was overexpressed in both colorectal cancer cell lines and tumor tissues and that FOXO6 knockdown inhibited the migration and invasion of colorectal cancer cells 43 ; another study showed that overexpression of FOXO6 promoted gastric cancer cell tumorigenicity through regulation of C-myc expression 44 ; likewise, FOXO6 was also shown to be highly overexpressed in both breast cell lines and tumor tissues 45 , indicating that FOXO6 may act as an oncogene in BC. Although one GWAS study identified a locus on HIVEP3 that shares the same location 1p34.2 as FOXO6 41 , FOXO6 rs61229336 is not in LD with the previously published HIVEP3 rs79724016 ( Figure S6b ) without functional analysis. In the present study, we found that rs61229336 T allele might down-regulate expression of the likely oncogenic FOXO6 , leading to a reduced BC risk. Nevertheless, we did not have additional experimental data that could explain how FOXO6 rs61229336 C>T influenced BC risk. It should be pointed out that there are several limitations in the present study. First of all, the small sample sizes in some of DRIVE GWAS studies, such as WAABCS with African Americans, were not sufficient for stratification analysis for ethnic difference in genetic background and therefore excluded from the analysis. Thus, the genotype data available for the analysis consisted of non-Hispanic whites only; therefore, the results may not be generalizable to the general population. Second, several known risk factors, such as physical activities, unhealthy lifestyle and reproductive factors, prevented us from performing complete adjustment in the analyses. Finally, the publically available data for the eQTL analysis was very limited, and gene expression analysis of biological samples from the participants in the DRIVE study was not performed, because relevant data from the samples of target breast tissues were not available. In summary, we analyzed the associations between genetic variants in 55 FOXO pathway-related genes and BC risk by using genotyping data from 14 previously published GWASs of 53,107 participants of European descendants in the DRIVE study. We identified three independent BC susceptibility loci in FOXO pathway genes (i.e., FBXO32 rs10093411 and rs62521280 at 8q24.13, and FOXO6 rs61229336 at 1p34.2), which may help with the identification of individuals at high-risk of developing BC, once these SNPs are further validated by other studies. Our findings offer some new evidence for genetic variants in the FOXO pathway as possible biomarkers that may provide additional insights into molecular biological mechanisms underlying the observed associations with BC risk.

Worldwide data show a continuing increase in the incidence of breast cancer (BC) 1 that has become the second leading cause of cancer deaths among women in the United States 2 . In 2020, about 279,100 cases were diagnosed with and 42,690 died from BC in the United States 3 . Therefore, additional biomarkers are needed to identify individuals who are at a greater risk of BC for early detection and prevention to reduce the incidence of BC. Several risk factors (e.g. physical activities, unhealthy lifestyle and reproductive factors) are known to contribute to BC risk, and gene mutations only explain approximately 9% ~ 13% of the heritability of BC 4 , 5 . It is well recognized that single nucleotide polymorphisms (SNPs), the most common form of genetic variation, are also contributed to BC risk, suggesting their importance as a molecular biological mechanism of carcinogenesis. In the post genome-wide association study (GWAS) era, it is possible to take more sophisticated analyses to identify cancer risk-associated functional SNPs in a biological pathway manner. With such a targeted pathway-based and hypothesis-driven approach, investigators may identify cancer risk-associated SNPs from previously published GWAS datasets by using available genotyping data with further evaluation of their biological functions. It is known that forkhead box class O (FOXO) transcription factors are important regulators of gene expression and play a pivotal role in regulating a variety of biological processes, including organismal development, cell signalling, cell metabolism and tumorigenesis 6 – 12 . FOXOs are involved in the regulation of the upstream signalling pathway of PD-L1 by transcriptional or post-translational manner. For example, FOXOs may inhibit PD-L1 expression in cancer cells and indirectly upregulate T cell response 13 . FOXOs are considered putative tumor suppressors, because the activation of FOXOs inhibits cell cycle and induces apoptosis in a variety types of tumor cells 14 – 17 . However, the role of FOXOs in carcinogenesis remains to be determined, because only few studies have investigated the effect of genetic variation in FOXO pathway genes on BC risk, such as an association of AKT1 gene mutation with breast cancer risk in the high-altitude Ecuadorian mestizo population 18 and the relationship between SIRT1 gene polymorphisms and breast cancer in Egyptians 19 . Considering the importance of the FOXO pathway in the biology of carcinogenesis, it is very likely that genetic variants in FOXO pathway genes are associated with BC risk. To test this hypothesis, we conducted a targeted pathway-based and hypothesis-driven approach to identify SNPs of FOXO pathway-related genes and examined their associations with BC risk by using available genotyping data from 14 previously published GWAS datasets of 53,107 BC case-control study participants in the DRIVE study.