Assessing the role of rare genetic variants in drug-resistant, non-lesional focal epilepsy.

OBJECTIVE: Resistance to antiseizure medications (ASMs) is one of the major concerns in the treatment of epilepsy. Despite the increasing number of ASMs available, the proportion of individuals with drug-resistant epilepsy remains unchanged. In this study, we aimed to investigate the role of rare genetic variants in ASM resistance.
METHODS: We performed exome sequencing of 1,128 individuals with non-familial non-acquired focal epilepsy (NAFE) (762 non-responders, 366 responders) and were provided with 1,734 healthy controls. We undertook replication in a cohort of 350 individuals with NAFE (165 non-responders, 185 responders). We performed gene-based and gene-set-based kernel association tests to investigate potential enrichment of rare variants in relation to drug response status and to risk for NAFE.
RESULTS: We found no gene or gene set that reached genome-wide significance. Yet, we identified several prospective candidate genes - among them DEPDC5, which showed a potential association with resistance to ASMs. We found some evidence for an enrichment of truncating variants in dominant familial NAFE genes in our cohort of non-familial NAFE and in association with drug-resistant NAFE.
INTERPRETATION: Our study identifies potential candidate genes for ASM resistance. Our results corroborate the role of rare variants for non-familial NAFE and imply their involvement in drug-resistant epilepsy. Future large-scale genetic research studies are needed to substantiate these findings.

Introduction

Epilepsy is one of the most frequent neurological disorders, with a lifetime prevalence of approximately 7 in 1000 individuals. A genetic component for many types of epilepsy has been established for many years. , Focal epilepsies of unknown etiology, also known as non‐acquired focal epilepsies (NAFE), are characterized by focal seizures, focal epileptiform EEG findings, and the absence of epileptogenic lesions on magnetic resonance imaging (MRI) except hippocampal sclerosis. They account for 20–40% of all epilepsies and harbor a significant genetic component.

Previous studies have identified common single nucleotide variants (SNPs) as significant predictors of temporal lobe epilepsy with hippocampal sclerosis, and of NAFE in general, but the associated effect sizes are low. The conceptualization of NAFE as a syndrome with a polygenic component is further corroborated by a recent finding that polygenic risk scores allow to differentiate healthy individuals from individuals with NAFE. , Besides common variants, rare variants associated with developmental and epileptic encephalopathy (DEE) and ultra‐rare truncating and deleterious missense variants are enriched in NAFE. The latter study demonstrated an enrichment of variants in a group of 19 genes encoding all GABAA receptors and in a group of 43 dominantly inherited known epilepsy genes, though not reaching exome‐wide significance for any single gene.

Resistance to antiseizure medications (ASMs) presents one of the major challenges in the treatment of individuals with epilepsy. Individuals are considered drug‐resistant when at least two tolerated and appropriate ASMs fail to achieve ongoing seizure freedom. Despite more than 20 available ASMs, these individuals are unlikely to become seizure‐free with further ASM changes or polytherapy. Although multiple new ASMs have been licensed in recent years, the proportion of people with epilepsy who are drug‐resistant has not significantly decreased.

Pharmacogenetic markers to identify early individuals likely to have broad pharmacoresistance could prove useful to streamline the management of people with drug‐resistant epilepsy, for example by directing them to alternative treatment approaches such as epilepsy surgery. Few studies have addressed this issue, and those have shown no or only marginal association of genetic markers with response to specific ASMs or broad pharmacoresistance. , , , Various theories have been proposed to explain drug resistance in epilepsy. The drug transporter hypothesis purports that genetic variation of transporter genes could influence the pharmacokinetics of ASMs. The target hypothesis claims that genetic variants in genes that encode target proteins for ASMs could cause drug resistance. , Yet, a considerable portion of individuals with epilepsy is resistant to multiple or any ASMs regardless of the drugs’ target proteins or kinetics. This is addressed by the intrinsic severity hypothesis. Individuals with frequent and severe seizures are more likely to develop resistance to treatment. , Thus, the same genetic/biologic factors that give rise to interindividual differences of epilepsy severity, despite similar etiology, could also influence drug resistance. Epigenetic modification of gene expression, via DNA methylation or histone acetylation, presents another viable theory for drug resistance.

