Identification of a rare coding variant in complement 3 associated with age-related macular degeneration.

Macular degeneration is a common cause of blindness in the elderly. To identify rare coding variants associated with a large increase in risk of age-related macular degeneration (AMD), we sequenced 2,335 cases and 789 controls in 10 candidate loci (57 genes). To increase power, we augmented our control set with ancestry-matched exome-sequenced controls. An analysis of coding variation in 2,268 AMD cases and 2,268 ancestry-matched controls identified 2 large-effect rare variants: previously described p.Arg1210Cys encoded in the CFH gene (case frequency (fcase) = 0.51%; control frequency (fcontrol) = 0.02%; odds ratio (OR) = 23.11) and newly identified p.Lys155Gln encoded in the C3 gene (fcase = 1.06%; fcontrol = 0.39%; OR = 2.68). The variants suggest decreased inhibition of C3 by complement factor H, resulting in increased activation of the alternative complement pathway, as a key component of disease biology.

Genetic and environmental factors contribute to age-related macular degeneration (AMD)[1,2], a major cause of vision loss in elderly individuals[3]. Pioneering discovery of association of AMD with complement factor H (CFH[4-6]) was quickly followed by the identification of additional susceptibility loci that now include ARMS2/HTRA1[7,8] and complement genes C3, C2/CFB and CFI[9-12]. Genome-wide association studies (GWAS) of AMD cases and controls have now revealed common susceptibility variants at ~20 different loci[13,14] and begun to uncover specific cellular pathways involved in AMD biology.

While common variants tag the associated genomic region, rare coding variants can provide more specific clues about the underlying disease mechanism[15]. For example, rare variant R1210C in the CFH gene was recently associated with a large increase in AMD risk using targeted sequencing of rare CFH risk haplotypes[16]. The resulting altered protein has decreased binding to C3b, C3d, heparin and endothelial cells[17-19]. A reduction in CFH’s ability to inactivate C3, leading to increased cell killing activity by the complement pathway, could contribute to AMD – a much more specific and testable hypothesis about disease mechanism than provided by common CFH variants whose mechanistic consequences are unclear.

To systematically identify rare, large-effect variants, we carried out targeted sequencing of eight AMD risk loci identified in GWAS[20] (near CFH, ARMS2, C3, C2/CFB, CFI, CETP, LIPC and TIMP3/SYN3) and two candidate regions (LPL and ABCA1) (Supplementary Table 1). We re-sequenced these regions in 3,124 individuals (2,335 cases and 789 controls) recruited in ophthalmology clinics at the University of Michigan and at the University of Pennsylvania and among Age-Related Eye Disease Study (AREDS) participants[20,21]_ENREF_17. Genomic targets were enriched using a set of 150-bp probes designed by Agilent Technologies, and sequence data was generated on Illumina Genome Analyzer and HiSeq instruments. The ten loci comprised 115,596 nucleotides of protein coding sequence and totaled 2,757,914 nucleotides overall. We designed probes to capture 111,592 protein coding nucleotides (96.5% of coding sequence) and 966,607 nucleotides overall (35.1 % of the locus sequence), generating an average of 123,221,974 mapped bases of on-target sequence per individual (127.5× average depth counting bases with quality >20 in reads with mapping quality >30, after duplicate read removal); 98.49% of sites with designed probes were covered at >10× depth. We applied variant calling tools and quality control filters similar to those used to analyze NHLBI Exome Sequencing Project data[22] (Supplementary Table 2). We identified an average of 1,714 non-reference sites in each sequenced individual. In total, this resulted in 31,527 single nucleotide variants of which 18,956 were not in dbSNP 135. Discovered sites included 834 synonymous variants, 1,379 nonsynonymous variants and 43 nonsense variants, most of which were extremely rare (see Supplementary Table 3). Among 13 samples sequenced in duplicate, genotype concordance was 99.82% (when depth >10×). Among 908 samples previously examined with GWAS arrays[20], sequence-based genotypes were 98.99% concordant with array-based calls (again, when depth >10×).

