The large-scale MS GWAS data set consisted of 9,772 MS cases from the IMSGC and 17,376 controls from the WTCCC2. And cases were gathered by 23 teams operating in 15 countries. All individuals were self-reported as being of European descent. The original research conducted sample and single-nucleotide polymorphism (SNP) QC to create SNP metric strainers. In sample QC, Bayesian clustering and principal components analysis were mainly used to identify qualified samples. Researchers developed and utilized two novel methods (automated cluster checking and beta-binomial model) in SNP QC to exclude SNPs identified because of incorrect genotype calling. Finally, 464,357 autosomal SNPs from the entire data set were obtained for further analysis (International Multiple Sclerosis Genetics Consortium et al., 2011).

The IS data set was summarized from the METASTROKE collaboration’s discovery stage. In this stage, genotyped data from 12 case–control GWASs correspond to 10,307 IS cases and 19,326 controls (ASGC: 1,162 IS cases and 1,244 controls; BRAINS: 371 IS cases and 2,640 controls; GASROS_affy: 485 IS cases and 3,030 controls; GASROS_illumina: 296 IS cases and 377 controls; GEOS: 448 IS cases and 498 controls; HPS: 588 IS cases and 571 controls; ISGS-SWISS: 1,014 IS cases and 1,370 controls; MILANO: 366 IS cases and 407 controls; VISP: 1,723 IS cases and 1,047 controls; WHI: 306 IS cases and 2,170 controls; WTCCC2-D: 1,174 IS cases and 797 controls; WTCCC2-UK: 2,374 IS cases and 5,175 controls; all individuals in the discovery stage were Caucasian). After QC by logistic regression analysis and meta-analysis, 9,541,572 SNPs were available for gene-based testing (Malik et al., 2016).

We utilized the VEGAS2 to execute more elastic gene-based testing for MS and IS GWAS. The chief distinguishing feature of VEGAS2 is that this approach is computationally feasible and can be applied to any GWAS experimental design by using the summary GWAS data (Mishra and Macgregor, 2015). After uploading the SNP information (rs-number and P value of both GWASs) to VEGAS2, we selected 1000G EUROPEAN as the population reference group. For considering both the physical position and LD profiles of SNPs with others in gene-based testing, we chose “0kbldbin outside gene and SNPs in LD r² > 0.8” (Christoforou et al., 2012). For each gene definition, P values of SNPs were first changed to upper tail χ² statistics with one degree of freedom (df) in the gene definition. And then, if SNPs are in linkage equilibrium, the statistics of gene-based tests would have a χ² distribution with n df under the null hypothesis. We used the reference population group of 1000 Genomes European and the model of Σ (n × n matrix of LD [r] values) for SNP correlation, because LD for the n SNPs occurs frequently. Significance was tested with reference to the aggregate χ² statistics for each gene to simulated repetition from a multivariate normal distribution in which mean = 0 and variance = Σ. The empirical P value of every gene was calculated via the formula, P = r + 1/m + 1, where r was the number of cases where the simulation statistics surpassed the observed statistics, and m was the number of simulations (Mishra and Macgregor, 2015). In this gene-based testing, we screened the shared genes with nominal significance respectively (P-_MS < 0.05, P-_IS < 0.05).

To combine P values calculated by VEGAS2 for MS and IS GWAS, we used Fisher’s method to perform meta-analysis for every overlapped gene. For a given gene, the formula for the statistic is as follows:

For the i_th study, P_i is the P value of the genes, and k is the entire count of studies. x² abides by a c² distribution with 2k degrees of freedom (Begum et al., 2012). The gene-based meta-analysis was performed with R software.¹

Following the meta-analysis, we performed Bonferroni correction to control type I error. The combined P value was less than 0.05/2n, where n was the number of shared genes with previous nominal P value (P < 0.05).

We further investigated the differential expression of shared genes in the MS and IS patients and healthy control subjects using gene expression data sets from the NCBI GEO database.²

Immune cells from whole blood are the typical samples for gene expression analysis in patients with MS. PBMCs from 12 MS female patients under Poser’s criteria and 15 unrelated female controls were isolated from whole blood. The extracted samples were tested on the Affymetrix Gene Chip Human Genome U133 Plus 2.0 Array (Kemppinen et al., 2011). The expression data were obtained from GEO series GSE21942. In GEO series GSE43591, samples were extracted from 10 relapsing–remitting MS patients diagnosed in accordance with the criteria of McDonald et al. (2001) and 10 age- and sex-matched controls. T cells were purified through CD14⁺ removal sorting and CD3⁺-positive selection from PBMCs. Transcriptional data were examined with Human Genome HG-U133 plus 2.0 arrays (Jernas et al., 2013).

