Association of Single Nucleotide Polymorphisms in Glycosylation Genes with Risk of Epithelial Ovarian Cancer

Epithelial ovarian cancer has a recognized genetic component. There is an autosomal dominant pattern of ovarian or breast-ovarian cancer in some families (1, 2) with a high lifetime risk (3-5). Mutation screening of BRCA1 and BRCA2 in population-based series of ovarian cancer cases has shown that mutations in these genes collectively account for only a small percentage of cases (3, 6). Ovarian cancer also occurs in Lynch syndrome II, an autosomal dominant predisposition to cancers of the colon, endometrium, ovary, and other sites (7, 8). The genetic defects underlying Lynch syndrome II include the mismatch repair genes MLH1, MSH2, PMS1, PMS2, and MSH6 (9, 10). The exact contribution of these genes to ovarian cancer risk is not known but is estimated to be lower than that of BRCA1 and BRCA2. Collectively, these high-risk, high-penetrance genes are rare in the population and are estimated to account for 10% to 15% of ovarian cancer (11). Studies to identify other chromosomal regions that might harbor major genes for ovarian cancer risk provide little evidence for the existence of additional high-risk genes for ovarian cancer susceptibility (12). Thus, there is emerging consensus that most of the genetic component of ovarian cancer risk is due to genetic polymorphisms that confer low to moderate risk.

A common approach to identify risk variants is to rely on known biology to identify plausible candidate genes. The present report was motivated by two converging lines of evidence from molecular biology and epidemiology. The mucin glycoprotein MUC1 is frequently overexpressed in epithelial ovarian cancer (13) and found in circulation (14). Mucin glycoproteins have a protein backbone composed of repetitive domains rich in peptides that can serve as sites for O-linked glycoslyation. Multiple O-linked oligosaccharides bind to these sites and represent 50% to 80% of the total glycoprotein mass. Normal epithelial cells express a highly glycosylated form of MUC1 at low levels, but ovarian cancer cells express high levels of a hypoglycosylated form (15). Two recent studies provide evidence that links circulating levels of anti-MUC1 antibodies with risk factors for ovarian cancer. The first report (16) observed that oral contraceptive use, breast mastitis, bone fracture or osteoporosis, pelvic surgeries, nonuse of talc in genital hygiene, and (to a lesser extent) current smoking predicted anti-MUC1 antibody levels among 705 women selected as controls for an ovarian cancer study. A second report by the same group (17) on 721 controls found that early age at first birth, menstrual cycles longer than 30 days, and lifetime number of ovulatory cycles were inversely associated with anti-MUC1 antibodies.

We hypothesized that polymorphisms in genes encoding glycosylation enzymes may be risk factors for epithelial ovarian cancer. The current study compared the frequencies of genotypes and selected haplotypes between women recently diagnosed with epithelial ovarian cancer and a group of control women without ovarian cancer. The findings may have important implications for our understanding of the pathogenesis of this important malignancy.

Details of the study design have been presented previously (18). Briefly, the subjects recruited for this research project were identified through two institutions: Duke University and Mayo Clinic. The study is ongoing, with the current results based on recruitment from June 1999 to March 2006. The protocol was approved by the institutional review board at each institution, and all study subjects provided written informed consent. Although there were several differences in the specific study designs used at each site, the overall inclusion-exclusion criteria were similar. Cases were Caucasian or African American women ascertained within 1 year of a diagnosis of histologically confirmed primary epithelial ovarian cancer, either borderline or invasive. The lower age limit was 20 years at both sites, but Mayo had no upper age limit, whereas Duke excluded women ages >74 years.

The Mayo ascertainment began in January 2000 and was clinic based. The catchment area was limited to a six-state region that represents >85% of all ovarian cancer cases seen there: Minnesota, Iowa, Wisconsin, Illinois, North Dakota, and South Dakota. We selected clinic-based controls from women seeking general medical evaluation, frequency matched to cases on age (5-year age category), race, and state of residence. Potential controls were excluded if they had a history of ovarian cancer or oophorectomy. Response rates for those invited to participate at the Mayo site were 83% for cases and 74% for controls.

