A colorectal cancer genome-wide association study in a Spanish cohort identifies two variants associated with colorectal cancer risk at 1p33 and 8p12

Abstract

Background

Colorectal cancer (CRC) is a disease of complex aetiology, with much of the expected inherited risk being due to several common low risk variants. Genome-Wide Association Studies (GWAS) have identified 20 CRC risk variants. Nevertheless, these have only been able to explain part of the missing heritability. Moreover, these signals have only been inspected in populations of Northern European origin.

Results

Thus, we followed the same approach in a Spanish cohort of 881 cases and 667 controls. Sixty-four variants at 24 loci were found to be associated with CRC at p-values <10-5. We therefore evaluated the 24 loci in another Spanish replication cohort (1481 cases and 1850 controls). Two of these SNPs, rs12080929 at 1p33 (P_replication=0.042; P_pooled=5.523x10^-03; OR (CI95%)=0.866(0.782-0.959)) and rs11987193 at 8p12 (P_replication=0.039; P_pooled=6.985x10^-5; OR (CI95%)=0.786(0.705-0.878)) were replicated in the second Phase, although they did not reach genome-wide statistical significance.

Conclusions

We have performed the first CRC GWAS in a Southern European population and by these means we were able to identify two new susceptibility variants at 1p33 and 8p12 loci. These two SNPs are located near the SLC5A9 and DUSP4 loci, respectively, which could be good functional candidates for the association signals. We therefore believe that these two markers constitute good candidates for CRC susceptibility loci and should be further evaluated in other larger datasets. Moreover, we highlight that were these two SNPs true susceptibility variants, they would constitute a decrease in the CRC missing heritability fraction.

Background

Even though genetic susceptibility is thought to be responsible for almost 35% of the variation in colorectal cancer (CRC) risk [1], high penetrance mutations in Mendelian predisposition genes, such as APC, the mismatch repair (MMR) genes, or MUTYH are only able to explain <5% of cases [2]. The recent advances in the field of genetic epidemiology have validated the hypothesis that at least part of that remaining missing susceptibility lies in the form of multiple common low-risk variants, each conferring a modest effect on disease risk.

Genome-wide association studies (GWAS) are one of the most widespread methodologies for the detection of such susceptibility loci. The procedure (in distinction to candidate-gene association studies) offers an untargeted strategy for the detection of new low-penetrance variants, for it does not assume any a priori hypothesis on the location of these loci. This advantage has been proved important, since so far this kind of survey has successfully identified 20 variants at 8q24.21, 8q23.3, 10p14, 11q23, 15q13.3, 18q21.1, 14q22.2, 16q22.1, 19q13.1, 20p12.3, 1q41, 3q26.2, 12q13.3, 20q13.33, 6p21, 11q13.3 and Xp22.2 [3–6]. The combined effect of these variants altogether is thought to explain ~7% of the familial cancer risk [7]. Still, there is a high proportion of genetic contribution to CRC risk that has not been identified.

In this study we have undertaken a new screen for CRC susceptibility variants using a GWAS approach on our cohort of 881 CRC cases and 667 controls from the Spanish population. The use of a Southern-European dataset is a novelty that could lead to the identification of new candidate loci, since all of the populations where GWAS analyses have been conducted so far have been of Northern European origin. Although this may provide additional confirmation of the relationship of the 20 described variants to CRC risk in Southern Europe, we must also consider the possibility that there may be differences, at these or other particular loci in the genome, between these populations.

Materials and methods

Study populations

Subjects in Phase I were 882 cases and 473 controls ascertained through the EPICOLON II Project and 194 additional controls from the Spanish National DNA bank. The EPICOLON Consortium comprises a prospective, multicentre and population-based epidemiology survey of the incidence and features of CRC in the Spanish population [8, 9]. Cases were selected as patients with de novo histologically confirmed diagnosis of colorectal adenocarcinoma. Patients with familial adenomatous polyposis, Lynch syndrome or inflammatory bowel disease-related CRC, and cases where patients or family refused to participate in the study were excluded.

