Genome-wide Association Study Identifies 2 New Loci Associated With Anti-NMDAR Encephalitis.

BACKGROUND AND OBJECTIVES: To investigate the genetic determinants of the most common type of antibody-mediated autoimmune encephalitis, anti-NMDA receptor (anti-NMDAR) encephalitis.
METHODS: We performed a genome-wide association study in 178 patients with anti-NMDAR encephalitis and 590 healthy controls, followed by a colocalization analysis to identify putatively causal genes.
RESULTS: We identified 2 independent risk loci harboring genome-wide significant variants (p < 5 × 10-8, OR ≥ 2.2), 1 on chromosome 15, harboring only the LRRK1 gene, and 1 on chromosome 11 centered on the ACP2 and NR1H3 genes in a larger region of high linkage disequilibrium. Colocalization signals with expression quantitative trait loci for different brain regions and immune cell types suggested ACP2, NR1H3, MADD, DDB2, and C11orf49 as putatively causal genes. The best candidate genes in each region are LRRK1, encoding leucine-rich repeat kinase 1, a protein involved in B-cell development, and NR1H3 liver X receptor alpha, a transcription factor whose activation inhibits inflammatory processes. DISCUSSION: This study provides evidence for relevant genetic determinants of antibody-mediated autoimmune encephalitides outside the human leukocyte antigen (HLA) region. The results suggest that future studies with larger sample sizes will successfully identify additional genetic determinants and contribute to the elucidation of the pathomechanism.

Antibody-mediated encephalitides are a recently discovered group of rare diseases caused by autoantibodies against CNS system antigens.[1] Subgroups are defined by the respective target antigens. The most common subgroup is caused by immunoglobulin G (IgG) class antibodies against the GluN1 subunit of the NMDA receptor (NMDAR), the most important excitatory neurotransmitter in the CNS. NMDAR antibodies cause internalization of surface NMDAR, thereby reducing signal transduction. Anti-NMDAR encephalitis affects children and adults with a female preponderance. The estimated prevalence is 0.6/100,000 population.[2] It manifests with behavioral changes, psychiatric symptoms, epileptic seizures, memory dysfunction, movement disorders, and loss of consciousness and often responds favorably to immunotherapy.[3] A definite diagnosis requires the detection of anti-NMDAR antibodies of the IgG class in serum and/or CSF. Known trigger factors include ovarian teratomas with ectopic expression of NMDAR and viral (mostly herpes simplex virus type 1 [HSV-1]) encephalitis.[1] In a first genome-wide association study (GWAS) of antibody-mediated encephalitides in a much smaller sample, we found no variants showing genome-wide significant association with anti-NMDAR encephalitis.[4] For this study, we doubled the sample size, modified the analysis parameters, and added colocalization analysis to identify putatively causal genes.

Methods

Study Population

In this case-control study, we analyzed 212 samples from patients with NMDAR antibodies (197 German and 15 Czech patients) collected in the years 2015–2020. Anti-NMDAR encephalitis was classified according to consensus criteria based on a compatible clinical syndrome together with detection of IgG-NMDAR antibodies in serum and/or CSF.[3] One hundred seven patients were already included in our previous GWAS.[4] The additional 105 individuals were either recruited via the German Network for Research on Autoimmune Encephalitis (GENERATE, n = 80) or specifically for this study (n = 25). For patient recruitment, we contacted the centers of the GENERATE network as well as other neurologic departments caring for patients with antibody-mediated encephalitides. All contributing scientists are listed in eTable 1, links.lww.com/NXI/A587. Healthy control samples (n = 1,219) were obtained from the PopGen study, a population-based biobank from northern Germany.[5,6]

Genotyping

Genomic DNA was isolated from blood (n = 150) or saliva (n = 62) using standard procedures. The samples were genotyped in 3 batches at the Institute of Clinical Molecular Biology, Kiel University on Infinium Global Screening Array (GSA; Illumina). Array version 2.0 was used for cases and version 1.0 for healthy controls. Genotypes were called using Illumina GenomeStudio 2.0 according to the manufacturer's instructions using in-house clusterfiles. We previously determined a > 99.8% genotype concordance for DNA isolated from blood and saliva genotyped on the GSA array in 8 individuals.

