The genetic architecture of membranous nephropathy and its potential to improve non-invasive diagnosis

Membranous Nephropathy (MN) is a rare autoimmune cause of kidney failure. Here we report a genome-wide association study (GWAS) for primary MN in 3,782 cases and 9,038 controls of East Asian and European ancestries. We discover two previously unreported loci, NFKB1 (﻿rs230540, OR = 1.25, P = 3.4 × 10−12) and IRF4 (﻿rs9405192, OR = 1.29, P = ﻿1.4 × 10−14), fine-map the PLA2R1 locus (﻿rs17831251, OR = 2.25, P = 4.7 × 10−103) and report ancestry-specific effects of three classical HLA alleles: DRB1*1501 in East Asians (OR = 3.81, P = 2.0 × 10−49), DQA1*0501 in Europeans (OR = 2.88, P = 5.7 × 10−93), and DRB1*0301 in both ethnicities (OR = 3.50, P = 9.2 × 10−23 and OR = 3.39, P = 5.2 × 10−82, respectively). GWAS loci explain 32% of disease risk in East Asians and 25% in Europeans, and correctly re-classify 20–37% of the cases in validation cohorts that are antibody-negative by the serum anti-PLA2R ELISA diagnostic test. Our findings highlight an unusual genetic architecture of MN, with four loci and their interactions accounting for nearly one-third of the disease risk.

M embranous Nephropathy (MN) is a rare cause of kidney failure, manifesting as nephrotic syndrome with a peak incidence between 30 and 50 years of age 1 . The landmark discoveries of pathogenic antibodies against neutral endopeptidase in antenatal MN 2 , and anti-phospholipase A2 receptor (PLA2R) antibodies in adult MN 3 have established MN as the disease of autoantibodies directed against podocyte antigens. Several studies have confirmed the presence of autoantibodies against PLA2R in~60-70% of cases of primary MN 4 , with another 3-5% potentially explained by antibodies against thrombospondin type 1 domain-containing 7A 5 .
Previous genome-wide association study (GWAS) for MN conducted in 75 French, 146 Dutch and 335 British cases genotyped with low resolution arrays identified impressively strong associations of the HLA region and the PLA2R1 locus encoding the dominant antigen in MN 6 . These findings suggest that genetic variation controls the immunogenicity and/or expression level of the PLA2R auto-antigen, as well as the production of anti-PLA2R autoantibodies in individuals with a permissive HLA haplotype. However, specific causal alleles underlying GWAS associations have not yet been mapped at high resolution. Moreover, prior GWAS was limited to Europeans, and the reported associations have not been examined comprehensively across different ethnicities. Lastly, because of small sample size, the prior study might have missed additional disease relevant loci.
Herein, we report a genetic study of primary MN involving 12,820 individuals (3782 biopsy-documented cases and 9038 ancestry-matched controls), across nine cohorts of East Asian and European ancestries. The composition of our cohorts reflects the demographics of the centres that have collected DNA samples for genetic studies of this rare disease over the past 15 years. By using high resolution arrays with genome-wide imputation and over 7fold increase in sample size compared to the prior GWAS, we discover two previously unreported genome-wide significant risk loci for MN and perform high resolution mapping and ethnicityspecific analyses of the known loci.
We describe an unusual genetic architecture of MN, with four loci and their genetic interactions accounting for nearly one-third of the disease risk. Our study implicates dysregulation of NFKB1 and IRF4 genes in the disease pathogenesis, providing genetic support for potential targeting of the NF-κB and interferon signalling pathways in primary MN. We also refine ethnicity-specific effects at the HLA locus, defining DRB1*1501 as a major risk allele in East Asians, DQA1*0501 in Europeans, and DRB1*0301 in both ethnicities. We describe a risk haplotype at the PLA2R1 locus that has a regulatory function and exhibits strong genetic interactions with the HLA-DRB1 risk alleles. Lastly, we calculate a genetic risk score (GRS) based on these findings which, when used in combination with a serum anti-PLA2R ELISA (a serologic test for MN currently in clinical use), shows superior performance in discriminating cases and controls than the ELISA or GRS alone. We validate the performance of this combined risk score (CRS) in external validation cohorts. Our results demonstrate that a combined serum-genetic test can potentially be used to establish a new diagnosis of primary MN, obviating the need for a high risk kidney biopsy procedure in the majority of cases.