Mean age for cases in Phase I was 71.20 years (SD±0.70). Hospital-based controls were recruited through the blood collection unit of each hospital, together with cases. All of the controls were confirmed to have no history of cancer or other neoplasm and no reported family history of CRC. Controls were randomly selected and matched with cases for hospital, sex and age (± 5 years). Population controls from the National DNA bank were also genotyped, to lessen the deficit of controls. They were matched for sex, age (± 10 years) and geographical origin of the sample with the remaining cases. Both cases and controls were of European ancestry and from Spain (stated, when possible, as all four grandparents being Spanish).

Samples in Phase II consisted of 1436 CRC patients and 1780 controls: samples from Hospital Sant Pau were 125 CRC patients from a previously described cohort [10]; the Hospital Gregorio Marañón dataset consisted of 104 CRC patients participating in a pharmagogenetic survey; Catalan Institute of Oncology (ICO) samples were 439 patients who belonged to the Bellvitge Colorectal Cancer Study; the CHUS hospital in Santiago de Compostela was a subsample of 153 participants included in a pharmacogenetic study; 105 CRC cases and 1330 controls came from the Spanish National DNA bank; and 510 CRC cases and 450 DNA controls belonged to the EPICOLON I Project [8]. Of the cases, 60.4% were male and 39.6% female. Controls were matched for gender. Age mean was 69.61 (SD±0.59) for cases and 52.00 (SD±0.58) for controls. Gender and hospital distribution of samples for case and control groups on both Phases is shown on (Additional file 1: Table S1).

DNA was obtained from frozen peripheral blood by standard extraction procedures for all samples. Cases and controls were extracted in mixed batches to avoid bias.

Ethical standards

The study was approved by the “Comité Ético de Investigación Clínica de Galicia”, and each of the institutional review boards of the participating hospitals. All samples were obtained with written informed consent reviewed by the ethical board of the corresponding hospital, in accordance with the tenets of the Declaration of Helsinki.

SNP genotyping and QC

Affymetrix array 6.0 (Affymetrix, CA, USA), which includes probes for almost 1M SNP markers, was chosen to obtain genome-wide coverage for Phase I genotyping. Genotype calling for Affymetrix 6.0 was performed with the Birdseed algorithm, included within Birdsuite v1.4 [11]. Samples were organized in 23 batches of 16<n<99 according to hospital of origin for computational purposes. We obtained valid genotypes for 909 622 SNPs by these means. Quality control of the data, performed mainly with PLINK v1.07 [12], included the removal of both SNPs and samples with genotyping success rates <99% (N=5984 for SNPs and N=0 for samples) and samples with discordant gender between clinical recorded data and Affymetrix-asigned sex (N=7). Hardy-Weinberg equilibrium (HWE) was evaluated and markers with P_HWE<1x10^-4 in controls were removed from further analyses (N=6984). SNPs with MAFs below 0.05 (N=221 799) were also eliminated due to low power to detect true signals and to avoid unnecessary noise. Finally, differential missingness between cases and controls was also accounted for by excluding markers with p-values below 1x10^-4 (N=137). This test compares genotyping error rates for the affected vs. unaffected groups in order to avoid an increase in false positive findings due to this bias. Finally, SNPs with poor clustering were also excluded after visualisation with Evoker (briefly, associated SNPs at a selected threshold were selected for comparison of the two intensity channels against each other to manually check the proper assignment of the genotype-calling algorithm) [13]. A total of 674 718 SNPs remained after this filtering.

To address the possibility of underlying population stratification, Principal Component Analysis (PCA) on a subset of 98 986 randomly chosen independent SNPs (pairwise r²<0.1) was also performed on the full Phase I cohort using the EIGENSOFT smartpca software [14]. Long-range LD regions, as described by Price et al.[15], were also removed from this analysis. Outliers, taken as samples spread on principal components 1 and 2 were removed from subsequent analyses, since they deviated from the main cloud. No evidence was found of population differences between cases and controls for the first 10 components of the PCA analysis, as stated in the Tracy-Widsom test (Figure 1a). Other potentially confounding variables, such as markers typed using the Nsp or Sty restriction endonucleases, hospital of collection, genotyping plate, or geographical origin of the samples were also checked for as sources for stratification (data not shown). All results were concordant with the original assumption of a single originating population except for hospital of origin. When considered as a confounding variable, the EPICOLON cohort clustered into three separate subgroups: samples from the Donostia hospital (VAS dataset), the only collection centre for the Basque Country region (North of Spain), samples from the Meixoeiro hospital (GAL dataset), the single collection point in Galicia (NW Spain), and all others (REST dataset) (Figure 1b). An additional PCA comprising the full Phase I cohort and the HapMap3 populations (all ancestries) was also performed to illustrate the clustering of these populations (Figure 1c) [16]. EPICOLON II samples that clustered away from the European end of the plot (showing evidence of non-European ancestry) were excluded from further analyses.

