Exploring SVA Insertion Polymorphisms in Shaping Differential Gene Expression in the Central Nervous System

1. Introduction

Amyotrophic lateral sclerosis (ALS) is a fatal neurodegenerative disease, originally termed by the French neurologist Jean-Martin Charcot in 1869 to describe muscular atrophy (amyotrophic) and tissue scarring and hardening of tissue within the lateral spinal cord [1,2]. ALS is the most common form of motor neuron disease (MND) and is characterized by the progressive deterioration of both upper and lower motor neurons within the brain and spinal cord [2,3]. The incidence of which is 1.75 – 3 per 100 000 people per year and an increased incidence of 4 – 8 per 100 000 people per year in the highest ALS risk age group (45 – 75 years old) [4]. Categorised as 2 groups, ALS can present as familial ALS where at least one family member of the affected individual has ALS, accounting for up to 10% of ALS cases, or sporadic ALS (sALS) where the affected individual has no prior family history which accounts for 90 – 95% of cases [5]. The clinical features of typical ALS patients consist of muscle spasticity, atrophy, muscle wasting, weakness and death due to respiratory failure, with the average survival after symptom onset lying between 3-5 years [1,4,6].

Neurodegenerative diseases are complex disorders, involving both environmental and genetic factor interactions. Therefore, attempts to elucidate disease mechanisms and pathogenic genetic variants are essential. Previous research has identified more than 30 genes associated with ALS, highlighting four genes SOD1, TARDBP (TDP-43), C90RF72 and FUS for harbouring pathogenic mutations which cause the greatest number of ALS cases [7,8,9,10]. Although these four genes have been identified as major ALS-associated genes, a 2017 meta-analysis study demonstrated that within European and Asian populations, these genes only contribute to 47.7% and 5.2% of familial and sporadic cases, respectively [11]. Following the identification of these pathogenic mutations, several pathological mechanisms have been implicated in ALS pathogenesis including oxidative stress, mitochondrial dysfunction, axonal transport, inflammation, toxic protein aggregation and RNA metabolism and toxicity [7]. In addition to identified genetic variants only explaining a small fraction of sporadic aetiology, twin studies have highlighted the importance of genetic risk factors within sporadic ALS, estimating a 61% heritability [12]. Therefore, as the exact causation of ALS is still undetermined, a better understanding of disease pathogenesis and identification of genetic biomarkers is essential. To further investigate the missing heritability of ALS, Theunissen et al., proposed structural variations (SVs) as an area of potential significance [13]. SVs are classified as insertions, inversions, deletions and microsatellites, usually of repetitive structure, that are predominantly present within non-coding DNA regions and contribute towards genomic variation [13,14]. Furthermore, SVs have demonstrated the ability to modulate gene expression and have already been implicated in neurodegenerative diseases, such as the C9orf72 repeat expansion in ALS and frontotemporal dementia (FTD) [15,16]. Hence, as 99% of the human genome is non-coding DNA, continued research into these regions is crucial to provide new insights into disease pathogenesis and to identify new potential targets for therapeutics.

Repetitive DNA is a major contributor to structural genomic variation. One form of repetitive DNA is a group of endogenous transposable elements (TEs), which can exist in both a static form and a mobile form. TEs are categorised into two classes known as DNA transposons and retrotransposons, whereby the later class possesses the ability to propagate throughout the genome via a ‘copy-and-paste’ mechanism, involving an RNA intermediate [17]. This results in the insertion of a new retrotransposon copy at a new locus within the host genome. Although, DNA transposons are capable of mobilisation via a ‘cut-and-paste’ mechanism, they are not currently active within the human genome [8,18]. Originally dismissed as “junk” DNA, TEs are known to drive genetic diversity not only by contributing to regulation and evolution of the genome but also by contributing to genetic instability and disease progression [18,19]. Previous research by Prudencio et al., highlighted the potential implication of retrotransposons in ALS through the analysis of repetitive element expression using RNA sequencing data from both healthy controls as well as C9orf72 expansion positive carriers and sporadic ALS patients (C9orf72-negative) [20]. This research revealed that repetitive element expression, including retrotransposons, was significantly increased in ALS patients with the C9orf72 expansion in comparison to C9orf72-negative patients and healthy controls, thus, suggesting retrotransposon involvement in ALS [20]. However, until recently these elements have been largely overlooked in relation to neurodegenerative diseases, even though TEs constitute to around 45% of the human genome [21].