Results
Study design. Our study involved nine case-control cohorts, including four East Asian cohorts of 4841 individuals (1632 primary MN cases and 3209 controls) and five European cohorts of 7979 individuals (2150 primary MN cases and 5829 controls). Eight cohorts were genotyped with high density SNP arrays, imputed using the latest whole genome sequence reference panels, and meta-analyzed genome-wide, and the top 46 loci selected based on P < 5 × 10 −5 were tested by targeted genotyping in the ninth cohort (Supplementary Table 1). The summary of study cohorts, genotyping methods, and ancestry-specific imputation panels is provided in Table 1.
All cases used in this study were defined by a kidney biopsy diagnosis of idiopathic MN and any suspected secondary cases due to drugs, malignancy, infection, or autoimmune disease were excluded. With the exception of the German Chronic Kidney Disease (GCKD) cohort, all controls used for discovery involved healthy population controls and any individuals with a known diagnosis of kidney disease were excluded. The GCKD cohort was drawn entirely from a prospective observational study of patients with CKD and consisted of biopsy-defined cases and controls for whom CKD etiology was clearly assigned to a non-MN cause, as previously described 7 .
All genome-wide significant loci (P < 5 × 10 −8 ) were refined by cohort-stratified stepwise conditional analyses to define independently associated haplotypes. We also analyzed classical HLA alleles and all common amino acid polymorphisms at class I and class II genes imputed at high resolution. We performed detailed genomic annotations and explored epistatic effects for significant loci. Based on significant GWAS loci, we designed a GRS for MN Table 1 Baseline characteristics of participants in the discovery and replication cohorts. ARTICLE and performed its validation in external cohorts, including the three previously published European GWAS cohorts 6 and the Nephrotic Syndrome Study Network (NEPTUNE) study 8 . The descriptions of all cohorts including ancestry analyses and details of statistical approaches are provided in Methods, Supplementary Methods, and Supplementary Figs. 1 and 2.
Genome-wide association. The results of combined genomewide meta-analyses are summarized in Fig. 1 and Table 2, with  more information provided in Supplementary Table 1    . We also confirmed strong and highly significant associations at the previously described loci, including chromosome 2q24.2 encoding PLA2R1 (rs17831251, OR = 2.25, Meta-analysis P = 4.7 × 10 −103 ) and 6p21.32 encoding HLA-DQA1/DRB1 genes (rs9271573, OR = 2.41, Metaanalysis P = 2.7 × 10 −154 ). Conditional analyses of the three non-HLA loci revealed that each signal is explained by a single SNP in each cohort, suggesting a single shared risk haplotype per locus in East Asian and European populations (Fig. 1b-d). To further test if the causal variants at these loci are likely shared between Europeans and East Asians, we performed 99% credible set analyses using summary statistics for each ancestry-defined subgroup and compared them with credible sets derived from the trans-ethnic meta-analysis. We confirmed that the predicted causal variants derived from the trans-ethnic analysis were largely overlapping with ancestry-specific results ( Supplementary Fig. 5). In contrast, stepwise conditional analyses of SNPs at the HLA region revealed a complex pattern of association, with at least three independently genome-wide significant SNPs explaining the signal across all cohorts (Supplementary Table 2).
Given the complexity of the association signal at the HLA locus and known differences in linkage disequilibrium (LD) patterns by ancestry, we performed additional analyses of this region separately in East Asians and Europeans. In the conditional analyses of the East Asian cohorts, only two independently associated SNPs explained the entire signal at this locus (rs9269027 and rs1974461). In Europeans, stepwise conditional analyses revealed three independently associated genome-wide significant SNPs (rs9271541, rs9265949, and rs2858309), suggesting a more complex pattern of association (Supplementary  Table 2). In both ethnicities, the top signal centred on HLA-DRB1 and DQA1 genes (Fig. 2a, b).
Classical HLA alleles and amino acid polymorphisms. We next imputed classical HLA alleles at two-and four-digit resolution using ethnicity-specific reference panels (see Methods). The first two digits specify a group of HLA alleles known as super-types as defined by older typing methodologies. The third through fourth digits specify nonsynonymous substitutions. Moreover, we imputed individual amino acid polymorphisms at class I (HLA-A, -B, and -C) and class II (HLA-DQB1, -DQA1 and -DRB1) genes.
In East Asian cohorts, stepwise conditioning on classical HLA alleles defined two independent risk alleles, DRB1*1501 (OR = 3.81, Wald test P = 2.0 × 10 −49 ) and DRB1*0301 (OR conditioned = 3.88, Wald test P = 4.5 × 10 −24 , Fig. 2c, Supplementary Table 3). In the analysis of polymorphic amino acid sites, genetic variation at only two codons encoding residues at positions 13 and 71 in DRβ, explained the entire HLA-DRB1 signal (Fig. 3a Table 5). Consistent with a prior study in Chinese patients 9 , these amino acids define the classical risk alleles DRB1*1501 and DRB1*0301 (Supplementary Table 6), and their side chains map adjacent to each other within the antigen-binding pocket of the β-chain of DR (Fig. 3c).
Given that our HLA imputation reference panels were considerably smaller for East Asians compared to Europeans, we sought additional validation of the observed classical HLA associations that were Asian-specific. We therefore created another reference panel based on the MHC sequence data from Zhou et al. 10 including 10,689 control individuals of Han Chinese ancestry. Using SNP2HLA software, we then re-imputed classical HLA alleles for our East Asian cohorts. We used the same quality control filters (MAF > 0.01 and imputation R 2 >0.8) and methods for association testing as described above. We observed no major differences in the association statistics for the two Asian risk alleles, DRB1*1501 (OR = 3.49, P = 3.85e−40) and DRB1*0301 (OR = 4.08, P = 6.3E−24), demonstrating that these effects do not represent artifacts of smaller imputation panels.
We next performed the analysis of HLA amino acid substitutions in Europeans. Consistent with the association analyses of classical alleles, five bi-allelic sites in DQA1 that correlate with the DQA1*0501 allele were most strongly associated with the risk of MN (75Ser-107Ile-156Leu-161Glu-163Ser, Wald test P = 5.7 × 10 −93 , Fig. 3b Tables 9 and 10). Notably, positions 74 (Europeans) and 71 (East Asians) are separated by a single turn along the α-helix, and their side chains are spatially close to that of position 13, located on the beta-sheet floor with its side chain oriented into the peptidebinding groove (Fig. 3c).
To confirm ethnicity-specific HLA effects, we repeated stepwise conditioning in a joint stratified analysis of all cohorts using biallelic tests of HLA alleles with formal tests of heterogeneity (Supplementary Table 11, Fig. 4a). The top classical allele supported by all cohorts regardless of ethnicity was DRB1*0301 PLA2R1 locus and its genetic interactions. Consistent with prior GWAS, the most significant non-HLA locus resided on chromosome 2q24.2 6 . The top SNP was in the first intron of PLA2R1, which encodes the main podocyte autoantigen in primary MN. This signal was supported by both ethnicities, but the effect appeared stronger in East Asians (OR = 2.81, Meta-analysis P = 3.5 × 10 −61 ) compared to Europeans (OR = 1.98, Meta-analysis P = 4.7 × 10 −48 , Table 2). After conditioning the association on the top SNP, rs17831251, there was no residual association at this locus, suggesting a common risk haplotype in both ethnicities (Fig. 1b).
We next refined the previously reported genetic interactions between the PLA2R1 locus and HLA risk haplotypes. The PLA2R1 risk genotype exhibited significant multiplicative interaction with both Asian and European HLA risk haplotypes (Fig. 4b, c)  (Interaction test P = 0.1). Similarly, there was no significant interaction between DQA1*0501 allele and PLA2R1 locus in East Asians. These analyses suggest that the PLA2R1 locus interactions are driven primarily by the DRB1 alleles. We annotated all SNPs in LD with rs17831251 for potential impact on the structure and/or transcriptional regulation of PLA2R1. We found two common missense variants in moderate LD, rs35771982 (p.H300D, r 2 = 0.69) and rs3749117 (p.M292V, r 2 = 0.68), but the effects of these variants were considerably weaker compared to rs17831251, suggesting that they are unlikely to represent causal variants (Supplementary Table 12). Our tissue-specific functional scoring method for non-coding variants based on the ENCODE and Roadmap Epigenetics data 11 prioritized another variant in intron 1, rs17241973 (r 2 = 0.93 with rs17831251) that intersects a putative enhancer element across multiple tissues ( Supplementary Fig. 6). Both rs17831251 and rs17241973 exhibit strong cis-eQTL effects on PLA2R1 expression wherein the MN risk alleles associate with lower mRNA expression of PLA2R1 across multiple tissues in GTEx 12 , but this effect appears reversed for the kidney tissue (Supplementary Fig. 7). To further confirm these kidney-specific effects, we used gene expression data from manually micro-dissected human kidney compartments of 166 NEPTUNE participants 13 . We detected suggestive glomerular eQTL effects that were weak, but direction-consistent with GTEx for rs17831251 (Wald test P = 0.055) and rs17241973 (Wald test P = 0.024), wherein MN risk allele were associated with increased glomerular PLA2R1 mRNA levels ( Supplementary Fig. 8).
Because kidney tissue compartments are not well represented in either ENCODE or Roadmap datasets, we next examined the genomic location of rs17831251 and rs17241973 in relationship to the recently published kidney compartment-specific chromatin landscape 14 (Supplementary Fig. 9); rs17241973 lies within intron 1 of PLA2R1 in open chromatin that is active in both glomerular and tubular compartments and is contiguous with the gene promoter. In contrast, rs17831251 lies within a broad region of increased chromatin accessibility in glomeruli, and is only 2.1-kb away from a glomerulus-specific DHS that contains a highconfidence NFKB1-binding motif. In glomerulus-specific chromatin conformation (Hi-C) data, both rs17831251 and rs17241973 make regional and distal contacts with other glomerular DHS, emphasizing a composite cis-regulatory module for PLA2R1 gene expression.
Novel loci encoding NFKB1 and IRF4. The 4q24 locus contains the NFKB1 gene, which encodes an active DNA binding subunit of the NF-κB transcriptional complex. The top SNP, rs230540 (OR = 1.25, Meta-analysis P = 3.4 × 10 −12 ), is an intronic variant predicted to have a functional effect specific to immune cells ( Supplementary Fig. 10). In agreement with our prediction, rs230540 has been associated with higher mRNA expression of NFKB1 in whole blood (P = 2.6 × 10 −11 ) 15 and in CD4 + T cells (P = 2.0 × 10 −9 ) 16 . Consistent with pro-inflammatory effects of NF-κB, the MN risk haplotype at this locus determined higher leukocyte counts 17 and increased risk of ulcerative colitis 18,19 and primary biliary cholangitis 20,21 (Supplementary Table 13, Fig. 5a). NFKB1 is also expressed in human podocytes 22 , as well as primary human glomerular and tubular epithelial cell cultures 14 , and rs230540 intersects an active glomerular DHS with several     Fig. 11). Notably, the MN risk allele has previously been associated with lower estimated glomerular filtration rate in GWAS of renal function 23,24 , thus this locus may be more broadly associated with the risk of kidney disease. The top SNP on chromosome 6p25.3, rs9405192 (OR = 1.29, Meta-analysis P = 1.4 × 10 −14 ), resides upstream of IRF4 gene, which belongs to the family of transcription factors regulating interferon-inducible genes. IRF4 is lymphocyte specific and negatively regulates Toll-like-receptor signalling that is central to the activation of innate immune system; this gene is known to be under the transcriptional control of the NF-kB complex [25][26][27] . Unlike PLA2R1 and NFKB1, IRF4 does not appear to be expressed in human kidney cells by single nuclei RNA-seq 28 ( Supplementary Fig. 12). We did not find functional or coding SNPs in LD with rs9405192, nor did we observe any cis-eQTL effects for this variant, thus the precise mechanism underlying this association remains unknown. However, the analysis of binding sites for individual components of the NF-κB complex in lymphocytes 29 suggested binding of the complex in close proximity of rs9405192 ( Supplementary Fig. 13). The risk allele at this locus is also in strong LD with variants previously  associated with increased risk of inflammatory bowel disease 19 , and in weaker LD with several risk variants for chronic lymphocytic leukemia (Supplementary Table 13, Fig. 5a), suggesting the pattern of pleiotropy that is similar to the NFKB1 locus. Nevertheless, we detected no statistically significant genetic interactions of IRF4 and NFKB1 loci.
In addition, we systematically annotated all other suggestive non-HLA loci defined by P < 5.0 × 10 −5 and these results are summarized in Supplementary Table 14. To enhance potential genetic discovery of novel podocyte antigens, we also repeated genome scans after conditioning for the PLA2R1 locus, but detected no additional suggestive loci.
SNP-based heritability and risk explained by GWAS. Using our genotype data and genome-based restricted maximum likelihood method (GREML) 30 , we estimated the overall SNP-based heritability of MN at 0.43 (SE = 0.039) in East Asians and 0.36 (SE = 0.0046) in Europeans. Remarkably, all genome-wide significant risk alleles exhibited unusually large effect sizes for GWAS. In order to quantify the fraction of disease variance cumulatively explained by genome-wide significant SNPs and their interactions, we performed ethnicity-specific GRS analyses (see Methods). Each GRS was expressed as a weighted sum of risk alleles with weights defined by their mutually adjusted effect estimates and included the 3 independent non-HLA SNPs (rs6707458, rs230540, rs9405192) as well as ethnicity-specific HLA risk alleles and their interactions. This included rs9269027, rs1974461, and rs9269027*rs6707458 interaction term for East Asians, and rs9271541, rs9265949, rs2858309 and rs9271541*rs6707458 interaction term for Europeans (Supplementary Table 15). The GRS calculated using this method explained 32% disease risk in East Asians, 25% in Europeans, and 29% of overall disease risk across all cohorts combined. Remarkably, the magnitude of the GRS effect was comparable to rare, highly penetrant mutations causing Mendelian forms of kidney disease, with individuals in the top decile of GRS having 30 to 40-fold higher disease risk compared to the lowest decile (Fig. 6a, b).
Clinical correlations of the GRS. For a subset of patients with available clinical data, we performed genetic correlation analyses with selected clinical features reflective of disease severity. The GRS was positively correlated with PLA2R antibody seropositivity (Wald test P = 9.0 × 10 −8 ), and in those with detectable antibodies, higher titers at the time of biopsy (Slope test P = 1.2 × 10 −9 , Fig. 5b). The GRS also predicted worse proteinuria at the time of biopsy, which represents the key marker of MN severity and prognosis (Slope test P = 1.3 × 10 −3 , Fig. 5c). Other clinical features, such as age at diagnosis, renal function, or serum albumin levels at the time of biopsy were not significantly correlated with the GRS after multivariate adjustment (Supplementary Table 16).
Potential diagnostic implications of the GRS. Although the diagnosis of MN is traditionally established by a kidney biopsy, the detection of circulating PLA2R antibodies by ELISA has recently emerged as a useful diagnostic modality 31 . In this study, we performed ELISA in sera obtained within 6 months of a diagnostic kidney biopsy in a total of 2331 individuals (1488 cases, 300 healthy controls, and 543 disease controls). In East Asians, we estimated that the standard ELISA cut-off of 20 U/mL provided 100% specificity and 60% sensitivity for the diagnosis of MN. In the analysis of Europeans, depending on the specific cohort, the same cut-off provided 99-100% specificity and 51-57% sensitivity (Supplementary Table 17). While the antibody level of 20 U/mL represents the manufacturer's recommended cut-off, levels 2-20 U/mL are frequently considered as borderlinenegative, and levels <2 U/mL as negative 31 . In our cohorts, the cut-off 2 U/mL had inadequate diagnostic specificity (range 73-92%). These results confirm the key limitation of the PLA2R antibody ELISA, which has high specificity (99-100%) but low sensitivity (51-60%) at the standard recommended cut-off point; while lowering the cut-off increases sensitivity, it results in inadequate specificity. Consequently, the levels in the borderlinenegative range (2-20 U/mL) are difficult to interpret clinically.
Given this limitation, we evaluated if the addition of genetic risk information can improve the performance of ELISA, especially in cases that fall in the borderline-negative range. First, we evaluated diagnostic properties of the GRS alone in our discovery cohorts. In East Asians, the genetic test had area under the receiver operating characteristics curve (AUROC) of 0.80 (95% CI: 0.78-0.82), while in Europeans the AUROC was 0.75 (95% CI: 0.74-0.77). Combining genetic and serologic tests in the form of a CRS provided superior case discrimination with AUROCs of 0.96 (95% CI: 0.95-0.98) in East Asians and 0.89 (95% CI: 0.87-0.91) in Europeans (Fig. 6, Supplementary Table 18).
We next tested the GRS performance in several external validation cohorts, including three independent GWAS cohorts of European ancestry as well as in the European-American NEPTUNE participants with incident nephrotic syndrome (Supplementary Table 18C-E). Overall, the effects and the diagnostic performance of the GRS were comparable between the European discovery and each validation cohort (Fig. 6c). When combined with the serum antibody titer, the CRS achieved AUROC of 0.96 (95% CI: 0.94-0.97) across all validation cohorts combined (Supplementary Table 18E, Fig. 6f).
In the subgroup analyses, we compared the diagnostic properties of GRS and CRS by the antibody status in all cohorts pooled by ancestry (Fig. 7). These analyses demonstrated that both GRS and CRS were predictive of case status even for antibody-negative MN. Importantly, the CRS continued to have excellent performance in classifying the borderline-negative cases (antibody level 2-20 U/mL range), with AUROCs of 0.98 (95% CI: 0.97-0.99) in East Asians and 0.95 (95% CI: 0.93-0.96) in Europeans. Notably, among all cases for which the ELISA test was either negative or inconclusive, adding genetic information in the form of CRS can establish the diagnosis in 20-37% cases with 99% specificity. The comparison of AUROCs between GRS, CRS, and serum anti-PLA2R Ab test by ancestry is provided in Supplementary Fig. 14, and the clinical implications of these findings are summarized in Supplementary Note 1 and Supplementary Table 19.
Lastly, we expanded our validation studies to non-European participants of the NEPTUNE study (Supplementary Tables 20  and 21). Although the European risk score performance was diminished in Hispanic Americans, the European GRS performed well in African Americans, and this is despite substantial differences in risk allele frequencies between Europeans and Africans (Supplementary Table 22). Similar to the European validation cohorts, the European GRS was superior compared to the trans-ethnic GRS when applied to the NEPTUNE minority populations, while the Asian GRS had relatively poor performance in both African-American and Hispanic/Latino cohorts (Supplementary Table 21, Supplementary Fig. 15).