Ischemic stroke case–control expression analyses were performed to identify differentially expressed genes in the GEO series. GEO series GSE16561 provided the raw mRNA expression data from peripheral whole blood of a study of 39 IS patients and 24 healthy controls (O’Connell et al., 2017). More information about the inclusion and exclusion criteria was provided in the original study (O’Connell et al., 2017). RNA data from the blood of 23 control samples and 69 cardioembolic stroke samples were analyzed in GEO series GSE58294. Cardioembolic stroke is a subtype of IS and has unique genetic traits (Stamova et al., 2014). In the original research, 69 cardioembolic stroke samples were collected from 23 patients at three time points to monitor transcriptional changes in the first 24 h. Twenty-three subjects without symptomatic vascular disease history were recruited as VRFCs (Stamova et al., 2014). And we used all the 69 cardioembolic stroke samples as a whole (including the samples of three time points) to compare with the samples from VRFCs.

To calculate differential expression between case–control samples, we utilized GEO2R and performed GEO query and limma R packages (Barrett et al., 2013). Among four differential expression data sets, those with the smallest P value for each shared gene in MS and IS were selected for further analysis.

We uploaded the 464,357 MS and 9,541,572 IS SNP rs IDs and associated P values for gene-based analysis. We sorted 1,353 genes with P-_MS < 0.05 from a total of 14,811 MS gene sets obtained by VEGAS2. After Bonferroni correction (P < 0.05/14,811 = 3.38E-06), 47 non-major histocompatibility complex (MHC) and MHC risk variants were sorted including known risk genes such as AHI1, RGS1, SP140, IQCB, IL2RA, TNFRSF1A, CLEC16A, GALC, TNFSF14, and CYP24A1. The good applicability of VEGAS2 for gene-based test was proven by the well-replicated genes above. Information of the gene-wide significant genes of MS is shown in Supplementary Table S1. In IS data sets, 1,290 genes with P-_IS < 0.05 were screened out from 21,913 IS genes via VEGAS2. However, no IS variant achieved the strict threshold (P-_IS < 0.05/21,913 = 2.28E-06). And we utilized a looser significance threshold (P < 1.00E-04) and identified two genes that satisfied this threshold set: ZYX (P < 2.60E-05) and NCR3LG1 (P < 8.10E-05).

After making an intersection between 1,353 MS genes and 1,290 IS genes (P-_MS < 0.05, P-_IS < 0.05), we obtained 108 nominal shared genes. We performed meta-analysis for the 108 shared genes and obtained the P_–combined for each gene. And then 24 overlapped significant genes were identified with P_–combined < 2.31E-04 (0.05/2^∗108), including LBH at 2p23.1, STAT4 at 2q32.2-q32.3, FOXP1 at 3p13, TEC at 4p12-p11, MICB/NFKBIL1 at 6p21.33, IL6 at 7p15.3, CAMK2G at 10q22.2, B4GALNT1 at 12q13.3, BRAP at 12q24.12, CLEC2D/CLECL1 at 12p13.31, HECTD4/NAA25 at 12q24.13, OS9 at 12q13.3-q14.1, NUDT14 at 14q32.33, PDE4C at 19p13.11, SLC44A2 at 19p13.2, and LIME1/RTEL1/RTEL1-TNFRSF6B/SLC2A4RG/ZBTB46/ZGPAT at 20q13.33. The detailed results are listed in Table 1.

We analyzed and integrated GEO profiles of MS and IS to investigate the differential expression levels of shared genes. Notably, 16 of 24 common genes were determined to have significant changes in at least one of the four GEO data sets with P < 2.08E-03 (0.05/24) (Table 2).

In MS expression profiles, the transcriptional levels of CAMK2G (P = 4.15E-04), CLEC2D (P = 9.11E-04), FOXP1 (P = 2.35E-09), and NUDT14 (P = 1.51E-04) in the GEO data set GSE21942 (Supplementary Table S2) and CAMK2G (P = 3.09E-06), LBH (P = 5.60E-04), SLC2A4RG (P = 1.29E-03), and ZGPAT (P = 6.93E-04) in MS patients were markedly changed compared to controls (Supplementary Table S3). In samples from peripheral whole blood, CLEC2D (P = 2.68E-07), LBH (P = 5.45E-06), LIME1 (P = 2.33E-07), OS9 (P = 1.41E-04), SLC2A4RG (P = 5.07E-05), and SLC44A2 (P = 4.73E-04) in the acute IS case–control set (GSE16561) (Supplementary Table S4) and BRAP (P = 1.71E-20), CAMK2G (P = 5.14E-07), CLECL1 (P = 6.30E-06), FOXP1 (P = 1.55E-05), HECTD4 (P = 8.83E-05), LIME1 (P = 1.16E-09), MICB (P = 2.45E-04), NAA25 (P = 3.92E-08), NFKBIL1 (P = 3.33E-06), OS9 (P = 1.24E-12), and SLC44A2 (P = 1.27E-03) in the cardioembolic stroke case–control set (GSE58294) (Supplementary Table S5) showed significantly different expression.