Retrotransposons are further subdivided into two groups dependent on the presence of long terminal repeats (LTRs), known as LTR and non-LTR retrotransposons [18]. SINE-VNTR-Alu (SVA) elements are a member of the non-LTR retrotransposon family, which are typically 0.7 – 4kb in length [22]. SVAs are classified by evolutionary age based on their SINE-R region into sub-family groups A-F, whereby subfamily SVA-F is the youngest in evolutionary history [23]. Full-length SVA elements contain a 5’ CT element, Alu-like region, GC-rich VNTR (variable number tandem repeat), SINE (short interspersed nuclear element)-R domain and a 3’ poly-A tail [23] (Figure 1b). SVAs are major contributors to genetic diversity through a variety of different mechanisms, including acting as transcriptional regulators by providing alternative splice sites, polyadenylation signals and promoters whilst harbouring sites for transcription factor (TF) binding, thus, modulating gene expression [24]. SVAs are not only polymorphic in structure but ongoing mobilisation has resulted in SVAs being polymorphic for presence or absence within the genome and thus are termed retrotransposon insertion polymorphisms (RIPs) [19]. This adds additional layers of complexity to gene expression dynamics, indicating that SVAs could be associated with predisposition to disease [25].

Previous work within our group has shown that seven SVAs polymorphic for their presence or absence were significantly associated with Parkinson’s disease progression and differential gene expression in Parkinson’s disease (PD) patients using the Parkinson’s Progression Markers Initiative (PPMI) cohort [22,26]. This involved an SVA named SVA_67 which is located at the MAPT locus, a locus which has been implicated in neurodegenerative disease risk including PD, FTD and Alzheimer’s disease [27,28,29]. By using a clustered regularly interspaced short palindromic repeats (CRISPR) cell line knock-out model for SVA_67, we demonstrated that this SVA was significantly associated with differential gene expression of three genes at the MAPT locus, associated with neurodegeneration [30].

In this study, we analysed SVAs which are present in the human reference genome identified as RIPs in the cohort analysed, herein termed reference SVA RIPs to investigate differential gene expression patterns associated with ALS. For this analysis, we used whole genome sequencing (WGS) and transcriptomic data obtained from the New York Genome Center (NYGC) ALS Consortium, to elucidate the role of SVA insertion polymorphisms on gene expression in central nervous system (CNS) tissues in healthy controls and ALS patients (Figure 1a). Our analysis demonstrated that reference SVA RIPs regulate single and multiple gene targets genome-wide within CNS tissues and display tissue-specific gene modulation by comparing spinal cord, cerebellum, motor cortex and frontal cortex tissue analyses. Furthermore, we discovered that SVA RIPs influence the expression of genes at loci (HLA and MAPT) previously associated with ALS, highlighting SVA regulation of genes at this locus could be a potential mechanism involved in ALS pathology.

2. Results

2.1. SVA RIPs act as eQTL genome wide in CNS tissues

To determine the impact of reference SVAs RIPs on gene expression profiles, we analysed WGS and transcriptomic data of both ALS patients (n=1544) and healthy controls (n=359) from the ALS Consortium cohort. For this analysis, we initially combined data from different CNS tissues, including motor cortex, frontal cortex, cerebellum, occipital cortex, temporal cortex, hippocampus, and spinal cord tissue (Figure 2).

For eQTL analysis we assessed changes at a gene-based level. We demonstrate that, 14,830 genomic loci were significantly differentially modulated by reference SVA insertion polymorphism (FDR p<0.05). Of these effects, 167 were cis regulated by SVA RIPs, whilst over 98% (14,663 targets) were trans-regulated, where trans-regulated is defined as genes regulated greater than 1Mb away from the SVA site (Figure 3a). This indicates that most regulatory effects observed by these SVA RIPs impact target genes >1Mb away from the SVA element itself, these could be both direct and indirect effects (Figure 1c). In addition, single SVA elements had the capacity to regulate either a single gene or multiple gene targets. Our analysis identified that 92 reference SVA RIPs were responsible for the modulation of the 14,830 genomic loci, with SVA_16 affecting the greatest number of genes (917) (Figure 3b). Of these 92 SVA elements polymorphic for their presence or absence, 26 SVAs affects on more than 200 genes in total, whilst the lowest number of targets for an SVA in this analysis was 6, indicating that all reference SVA RIPs analysed in this study had the ability to modulate expression profiles of multiple gene targets in the CNS (Figure 3b).