In an initial comparison of AMD cases and controls (see Supplementary Table 4), no rare coding variants with frequency <1% reached experiment wide significance (p < 0.05 / 31,527 = 1.6×10−6, including all discovered variants, or p < 0.05 / 1,422 = 3.5×10−5 considering only protein altering variants), although several showed encouraging patterns. For example, rare variant R1210C in the CFH gene was observed in 23 of the 2,335 sequenced cases, but in none of the 789 sequenced controls (exact test p=0.0025). Common variants in several loci exhibited strong evidence of association, including in CFH (peak variant rs9427642 with case frequency fcase = 12%, control frequency fcontrol = 27%, P-value = 2.52×10−48), ARMS2 (rs10490924, fcase = 33%, fcontrol = 18%, P-value = 5.48×10−27), C3 (rs2230199, fcase = 25%, fcontrol = 17%, P-value = 3.94×10−9) and C2/CFB (rs556679, fcase = 7%, fcontrol = 12%, P-value = 1.32×10−10).

A key requirement for establishing significance of rare disease associated variants is the availability of sufficient numbers of control samples. To increase power, we sought to identify additional controls and focused on samples from the NHLBI Exome Sequencing Project (ESP)[23], which sequenced 15,336 genes across 6,515 individuals. Sequence data for our samples and the NHLBI Exome Sequencing Project samples were analyzed with the same analysis pipeline, which minimized potential differences due to heterogeneity in analysis tools and parameters. To further avoid sequencing and variant calling artifacts, we restricted our analysis to sites within regions targeted in both sequencing experiments, genotyped and covered with >10 reads in >90% of the samples examined in each project, and >5-bp away from insertion/deletion polymorphisms catalogued by the 1000 Genomes Project[24]. Since careful matching of genetic ancestry is critical for rare variant association studies[24,25], we selected an ancestry-matched subset of our samples and of samples from the NHLBI Exome Sequencing Project. We used principal component analysis to construct a genetic ancestry map of the world with samples from the Human Genome Diversity Project, each genotyped at 632,958 SNPs[26]. If GWAS array genotypes were available for our samples and for the NHLBI Exome Sequencing Project samples, it would be straightforward to place them directly in this genetic ancestry map. Using targeted sequence data, however, the analysis is more challenging: targeted regions include too few variants to accurately represent global ancestry and off-target regions are covered too poorly, precluding estimation of the accurate genotypes needed for standard principal component analysis. Thus, we relied on the new LASER algorithm (Wang and Abecasis, personal communication) to place each sequenced sample in a pre-defined genetic ancestry map of the world. The method can accurately place individuals on this worldwide ancestry map with <0.05× average coverage of the genome and is thus ideal for targeted sequence data, such as ours and the NHLBI Exome Sequence data, which have average off-target coverage of ~0.23× and ~0.90×, respectively (see Supplementary Figures 1A, 1B, 1E and 1F, which show that PCA coordinates inferred using 0.10× genome coverage or using GWAS array genotypes are highly similar). We focused on samples where PCA coordinates could be estimated confidently (Procrustes similarity larger than 0.95; see Online Methods) and used a greedy algorithm to match cases and controls based on estimated genetic ancestry. As shown in the Online Methods, alternative matching algorithms do not alter our conclusions. After matching, we focused on a set of 2,268 AMD cases and 2,268 controls, ancestry-matched one-to-one (Supplementary Figure 1C and 1G). Since AMD phenotype information was not available for most controls, we expect that a small proportion may eventually develop disease; however, this should not impact power substantially[27]. After matching case-control samples, we excluded 1 variant with Hardy-Weinberg Equilibrium test p-value <10−6 and focused our analysis on 430 protein changing variants in regions that were targeted and deeply sequenced in both experiments as well as far away from insertion deletion polymorphisms.

