Multi-ancestry genetic study in 5,876 patients identifies an association between excitotoxic genes and early outcomes after acute ischemic stroke

Study Design

A detailed description of the acute ischemic stroke patients recruited from 21 sites from seven countries throughout the world, has been published elsewhere.⁶ Briefly, adult acute ischemic stroke patient with measurable deficit on the NIHSS that presented within 6 hours of stroke onset (or last known normal) were enrolled in the study after obtaining informed consent, including patients treated with tPA. Patients who underwent a thrombectomy, were enrolled in other treatment trials, or for whom consent and/or a blood sample could not be obtained were excluded. Demographics, co-morbidities, acute treatment variables, imaging data and TOAST classification were collected.

To accommodate the difference in the genetic architecture intrinsic to the country of origin, we performed a three-stage analysis (Figure 1A). First, we used an additive model to perform a GWAS in each country individually, except for the United States, where the population was stratified into European and African ancestry cohorts. We then performed a fixed effects meta-analyses within the same ethnic cohorts. Finally, we used a multi-ancestry Bayesian meta-analysis to collapse all the ethnic backgrounds. Unlike a fixed effect meta-analysis, the Bayesian approach is able to account for population structure differences.¹⁹ Genetic loci that passed multiple test correction, a threshold set at Log Bayes Factor (LBF) > 5, were annotated using bioinformatics tools to identify the gene driving the genetic signal (Figure 1B). We used functional annotation, multi-tissue expression quantitative trait loci (eQTL) data, and summary-data-based Mendelian randomization (SMR) to map the genome-wide to specific genes. Single nuclei-RNA-seq data derived from cortex samples was used to determine potential correlation between the transcripts of the identified genes and determine in which brain cell types the genes are expressed.²⁰

• Download figure
• Open in new tab

Figure 1. Study design.

Summarized description of A. the multi-step approach used to account for the genetic heterogeneity intrinsic to the multi-ancestry nature of the GENISIS study. We performed single variant analysis in each of the participating countries separately. Then we meta-analyzes all the non-Hispanic whites (blue) and Hispanic (green) ethnicities. Finally, we analyzed the non-Hispanic whites, Hispanics, Korea (orange) and US participants with African descent (US AfA – yellow) using a Bayesian model. The variants with genome-wide significant or suggestive results were annotated using B. sequential steps to elucidate the gene driving the association. We performed gene-based and pathway analyses, we collected the information available in publicly available datasets and we performed Mendelian randomization. We also performed genetic architecture overlap tests to test if there was any overlap with known risk factors.

The study was approved by the Institutional Review Boards at every participating site. Written informed consent was obtained from all participants or their family members. All research was performed according to the approved protocols and consents.

Genotyping

All participants were genotyped using Illumina SNP array technology. Samples were genotyped in seven batches during the GENISIS recruitment (see Supplementary Methods). Genotyping quality control and imputation were performed separately for each genotyping round using SHAPEIT²¹ and IMPUTE2.²² For each genotyping batch, SNPs with a call rate lower than 98% and autosomal SNPs that were not in Hardy-Weinberg equilibrium (P<1×10⁻⁰⁶) were removed from the dataset. The X chromosome SNPs were used to determine sex based on heterozygosity rates, and samples with discordant inferred sex and reported sex were removed. Only samples with call rate greater than 98% were considered to pass quality control. Finally, the genotype batches were merged in a single file to perform the analyses.

Additional QC was performed in the merged dataset. We tested pairwise genome-wide estimates of proportion identity-by-descent, the presence of unexpected duplicates, and cryptically relatedness (PI-HAT>0.30). Of the pairs of these samples flagged, the sample with higher genotyping rate was kept for downstream analysis. Principal component analysis (PCA) was performed using HapMap as an anchor to remove ethnic outliers and keep the populations as homogeneous as possible for each of the participant countries. Principal components were also used to cluster and identify ancestry populations for US participants with European descent (EuA) and African-American descent (AfA). Samples outside two standard deviations from the center of the Non-Hispanic White or the Asian cluster were considered outliers for Spain, Finland and Poland. We confirmed the ethnicity of the AfA and Hispanic populations, however, due to the genetic heterogeneity present in these populations we did not remove the samples outside two standard deviations from the mean.

Analysis of variance

We used genome-wide complex traits analysis (GCTA) to determine the heritability of ΔNIHSS.²³ GCTA estimates the amount of phenotypic variance in a given complex trait explained by all the SNPs and fits the effects of these SNPs as random effects in a linear mixed model. Because it relies on a large, homogeneous populations for accurate results, we only included the individuals with non-Hispanic White ancestry.