2.2. Mitochondrial genes are significantly modulated by SVA RIPs

To further investigate the capacity of reference SVA elements to act as eQTL, we analysed the beta values obtained from the eQTL analysis, to determine the magnitude of the measured effects induced by SVA RIPs. The ten reference SVAs RIPs with the greatest effect size on gene upregulation (highest beta value) and gene downregulation (lowest beta value) and the affected gene targets are displayed in Figure 4a,b. All top ten hits for the highest and lowest beta values are trans-regulatory effects. SVA_55 demonstrated the greatest increase in activation of the MBP gene with a beta coefficient 272,847, whilst SVA_70 demonstrated the greatest repressive effect of the MT-ND1 gene with a beta coefficient of -36441 (FDR p<0.05) (Figure 4a,b). Interestingly, it was noted that the only two SVAs (SVA_55 and SVA_15) were responsible for the top ten hits regarding the effect size for gene activation. Furthermore, although the top hit for effect size for gene repression was SVA_70 modulating MT-ND1 expression, our results also demonstrated that the presence of both SVA_15 and SVA_55 showed the opposite effect by upregulating MT-ND1 expression (Figure 4a,b). To determine the impact of SVA presence on gene modulation in more detail, we examined the top hits for SVAs with the greatest effect on gene upregulation (SVA_55) and downregulation (SVA_70) independently. For MBP, the presence of two copies/alleles (PP genotype) of SVA_55 significantly upregulated MBP gene expression compared to the presence of one copy of SVA_55 (PA genotype) or complete absence (AA genotype) (p=0.0218) (Figure 4c). In comparison to the AA and PA genotypes, individuals with the PP genotype displayed a 9.5-fold and 19.2-fold increase in MBP gene expression, respectively. In contrast, individuals homozygous present for SVA_70 demonstrated significant repression of the MT-ND1 gene in comparison to individuals heterozygous for SVA_70, namely down-regulating gene expression of MT-ND1 by 55% (p=0.0098) (Figure 4d). No statistical significance for differential gene expression between the PP and AA genotype was obtained for MT-ND1 expression using the Wilcoxon pairwise comparison with FDR adjusted p-values (FDR p<0.05).

We next analysed the top ten most significant SVA effects based on FDR p-value. The most significant effects observed (lowest FDR p-value) were SVA_67 acting on the genes MAPK8IP1P2 and ENSG00000285668.1, displaying an FDR p-value of 1.93E-303 (Table 1). Followed by SVA_67 acting on the gene LRRC37A, with an FDR p-value of 5.12E-299 (Table 1). Intriguingly, 6/10 hits for highest significance were cis-specific effects on gene upregulation by SVA_24, SVA_33 and SVA_58 and gene downregulation by SVA_67. This is different to our previous results as the top ten reference SVAs with the greatest effect size on upregulation and downregulation were all trans-regulatory effects. Table 1 also shows that the four most significant trans-specific effects induced by SVA_15, SVA_84, SVA_87, and SVA_93 all influenced the same mitochondrial gene MTND4P24, whereby SVA_84, SVA_87, and SVA_93 increased gene expression and SVA_15 repressed gene expression. This demonstrates that multiple SVA RIPs can act on one gene influencing gene regulation, simultaneously activating and repressing gene expression.

To illustrate the impact of SVA_67 on target gene expression, we investigated the specific influence of SVA_67 allele dosage on transcriptomic profiles of the respective target gene (Figure 5). As an example, this analysis is shown for the genes MAPK8IP1P2 and LRRC37A, whereby gene expression was stratified by SVA_67 genotype (PP, PA, and AA). Analysis using the Wilcoxon pairwise comparison with FDR adjusted p-values (FDR<0.05) demonstrated that SVA_67 presence significantly downregulated gene expression of both MAPK8IP1P2 and LRRC37A (p<0.001) (Figure 5). Significant differences were observed between all genotype groups PP vs PA (MAPK8IP1P2 p=7.33E-237, LRRC37A p=1.95E-204), PP vs AA (MAPK8IP1P2 p=4.26E-55, LRRC37A p=1.06E-39), and PA vs AA (MAPK8IP1P2 p=1.33E-07, LRRC37A p=2.13E-12) for both genes (Figure 5). For individuals homozygous absent for SVA_67, a change of 262-fold and 6-fold was observed for gene expression of MAPK8IP1P2 and LRRC37A respectively, in comparison to individuals homozygous present for SVA_67 (Figure 5).