Figure 1

PCA analysis on the EPICOLON cohort. (a) filtered data by case/control; (b) filtered data by hospital of origin; (c) EPICOLON and HapMap3 populations. Significant differences may be seen in section b with Meixoeiro and Donosti hospitals deviating from the main cloud

The final case dataset comprised 1477 samples (848 cases and 629 controls). The total count per subgroup was 167 for VAS, 366 for GAL and 944 for REST.

Genotyping in Phase II was conducted by Sequenom MassARRAY technology (Sequenom Inc. San Diego, CA, USA). rs7087402 at 10q23.31 could not be included in the analyses for genotyping design reasons. Quality control was performed with PLINK similarly to Phase I. Genotyping for both Phases was performed at the Santiago de Compostela node of the Spanish Genotyping Centre.

Statistical analysis

Association analysis was assessed as a 1°-of-freedom χ² allelic test for each of the three subgroups independently for Phase I, and for second Phase replication, with PLINK [12]. The adequacy of the distribution of p-values was evaluated using quantile-quantile (Q-Q) plots of test statistics. Meta-analysis was also conducted using PLINK in Phase I. The method is based on a Mantel-Haenszel approach for data pooling. Cochran´s Q statistic and the I² heterogeneity index were also estimated to account for inter-population heterogeneity between groups, which was defined as I²>75% [17, 18]. For markers above this threshold, a random-effects model was considered, whereas fixed-effect results were otherwise reported. Risks (odds ratios, ORs) with 95% confidence intervals (CIs) associated with each marker were then estimated assuming the appropriate model. Phase II analyses were adjusted by age, given the mean differences between the case and control populations for these cohorts. Pooled analysis was performed by logistic regression, considering genotyping Phase and population subgroup as covariates. Associations by phenotype (age at diagnosis, MSI status, tumour location, presence of previous adenomas, family history of CRC and sex) were examined by logistic regression in case-only analyses for the two associated SNPs (Additional file 1: Table S2). Additional statistical calculations and plots were performed using R [19].

Imputation

Imputation between the two recombination hotspots encompassing each of the 24 loci that showed evidence of association in Phase I was accomplished with Impute v2 using two reference panels: 1000 Genomes Project (b36) for wide coverage, and HapMap3 (r2 b36) for deep coverage [16, 20, 21]. Results from the imputation were later tested for association with SNPTEST [20]. Imputation results were filtered by minor allele frequency (MAF) of the markers (SNPs with MAFs<5% were excluded, since the procedure generates genotypes for a high number of rare variants that could give spurious association results), by missing data proportion (set to a 5% max), and the frequentist-add-proper-info column of the output. This latter statistic is the ratio of the empirically observed variance of the allele dosage to the expected binomial variance p(1-p) at HWE, where p is the observed allele frequency from HapMap [22]. Optimal values should be within the (0.4-1) range and provide a measure for quality and accuracy of the imputation. Since the proportion of cases and controls deviates significantly from the standard 1:1, we also considered the possibility that the imputed genotype probabilities for each marker were different in both subsets. Since IMPUTE gives back the imputation results as a probability for each genotype, we decided to filter out SNPs for which the probability of two out of the three genotypes was ≥25% (i.e. the genotype in this sample was not clear) in at least 5% of the cases or the controls. This way we can eliminate samples for which genotype imputation has been inaccurate. Imputation results were plotted with the SNAP on-line tool [23].

Results

Stratification within the EPICOLON cohort