Discussion
Our study provides important insights into an autoimmune disease and the genetic architecture of MN. First, we discover novel genome-wide significant risk loci for MN with large effects encoding two transcriptional master regulators of inflammation, NFKB1 and IRF4. The association at the NFKB1 locus highlights the role of the canonical NF-κB pathway in primary MN. Upon activation by pro-inflammatory signals, NFKB1 undergoes proteasome processing to p50, an active DNA binding subunit of the NF-κB complex. Inappropriate activation of this pathway has previously been studied in progressive diabetic nephropathy 32 , and in the context of inflammatory diseases, including IBD 33,34 and MN [35][36][37] . Importantly, the MN risk allele at this locus has a concordant effect on the risk of ulcerative colitis 18,19 and primary biliary cirrhosis 20,21 . It has also been associated with increased mRNA expression of NFKB1 in cis, and reduced DNA  methylation in trans across >400 CpGs that overlap with NF-κBbinding sites, suggesting enhanced baseline activity of NF-κB 38 .
The NF-κB complex is known to up-regulate IRF4 expression with cross-regulatory feedback loops between NFKB1 and IRF4 described in several prior studies [25][26][27] . Taken together, NFKB1 and IRF4 loci participate in a common regulatory pathway in immune cells, and our genetic findings clearly establish a critical role of this pathway in the pathogenesis of MN. Second, due to the bi-ethnic composition of our cohorts, we were able to refine ethnicity-specific effects at the HLA locus, defining DRB1*1501 as a major risk allele in East Asians, DQA1*0501 in Europeans, and DRB1*0301 in both ethnicities.
These findings suggest that different epitopes are likely presented to T cells to initiate the anti-PLA2R response in East Asians and Europeans. We also identified specific high-risk amino acid substitutions, at positions 13, 71, and 74, mapping to the P4 pocket of DRβ1. Although the same positions contribute to the risk of T1D 39 and rheumatoid arthritis 40 , the effects of individual residues at each position are discordant, likely reflecting differences in target epitopes.
Third, we confirm that a single haplotype at the PLA2R1 locus conveys the disease risk in both East Asians and Europeans, and exhibits genetic interactions with HLA-DRB1 risk alleles. Our analysis supports a regulatory function of the PLA2R1 risk  haplotype. The candidate causal variant resides in the first intron of PLA2R1 and intersects a predicted enhancer element. While this variant is normally associated with suppressed PLA2R1 transcription across multiple tissues, it appears to increase expression of PLA2R1 in the kidney. This finding highlights the importance of studying target tissues and is consistent with the findings that among CKD loci that are transcriptionally active in renal tissue, 15.8% of effects are kidney-specific 41 . Notably, the top variants at the PLA2R1 locus also intersect a putative NF-κB binding site in lymphocytes, although no similar data is presently available for podocytes. Further experimental work is thus needed to test if the glomerular-specific eQTL effect is under the transcriptional control of NF-κB. Moreover, larger glomerular compartment-specific datasets will be needed to confirm the observed eQTL effects.
Another observation is that all four genome-wide significant risk loci (PLA2R1, IRF4, NFKB1, and HLA) exhibit highly pleiotropic effects and all four lead SNPs have a concordant effect on the risk of inflammatory bowel disease (IBD). This observation suggests shared pathogenic mechanism between IBD and MN. Considering that MN is an orphan disease without a targeted treatment, there may now be opportunities for drug repositioning approaches from IBD, where several new antiinflammatory agents are currently under development. Our study suggests that the NF-κB and interferon pathways may represent particularly attractive drug targets.
Remarkably, our GWAS loci are highly predictive of the disease status and jointly explain up to one third of disease risk, an exceptionally large fraction for common alleles. This may be partially explained by the fact that MN frequently occurs after the peak reproductive age, allowing the risk alleles to escape purifying selection. Moreover, even though the risk alleles are common, our interaction analysis demonstrates that specific high-risk genotype combinations are relatively rare in the general population, potentially explaining the low overall prevalence of MN 42 . The alternative hypothesis is that of balancing selection. NF-kB and IRF4 are both involved in immune defenses against common pathogens and some Phospholipase A2 ligands for PLA2R1 represent downstream NF-kB targets with antibacterial properties 43,44 . Therefore, the observed high frequencies of MN risk alleles could be explained by their protective effects against common infections.
Finally, a simple GRS based on our GWAS loci has excellent discriminant properties when combines with anti-PLA2R ELISA test. Importantly, the combined genetic-serum test has superior diagnostic properties compared to serologic test alone, mitigating the key issue of low sensitivity. The GRS provides complementary information to the serum test and correctly reclassifies 20-37% of antibody-negative cases, potentially sparing the need for a kidney biopsy in this large subgroup of patients. In the clinical settings where neither a serum test nor a kidney biopsy is possible, the GRS itself can establish a diagnosis of MN with 99% specificity in 13-15% of cases. The practical advantage of this approach is that the GRS can be readily determined at any time after birth and, unlike the serum test, it does not fluctuate with time or in relationship to the disease onset, activity, or treatment. One important limitation, however, is that genetic effects may be population-specific and may not be generalizable to populations not represented in our GWAS. The performance of the GRS is remarkably consistent in our discovery and validation cohorts, including African-Americans, but it appears to be lower in self-reported Latino/ Hispanics. Therefore, future efforts extending GWAS for MN to more diverse populations will be important.
In summary, we described a highly unusual genetic architecture of MN, including large effect sizes for a small number of common alleles and a strong evidence for ethnicity-specific genetic interactions. These insights enabled formulation a powerful genetic disease predictor that provides means to enhance a non-invasive diagnosis of MN, and can be especially useful in the settings where kidney biopsy represents too great of a risk or is not readily available.