In this expanded analysis (see Table 1), common variant signals at all loci increased in significance (in comparison to Supplementary Table 4). In addition, two rare coding variants exhibited association with p < 0.01. The first was R1210C in the CFH gene (observed in one control and 23 cases, OR = 23.11, pexact = 2.9×10−6), providing strong support for the original report[16]. The second variant was K155Q in the C3 gene (18 controls, 48 cases, OR = 2.68, pexact = 2.7×10−4; Supplementary Figure 1D and 1H for carrier ancestry distribution). When controlling for a previously described common variant signal nearby, rs2230199 (fcontrol = 20.63%, fcase = 25.26%, marginal pexact = 1.8×10−7, OR = 1.31), the evidence for association with K155Q increased slightly (conditional OR = 2.91, pexact = 2.8×10−5). Inspection of the raw read data shows the variant is well supported and is unlikely to be a sequencing or alignment artifact, a result further confirmed by Sanger sequencing (see Supplementary Figures 2, 3 and 4). Finally, in an examination of our sequenced samples and available whole genome sequences (Online Methods), we observed no additional variants in strong linkage disequilibrium with K155Q that might account for the association signal. Analysis with burden tests, which jointly evaluate evidence for association with rare variants at each gene, identified no significant association signals (Supplementary Figure 5)[28-30].

Table 1

Summary association results for 2,268 sequenced AMD cases and 2,268 sequenced controls. Samples in this expanded analysis include our sequenced AMD samples and genetically matched controls, sequenced by us or by the NHLBI Exome Sequencing Project. The top coding variant in each locus is included in this table when p < 10−6. Rare coding variants are included when the corresponding p-values for the conditional or marginal analysis are less than 1×10−4.

SNP	Chromosome	Position (bp)	Nearest Gene	Consequence	Alleles (ref/alt)	Frequency (alt allele)		OR	P-value	Conditional P-value*
						Cases	Controls
Common variant hits
rs1061170	1	196,659,237	CFH	H402Y	C/T	0.478	0.623	0.555	1.01 ×10−43
rs438999	6	31,928,306	SKIV2L	Q151R	A/G	0.058	0.098	0.566	1.26 ×10−12
rs10490924	10	124,214,448	ARMS2	A69S	G/T	0.329	0.197	1.990	1.04 ×10−45
rs2230199	19	6,718,387	C3	R102G	G/C	0.253	0.206	1.300	1.58 ×10−7
Rare variant hits MAF < 1%, marginal and conditional P <.01 (after conditioning on nearby common variants)
rs121913059	1	196,716,375	CFH	R1210C	C/T	0.005	0.000	23.11	2.9×10−6	6.0×10−4 (rs1061170)
rs147859257	19	6,718,146	C3	K155Q	T/G	0.011	0.004	2.68	2.7×10−4	2.8 ×10−5 (rs2230199)

All p-values were calculated using exact logistic regression. For rare variants, we re-evaluated statistical significance after adjusting for the top common variant in the locus, to avoid shadow signals driven by linkage disequilibrium. The variant used for conditioning is named in (parenthesis).

To confirm the K155Q signal, we genotyped additional samples totaling 4,526 cases and 3,787 controls and, again, observed strong association (fcontrol = 0.5%, fcase = 1.3%, pfollow-up = 7.7×10−7, pcombined = 1.1×10−9, Table 2). In addition, we genotyped 471 families with multiple AMD cases to identify 18 nuclear families where K155Q segregates. These families included 49 affected individuals, where at least one individual carries K155Q and, adjusting for ascertainment, we estimate that 75% of first degree relatives of a K155Q carrier who also have AMD will carry the variant, consistent with an OR of ~3 (Supplementary Table 5 and Online Material). Further strong evidence for association of this variant with macular degeneration is provided in independent work by deCODE Genetics[31], examining 1,143 Icelandic macular degeneration cases and 51,435 Icelandic controls (control frequency 0.55%, OR = 3.45, pdeCODE = 1.1×10−7, pcombined = 1.6×10−15). In 1,606 directly genotyped cases of macular degeneration from the Age Related Disease Study II[32] the variant has frequency 1.77%, similar to our sequenced AMD cases (frequency 1.10%) and our follow-up AMD cases (1.30%) and is notably higher than in our sequenced controls (0.30%), our genotyped controls (0.50%), in NHLBI Exome Sequencing Project participants with primarily European Ancestry (0.40%) and in deCODE controls (0.55%). We found no evidence of the K155Q variant in a small sample of patients with atypical haemolytic-uremic syndrome (aHUS, n=53), a rare disorder whose genetic risk factors partially overlap with macular degeneration.