Single Variant Analyses

To mitigate the effects of genetic heterogeneity due to the diverse ancestry of participants enrolled in the GENISIS study, we used a multi-step study design (Figure 1A). First, we performed single variant analyses each participant country separately. We tested the association of SNPs across the genome with ΔNIHSS using an additive linear model with PLINK 1.9.²⁴ Sex, age, and the two Principal Components calculated for each population were included in the model. Additional covariates include the SNP genotyping batch, TOAST classification (using dummy variables to incorporate all subtypes), and baseline NIHSS to adjust for stroke severity. Although baseline NIHSS was used to calculate ΔNIHSS, it does not fully explain the observed variance in ΔNIHSS; further, there is no multicollinearity between these two variables, permitting their inclusion in the model.²⁵ Fifty four percent of study participants were treated with tPA, however treatment was not included as covariate. While tPA influences ΔNIHSS⁶, it has little impact on the genetic architecture of ΔNIHSS. The inclusion of tPA as covariate lead to highly correlated results when compared to a model without tPA (r²=0.96, Supplementary Figure 1). Second, we meta-analyzed the populations with similar ethnic backgrounds using with fixed effect meta-analyses using METAL.²⁶ We performed two meta-analyses, one for the non-Hispanic Whites (Spain, Finland, Poland, and United States EuA) and one for the Hispanics (Costa Rica and Mexico). Finally, we analyzed the four available ethnicities non-Hispanic Whites (meta-analysis), Hispanics (meta-analysis), Asians (Korea) and African Americans (United States AfA) using MANTRA, a Bayesian-based multi-ancestry meta-analyses.¹⁹ Log Bayes Factor (LBF) greater than 5 was considered to be genome wide significant after multiple test correction.

Functional Annotation

We annotated all the variants with suggestive associations (LBF>4) with ANNOVAR²⁷ and SnpEff²⁸ to identify the nearest gene and to determine if any variant is predicted to change protein sequence (non-synonymous variants) or could affect expression. We also confirmed if any of the SNPs were possible regulatory elements or DNA features using RegulomeDB.²⁹