The Duke study is population based with a rapid case ascertainment network covering a 48-county region of North Carolina. Recruitment has been ongoing since May 1, 1999. List-assisted random digit dialing and Health Care Financing Administration roster methods were used to identify control subjects. Controls were frequency matched to the cases on the basis of race (Black versus non-Black), age (5-year age categories), and residence. The response rate was 71% among eligible cases (87% among cases who could be contacted) and 64% among the controls.

Information on known and suspected ovarian cancer risk factors and demographic data were collected through in-person interviews. Mayo cases were interviewed by telephone if discharged from the hospital before contact about the study. Similar questionnaires were used at each institution. Information collected included race/ethnicity, menstrual and reproductive history, use of exogenous hormones, medical and surgical history, height and weight 1 year before the interview, use of tobacco, education level, and family history of breast or ovarian cancer in first- or second-degree relatives.

Genomic DNA was obtained from cases and controls in one of two ways. At Duke, venipuncture was done at the conclusion of the interview. At Mayo, participants had an extra vial of blood drawn in the course of their scheduled medical care. DNA was extracted from fresh peripheral blood using the Gentra AutoPure LS Purgene salting out methodology (Gentra). Due to limited quantity of available DNA from Duke subjects, we did whole-genome amplification (WGA) on all Duke samples (534 cases and 568 controls) using the REPLI-G protocol (Qiagen), with 200 ng genomic DNA. DNA concentrations were adjusted to 50 ng/μL and verified using PicoGreen dsDNA Quantitation kit (Molecular Probes). The samples were bar-coded to ensure accurate and reliable sample processing and storage.

Genes involved in the glycosylation process were identified through several sources, including peer-reviewed published literature (14, 19) and the Cancer Genome Anatomy Project Biokarta and Kegg pathway databases. A total of 26 genes were selected: FUT10, FUT3, FUT5, FUT6, FUT7, FUT8, FUT9, GALNT1, GALNT14, GALNT2, GALNT3, GALNT5, GALNT6, GALNT7, MGAT2, MGAT3, MGAT5, POFUT1, SIAT1, SIAT4A, ST3GAL3, ST3GAL4, ST3GAL5, ST6GALNAC5, ST8SIA2, and ST8SIA5.

For each gene, chromosome and protein attributes were selected and the data mined from the Ensembl database version 34 (Biomart) using the gene reference sequence identification number (RefSeq ID) and the approved gene symbol from HUGO or Entrez Gene. The chromosomal location on build 35 and strand (forward or reverse) were provided to Illumina. Illumina verified chromosomal coordinates. We requested all single nucleotide polymorphisms (SNP) within each gene as well as up to 10 kb in the 5′ and 3′ flanking regions and all nonsynonymous SNPs with a minor allele frequency (MAF) > 0.05 and Illumina Design Score > 0.6. The Illumina Assay Design Tool database includes all SNP data contained in the public domain, filtering out SNPs that are not suitable for the Illumina platform, such as insertions/deletions, tri- and tetra-allelic SNPs, and SNPs that are not uniquely localized.

Ninety-six glycosylation SNPs were included in an Illumina GoldenGate assay for 1,536 SNPs. Ninety-three percent of the SNPs had SNP_scores greater than 0.6 and none had scores less than 0.4. 2368 samples (1,086 genomic and 1,282 WGA; 250 ng each) were genotyped on 1,967 unique subjects following the Illumina protocol, including 930 cases (396 genomic and 534 WGA) and 1,037 controls (469 genomic and 568 WGA). Genotype calls were made using the Genotyping module of BeadStudio 2 software. Genomic and WGA DNA were analyzed separately, as WGA DNA clustered differently from genomic DNA. For this project, samples with GenCall scores below 0.25 and/or call rates below 90%, and SNPs with GenCall scores below 0.4 or call rates below 90%, were failed. For genomic DNA, 1,492 SNPs passed this initial quality assurance cutoff, whereas 1,435 passed in the WGA DNA set.

