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

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.

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].

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

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.

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.

Previous susceptibility loci

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].