Quality Control and Imputation

We used PLINK v.1.9,[7] R v.3.6.3,[8] and the Illumina GenomeStudio for genotype quality control. First, we excluded all nonoverlapping variants between the 2 different versions of the GSA chip, variants with multicharacter allele codes, insertions, deletions, duplicated markers, and ambiguous A/T and G/C variants. We determined genotyping sex by the X-chromosome inbreeding coefficients, with F < 0.2 being female and F > 0.8 being male, and excluded samples with discordance between reported and imputed sex. After that, we filtered first variants and then individuals with a relaxed threshold for a call rate of less than 85%, followed by a stringent threshold of 98%. We applied a minor allele frequency (MAF) filter of 1%, as well as filters for significant deviation from Hardy-Weinberg equilibrium (HWE; p < 1 × 10−6) in controls, informative missingness (p < 1 × 10−5), and outlying heterozygosity rate (mean ± 3 SDs). To determine duplicated or cryptically related individuals, we used pairwise genome-wide estimates of the proportion of identity by descent (IBD) on a pruned data set containing only markers in low linkage disequilibrium (LD) regions (pairwise r2 < 0.2) and excluded those more closely related than third-degree relatives (IBD > 0.125). Of each identified sample pair, we kept the individual with a higher call rate. To identify ethnic outliers, we used a procedure similar to the one suggested in the R package plinkQC[9]: we combined the genotype data with the samples of the publicly available 1000 Genomes Project[10] and performed a principal component (PC) analysis on the merged data set. A European center was determined by the first 2 PCs of known European samples, and the Euclidean distance from this center determined the ethnical assignment with samples more than 1.5 times the maximal European Euclidean distance away from the center being excluded. The remaining individuals were used for preliminary association analysis based on which we visually inspected the cluster plots of all variants with a p value < 10–4 and discarded variants without adequate cluster separation. To overcome issues with population stratification, we matched controls by ancestry and sex to cases with the R package PCAmatchR,[11] leading to 590 control samples for the analysis and approximately 3 controls per case. An exact match on sex was used because there were significantly more female samples in the case samples than in the control samples.

Imputation was performed on the quality-assured data set, containing 768 individuals (590 controls and 178 cases) and 446,353 variants. Subsequently, 26,356,529 variants were imputed based on the Trans-Omics for Precision Medicine (TOPmed) r2 panel[12] using the TOPMed Imputation Server,[13] which uses (mini Markov-chain haplotyper 4) for imputation.[14] A quality check was performed, including variants with an MAF > 1%, an imputation quality score R2 > 0.7, and no significant deviation from HWE (p < 1 × 10−6) in controls, resulting in 8,073,349 variants.

Association Analysis

We conducted an association analysis on the whole data set using a genome-wide significance threshold of p < 5 × 10−8. We applied an additive logistic regression model, including sex and PCs, to estimate the association of each single nucleotide polymorphism (SNP) with the disease status. The number of PCs was chosen using scree plot analysis.[15] Population stratification was examined using the inflation factor λ and the visual inspection of quantile-quantile plots. To further distinguish between confounding factors like population stratification and polygenicity of the anti-NMDAR encephalitis trait, we performed LD score regression (LDSC) using the LDHub web interface.[16] Conditional analyses in which successively each genome-wide significant variant was included as a covariate were conducted to identify adjacent independent signals. We used 7,122 genotyped and quality controlled variants from the human major histocompatibility complex region on chromosome 6 to impute four-digit human leukocyte antigen (HLA) alleles using the R package HLA imputation using attribute bagging.[17] It uses attribute bagging to impute genotypes, and we chose a prediction model specifically for European ancestry and the Illumina GSA chip. We performed the association analysis with python HLA[18] using an additive, logistic model including sex and the first PC as covariates and adjusted p values with the Benjamini-Hochberg false discovery rate step-up method.