Table 2

Follow-up Genotyping Summary and Meta-Analysis Summary.

SAMPLE SET	CONTROLS		CASES		P-value
	N	MAF	N	MAF
Discovery Sample
Sequenced Samples (N=4,536)	2,268	.004	2,268	.011	2.7×10−4
Follow-up Samples
Germany / University of Regensburg (N=2,976)	1,147	.006	1,829	.016	1.7×10−3
United States / Vanderbilt/Miami (N=1,819)	726	.004	1,093	.007	3.5×10−1
Netherlands / Rotterdam Study (N=1,409)	1,280	.005	129	.031	1.5×10−4
United Kingdom / Cambridge AMD Study(N=1,279)	423	.006	856	.015	6.2×10−2
United States / UCLA/University of Pittsburgh (N=830)	211	.004	619	.017	8.3×10−4
deCODE Study
deCODE Discovery Sample (N=52,578)	51,435	.005	1,143	*	1.1×10−7
Meta-Analysis
All Follow-Up Samples (N=8,313)	3,787	.005	4,526	.013	7.7×10−7
Discovery and All Follow-Up Samples (N=12,849)	6,055	.005	6,794	.013	1.1×10−9
Discovery, All Follow-Up Samples, and deCODE (N=65,427)	57,490	.005	7,937	*	1.6×10−15

The table includes the number of cases and controls in each comparison, the corresponding allele frequency for the K155Q allele in each set of samples, and the p-value for a comparison of allele frequencies in the two groups. Meta-analysis p-values are calculated using Stouffer’s method.

We next investigated the potential functional consequences of the K155Q variant in silico. Based on protein crystallography, the model in Figure 1 shows that CFH variant R1210C (OR=23.11), C3 variant K155Q (OR=2.91) and C3 variant R102G (OR=1.31) all map near the surface where CFH and C3b interact and suggests they might affect binding of complement factor H to C3b. Factor H inhibits C3b and limits immune responses mediated by the alternative complement pathway. We hypothesize that K155Q and R102G affect binding of the first macro-globular domain of C3 to CFH and thus interferes with inactivation of the alternative complement pathway_ENREF_31, a hypothesis that must be confirmed experimentally[33]. Interestingly, the three variants (R102G and K155Q in C3 and R1210C in CFH) all are associated with replacement of a positively charged residue.

Figure 1

C3 variants R102G and K155Q and CFH variant R1210C are in the interaction domains of the first alpha-macro-globular domain of C3b and CFH, respectively. The fragment of the crystal structure of the four Sushi domains (purple in figure, one not shown for clarity) of CFH in a complex with complement fragment C3b (PDB file: 2wii) was used to explore the effect of disease associated nonsynonymous changes. The CFH residues 987–1230 were used to generate the structure using the first four Sushi domains from 2wii as a structural template (shown in pink, with cysteine residue side chains in yellow). The C-terminal Sushi domains were docked to the binding site in C3b. The first two alpha-macro-globulin domains of C3b, MG-1 and MG-2, are shown in green and cyan, respectively. The location of mutations R102G, K155Q, and R1210C are marked in red.