Several quality-control procedures were established. For the genomic DNA, we included eight replicates of a CEPH family trio from the Coriell Institute (mother, father, and child) and replicates of an additional three standard DNAs in each 96-well plate. The replicate and inheritance data were used to review and refine clustering. In addition, two samples per 96-well plate were blindly duplicated (n = 20). Among the WGA samples, there were 88 replicates (same WGA preparation) and 15 replicate WGA (separate WGA preparation) samples genotyped. In addition, 124 of the WGA samples had sufficient genomic DNA to allow genotyping and thus a means to monitor the concordance between genomic and WGA DNA genotype calls. Of these, two had an unusually high number of discrepancies between the genome and WGA'd samples due to sampling errors and were therefore excluded. These genomic DNAs were analyzed along with the 1,086 genomic DNAs to generate genotype calls.

We attempted to genotype 2,058 subjects on 96 SNPs in 26 glycosylation genes. Of these subjects, we excluded 87 with poor clustering and therefore failed for every SNP, 2 with low call rates (<95% of SNPs successfully genotyped) and 2 with an unusually high number of discrepancies between genomic DNA genotype results and WGA results. Of the 1,536 SNPs targeted for genotyping, we excluded 44, which failed completely (or nearly completely), 6 because of low call rates (<95%), and 7 with no genetic variability (that is, monomorphic), resulting in a total of 1,480 SNPs. Of the 57 SNPs that failed, 3 were in glycosylation genes. The exact number of SNPs for analysis varied slightly per study site: 1,421 had results of sufficient quality for both Mayo and Duke subjects, 51 were successfully genotyped for Mayo subjects but failed for Duke subjects, and 7 were successfully genotyped for Duke subjects but had either low call rates or no genetic variability for the Mayo samples. Genotypes at these 58 problematic loci were coded as missing. This resulted in a final sample size of 1,967 subjects and 93 SNPs in glycosylation genes.

Before analysis, we determined descriptive statistics using frequencies and percents for categorical variables and means and SDs for continuous variables. The distributions of covariates were compared across study site and case status using ANOVA methods for continuous variables and χ² tests for categorical variables. SNP genotype frequencies among the controls were tested for Hardy-Weinberg equilibrium (HWE) using χ² goodness-of-fit tests.

