Genome-wide association study of febrile seizures implicates fever response and neuronal excitability genes.

Febrile seizures represent the most common type of pathological brain activity in young children and are influenced by genetic, environmental and developmental factors. In a minority of cases, febrile seizures precede later development of epilepsy. We conducted a genome-wide association study of febrile seizures in 7635 cases and 83 966 controls identifying and replicating seven new loci, all with P < 5 × 10-10. Variants at two loci were functionally related to altered expression of the fever response genes PTGER3 and IL10, and four other loci harboured genes (BSN, ERC2, GABRG2, HERC1) influencing neuronal excitability by regulating neurotransmitter release and binding, vesicular transport or membrane trafficking at the synapse. Four previously reported loci (SCN1A, SCN2A, ANO3 and 12q21.33) were all confirmed. Collectively, the seven novel and four previously reported loci explained 2.8% of the variance in liability to febrile seizures, and the single nucleotide polymorphism heritability based on all common autosomal single nucleotide polymorphisms was 10.8%. GABRG2, SCN1A and SCN2A are well-established epilepsy genes and, overall, we found positive genetic correlations with epilepsies (rg = 0.39, P = 1.68 × 10-4). Further, we found that higher polygenic risk scores for febrile seizures were associated with epilepsy and with history of hospital admission for febrile seizures. Finally, we found that polygenic risk of febrile seizures was lower in febrile seizure patients with neuropsychiatric disease compared to febrile seizure patients in a general population sample. In conclusion, this largest genetic investigation of febrile seizures to date implicates central fever response genes as well as genes affecting neuronal excitability, including several known epilepsy genes. Further functional and genetic studies based on these findings will provide important insights into the complex pathophysiological processes of seizures with and without fever.