This study aimed to identify the role of rare genetic variants for ASM resistance and as a predictor of NAFE in a cohort of 1,128 individuals with NAFE (762 non‐responders, 366 responders) and 1,734 healthy controls.

Patients and Methods

Main cohort

The epilepsy cohort is derived from the EpiPGX Consortium and from the Canadian Epilepsy Network (CENet). 1,128 individuals with non‐familial NAFE (i.e. no 1st or 2nd degree relatives with reported epilepsy) underwent whole‐exome sequencing (WES), including 762 non‐responders (NR) (396 women) and 366 responders (R) (171 women). 975 individuals were of Non‐Finnish‐European descent, the remainder was of French‐Canadian origin. Recruitment sites and their respective sample contribution are listed in Table S1. All individuals gave written informed consent to participate. The study was approved by local institutional review boards at each recruitment site.

Epilepsy syndrome classification was based on the current guidelines of the International League against Epilepsy (ILAE). Only individuals that fulfilled criteria for NAFE were included: history of focal seizures, non‐lesional MRI with the exception of hippocampal sclerosis, and EEG findings compatible with focal epilepsy. Individuals with reported moderate to severe intellectual impairment were excluded to avoid overlap with DEE.

We classified individuals as drug‐responsive if they achieved 12 consecutive months of seizure remission to the first tolerated and appropriate ASM in monotherapy, starting within two years of the institution of treatment. Cases with known relapse after the initial 12‐months remission could be included. We classified individuals as non‐responders if they experienced recurring seizures at a frequency of ≥4/ year for 12 months prior to the latest recorded visit, despite adequate trials of at least two appropriate and tolerated ASM trials; individuals that met this definition of non‐response in the past but achieved seizure control owing to surgery or alternative treatments (e.g. vagus nerve stimulator) were included.

Control cohort

For the case–control study, we were granted access to bam files of 10 UK10K WES datasets from the European Genome‐Phenome Archive (EGA) including 1,734 healthy individuals: EGAD00001000417, EGAD00001000418, EGAD00001000419, EGAD00001000420, EGAD00001000429, EGAD00001000431, EGAD00001000433, EGAD00001000438, EGAD00001000440 and EGAD00001000442.

Replication cohort