To examine the origin of the variant-trait association signals more closely, we analyzed the subsamples of patients with early/late disease onset (< or ≥25 years) and patients with/without a tumor. To prioritize genes putatively involved in the disease etiology, we investigated the overlap of expression quantitative trait loci (eQTL) from the Genotype-Tissue Expression (GTEx) project[19] as well as immune cell eQTL from the BLUEPRINT (A BLUEPRINT of haematopoietic epigenomes) project[20] and variants in the risk loci identified by this GWAS. We investigated whether these 2 independent signals might stem from the same causal variant by colocalization analysis conducted with coloc.[21] Coloc uses approximate Bayes factors to estimate posterior probabilities (PP) for common variants causal in the GWAS as well as the eQTL study. We studied all variants present in the GWAS results as well as in GTEx V7 for the 13 available brain tissues or present in the BLUEPRINT immune cell eQTL data within 100 kb up- and downstream of each gene in the 2 encephalitis risk loci. Coloc estimates PPs for 4 different scenarios. PP4 denotes the posterior probability that both traits—the disease association and the eQTL—are caused by the same variant. A PP4 over 70% was considered as evidence for colocalization. We used LocusZoom[22] and R to visualize the association results. All analyses and the presentation of the results in this article are based on the human genome version 38 (GRCh38/hg38).

Standard Protocol Approvals, Registrations, and Patient Consents

All participants gave written informed consent. Institutional review board approval was obtained from the ethical advisory boards of the Universities of Kiel and Luebeck (B337/13;13-162).

Data Availability

Summary level genetic data for all variants with p values < 1 × 10−4 are available from the corresponding author on reasonable request to any qualified investigator.

Results

Table 1 summarizes the clinical features of patients and control individuals demonstrating that patients included in our first GWAS[4] are comparable to the additional patients in this study regarding age, sex, and clinical features. However, the control individuals were much older than the patients. Genotype data for 212 individuals with anti-NMDAR encephalitis and 1,219 controls were available for analysis. After quality control procedures and control matching, 178 cases and 590 healthy controls remained (eTable 2, links.lww.com/NXI/A588). Imputation resulted in 8,073,349 quality-assured variants with an MAF > 1%. We incorporated sex and the first dimension of the PC analysis as indicated by scree plot analysis as covariates. In contrast to our first GWAS of anti-NMDAR encephalitis, we did not include age as a covariate (for rationale, see Discussion). The genomic inflation factor of λ = 1.03 indicated a low degree of population stratification (Figure 1A). The LDSC intercept was 1.01 (standard error = 0.01), with 44% of the genomic inflation attributable to confounding bias including population stratification and cryptic relatedness. This indicates that the majority of inflation is caused by polygenicity. We found 13 genetic variants in 2 distinct loci below the threshold of p = 5 × 10−8 for genome-wide significance (Figure 1, B-D, Table 2), with leading variants rs10902588 on chromosome 15 (OR = 2.24 [95% CI = 1.70–2.95], p = 1.78 × 10−8) and rs75393320 on chromosome 11 (OR = 2.20 [1.66–2.92], p = 3.78 × 10−8) as well as only 14 Kb further downstream rs11039155 with the same p value and OR. Conditional analysis including the top variants at each locus argues against the presence of any independent secondary signals (eFigure 1, links.lww.com/NXI/A586). The significant variants on chromosome 15 are located in the leucine rich repeat region 1 (LRRK2) gene. In the chromosome 11 locus, rs75393320 lies in the lysosomal acid phosphatase 2 (ACP2) gene, and rs11039155 is located in the nuclear receptor subfamily 1 group H member 3 (NR1H3) gene.

Table 1

Sample Characteristics

Figure 1

Association Plots for Anti-NMDAR Encephalitis

(A) Quantile-quantile plot of association analysis for 8,073,349 variants. The plot shows deviation from the null distribution in the upper tail, which corresponds to variants with the strongest evidence for association. (B) Manhattan plot of the association results. The plot shows −log10 marker-wise p values against their genomic base pair position. The red line indicates the genome-wide significance threshold of 5 × 10−8. (C) LocusZoom plot for the association between anti-NMDA receptor encephalitis and variants on chromosome 11 in the genomic region from 46.6 to 48.2 Mb. A circle represents a genotyped and a plus symbol an imputed variant. The r2 metric displays the pairwise LD between the leading and the respective variant. Gene positions are present in the bottom part. (D) LocusZoom plot for associations on chromosome 15 in the genomic region from 100.9 to 101.1 Mb.