We first did SNP-specific analyses to examine the main effects of the SNPs on risk of ovarian cancer. We fitted log additive, dominant, and recessive models for each individual SNP. Unconditional logistic regression models were fitted using the R function glm (http://www.r-project.org) to estimate odds ratios (OR) and corresponding 95% confidence intervals (95% CI) between genotypes and case status. Wald χ² tests were calculated using the R function ANOVA to obtain P values for SNP effects. Due to differences in genotype frequencies by race, separate logistic regression models were fitted for Caucasian and non-Caucasian subjects. All models were adjusted for the design variables of geographic area/study site and age group. We also adjusted for several nongenetic risk factors: body mass index (BMI), months of hormone replacement therapy (HRT) use, months of oral contraceptive use, and parity/age first birth combination. Although none were traditional confounding factors, results were similar regardless whether they were in the model. In addition, we carried out sensitivity analyses by excluding the borderline cases in the Caucasian samples. Individual SNP associations with risk were compared before and after exclusion of borderline cases.

We next examined the association of haplotypes in each glycosylation gene with ovarian cancer status. These analyses were restricted to Caucasian subjects due to the relatively low number of enrolled non-Caucasian subjects and done using the statistical program Haplo.stats (20). For each gene with multiple SNPs, we first used the function Haplo.score to estimate haplotypes and test global significance, adjusting for the covariates mentioned above. If the global P value for haplotype association was smaller than 0.05, we examined individual haplotypes by comparing the risk of ovarian cancer associated with each inferred haplotype to the risk associated with the highest estimated frequency haplotype. The estimated OR and the corresponding P value were obtained using the function Haplo.glm with adjustment for the covariates. Rare haplotypes were pooled into a single category to reduce the burden of sparse table cells.

In addition to performing single SNP and haplotype analyses, we assessed potential modifying effects of selected demographic and clinical variables by fitting a series of interaction models. The modifying effects of the following variables were considered based on published association with MUC1 antibody levels: ever smoked, BMI (median split), live birth (0 versus ≥1), ever used oral contraceptive, ever used HRT, and menopausal status (premenopausal and postmenopausal). These analyses were also restricted to Caucasian subjects. To avoid further the possibility of sparse table cells, all the SNPs were first modeled as dichotomous variables based on the presence (one or two copies) or absence (zero copies) of the variant alleles. Two SNPs (rs3828139 and rs37460) with significant recessive main effect and high MAF were also modeled as dichotomous variables based on the presence or absence of the common alleles. Similarly, each environmental variable was dichotomized based on either a median split of the distribution among controls or a pooling of categories. For each interaction between the nongenetic variables and SNPs, a logistic regression model was fitted with the corresponding SNP, the corresponding environmental variable, the interaction term, and other potentially confounding covariates. Based on the models, P values for testing the interaction effects were obtained using Wald χ² tests. For each interaction with P < 0.05, we obtained OR and 95% CI estimates for individual strata of SNP and environmental variable level.

The distribution of nongenetic risk factors by study site and case-control status is summarized in Table 1. Cases tended to be more obese, less likely to use oral contraceptive pills, have lower parity, and more likely to report a family history of ovarian cancer. Case-control differences were generally similar across sites, with the exception of HRT use, with Duke cases having higher exposure levels than the other groups.

Differences in genotype frequencies between cases and controls were compared by race and study site for each of the 93 SNPs in the 26 genes investigated. Caucasian and non-Caucasian subjects had significantly different allele frequencies for some of the SNPs, so all analyses were stratified by race. There were no differences in allele frequency by study site, so data from Mayo and Duke were combined for analysis. After fitting dominant, recessive, and additive models and adjusting for age and the covariates listed as a footnote in Table 2, statistically significant main effects before adjusting for multiple testing were identified in 8 of the 26 genes studied and 9 of the 93 SNPs (Table 2). One gene (ST3GAL3) had 2 SNPs associated with risk. It is interesting to note that 7 of the 9 SNPs were inversely associated with risk and that only 1 of the 9 was based on a dominant genetic model. Only one SNP (rs17647532 in GALNT1) was not in HWE (P = 0.002), but this was limited to controls from Mayo (P = 0.0003), not Duke (P = 0.35), although the estimate of the MAF was 0.10 at both sites. Table 2 also presents the results for Mayo and Duke samples separately.

Given the large number of statistical tests (93 SNPs ×3 genetic models), one needs to be cautious about interpretation. Therefore, we used the Bonferroni method to adjust for multiple tests and keep the overall type I error rate under 5%. Only the GALNT1 SNP (rs17647532) remained statistically significant. This SNP was rare (MAF of 0.10); therefore, the results are based on few homozygous carriers (1 case and 20 controls). Given the lack of HWE for this SNP among the Mayo controls, we did post hoc analyses by study site to see if this influenced the findings. The only case homozygous for the minor allele was ascertained at Mayo, which meant the association could not be estimated among the Duke subjects. However, the combined OR (95% CI) of 0.07 was somewhat attenuated among the Mayo subset [0.14 (0.02-1.14)]. This suggests that the overall result was not affected by the absence of HWE at Mayo. The SNP is located in the promoter/regulatory region in GALNT1. We applied bioinformatic tools to help interpret the plausibility of this finding. In particular, we ran the program FASTSNP (21) to predict the potential functional significance of the polymorphism. The SNP was categorized to have low to medium risk, indicating that it may affect the level, location, or timing of gene expression. Because there is some evidence that borderline ovarian cancers may be etiologically different from invasive epithelial cancers, we did post hoc sensitivity analyses after excluding cases with borderline tumors. Results were unchanged with one exception: the two SNPs in ST3GAL3 both became more strongly and significantly inversely associated with risk (rs3828139 OR, 0.72; P without multiplicity adjustment = 0.010; rs37460 OR, 0.71; P without multiplicity adjustment = 0.0069).

The analysis of Caucasian subjects continued with haplotype analyses for genes that had more than one SNP genotyped. Based on the global score statistic, only GALNT1 was significantly associated with risk (P = 0.038). There were five measured SNPs in GALNT1 (rs6507133, rs11663626, rs556736, rs607498, and rs17647532) and they were in LD with r² values ranging from 0.06 to 1.0. Although 32 haplotypes were possible, we observed only 5 common haplotypes (frequency > 0.05); their frequencies and individual haplotype effects are presented in Table 3. Because the individual SNP analysis revealed a SNP (rs17647573) with significant recessive effect, we fitted the haplotype genetic model to be a recessive model. The result of the global score test was largely the reflection of the ACTAC haplotype (0.10 frequency) that was associated with decreased risk (OR, 0.07; P = 0.01) compared with the most common haplotype TTTAT (0.39 frequency), consistent with an untyped causal allele or the haplotype.

SNP-specific results based on the African American subjects are shown in Table 4. A total of 10 SNPs in 7 genes had multivariate-adjusted P values smaller than 0.05 from any of the three genetic models. All of the SNPs were in HWE. Three of the 7 genes were also associated with risk among the Caucasians (GALNT1, MGAT5, and ST3GAL3), but only 1 of the SNPs was in common: rs1257189 in MGAT5. Most of the SNPs (7 of 10) were associated with increased risk, and most often (6 of 10) this was under an additive genetic model. None of these associations were statistically significant after adjustment for multiple comparisons.

Exploratory models were then fitted to the Caucasian subjects to test gene × environment interactions (technically SNP × environment interactions). We included a fixed set of nongenetic variables based on the epidemiologic literature linking them with MUC1 antibodies (16, 17): BMI, oral contraceptive use, HRT use, menopause status, parity, and smoking. Rather than fitting interactions under all genetic models, this exercise was limited to a dominant genetic effect for all 93 SNPs plus 2 recessive models (based on main effect results). There were 18 interactions significant at the 0.05 level (very close to the expected value of 20 based on chance alone). Because we conducted 570 tests (93 SNPs × 6 environmental factors + 2 SNPs under other genetic models × 6 environmental factors), the P value cutoff based on the Bonferroni method is 0.05/570 or 0.000088. None of the interactions exceeded this threshold (data not shown).

This study sought to evaluate whether inherited variation in genes involved in the glycosylation of MUC1 were potential risk factors for ovarian cancer. The analysis focused on 93 SNPs in 26 genes and considered different genetic models, haplotypes, and gene × environment interaction. Before adjustment for multiple comparisons, statistically significant main effects were identified among Caucasians in 6 of the 26 genes and 9 of the 93 SNPs and another 10 genes had at least one statistically significant interaction effect with environment. After adjustment for multiple comparisons, one SNP in GALNT1 remained statistically significant and inversely associated with risk. Seven genes were significantly associated with risk among African Americans: three genes (GALNT1, MGAT5, and ST3GAL3) and one SNP (rs1257189 in MGAT5) were consistent with the results among the Caucasians, but the results did not hold after correction for multiple testing. To our knowledge, this is the first study of glycosylation genes as risk factors for ovarian cancer and the early results merit further investigation.

Interpretation of the current results is limited by incomplete information regarding the biological significance of the SNPs studied. During the selection of the genes and SNPs, we targeted exons and regulatory regions but in some instances necessarily included intronic SNPs to ensure good coverage of the candidate gene. We also selected every cSNP reported in the public databases, but unfortunately none of the 23 cSNPs passed the GoldenGate assay review at Illumina and were dropped. Thus, we were left with many SNPs that had no known biological function. Another issue is the absence of HWE among the Mayo controls for the candidate gene found to show the strongest association with ovarian cancer risk (GALNT1).

The results were consistent when analyses were restricted to the Duke subjects, careful inspection of the Illumina results provided no indication of genotype errors, and the number of SNPs that violated HWE (n = 7) is consistent with expectations based on chance. Moreover, one must also consider the biological importance of O-linked glycosylation of proteins begins with the addition of a single GalNAc monosaccharide to a serine or threonine residue on the polypeptide. Attachment is catalyzed by a UDP-N-α-d-galactosamine: polypeptide N-acetylgalactosaminyltransferase (ppGalNAc-T). During glycosylation, additional glycosyltransferases are responsible for catalyzing other types of glycosidic linkage—an extensive family of up to 24 ppGalNAc-T’s has been described (19).

Despite the incomplete data on functional significance of the actual SNPs studied, the possibility that genes involved in glycosylation may be risk factors for ovarian cancer could nonetheless contribute to our understanding of the biology of this disease. Although there is a substantial literature that expression of mucins in epithelial ovarian cancer have diagnostic and prognostic value (22, 23) and may represent therapeutic targets (24), the basis for their altered expression remains unknown. There is now emerging evidence that glycans in general (25), and MUC1 in particular (22) may play a role in host defense against pathogenic molecules. Because GALNT1 is known to participate in the glycoslyation of MUC1 (26), and because MUC1 is involved in defense against pathogenic agents, it is tempting to speculate that genetic polymorphisms that encode enzymes involved in glycosylation may reflect interindividual differences to respond to an infectious agent. This may partly explain the observation that tubal ligation is inversely associated with ovarian cancer risk (27-30). Although this association has been hypothesized to reflect tubal ligation effects on hormonal mechanisms (31) or inflammation (32), the current findings lend some credence to the hypothesis raised by Walhberg (33) of an infectious etiology to ovarian cancer.

A critical issue for studies of this type for which there is no consistent agreement among the scientific community is adjustment for multiple testing. Exacerbating the debate for the current study was the a priori decision to fit three genetic models to the data. This turned out to be potentially important, as among the Caucasians only 1 of the 9 statistically significant SNPs (before adjustment for multiple comparisons) was based on a dominant genetic model and 6 of 10 statistically significant SNPs among the African Americans were based on an additive model. Regardless, there were several strategies and methods one could have adopted that merit some discussion. For example, adjustment could be based at a gene level or at a SNP level, with the former being a more relaxed stringency in which fewer degrees of freedom are spent. In this report, we opted for adjustment based on the number of SNPs examined (n = 93), with further correction for the fact that we examined the association of these SNPs with ovarian cancer risk as dominant, additive, or recessive effects. This led to an even greater stringency on the significance level. Although this is the most conservative approach to protect against false-positive results, there may be some true-positive findings that are being downplayed but still merit consideration in future studies.

Strengths of the study include the high participation rates, relatively large sample size, the careful inclusion/exclusion criteria, the tight quality control on the genotype data, and data on nongenetic risk factors that were comparable across sites and available for adjustment in statistical models. Despite these strengths, interpretation of the findings should consider that there were too few African American subjects for stable estimates of effect. Genotyping included a mixture of genomic and amplified DNA. Although stringent call rates were applied, some SNPs and subjects did drop out from analysis.

This decreased sample size slightly, but the dropout rates did not vary by case-control status and the concordance of native to WGA DNA was 99.16%, so there is little reason to believe that this introduced any sort of bias. The Mayo ascertainment was not population based, but the proportion of in-state cases was similar for cases and controls (52% and 58%, respectively). Moreover, residence was adjusted for in the multivariate analyses, so any imbalances were accounted for statistically. Finally, although the results may provide support for an infectious component of ovarian cancer, we did not have data on tubal ligation, which could have been instructive. Identification of the genes involved in glycosylation relied on a review of the existing literature and consultation with experts on MUC1 biology. Although we were thorough in our approach, there may be additional genes involved in glycosylation that are relevant to MUC1.

In summary, this large case-control study of epithelial ovarian cancer provides some evidence to support the role of interindividual differences in glycosylation genes as risk factors for the disease. More detailed investigation is warranted with a larger and more ethnically diverse study sample, greater coverage in genetic variation, stratification by histologic subtype, and inclusion of additional nongenetic covariates that might help further elucidate the actual pathogenesis.