It was remarkable that there were five overlapped genes with a significant expressed difference in at least one data set of each disease GEO data sets (P < 2.08E-03) (Table 3). The expression difference of every individual gene presented the same tendency between two GEO data sets in MS, and CAMK2G was also upregulated in IS GEO data sets, whereas three of five genes (CLEC2D, LBH, and SLC2A4RG) had opposite alterations in GSE16561 and GSE58294, which might be due to the heterogeneity of IS and its subtypes. FOXP1 expression in patients was remarkably increased in every data set compared to controls [GSE21942: log fold change (FC) = 1.21, GSE43591: logFC = 0.528, GSE16561: logFC = 0.15, and GSE58294: logFC = 0.34].

Forkhead box protein P1 protein is a member of the forkhead box (FOX) family. The forkhead TF plays vital roles in diverse cell and tissue processes including development, aging, metabolism, and cancer (Coffer and Burgering, 2004; Wijchers et al., 2006). FOXP1 is also critical in the development and function of immune cells. In peritonitis models, Shi et al. (2008) illustrated that FOXP1 upregulation blunted monocyte development and macrophage biological activities by inhibiting the production of c-Fms/macrophage colony-stimulating factor receptor. Via Smad2/Smad3 and transforming growth factor β (TGF-β) signaling, increased FOXP1 expression restrained CD8⁺ T cells from proliferation and activation in cancer (Stephen et al., 2014). Combined with FOXP3, FOXP1 was demonstrated to play an essential role in maintaining expression of CD25, SATB1, and CTLA-4 and responsiveness to interleukin 2 in regulatory T cells (Konopacki et al., 2019). FOXP1 was also reported as a crucial negative transcriptional modulator in the differentiation of CD4⁺ follicular helper T cells (Wang et al., 2014). Similarly, FOXP1 is a fundamental TF in the early development of B cells (Hu et al., 2006), and abnormal FOXP1 upregulation leads to a reduction and irregular distribution of B cells in germinal centers (Sagardoy et al., 2013). FOXP1 is associated with both CNS development and abnormalities. Braccioli et al. (2017) reported that FOXP1 can adjust neural stem cells (NSCs) neurogenesis via Notch signaling and foster the differentiation of embryonic NSCs to neurons and astrocytes in vitro. Li et al. (2015) demonstrated that FOXP1 has vital influences on modulating neuronal migration and morphogenesis in cortical regions. FOXP1 mutations reportedly contribute to nervous system disorders including Huntington disease (Tang et al., 2012), autism (O’Roak et al., 2011; Chien et al., 2013), and epilepsy (Jay et al., 2019).

In this study, FOXP1 exhibited an elevated tendency of expression in both MS and IS patients compared to controls (GSE21942: logFC = 1.21; GSE43591: logFC = 0.528; GSE16561: logFC = 0.15, and GSE58294: logFC = 0.34). FOXP1 was significantly more than twofold upregulated in the PBMC of MS (P = 2.35E-09, logFC = 1.21) and significantly 1.3-fold upregulated in the blood cells of cardioembolic stroke (p = 1.55E-05, logFC = 0.34). The IMSGC assessed two sets of GWAS data and identified rs9828629/FOXP1 (P = 1.9E-10, odds ratio [OR] = 1.08) as an immune-related risk variant linked to MS (International Multiple Sclerosis Genetics Consortium et al., 2013). In research on atherosclerosis, Bot et al. (2011) illustrated that FOXP1 was expressed in diverse cell types and was related with stable plaques through the TGF-β pathway. As atherosclerosis is one cause of large vessel disease, a major subtype of IS, there may be more associations between FOXP1 and IS. In summary, FOXP1 could be a potential target in terms of its multiple functions in the immune and nervous systems.