We observed using PCA that there was a batch effect due to differences by hospital at which sample had been collected, thereby dividing the EPICOLON cohort into three separate subgroups. We thus proceeded on the basis that each cluster - the GAL, VAS and REST case-control groups - was a separate sample set. Association results were then obtained for each of the subpopulations separately and then meta-analysed. Q-Q plots for the subgroups (after QC) showed an improvement in the systematic inflation for the distribution of the association p-values for the GAL, REST and VAS subgroups (Additional file 2: Figures S1A, 1B and 1C. respectively). Lambda genomic factor calculations (1.04192, 1.02323 and 1.02292 for GAL, VAS and REST, respectively) were consistent with no evidences of an increased false discovery rate. Additional file 2: Figure S1D represents the Q-Q plot distribution after meta-analysis of the three subpopulations. Some SNPs still seem to deviate from the expected distribution. These were later discovered to be artefacts of the calling procedure and were further removed with Evoker.

Association analyses

We found 64 SNPs at 24 genomic loci SNPs associated with CRC risk at P≤10^-4 (lowest p-value=9,7x10^-8 for rs11996339 at 8p12) (Table 1). Notably, I² heterogeneity values for the three subgroups (GAL, VAS, REST) were all 0 for these markers, thereby reflecting homogeneity in these associated SNPs.

Table 1 Association results for phase I

Using imputation, we examined the LD blocks defined by recombination hotspots (as obtained from Haploview in the CEU+TSI HapMap3-r2 populations) around the 64 SNPs for evidence of stronger signals of association. Additional file 1: Table S3 provides a summary of the loci and extent of the imputed regions. This analysis improved the association at 4 of the 24 loci: 1p33 (best SNP rs12060081); 14q31.3 (rs2057115); 15q21.3 (rs7176932); and 22q12.3 (rs17725348) (Additional file 3: Figure S2).

Replication

The ±1Mb regions flanking all 24 loci with evidence of association in Phase I were examined in the United Kingdom CORGI GWAS cohort through proxy SNP assessment [24]. Only 5 of these locations showed to have some CORGI associated SNP at an established p value threshold of P≤10^-4. However, LD analysis (in HapMap3-r2 CEU+TSI available data) showed them all to be independent signals (data not shown). Therefore, our attempt to replicate the association signals in silico was not successful for any of the variants.

Since there has been extensive literature on the differences amongst Northern and Southern European populations [25–28], we decided to perform a further PCA on 15,000 independent markers in order to compare allelic frequencies among the EPICOLON controls, the HapMap3 CEU and TSI populations, and the Wellcome Trust Case Control Consortium (WTCCC2) control cohorts [29]. This analysis effectively separated the Northern and Southern-European populations (Figure 2). Given this evidence, we decided to attempt replication of the best-associated markers (whether directly genotyped or imputed) at these 24 loci in an independent Spanish cohort (Phase II). Thirty-two SNPs were finally selected to be genotyped at this Phase II according to LD measures at the 24 loci and experimental design. The markers genotyped at each locus and their association values in the replication phase are described in Table 2.

Figure 2

PCA analysis on the WTCCC (Affymetrix 6.0 data), HapMap3 CEU and TSI and EPICOLON populations. A set of 15 000 independent markers was used to perform the analysis. The first eigenvector separates the Northern and Southern European populations.

Table 2 Association results for phase II and pooled analysis

Two out of the replicated 32 SNPs, rs12080929 (chromosome 1p33) and rs11987193 (chromosome 8p12) were successfully replicated at a nominal level of p<0.05 in this second Phase (P=0.044, OR=0.867 (0.722-0.994) and P=0.039, OR=0.847 (0.724-0.992), respectively). Although the association signals were very modest, pooled analysis of the data from both phases was consistent with the presence of a potential CRC susceptibility variant in these locations (pooled p values P=5.523x10^-3 and P=6.985x10^-5, respectively), and the signals remained significant across population subgroups (with the exception of the smallest VAS dataset) (Figure 3). Another two variants at locus 14q31.3 were significant in Phase II, but their OR were in different directions for each of the Phases, thereby reflecting a potential false positive event in Phase I. Given the marker sizes, this finding is entirely compatible with a random positive finding at a significance level of 0.05. With regards to phenotype analysis, rs12080929 seemed to be slightly overrepresented in males vs. females (P=0.042, OR=0.771 (0.600-0.991)), whereas rs11987193 was more prevalent in rectal cancers (P=0.028, OR=1.327 (1.031-1.707)). None of the other variables used in the subgroup analyses provided statistically significant results (Additional file 1: Table S2). It is also remarkable that none of the remaining three SNPs genotyped at these two loci (rs12080061 at 1p33, and rs11996339 and rs12548021 at 8p12) appeared to be replicated in this second phase, although rs12080061 showed a borderline p value in Phase II (P=0.087); in the case of rs11996339 and rs12548021 at 8p12 r-squared LD values among the three SNPs seem to show that the three markers are independent.