In summary, our work and the companion paper identify K155Q as a rare C3 variant associated with a ~2.91-fold increased risk of macular degeneration. Together with rare CFH variant R1210C and previously described common C3 variant R102G, K155Q may reduce binding of CFH to C3b, inhibiting the ability of Factor H to inactivate the alternative complement pathway. Clarifying the mechanistic impact of K155Q is likely to be challenging, as illustrated by contradictory results of previous functional follow-up of AMD loci[34-36], but functional studies of complement activity suggest potential next steps[33,37]. Our work relied on targeted sequencing of GWAS loci, genetic ancestry matching of our sequenced samples to additional sequenced controls analyzed with the same variant calling and filtering tools, focused analysis of regions deeply sequenced in both our project and previously sequenced controls, and avoidance of common calling artifacts near insertion/deletion polymorphisms. The use of publicly available samples to augment control sets may be useful to many targeted sequencing studies, but the strictness of matching and variant filtering required for preventing false-positive findings due to population stratification and/or sequence analysis artifacts are areas deserving of further study. As the number of sequenced human genomes and exomes grows, we expect that the utility of the approach will grow – making it possible to match multiple controls to each case and to focus on progressively finer ancestry matches. Our results also emphasize the challenges and the large sample sizes will be required for rare variant studies of complex human traits, as well as the promise of these studies to highlight disease biology, as illustrated by the interaction between Factor H and C3b that is suggested as a key factor in AMD biology here.

Online Methods

Study samples

Macular degeneration cases and controls were recruited at Ophthalmology clinics at the University of Michigan and the University of Pennsylvania and through the Age Related Eye Diseases Study, as previously described. For replication, we contacted members of the International AMD Genetics Consortium; their samples are described in Fritsche et al[13]. All participants provided informed consent allowing for collection of genetic data and all data contributors obtained approval from their local Institutional Review Boards before generating genetic data. Our discovery sample, with ~2350 sequenced cases and ~750 sequenced controls, provides 90% power to discover variants with a frequency of 0.1% and an associated relative risk of 19.2 or greater (similar to CFH R1210C) at significance level alpha = 0.00005, which corresponds to an adjustment for analysis of 1,000 independent coding variants.

Sequence production and quality control

Illumina multiplexed libraries were constructed according to the manufacturer’s protocol (Illumina Inc, San Diego, CA) with modifications: 1) DNA was fragmented using a Covaris E220 DNA Sonicator (Covaris, Inc. Woburn, MA) to range in size between 100 and 400bp. 2) Illumina adapter-ligated library fragments were amplified in four 50μL PCR reactions for eighteen cycles. 3) Solid Phase Reversible Immobilization (SPRI) bead cleanup was used for enzymatic purification and final library size selection targeting 300–500bp fragments. Samples were pooled in groups of 4–24 before hybridization. A custom targeted probe set of 150bp probes was designed (Agilent Technologies, Santa Clara, CA) and captured 0.97 Mb of sequence. The concentration of each captured library pool was determined through qPCR (Kapa Biosystems, Inc, Woburn, MA) to produce cluster counts appropriate for the Illumina GAIIx and HiSeq 2000 platforms. We generated approximately 1.7Gb of sequence per sample, covering 80% of the targeted space at depth >20×. Reads were aligned to the NCBI37/hg19 reference sequence using BWA[38]. Where pre-existing genotype information was available, sample identity was confirmed by comparing sequence data with pre-existing array data.

Quality control and variant calling

Quality control steps for all BAM files included: removal of duplicated reads; recalibration of base qualities[39]; generation of diagnostic graphs and evaluation of sequencing quality[40]; checks for DNA contamination[41]. After removing samples with high contamination, unexpected relatedness or with high discordance rate, we retained 2,335 cases and 789 controls for an initial round of analysis. We calculated the sequencing depth using reads with mapping quality >30 and bases with quality >20. Across the 966,607 base pair target region, we retained an average 123,221,974 bases per individual (127.5× average coverage). Within targeted regions, 98.49% of the protein coding exons had coverage >10×.

We performed variant calling step using UMAKE[23]. Genotype calling and polymorphism discovery were attempted across the original target +/− 50 basepairs. To remove low quality variants, we excluded: 1) sites with average depth <0.5 or >500; 2) sites with evidence of strand bias or cycle bias; 3) sites within 5 basepairs of a 1000 Genomes Project indel; 4) sites with excess heterozygosity. These filters excluded 15,219 low quality variants. The transition-transversion ratio (Ts/Tv) for the remaining 31,527 site was 2.10. Concordance rates between sequence-based genotypes in 13 duplicates were 99.82% when depth >10×. Concordance with array-based genotypes[20] was 98.99% when depth >10×.