2.3. SVA RIPs demonstrate cis effects on HLA and MAPT loci

To elucidate the impact of cis-acting SVA RIPs, we explored which SVA elements had the largest significant effects on genes regulated in cis, including gene upregulation (positive beta value) and downregulation (negative beta value) using beta values from matrix eQTL analysis (FDR p<0.05) (Figure 6). Here, we highlight that SVA_67 largely influences gene regulation in cis of six different genes, including the upregulation of MAPT and LRRC37A4P and downregulation of MAPK8IP1P2, LRRC37A2, LRRC37A, and KANSL1 (Figure 6a). SVA_67 was responsible for the two greatest increases in gene activation of the MAPT and LRRC37A4P genes with beta coefficients of 1384 and 1101, respectively (Figure 6a). Interestingly, six of the ten most positive hits were effects induced on HLA genes by SVA_24, SVA_25, SVA_27 and SVA_88 (Figure 6a). Furthermore, we demonstrate that SVA_73 had the greatest effect out of all cis-regulatory effects with a beta coefficient of -2843 for the gene FCGBP, indicating a large inhibitory effect of FCGBP gene expression (Figure 6b). FCGBP gene expression is further repressed by SVA_72, the second highest hit in our analysis for the top ten SVAs with the greatest effect size on gene downregulation (Figure 6b).

To determine the influence of cis-acting SVA RIPs on gene regulation, we next analysed gene expression patterns of the top two hits for the largest effects on gene upregulation and downregulation (Figure S1). The presence of either one (PA) or two (PP) SVA_67 alleles significantly upregulated both LRRC37A4P and MAPT gene expression when tested using the Wilcoxon pairwise comparison with FDR adjusted p-values (FDR p<0.05) (Figure S1a,b). In comparison to individuals homozygous absent for SVA_67, homozygous present and heterozygous individuals displayed a 447.3-fold (p=1.04E-40) and a 214.4-fold (p=2.55E-38) increase in LRRC37A4P expression and a 1.4-fold (p=2.53E-09) and 1.2-fold (p=2.51E-04) increase in MAPT gene expression, respectively (Figure S1a,b). Analysis of the top two hits for the greatest effect on gene downregulation displayed differing gene expression patterns between the influence of SVA_73 and SVA_72 on FCGBP gene expression (Figure S1c,d). When comparing to the complete absence (AA) of SVA_73, the PA and PP genotype groups displayed a significant 3.7-fold (p=3.61E-02) and 10-fold (p=2.61E-03) increase in FCGBP gene expression (Figure S1c). Although our analysis displayed that the presence of SVA_73 upregulates the FCGBP gene in comparison to AA individuals, FCGBP gene expression is significantly repressed by 63% (p=3.29E-08) in the PP genotype group compared to the PA genotype group (Figure S1c). SVA_72 analysis highlighted that homozygous presence of the SVA_72 allele significantly downregulated FCGBP gene expression by 23% (p=5.74E-03) in comparison to homozygous absence (Figure S1d). PP genotype individuals also displayed a significant reduction of 32% (p=1.85E-05) in FCGBP gene expression in comparison to the PA genotype (Figure S1d). Similar to SVA_73 gene expression analysis, individuals heterozygous for SVA_72 indicated a 1.2-fold (p=4.23E-01) increase in FCGBP expression, however statistical significance was not achieved (Figure S1d).

As a recent genome-wide association study (GWAS) reported the association of the HLA locus with ALS risk [31,32], we investigated the influence that SVA RIPs had on cis regulated HLA genes (Figures S2 and S3). This analysis demonstrated that the top ten hits for greatest effect of SVA RIPs regarding upregulation was for HLA genes. We discovered a common gene expression pattern across all HLA genes analysed, demonstrating that individuals homozygous present for SVA RIPs (SVA_24, SVA_25, SVA_27, SVA_88) significantly upregulate HLA genes in comparison to homozygous absent individuals (p<0.01) (Figures S2 and S3).