Table 2

Identified Associations With a p Value < 5 × 10−8

Additional analyses of subpopulations defined by tumor status and age at onset yielded no genome-wide significant associations. In the previous GWAS, we observed a weak association between patients with anti-NMDAR encephalitis and the HLA-B*07:02 allele, preferentially in the patients with late disease onset. We did not confirm this association in the current GWAS. We were not able to analyze the recently reported association with HLA-DRB1*16:02 in a Chinese population[23] because the frequency of this allele is very low in the German population and no patient in our sample and only 1 control individual carried this allele. We did not detect any novel significant HLA associations. We performed an additional analysis, splitting the patient sample into the samples already included in our prior GWAS[4] and the newly acquired samples only (eTable 3, links.lww.com/NXI/A589). For all genome-wide significant variants in the complete sample, we found an identical direction and comparable magnitude of the ORs as well as p values between 1.37 × 10−7 and 1.13 × 10−3 demonstrating that both subsamples contributed to the final result. As readily apparent in Figure 1C, the significant variants on chromosome 11 are located in a gene-rich area with numerous further variants in high LD with the leading variant and p values less than 1 × 10−5. Coloc analysis showed colocalization with a PP4 > 0.7 between the sum of GWAS variants and GTEx eQTLs for the 3 genes NR1H3, ACP2, and mitogen-activated protein (MAP) kinase activating death domain (MADD) on chromosome 11 in brain tissues (Figure 2) and for the 4 genes NR1H3, ACP2, damage-specific DNA binding protein 2 (DDB2), and chromosome 11 open reading frame 49 (C11orf49) in various immune cells (Figure 3). We did not identify any single variant with a PP4 > 0.7. In contrast, we found no colocalizing eQTL signals for the GWAS signal on chromosome 15.

Figure 2

Colocalization Results for Brain Tissues

Gene- and SNP-wise results of the colocalization analysis for brain tissues represented in Genotype-Tissue Expression types. Only genes with a PP4 > 0.7 and variants with a p value < 10−5 are shown. ACP2 = acid phosphatase 2, lysosomal; GWAS = genome-wide association study; MADD = mitogen-activated protein kinase activating death domain; NR1H3 = nuclear receptor subfamily 1 group H member 3.

Figure 3

Colocalization Results for Immune Cells

Gene- and SNP-wise results of the colocalization for immune cells represented in the BLUEPRINT data set. Only genes with a PP4 > 0.7 and variants with a p value < 10−5 are shown. ACP2 = acid phosphatase 2, lysosomal; C11orf49 = chromosome 11 open reading frame 49; DDB2 = damage-specific DNA binding protein 2; GWAS = genome-wide association study; NR1H3 = nuclear receptor subfamily 1 group H member 3.

Discussion

Except for the HLA complex, the genetic determinants of antibody-mediated encephalitides are unknown. The collection of sufficiently large sample sizes for genetic analyses is hampered by the low disease prevalence of anti-NMDAR encephalitis, which is estimated to be around 0.6/100,000 population.[2] Despite the small sample size, we were able to find 2 distinct genomic regions on chromosomes 11 and 15 harboring genome-wide significant disease-associated variants. We did not identify additional loci containing variants with suggestive p values very close to genome-wide significance. The locus on chromosome 15 encompasses ∼100.000 bp and contains only 1 gene, LRRK1. Although genetic variants in the homolog leucine rich repeat region 2 (LRRK2) are the most common cause of autosomal dominant Parkinson disease, no neurologic diseases are currently linked to LRRK1. In mice, LRRK1 and LRRK2 complement each other at least partially in the nervous system because only deficiency of both proteins causes a neurodegenerative phenotype and both proteins regulate autophagy.[24] LRRK1 is expressed in B cells and monocytes, suggesting a role in the immune system.[25] Indeed, LRRK1-deficient mice show alterations of B-cell development, failure to produce IgG3 class antibodies in response to non–T-cell dependent antigens, and a proliferation and survival defect on B-cell receptor stimulation.[26] Yet, there is currently no known connection between LRRK1 and autoimmunity. However, it is intriguing to speculate that LRRK1-mediated control of non–T-cell dependent B-cell activation could be dysregulated in patients with anti-NMDAR encephalitis. Indeed, this could explain the observation of frequent nonmutated, germ-line encoded NMDAR antibodies in patients,[27] the childhood and early adult manifestation, and the coexistence of additional autoantibodies.[28]

