Wen Canhong, LI Jiarui, Cai Jiahui, Tan Haizhu*
1. School of Data Science, University of Science and Technology of China, Hefei 230027, China;2. Department of Statistics and Finance, School of Management, University of Science and Technology of China, Hefei 230026, China;3. Department of Preventive Medicine, Shantou University Medical College, Shantou 515031, China
Abstract: Alzheimer’s disease(AD) is a common neurodegenerative disorder with high heritability. An increasing number of common variants have been found to be associated with AD, but these common variants can only explain a small proportion of the heritability. Theory and practice have shown that rare variants can explain the remaining heritability. We explored rare functional variants that altering susceptibility to AD among 600470 variants in 389 individuals (175 with AD and 214 with cognitively normal). Firstly, after imputing the missing genotypes on the Michigan imputation server, quality control and gene-based annotation were carried out. Secondly, the efficient resampling sequence kernel association test was performed on 311 annotated exonic variants. Finally, the underlying biological interpretations of the identified risk gene were predicted through several bioinformatics tools. The results showed that under the Bonferroni correction, the rare missense variant rs2275303 in SIPA1L2 gene was significantly associated with AD (P=6.00E-04), and its pathogenicity was verified by bioinformatics analysis. SIPA1L2 gene is expected to play an important role in the prevention, diagnosis, prognosis and treatment of AD.
Keywords: Alzheimer’s disease; rare variants; association analysis; ER-SKAT
Alzheimer’s disease (AD) is a chronic, incurable, common neurodegenerative disorder with high heritability (estimated heritability up to 80%)[1,2]. Linkage analysis revealed that the apolipoprotein E (APOE) gene was the first confirmed susceptibility gene for AD[3]. However, APOE can account for at most 50% of the total genetic effect in AD[4]. A large number of genome-wide association studies (GWASs) have successfully identified other genomic risk factors for AD. For example, two extensive GWAS, of 12000 probable AD cases and 18000 age-matched non-demented controls, revealed three new candidates for the genetic risk of late-onset or sporadic AD: CLU, CR1 and PICALM[5]. A staged association study and testing of suggestive loci revealed variants in ABCA7, MS4A6A/MS4A4E, EPHA1, CD33 and CD2AP associated with AD[6]. Despite these discoveries, these risk common variants (minor allele frequency, MAF>0.05) identified by GWAS account for less than half of the heritability to date[7,8]. The reasons might be the common variant hypothesis for GWAS, that allelic variants exist in more than 1% to 5% of the population, which is only partly implicated in usual diseases[9,10].
Hence, there are many explanations for the remaining unexplained heritability (also called missing heritability), such as rare variants (0.005 How to detect rare variants in association studies is a challenge. The statistical power of single-variant association study is low for assessing the role of rare and low-frequency variants unless the sample size is huge. Multiple comparison adjustments are also required[16]. The main multiple-comparison problem is that type I error increases with each additional test. Hence, numerous statistical methods specifically designed to increase the power for the rare variant association analysis have been proposed. For example, the burden tests method assumes that all rare variants in the target region have effects on the phenotype in the same direction with the similar magnitude. The Sequence Kernel Association Test (SKAT), an extension of the C-alpha test, provides a robust test that is particularly powerful in the presence of protective and deleterious variants and null variants[17]. The Efficient Resampling-SKAT (ER-SKAT), calibrating single and gene-based rare variant association analyses, can improve computational efficiency over conventional resampling with low minor allele count (MAC) variant sets and can ameliorate the conservativeness of results in case-control studies[18]. To identify the risk rare variants related to AD, we explored rare functional variants altering susceptibility to AD in 389 individuals with 600470 variants (175 individuals with AD and 214 cognitively normal(CN) controls) from the Alzheimer’s disease neuroimaging initiative (ADNI) database (www.loni.ucla.edu/ADNI). After imputing missing genotypes on the Michigan imputation server, we performed quality control (QC) and gene-based annotation. Subsequently, we performed a gene-based association study of 311 selected annotated exonic variants in 254 genes by using ER-SKAT. The identified significant rare variant could explain some proportions of genetic heritability of AD. Data used in this study were downloaded from the ADNI-1 database. In total, 819 individuals were enrolled at the baseline, including 175 with AD, 398 with mild cognitive impairment (MCI), and 229 CN controls. The present study involved 175 individuals with AD and 214 CN controls with 600470 variants located on chromosomes 1-22 in PLINK data format. Information for the participants, including selection criteria, demographic data, and baseline assessments, are available on the ADNI website. Genotype imputation can increase genome coverage and improve power in GWAS. It greatly boosts the number of single nucleotide polymorphisms (SNPs) tested in association studies and helps in discovering related loci[19]. Moreover, rare GWAS variants to be imputed are easier to be tagged than are common GWAS variants[20]and could increase the imputation quality for rare variants. In our study, pre-imputation QC was conducted before phasing and imputation. We eliminated SNPs with genotype call rate>10%, 0.05 Although the human exome represents less than 2% of the genome, it contains about 85% of the known disease-related variants[22]. The exon-transcribed sequence remains in mature RNA after the intron is removed by RNA splicing and can be expressed as a protein during protein biosynthesis[23]. To identify whether the variants in the exonic or splice region might cause protein coding changes and the amino acids that are affected by the mutations[24], we used ANNOVAR (http://annovar.openbioinformatics.org/en/latest/) for gene-based annotation. Therefore, only exonic and splicing (variant is within 2 bp of a splicing junction) variants were chosen from the results of ANNOVAR. The gene-based association study based on the selected exonic rare variants involved using the SKAT R-package v1.3.2. In this procedure, small-sample-adjusted SKAT null model was used for a binary phenotype correcting for the gender, APOE4 allele and age. Finally, significant SNP sets were obtained after controlling by family-wise error rate (FWER) = 0.05. We used the Minimum Achievablep-values (MAP), the lower limit ofp-values for each variant set, to estimate the effective number of tests. MAP-adjusted and unadjusted quantile-quantile (Q-Q) plots were drawn. Several prediction tools were used to show the protein structure and verify the underlying biological interpretations of the identified gene. These bioinformatics databases and resources include NCBI (https://www.ncbi.nlm.nih.gov/gene/), genecards (https://www.genecards.org/), targetvalidation (https://www.targetvalidation.org/), varsome (https://varsome.com/), genepine (http://grch37.genepipe.ncgm.sinica.edu.tw/variowatch/), Noncode (http://www.noncode.org/), hosphosite (https://www.phosphosite.org/homeAction),gnomad (https://gnomad.broadinstitute.org/), UCSC Genome Browser on Human (https://genome.ucsc.edu/), and Online Mendelian Inheritance in man (OMIM) (https://omim.org/entry/). Table 1 describes the characteristics of the AD and CN groups. 93 (53.1%) participants were male in the AD group, and 115 (54%) in the CN group. The average age in AD and CN groups was 75.4 (SD =7.4) and 75.7 (SD=4.9). Approximately 67% and 27% of the AD and CN groups carried at least one APOE4 allele. Table 1. Demographics and clinical indicators of participants with Alzheimer’s disease (AD) and cognitively normal (CN) controls. After pre-imputation QC, 8,121 variants and 389 individuals were retained. We chose an Rsq threshold of 0.8 for SNPs with 0.005 Figure 1. The x-axis is the MAP-adjusted or unadjusted expected quantile of -log10 p-values, and the y-axis is quantiles of -log10 p-values. Observed p-values are plotted against the MAP-adjusted expected quantiles (black dots) and unadjusted expected quantiles (gray dots). The dashed line represents a 95% confidence interval based on 1000 random draws from the MAP-based mixture distribution. The protein primary structure of the SIPA1L2 gene is shown in Figure 2. The accession number for SIPA1l2 protein is NP_065859.3 (for human SIPA1l2). Structural analysis revealed a conserved RapGAP domain followed by a PDZ signaling domain. The PDZ domain may be responsible for specific protein-protein interactions, as most PDZ domains bind C-terminal polypeptides and bind internals (non-C-terminal). SPAR_C region is the C-terminal domain of SPAR protein. Figure 2. The protein primary structure of the SIPA1L2 gene. We also found that the identified variant in SIPA1L2 is a missense mutation. The Gene Ontology annotations related to SIPA1L2 included GTPase activator activity. SIPA1L2 is suggested to be most abundant in granule cells of the dentate gyrus and cerebellum and show RapGAP activity for the small GTPases Rap1 and Rap2. The gene is involved in two related pathways in the Ras signaling pathway including the Ras and Rap1 signaling pathways (one KEGG pathway is the Rap1 signaling pathway, which belongs to Ras signaling pathway). In mammalian models, some evidence has implicated components of the Ras signaling pathway in the aging and metabolic regulation[25]. Cellular senescence is a feature of age-associated diseases, such as AD and osteoarthritis[26]. Oncogenic mutations in Ras result in its hyperactivation to trigger senescence[27]. From OMIM, SIPA1L2, a 1514-amino acid protein, shares significant similarity with SIPA1L1, and SIPA1L2 is moderately expressed in all adult and fetal tissues and specific adult brain regions. Target validation showed that on microarray analysis of six brain areas from AD patients and normal individuals, RNA expression was decreased in the posterior cingulate cortex of AD patients. Alias for SIPA1L2 Gene is Spine-associated RapGAP 2 (SPAR2). SPAR2 is a novel GTPase activating protein (GAP) for the small GTPase Rap that shows significant sequence homology to SPAR.SPAR is a synaptic RapGAP, which has been reported to regulate the spine morphology in hippocampal neurons[28]. Overexpression of Rpph1 increased the density of dendritic spines in primary cultured hippocampal pyramidal neurons, while the knockdown of Rpph1 had the reverse effect, which may represent a compensatory mechanism at the early stage of the AD pathogenesis[29]. Like SPAR, SPAR2 interacts with the recently described synaptic scaffolding protein ProSAP interacting Protein (ProSAPiP), which in turn binds to the PDZ domain of ProSAP/Shank post-synaptic density protein. SPAR2 transcripts were mainly expressed in cerebellar and hippocampal granule cells. In a recent study, SPAR2 was identified as a brain-derived neurotrophic factor (BDNF) responder gene that is upregulated during cerebellar granulosa cell development[30]. In neurons, BDNF is transported in amphisomes, which signal locally at presynaptic boutons during retrograde transport to the soma. In hippocampal neurons, TrkB-signaling endosomes are actually amphisomes with local signaling capacity in the context of presynaptic plasticity during retrograde transport. The autophagosomal protein LC3 regulates RapGAP activity of SIPA1L2 and controls the retrograde transport and local signaling of TrkB[31]. Hence, SIPA1L2 is a potential rare risk variant for AD. After a series of bioinformatics analysis, the functions of SIPA1L2 gene in biomedicine were obtained, see Table 2. Table 2. Functions of SIPA1L2 gene. Although genotype imputation better enables the imputation of genotypes of rare variants and low-frequency variants from existing GWAS datasets, we analyzed rare variants with only 0.005 An increasing number of rare-variant association methods have been built depending on the underlying genetic architectures of complex traits. We used the ER-SKAT method to analyze the rare-variant association with AD based on integrating genomic and demographic data in our study. ER-SKAT aggregates variants within a specified region without considering the direction of the effect for individual variants, and it successfully solves the problem of linkage disequilibrium between genetic variations. Therefore, ER-SKAT is more suitable to detect associations in which both risk and protective variants or numerous non-causal variants are present. Lee et al. proposed an efficient resampling method based on ER-SKAT to calibrate single and gene-based rare-variant association analyses in case-control studies. This method can improve computational efficiency >1000-fold and significantly shorten the calculation time over conventional resampling for low minor-allele-count variant sets. It ameliorates the conservativeness of results by using the mid-p-value and the estimated minimum achievablep-value for each test. By applying this method, we obtained the calibrated QQ plots and the number of effective tests to judge the significance and identified the significant gene SIPA1L2. However, our research has some limitations. First, the identified gene has not been verified with subsequent molecular genetics experiments. Second, the number of samples in the ADNI-1 dataset is limited for the rare-variant association research. In conclusion, we found that SIPA1L2 is significantly associated with AD. Our research has demonstrated that rare variants may disrupt the gene function of AD as well. Hence, determining the gene function of rare variants with full-genome resequencing is a promising area for the AD genomic research. Availability of data and material The ADNI was initiated in 2003 by the U.S. National Institute on Aging, the U.S. National Institute of Biomedical Imaging and Bioengineering, the U.S. Food and Drug Administration (FDA), private pharmaceutical companies, and nonprofit organizations as a public-private partnership, led by the principal investigator Michael W. Weiner, MD. The primary goal of ADNI is to test whether serial MRI, positron emission tomography, other biological markers, and clinical and neuropsychological assessment can be combined to measure the progression of the mild cognitive impairment and early Alzheimer’s diseases. Code availability Our codes and software lists are available for free at https://github.com/ruijiali/SIPA1L2. Acknowledgments The work is supported by the National Natural Science Foundation of China (11801540), the Natural Science Foundation of Anhui (BJ2040170017), the Fundamental Research Funds for the Central Universities (WK2040000016) and the Science and Technology Program of Guangzhou (202002030129); the Science and Technology Planning Project of Guangdong Province (2017A010101030),the National Key Research and Development Program of China (2018YFC1315400) and the Third Medical Technology Projects of Shantou City in 2018 (ShanFuKe【2018】128-08). Conflictofinterest The authors declare no conflict of interest. Authorinformation WenCanhongreceived the BS degree in 2009 from Hefei University of Technology and the PhD degree in 2015 from Sun Yat-sen University. Currently, she is an associate professor in the School of Data Science, and the Department of Statistics and Finance in the School of Management at the University of Science and Technology of China. Her main research interests include statistical learning, biostatistics, and image genetics. TanHaizhu(correspongding author) received the BS degree in 1999 from Shantou University of Computer Technology and the PhD degree in 2014 from Sun Yat-sen University. Currently, she is an associate professor in the Department of Preventive Medicine, Shantou University Medical College,China. Her main research interests include biostatistics, bioinformatics and image genetics.2 Materials and methods
2.1 Participants
2.2 Statistical analysis
2.3 Bioinformatics analysis
3 Results
3.1 Basic demographic features

3.2 Significant gene identified by ER-SKAT

3.3 Protein structure and functions of SIPA1L2 derived from bioinformatics analysis


4 Discussion
5 Declarations