2.4. SVA RIPs display tissue-specific modulation if gene expression

To investigate tissue-specific influences of reference SVA polymorphisms, we conducted eQTL analysis on the CNS tissue types spinal cord, cerebellum, motor cortex and frontal cortex individually (Table 2 and Table 3). The tissues temporal cortex, occipital cortex and hippocampus were excluded in the analysis due to the low n numbers. The results tables display the eQTL analysis from each individual tissue to allow determination of the top 40 hits for greatest effect size on gene upregulation (positive beta values) (Table 2) and downregulation (negative beta values) (Table 3). From this analysis, we established that a large proportion of the gene regulation effects were present within spinal cord tissue, displaying 27/40 and 26/40 tissue-specific hits for both gene upregulation and downregulation, respectively (Table 2 and Table 3). Furthermore, our analysis revealed that five SVAs (SVA_55, SVA_15, SVA_37, SVA_85, and SVA_4) were responsible for the 40 most positive beta values and eight SVAs (SVA_90, SVA_5, SVA_87, SVA_93, SVA_30, SVA_70, SVA_84, SVA_91, SVA_16) were responsible for the 40 most negative beta values, all effects of which were in the trans position (Table 2 and Table 3). This data demonstrates the capability of SVAs to modulate the expression of multiple gene targets.

Moreover, the MBP gene located on chromosome 18 and the PLP1 gene located on the X chromosome, were the most significantly upregulated and downregulated genes simultaneously (Table 2 and Table 3). Here, the gene expression of both genes was modulated in spinal cord tissue (Table 2 and Table 3). In addition, the genes MTURN and MOBP, and mitochondrial genes (MTCO1P12, MT-ND1, MT-ND2, MT-ND3, MT-ND4, MT-ND5, MT-CO2, MT-CO3, and MT-CYB) appear multiple times within both gene upregulation and downregulation hits (Table 2 and Table 3). Upon analysis of the whole tissue-specific dataset, the genes MBP, PLP1, MTURN, and MOBP were only regulated within spinal cord tissue, whilst the expression of mitochondrial genes was modulated within spinal cord, cerebellum, and motor cortex while MTCO1P12 gene expression was regulated in all four tissue types (Table 2 and Table 3). Thus, this illustrates that SVAs possess the capacity to influence tissue-specific gene expression.

3. Discussion

In this study, we evaluated the role of reference SVAs polymorphic for their presence in the human genome, to modulate gene expression within CNS tissues of ALS patients and healthy controls. Analysis of the NYGC ALS consortium dataset demonstrated that SVA RIPs significantly regulate gene expression genome-wide and in a tissue-specific manner. This study continues to expand on our previous findings, demonstrating the capability of SVAs to differentially regulate gene expression [24,33]. We have previously illustrated the capacity of SVAs to influence gene expression within Parkinson’s disease, highlighting the role of SVA_67 to modulate genes at the MAPT locus within the PPMI cohort and a CRISPR deletion model [22,26,30]. This analysis further validates our previous research, illustrating the functional capacity of SVA_67 within neurodegenerative disease and potentially expanding the importance of SVA_67 and the MAPT locus to ALS (Figure 5 and Figure 6). Therefore, this study not only emphasizes the correlation between SVA presence or absence and differential gene expression but also the involvement of SVAs within disease pathology.

Using WGS data from the ALS consortium for matrix eQTL analysis, we identified that polymorphic SVA RIPs possess the ability to modify expression of multiple target genes, but also many SVAs can regulate gene expression of a single target. Of these SVA RIPs, a greater proportion of trans regulatory effects were displayed in comparison to cis. Similarly, eQTL studies by both Wang et al. (2017) and Koks et al. (2021) analysing RNA-seq data from the 1000 genome project and the PPMI cohort, respectively, identified the impact of TEs on gene expression. In line with our analysis, these studies highlighted that several TEs simultaneously modulate gene expression of a single gene, an individual TE can regulate the expression of multiple genes, and that a greater number of TEs analysed were in the trans position [33,34]. A potential mechanism for trans-acting eQTLs could be through the binding of CCCTC-binding factor (CTCF) to SVAs. Through the cooperation with protein complex cohesion, CTCF plays a key role in three-dimensional chromatin regulation and chromatin looping, bringing promoters and regulatory elements within close proximity to activate or repress gene expression [35,36]. An additional mechanism could be through indirect transcription factor (TF) mediated associations. This suggests that SVAs influencing the expression of TF could indirectly regulate one or multiple TF gene targets, ultimately activating or repressing TF target gene expression [34].

