For this study, we used existing RNA-seq data for laser-captured dentate granule cells of post-mortem hippocampal tissues from the mesial temporal lobe epilepsy (MTLE) patients with and without hippocampal sclerosis (HS) generated by Griffin et al. (2016) (NCBI BioProject accession: PRJNA290212). The subject number for MTLE without hippocampal sclerosis (TLE) and with hippocampal sclerosis (TLE-HS) is 7 and 5, respectively, with the sample size being 14 and 8, respectively, due to replicated samples taken from some subjects. We also included the RNA-seq samples of bulk neocortex tissues from MTLE patients (TLE-NC) as a type of negative tissue control at a sample size of 17 from the datasets of Kjær et al. (2019) (NCBI BioProject accession: PRJNA556159). Another RNA-seq dataset for laser-captured dentate granule cells of post-mortem hippocampal tissues from healthy individuals was chosen as the normal control at a sample size of 51 (NCBI BioProject accession: PRJNA600414). Based on the original descriptions, all laser-captured hippocampal tissues were originated from dentate gyrus with the RNA-seq data generated using Illumina HiSeq 2000/2500 platforms in paired-end reads at 150 bp. The list of samples used in this study is summarized in Table 1 with the detailed NCBI sample and sequence accessions listed in Supplementary Table S1. The human reference genome (GRCh38/UCSChg38, Dec 2013) and the corresponding gene annotation (refGene.txt) files were downloaded from the UCSC website (http://genome.ucsc.edu).

The RNA-seq data for the samples listed in Table 1 and Supplementary Table S1 were downloaded from the NCBI SRA database (https://www.ncbi.nlm.nih.gov/sra) to the local high-performance computing server for de novo transcriptome assembly, identification of ME sequences, and position correlation with annotated genes. De novo assembly was performed using SOAPdenovo-Trans (Xie et al., 2014) with the default parameter settings keeping only transcripts that are 300 bp or longer.

To identify transcripts containing ME sequences, all assembled transcripts from the above step were subjected to repeat annotation using RepeatMasker (version 4.1.0, http://www.repeatmasker.org). Any transcripts for known genes containing one or more MEs are referred to as ME-transcripts. The individual sequences of all ME-transcripts were extracted from the transcript assemblies using an in-house tool, fatools (available at https://github.com/pliang64/fatools). Extracted ME-transcript sequences for each sample were then aligned to the human reference genome (GRCh38/UCSChg38, Dec 2013) using pblat (Kent, 2002; Wang and Kong, 2019) to obtain their aligned positions in the reference genomes. The blat output was processed using an in-house Perl script to generate a delimited text file, which provides the chromosomal position, classification of the MEs in the ME-transcripts, the symbols and positions of the genes in which the ME(s) reside, and the gene regions to which MEs correspond. This was performed based on the position of MEs in the transcripts and the position of the transcripts in the reference genome and the location of gene exons in the reference genome. The human reference genome sequence and the corresponding gene annotation (refGene.txt) files were downloaded from the UCSC website (http://genome.ucsc.edu).

For protein coding genes, the distribution of MEs in transcripts is also further broken down into 5′ UTR, CDS, and 3′ UTR regions. The frequency of gene transcripts containing ME sequences (ME-transcripts) in each sample was normalized as the number of ME-transcripts per million transcripts (TPM). TPM values were also calculated for each major ME class (SINE, LINE, LTR, SVA, and DNA transposons) and for the subfamilies. More related details are provided in Supplemental Methods, while the statistical analysis was performed as described in section 2.5.

To find out whether the presence of the ME-transcripts shows any bias toward certain types of genes, we performed enrichment analysis for genes associated with ME-transcripts for each sample group individually using the DAVID tool (version 6.8) (https://david.ncifcrf.gov) (Huang et al., 2007). Specifically, a nonredundant list of genes associated with ME-transcripts in all individual samples from a group was generated and used for enrichment analysis. For enrichment analysis, gene ontology (GO terms), KEGG pathways, and additional functional categories including UP_KEYWORDS, UP_SEQ_FEATURE, and UP_TISSUES provided by the DAVID tool were included. We collected all categories showing statistically significant enrichment with a Benjamini adjusted p-value below 0.05.

We collected a list of 977 genes which are known to be associated with epilepsy compiled by Wang et al. (Wang et al., 2017) and a list of 2,449 genes associated with neurogenesis and neural system development from the MGI database (http://www.informatics.jax.org). The numbers of epilepsy and neural genes associated with the ME-transcripts were collected for each sample by cross-mapping between the list of genes associated with the ME-transcripts and the epilepsy associated and neural gene lists for comparison among groups. For a selected number of ME-transcripts associated with epilepsy genes, detailed sequence analysis was performed to predict the functional impact of the ME sequence on the host genes. For the SCN1A gene, the expression level of the ME-form in comparison with the normal splicing forms was analyzed and compared across sample groups. A more detailed method description for this part is provided in the Supplementary Methods.

For a selected list of epilepsy-associated genes involving MEs in CDS, we performed detailed sequence analysis of the ME-transcripts manually using the resources and utilities on the UCSC genome browser (https://genome.ucsc.edu). Specifically, the entire ME-transcript sequences in the fasta format were retrieved from the transcript assemblies of the corresponding samples and used to search the reference genome and locate the MEs in CDS exons as a way of validation of the predicted CDS location of the MEs and predicting the impact on protein coding. Furthermore, for each of such ME-transcripts, we extracted the 20 bp sequences flanking the junction point (10 bp on each side), at which the normal transcript(s) and the ME-transcript differentiate due to the presence of an ME for both the normal and ME-transcripts. This pair of 20 bp sequences was then used to search the entire set of RNA-seq raw reads to count the number of reads supporting each of the two forms of transcript sequences at that specific location. These counts were collected for each sample to calculate the expression level of the gene in the ME-transcript and normal transcript as reads per million reads (RPM) and the rate of ME-transcript as its percentage in the total reads for the gene at the ME junction position. The RMP values and ME-transcript rate for a gene were compared across sample groups with statistical testing performed using an unpaired 2-tailed t-test.

Most of the analysis for RNA-seq data was performed on Compute Canada high-performance computing facilities (https://computecanada.ca). Downstream data analyses and figure plotting were performed by using a combination of R package, GraphPad Prism, and Microsoft Excel. Statistics was performed using SAS 9.4 (SAS Inc, Cary, NY, United States). Data were reported using descriptive statistics. Variables were presented as the arithmetic mean plus standard deviation (mean ± STD). To assess the statistical significance of the variables across different sample groups, one-way analysis of variance (ANOVA), followed by Student-Newman-Keuls post hoc tests (SNK-q test), was used for those showing normal distribution, while the Kruskal–Wallis rank sum test, followed by the Nemenyi test or Wilcoxon rank sum test, was used for those showing non-normal distribution, and the Chi-square test was used for categorical variables. All p-values were from two-sided tests, and the results were deemed statistically significant at different levels of p < 0.05 ( _* ), p < 0.01 ( _** ), or p < 0.001 ( _*** ).

As shown in Table 2 and Figure 1, four classes of MEs, including SINE, LINE, LTR, and DNA transposon, were observed in the gene transcripts from all four sample groups (control, TLE-NC, TLE, and TLE-HS). It is interesting to notice that SVA, as the youngest ME class in the human genome, is completely missing from the ME-transcripts. As shown in Figures 1A–D, except for SINEs, all other three classes of MEs showed a gradient of the ME frequency from low to high in the order of the normal control, TLE-NC, TLE, and TLE-HS. For SINEs, the only difference from the other ME classes is that the frequency in TLE-NC is higher than that of the TLE. However, for all ME classes, it is always that the TLE-HS group has the highest frequency of ME-transcripts, while the normal control group has the lowest, when all transcripts (including noncoding transcripts) were considered. The differences were statistically significant (p < 0.05) for the TLE and TLE-HS groups in comparison with the normal control group for all four ME classes (Figures 1A–D). It is worth noting that except for SINEs, the TLE-NC group, which was included as a control brain tissue from the epilepsy patients, showed a ME-transcript frequency that is lower than that of TLE and TLE-HS and higher than (except for DNA transposons) that of the normal control group, but without significant difference (Figures 1B–D).

We also examined the distribution of MEs across the three sub-exon regions of protein coding genes. As shown in Table 2 and Figures 1E–H, the 3′ UTR region showed the highest frequency of ME-transcripts, being more than 10 times higher than that in the CDS and 5′ UTR regions. The differential pattern among the sample groups in the 3′ UTR region (Figure 1F) is similar to that of all exon regions (for all genes) (Figure 1E). For the CDS regions (Figure 1G), the TLE-HS group has the highest ME-transcript number when all ME classes are considered, being significantly higher than that of TLE-NC, with the latter being the lowest among the four sample groups. The 5′ UTR regions had the lowest number of MEs among the three exon regions, likely due to their much shorter sequences. In contrast to the other two regions and all regions, the ME number in 5′ UTR regions was shown to be the lowest in the TLE-HS group, but with no statistical significance between any sample groups (Figure 1H).

More detailed comparisons were made among the three exon regions with MEs broken down to ME classes (Table 2). Overall, the frequency of SINEs was always shown to be the highest in the entire transcripts and for the CDS and 3′ UTR regions. In comparison, the frequency of LINEs was shown to be the highest in TLE-HS only for 3′ UTR. For LTRs, the highest frequency was seen in 3’ UTR regions of TLE-HS, but for CDS, the frequency in TLE-HS is the lowest among all sample groups. For DNA transposons, the highest frequency was seen in TLE-HS for all three exon regions (Table 2). These results indicate that the distribution patterns in the transcripts by gene context shared some similarities but alsoshowed differences among different ME classes.

Among the ME classes, SINE is the most frequent type in all sample groups, ranging from 3,734 TPM in the control group as the lowest to 9,983 TPM in the TLE-HS group as the highest, when all exon regions were considered (Table 2 and Supplementary Figure S1A). Following SINE are LINE, DNA, and LTR in the order from high to low with the rates of LTRs being more than 4 times lower than that of SINE. This order of abundance among ME classes is the same for all sample groups and in all three exon regions of protein coding genes, except for the TLE-NC group, which showed a slightly higher rate of LTRs than DNAs in the 5′ UTR and CDS regions (Table 2). The higher rates for SINEs and LINEs can be explained by their higher copy numbers in the human genome [1,779,233 and 1,516,226, respectively (Tang et al., 2018b)]. However, DNA transposons showed a higher rate in ME-transcripts than LTRs despite having a lower copy number in the genome than LTRs [483,994 vs 720,177 (Tang et al., 2018b)] (Table 2 and Supplementary Figure S1).

To assess and compare the degree of involvement of MEs in transcripts across ME classes proportionally by their abundance in the genome, we normalized the TPM values in Table 2 as TPM per million MEs (TPMPM) from the ME class in the entire human genome. As shown in Supplementary Table S2 and Supplementary Figure S1A, after this normalization for all exon regions, while the relative portion dropped, SINE was still shown to be the most frequent ME class, ranging from 3,162.2 TPMPM in the normal group as the lowest to 8,452.8 TPMPM in the TLE-HS group as the highest. The relative portion of LINEs and LTRs did not seem to change much by the normalization (Supplementary Figure S1A). However, interestingly, the frequency of DNA transposons in TPMPM increased in all sample groups from their TPM values and became the second highest among the ME classes, being higher than that of LINEs and LTRs in all sample groups and being very close to that of SINEs in the TLE group (Supplementary Table S2 and Supplementary Figure S1A). This indicates that proportionally, DNA transposons were shown to be more active contributors to ME-transcripts than LINEs and LTRs, with the largest increase by normalization seen in the TLE-HS group (Supplementary Table S2 and Supplementary Figure S1A). We also examined the situation in the CDS regions and observed a similar increase for DNA transposons when normalized as TPMPM (Supplementary Figure S1B). We would expect this also to be the case for the 3′ UTR regions (data not shown).

We further analyzed the frequency of ME subfamilies within each class and identified subfamilies contributing to differential ME subfamilies. Among the 36 ME subfamilies detected in epileptic transcriptomes, 10 subfamilies showed significant difference in their frequencies in ME-transcripts at least for one pairwise group comparison. As shown in Supplementary Table S3, the 10 subfamilies include Alu and MIR from SINEs; L1, L2, and CR1 from LINEs; ERVL-MaLR, ERV1, and ERVK from LTRs; and hAT-Charlie, TcMar-Tigger, and hAT-Tip100 from DNA transposons in order from high to low within their respective ME classes.

To better understand the contributing factor(s) related to the degree of MEs’ involvement in gene transcripts, we performed age profiling of these MEs based on sequence divergences from their respective consensus sequences. In this case, the higher the divergence, the older the age of the ME in the genome.

Overall, as expected, the age profile is different among ME classes with SINEs showing the youngest peak among all ME classes, followed by DNA transposons and LTRs and then by LINEs as the oldest (Supplementary Figure S2). For SINEs, the age profile is pretty much identical among the three datasets (the genome, the normal control group, and the TLE-HS group), indicating that their degree of involvement in ME-transcripts by the age profile is simply a reflection of what is in the genome and is not impacted by the disease status. For DNA transposons, the two sample groups showed a younger age profile than the genome-wide profile, indicating a bias for younger elements in ME-transcripts. LTRs showed a very wide and flat peak in the genome, likely reflecting a continuous LTR domestication process and/or steady propagation during the human genome evolution. However, LTRs showed a single main peak in the control sample group somewhere in the middle for the flat peak for the genome but two peaks in the LTE-HS group with one being younger and the other being older than the main peak in the control group (Supplementary Figure S2), indicating a unique age profile for those involved in the LTE-HS transcriptomes. For LINEs, the genome profile showed a clear peak toward the old age, while the control group showed two peaks, both being younger than that of the genome profile. However, the TLE-HS group showed a main peak aligned to the older peak in the normal group and a secondary peak at a much younger age, which happens to be of the same age as the main SINE peak. This indicates that LINEs in ME-transcripts in the TLE-HS group have a higher proportion contributed from young members than the control (Supplementary Figure S2). Overall, as the general pattern, while TLE-HS seemed to involve more younger MEs, older MEs are the main players in ME-transcripts, agreeing also with the lack of SVAs in these ME-transcripts.

To find out whether the impact of MEs on RNA splicing is completely random or selective for a certain type of genes, we performed enrichment analysis of genes associated with ME-transcripts for gene function categories including GO terms, KEGG pathways, and other category terms provided by the DAVID tool (Huang et al., 2007). The list of categories showing statistically significant enrichment in at least one of the sample groups was collected and is shown in Figure 2A. It appears that a common theme among the enriched function categories is the mRNA splicing process, involving a total of 206 genes from all sample groups. In addition to the first four categories with descriptions directly associated with RNA splicing, all remaining categories share most of their genes with these categories, thus making the RNA-splicing process as the main theme of the enriched function category. The other three enriched keywords (phosphoprotein, zinc finger region: C2H2-type 2, and acetylation) may reflect other specific characteristics of these genes. No extra categories of enrichment were observed in LTE or LTE-HS over the normal control, despite their specific gene lists being similar but not identical. However, interestingly, the TLE-NC (neocortex) group showed the least enrichment among the sample groups for these function categories with none being statistically significant, indicating that neural cells in the neocortex are quite different from those in the hippocampal tissue at least for mRNA splicing. The enrichment for brain tissues in all sample groups can be considered expected (Figure 2A) since these genes were detected to be expressed in brain tissues as a starting point of the analysis.

A list of epilepsy-associated and neural genes and those impacted by MEs in the CDS regions in the epileptic sample groups was generated. In total, 15 genes were detected to harbor MEs among the epilepsy-associated genes, and genes which are more frequent to harbor MEs in TLE include BRD2, DOCK7, HERC2, IGSF8, NF1, NIN, PCLO, RFX3, SCN1A, SCNM1, and TBCK (Supplementary Table S4).

We further compared genes associated with ME-transcripts against the lists of epilepsy-associated genes and neural genes to see whether there are any differences between the control and epileptic groups by frequency. Specifically, we calculated the number of genes from each of the two gene lists (2,449 and 977 for neural and epilepsy-associated genes, respectively) in each sample in a group and compared among the groups. In this case, we only included genes with MEs in the 3′ UTR and CDS regions of the transcripts. A nonredundant list of genes associated with ME-transcripts in the CDS regions for the epileptic groups is provided in Supplementary Table S4.

As shown in Supplementary Table S5 and Figures 2C–E, for both neural and epilepsy-associated genes, the TLE-HS group always showed a higher number of genes impacted by MEs than all other groups in both the 3′ UTR and CDS regions. Statistically significant differences were observed for neural genes in CDS regions and for both neural and epilepsy-associated genes in the 3′ UTR regions in comparison with the normal control and/or TLE-NC group (Figures 2B,D,E). It is interesting to note that the TLE-NC group showed the lowest number of genes impacted by MEs in CDS regions from both gene lists with statistical significance for its differences from the two TLE groups (Figures 2B,C).

To understand the functional impact of the MEs in gene transcripts, we took the ME-transcripts for genes known to be associated with epilepsy as examples and examined in detail with the transcript sequences. In almost all cases we examined, the MEs in the transcripts represent those normally located in the intron regions but were retained to generate extended exons (exon extension) as shown in Figure 3A–C for the SCN1A gene as an example. The outcome of such abnormal splicing involving MEs in the CDS region is the disruption of the normal open reading frame via generation of a premature stop codon (Figure 3D). The predicted functional impact for SCN1A is a truncated protein product, which misses more than 1,000 amino acids at the C-terminus (Figure 3E), and/or mRNA degradation by the (NMD) machinery (Chang et al., 2007), both expected to result in loss of functions for the gene.

We further examined to see whether, for genes impacted by abnormal ME-transcripts, the rate of the abnormal splicing form is differentially expressed among the sample groups. As shown in Figure 3F, while the ME-transcripts of the SCN1A gene were observed in all sample groups, the rate measured as the percentage of ME-transcripts in all transcripts for the gene was shown to be more than 2 times higher in the TLE-HS and TLE groups than in the control group (TLE at 27.0%, TLE-HS at 26.3% vs control at 10.9%) despite not having statistical significance, possibly due to the small sample size for the TLE and TLE-HS groups and/or the large relatively variations within the groups. Notably, the rate of the ME-form in TLE-NC (3.0%) was shown to be statistically lower than that of all other three groups (p < 0.05) (Figure 3F). What is also interesting is that the TLE and TLE-HS groups showed a significantly lower expression level of the normal SCN1A transcript form than the control and TLE-NC groups (Figure 3G). This is also true for the total level of expression (Figure 3H), likely due to the high ratio of these abnormal splice variants, which are likely targeted by NMD. For this reason, we suspect that the observed rate of the ME-form in this example might be an under-representation of what actually was happening, and this could explain the overall lower expression of the genes in the epileptic groups.