Figure 3

Forest Plots for rs12080929 and rs11987193. The figure represents the odds ratios and 95% confidence intervals for the two markers in all 3 populations subgroups from Phase I, the Phase II replication dataset and the pooled analysis of both Phases

Previous susceptibility loci

In addition to the search of new susceptibility variants, we also investigated the association signals for 19 of the known CRC susceptibility variants [3–5]. rs5934683 on Xp22.2 could not be evaluated due to the fact that sexual chromosome data need to be processed differently. Considering these markers were described from Illumina array tagSNP panels, most of them were not directly genotyped in our chip; therefore we proceeded with the evaluation of the association signals by considering the closest related proxy SNP (Table 3). Direct evidence of replication (taken as the presence of an associated SNP with one-tailed P<0.05 in the same LD block as the described tagSNP) was found for 2 of the SNPs (Table 3). The remaining loci, although not significant, showed ORs in the same directions as those described in the literature. We must however highlight that some of these SNPs (rs4444235 and rs1957636 at 14q22 and rs961253 andrs4813802 at 20p12) have already been found to present differences in Northern and Southern European populations, a fact that is consistent with our study not being able to replicate the association signals at these loci [28]. Imputation of the LD regions around these associated loci was conducted to search for an enhancing of the signals. No significant improvements were found, except for locus 15q13, for which an imputed SNP, rs16970016, 15kb upstream the GREM1 gene, scored the best p value in our dataset (P=9.847x10^-5). This SNP has a good pairwise relationship with the formerly described rs4779584 (r²=0.882) (data not shown).

Table 3 Replication results for the already-described loci

Discussion

Genome-wide association studies have so far successfully identified 20 CRC susceptibility SNPs [3–6]. Although this has been a significant improvement in the unravelling of the genetic basis of the disease, these variants alone do not completely explain all the inherited variation that has been attributed to CRC.

Following the lead of the previous studies, we addressed the issue of trying to detect new colorectal cancer susceptibility variants through the performance of a GWAS in a Spanish cohort. This was the first attempt to perform a CRC GWAS in a Southern European population. By these means, we were able to positively identify two new candidate variants that have shown good evidence of association with CRC risk: rs12080929 at 1p33 and rs11987193, at 8p12.

Previous susceptibility loci