Upon the examination of beta values from the eQTL analysis, we determined that within the combined analysis of CNS tissues SVA_55 demonstrated the greatest effect on myelin basic protein (MBP) gene upregulation (Figure 4). Following further analysis, we identified that through individual tissue analysis the top hits for gene upregulation and downregulation, were SVAs influencing MBP gene expression (Table 2 and Table 3). However, MBP expression was only significantly modulated within spinal cord tissue. MBP is a key protein involved in the myelination process, whereby myelin sheaths are formed around CNS axons by oligodendrocytes [37,38]. As oligodendrocyte loss and myelin dysfunction has recently been emphasized in neurodegenerative disease, including ALS, it is essential to investigate this relationship between SVAs and MBP differential expression [39,40]. Lorente Pons et al. demonstrated the potential significance of MBP in ALS through post-mortem analysis of both sporadic ALS and C9orf72-related ALS cases, identifying a significant reduction in MBP protein abundance when normalised to proteolipid protein (PLP) in the spinal cord corticospinal tracts in ALS cases in comparison to controls [41]. As MBP mRNA is transported to the myelin compartment by the RNA transport granule and PLP is transported as a protein, this suggests that the reduction in MBP could be due to impaired mRNA transport [41]. Our data suggests that certain SVAs (SVA_90 and SVA_5) act to downregulate MBP gene expression, therefore SVAs could play a role in this mechanism leading to the reduction of MBP protein levels observed in ALS patients compared to controls.

Our analysis also demonstrated that SVAs can activate (SVA_37) and repress (SVA_87 and SVA_93) PLP1, a form of PLP. PLP1 has previously been implicated in Pelizaeus-Merzbacher disease (PMD), an X-linked neurodegenerative disease whereby mutations within this gene inhibit CNS myelination [42,43]. In addition, we demonstrate that the myelin-associated oligodendrocyte basic protein (MOBP) gene, the locus of which has been highlighted for ALS risk, is again potentially regulated by multiple SVAs (SVA_5, 15, 37, 55, 84, 85, 87, 91 and 93) only within the spinal cord [31]. Hence, SVA regulation could be a potential mechanism involved in ALS risk at the MOBP locus. Furthermore, our tissue specific analysis displayed that multiple mitochondrial genes (MT-ND1,2,3,4,5, MT-CYB, MT-CO2,3 and pseudogene MTCO1P12) were largely activated and repressed within all four tissues analysed. SVA_30 and SVA_70 displayed the greatest regulatory effects on mitochondrial genes modulating five and four gene targets, respectively. As mitochondrial dysfunction is known to be implicated in other neurodegenerative diseases including PD and Alzheimer’s disease, as well as ALS, SVA regulation resulting in differential expression could be an underlying mechanism involved in disease pathology [44,45].

Upon analysis of CNS tissues from both ALS individuals and healthy controls, seven of the top ten cis-acting reference SVA RIPs imposed the greatest effects on the upregulation of genes at the human leukocyte antigen (HLA) locus (Figure 6A). Four SVAs are responsible for the effects on HLA gene expression, whereby SVA_24 influences the expression of one gene (HLA-A), SVA_25 two genes (HLA-B and HLA-C), SVA_27 two genes (HLA-DRB1 and HLA-DRB5) and SVA_88 one gene (HLA-DQB1). HLA also referred to as the major histocompatibility complex (MHC), acts to regulate both innate and adaptive immunity involved in the human immune response [46]. Since, the involvement of the immune response in neurological disease has been recognised the HLA locus has been highlighted as a region of importance in numerous neurodegenerative disease, including ALS [46]. Various studies have investigated the significance and mechanism of HLA in ALS, demonstrating increased frequencies of HLA-A, HLA-B and HLA-C alleles in ALS cases compared to controls [47,48,49,50]. A recent large-scale GWAS conducted by Van Rheenen et al. has identified the HLA region as a locus significantly associated with ALS, further highlighting the importance of this locus [31]. Previous studies within our group have highlighted the capability of SVAs to modulate HLA gene expression; analysis of whole genome sequencing and transcriptomic data obtained from the whole blood of individuals within the PPMI cohort discovered that SVA_24, SVA_25 and SVA_27 modulate the expression of HLA-A, HLA-B and HLA-C, and HLA-DRB1 and HLA-DRB5, respectively [33,51]. This suggests that the modulation of HLA genes by SVAs could be a common mechanism within neurodegenerative disease.