The borders of the genomic region on chromosome 11 harboring the second association signal are less defined. The region is much larger, exceeding 1 Mb, and comprises multiple genes. To generate a hypothesis concerning putatively causal genes in this region, we used colocalization analysis between eQTL data from GTEx for different brain regions as well as immune cell eQTL from the BLUEPRINT project. In brain tissues, we found evidence for colocalization between the genes ACP2 and MADD with eQTL for cerebellum and NR1H3 with eQTL for the hypothalamus. Although it is well known that the hippocampus is a prime target of anti-NMDAR encephalitis, the ubiquitous expression of NMDA receptors containing the GluN1 subunit in the brain, the manifold symptoms of anti-NMDAR encephalitis, and pathologic studies suggest an involvement of most if not all brain regions.[29] Therefore, we think that the cerebellum and hypothalamus are valid target regions. In immune cells, we detected colocalization of the genes NR1H3, ACP2, DDB2, and C11orf49 with eQTL in various immune cells including T-lymphocytes. NR1H3 and ACP2 show evidence for colocalization in both brain and immune cells. Unfortunately, B-lymphocytes/plasma cells, the producers of antibodies, are not represented in the BLUEPRINT data. Of the genes identified in the colocalization analysis, NR1H3 encoding the liver X receptor alpha (LXRα) is the best functional candidate. LXRα is a transcription factor whose activation inhibits inflammatory processes.[30] In the CNS, LXRα agonists inhibit proinflammatory cytokine production by microglia and astrocytes.[31] Knockout of LXRα in brain endothelial cells led to blood-brain barrier dysfunction, inflammation, and increased transendothelial mononuclear cell migration.[32] ACP2 is a lysosomal acid phosphatase used in lysosomal protein degradation, MADD is an adaptor protein involved in transmitting tumor necrosis factor alpha (TNF-α)-induced apoptotic signals, DDB2 is involved in DNA repair, e.g., after ultraviolet light damage, and C11orf49 encodes a protein of unknown function.

In contrast to our first GWAS of antibody-mediated encephalitis, we identified 2 independent genome-wide significant associations in this study. There are 3 important differences between our previous GWAS and the current one. First, doubling of the sample size led to larger statistical power; second, we carefully removed population outliers; and third, we chose a different set of covariates in the logistic regression model. In contrast to the first GWAS, we included only sex and the first PC in the current analysis. Scree plot analysis suggested using only 1 PC, which might in part be due to the stringent exclusion of ethnic outliers and careful control matching. Another difference to the first GWAS is the exclusion of age as a covariate. Genetic variants are stable throughout life. For common late-onset diseases, significantly younger controls than patients warrant inclusion of age as a covariate because many controls will still develop the disease during their lifetime. However, in this study, the disease is rare, and the controls are significantly older than the patients. Including age as a covariate led to partial masking of the signals contributing to the effects in this GWAS. The major shortcoming of this study is its small sample size, which on the one hand limits the power to detect true variant-disease associations and on the other hand did not allow to include an independent replication sample, thereby increasing the likelihood of false positives. In our opinion, increasing the sample size by international cooperation as well as locus fine-mapping by sequencing and analysis of multiethnic samples will be key in future studies. The history of GWAS has shown that in most diseases, increasing sample size was more important than more detailed phenotyping. This might be particularly true in antibody-mediated encephalitis because the antigen specificity itself leads to a rather homogenous biologic disease entity compared with most other neurologic diseases, e.g., polygenic neurodegenerative disorders. In summary, we performed a GWAS of anti-NMDA receptor encephalitis and identified 2 independent genome-wide significant association signals. Both genomic regions contain putative functional candidate genes. In addition, eQTL for 5 genes show significant colocalization with the association signal on chromosome 11.