During the analysis, we were faced with the fact that, although there were no differences between case and control populations, there was a significant stratification issue determined by the hospital of origin of the samples. Because of this, the analyses had to be modified to match our case scenario without losing significant power. Nevertheless, the substructure in our cohort did not seem to greatly affect outcome quality. The evaluation on the 19 out of 20 already-described signals achieved direct replication for 2 of the loci (11q23 and 18q21). The other 17 markers did not show evidence of association, probably due to the lack of power in our cohort to detect such moderate effects. Nevertheless, OR directions were consistent with those previously published. We must highlight at this point that the best-associated markers for these regions did not always match with the best proxy for the already described SNPs. This would make sense if we consider that any given GWAS relies on an indirect approach, and we would expect the associated SNPs to only be tagging the real causative variant. Results for allele frequencies and ORs seem consistent with the bibliography [3–5]. We consider the replication of these loci an important achievement, since the majority of these association signals had not been previously evaluated in Southern European cohorts (with the notable exception of rs16892766 at 8q23.3, rs10795668 at 10p14, rs3802842 at 11q23, rs4779584 at 15q13, rs4444235 and rs1957636 at 14q22, and rs961253 and rs4813802 at 20p12 [5, 28, 30, 31].

Spanish GWAS results

The association analysis in itself provided with positive results for 24 different genomic loci at a p-value <0.0001. A first attempt at replication was aimed by inspection of these association signals on the British CORGI cohort [24]. However, none of the signals seemed to be shared between datasets. This lack of replication could be due to both false positive findings and artefacts from the calling algorithm, or to real differences between both populations leading to dissimilar abilities to tag the real causative variant [32]. The latter option has been recently proven to be possible, since differences in MAFs between Northern and Southern European populations have been described for SNPs rs4444235, rs1957636, rs961253 and rs4813802 [28].

A PCA analysis on the EPICOLON samples compared to the WTCCC control cohorts and the data from the HapMap3 CEU and TSI populations showed clear differentiation between the Northern and Southern European populations. Although not significant, SNP loadings also evidenced principal component 3 to be exclusively driven by a region of chromosome 8 (7.2-12Mb) where a common inversion is known to occur [33, 34], whereas Eigenvectors 4-7 were driven by the HLA-A locus in the 6q21.2-21.3 region of chromosome 6, which has been also described as highly variable between populations [35]. Given this evidence, we proceeded on to replicate these 24 loci in an independent Spanish cohort.

SNPs rs12080929 and rs11987193 were successfully replicated in Phase II analyses. The former SNP, rs12080929 is located on an intronic position within hypothetical locus LOC388630 on 1p33. This predicted gene is believed to code for a single-pass type I membrane protein. Moreover, it lies 252kb upstream the SLC5A9 gene, a member of the solute carrier family, which could be a feasible regulation target. SLC proteins constitute good candidates to harbour CRC susceptibility loci, since some family members have been proven to act as tumour suppressors. In fact, the SLC5A8 gene is properly expressed in normal colon, but silenced in aberrant crypt foci through gene methylation [36].

The rs11987193 SNP is located in the 8p12 locus, 128kb downstream DUSP4. This gene is a member of the dual kinase phosphatase family, which are well-known tumour suppressors too [37]. They act through the downregulation of MAP kinases, thus preventing cellular proliferation and differentiation. Deletions in this gene have already been described to happen in other types of cancers, such as those of the breast and lung. In the case of CRC, DUSP4 expression appears to be modulated by KRAS mutations [38–40]. Moreover, it has recently been described that DUSP4 expression is associated with microsatellite instability in CRC and causes increased cell proliferation [41].

The fact that these two SNPs were not replicated during the initial assessment of the association signals in the CORGI cohort, together with the North-South discrepancies seen in the PCA analysis, could be a sign of differences in the tagging of the real causative variant amongst populations. Even when Europeans are presumed to be genetically homogeneous, it is not unrealistic to believe that punctual LD variations may be actually happening within populations, and that these may constitute a certain impediment in our ability to replicate association signals [25, 26].

Although none of the markers reached a final genome-wide significant p value, two of the SNPs in our study (rs12080929 and rs11987193) were favourably replicated in the second Phase and pooled analysis. However the limitations of this study (namely the modest sample sizes, further power restrictions also derived from the case-control age bias in Phase II and the potential differences between northern and southern European populations), we believe these two regions are good candidates for CRC susceptibility loci. The peculiarities of these loci, particularly those relating to potential Northern-Southern European differences, may have important repercussions on subsequent analysis. For this reason, the eventual identification of the functional variant is of uttermost importance. Finer mapping of the locus, coupled with additional replication efforts in larger cohorts will be needed to fully ascertain the relationship between these variants and disease. Moreover, it is important to highlight that, were these two SNPs true susceptibility variants, they would constitute a decrease in the CRC missing heritability fraction. This is an essential point in our road towards the identification of high-risk individuals within populations [7].

Conclusions

The novelty of this study was the use of a Southern-European dataset to perform a CRC GWAS study that led to the identification of two new candidate variants: rs12080929 at 1p33 and rs11987193 at 8p12. These two SNPs are located near the SLC5A9 and DUSP4 loci, respectively. Some family members of the SLC proteins, as well as the DUSP gene, have been proven to act as tumour suppressors. Therefore, both of them could be good functional candidates for the association signals. Finer mapping and further replication in larger cohorts will be needed to ascertain their relationship with CRC susceptibility. The peculiarities of these loci, particularly those relating to potential Northern-Southern European differences, may have important repercussions on subsequent analyses.