See Olson and Poduri (.

Introduction

Fever is a systemic response to infection characterized by elevated core body temperature, orchestrated by a complex interplay between the immune system and the neuronal circuitry.[1] Although generally conferring a survival benefit, the fever response may occasionally induce adverse events. In young children, fever lowers the seizure threshold and 2–5% of children of European ancestries experience febrile seizures before 5 years of age.[2] Most febrile seizures are benign and self-limiting with no recurrence or further sequelae.[3] However, febrile seizures may herald the onset of various epilepsy syndromes, including very severe ones, with a long-term study indicating 7% of children with febrile seizures subsequently develop epilepsy.[4] Seizure disorders are also in varying degree associated with neurological and psychiatric comorbidities.[5]

The mechanisms by which fever can result in neuronal hyperexcitability and seizures are still incompletely understood. Increased brain temperature is known to directly affect the sensitivity of ion channels, which might be sufficient to generate febrile seizures.[6,7] Other suggested mechanisms include alkalosis from fever-induced hyperventilation[8] and effects of inflammatory mediators such as cytokines on neuronal excitability.[9] Although the exact mechanisms underlying febrile seizures remain elusive, it is well recognized that the exposures that precipitate febrile seizures act together with developmental and genetic factors.[10]

Most genetic studies of febrile seizures have focused on rare epilepsy syndromes that include febrile seizures as part of the clinical presentation.[11-14] Less is known about the genetic background of common febrile seizures. We previously conducted a genome-wide association study (GWAS) of febrile seizures and identified four susceptibility loci for febrile seizures in general, implicating the sodium channel genes SCN1A and SCN2A, a TMEM16 transmembrane protein family gene (ANO3), and a region on 12q21.33 associated with magnesium levels.[15] We further identified two loci distinctly associated with febrile seizures as an adverse event following the measles, mumps and rubella (MMR) vaccination, implicating the innate immune system genes IFI44L and CD46.[15] A few other loci have also been reported in candidate gene studies of smaller febrile seizure cohorts, including variants in the synaptic zinc transporter gene ZNT3,[16] and the immunoregulatory cytokine gene IL6.[17] Despite this progress, the identified variants only explain a small fraction of the variance in disease liability.

The present study is based on a strict GWAS meta-analysis design, including 7635 febrile seizure cases and >80 000 controls. Here, our aims are to identify novel robustly associated genetic loci for febrile seizures, to investigate potential underlying genetic mechanisms, to analyse shared genetic susceptibility with epilepsies and neuropsychiatric diseases and to study the association between polygenic risk for febrile seizures and an individual’s history of hospital admission with a febrile seizure diagnosis.

Materials and methods

Participants

The discovery stage included two cohorts of Danish ancestry with existing GWAS data and phenotypic information from the Danish National Patient Register.[18] The register includes information from all inpatient admissions in Denmark since 1977 and all emergency and outpatient hospital contacts since 1995. Diagnoses are coded according to the International Classification of Diseases (ICD) using the 8th revision (ICD-8) from 1977 to 1993 and the 10th revision (ICD-10) from 1994 onwards. Eligible febrile seizure cases were defined as children registered with ICD-8 code 78021 or ICD-10 code R560 before the age of 7 years, and controls were required not to have any febrile seizures or epilepsy diagnosis codes (ICD-8 codes 345, 78020, 78029 and ICD-10 codes G40-G41 or R568) in the National Patient Register. The first discovery cohort included 1991 febrile seizure cases and 4097 controls from our previous GWAS of febrile seizures conducted at Statens Serum Institut (SSI).[15] This study had a focus on febrile seizures occurring as an adverse event following MMR vaccination, and 925 of the cases included in the current study had a febrile seizure episode that occurred in a risk window of 9 to 14 days after the date of their MMR vaccination. These are referred to as MMR-related febrile seizure cases. The remaining 1066 cases included in the current study experienced febrile seizures with no temporal relation to MMR vaccination.[15] The second discovery cohort was composed of 2511 febrile seizure cases and 46 952 controls from the Initiative for Integrative Psychiatric Research (iPSYCH) study.[19] The iPSYCH sample included patient groups of six neuropsychiatric disorders [autism spectrum disorder (ASD), attention deficit/hyperactivity disorder (ADHD), schizophrenia, affective disorder, bipolar disease and anorexia] and a random population sample.

The replication stage included three cohorts with genome-wide imputed data available or custom genotyping. First, from SSI, we included an additional 2004 cases (411 MMR-related febrile seizures and 1593 MMR-unrelated) and 3208 controls. Second, from the Danish Blood Donor Study (DBDS), we analysed 928 cases and 28 862 controls. Third, we included an Australian cohort of 201 febrile seizure cases and 847 controls. In the Danish replication cohorts, case and control definitions were the same as in the discovery stage. In the Australian cohort, a comprehensive assessment of each patient for diagnosis of febrile seizure was obtained using a validated seizure questionnaire, clinical evaluation by an experienced epileptologist and review of relevant medical records and clinical investigations as described previously.[16] Cases were unrelated patients of European ancestries recruited from Australia diagnosed with febrile seizures and controls were individuals without febrile seizures from the same population, a proportion of which were reported in a previous study.[16] The study design is schematically illustrated in Fig. 1 and Supplementary Fig. 1.

Figure 1

Overall study design.

Ethics statement

The SSI study was approved by the Scientific Ethics Committee for the Capital City Region (Copenhagen) and the Danish Data Protection Agency. The Scientific Ethics Committee also granted exemption from obtaining informed consent from participants (H-3–2010-003) as the study was based on biobank material. The Danish Scientific Ethics Committee, the Danish Data Protection Agency and the Danish Neonatal Screening Biobank Steering Committee approved the iPSYCH study. DBDS was approved by the Central Denmark Regional Committee on Health Research Ethics (M-20090237) and the Danish Data Protection Agency, Copenhagen (2012–58-0004, RH-30–0444/I-suite no.: 00922).

For the Australian replication cohort, all experimental protocols were approved by the Human Research Ethics Committee of Austin Health, Melbourne, Australia. All methods were carried out in accordance with the approved guidelines of the Human Research Ethics Committee of Austin Health, Melbourne, Australia. Informed consent was obtained from all human participants.

Genotyping, data cleaning and imputation

The SSI discovery cohort had existing GWAS data available from genotyping with Illumina Omni Bead Arrays and GenomeStudio software.[15] Data cleaning included filtering out individuals that (i) had >3% missing genotypes; (ii) had an autosomal heterozygosity rate deviating by >2.5 standard deviations (SD) from the mean; (iii) had discordant sex information; or (iv) were >6 SD away from the mean of any of the first five principal components in a principal component analysis. We also estimated pairwise identity by descent based on the GWAS data for all possible pairs of individuals. We then excluded one individual from each pair with an estimated identity by descent proportion >0.1875, corresponding to the point halfway between second- and third-degree relatives. Next, we excluded single nucleotide polymorphisms (SNPs) based on a missing rate of >2%, a minor allele frequency (MAF) of <0.01, deviations from Hardy–Weinberg equilibrium (P < 1 × 10−6), and genotyping batch differences in call rate or allele frequencies. The remaining 540 064 SNPs were used for imputation, which was done remotely at the Michigan Imputation Server (see ‘Web resources’ section) using the 1000 Genomes Project Phase 3 reference panel. After imputation, 8 012 441 autosomal variants with MAF ≥ 0.01 and imputation INFO > 0.8 were available for analysis in 1991 cases and 4097 controls.

The iPSYCH discovery cohort comprised 78 050 individuals with existing GWAS data available from genotyping with the Illumina Infinium PsychChip v.1.0.[19] Data cleaning involved filtering out SNPs that had a MAF of <0.01 or deviated from Hardy–Weinberg equilibrium (P < 1 × 10−6). Structural variants, tri-allelic SNPs, non-autosomal SNPs and SNPs that were strand ambiguous (A/T or C/G) or did not uniquely align to the genome were also excluded, resulting in an imputation backbone of 246 369 autosomal SNPs. A total of 364 samples failing basic quality control (discordant sex information, >1% missing genotypes, abnormal heterozygosity, duplicate sample discordance) were also excluded. Based on the good quality SNPs, all individuals were phased in a single batch using SHAPEIT3[20] and imputed in 10 batches using Impute2[21] with reference haplotypes from the 1000 genomes project phase 3. Before analysis, we filtered out participants who were: (i) of non-European ancestries (outliers in principal component analysis based on genotyped SNPs); (ii) related to another participant in the sample (corresponding to an identity by descent proportion >0.1875); (iii) were included in the SSI discovery or replication sample; or (iv) had an epilepsy or non-febrile seizure diagnosis. We also balanced the ratio of febrile seizure cases to controls, so that it was the same within each of the participant groups of the iPSYCH study as illustrated in Supplementary Fig. 1B. After these steps, 7 264 301 autosomal variants with MAF ≥ 0.01 and imputation INFO >0.8 were available for analysis in 2511 febrile seizure cases and 46 952 controls.

GWAS meta-analysis

In each discovery cohort, the data were analysed by logistic regression under an additive genetic model, including sex and four principal components as covariates. Variant positions were based on GRCh37/hg19 and alleles were labelled on the positive strand of the reference genome. A total of 6 799 883 variants were available in both discovery cohorts and were included for further analysis. We carried out combined analysis of the two discovery cohorts using fixed-effects inverse variance-weighted meta-analysis as implemented in METAL.[22] Heterogeneity between studies was assessed using the I2 statistic and Cochrane’s Q-test.[23] To assess possible confounding effects on the distribution of test statistics, we calculated the genomic inflation factor[24] and the linkage disequilibrium (LD) score regression intercept.[25] Loci with either one or more variants associated at P < 5 × 10−8 or two or more variants associated at 5 × 10−8 < P < 1 × 10−6 were identified and one top SNP from each locus was taken forward to the replication stage. To search for additional independent associated variants, we conditioned on the imputed allelic dosage of the top SNP at each locus in each discovery stage cohort and meta-analysed the results. Any SNPs reaching P < 1 × 10−6 in these analyses were also considered eligible for the replication stage. In common variant GWAS, relaxing the P-value threshold of 5 × 10−8 somewhat in the discovery stage may substantially increase the number of true positive association signals that can be confirmed in replication analyses.[26]

Replication

For the SSI replication cohort, genotyping of the selected SNPs was performed using competitive allele-specific PCR (KASP) chemistry (LGC Genomics). In the Australian cohort, genomic DNA was extracted from whole blood (Qiagen QIAamp DNA Maxi Kit) and selected SNPs were genotyped by PCR amplification (Veriti Thermal Cucler, Applied Biosystems) and bidirectional sanger sequencing (BigDyeTM v.3.1 Terminator Cycle Sequencing Kit and 3730xl DNA Analyzer, Applied Biosystems). The DBDS replication cohort had existing GWAS data available from genotyping with Illumina Global Screening Array and subsequent imputation, as described previously.[27] Association testing was done by logistic regression under an additive genetic model adjusting for sex. With the availability of GWAS data, the replication results from DBDS were also adjusted for the first four principal components. Combined analysis of the discovery and replication stage data was done by fixed-effects inverse variance-weighted meta-analysis. For the 1p31.1 locus, rs3001038 failed in genotyping of the Australian cohort and rs998550 was used as a proxy instead. For the 15q22.31 locus, rs62023078 failed in the replication cohorts from SSI and Australia, and rs13379670 was used as a proxy instead. We considered SNPs with P < 0.05 in the replication stage and P < 5 × 10−8 in the combined analysis to indicate robust evidence of association.

Power analysis

We assessed the statistical power of our study by computer simulations in R.[28] Disease state was simulated from a logistic regression model allowing the odds ratio for a log-additive genetic effect and the frequency of the effect allele to vary and assuming a 3.6% population incidence of febrile seizures before the age of 5 years.[29] For each combination of effect size and effect allele frequency, we simulated 1000 datasets using the study sample size (4502 cases and 51 049 controls). We then conducted association tests on the simulated datasets and calculated power as the proportion of tests with a P-value lower than 5 × 10−8. To compare with statistical power in the previous GWAS of febrile seizures, we also performed simulations using the sample size (1999 cases and 4118 controls) of that study.[15]

Functional annotation

For each of the seven novel and four previously known genome-wide significant loci, we annotated all variants with discovery stage association P < 1 × 10−4 within 500 kb of the lead SNP at the locus using ANNOVAR.[30] Gene-based annotation was performed with gene definitions from the NCBI RefSeq database,[31] and SIFT,[32] PolyPhen-2[33] and FATHMM[34] were used to predict the pathogenicity of missense mutations. Using the FUMA web application,[35] we searched the GWAS Catalog for previously reported associations with P < 5 × 10−8 for SNPs at the febrile seizure loci. We also used FUMA to perform expression quantitative trait locus (eQTL) annotation based on the BRAINEAC, DICE, eQTL Catalogue, eQTLGen and GTEx/v8 databases (see ‘Web resources’ section).

Exome analyses

Exome sequencing data were available for a subset of samples in the iPSYCH cohort. In total, n = 10 967 individuals, sampled from ADHD (n = 2382), ASD (n = 3612), affective disorder (n = 1), bipolar disease (n = 652), schizophrenia (n = 1471) or the population representative sample (n = 2849) were analysed. Of these, 589 were febrile seizure cases. There were more potential febrile seizures controls with exome data available, but within each iPSYCH subset the ratio of febrile seizure cases to controls was approximately balanced at 1:17.6 to avoid spurious associations driven by the iPSYCH diseases. For these samples, and for each of the novel febrile seizure loci, we considered all genes within a 1-MB region of the lead variant and performed gene-based tests for rare-variant association, using the optimal sequence kernel association test (SKAT-O)[36] approach as implemented in EPACTS v.3.2.6, using default settings.

Variance explained and genetic correlation analyses

Assuming a multifactorial liability threshold model, we used risk allele frequencies and estimated odds rations (ORs) from the discovery meta-analysis for the lead SNPs at all 11 genome-wide significant loci to estimate the proportion of variance in the liability to febrile seizures explained by these variants.[37] The estimate of variance explained by all common (MAF >1%) autosomal SNPs (also known as SNP heritability) was calculated based on the discovery stage meta-analysis results using LD score regression.[25] In the calculations of variance explained, we assumed a 3.6% population incidence of febrile seizures before the age of 5 years.[29]

We retrieved GWAS summary statistics for focal epilepsy, genetic generalized epilepsy and all epilepsies combined from the International League Against Epilepsy (ILAE) Consortium on Complex Epilepsies,[38] and used LD score regression to estimate genetic correlations between febrile seizures and epilepsies.

Polygenic risk score analyses

To examine associations between general susceptibility to febrile seizures with other outcomes of interest, we constructed polygenic risk scores (PRS) by summing over risk allele counts/dosages at sets of relevant SNPs, with weights for each SNP corresponding to its estimated effect size. Two sets of scores were constructed. First, we included only the lead variants at each of the seven novel and four known genome-wide significant and robustly replicated loci.

Second, we constructed scores for the individuals in the iPSYCH cohort based on the GWAS summary statistics for febrile seizures from the SSI discovery cohort alone. We generated PRS applying an LD-clumping and P-value thresholding approach with nine different P-value thresholds (in −log10 scale: 0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4) as well as LDpred[39] infinitesimal model and LDpred point-normal mixture with seven different proportions of causal variants (1, 0.3, 0.1, 0.03, 0.01, 0.003, 0.001), yielding a total of 17 sets of PRSs. We then performed out-of-sample validation of how well each set of PRS predicted febrile seizures in the iPSYCH discovery cohort. The best predictive performance (assessed with Nagelkerke R measure) was found for LDpred with a proportion of 0.3% causal variants in the model.

We used logistic regression models to assess the association between PRS for febrile seizures (based on the SSI cohort GWAS results) and risk of epilepsy in the iPSYCH cohort with and without stratification by febrile seizure diagnosis. In these analyses, the larger iPSYCH cohort was divided into individuals with epilepsy (n = 1978; excluded from the main febrile seizures discovery analysis) and a comparison group of individuals without any epilepsy record, who were unrelated, of European ancestry, and not included in the SSI cohort (n = 62 906). The comparison group of individuals without epilepsy included the febrile seizure cases and controls from the main iPSYCH febrile seizures discovery analysis.

To investigate relations between iPSYCH disease groups and genetic susceptibility to febrile seizures, we used a multiple regression model with febrile seizures PRS (calculated for the iPSYCH individuals based on the SSI cohort GWAS results) as the outcome and 14 factors defined by iPSYCH group (ADHD, ASD, schizophrenia, affective disorder, bipolar disease, anorexia, and the population representative sample) and febrile seizure diagnosis (yes/no), and also adjusting for 10 principal components. We used Wald tests for pairwise comparisons of the 14 different groups’ effects on the PRS.

For the hospitalization data for febrile seizure patients (available in the iPSYCH and DBDS cohorts), we used a negative binomial model to analyse the effect of febrile seizures PRS (based on meta-analysis results for the 11 robustly associated loci) on the number of febrile seizures hospital admissions, including sex and four principal components as covariates in the model. We also divided febrile seizure patients into three groups by PRS, corresponding to the first quintile, the second to fourth quintile and the fifth quintile, and conducted time-to-event analysis of age at first febrile seizures hospital admission by PRS group by the Kaplan–Meier method. We used a log-rank test to test for differences between the three PRS groups.

Data availability

GWAS summary statistics from this study will be made available via the Danish National Biobank website (https://www.danishnationalbiobank.com/gwas, accessed 18 November 2021) on publication of the study.

Results

Our study design is illustrated in Fig. 1 and Supplementary Fig. 1. The discovery stage GWAS meta-analysis included 4502 cases and 51 049 controls of European ancestries. Statistical power simulations suggested that our discovery analysis had >80% power to detect variants with ORs down to 1.15–1.25 for common variants (Supplementary Fig. 2). After data cleaning and imputation based on the 1000 genomes project phase 3 reference panel,[40] 6.8 million autosomal variants were analysed for association with febrile seizures. The genomic inflation factor was 1.08 and the LD score regression intercept was 1.03, indicating minimal population stratification. Quantile–quantile and Manhattan plots are shown in Fig. 2. Eight loci reached genome-wide significance (P < 5 × 10−8), including four previously described loci (SCN1A, SCN2A, ANO3, 12q21.33)[15] and four novel loci (PTGER3, IL10, BSN, ERC2; Fig. 2 and Supplementary Fig. 3A–D). Three additional novel loci (GABRG2, MAP3K9, HERC1) had multiple variants associated at 5 × 10−8 < P < 1 × 10−6 and were considered suggestive in the discovery stage (Fig. 2 and Supplementary Fig. 3E–G). For each of the seven novel genome-wide significant or suggestive loci, we selected one SNP for replication stage genotyping in 3133 cases and 32 917 controls from Danish and Australian cohorts (Fig. 1). Furthermore, we performed analyses conditioning on the lead SNP at each locus, but did not identify any additional SNPs fulfilling the criteria for replication genotyping. All seven novel genetic loci were replicated and reached genome-wide significance (P < 5 × 10−8) in the combined analysis (Table 1 and Supplementary Fig. 4). There was no interaction between genotype and sex (Supplementary Table 1) and the additive genetic model fitted the data well (Supplementary Table 2). Also, a sensitivity analysis excluding 57 cases with first hospital admission with febrile seizures at age 5 or 6 years, and 92 cases with later epilepsy or non-febrile seizures in additional follow-up data produced very similar results in comparison to the main analysis (Supplementary Table 3).

Figure 2

Febrile seizures genome-wide association analysis. Manhattan plot of −log10 P-values across the chromosomes (left) and corresponding quantile–quantile plot of observed versus expected −log10 P-values (right). Robustly associated loci are named; previously known loci are indicated in black; novel loci in red. The names above each locus represent the nearest protein coding gene at the locus or the chromosome band if the locus is in an intergenic region.

Table 1

Discovery, replication and combined results for seven novel loci associated with febrile seizures

ChromosomePosition (bp)SNP			Effect allele	Alt. allele	Gene	Sample set	EAF cases	EAF controls	OR (95% CI)	P	I 2	P het
Locus 1p31.1
1 71163715 rs3001038			A	G	PTGER3	SSI	0.527	0.486	1.19 (1.10–1.28)	1.42 × 10−5
						iPSYCH	0.524	0.488	1.15 (1.09–1.22)	8.59 × 10−7
						Combined discovery			1.17 (1.11–1.22)	6.54 × 10−11	0	0.58
						Replication			1.10 (1.03–1.16)	0.00214	37.8	0.2
						All combined			1.14 (1.10–1.18)	1.84 × 10−12	58.9	0.12
Locus 1q32.1
1 206944645 rs1518111			T	C	IL10	SSI	0.215	0.194	1.15 (1.04–1.26)	0.00468
						iPSYCH	0.212	0.186	1.18 (1.10–1.27)	3.65 × 10−6
						Combined discovery			1.17 (1.11–1.24)	4.62 × 10−8	0	0.62
						Replication			1.14 (1.06–1.23)	0.000448	0	0.76
						All combined			1.16 (1.11–1.21)	1.02 × 10−10	0	0.59
Locus 3p21.31
3 49644012 rs11917431			T	C	BSN	SSI	0.308	0.282	1.13 (1.04–1.23)	0.00284
						iPSYCH	0.319	0.286	1.17 (1.10–1.25)	4.6 × 10−7
						Combined discovery			1.16 (1.10–1.22)	4.42 × 10−9	0	0.53
						Replication			1.13 (1.06–1.20)	0.000235	0	0.95
						All combined			1.15 (1.10–1.19)	5.63 × 10−12	0	0.52
Locus 3p14.3
3 56298138 rs4974130			C	T	ERC2	SSI	0.076	0.058	1.34 (1.15–1.55)	0.00019
						iPSYCH	0.072	0.055	1.33 (1.19–1.48)	1.5 × 10−6
						Combined discovery			1.33 (1.22–1.45)	4.24 × 10−10	0	0.95
						Replication			1.27 (1.13–1.44)	6.68 × 10−5	21.1	0.28
						All combined			1.31 (1.22–1.41)	1.52 × 10−13	0	0.57
Locus 5q34
5 161508503 rs1997583			T	A	GABRG2	SSI	0.108	0.092	1.20 (1.06–1.36)	0.00496
						iPSYCH	0.101	0.085	1.24 (1.12–1.37)	4.65 × 10−5
						Combined discovery			1.22 (1.13–1.32)	5.11 × 10−7	0	0.7
						Replication			1.33 (1.21–1.46)	8.8 × 10−9	64.8	0.06
						All combined			1.26 (1.19–1.34)	5.48 × 10−14	42.4	0.19
Locus 14q24.2
14 71323249 rs244620			C	T	MAP3K9	SSI	0.469	0.447	1.10 (1.02–1.19)	0.0151
						iPSYCH	0.485	0.451	1.15 (1.08–1.22)	2.09 × 10−6
						Combined discovery			1.13 (1.08–1.19)	1.33 × 10−7	0	0.42
						Replication			1.11 (1.05–1.18)	0.000384	0	0.54
						All combined			1.13 (1.09–1.17)	2.35 × 10−10	0	0.68
Locus 15q22.31
15 64132818 rs62023078			G	T	HERC1	SSI	0.848	0.823	1.21 (1.09–1.34)	0.000346
						iPSYCH	0.841	0.821	1.15 (1.06–1.24)	0.000337
						Combined discovery			1.17 (1.10–1.25)	7.43 × 10−7	0	0.46
						Replication			1.17 (1.08–1.27)	0.000161	55.8	0.1
						All combined			1.17 (1.11–1.23)	4.97 × 10−10	0	0.99

EAF = effect allele frequency; I2 = heterogeneity estimate; P = association P-value for febrile seizures; Phet = P-value from the Cochran Q-test of heterogeneity.

Two novel loci harbouring fever response genes

At the first locus for febrile seizures on chromosome 1p31.1, rs3001038 yielded the lowest P-value (combined OR = 1.14, 95% CI = 1.10–1.18; P = 1.84 × 10−12). The lead variant rs3001038 was located in an LD block containing the important fever response gene PTGER3[1] (Supplementary Fig. 3A). In search of possible functional mechanisms underlying the association signal, we annotated the 109 variants within 500 kb of the lead SNP that had discovery stage association P < 1 × 10−4 (Supplementary Table 4). Of note, many associated SNPs at the locus were reported as eQTLs for PTGER3 in various brain tissues in the Genotype Tissue Expression (GTEx) project with risk-increasing alleles corresponding to increased PTGER3 expression (Supplementary Table 5).

The association signal at the second locus on chromosome 1q32.1 was sharply defined by a few SNPs intronic or upstream of IL10, which encodes the antipyretic cytokine interleukin 10 (IL10)[41] (Supplementary Fig. 3B). The SNP rs1518111 was taken forward to the replication stage and confirmed (combined OR = 1.16, 95% CI = 1.11–1.21; P = 1.02 × 10−10). Of the six variants with P < 1 × 10−4, two have been reported to associate with the inflammatory disorder Behçet’s disease (P < 5 × 10−8), with the risk allele also corresponding to an increased risk of febrile seizures in our data (Supplementary Table 6). Notably, rs1518111 and other associated SNPs at the locus are all strong eQTLs for IL10 in whole blood with the risk alleles for febrile seizures corresponding to lower IL10 expression (Supplementary Table 5).

Five additional novel loci including genes affecting neuronal excitability

At the third novel locus on chromosome 3p21.31, the SNP rs11917431, which is intronic in BSN yielded the lowest P-value (combined OR = 1.15, 95% CI = 1.10–1.19; P = 5.63 × 10−12). BSN encodes bassoon, a scaffolding protein of the presynaptic active zone.[42] The locus is characterized by extensive LD with SNPs highly correlated to the lead SNP spanning a region of almost 500 kb (Supplementary Fig. 3C). Fourteen of the SNPs at the locus with P < 1 × 10−4 were exonic, including five non-synonymous variants resulting in amino acid changes of USP4, GPX1, BSN, MST1 and RNF123, respectively (Supplementary Table 4). There were also a large number of previously reported GWAS associations at the locus, particularly for cognitive traits, chronic inflammatory diseases and blood protein levels (Supplementary Table 6). The risk alleles for febrile seizures SNPs at this locus were previously positively associated with math ability and other cognitive measures (Supplementary Table 6). A massive number of eQTL associations involving a wide array of genes and tissues have been reported for SNPs at the locus (Supplementary Table 5).

The lead SNP at the fourth locus on chromosome 3p14.3, rs4974130 (combined OR = 1.31, 95% CI = 1.22–1.41; P = 1.52 × 10−13) was intronic in ERC2 with a swathe of highly correlated SNPs within the gene having similar P-values (Supplementary Fig. 3D). ERC2 encodes CAST, another protein of the cytomatrix of the active zone.[43] Four of these SNPs were exonic, including a missense variant in ERC2 (Supplementary Table 4). Furthermore, the lead SNP and many other associated SNPs at the locus were reported as eQTLs for ERC2 in cerebellum in the GTEx project with risk-increasing alleles corresponding to increased ERC2 expression, and a few GWAS Catalog associations were reported (Supplementary Tables 5 and 6).

At the fifth locus for febrile seizures on chromosome 5q34, the lead SNP, rs1997583 (combined OR = 1.26, 95% CI = 1.19–1.34; P = 5.48 × 10−14) was intronic in GABRG2, which encodes the γ2 subunit of gamma-aminobutyric type A (GABAA) receptor[44] (Supplementary Fig. 3E), and where pathogenic missense variants cause febrile seizures in rare dominant families.[11,12] Other associated SNPs at the locus were located within or near the gene; one was exonic (synonymous) and none were reported as eQTLs or in the GWAS Catalog (Supplementary Tables 4–6).

The sixth locus comprised associated SNPs in an LD block on chromosome 14q24.2 containing the gene MAP3K9 encoding mitogen-activated kinase kinase kinase 9 (Supplementary Fig. 3F). The lead SNP, rs244620 (combined OR = 1.13, 95% CI = 1.09–1.17; P = 2.35 × 10−10) has been reported as an eQTL for MAP3K9 in a range of different GTEx tissues (Supplementary Table 5) with the risk allele rs244620-C corresponding to increased expression in most tissues. Associated SNPs at the locus also showed strong cis-eQTL effects for MAP3K9, TTC9 and MED6 in whole blood samples from the eQTLGen Consortium also with the febrile seizures risk alleles corresponding to increased expression (Supplementary Table 5).

Finally, at the seventh locus on chromosome 15q22.31, rs62023078 (combined OR = 1.17, 95% CI = 1.11–1.23; P = 4.97 × 10−10) was intronic in HERC1, which encodes HERC1, an E3 ubiquitin ligase located at the presynaptic terminal[45] (Supplementary Fig. 3G). Among SNPs at the locus with P < 1 × 10−4, many have been reported as cis-eQTLs in whole blood and other tissues with febrile seizures risk alleles corresponding to decreased expression of HERC1, USP3, DAPK2 and PLEKHO2, and increased expression of FBXL22, APH1B and RBPMS2 (Supplementary Table 5).

Exome sequencing data were available in the iPSYCH cohort for 589 febrile seizure cases and 10 378 controls. At each locus and for all genes within 1 Mb of the top variant, we used these data to conduct gene-based tests for aggregated effects of rare coding variants (optimal sequence kernel association test, SKAT-O; see the ‘Materials and methods’ section). Among the 159 genes tested, we were not able to detect any locus-wide significant gene-based rare-variant association to febrile seizures, except for an association for SNX22 at the HERC1 locus, which was, however, driven by only 0.36% of individuals carrying rare variants in the gene (Supplementary Table 7).

Previously reported loci

Four loci associated with febrile seizures in general (SCN1A, SCN2A, ANO3 and 12q21.33) are already reported based on the SSI cohort included in this meta-analysis.[15] These were all corroborated by analysis of the independent iPSYCH cohort (Supplementary Table 4) and functional annotation revealed many entries in the GWAS Catalog as well as eQTLs for several genes, including the voltage gated Na+ channel genes SCN1A and SCN2A (Supplementary Tables 5 and 6). A comparison of effect sizes against MAFs of the seven novel and four known genome-wide significant loci reflects theoretical expectations from simulations of statistical power based on sample sizes of the previous and the current GWAS (Fig. 3). At a given MAF, the effects of the novel loci are smaller than those of the known loci, and further undetected loci are likely to have even smaller effects. An exception to this pattern may be loci for febrile seizures elicited by a specific exposure. In our previous study, we detected two loci (harbouring the interferon-stimulated gene IFI44L and the measles virus receptor gene CD46) that were distinctly associated with febrile seizures occurring as a rare adverse event after MMR vaccination.[15] Vaccination data were not available for cohorts other than SSI, but most febrile seizure episodes have no temporal connection to MMR vaccination,[46] and analysis of the IFI44L and CD46 variants in the iPSYCH and DBDS cohorts showed little evidence of association (Supplementary Table 8). Furthermore, none of the seven novel loci differed between MMR-related febrile seizures and MMR-unrelated febrile seizures (Supplementary Table 9).

Figure 3

MAFs and effect sizes (OR in discovery stage) for 11 robustly associated febrile seizure loci. Previously known loci are indicated in blue; novel loci in red. The red dotted lines represent ORs corresponding to 80% statistical power based on simulations with the sample size of the current study. The blue lines represent ORs needed for 80% power in simulations with the sample size of the previous febrile seizures GWAS.[15]

Heritability and shared susceptibility with epilepsy

To assess the combined effects of loci, we used a multifactorial liability threshold model[37] and found that the lead SNPs at the 11 robustly associated loci collectively explained 2.8% of the variance in liability to febrile seizures. Looking more broadly across the genome, the estimated proportion of variance in febrile seizures liability explained by all autosomal SNPs with a MAF >1% (SNP heritability) was 10.8% (SE 1.9%).

Given the fact that rare variants in GABRG2, SCN1A and SCN2A at three of the febrile seizure loci have been linked to a range of epilepsy syndromes[47,48] and the epidemiological links between febrile seizures and epilepsy, we wanted to examine potential shared genetic susceptibility. First, we performed genetic correlation analyses between our febrile seizures results and summary statistics from the latest GWAS mega-analysis from the ILAE Consortium on Complex Epilepsies.[38] These analyses revealed positive genetic correlations with focal epilepsy (r = 0.59, SE = 0.24, P = 0.01), genetic generalized epilepsy (r = 0.23, SE = 0.07, P = 0.002) and all epilepsies combined (r = 0.39, SE = 0.10, P = 1.68 × 10−4).

Next, we included an additional subset of the iPSYCH cohort of patients diagnosed with epilepsy. These epilepsy patients (n = 1978) were excluded from the main discovery analyses of the iPSYCH cohort. The comparison group comprised unrelated, European ancestry iPSYCH individuals without epilepsy, who were also not in the SSI cohort (n = 62 906). Using these data, we assessed the association of epilepsy as outcome to a PRS for febrile seizures based on the SSI discovery cohort alone. One unit increase in the standardized PRS was associated with an increased risk of epilepsy (OR 1.07, 95% CI 1.02–1.12, P = 2.64 × 10−3; nno_epilepsy = 62 906, nepilepsy = 1978), which remained significant when correcting for febrile seizures case status (OR 1.05, 95% CI 1.00–1.10, P = 0.045). Among the 1978 epilepsy patients, 288 also had a history of febrile seizures. Stratification of the data into individuals with and without a history of febrile seizures revealed no difference in effect of PRS on risk of epilepsy (P = 0.93).

Neuropsychiatric outcomes

iPSYCH is designed as a case-cohort study and consists of a random population sample, and groups of patients with neuropsychiatric diagnoses, including ADHD, affective disorder, anorexia, ASD, bipolar disease and schizophrenia. We tested whether the febrile seizure association of each lead variant in the seven novel and four known loci was constant across these different iPSYCH subgroups and when correcting for the 11 loci tested, no statistically significant heterogeneity could be claimed. To generalize this analysis, we used a PRS based on the SSI discovery cohort alone and regressed the PRS on 14 factors defined by iPSYCH group and febrile seizure diagnosis. Figure 4A shows that the adjusted mean PRS for febrile seizures among iPSYCH individuals without febrile seizures was close to the overall adjusted mean PRS. For iPSYCH individuals with a febrile seizure diagnosis, the adjusted mean PRS was elevated in five of seven groups. The largest effect was seen for the iPSYCH population representative sample with febrile seizures, which in pairwise comparisons had a higher mean adjusted PRS than all other groups except the group defined by anorexia and febrile seizures (Fig. 4B).

Figure 4

Multiple regression of PRS for febrile seizures on iPSYCH group and febrile seizure diagnosis. (A) Effect estimates of adjusted mean PRS for 14 factors defined by iPSYCH group and febrile seizure diagnosis with 95% CI. P-values are from tests of difference from overall mean adjusted PRS. (B) P-values from Wald tests of the equal group effect in pairwise comparisons of the 14 groups. ADHD = attention deficit hyperactivity disorder; ASD = autism spectrum disorders; Schizo = schizophrenia.

Analysis of hospitalization data

Next, we assessed the combined impact of the identified febrile seizure loci on (i) number of hospital admissions due to febrile seizures; and (ii) the age at first admission with febrile seizures. We constructed a PRS for febrile seizures based on the 11 robustly associated loci and analysed the association between the PRS and the number of additional febrile seizures in cases using negative binomial regression (see the ‘Materials and methods’ section). In the iPSYCH discovery cohort each unit increase in standardized PRS was associated with a 1.17-fold increase in expected additional febrile seizures hospital admissions (negative binomial test, P = 6.47 × 10−4; n = 2511 cases). This analysis could not be done in the SSI discovery cohort, due to the cohort’s special ascertainment scheme, but a similar estimate (1.16-fold increase in expected additional febrile seizures per PRS unit) was found in the DBDS replication cohort (negative binomial test, P = 0.07; n = 934 cases). The distribution of number of febrile seizures hospital admissions among cases and PRS by number of hospital admissions is seen in Fig. 5, which also shows that febrile seizure patients with a high polygenic risk were more likely to be hospitalized with febrile seizures at a younger age.

Figure 5

Association between PRS for febrile seizures and history of hospital admissions with a febrile seizure diagnosis among 2511 febrile seizure cases in the iPSYCH cohort. (A) Number of separate admissions to hospital with a febrile seizure diagnosis among the 2511 cases. (B) Box plots of PRS for febrile seizures based on 11 genome-wide significant loci by number of febrile seizures hospital admissions. The boxes define the first and third quartiles, with a black line representing the median. Whiskers extend to the smallest and largest values no further than 1.5 times the interquartile range from the first and third quartiles, respectively. The mean PRS is marked with an ‘x’. (C) Kaplan–Meier curves of time-to-event analysis for first febrile seizures hospital admission by groupings of PRS. A log-rank test was used to test for difference between the three curves.

Discussion

In this genome-wide meta-analysis of febrile seizures including a total of 7635 cases and >80 000 controls, we identified seven novel robustly associated loci and confirmed four previously reported loci. Our results implicate central fever response genes as well as loci related to ion channel function, neurotransmitter release and binding, and vesicular transport and membrane trafficking at the synapse. We found positive genetic correlations between febrile seizures and epilepsy syndromes, and we found that higher PRSs for febrile seizures were associated with epilepsy and with history of hospital admission for febrile seizures. We also found that polygenic risk of febrile seizures was lower in febrile seizure patients with neuropsychiatric disease compared to febrile seizure patients in a general population sample.

Febrile seizures are an outcome of intricate aetiology involving feedback loops between the body’s fever response and the neuronal circuitry. It has been suggested that the regulation of body temperature may be different in children susceptible to febrile seizures.[49,50] Our association findings at the PTGER3 and IL10 loci are consistent with this hypothesis. PTGER3 encodes EP3, one of four receptors for prostaglandin E2 (PGE2), which is considered to be a major pyrogenic mediator of fever.[1] It has been shown that mice lacking the EP3 receptor fail to show a febrile response to direct stimulation with PGE2 and to exposure to endogenous (IL-1β) and exogenous (lipopolysaccharide, LPS) pyrogens.[51] Additional experiments established that selective disruption of Ptger3 in the median preoptic nucleus of the mouse brain was sufficient to block fever in response to PGE2 and LPS.[52] While the median preoptic nucleus is not assayed in GTEx, we found that risk alleles for the SNPs at the locus corresponded to higher expression of PTGER3 in several other brain regions. We hypothesize that the genetic upregulation of PTGER3 leads to a more pronounced fever response, which in turn increases the susceptibility to febrile seizures. Furthermore, we found that febrile seizure risk alleles for SNPs at the IL10 locus corresponded to lower IL10 expression in blood. IL10 encodes the anti-inflammatory cytokine IL10, which in the complex landscape of cytokine signalling acts as a central endogenous antipyretic.[41] Thus, the genetic downregulation of IL10 could conceivably lead to a more pronounced fever response and thereby confer increased risk of febrile seizures.

Other novel loci harboured genes such as BSN, ERC2 and HERC1, which have functions related to neurotransmitter release, vesicular transport and membrane trafficking. Presynaptic plasma membranes contain specialized regions called active zones, where synaptic vesicles release neurotransmitters, such as glutamate or GABA, depending on the type of neuron.[42] The release sites are characterized by a dense cytomatrix at the active zone (CAZ) organized by a number of proteins, including the large scaffolding protein bassoon encoded by BSN. Disruption of bassoon function in mice is associated with pronounced spontaneous seizures, and is used as an experimental model for epilepsy.[53,54] Bassoon interacts with multiple other CAZ proteins, including CAST (also known as ELKS2) encoded by ERC2.[43] Uncovering the complex roles of CAST, bassoon and other proteins in the assembly of the CAZ and in synaptic transmission in different types of neurons is an active field of research.[42,43,55] Post-translational modifications such as ubiquitination add to the complexity of presynaptic dynamics, and recent mouse studies show the importance of the HERC1 E3 ubiquitin ligase in the regulation of vesicular transport and membrane trafficking at the synapse.[45,56] While these observations indicate BSN, ERC2 and HERC1 are potential candidates at the loci on chromosomes 3p21.31, 3p14.3 and 15q22.31, follow-up studies are required to establish causal links and to illuminate mechanisms related to febrile seizures.

Our findings reinforce the notion that genetic studies of isolated febrile seizures can be informative for understanding epilepsy and seizure genesis in general. Previously, rare familial missense mutations in GABRG2 have been found to segregate with febrile seizures and childhood absence epilepsy or genetic epilepsy with febrile seizures plus (GEFS+),[11,12] and nonsense mutations of the gene have been found in more severe epilepsy phenotypes.[48] GABA is the major inhibitory neurotransmitter acting at GABAA receptors, responsible for fast synaptic inhibition in the brain,[44] and it is thought that the variability of epilepsy phenotypes associated with GABRG2 mutations is related to the extent of GABA signalling impairment.[48] Extending this line of thinking, we hypothesize that the common GABRG2 variants identified here increase febrile seizures risk through subtle effects on GABA signalling. We also confirmed common variant associations for febrile seizures in the well-known epilepsy genes SCN1A and SCN2A. Previous studies have linked rare missense mutations in these voltage gated Na+ channel genes to a range of epilepsy syndromes, some involving febrile seizures.[47] Looking across the genome, we found positive genetic correlations between febrile seizures and focal epilepsy, genetic generalized epilepsy and all epilepsies combined. Delineating genetic differences and similarities between febrile seizures and epilepsies for common variants as well as using approaches focused on rare genetic variation[57] is an important area for future research.

From a clinical perspective, an enhanced understanding of potential etiological relations between febrile seizures and epilepsies is highly warranted. A central unresolved question is whether febrile seizures may in some cases lead to hippocampal damage, provoking a chronic process of epileptogenesis or whether shared genetic and developmental factors increase the susceptibility to febrile seizures as well as epilepsy.[10,58] Population registers and large biobanks, such as those available in Denmark, offer opportunities for carefully designed massive-scale studies addressing this question through genetic analyses of febrile seizures and epilepsy patients with long-term follow-up. Such studies may also enable prediction of individuals at high risk of febrile seizures and/or epilepsy, and possibly suggest preventive strategies. While the current study was not sufficiently powered to develop prediction tools, the fact that a score based on all genome-wide significant loci was associated with number of hospital admissions for febrile seizures and age at first admission, suggests future potential clinical utility of a deeper understanding of the genetic underpinnings of febrile seizures.

Epidemiological studies have reported associations between epilepsy and, to a lesser degree, febrile seizures and risk of psychiatric disorders.[5,59-61] However, somewhat surprisingly, an adjusted PRS for febrile seizures based on the SSI discovery cohort was on average lower in all iPSYCH disease groups except anorexia, compared to the iPSYCH population representative sample. While this question should be revisited in a later study with a larger sample for PRS construction, it raises the possibility that shared environmental or developmental risk factors or rare genetic variants may underlie epidemiological associations between febrile seizures and psychiatric disorders rather than shared common genetic variation. Interestingly, a large-scale genetic correlation study also found limited sharing of common genetic variants between epilepsy and psychiatric disorders, despite strong epidemiological associations.[62]

Our study had limitations, including the restriction to individuals of European ancestries. Reported incidence rates of febrile seizures vary considerably around the world,[2] and differences in allele frequencies and LD structure at associated loci also makes it inappropriate to generalize our findings to other ancestries. Thus, further studies are warranted to characterize genetic susceptibility to febrile seizures in non-European populations. A second limitation was the granularity of the outcome definition. With a positive predictive value of a diagnosis of febrile seizures (ICD-8 and ICD-10) of 93%, the Danish National Patient Register is an excellent resource for research,[63] but the register does not contain information about complex febrile seizures including characteristics such as the duration and type (focal or generalized) of febrile seizures, postictal neurological abnormalities and possible recurrence within 24 h. Further studies are required to understand effects of variants at the identified loci across a phenotypic spectrum of thoroughly characterized febrile seizure patients. Information about the specific infectious exposure or inflammation leading to fever and febrile seizures would also be revealing, but is usually difficult to ascertain.

In conclusion, this largest genetic investigation of febrile seizures to date implicates central fever response genes as well as genes affecting neuronal excitability, including several known epilepsy genes. This adds to knowledge from clinical, epidemiological and laboratory studies that febrile seizures are the manifestation of a highly complex physiological process involving environmental triggers, neurodevelopmental vulnerability and genetic susceptibility. Furthermore, the positive genetic correlation between febrile seizures and both focal and generalized epilepsies underscores that carefully conducted genetic studies of febrile seizures can yield fundamental insights into the genetic regulation of neuronal excitability highly relevant to our understanding of epileptogenesis. Clearly, additional functional and genetic studies, including of rare as well as common variants at the new loci we have identified, are required to enhance our understanding of the pathophysiological pathways and triggers of seizures with and without fever.

Web resources

All URLs were accessed 18 November 2021: 1000 Genomes Project, http://www.1000genomes.org/; ANNOVAR, http://annovar.openbioinformatics.org/; Braineac, http://www.braineac.org/; DICE, https://dice-database.org/; eQTL Catalogue, https://www.ebi.ac.uk/eqtl/; eQTLGen Consortium, https://www.eqtlgen.org/; FUMA, https://fuma.ctglab.nl/; GTEx, https://www.gtexportal.org/; iPSYCH, http://ipsych.au.dk/about-ipsych/; LD score regression, https://github.com/bulik/ldsc/; METAL, http://www.sph.umich.edu/csg/abecasis/metal/; Michigan Imputation Server, https://imputationserver.sph.umich.edu; NHGRI-EBI GWAS Catalog, https://www.ebi.ac.uk/gwas/; R software, http://www.r-project.org/; SNPTEST, https://mathgen.stats.ox.ac.uk/genetics_software/snptest/snptest.html.

Funding

The study was supported by grants from the Danish Medical Research Council (0602-01818B), the Oak Foundation (OCAY-18–598), the US National Institutes of Health (NIH)/National Institute of Allergy and Infectious Diseases (R01AI093697), the Novo Nordisk Foundation Challenge programme (NNF17OC0027594), a Lundbeck Foundation Ascending Investigator grant (R313-2019–554) to B.F., and National Health and Medical Research Council (NHMRC) Program Grant (1091593) to S.F.B. and I.E.S, Practitioner Fellowship (1006110) to I.E.S., and R.D Wright Career Development Fellowship (1063799) to M.S.H., and a Novo Nordisk Foundation Hallas-Møller grant to A.H. The Danish National Biobank was established with the support of major grants from the Novo Nordisk Foundation, the Danish Medical Research Council and the Lundbeck Foundation. L.S. received support from a Carlsberg Foundation postdoctoral fellowship (CF15-0899); X.L. reports funding from the Nordic Center of Excellence in Health-Related e-Sciences; K.B. received support from the Novo Nordisk Foundation (NNF14CC0001, NNF17OC0027594); D.W. reports funding from the Novo Nordisk Foundation (NNF18SA0034956, NNF14CC0001, NNF17OC0027594); T.H.P. acknowledges the Novo Nordisk Foundation (NNF18CC0034900) and the Lundbeck Foundation (R190-2014–3904); J.C. and J.W.D. report funding from the Novo Nordisk Foundation (NNF16OC0019126), The Danish Epilepsy Association and the Central Denmark Region; C.A. is supported by The Danish National Research Foundation (Niels Bohr Professorship to John McGrath). The content is solely the responsibility of the authors and does not necessarily represent the official views of the funders.

Competing interests

The authors report no competing interests.

Click here for additional data file.