Identification of a Rare Coding Variant in Complement 3 Associated with Age-related Macular Degeneration

Genetic and environmental factors contribute to age-related macular degeneration (AMD)^1,2, a major cause of vision loss in elderly individuals³. 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¹⁵. 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¹⁶. 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²⁰ (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²² (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²⁰, 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⁻⁶, including all discovered variants, or p < 0.05 / 1,422 = 3.5×10⁻⁵ 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 f_case = 12%, control frequency f_control = 27%, P-value = 2.52×10⁻⁴⁸), ARMS2 (rs10490924, f_case = 33%, f_control = 18%, P-value = 5.48×10⁻²⁷), C3 (rs2230199, f_case = 25%, f_control = 17%, P-value = 3.94×10⁻⁹) and C2/CFB (rs556679, f_case = 7%, f_control = 12%, P-value = 1.32×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)²³, 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²⁴. 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²⁶. 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²⁷. After matching case-control samples, we excluded 1 variant with Hardy-Weinberg Equilibrium test p-value <10⁻⁶ 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, p_exact = 2.9×10⁻⁶), providing strong support for the original report¹⁶. The second variant was K155Q in the C3 gene (18 controls, 48 cases, OR = 2.68, p_exact = 2.7×10⁻⁴; Supplementary Figure 1D and 1H for carrier ancestry distribution). When controlling for a previously described common variant signal nearby, rs2230199 (f_control = 20.63%, f_case = 25.26%, marginal p_exact = 1.8×10⁻⁷, OR = 1.31), the evidence for association with K155Q increased slightly (conditional OR = 2.91, p_exact = 2.8×10⁻⁵). 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.

To confirm the K155Q signal, we genotyped additional samples totaling 4,526 cases and 3,787 controls and, again, observed strong association (f_control = 0.5%, f_case = 1.3%, p_follow-up = 7.7×10⁻⁷, p_combined = 1.1×10⁻⁹, 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³¹, examining 1,143 Icelandic macular degeneration cases and 51,435 Icelandic controls (control frequency 0.55%, OR = 3.45, p_deCODE = 1.1×10⁻⁷, p_combined = 1.6×10⁻¹⁵). In 1,606 directly genotyped cases of macular degeneration from the Age Related Disease Study II³² 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.

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³³. Interestingly, the three variants (R102G and K155Q in C3 and R1210C in CFH) all are associated with replacement of a positively charged residue.

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.