We were kindly provided by the Epi4k group (https://www.epi4k.org/) with a replication cohort of 350 individuals with WES data, 165 NR (101 women) and 185 R (89 women). All individuals were of European descent. Phenotype definitions were equivalent to the main cohort.

Bioinformatics

Genomic data from individuals with epilepsy was generated by the Canadian Epilepsy Network (CENet). Sequencing of whole exomes was performed at Genome Quebec Innovation Center (http://gqinnovationcenter.com/index.aspx?l=e). Genomic DNA was quantified using the Quant‐iT™ PicoGreen® dsDNA Assay Kit (Life Technologies). Libraries were generated robotically on a Sciclone (PerkinElmer) using the KAPA HTP Library Preparation Kit Illumina® platforms (Kapa Biosystems) as per the manufacturer’s recommendations. TruSeq adapters and PCR primers were purchased from IDT. Libraries were quantified using the Kapa Illumina GA with Revised Primers‐SYBR Fast Universal kit (Kapa Biosystems). The average size fragment was determined using a LabChip GX (PerkinElmer) instrument. Two hundred and fifty ng of 4 libraries were pooled together (total of 1000 ng per capture) prior to proceeding with the enrichment of the targeted regions using the Roche Nimblegen EZ Choice custom baits. Captures were performed robotically according to the manufacturer’s recommendations. Final libraries were quantified using the Quant‐iT™ PicoGreen® dsDNA Assay Kit (Life Technologies) and the Kapa Illumina GA with Revised Primers‐SYBR Fast Universal kit (Kapa Biosystems). The average size fragment was determined using a LabChip GX (PerkinElmer) instrument. For cases only analysis, the Illumina control software was HCS 2.2.58, the real‐time analysis program was RTA v. 1.18.64. Program bcl2fastq v1.8.4 was used to demultiplex samples and generate fastq reads. The filtered reads were aligned to reference Homo_sapiens assembly b37. Each readset was aligned to create a Binary Alignment Map file (.bam) and then a gvcf using the MUGQIC pipeline for DNAseq (https://bitbucket.org/mugqic/mugqic_pipelines#markdown‐header‐dna‐seq‐pipeline). For cases and controls we performed a coverage analysis for each sample to eliminate 1) individuals with less than 85% of sites with coverage between 10 and 300, 2) sites with less than 90% of the samples with coverage between 10 and 300 leaving 29. 3Mb for our case cohort alone and 23.5Mb for both cohorts merged. Then we performed joint calling of gvcfs that were merged into a single vcf using GATK version 3.7‐0 (https://software.broadinstitute.org/gatk/). The vcf was recalibrated, filtered, and annotated following the GATK best practice guideline. VEP software version 84 (https://useast.ensembl.org/info/docs/tools/vep/index.html) was used for variant effect prediction. Joint calling was performed for 1) cases only, 2) cases and controls 3) the replication cohort only, and 4) the replication cohort and controls. The further filtering steps included: selection of biallelic sites present on the consensus coding sequence (CCDS), exclusion of indel variants, selection of sites with a genotyping rate of at least 98% overall samples, and with a Hardy‐Weinberg equilibrium (HWE) of greater than 0.001 using Plink version 1.9 (https://www.cog‐genomics.org/plink2).

After joint calling of both our cases and controls, and of our replication cases and controls, we performed additional cleaning steps to minimize batch effects. We performed a logistic regression of base quality (Q) as the dependent variable and the genotype as the independent variable to identify variants that were associated with a low Q. We set a p‐value threshold of 0.01. Base quality was determined from bam‐files using samtools (http://samtools.sourceforge.net).

Analyses

We stored all data and performed all analyses on Compute Canada’s systems (https://www.computecanada.ca/). We used R version 3.4.0 (https://cran.r‐project.org/) to create all plots. We performed PCA using smartpca from Eigensoft package version 3 with SNPs at 0.01 frequency or more. SNPs were pruned using plink version 1.9 (‐‐indep‐pairwise 50 5 0.2).

Gene‐based association tests

We performed SNP set kernel association tests using the SKAT‐O function from the SKAT R package (https://cran.r‐project.org/web/packages/SKAT/index.html). The first 10 principal components were used as covariates to account for potential differences in populations structure as well as sex. We used Annovar to annotate variants for effect (synonymous, nonsynonymous, or truncating) and frequency. We defined genes as SNP sets. Default beta weights (1,25) were used to put more weight on rare SNPs. Bonferroni correction was applied for single gene testing for a given significance level of p = 0.05. The number of genes and the resulting Bonferroni‐corrected p‐value thresholds are shown in Table S2. We defined ultra‐rare variants (URVs) as MAF ≤ 0.001 in gnomAD. We performed association tests for all variants and URVs independently of variant effect and separately for nonsynonymous and truncating variants.

Study power

We calculated the necessary sample size to achieve 80% power using the SKAT package. For URVs, we assumed a MAF ≤ 0.001, for all variants a MAF ≤0.4, given a prevalence of NAFE of 0.01, a prevalence of non‐responders versus responders of 0.3, and alpha levels according to the respective p‐value thresholds. Necessary sample sizes are shown in Table S2.

Gene set‐based association tests

We compiled 4 gene sets: ADME (absorption, distribution, metabolism, excretion) genes, ASM target genes, epilepsy genes, and NAFE genes). They were based on proposed mechanisms of drug resistance or association with epilepsy (Table S3) and in analogy to previous studies. , We analyzed the gene sets in a similar approach to the gene‐based association tests using the SKAT‐O function. Since the gene sets were not entirely independent, we chose a false discovery rate (FDR) correction to account for multiple testing. A significant enrichment was defined at an FDR < 0.05.

Results

Cohort description

In total 1,128 individuals with non‐familial NAFE (762 non‐responders, 366 responders) and 1,734 healthy controls satisfied our inclusion criteria in the main cohort. Our replication cohort comprised 350 individuals with NAFE (165 non‐responders, 185 responders). PCA showed that population structure was similar in the main cohort and control cohort (Fig. S1A), as well as in the replication cohort and control cohort (not shown).

Assessing enrichment of SNPs in non‐responders with responders

After quality control and filtering, 377,416 variants remained in the analysis. In order to determine if rare genetic variants were predictors of resistance to ASMs, we performed gene‐based enrichment analyses in responders versus controls for all variants in the dataset (Fig. 1A), for ultra‐rare variants (Fig. 1B), and for nonsynonymous and truncating variants (Fig. S2). No gene was genome‐wide significantly associated with drug resistance after adjusting the p‐value threshold by Bonferroni correction. The most strongly associated genes are depicted in Table 1.

FIGURE 1

(A) SKAT‐O Manhattan plots of non‐responder epilepsy cases and responder epilepsy cases using WES variants with all types of effects from 18,248 genes and (B) using only URVs (nonsynonymous and truncating) from 16,580 genes. Red line represents the 0.05 significance threshold after Bonferroni correction on the number of genes.

Table 1

Overview of genes with strongest association in gene‐based analyses.

Variant types	Reporting threshold (P‐value)	Non‐Responders/ Responders (P‐value)	Responders/Controls (p‐value)	Non‐Responders/ Controls (p‐value)	All NAFE/ Controls (p‐value)
All	<10‐4	SIRT4 (2.7x10‐5)	RAB40AL (1.5 × 10‐5) SPOUT1 (2.8 × 10‐5) SERAC1 (5.9 × 10‐5)	ZKSCAN4 (2.8 × 10‐5) OR4Q3 (8.7 × 10‐5)	ZKSCAN4 (1.8 × 10‐5) DUPD1 (5.1 × 10‐5) CCND3 (9 × 10‐5)
All nonsynonymous	<10‐4	FAM46D (4.2 × 10‐5) DSCR3 (7.5 × 10‐5)	LRRTM3 (1.1 × 10‐5) RAB40AL (4.1 × 10‐5) RCVRN (8.7 × 10‐5) SERAC1 (8.9 × 10‐5) GRIN2B (8.9 × 10‐5)	‐‐	DUPD1 (6.4 × 10‐5)
All truncating	<10‐3	IFNA5 (4.2 × 10‐4) SDCBP2 (8.9 × 10‐4)	‐‐	‐‐
URVs	<10‐4	‐‐	LRRTM3 (9.9 × 10‐6) IFNW1 (5.3 × 10‐5) CFAP45 (9.8 × 10‐5)	UHRF1BP1 (5.6 × 10‐5) OR12D2 (6.3 × 10‐5)	‐‐
URVs nonsynonymous	<10‐4	‐‐	LRRTM3 (9.9 × 10‐6) IFNW1 (5.3 × 10‐5) CFAP45 (9.8 × 10‐5)	OR12D2 (6.3 × 10‐5) UHRF1BP1 (7.6 × 10‐5)	‐‐
URVs truncating	<5 × 10‐3	‐‐	‐‐	DEPDC5 (3.8 × 10‐3)	‐‐

Genes with the strongest association in the gene‐based SKAT‐O analyses for the four comparison groups and the six variant types. No gene reached genome‐wide significance. Reporting P‐value‐threshold has been adapted to take into account the number of variants included in the respective analyses.

URV nonsynonymous = ultra‐rare variants (MAF ≤ 0.001 in gnomAD, nonsynonymous), URV truncating (minor allele frequency ≤ 0.001 in gnomAD, ultra‐rare truncating variants).

(1.5 × 10

SPOUT1 (2.8 × 10

SERAC1 (5.9 × 10

(2.8 × 10

OR4Q3 (8.7 × 10

(1.8 × 10

DUPD1 (5.1 × 10

CCND3 (9 × 10

(1.1 × 10

(4.1 × 10

RCVRN (8.7 × 10

SERAC1 (8.9 × 10

GRIN2B (8.9 × 10

DUPD1 (6.4 × 10

IFNA5 (4.2 × 10

SDCBP2 (8.9 × 10

(9.9 × 10

IFNW1 (5.3 × 10

CFAP45 (9.8 × 10

UHRF1BP1 (5.6 × 10

OR12D2 (6.3 × 10

(9.9 × 10

IFNW1 (5.3 × 10

CFAP45 (9.8 × 10

OR12D2 (6.3 × 10

UHRF1BP1 (7.6 × 10

DEPDC5 (3.8 × 10

To determine whether rare variant enrichment in groups of functionally related candidate genes could predict resistance to ASMs, we performed gene‐set‐based enrichment analysis in four gene groups for ultra‐rare truncating and missense variants (Table 2). We found no significant enrichment for any gene group.

Table 2

Result of gene set analyses.

Variant type	Gene Sets (number of genes)
	ADME (406)	Target (76)	Epilepsy (80)	NAFE (20)
Non‐Responders vs Responders
URV nonsynonymous	0.91 (0.69)	0.91 (0.59)	0.91 (0.55)	0.55 (0.06)
URV truncating	0.91 (0.36)	0.91 (0.93)	0.93 (0.79)	0.91 (0.47)
Non‐Responders vs Controls
URV nonsynonymous	0.23 (0.05)	0.31 (0.12)	0.76 (0.73)	0.45 (0.34)
URV truncating	0.46 (0.32)	0.75 (0.75)	0.46 (0.26)	0.03 (0.004)
Responders vs Controls
URV nonsynonymous	0.25 (0.05)	0.60 (0.51)	0.60 (0.45)	0.55 (0.34)
URV truncating	0.54 (0.30)	0.81 (0.81)	0.30 (0.11)	0.25 (0.06)
NAFE cases vs Controls
URV nonsynonymous			0.54 (0.54)	0.40 (0.30)
URV truncating			0.31 (0.15)	0.05 (0.01)

Gene set‐based SKAT‐O results of four different gene sets and two SNP sets with different functional effects. The table shows the FDR‐adjusted p‐values and the raw P‐values in parentheses. After correction for multiple testing, the NAFE set showed a borderline significant enrichment of truncating variants in individuals with epilepsy versus controls, and in non‐responders versus controls. Significant results are depicted in bold.

URV nonsynonymous = ultra‐rare variants (MAF ≤ 0.001 in gnomAD, nonsynonymous), URV truncating (minor allele frequency ≤ 0.001 in gnomAD, ultra‐rare truncating variants).

Assessing enrichment of SNPs in responders with controls

After quality control and filtering, 477,200 variants remained in the analysis. In order to determine if rare genetic variants could predict responders to ASMs versus healthy controls, we performed gene‐based enrichment analyses for all variants in the dataset (Fig. 2A), for ultra‐rare variants (Fig. 2B), and for nonsynonymous and truncating variants (Fig. S3). No gene showed a genome‐wide significant association with drug response after adjustment by Bonferroni correction. Among the most strongly associated, yet not significant genes (Table 1), we saw an enrichment of nonsynonymous variants in LRRTM3 and GRIN2B.

FIGURE 2

(A) SKAT‐O Manhattan plots of responder epilepsy cases and controls using WES variants with all types of effects from 17,934 genes and (B) using only URVs (nonsynonymous and truncating) from 16,800 genes. Red line represents the 0.05 significance threshold after Bonferroni correction on the number of genes.

To determine whether rare variant enrichment in groups of functionally related candidate genes could predict response to ASMs, we performed gene‐set‐based enrichment analysis in four gene groups for ultra‐rare truncating and missense variants (Table 2). We found no significant enrichment for any gene group.

Assessing enrichment of SNPs in non‐responders with controls

After quality control and filtering, 477,200 variants remained in the analysis. In order to determine if rare genetic variants could predict non‐responders to ASMs versus healthy controls, we performed gene‐based enrichment analyses for all variants in the dataset (Fig. 3A), for ultra‐rare variants (Fig. S3), and for nonsynonymous and truncating variants respectively (Fig. 3B, Fig. S3). No gene showed a genome‐wide significant association with drug response after adjustment by Bonferroni correction. Among the most strongly associated, yet not significant genes (Table 1), we saw an enrichment of ultra‐rare truncating variants in DEPDC5.

FIGURE 3

(A) SKAT‐O Manhattan plots of non‐responder epilepsy cases and controls using WES variants with all types of effects from 17,934 genes and (B) using only truncating URVs from 2,656 genes. Red line represents the 0.05 significance threshold after Bonferroni correction on the number of genes.

To determine whether rare variant enrichment in groups of functionally related candidate genes could predict response to ASMs, we performed gene‐set‐based enrichment analysis in four gene groups for ultra‐rare truncating and missense variants (Table 2). We found a significant enrichment of truncating variants in the NAFE gene group in association with drug‐resistant epilepsy.

Assessing enrichment of SNPs in NAFE cases with controls

After quality control and filtering, 477,200 variants remained in the analysis. In order to determine if rare genetic variants could predict NAFE versus healthy controls, we performed gene‐based enrichment analyses for all variants in the dataset (Fig. 4A), for ultra‐rare variants (Fig. S3), and for nonsynonymous and truncating variants respectively (Fig. 4B, Fig. S5). No gene showed a genome‐wide significant association with drug response after adjustment by Bonferroni correction. Among the most strongly associated, yet not significant genes, we also found an enrichment of ultra‐rare truncating variants in DEPDC5.

FIGURE 4

(A) SKAT‐O Manhattan plots of all epilepsy cases and controls using WES variants with all types of effects from 17,934 genes and (B) using only truncating URVs from 2,656 genes. Red line represents the 0.05 significance threshold after Bonferroni correction on the number of genes.

To determine whether rare variant enrichment in groups of functionally related candidate genes could predict the risk of NAFE, we performed gene‐set‐based enrichment analysis in two gene groups for ultra‐rare truncating and missense variants (Table 2). We found a marginally significant enrichment of truncating variants in the set of NAFE‐genes.

Replication analysis

To test whether we could reproduce the sub‐threshold associations observed in the main cohort, we replicated the previous analysis steps in our replication cohort. After quality control and filtering, 328,145 and 327,901 variants remained in the responder/control and non‐responder/control analysis respectively. 1,109,232 variants remained in the non‐responder/responder analysis.

None of the analyses yielded genome‐wide significant results. The sub‐threshold associations of the aforementioned genes could not be reproduced (Fig. S6, Table S4). Among the strongest sub‐threshold signals were the gene PABPC3 and NPEPPS.

Discussion

In this exome‐based study of individuals with NAFE, we aimed to identify rare genetic variants as potential risk factor for drug‐resistant epilepsy. We also strived to corroborate the role of rare genetic variant for the risk of NAFE. We formed subgroups for responders and non‐responders and analyzed them against each other and against controls. No single gene reached exome‐wide significance, but we identified some potential candidate genes. However, we found an enrichment of rare truncating variants in known NAFE‐genes in association with drug‐resistance and as a predictor of NAFE.

Our finding of an enrichment of rare variants in NAFE replicates previous results from large‐scale sequencing studies that showed an enrichment of truncating variants in NAFE. In analogy to genetic generalized epilepsy, these findings corroborate previous studies that showed a polygenetic background for NAFE, however to a lesser degree , than for GGE. This could either imply that non‐genetic factors are more influential for the etiology of NAFE, but could also mean that the NAFE group is more heterogenous. Possibly, the presence of individuals in our cohort, whose epilepsy is not genetically determined at all, but due to undetected acquired inflammatory or structural changes, could have diluted a more robust effect. Moreover, unlike NAFE, the definition for GGE is cut more clearly and thus facilitates the assembly of more homogenous cohorts. Interestingly, we also found that enrichment of rare truncating variants in NAFE genes was associated with drug‐resistance, implying that the presence of rare variants promotes a more severe phenotype. This finding is in accordance with previous studies that showed the association of rare variants with drug‐resistance to specific ASMs: levetiracetam and valproic acid. The enrichment of rare variants in NAFE genes could not be shown in a direct comparison of non‐responders and responders. This was probably due to a lack of sufficient power owed to the much smaller sample size.

On the level of single genes, none proved to be exome‐wide significant. Given the limited power of the analyses this study was not able to detect exome‐wide significant loci. Yet, among the strongest associations for drug‐resistance was the gene DEPDC5, which was also part of the NAFE gene set and suggestively one of the main drivers of the aforementioned association. Variants in DEPDC5 have been identified in various MRI‐negative familial forms of NAFE, , , but also in individuals with cortical malformations. , Recent large‐scale sequencing studies identified DEPDC5 among the genes with the strongest association with familial and non‐familial NAFE. , It is therefore not surprising that DEPDC5 showed one of the strongest signals in our analysis of all individuals with epilepsy versus controls. However, so far it has not been analyzed whether DEPDC5 was associated with any ASM response profile. We found that DEPDC5‐variants were only enriched in non‐responders, but not responders. Potentially, DEPDC5 variant carriers feature less responsive forms of epilepsy. It is not known that DEPDC5 directly affects target proteins or kinetic pathways of current ASMs. However, DEPDC5‐associated ASM resistance could be related to potential subtle cortical malformations that evade detection by standard MRI. It has been well established that cortical malformations are associated with drug resistant seizures, and that overall patients with cortical malformations fare worse than those with other epilepsy‐associated entities after epilepsy surgery. Notwithstanding the epilepsy syndrome, DEPDC5‐related epilepsies have a rate of >50% of drug‐resistant cases and about only 10% responder rate to the first ASM. On the other hand, individuals with DEPDC5‐relateted malformations and focal epilepsy show a favorable outcome after epilepsy surgery. Thus, the identification of DEPDC5 variations could be a promising predictor for drug resistance in NAFE and could be useful to fast‐track ASM resistant individuals for surgery evaluation.

In the responder‐control analysis, we found a non‐significant enrichment of nonsynonymous variants in LRRTM3 in ASM responders. LRRTM3 is a regulator of excitatory synapse development. LRRTM3 regulates excitatory synapse density and also controls AMPA receptor surface expression in the dentate gyrus, one of the pivotal regions of epileptogenicity. Although there is no established link of LRRTM3 with epilepsy, variants in LRRTM3 could lead to reduced synaptic excitability. The likelihood of response to ASMs could thus be increased. Another finding in this analysis was a non‐significant enrichment of variants in GRIN2B, which encodes the beta‐2‐subunit (NR2B) of the glutamate‐activated N‐methyl‐D‐aspartate (NMDA) receptor. Gain‐of‐function GRIN2B variants have been described as the cause of a form of DEE, , probably as a result of increased neuronal hyperexcitability. Possibly loss‐of‐function variants could promote the opposite effect – a decrease in neuronal excitability, favorable to ASM response.

We strived to reproduce our results in a second cohort, even though we did not find any exome‐wide significant genes with our main cohort. We identified several loci in a p‐value range comparable with our main analysis. Yet, the top hits did not match the findings of the main cohort. One explanation could be differences in the population structure of the main and the replication cohort. About 15% of our main cohort were of French‐Canadian origin – a population known to harbor specific genetic characteristics. Nonetheless, we would not expect this to profusely affect the burden of URVs. It is far more likely that the discrepancies were related to the very limited power for gene‐based analyses. Our power calculations show that especially for URVs very large sample size are needed to achieve exome‐wide significance. The very limited sample size of the replication cohort makes it even improbable to generate the same results. To exemplify this, we could consider the analysis of URVs (nonsynonymous and truncating) for NAFE versus controls. This analysis, after all filtering steps, contained 161,701 variants in 18,697 genes, i.e. ~8.6 variants per gene. Assuming an even distribution of variants and a share of NAFE of ~17% this amounts to 1.5 variants per gene in the 350 NAFE individuals. Thus, the odds to replicate the observed sub‐threshold associations of our main cohort seem small.

It is likely that additional factors are involved in drug‐resistant epilepsy. There is some evidence that rare variants are associated with drug resistance to specific ASMs. These effects are likely to remain undetected in a cohort with broad drug resistance. By design, this study did not assess the role of intergenic and non‐exonic regions. For instance, enhancer regions that can be found thousands of base pairs away upstream or downstream of the gene could present interesting targets for future research. Somatic mutations in MRI‐negative, subtle cortical malformations also have to be considered. Epigenetic mechanisms such as DNA methylation and histone deacetylation are another possible factor in drug‐resistant epilepsy and could be a relevant factor in cortical malformations.

To corroborate the role of rare genetic variants in drug‐resistant epilepsy larger, preferably genome sequenced cohorts of patients will be necessary. Nowadays, the limiting factor is not the sequencing, but the deep phenotyping of large cohorts that require expertise and manpower. Thus, our study points out potential candidate genes, whose role has to be substantiated by future studies.

Conflicts of Interest

The authors declare that they have no conflicts of interest related to the present study and publication.

Author Contributions

SW, CM, SG, and PC designed the study. SW, CM, MM, and SG analyzed the data. SW, MK, SB, GLC, ND, CD, MRJ, BPCK, WSK, HL, AGM, TJO, SP, JWS, GJS, PS, FZ, SMS, PC, and the EpiPGX Consortium collected data. RK, SG and PC provided the computational infrastructure. SG and PC supervised the study. SW, CM, SG, GPC, HL, SMS, SG, and PC interpreted the analyses. SW, CM, and SG wrote the manuscript. All authors interpreted the data and revised the manuscript.

Figure S1. Principal component analysis. Depiction of the first and second principal component for non‐responders vs responders (A) and all epilepsy cases vs controls (B). Colouring discriminates by responder status in A, and case–control status in B.

Figure S2. SKAT‐O Manhattan plots for non‐responders with responders for different variant groups: (A): Truncating URV. (B): Nonsynonymous URV. (C): All truncating variants. (D): All nonsynonymous variants.

Figure S3. SKAT‐O Manhattan plots for responders with controls for different variant groups: (A): Truncating URV. (B): Nonsynonymous URV. (C): All truncating variants. (D): All nonsynonymous variants.

Figure S4. SKAT‐O Manhattan plots for non‐responders with controls for different variant groups: (A): all URV (truncating, nonsynonymous). (B): Nonsynonymous URV. (C): All truncating variants. (D): All nonsynonymous variants.

Figure S5. SKAT‐O Manhattan plots for all epilepsy cases with controls for different variant groups: (A): all URV (truncating, nonsynonymous). (B): Nonsynonymous URV. (C): All truncating variants. (D): All nonsynonymous variants.

Figure S6. SKAT‐O Manhattan plots for replication analysis: (A, B): Non‐responders with controls for all variants and URVs. (C, D): Responders with controls for all variants and URVs. (E, F): Non‐responders with responders for all variants and URVs. (G, H): All epilepsy cases with controls for all variants and URVs.

Click here for additional data file.

Table S1. Recruiting site contributions. Abbreviations: EKUT=University Hospital Tübingen; IGG=Insituto Gaslini Genova; RCSI=Royal College of Surgeons in Ireland; UCL=University College London; UKB=University Hospital Bonn; ULB=Université Libre de Bruxelles; ULIV=University of Liverpool; UMCU=University Medical Centre Utrecht; UV=University Hospital Vienna; CHUM=Centre Hospitalier de l’Université de Montréal.

Table S2. Number of genes for each respective gene‐based analysis, Bonferroni corrected p‐value threshold for significant findings, given a significance level of p = 0.05, and estimated necessary sample size to achieve 80% power given the respective alpha (corrected p‐value). Sample size was calculated based on a prevalence of non‐responders among all epilepsy patients of 0.3, a prevalence of NAFE of 0.01, of NAFE responders of 0.07, and of NAFE non‐responders of 0.03. For URV analyses a MAF ≤0.001 was selected, for the analysis of all variant a MAF of ≤0.4.

Table S3. Gene set compositions: Genes included in the gene sets for the gene set‐based analyses.

Table S4. Overview of genes with strongest association in gene‐based replication analyses. Genes with the strongest association in the gene‐based replication SKAT‐O analyses for the four comparison groups and four variant types. No gene reached genome‐wide significance. Reporting P‐value‐threshold has been adapted to take into account the number of variants included in the respective analyses. URV nonsynonymous = ultra‐rare variants (MAF ≤ 0.001 in gnomAD, nonsynonymous), URV truncating (minor allele frequency ≤ 0.001 in gnomAD, ultra‐rare truncating variants).

Click here for additional data file.