DEPICT³⁰ and FUMA³¹ were used to perform gene ontology and pathways analyses. We also leveraged brain single nuclei RNA expression data (http://ngi.pub/snuclRNA-seq/)²⁰, to determine if the gene expression of the genes located in each identified loci was expressed in brain. For the ones expressed in brain, we also investigated if they were expressed in any specific brain cells (Figure 1B). Finally, we accessed blood RNA expression data taken at different times after stroke onset (3h, 5h and 24h) from the CLEAR trial³² (NCT00250991 at www.Clinical-Trials.gov) to test if the expression of genes located in the identified loci were associated with NIHSS or ΔNIHSS (NIHSS_5h-NIHSS_24h). We extracted the correlation between ΔNIHSS and gene expression (meausered using Affymetrix U133 Plus 2.0 array).³³

eQTL mapping, Mendelian Randomization and Colocalization

To identify the most likely functional gene, we accessed available expression quantitative trait (eQTL) datasets: the Genotype-Tissue Expression (GTEx) Project V8 (accessed on 10/10/2019), the Brain eQTL Almanac (Braineac) and an in-house dataset that includes brain expression data for 613 brains³⁴. We used the summary-data-based Mendelian Randomization (SMR)³⁵ and colocalization³⁶ to test for pleiotropic association between the expression level of a gene and a complex trait to evaluate if the effect size of a genetic variant on the phenotype is mediated by gene expression (Figure 1B). We tested GWAS-significant and -suggestive loci from the ΔNIHSS analysis in two datasets: selected GTEx tissues (brain anterior cortex, cerebellum, brain cerebellar hemisphere, substantia nigra, hippocampus, frontal cortex, and putamen) and the Westra et al. dataset³⁷ derived from whole blood. Both SMR and colocalization require effect sizes and the respective standard error to test the causal relationship, but MANTRA does not provide effect sizes. As a consequence, we used the summary statistics from the joint analysis for all populations to perform these analyses that are correlated with the results from MANTRA (r=-0.57; p<1.07×10⁻⁰⁵ – data not shown). To complement the Mendelian randomization analyses with the posterior probability of a variant being causal in both GWAS and eQTL studies accounting for the genetic heterogeneity and linkage disequilibrium (LD), we used eCAVIAR³⁸ which will consider several variants within the GWAS significant loci to perform the test.

Genetic Correlation

We examined similarities in the genetic architectures of stroke early outcomes (ΔNIHSS) and stroke risk¹⁵ using PRSice³⁹, LDSC⁴⁰ and GNOVA⁴¹ (Figure 1B). Briefly, PRSice calculates polygenic risk scores at different p value thresholds by weighting each SNP by their effect size estimates. SNPs present in one dataset, ambiguous SNPs (A/T or C/G) and all SNPs in LD are removed prior to polygenic risk score calculation. LDSC and GNOVA estimate the genetic covariance and the variant-based heritability for two sets of summary statistics, each one corresponding to one trait of interest. These two parameters are used to calculate the genetic correlation and covariance respectively between the two traits. We limited our comparisons to the non-Hispanic White population to keep the population genetically homogeneous and use the 1000 Genomes European population-derived reference dataset. We calculated the genetic correlation between the European ischemic stroke summary statistics of the MEGASTROKE¹⁵ study and the non-Hispanic Whites meta-analysis summary statistics from the GENISIS study. We also determined if traits related to cardiovascular and general health (age at death⁴², lipid levels⁴³ and body mass index (BMI)⁴⁴) are genetically correlated to ΔNIHSS.

Summary statistics of the GENISIS dataset used for these analyses are available upon request. Individual data for the full GENISIS dataset will be uploaded to dbGAP titled: “Genetics of Early Neurological Instability After Ischemic Stroke (GENISIS)”.

Identification of novel loci associated with stroke early outcomes

We performed single variant analyses for each individual cohort separately; then we combined cohorts with similar ethnic backgrounds; finally, we performed a multi-ancestry meta-analysis with the four ethnic groups available in this study (Non-Hispanic Whites, Hispanics, Asians and African Americans) (Figure 1A). We identified seven GWAS significant loci (Figure 2A and Table 2) associated with ΔNIHSS.

Three independent loci were identified in chr2. The first locus, tagged by rs58763243 (MAF_G=0.07; LBF=6.51), was located in a region comprised by several long non-coding RNAs and microRNAs (Supplementary Figure 3A). For this locus, all of the populations contributed to the association with negative betas, indicating that the minor allele was associated with lower (or more negative) ΔNIHSS. In addition this locus reached genome-wide significance in the US AfA population and was nominally significant in the Finnish population (Supplementary Figure 3B).

• Download figure
• Open in new tab

Figure 3. Gene prioritizing summary.

Summary table showing the seven genome-wide significant loci from the multi-ancestry analysis (first column), the total number of genes identified in each of the locus (second column) and gene name for genes for which we have found some kind of evidence (third column). We have included the results from the gene-based analyses, the presence of any eQTL in GTEx portal or Braineac for any of the genome-wide or suggestive variants, if the gene is differentially expressed in neurons according to the snRNAseq data and the results from Mendelian randomization using Westra dataset (whole blood) or GTEx portal (all tissues). Black dots indicate that the gene was not found, red is that it was found but was not significant, yellow it was moderately significant (0.05<p<1×10⁻⁰³) and green shows a significant association (p<1×10⁻⁰³). Shape indicates the direction of effect, round positive, diamond negative and hexagon unknown.

The second locus, rs13403787 (MAF_A=0.16; LBF=5.57), was also located on chr2 in a region with more than 20 genes between two recombination sites (Supplementary Figure 3C). The minor allele was associated with higher (or more positive) ΔNIHSS in all cohorts (Supplementary Figure 3D). The last genome wide significant locus in chr2 was rs72958644 (MAF_T=0.04; LBF=6.34), located in a region that includes ADAM23, CREB1, DYTN, NRP2, MDH1B, among many others (Figure 2C). The signal is driven by the Non-Hispanic Whites (meta-analysis p=3.44×10⁻⁰⁸), but virtually all ethnic groups contributed to this association, as the directionality was consistent across Hispanic, non-Hispanic White and AfA ethnic groups (Figure 2D and Supplementary Table 1). However, the SNPs in this locus were monomorphic in the Asian population. We observed a trend for positive correlation between ΔNIHSS with the genotype in this locus (R²=0.005, p=0.04; Supplementary Table 2).

Four additional loci were identified outside chr2. One locus at 4q34.3, rs12641856 (MAF_T=0.05; LBF=5.50), in a region that includes four genes (LINC00290, MGC45800, MIR1205 and TENM3; Supplementary Figure 3E). The signal was driven mainly by the Korean cohort (p=1.38×10⁻⁰⁷) (Supplementary Figure 3F). The locus identified on chr5 is located on a region containing nine genes (LOC101927134, GRIA1, FAM114A2, SAP30L, SAP30L-AS1, MFAP3, GALNT10, HAND1 and MIR3141; Figure 2F). The minor allele for the top hit in this locus, rs114248865 (MAF_T=0.05; LBF=5.29) was associated with greater (more positive) ΔNIHSS across all cohorts, and was significant in the Spanish (p=8.41×10⁻⁰⁷) and Finnish cohorts (p=0.03; Figure 2G). Another locus on chr6, tagged by rs6930598 (MAF_T=0.06; LBF=5.30), was located in a region with the genes PARK2, PACRG, MAP3K4, AGPAT4, AGPAT4-IT1 and LOC101929239 (Supplementary Figure 3G). The variants in the region were significant in the Mexican and Korean cohorts but the direction of effect was not consistent (Supplementary Figure 3H). Moreover, the MAF for these variants ranged between 16% in the AfA population to 6% in the Asian and Finnish populations, suggesting that the region is very polymorphic depending on ethnicity. Thus, even though the locus is important for ΔNIHSS, it is possible that it is not the causal variant. Finally, we identified a locus on chr7, tagged by the variant rs10807797 (MAF_G=0.34; LBF=5.70), and located in a gene rich region with 15 genes, including TWISTNB, MACC1, TMEM196, ABCB5, RPL23P8 (Supplementary Figure 3I). This locus is tightly encompassed by two recombination sites. The top signal was significant or suggestive in all populations except the Polish and Mexican cohorts. Consistently, the direction of effect was the same in all cohorts except the Mexican cohort (Supplementary Figure 3J, Supplementary Table 1).

Genetic contribution to early outcomes after ischemic stroke

We used GCTA to quantify the phenotypic variance explained by common SNPs. Because GCTA exploits LD patterns to calculate the explained variance, we restricted our analysis to non-Hispanic Whites. Due to founder effects present in the Finnish population, we also removed this cohort from the variance calculation (Final N=4,573). GTCA revealed that common genetic variants explained 8.7% of the variance of ΔNIHSS (p=0.001), confirming that genetic variants and genes are implicated on stroke outcomes. Next we determine what proportion of the genetic component is explained by the GWAS signals.

The SNPs comprised within the 7 genome-wide significant loci, defined as 500 bp upstream or downstream of the top signal, explained 2.1% of the total variance (p=1.80×10⁻⁰⁶) of ΔNIHSS, or just 24.1% of the genetic component of ΔNIHSS. This suggests that there are additional loci associated with ΔNIHSS yet to be discovered. Thus, studies with larger sample size and more statistical power are needed to identify these additional loci.

Functional Annotation of the Genome-Wide Significant Loci

Identifying the likely causal gene from each loci driving the association is a multi-step process (Figure 1B). We first annotated the suggestive variants (LBF>4), but none of them were predicted to change the protein sequence, comprise a regulatory element, or affect the chromatin architecture. Next, we explored publicly available datasets to investigate if any of the SNPs with suggestive LBFs were eQTLs (Supplementary Table 3). We performed gene-based analyses (Supplementary Table 4) and Mendelian Randomization (MR) analyses to identify possible causal relationships between gene expression and ΔNIHSS (Supplementary Table 5). Summary results can be found in Figure 3.

Gene-based analyses using FUMA suggested that DYTN (p=2.86×10⁻⁰⁵, Z=4.02) and ADAM23 (p=1.59×10⁻⁰⁴, Z=3.60) were the genes driving the association at 2q33.3. Several variants in the ADAM23 region were strong eQTLs for this gene in multiple tissues, based on the GTEx data (esophagus mucosa: p=1.90×10⁻⁰⁶; Cultured Fibroblasts: p=3.70×10⁻⁰⁵). Several SNPs in this region were also eQTLs for NDFUS1 (thyroid p=1.20×10⁻⁰⁴) although with less significant p-values than those for ADAM23. The variant was also significant in an independent dataset from the Braineac eQTL study, where rs13422013 (LFB=4.05) was also associated with ADAM23 gene expression (thalamus p=7.60×10⁻⁰⁴) and NRP2 (cerebellum p=4.40×10⁻⁰⁴). MR analyses indicated that ADAM23 (p=0.05) and not NRP2 (p>0.05) was the gene driving the association in this locus. Human brain single nuclei RNA-seq data indicate that ADAM23 expression is enriched in neurons (p<2.20×10⁻¹⁶); compared to all the other brain cell types. More specifically, its expression is enriched in excitatory neurons (Figure 2E).

Gene-based analyses using FUMA revealed that GRIA1 located in 5q33.2 was the gene most likely driving the association in that region (p=0.03, Z=1.83). However, Braineac identified several eQTLs for GALNT10 (p=3.60×10⁻⁰⁴) in the occipital cortex, but GALNT10 did not reach statistical significance in the gene-based analysis. GTEx portal and the protein atlas reveals that GRIA1 is mainly expressed in brain tissue. While GALNT10 is also expressed in the brain, it has higher expression in other tissues. The human brain single nuclei RNA-seq data confirmed that both GRIA1 and GALNT10 are expressed in divergent brain cell types (Figures 2H and Supplementary 4A). GRIA1 is highly-expressed in neurons compared to other cell types (p<2.20×10⁻¹⁶), but not expressed in oligodendrocytes (p<2.20×10⁻¹⁶) or astrocytes (p<2.20×10⁻¹⁶). In contrast, GALNT10 is expressed in microglia, oligodendrocytes and astrocytes, but expression in neurons is low (p<2.20×10⁻¹⁶). GRIA1 expression was also nominally associated with ΔNIHSS in the CLEAR trial dataset (p=0.002, r2=0.22).

Of the remaining five loci, we were able to map four (Supplementary Results). Briefly, eQTL analysis, revealed that 2q31.2 was likely to be driven by DFNB59, and 4q34.3 by MGC45800. No eQTLs were identified in 6q26, however, the locus falls within the boundaries of PARK2. Finally, 7p21.1 contains several eQTLs for TWISTNB and ABCB5 (Supplementary Results).

Pathway analyses

Gene ontology and pathway analyses using DEPICT and summary statistics for ΔNIHSS revealed consistent suggestive associations with functions relating to the brain and central nervous system. The top tissue enrichment from DEPICT identified the central nervous system (p=4.9×10⁻⁰³), including the brain (p=4.9×10⁻⁰³) and some brain regions: occipital lobe (p=1.52×10⁻⁰³), parietal lobe (p=5.31×10⁻⁰³), and hippocampus (p=5.07×10⁻⁰³; Supplementary Table 6). The most significant pathways in the gene-set enrichment were the AKAP5 protein-protein interaction subnetwork (p=5.11×10⁻⁰⁴), the CAMK1D protein-protein interaction subnetwork (p=5.27×10⁻⁰⁴), the RALA protein-protein interaction subnetwork (p=6.32×10⁻⁰⁴) and neurotransmitter uptake (p=6.74×10⁻⁰⁴; Supplementary Table 7). Several genome-wide significant candidate genes fell within these networks, including AGPS (2q31.2, LBF=5.57) in the RALA subnetwork and GPRS1 (2q33.3, LBF=6.34) in the CAMK1D subnetwork. AKAP5, expressed almost uniquely in brain (GTEx portal and the protein atlas.), interacts with post-synaptic density protein 95 (PSD-95) in dendritic spines and modulates synaptic transmission.^45-47 Both CAMK1D and RALA, expressed in most tissues (CAMK1D demonstrates higher expression in brain), are involved in signal transduction functions. MAGMA gene-set analyses revealed the gene-set “meissner_brain_hcp_with_h3k4me3_and_h3k27me3” (p=9.41×10⁻⁰⁶, Supplementary Table 8). Unfortunately, there is no function associated with this gene set that has been described in animal models. Together these results indicate that our unbiased genetic screening has identified genes that are largely expressed in the brain and in some cases involved in pathways associated with neuronal and synaptic function.

Unique genetic architecture of early outcomes after stroke

We examined the genetic architecture of ΔNIHSS for shared genetic variation and variance with other cardiovascular and aging-related traits, including stroke risk, age at death, plasma lipid levels and body mass index using PRSice (Supplementary Table 9 and Figure 5), LDSC (Supplementary Table 10) and GNOVA (Supplementary Table 11). Supported by a previous report¹⁸, PRSice did not reveal any genetic overlap between the genetic architecture of ΔNIHSS and stroke risk. Similarly, no overlap with age at death, lipid levels, or BMI was found. LDSC was unable to calculate the heritability estimate for ΔNIHSS. GNOVA, was successful at estimating the heritability for ΔNIHSS, and revealed significant covariance with stroke risk, but only when correcting for sample overlap (p=4.67×10⁻⁰⁵, correlation=0.80) and age at death (p=9.26×10⁻⁰⁴, correlation=1). However, the heritability estimate for ΔNIHSS for overlap with lipid levels and body mass index was negative, likely due to the low number of variants included in the analyses. Because both GNOVA and LDSC require larger sample sizes, the results of these analyses were inconclusive.