In conclusion, we show that SVAs demonstrate a significant impact on the expression of individual or numerous genes including those previously associated with neurodegenerative diseases, such as ALS. Ultimately, the ability of SVAs to act as a regulatory domain could highlight the importance of TEs in the missing heritability of neurodegenerative disease. However, due to limitations such as low n numbers for some SVA genotypes and CNS tissue types (occipital cortex, temporal cortex, and hippocampus) further research investigating the involvement of TEs in the pathogenesis of neurodegenerative disease, specifically ALS is crucial. In addition, due to low n numbers in the ALS and control groups, individual group analysis was not possible. Therefore, future experiments investigating the influence of SVAs in ALS patient and control groups individually is required. Furthermore, although this analysis continues to demonstrate the potential role of SVAs in neurodegenerative disease further experiments such as CRISPR are essential to validate SVA-specific influences. For example, we have previously shown that SVA_67 deletion in a CRISPR model resulted in a significant increase in MAPT and LRRC37A gene expression [30].

4. Materials and Methods

4.1. Genotyping reference SVAs polymorphic for their presence/absence and disease association from whole genome sequencing data from ALS Consortium dataset

Whole genome sequencing (WGS) data from the New York Genome Center (NYGC) as part of the ALS consortium dataset, were obtained in cram file format and aligned to hg38. The ALS consortium dataset contains data from individuals diagnosed with ALS spectrum MND, other neurological disorder (Alzheimer’s, Parkinson’s disease etc), other MND and ALS with other neurological disorders, along with non-neurological controls (healthy controls). The structural variant caller Delly2 (https://github.com/dellytools/delly), with default settings, was used to genotype the ALS consortium cohort (4403 individuals) for previously identified reference SVA sites, as outlined in [52]. In total 92 reference SVAs polymorphic for presence/absence [52] were used for matrix eQTL and differential gene expression analysis.

4.2. RNA-seq differential gene expression analysis

To elucidate the influence of SVA RIP genotypes on differential gene expression, RNA-seq data from CNS tissues obtained from the NYGC ALS consortium were analysed (https://www.nygenome.org/als-consortium/). For this analysis both ALS cases and healthy controls were combined. The Salmon quantification tool was used to quantify gene-specific expression levels from previously obtained FASTQ files (https://salmon.readthedocs.io/en/latest/). Salmon-generated quant files were imported using the tximport function in R and raw counts were extracted using the DESeqDataSetFromTximport. The DESeq2 R package was implemented to normalise the raw counts through the median-of-ratios method. DESeq2 was also used to detect statistically significant differential gene expression across SVA genotype groups (AA, PA, PP), whereby P represents SVA presence and A represents SVA absence. Results of the analysis were visualised using the ggplot2 package in R. The Kruskal-Wallis multiple comparison test (non-parametric) was applied to determine statistical significance between genotypes. The Wilcoxon test was used for pairwise comparison to assess significant differences between genotypes and obtained p-values were corrected for multiple comparisons using false discovery rate (FDR).

4.3. eQTL (expression quantitative trait loci) analysis

Matrix eQTL was applied to evaluate the genetic loci (SVA elements) regulating expression of nearby (cis) and distant (trans, >1 Mbp away) genes, using the matrixEQTL package in R. For this analysis, the additive linear model was applied with normalised gene expression levels, chromosomal positions of SVAs (hg38) and covariates (sex and age). Effect size estimates were reported as beta values from matrix eQTL. Obtained p-values were corrected for multiple testing and the threshold was set to 0.05 for FDR corrected p-values.