Methods
Study design overview. We performed a genome-wide meta-analysis of eight discovery cohorts of East Asian and European ancestry (2956 cases and 7799 controls), all genotyped with high resolution arrays and imputed to~7 million common high-quality markers using ancestry-matched reference panels. The top signals from the meta-analysis (P < 5 × 10 −5 ) were typed in the additional East Asian replication cohort of 826 cases and 1239 controls. Subsequently, all cohorts (3782 cases and 9038 controls) were analyzed jointly to define genome-wide significant signals. All subjects provided informed consent to participate in genetic studies, and the Institutional Review Board of Columbia University as well as local ethics review committees for each of the individual cohorts approved our study protocol. The individual cohorts, genotyping methods, and quality control analyses are described in the Supplementary Methods.
Primary association analyses and genome-wide meta-analyses. Within each cohort, primary association scans were performed for markers that were common (MAF > 0.01) and imputed at high quality (r 2 > 0.8) using logistic regression under additive coding of dosage genotypes, and with adjustment for cohort-specific significant principal components (PCs) of ancestry. To quantify potential inflation of type I error due to stratification or technical artifacts, we estimated genomic inflation factors 45 for each genome-wide scan after excluding HLA and PLA2R loci. No substantial inflation was noted in any individual scan (lambda consistently <1.05 for each individual cohort). Subsequently, a fixed effects meta-analysis was performed to combine the results of the eight discovery cohorts using METAL 46 . Genome-wide distributions of P-values were examined visually using quantilequantile plots for each individual cohort as well as for the combined analysis. The final meta-analysis quantile-quantile plot showed no global departures from the expected null distribution (Supplementary Fig. 3), with the genomic inflation factor estimated at 1.03 for the overall meta-analysis. Suggestive signals were defined by P-value < 5.0 × 10 −5 . To declare genome-wide significance of a novel locus, we used the generally accepted P-value threshold of 5.0 × 10 −8 .
Conditional analyses. To detect additional independent SNPs at genome-wide significant loci, we performed stepwise conditional analyses of each locus using logistic regression. This was done by including the genotype of conditioning SNP (s) under additive coding as covariate(s) in the outcome model. The conditional analyses were performed individually within each cohort and with adjustments for cohort-specific ancestry PCs. Subsequently, the conditioned summary statistics were combined across cohorts using fixed effects meta-analysis, similar to our primary association analyses.
Credible set analyses. For each of the three genome-wide significant non-HLA loci, we derived 99% credible sets using the trans-ethnic and ethnicity-specific meta-analysis results. First, we derived approximate Bayes factors from GWAS association statistics using Wakefield's formula, as implemented in the R package gtx 47 . Using CAVIAR software 48 , we next calculated the posterior probability (PP) for each SNP driving the association signal at each locus. We assumed there was only a single causal variant at each locus, since no additional independent SNPs were detected on stepwise conditioning analyses. We derived both trans-ethnic as well as ethnicity-specific 99% credible sets based on ranking the variants by their PPs and adding the variants to the set until cumulative PP > 99% was reached for each region. The overlaps between ethnicity-specific and trans-ethnic analyses were visualized in Supplementary Fig. 5.
HLA imputation. Six discovery cohorts (Chinese, South Korean, Japanese, European-1, European-2, and Turkish) had primary genotype data available for HLA imputation and association testing. For each of these cohorts, we imputed classical HLA alleles at two-and four-digit resolution, as well as individual amino acid polymorphisms at class I (HLA-A, -B, and -C) and class II (HLA-DQB1, HLA-DQA1 and HLA-DRB1) loci using SNP2HLAsoftware 49 . The European cohorts and the East Asian cohorts were imputed separately, using ethnicity-specific reference panels. For European reference, we used the pre-phased HLA reference dataset generated by the Type 1 Diabetes Genetics Consortium (T1DGC, 5,225 individuals) 49 . For our East Asian cohorts, we used the Pan-Asian HLA Reference Panel (268 individuals) 50 . For validation of classical HLA association results in East Asians, we built additional East Asian reference panel based on the MHC sequence data from Zhou et al. (10,689 Han Chinese) 10 . In the association analyses, we included only common HLA alleles (MAF > 0.01) that were imputed with high certainty (R 2 > 0.8).
Statistical framework for HLA association testing. Given that the frequency of HLA alleles can vary by ethnicity, we performed HLA association testing Europeans and East Asians separately. We used logistic regression models to test the additive effects of HLA allele dosages with adjustment for significant PCs of ancestry. For multi-allelic loci, we used the following logistic regression model: where m indicates a total number of alleles at a specific multi-allelic locus, j indicates a specific allele being tested, and x j,i is the imputed dosage for allele j for individual i; β 0 represents the intercept and β j represents the additive effect of an allele j; P k,i denotes the value for kth PC of individual i, n is the total number of significant PCs in the dataset, β k is the effect size of principal component k. We compared log-likelihoods of two nested models: the full model containing the test locus and relevant covariates with the reduced model without the test locus, but with the same set of covariates. The deviance was defined as −2×log likelihood ratio, which follows a X 2 -distribution with m − 1 degrees of freedom, from which we calculated P-values. In addition to multi-allelic tests, we also performed biallelic tests of association for all individual SNPs, classical HLA alleles, and individual amino acid residues in HLA molecules. All analyses were performed using dosage method under additive coding. Stepwise conditioning analyses across the HLA region were performed using both multi-allelic and bi-allelic coding of HLA variants. In each round of stepwise conditioning, we first included the most significant variant as the covariate in the logistic regression model. If additional independently associated markers are detected, they are included as covariates in our subsequent models. We repeated these analyses until no residual associations across the entire locus were observed.
Analysis of polymorphic amino-acid sites in HLA genes. To test the effects of individual amino acid substitution sites within the HLA-DRB1 and HLA-DQA1 genes, we applied a conditional haplotype analysis using fully phased haplotypes across the HLA region. We tested each single amino acid position by first identifying the m possible amino acid residues occurring at that position and then using m − 1 degrees of freedom test to derive P-values, with a single amino acid residue arbitrarily selected as a reference. For conditioning on individual amino acid sites, we used the following procedure: by adding a new amino acid position to the model, a total of k additional unique haplotypes were generated and tested over the null model (without a new amino acid site) using the likelihood ratio test with k degrees of freedom. If the new position was independently significant, we further updated the null model to include all unique haplotypes created by all amino acid residues at both positons to identify another independent position. The procedure was repeated until no statistically significant positions were observed.
Testing for pairwise epistasis. Multiplicative interactive effects were tested using logistic regression model; SNPs were coded under additive genotype coding (0, 1, 2), and interaction terms were defined as simple products of genotypes. To screen for interactions, we tested a lead SNP at each of the three non-HLA loci against each of the five independent HLA SNPs (three in Europeans and two in East Asians) resulting in a total of 15 independent tests. We additionally tested for all pairwise interactions between the three non-HLA loci shared between both ethnicities resulting in three additional tests. In each case, we used a likelihood ratio test comparing two nested models: the full model with both main effects and the interaction term to the reduced model with main effects only. Given a total of 18 pairwise interaction tests, we used a Bonferroni-corrected significance threshold of 0.05/18 = 2.8 × 10 −3 . In secondary analyses, we explored if significant HLA risk haplotypes interact with the PLA2R1 risk allele in the six cohorts with fully imputed classical HLA alleles. This included a total of 2759 East Asians (803 cases and 1956 controls) and 4507 Europeans (1880 cases and 2627 controls).
Functional annotations of GWAS loci. We used several different approaches to perform functional annotations of our significant loci. We first defined the region of each locus as +/−400 kb of the index SNP. Using ANNOVAR software 51 , we identified functional variants within each region that were in strong linkage disequilibrium (r 2 > 0.8) with the top SNP, including all known coding, splicing, 3′UTR and 5′UTR variants (Supplementary Table 12). To assess for potential functional variation in non-coding regions, we used our recently proposed tissue-specific functional scoring method (FUN-LDA) 11  To test for eQTL effects in kidney tissue, we used gene expression data from the NEPTUNE study 13 . This dataset is comprised of whole genome DNA sequence data and genome-wide transcriptome data (Affymetrix 2.1 ST chips) performed on microdissected glomerular (N = 136) and tubulointerstitial (N = 166) tissue compartments from kidney biopsies of patients with nephrotic syndrome. For each locus, we tested the index SNP and its high LD SNPs (r 2 > 0.8). Testing for cis-eQTL effects involved all transcripts within 1-Mb region centred on each SNP using the additive linear regression with adjustments for age, sex, PEER factors and first 4 PCs of ancestry as described previously 13 . In addition to kidney tissue, all loci were similarly interrogated for eQTL effects in the GTEx database version 8. Because NFKB1 and IRF4 both encode transcription factors, we also explored their potential binding in close proximity of each other or PLA2R1 gene. Based on the Chip-seq data for all five subunits of NFκB complex in immortalized lymphocytes (GSE55105) 29 , we found that the top SNPs at PLA2R1 and IRF4 loci intersect potential NFκB complex binding site ( Supplementary Figs. 8 and  13). In addition, rs230492, a variant in strong LD (r2=0.94) with the top SNP at the NFKB1 locus intersects a potential IRF4 binding site based on the Chip-seq data of IRF4 (GEO: GSM803390).
Pleiotropy analysis. We used the latest GWAS catalogue data to perform systematic cross-annotation of our top risk alleles against all other published GWAS findings. We first identified all genome-wide significant SNPs (P < 5 × 10 −8 ) reported in the catalogue that resided within the genomic regions of association with MN. We assessed the extent of linkage disequilibrium (r 2 ) between these SNPs and the top MN risk alleles based on the combined European and East Asian sequence data from 1000 Genomes (phase 3). We next defined the directionality of pleiotropic effects as either concordant or opposed in relationship to the MN risk alleles. In addition, we queried each qualifying SNP from the catalogue against our genome-wide summary statistics to extract the odds ratios and P-values for associations with MN. We defined overlapping susceptibility alleles if r 2 exceeded 0.2 (Supplementary Table 13). Lastly, we constructed a susceptibility overlap map that connects each of the MN loci to the previously associated GWAS traits and highlights associations with SNPs in high LD with the top MN signals (Fig. 5a).
The map was visualized with Cytoscape v.3.6 software.
Podocyte gene annotations of GWAS loci. To search for potential novel podocyte antigens encoded by our suggestive loci, we interrogated each locus against a podocyte-specific gene list predicted with in silico "nanodissection" approach (Supplementary Table 14). This computational approach used Affymetrix gene expression data from micro-dissected glomerular and tubulo-interstitial compartments of 452 renal biopsies 52 . In addition, we cross-annotated our positional candidates against the list of native podocyte proteins discovered by proteomic profiling of mouse podocytes 53 . All suggestive loci were additionally tested for multiplicative interactions with classical alleles, and the top most significant interactions were summarized in Supplementary Table 14. Lastly, we annotated all positional candidate genes using manual PubMed literature searches to prioritize genes with a previously established role in podocyte biology.
GR analysis. To estimate the cumulative effect of independently significant GWAS loci, we used a GRS approach. We first identified SNPs with independent contributions to MN risk at each locus using stepwise conditional analyses. Next, we tested for multiplicative interactions among those independently associated SNPs. We then built a logistic regression model including all independent SNPs along with their significant interaction terms to derive mutually adjusted effect sizes (Supplementary Table 15). Because we observed ethnicity-specific signals at the HLA locus, we generated ethnicity-specific for East Asians and Europeans separately. The ethnicity-specific risk scores were defined as a weighted sum of independent risk alleles and their significant interaction terms, weighted by their mutually adjusted effect sizes. The GRS was standardized using a Z-score transformation based on the mean and standard deviation for the distribution of ethnically matched controls, so that the standardized GRS was reflective of the distance between the raw score and the control mean in units of standard deviation. The final formulation of the GRS is as follows: where GRS i = genetic risk score for individual i, d i = dosage of risk allele (from 0 to 2) for individual i, East Asian Control Mean GRS = 1.6804, East Asian Control SD GRS = 1.003 where GRS i = genetic risk score for individual i, d i = dosage of risk allele (from 0 to 2) for individual i, European GRS control mean = 1.5089, European GRS control SD = 0.8202. The percentage of the total variance in disease risk explained was estimated using Nagelkerke's pseudo R 2 from the logistic regression model with the standardized GRS as a predictor and case-control status as an outcome. The performances of the ethnicity-specific GRS were estimated by the AUROC. We performed detailed GRS cut-off analyses in the discovery cohorts by selecting cutoff points on the ROC curve that provide specificities in the range from 95 to 100%. For each cut-off point, we calculate sensitivity, specificity, positive likelihood ratio (LR+), and negative likelihood ratio (LR-). All GRS analyses were implemented in R version 3.3.2.
CRS formulation. The CRS was formulated as a weighted sum of the GRS and serum anti-PLA2R antibody levels in U/mL. The weight was determined using logistic regression model, with standardized ethnicity-specific GRS (formulated as described above) and natural log-transformed anti-PLA2R antibody levels as two predictors and case-control status as an outcome. This resulted in the following model: where GRS i indicates the Z-transformed ethnicity-specific GRS for individual i, β 0 indicates the intercept, β 1 indicates the effect size of the GRS estimated based on discovery cohorts; αPLA2R i is the serum level of PLA2R antibodies in U/mL; the constant 0.001 is added to enable log transformation of undetectable (zero) levels; β 2 represents the effect size for anti-PLA2R antibody positivity estimated based on discovery cohorts. The weight for the In(αPLA2R i + 0.001) term was then defined as follows: Consequently, the weights for antibody levels were calculated for East Asian and European separately, which resulted in the following Crude CRS formulation: In the final step, the Crude CRS was Z-transformed using the mean and standard deviation for ethnicity-matched healthy controls: where European Control Mean CRS = −1.4982, European Control SD = 1.4354, E. Asian Control Mean CRS = 0.3724, E. Asian Control SD = 2.7503. The CRS performance was estimated by the area under the receiver operating curve (AUROC). Similar to GRS, we explored different CRS cut-offs that maximized specificity in the range 95 to 100%. For each cut-off, we calculated sensitivity, specificity, positive likelihood ratio (LR+), and negative likelihood ratio (LR−). The integrated discrimination improvement (IDI) and net reclassification improvement (NRI) 54 were calculated comparing the CRS test to the serum PLA2R antibody test alone. All CRS analyses were implemented in R version 3.3.2 (CRAN).
European GWAS validation cohorts. For the purpose of GRS validation, we utilized three previously published external GWAS cohorts of European ancestry: the UK, French, and Dutch Validation Cohorts 6 . These cohorts were composed of biopsy-documented cases of primary MN and ethnicity-matched healthy controls, totalling 2,887 individuals (550 cases and 2337 controls, see Supplementary Methods for details). For all three cohorts, we obtained primary genotype data after quality control analysis as published previously 6 and performed phasing and imputation using Eagle v.2.3 and Minimac3 with European populations of Phase 3 1000 Genome Project as a reference. The cases and controls were imputed jointly. The GRS for each individual was determined using genotype dosages for the SNPs included in the score. The performance of GRS was then analyzed individually in each cohort, and in all cohorts combined (Supplementary Tables 17-19). PLA2R antibody testing. In total, we determined serum antibody levels in N = 2331 study participants with genetic data (N = 1488 cases, N = 300 healthy controls, and N = 543 disease controls) across all cohorts and ethnicities. The ancestrymatched diseased controls were recruited among patients commonly presenting with nephrotic syndrome, including FSGS, MCD, IgAN. For all MN cases and disease controls, serum samples were obtained near or at the time of kidney biopsy and any samples obtained more than six months after the biopsy were excluded from the analysis. All individuals that underwent antibody testing had matched genetic data, enabling derivation and diagnostic testing of the CRS.
We performed a standardized measurement of serum anti-PLA2R Ab levels using the anti-PLA2R ELISA (IgG) test kit (EUROIMMUN Medizinische Labordiagnostika AG), which employs the indirect ELISA methodology. The kit includes a 96-well microplate pre-coated with PLA2R, 5 calibrators (2, 20, 100, 500, and 1500 U/mL respectively), positive and negative control samples, peroxidaselabelled anti-human IgG (rabbit) enzyme conjugate, kit specific sample and wash buffers, Chromogen/substrate solution (TMB/H 2 O 2 ), and stop solution. The assay was run as per the protocol included with the kit. A 5-point calibrated analysis was used to calculate the results for each assay performed. A standard curve was generated based on the spectrophotometric reading of the five calibrators included on each microplate. As recommended by EUROIMMUN, the sample was called positive if antibody level was ≥20 U/mL.
In the discovery stage, we analyzed sera for a total of N = 459 East Asian and N = 1034 European participants. The European cohorts included 810 cases with MN, 99 healthy controls, and 125 disease controls (37 FSGS and 88 IgAN). The East Asian cohorts included 304 cases with MN, 56 healthy controls, and 99 disease controls (52 FSGS and 47 IgAN). The disease controls were not included in the GWAS discovery analysis, but were added to test for disease specificity of the serologic test. We note that one East Asian disease control with a clinical diagnosis of IgAN tested anti-PLA2R antibody positive at high titer; follow-up pathology review found evidence of previously unrecognized sub-epithelial deposits diagnostic of MN in addition to IgA-dominant mesangial deposits; this individual was subsequently removed from the analysis.
In the validation phase, we analyzed sera for a total of N = 540 individuals including European case-control cohorts (N = 248 cases and N = 145 healthy controls) and a total of 147 European NEPTUNE participants (N = 36 cases and N = 111 disease controls). In secondary analysis, we extended testing to additional 180 NEPTUNE participants of non-European ancestry (N = 28 cases and N = 152 controls). This included 103 NEPTUNE participants of African American ancestry (16 cases and 87 disease controls) and 77 participants of Hispanic/Latino ancestry (12 cases and 65 disease controls). These numbers are smaller compared to the GRS validation cohorts, since not all NEPTUNE participants had sera sampled within 6 months of biopsy.
Testing for phenotypic correlations of the GRS. Using extensive clinical information available for our discovery cohorts, we investigated correlations between GRS and clinical traits from the time of kidney biopsy including age at diagnosis, 24-h proteinuria (P24), estimated glomerular filtration rate (eGFR), serum albumin (Alb) and serum anti-PLA2R1 antibody level (Supplementary Table 16). The eGFR was estimated based on serum creatinine level using the Modification of Diet in Renal Disease (MDRD) equation 55 . The proteinuria was quantified using spot urine protein-to-creatinine ratios. The values of proteinuria and eGFR were normalized by natural log-transformation. The serum albumin and anti-PLA2R1 antibody levels also required natural log-transformation. For each quantitative trait, we built a linear regression model with GRS as a predictor and each corresponding trait as an outcome. For dichotomous traits (anti-PLA2R seropositivity or presence of nephrotic range proteinuria), we used logistic regression with a GRS as a predictor and each binary trait as an outcome. The association analysis for age at diagnosis (biopsy) was performed before and after adjustment for sex and ancestry.
The tests of proteinuria, eGFR, serum albumin and serum anti-PLA2R1 antibody levels were carried out before and after controlling for age, sex and ancestry. Statistical analyses were implemented in R version 3.3.2.
SNP-based heritability. We estimated the SNP-based heritability of MN in East Asians and Europeans using the GCTA-GREML algorithm 30,56 . For this analysis, we included three European cohorts (European discovery 1, European discovery 2 and Turkish discovery) and three East Asian cohorts (Chinese, Korean and Japanese discovery) that had primary genotype data available for joint heritability analysis. We first estimated pairwise genetic relationship matrix between all individuals using autosomal SNPs. With the GCTA software 57 , we next estimated the disease variance explained by all autosomal SNPs. We transformed the estimate assuming an underlying liability scale and disease prevalence of 0.001. We derived a standard error (SE) and 95% confidence interval for each estimate.
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
All genome-wide summary statistics, including those presented in Fig. 1, are freely available for download on our lab website: www.columbiamedicine.org/divisions/kiryluk/ resources.php. The calculations of genetic risk score (GRS) and combined risk score (CRS) are implemented in the form of an online risk calculator, which is also freely available on our lab website. The PAGE consortium control genotype data is available on dbGAP under accession number phs000356.v2.p1. Primary genotype data for the European-1 discovery cohort is available under dbGAP accession number phs001984.v1. p1. Our IRB determined that the use of this dataset is restricted to genetic studies of kidney disease. Because of consent restrictions and/or country-specific privacy laws, we are unable to share primary genotype data on dbGAP for other cohorts. All data and summary statistics are available from the corresponding authors upon reasonable request.