59.8% of discovered variants are novel (versus dbSNP 135 and the 1000 Genomes Project). On average, each sample carried 40 synonymous variants, 34 nonsynonymous variants and 1 nonsense variant.

Initial analyses

We first performed single variant association tests using Fisher’s exact test. This analysis confirmed strong association for common variants near CFH, C2, ARMS2 and C3 genes. An initial examination of rare variants suggested some signals were shadows of common variants with larger effects, so we focused on those where association remained significant after accounting for nearby common variants. Conditional signals were evaluated by exact logistic regression[42,43]. Three coding variants had conditional exact P-value less than 0.01 (all also had marginal p-values < 0.01).

Augmenting our sample

We sought ancestry matched controls among samples sequenced in the ESP project. First, we used genome-wide reads to infer sample ancestries on a worldwide population map. Briefly, we first generated a genetic ancestry PCA space using genotyped reference samples (such as those from the Human Genome Diversity Panel). Then, we generated a series of sample specific genetic ancestry PCA that are calibrated to the exact sequencing depth and coverage pattern of each sample and include the reference samples together with a single sequenced sample. Finally, we transformed sample specific PCA coordinates to the original map using Procrustes analysis. This procedure generates a metric (the Procrustes similarity) that summarizes similarity of reference sample placements using array genotypes to placements using sequence data and we only considered samples where this metric was >0.95 as candidates for matching. Second, we used a procedure inspired on propensity score matching to pair cases and controls[44]. Briefly, this procedure uses logistic regression to predict the probability that an individual is a case using the four principal components of ancestry as predictors and disease status as the outcome. This estimated probability of being a case for each sample is a propensity score and can be used to match cases and controls. For matching, we used a greedy algorithm to match cases and controls; allowing matches when the respective propensity scores differed by <.0001. An alternative matching algorithm that matched cases and controls mapping close together in principal component space according to the Euclidean distance between them gave similar results (association at K155Q had OR=2.68, exact p-value 4.5×10−5 using Fisher’s exact test).

To avoid variant calling artifacts, we applied very stringent filters to both the AMD study and ESP study call sets. For both studies, we examined only sites with call rates >90%, Phred-scaled variant quality scores >30, passing all study specific quality control filters, with depth >10× for >90% of the samples in the AMD or ESP callsets, and >5-bp from a 1000 Genomes Project indel. Primers used to confirm the presence of K155Q by Sanger Sequencing are given in Supplementary Table 6.

Analyses using the combined AMD and ESP data set

Similar to our initial analysis, we first applied Fisher’s exact test for association to all variants. Next, we examined variants with frequency <1% for which signal remained significant after adjusting for common variants. This analysis highlighted R1210C in CFH and K155Q in C3 (Figure 1).

Linkage disequilibrium analysis

To search for variants that might explain the signal at K155Q, we evaluated linkage disequilibrium between K155Q and all variants within 1 Mb both within the samples sequenced for this experiment and also in preliminary whole genome sequence data for 600 individuals (300 macular degeneration cases, 300 controls; Swaroop, Stambolian and Abecasis, personal communication). This analysis did not find variants in strong linkage disequilibrium in the nearby region. The variant is only present in one 1000 Genomes Project sample, which does not allow for reliable estimates of linkage disequilibrium.

Segregation analysis

In a segregation analysis, one identifies probands who carry K155Q and then evaluates the probability that they transmit the variant to affected relatives (under the null, we would expect to find the variant in 50% of first degree relatives of a carrier). We genotyped 471 pedigrees with multiple affected individuals. In each pedigree where K155Q was found in >1 affected individual, we selected the nuclear family with the largest number of affected individuals. We recorded the number of affected individuals (N) and the number of K155Q carriers (C). Then, to average over possible choices of proband, we assigned each family a weight of C/N (this is the probability that a randomly selected proband in the family carries K155Q) and then scored the number of affected first degree relatives (N-1) and of carriers among those (C-1). The estimated fraction of carriers among affected first degree relatives of a proband is then calculated by summing C/N * (C-1) and C/N * (N-1) over families and taking the ratio of the two quantities.