CAMK2G encodes CaMKIIγ, an isoform of calcium (Ca²⁺)/calmodulin-dependent protein kinase II (CaMKII) that participates in Ca²⁺-related biological activities (Hudmon et al., 2001; Rokita and Anderson, 2012). In a model of acute ischemia/reperfusion, CAMK2G/CaMKIIγ improved neuronal survival by activating nuclear factor κB signaling (Ye et al., 2019). CAMK2G was identified as an enhancer gene for coronary artery disease in a GWAS meta-analysis (Gong et al., 2018). CAMKIIγ was also a key modulator to repress macrophage phagocytosis and enhance the necrosis of atherosclerotic plaques (Doran et al., 2017). We found that CAMK2G (P = 5.14E-07) was 1.3-fold upregulated in cardioembolic stroke compared to controls. In addition, CAMK2G could act through classical inflammatory pathways such as mTORC1 and STAT3, both of which play vital roles in MS (Ma et al., 2017; Meng et al., 2017).

CLEC2D encodes lectin-like transcript 1 (LLT1) and is located next to the NK gene complex (Yokoyama and Plougastel, 2003; Germain et al., 2010). LLT1 was identified as a negative ligand for NKR-P1A (CD161) receptor in humans (Rosen et al., 2005). The interaction of CD161 on NK cells with matched LLT1 can reportedly repress NK cell-mediated cytotoxicity (Aldemir et al., 2005). Llibre et al. (2016) illustrated that LLT1/CD161 interaction influences B-cell development in human germinal centers. Toll-like receptor (TLR)-activated dendritic cells and TLR or B-cell receptor-activated B cells express LLT1 (Rosen et al., 2008). LLT1 is also expressed by monocyte/macrophage in RA joints, and soluble LLT1 is elevated in patient serum (Chalan et al., 2015). Here, we detected that the expression of CLEC2D (P = 9.11E-04, logFC = 0.49) was significantly altered in PBMCs from MS patients. However, the expression of CLEC2D in IS cases was downregulated by 0.62-fold compared to controls. The distinct alterations of CLEC2D expression in MS and IS need further study.

SLC2A4 regulator encoded by SLC2A4RG was identified as TF binding to the SLC2A4 promoter district where it regulates expression (Oshel et al., 2000; Knight et al., 2003). In a study of Huntington disease, SLC2A4RG was found to cooperate with a key cis-element and shuttle to and from the nucleus (Tanaka et al., 2004). A recent large GWAS demonstrated that rs2256814/SLC2A4RG (P = 3.5E-9, OR = 1.08) was a novel MS susceptibility gene with immune function (International Multiple Sclerosis Genetics Consortium et al., 2013). Furthermore, Dhaouadi et al. (2014) stated that SLC2A4RG could slightly induce TGFB1 expression in atherosclerosis. Notably, SLC2A4RG showed differential expression in MS and IS GEO data sets (GSE43591, P = 1.29E-03; GSE16561, P = 5.07E-05).

Limb-bud and heart (LBH) encodes a transcription modulator that regulates cell development in multiple tissues (Ai et al., 2008; Al-Ali et al., 2010; Ekwall et al., 2015). LBH and related variants are reportedly associated with autoimmune diseases such as SLE, RA, and celiac disease (Zhernakova et al., 2011; Yu et al., 2013; Chang et al., 2016). Wnt signaling regulates the expression of LBH in tumors and the epithelium (Rieger et al., 2010). In the CNS, the Wnt pathway is related with male-specific genes of IS and positively impacts the regenerative process after injury (Tian et al., 2012; Lambert et al., 2016). Researchers integrated Alzheimer’s disease-GWAS data and identified LBH as a pathogenic gene of amyloid β accumulation that is tightly linked to immune system (Mhatre et al., 2015; Yamaguchi-Kabata et al., 2018). In our study, LBH was significantly downregulated both in T cells from PBMCs of patients with MS (P = 5.60E-04, logFC = −0.55) and blood cells from patients with IS (P = 5.45E-06, logFC = −0.55). The potential immunomodulatory function of LBH can assist us in uncovering common mechanisms between MS and IS.

Although we conducted powerful gene-based tests that are effective extensions of traditional GWASs to discover shared genes, our study has several limitations. First, subthreshold variants will undermine the significance of the causal SNPs in the same gene region. Next, IS consists of several subtypes with genetic heterogeneity. More refined data of IS subtypes should be excavated for further analyses. Furthermore, the single sample type and limited case numbers are two major disadvantages of our GEO data. All GEO sample sources are from peripheral blood cells because they are safe and accessible in clinical settings. If possible, we will integrate brain-derived data that can directly reflect the differential expression of genes in neurological diseases.