Exploratory Analysis of L1 retrotransposons expression in Autism

doi:10.21203/rs.3.rs-2233541/v1

Download PDF

Research Article

Exploratory Analysis of L1 retrotransposons expression in Autism

https://doi.org/10.21203/rs.3.rs-2233541/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 28 Jun, 2023

Read the published version in Molecular Autism →

You are reading this latest preprint version

Background

Autism Spectrum Disorder (ASD) is a set of highly heterogeneous neurodevelopmental diseases whose genetic etiology is not completely understood. Several investigations have relied on transcriptome analysis from peripheral tissues to dissect ASD into homogenous molecular phenotypes. Recently, alterations of gene expression from post-mortem brain tissues have identified sets of genes with altered expression, mapping on pathways known to be involved in ASD etiology. In addition to protein-coding transcripts, the human transcriptome is composed by a large set of non-coding RNAs and transposable elements (TEs). Advancements in sequencing technologies have proved that TEs can be transcribed in a regulated fashion, and their dysregulation have a role in neural diseases.

Methods

We exploited published datasets comprising RNA-seq data from 1) post-mortem brain of ASD subjects, 2) in vitro cell cultures where ten different ASD-relevant genes were knocked out and 3) blood of discordant siblings. We measured the expression levels of evolutionarily young full-length transposable L1 elements and characterized the genomic location of deregulated L1s assessing their potential impact on the transcription of ASD-relevant genes.

Results

We detected a strong upregulation of intronic full length L1s exclusively in a subset of samples where the number of deregulated genes correlates with L1 upregulation. In brain tissue, specific upregulated L1s overlap ASD-relevant downregulated genes, suggesting the existence of a negative effect of L1 transcription on host transcripts for a subset of ASD genes.

Limitations

Our analyses are limited by the low number of samples and by the lack of replicates among post-mortem brain samples. Measuring the transcription of locus-specific TEs is complicated by the repetitive nature of their sequence, which reduces the accuracy in mapping sequencing reads to the correct genomic locus.

Conclusions

We observed L1s upregulation in a subset of ASD samples. In these cases, L1 upregulation might directly impair the expression of ASD-relevant genes by a mechanism not yet understood. L1s upregulation in autism appears to be limited to a subset of subjects probably harbouring mutations in genes related to DNA/chromatin modification.

Autism

Transcriptomics

L1 retrotransposons expression.

Autism Spectrum Disorders (ASD) are a set of heterogeneous neurodevelopmental conditions mainly involving impaired communication and repetitive behaviors(1). ASD is the most common neurodevelopmental condition in human: the median worldwide prevalence of autism is around 1% according to the latest world-wide surveys(1). The etiology of autism presents a strong genetic component: twin and sibling studies have consistently shown that ASD is one of the most highly heritable complex disorders in humans(2). In the past decade several studies exploiting whole-exome sequencing (WES) data highlighted both de novo and inherited deleterious mutations which are either causative or contributing to the autistic phenotype(3–9). The increasing scope of WES studies and the ever-growing cohort sizes expanded the discovery rate of ASD genes, to the point where up to 1,000 genes are estimated to contribute to a different degree to ASD etiology(10,11). Yet, none of these genes accounts for more than 1% of idiopathic ASD cases(12). This extreme variability may be at the basis of the considerable phenotypic heterogeneity characteristic of ASD. Hence, the identification of subgroups with a more homogeneous molecular asset is essential to comprehend ASD etiology and program personalized treatments. Understanding how the large number of genes implicated in ASD susceptibility may converge to affect human brain development is critical(12). Similar to other neuropsychiatric disorders, most of the genes involved in ASD encode for neuronal components crucial for brain function(10,11). In addition, a relevant portion of genes are involved in general transcriptional regulation and/or chromatin remodeling(10,13). The link between deleterious mutations affecting a chromatin-regulatory gene and those involved in synaptic transmission and brain activity remains unclear, although a reasonable hypothesis suggests chromatin-regulatory gene mutations may affect transcriptional programs impacting genes involved in synaptic transmission and brain activity. Several recent works on post-mortem brain tissues were able to reveal shared abnormalities in gene expression in a large subset of autism cases(14–16). Namely, two main types of gene co-expression modules were consistently deregulated: 1) down-regulated genes involved in synaptic transmission, encoding for neuronal markers and enriched for ASD-associated genes; and 2) modules of up-regulated genes involved in immune and inflammatory responses, enriched for markers of microglia and astrocytes, but generally not for genes associated with ASD(14–18). Multiple independent studies observed perturbations of epigenetic marks distribution in post-mortem brain tissue of individuals with ASD(19–21). Both Wong et al(16)^⁠and Nardone et al(19)^⁠ reported an altered DNA methylation landscape among multiple brain regions of ASD individuals. Corley et al(21) showed that epigenetic alterations detected in ASD were preferentially directed at intragenic and bivalently modified chromatin domains of genes predominately involved in neurodevelopment. Interestingly, the methylation landscape in adult neurons affected by ASD closely resembled the pattern of earlier time points in fetal brain development(21). These findings suggest that a delay in the epigenetic program can contribute to deleterious transcriptional programs and to the establishment of ASD phenotypes. This model is supported by the identification of mutations in genes involved in the regulation of both DNA and histone methylation during brain development(13). A recent work by Wong et al(16)^⁠ integrated data from gene expression, DNA methylation and histone acetylation from ASD and healthy individuals, proposing the existence of two major subgroups of ASDs. The first subgroup recapitulated all known molecular changes typical of ASD(16). The second one was indistinguishable from control samples in terms of transcriptional and epigenetic alterations. It is therefore tempting to stratify deleterious genomic variants by classifying their molecular effects into subgroups at higher resolution. In this scenario, a recent publication identified a set of human regulatory regions evolved after the separation from old world monkeys(22) and highlighted how these were enriched within genomic regulatory regions altered in ASD(22). Indeed, genomic loci relative to hominoid-specific regulatory gains showed a significant overlap with regions that lose the H3K27ac histone mark, typical of active enhancers in ASD subjects(22). This evidence suggests that ASD-related epigenetic defects may be caused by altered activity of evolutionarily young regulatory regions.

Transposable elements (TEs) are genomic sequences able to move within the genome(23) and to shape the regulatory landscape of the host genome(24–27). Long Interspersed Nuclear Elements 1 (L1s) are the most abundant TEs (~17% of the human genome) and the only transposon class known to retain the ability to mobilize autonomously in human(28). Full-length (FL) L1 elements are about 6 Kbp long and take advantage of a “copy and paste” mechanism where a full-length sequence gives rise to a L1 RNA intermediate that is reverse-transcribed into a new genomic locus. FL L1s contain two ORFs (ORF1 and ORF2), encoding for respectively a nucleic acid chaperone and a protein with endonuclease and reverse transcriptase activity, which together mediate retrotransposition(28). Retrotranspositional events mostly produce 5′-truncated L1s that are unable to re-mobilize. Only around 100 of the more than 10,000 full-length L1s found in the human genome(29) are potentially active(30). To prevent potential deleterious effects of L1 abnormal activity, cells have developed several mechanisms to safeguard and fine-tune L1 retrotransposition. These include DNA methylation, transcriptional repression and L1s RNA degradation though the activity of the PIWI/piRNA pathway(31). Nevertheless, TEs are known to escape silencing at embryonic stages(32), affecting early human development by regulating nearby protein-coding genes. Waves of hypomethylation during embryogenesis are linked to higher rates of transcription and retrotransposition of L1 RNAs(33). Somatic L1 retrotransposition has been observed in the neuronal lineage leading to brain mosaicism and it has been shown deregulated in a plethora of neurodegenerative and neurodevelopmental diseases(33–35) although its extent and functional significance remain unclear(33,36,37). Recently, attention has been focused on the functional role of L1s independently from retrotransposition(38). The majority of L1 RNAs are retained in the nucleus, and they can function as regulatory long non-coding RNAs (lncRNAs), controlling transcriptional and chromatin landscapes. For example L1 RNAs are required for mouse embryonic stem cell (mESC) self-renewal and pre-implantation during development(39). L1 RNA expression and mobilization must therefore be considered two independent events under distinct regulatory pathways and with different functional outcomes(40,41).

Altered DNA methylation levels within L1 sequences have been shown in multiple neurodevelopmental diseases, including ASD(42). The reduction of methylation and an increase in L1 expression was reported in ASD post-mortem brains(42,43). In this tissue, trimethylation of histone H3K9 (H3K9me3), which is responsible for the formation of condensed heterochromatin and prevents L1 activation, was significantly reduced at L1 ORF1 and ORF2 sequences but not at the 5′-UTRs in ASD samples(42). Furthermore, Tangsuwansri et al(43) demonstrated that, in lymphoblastoid cell lines derived from a subset of ASD subjects with severe language impairment, the overall methylation level of L1 elements was decreased compared to controls, and this was inversely correlated with the level of expression of L1-containing genes(43).

In this study, we aimed at quantifying the expression of young and transcriptionally active L1 elements in ASD. Specifically, we wanted to assess whether specific individuals affected by ASD or in vitro models of ASD show an altered L1s expression. Furthermore, we wanted to characterize and explore the transcriptional dynamics of L1 elements and the possible impact of their expression on ASD-relevant genes.

Evolutionarily young L1 subfamilies are expressed in post-mortem brains, in in vitro models and in blood

L1 transcriptional deregulation is observed in multiple human disorders including ASD(42,43), through mechanisms that are not well understood. To explore the potential impact of L1 elements in ASD, we assessed whether and which L1 subfamilies are expressed in samples and experimental models of ASD. We retrieved public bulk RNA-seq data from post-mortem ASD brains(44). The dataset consists of a total of 31 individuals (from 15 ASD donors and 16 controls) with expression profiling of anterior cingulate cortex (ACC) and pre-frontal cortex (PFC). Data from both tissues was available for 10 individuals (7 ASD and 3 controls(45)), while for the remaining 21 individuals (8 ASD and 13 controls) expression data were produced from either ACC or PFC. Therefore, analyses were performed on a total of 41 RNA-seq derived samples, 18 for ACC (9 ASD and 9 controls) and 23 for PFC (13 ASD and 10 controls). We also retrieved RNA-seq data from induced pluripotent stem cell (iPSC) and neuron cells knockout (KO) for ten ASD-relevant genes (AFF2, ANOS1, ASTN2, ATRX, CACNA1C, CHD8, DLGAP2, KCNQ2, SCN2A, TENM1)(46). The dataset comprises a total of 87 RNA-seq samples, 42 from iPSCs and 45 from neurons. KOs were performed for each ASD-relevant gene through a CRISPR gene editing method to insert premature termination sites(46). Furthermore, we obtained RNA-seq data from whole blood of 73 pairs of discordant ASD siblings for a total of 146 samples.

To quantify L1 expression two software were used: TEspeX(45)andSQuIRE(47). TEspeX measures TE expression not counting reads possibly generated from TE fragments embedded in the exons of annotated transcripts, and cumulatively quantifies the expression of TE subfamilies annotated in RepeatMasker(29) by counting the sequencing reads mapped to the consensus sequence of each TE subfamily. SQuIRE(47) allows to count reads mapped to the genomic loci where annotated TEs are found. We quantified the expression of a total of 117 human L1 subfamilies. Firstly, we used TEspeX to determine which L1 subfamilies are more highly expressed independently from the host transcripts. We then exploited SQuIRE to pinpoint which specific L1 fragments belonging to the most expressed subfamilies are transcribed at the locus-specific level. For each dataset, we defined as expressed all L1s fragments with an average of at least 200 reads mapped among all samples. According to our analysis, L1HS is the L1 subfamily with the highest expression level in all datasets, followed by several L1PA subfamilies (Additional file 1). In ASD brain samples, KO cell lines and all their controls, L1HS and L1PA subfamilies make up ~80-90% of the total normalized L1 expression, while all other L1 subfamilies cumulatively show a modest overall expression (Figure 1A). On the other hand, L1HS and L1PA subfamilies constitute only about 35-40% of the total normalized L1 expression in the blood dataset, for both ASD and controls (Figure 1A).

Evolutionarily young TEs would be expected to show a lower number of mutations which might increase their probability to be transcribed compared to older TEs. A measure of the evolutionary age of TEs is measured as base mismatches in parts per thousand (MilliDiv). It is proportional to the sequence divergence of a given element compared to its consensus sequence. We retrieved MilliDiv values for each L1 fragment annotated in the human genome (version hg19) from RepeatMasker(29) and calculated the average MilliDiv for each subfamily. As expected, the average MilliDiv for L1HS and L1PA fragments is significantly lower compared to that of the fragments of other L1 subfamilies and, interestingly, the average MilliDiv for the group of expressed L1 elements is lower than the corresponding group containing also the non-expressed ones (Figure 1B). We therefore decided to focus our analyses on L1HS and 20 L1PA subfamilies.

L1 elements result upregulated in post-mortem brains of specific ASD individuals and ASD model samples

We wanted to assess whether any of our ASD samples would show altered expression of young L1 elements compared to controls. We therefore measured the expression of L1s annotated in the human genome (version hg19) with SQuIRE(47) in all samples, and performed L1-related differential expression analysis independently comparing ASD samples versus controls. To reduce the risk of quantifying the expression of small L1 fragments whose expression may be strictly linked to their host genes, we limited our analysis on likely full-length (>5,000 nts) fragments showing an average of at least 200 reads mapped among samples of each dataset and belonging to L1HS and L1PA subfamilies. A total of 128,506 L1HS and L1PA fragments are annotated in RepeatMasker(29) (hg19), of these, 9,077 are more than 5,000 bp long. From SQuIRE analysis we classified 100 L1s for Velmeshev et al and 175 L1s for Deneault et al as expressed (>200 average reads mapped among all samples) (Additional File 2). Of them 36 were upregulated in ASD brain samples or KO cell lines compared to controls in both datasets. We then performed L1 expression analyses on each single individual (see methods). Briefly, the number of reads mapped to each L1 fragment was first normalized with DEseq2(48). To infer for differential expression in each single individual in the absence of replicates we calculated a Z-score for each L1 fragment comparing its expression level in each sample against the distribution of its expression levels in the group of controls. For each sample, we classified as potentially upregulated all L1 fragments characterized by a Z-score > 3, and as potentially downregulated all L1 fragments characterized by a Z-score < -3. We then calculated the net number of expressed fragments by subtracting the number of downregulated L1s to the number of upregulated L1s for each sample. Overall, ASD brain and KO samples always show a positive net number of expressed L1s compared to controls. In most control samples the net number of expressed L1s is either null or very low but positive, except for two control samples of the Velmeshev et al ACC dataset, which show a moderately negative value (Figure 1c). In the case of the Velmeshev et al dataset, only three ASD samples in the ACC showed a considerably higher amount of L1 net expression compared to all other samples, showing a net number of about 60 upregulated L1s, while all other samples do not show more than 25 potentially upregulated L1s (Figure 1C). On the other hand, no PFC sample shows a substantially increased number of upregulated L1s compared to all other samples. Furthermore, the three samples characterized by strong L1 upregulation compared to all other samples in ACC (SRR9292614, SRR9292620, SRR9292621) do not show the same trend in their PFC counterpart (SRR9292613, SRR9292619, SRR9292623). Therefore, the observed L1 upregulation probably occurs only in specific areas of the brain. In the Deneault et al dataset, both iPSC and neurons KO for ATRX show a substantially higher mean net number of expressed L1s (about 60) compared to the other KOs, which do not exhibit an average net number of expressed L1s higher than 30 (Figure 1C). It is interesting to note how in all KO samples the average net number of expressed L1 compared to controls is always positive. We consider this result as an indication of L1 upregulation, suggesting that the loss-of-function of ASD-related genes might lead to an increase in L1 transcriptional activity.

We performed the same L1 expression analyses on a total of 146 RNA-seq whole blood samples derived from individuals affected by ASD and their unaffected siblings. However, no L1 fragment showed a Z-score above 3 in any ASD or control sample. We then plotted the cumulative expression of L1 subfamilies quantified with TEspeX. Overall, all blood samples (ASD and controls) show similar total L1HS/L1PA expression levels (Supplementary figure 1A). Furthermore, the distribution of the ratio of the total L1HS/L1PA normalized expression between each pair of discordant siblings (ASD/Control) show the absence of alterations in L1 levels in the blood of discordant siblings (Supplementary figure 1B).

Taken together, these results suggest that only a fraction of ASD samples display a strong increase in the transcriptional activity of young full-length L1 elements probably specific to specific regions of the brain.

Evaluation of L1 expression in the context of mutated genes

iPSCs and neurons KO for ATRX display a pervasive upregulation of young FL L1 elements, while the KO of nine other ASD-relevant genes leads to a more limited increase of L1 expression. Hence, we speculate that the loss-of-function of specific genes may strongly affect L1 transcriptional regulation. ATRX encodes for a protein which contains an ATPase/helicase domain belonging to the SWI/SNF family of chromatin remodeling proteins(49,50). It has been demonstrated that ATRX knockout mouse ES cells exhibit increased chromatin accessibility at genomic loci occupied by retrotransposons(49)⁠ suggesting that this gene has a direct role in the establishment and maintenance of heterochromatin also at the level of L1s retrotransposons in human. To validate this hypothesis, we retrieved ChIP-seq data from a recent publication(50), where the authors performed a comprehensive ChIP-sequencing on ATRX in human cell lines. To assess whether ATRX may specifically act on genomic regions occupied by upregulated FL L1s, we computed a profile of ATRX ChIP-seq enrichment separately on the genomic coordinates of total and upregulated loci compared to controls L1HS/L1PA, as well as of 1,000 of 6,000 bp long genomic coordinates randomly selected within introns as controls. The profile was computed with DeepTools(51)⁠ computeMatrix utility and plotted with plotProfiles. In both iPSC and neurons, total and upregulated L1HS/L1PA show a strong enrichment compared to both their up/down-stream boundaries and to random intronic segments (Figure 1D). Together with the strong L1 upregulation observed in ATRX KO cells, these results confirm that ATRX loss-of-function may directly lead to an increased transcription of young FL L1 elements.

We then asked whether ASD - ACC post-mortem samples characterized by strong L1 upregulation might be affected by disruptive mutations within genes having a role in repressing TEs transcription. We took advantage of the public blood-derived whole exome data samples of ASD individuals from Velmeshev et al(44)⁠ to search for putative deleterious mutations within ASD-relevant genes. Upon aligning raw fastq RNASseq reads to the reference genome (version hg19), the GATK4 variant discovery pipeline was used to detect all non-reference short variants, following GATK4 best practices(52) as described in methods. The variant discovery phase resulted in the detection of ~35,000 total non-reference variants for each sample. A filtering strategy was then employed to select a total of about 20 candidate strong deleterious variants per sample. We then selected only the variants within SFARI genes (SFARI score 1, 2, 3 or S)(11)⁠finding 16 potentially deleterious variants within SFARI genes in a total of 10 samples (Additional file 3). While we were not able to find a suggestive candidate variant for all ASD individuals, the three samples from the Velmeshev et al dataset characterized by a strong L1 upregulation bear candidate variants within genes DNAH3, PRUNE2, MBD1 and LAMA1. Although none of these genes has yet been specifically linked to the transcriptional control of TEs, further studies are needed to fully assess a possible link.

Transcriptome-wide level of deregulation in ASD brain and in vitro neurons correlates with L1 upregulation

Altered gene expression is currently considered as one of the main molecular manifestations of ASD(14–16,18). Since specific ASD samples show strong L1 upregulation compared to controls in independent datasets, we wanted to investigate whether altered L1 expression is associated to gene expression alteration. To define differentially expressed (DE) genes, we used the same approach applied to FL L1HS/L1PA fragments, which allowed us to characterize each sample individually in comparison with controls in absence of replicates. We applied this method only to genes associated to an average number of mapped reads above 200 across samples of each dataset. We refer to this subset of genes as expressed genes. Genes were considered DE if associated to a |Z-score| > 3. Overall, the number of DE genes is not homogeneous among samples, with some ASD samples showing a number of DE genes very close to the controls (Figure 2A). It is interesting to note that samples which showed stronger L1 upregulation are among the samples with the highest number of DE genes, especially in post-mortem brain and in vitro neurons (Figure 2A). For example, sample SRR9292620 is associated with both the highest number of net expressed L1s and with the highest number of DE genes. Furthermore, the top six samples of the Velmeshev et al - ACC dataset with the highest L1 net expression are at the same time the six samples with the highest number of DE genes. Similarly, in vitro neurons KO for ATRX, TENM1, KCNQ2 and DLGAP2 show the highest L1 net expression and gene dysregulation.

In order to further explore this result, we performed gene DE analysis with DEseq2³⁹ by exploiting the technical replicates of KO iPSC and neuronal cell lines. The number of DE genes was plotted for all KOs (FDR < 0.05, |Log2FoldChange| > 0.5). While we observed transcriptional deregulation for all KOs in iPSCs, only KO samples for ATRX show a considerable number of DE genes also in neurons (Figure 2B). Given that the very same samples showed the strongest L1 RNA upregulation, consistently with what observed in the Z-score analysis, these results suggest that FL L1HS and L1PA elements might be associated to a pervasive gene deregulation in mature neurons of ASD subjects.

Upregulated L1 elements are expressed within introns of actively transcribed genes

Given the original observation that gene deregulation is generally stronger in post-mortem brain tissue and in vitro neurons showing L1 upregulation, we wanted to explore the genomic context of upregulated L1s and assess whether their transcription might be independent from the host transcripts. First, we analyzed ChIP-seq data from the NIH Roadmap epigenomics mapping consortium(53), focusing our attention on six major histone modifications (H3K4me1, H3K4me3, H3K9me3, H3K36me3, H3K27me3 and H3K27ac) from post-mortem mid frontal lobe of a control individual. Upon retrieving ChIP-seq data in wig format, and converting them to bigwig using USCS wigToBigWig utility, we computed a profile of histone mark enrichment on a bed file containing the genomic coordinates of a unique list of upregulated L1s in all samples and of all FL L1HS/L1PA annotated in the human genome separately. The profile was computed with DeepTools(51)⁠ computeMatrix utility and plotted with plotProfiles, as described in methods. Overall, all histone marks profiles show a substantial drop at the level of L1 genomic coordinates, which may be due to the typically lower mappability of repeated regions. However, upregulated L1s show a slightly weaker enrichment for transcriptional repressive marks H3K9me3 and H3K27me3 compared to total FL L1HS/L1PA, both within the boundaries of the L1 fragment and their flanking upstream and downstream regions. Interestingly, upregulated L1s exhibit a considerably stronger enrichment for activator marks (H3K36me3, H3K4me1 and H3K27ac) at the level of their flanking regions (Figure 3A).

Furthermore, we compared the percentage of upregulated, expressed and total FL L1PA/L1HS which overlap expressed genes. In both datasets the percentage of upregulated and expressed L1s found within expressed genes greatly exceeds the percentage of total L1s. However, upregulated L1s are neither enriched nor depleted within expressed genes compared to non-deregulated and expressed L1s (Supplementary figure 2A). Subsequently, we overlapped the genomic coordinates of upregulated L1s with the genomic coordinates of exons and introns retrieved from UCSC. All L1s not overlapping with either introns or exons were considered intergenic. For the vast majority of samples in all datasets, the number of upregulated intronic L1s greatly exceeds the number of exonic and intergenic L1s, especially in the samples with the highest number of upregulated L1s (Figure 3B). These results suggest that upregulated L1s may be enriched within introns of actively transcribed genes.

Given that most upregulated L1s are intronic, we wondered whether their upregulation might be due to an increased intron retention. We therefore performed a differential intron retention analysis focused on the groups of samples characterized by L1 upregulation with IRFinder(54). While we detected only a few differentially retained introns in iPSC and neurons KO for ATRX, we detected 5 downregulated and 806 upregulated introns in the group comprising samples SRR292614, SRR292620 and SRR292621 (Figure 3C). However, only one of these overlapped an upregulated L1, suggesting that intron retention is not at the basis of the observed L1s upregulation. Subsequently, we wanted to test whether intronic upregulated L1s would be spliced into novel transcripts together with nearby annotated exons. To this aim we counted, for all groups of samples characterized by L1 upregulation, the percentage of mapped reads whose fragments are shared by pairs of genomic segments comprising: a) a random expressed exons and its closest expressed exon, b) an upregulated L1 and its closest exon, c) a random 6,000 bp long sequence sampled from introns of expressed genes and its closest exon. The distribution of the percentage of shared reads between upregulated L1s and their closest exon is similar to that of the random intronic segments, and significantly lower than the average percentage of shared reads between couples of neighboring exons (Figure 3D). Taken together, these results suggest that upregulated L1s are mostly located within introns of expressed genes and likely expressed as independent transcriptional units.

Upregulated L1 elements may negatively impact the expression of specific ASD-related host genes

Starting from the observation that intronic upregulated L1s are mostly expressed as independent transcriptional units, we wanted to assess how they relate with the expression of the host genes. We limited our analysis to the samples with the highest L1 upregulation: samples SRR292614, SRR292620 and SRR292621 from post-mortem ACC (Velmeshev et al), and iPSC and neurons KO for ATRX. First, we counted the number of upregulated L1s found within the genomic loci of DE genes for each sample. We considered as DE those genes and L1s featuring a |Z-score| above 3. Overall, most upregulated L1s are not found within DE genes (Figure 4A). However, focusing on L1s mapping within DE genes, we found that 5 out of 11 upregulated L1s are located inside down-regulated genes in post-mortem ACC (Figure 4B, Table 1), suggesting that L1 upregulation might have a negative impact on the expression of a subset of specific host genes in the brain. To further explore this result, we selected all intronic upregulated L1s (Z-score > 3) in each of the three ACC sample displaying strong L1 upregulation (SRR292614, SRR292620, SRR292621) separately. For each of the three samples, we then associated the Z-scores of each intronic upregulated L1 with the Z-score of their overlapping genes. In all the three samples, Z-score associated to genes overlapping upregulated L1s are mostly negative (Additional file 4, Figure 4C), while genes overlapping non-upregulated L1s (|Z-score| < 2) do not show either a preferentially positive or negative Z-score (Additional file 4, Figure 4C). This result suggests that L1 upregulation might have a negative impact on the transcription of host genes for a subset of them.

A recent publication(55)⁠ revealed that, in mouse ESCs, L1-enriched genes are transcriptionally silenced by a not yet well understood mechanism which involves the transcription of the L1 element and the subsequent binding of the L1 RNA to the DNA at the level of specific loci. Indeed, through chromatin isolation by RNA purification (ChIRP) experiments, the authors highlighted ~25,000 genomic loci where the L1 RNA binds genomic DNA(55). Around 2,400 genes are found in proximity to these genomic loci (distance < 5 kbp). Interestingly, compared to a random distribution computed with equally sized sets of expressed genes, genes overlapping upregulated L1s and associated to a negative Z-score are enriched within the human orthologs of genes nearby L1-ChIRP peaks in all three ACC samples characterized by a strong L1 upregulation (SRR292614, SRR292620, SRR292621) (Figure 4D, Z-score > 5). This result reinforces the idea that the L1 upregulation occurring in some ASD samples might lead to the dysregulation of a physiological mechanism repressing the expression of overlapping intronic L1 by the host gene.

We then wanted to test whether genes associated to a negative Z-score, but overlapping upregulated L1s (negatively correlated genes), may be relevant for ASD. For each of the three tested samples (SRR292614, SRR292620, SRR292621), negatively correlated genes are enriched for SFARI genes, compared to distributions computed with equally sized sets of randomly expressed genes (Figure 4E). Furthermore, these genes are not significantly enriched for those highly expressed in adult brain(56) (Figure 4E). Interestingly, most SFARI negatively correlated genes are not unique to one of the three samples analysed, since 5 genes (1/3 of the total) are shared among all three samples (CADM2, CSMD1, DLG2, DLGAP1, GPHN) (Supplementary Figure 3A). Notably, the main function of these genes lies in synapse organization and function (Table 2). This evidence suggests that genes whose expression is potentially altered by L1 upregulation have similar roles crucial for neuronal function. One example of an L1 upregulated in all 3 samples, found within a SFARI gene (DLGAP1) and associated to a negative Z-score is displayed in figure 4F. These results may point to a subset of genes expressed in adult brain whose expression is negatively influenced by the expression of the overlapping intronic FL L1s.

To validate this expression pattern in an independent dataset, we exploited data from a recent study(57), where the authors calculated the coefficient of correlation between the expression of intronic and inter-genic L1s and the expression of the overlapping genes in 49 human tissues. The authors stratified L1s into elements expressed in all tissues and in those with a particularly high expression (>1.5 fold higher) in one tissue compared to all others. To ease data representation, we grouped the 49 tissues into 11 classes showing that the average coefficient of correlation between intronic tissue-specific L1s and overlapping genes was substantially lower for brain-derived tissues compared to all other tissues (Figure 4G). This adds further support to the existence of a transcriptional inverse correlation between intronic L1 and the host genes in healthy neurons for a subset of genes.

Taken together, these results suggest that L1 upregulation might, in some instances, lead to the downregulation of host genes. This could happen more frequently in genes important for ASD, with a mechanism not yet well understood.

ASD is a highly heterogeneous group of neurodevelopmental disorders. The identification of common molecular targets is therefore instrumental in defining homogeneous groups of individuals affected by ASD for clinical diagnosis and personalized medicine. Several works determined that transcriptional deregulation affecting both coding and non-coding gene expression occurs in ASD(14–18). However, the transcription of TEs has often been overlooked. Here we devoted special attention to the pattern of expression of evolutionarily young FL L1s, since they are pervasively transcribed and are both a controller and controlled by the epigenetic status of a cell(32,36,58,59). The main aim of this work was to assess whether L1 expression is altered in post-mortem ASD brains, in an in vitro model of iPSC and neurons KOs for several genes known to be directly involved in the etiology of ASD, and in the blood of discordant siblings. Moreover, we aimed at evaluating the impact of L1 expression on the transcription of protein coding genes.

Firstly, all ASD/KO samples show a moderate positive net number of upregulated L1s. However, a massive increase was evident in three samples of post-mortem ASD ACC and in iPSC and neurons KO for ATRX. These samples show instances of pervasive L1 upregulation, with 30-50% of analyzed L1s showing significantly higher expression levels compared to controls. ATRX is a SFARI level 1 gene encoding for a protein which contains an ATPase/helicase domain belonging to the SWI/SNF family of chromatin remodeling proteins(49,50). ATRX KO cells were previously shown to present increased chromatin accessibility at genomic loci occupied by retrotransposons(49). By the analysis of public ChIP-seq data for ATRX in human cell lines, we show a strong enrichment of upregulated L1HS/L1PA sequences in both iPSC and neurons, confirming that ATRX loss-of-function may directly lead to an increased transcription of young FL L1 elements.

Large scale genome investigations have contributed to the identification of almost one thousand genes putatively involved in ASD(11). These genes may be arbitrarily divided into two large functional groups: a) genes which exert a crucial role in synaptic function and b) genes involved in transcription regulation and/or chromatin remodeling, including ATRX. It is therefore tempting to speculate that mutations in a subset of the latter may directly influence L1 transcriptional regulation, since they are better positioned to have a genome-wide impact on the epigenetic status of chromatin and therefore a widespread effect on the transcriptional landscape of cells. To identify genes altered in ASD that may account for the strong L1 upregulation in post-mortem brains, we analyzed WES data relative to the Velmeshev et al dataset but could not identify candidate variants in known regulators of L1 transcription. However, further analysis is required to fully assess this intriguing hypothesis.

Upregulation seems to occur in a cell-type specific manner since individuals characterized by the strongest L1 increase in ACC do not present the same pattern of expression in FPC. Furthermore, no changes are observed in the blood of ASD subjects compared with their healthy siblings. However, this dataset is characterized by a much lower coverage (~10-15 million 50 bp reads per sample) compared to RNA-seq data retrieved from Velmeshev et al and Deneault et al (~100-150 150 bp million reads per sample). Further studies are therefore needed to confront this important issue. It is interesting to note that despite all the KO of ten ASD-related genes cause a certain extent of gene deregulation in iPSC, only neurons KO for ATRX showed a high number of DE genes, the only neuronal sample with strong L1 upregulation. This result is consistent with the observation that brain samples characterized by the strongest L1 upregulation are also among the ones with the highest number of DE genes. An increased level of L1 RNA expression in the brain seems thus a biological marker which can be associated only to a subset of ASD cases characterized by a large number of DE genes. These results are in line with the recent work by Wong et al(16)^⁠that suggested the existence of two major subgroups of ASDs. While the first subgroup recapitulated all known molecular changes typical of ASD(16), the second one was indistinguishable from control samples in terms of transcriptional and epigenetic alterations.

Most upregulated L1s are intronic and seem to be transcribed independently from their host transcripts. In the three post-mortem ACC samples characterized by L1 upregulation, a portion of upregulated L1s were hosted in downregulated protein coding genes: a single gene/L1 pair presented an upregulation of L1 RNA hosted within a downregulated mRNA of the protein coding gene. The analysis of an independent dataset(57) confirmed the existence of this inverse L1/gene expression relationship and that this pattern was a feature of genes expressed in the brain and with neuronal function. Importantly, they are enriched in SFARI genes further suggesting a causative relationship between the deregulation of specific intronic L1 elements and the establishment of neurodevelopmental phenotypes in ASD. Several models can be proposed on how L1s can regulate the expression of their host genes, influence the differentiation and homeostasis of neurons and be a crucial player in triggering neuronal dysfunction. The function of TEs as regulatory non-coding RNAs is currently under intense investigation. The identification of a large number of DNA/RNA hybrids at L1s loci suggest that TEs might exert their function in cis(55). They can also act by organizing locally chromatin domains within the same locus or by recruiting other sequences belonging to different chromosomes. Being capable of interacting with different proteins, they can recruit complexes to specific regions of the genomes. In this context, a recent study showed that L1 RNAs were functionally crucial for binding of Nucleolin-KAP1 complex to its target chromatin, allowing for ESC self-renewal and promoting rRNA synthesis(60). According to this model, expressed L1 sequences can control and can be controlled by the deposition of epigenetic marks and may promote the silencing and subsequent re-activation of specific sets of genes during development. Within a gene/L1 pair, L1 expression may give rise to transcriptional interference on their host gene or may guide the establishment of new epigenetic marks. Given the increasing evidence that epigenetic alteration occurs in ASD(16,20,21,43) ⁠and may directly occur at L1s genomic sequences, a direct link between epigenetic control of L1s RNA by regulatory genes mutated in some ASD cases and the expression of host genes may be hypothesized.

Our study focuses on L1s fragments longer than 5,000 nts, including FL L1s that maintain their capacity to retrotranspose. It remains therefore open the possibility that the increased expression of L1s gives rise to uncontrolled retrotransposition and therefore to somatic mosaicism in ASD brains(33,34,36). Recently, ASD brains have been shown to present somatic single nucleotide mutations(61,62). A deep analysis of genomic sequences in ASD post-mortem brains is therefore required to have a full understanding of the molecular consequences of L1s RNA upregulation. Our results also have relevant clinical implications. If deleterious mutations within ATRX or a defined set of genes with similar role in L1s transcriptional control are indeed at the basis of the molecular etiology of a subset of ASD subjects, WES-derived data may be used for personalized medicine. Effective medications for the treatment of ASD core symptoms are still lacking. Most of the drug currently in development for ASD derive from knowledge of genes implicated in monogenic disorders associated with altered neurodevelopmental trajectories and autistic symptoms such as Fragile X, Landau–Kleffner and Rett syndromes(63,64). As a consequence, the therapeutic approach to ASD subjects has traditionally focused on associated conditions(63), with poor impact for its core symptoms. These drugs typically target genes associated to synaptic pathways such as dopaminergic and glutamatergic receptors(63). Interestingly, a recent study revealed that treatment with a low dose of romidepsin restored social deficits in animal models of autism(65). Romidepsin inhibits the activity of the enzyme histone deacetylase(65), thus restoring the expression of genes involved in neuronal signaling and downregulated in ASD. Our results may provide a basis to stratify ASD cases for clinical treatments with drugs modifying the epigenetic status of cells or by interfering with L1 RNA expression. Recent observations on the potential therapeutic use of manipulating TEs in disease conditions in other tissue districts(66), may suggest that molecular tools to interfere with TEs expression could represent a new strategy for personalized treatment of neurodevelopmental disorders.

Limitations

Our analyses on post-mortem brain tissue are mostly limited by the lack of replicates among samples and the small sample size. This is the reason why our analysis should be considered exploratory and needs to be validated in a larger number of samples. However, we were able to detect characteristics specific to single samples by restricting the scope of our work to highly expressed L1 elements and by assisting our analyses on tissue-derived samples with in vitro models comprising biological replicates. The accurate quantification of TEs expression is a challenging task. This is mainly due to the repetitive nature of their sequence, which complicates the process of mapping sequencing reads to specific loci. Furthermore, TEs are often embedded in host transcripts; therefore, in most cases, they are unlikely expressed as independent transcriptional units. We aimed at overcoming these challenges by focusing our analyses on only full length evolutionarily young L1 elements. We detected a strong upregulation of L1 elements in a subset of ASD cases and in vitro models of ASD, which may impact the expression of ASD-relevant genes. However, we cannot be sure whether this molecular phenotype is strictly linked to pathogenic outcomes or is a result of underling alterations at the level of chromatin accessibility. Nevertheless, we believe that a pervasive increase in L1 RNA content is unlikely not to cause any negative impact on neurons vitality.

Datasets used

RNA-seq and WES datasets

We took advantage of public data provided by Velmeshelv et al(44)(PRJNA434002). This dataset comprises post-mortem whole-tissue bulk RNA-seq data from anterior cingulate cortex (ACC) and pre-frontal cortex (PFC) and blood-derived WES data of 15 ASD cases and 16 matched controls. Since RNA-seq was performed for both tissues for only 7 out of 15 ASD individuals, while for the remaining individuals only sequencing of either ACC or PFC was performed, the two datasets are only partially overlapping. Furthermore, whole exome sequencing (WES) data were available for the 15 ASD samples. Additionally, we retrieved RNA-seq data from induced pluripotent stem cell (iPSC) and neuron cells knockout (KO) for ten ASD-relevant genes (AFF2, ANOS1, ASTN2, ATRX, CACNA1C, CHD8, DLGAP2, KCNQ2, SCN2A, TENM1)(46). The dataset comprises a total of 86 RNA-seq samples, 42 from iPSCs and 44 from neurons. KOs were performed for each one of the ten genes through a CRISPR gene editing method to insert premature termination sites(46).

We analyzed RNA-seq data derived from blood of 73 pairs of discordant ASD siblings for a total of 146 samples, generated by Poly-A RNA sequencing on Illumina RNASeq Platform as previously described (67).

Chromatin markers

We retrieved ChIP-seq data in BigWig format from the NIH Roadmap epigenomics mapping consortium (GSM772833, GSM772834, GSM773012, GSM773013, GSM773014, GSM773015)(53). Retrieved data was relative to ChIP-seq experiments for six major histone modification (H3K4me1, H3K4me3, H3K9me3, H3K9ac, H3K27me3 and H3K27ac) from post-mortem mid frontal lobe of a healthy individual. We used this data to calculate the enrichment for different histonic markers within upregulated and total FL L1HS/L1PA.

ATRX ChIP-seq

ATRX ChIP-seq data in BigWig format was retrieved from a recent publication (PRJNA289924)(50)⁠. XL-MNase ChIP-seq of ATRX was performed in the erythroleukemic cell line K562, two independent replicates using different ATRX antibodies were performed. Retrieved data are relative to the enrichment of the two replicates over the input used as a control.

Variant discovery from whole exome sequencing data

Variant discovery was performed on whole exome sequencing (WES) data relative to ASD samples from Velmeshev et al(44)⁠. The exome sequencing raw data were used as input for an automatized pipeline with the purpose to identify and annotate all loci differing from the human reference genome, that we define as variants. The pipeline workflow could be summarized in three steps: 1) alignment, where reads from the sequencing are mapped to the reference genome; 2) variant discovery, where the alignments are examined to detect variations with respect to the reference human genome; 3) annotation, where a series of functional information are associated with each detected variation. In general, the GATK4 best practices(52) for the discovery of germinal short variants were followed. Alignment of raw paired-end reads to the reference genome (version hg19) was performed with bwa mem(68) (version 0.7.17) with standard parameters. Subsequently, duplicated reads were marked with picard (version 2.18.20). Variant discovery was then performed with the GATK4 utility HaplotypeCaller(52), using the appropriate file containing the coordinates of the sequences targeted by the exome sequencing. At the end of this phase a single vcf file was produced. Finally, all variants were annotated using annovar(69) (databases updated to 27/05/2019). ASD-relevant variants were defined as single-nucleotide variants included all the following characteristics: 1) Unannotated or annotated as pathogenic in the database of pathogenic variants ClinVar(70), 2) Minimum of 15 supporting reads, 3) Non-synonymous or splicing effect on the transcript, 4) CADD phred score > 20, 5) SIFT pred. = D, 6) Within a SFARI gene (score 1, 2, 3 or S). The SFARI gene database release 2022 Q3 was used(11).

TEs and gene quantification and differential expression

Raw RNA-seq reads were aligned to the reference genome (version hg19) with STAR(71)⁠ (version 2.6, standard parameters). Gene expression was quantified with HTseq-count (version 0.11.3; parameters: -f bam -r pos -t exon -i gene_name -a 10) on the BAM files obtained with STAR. The number of reads mapping to each genomic feature defined as ‘exon’ in the GTF file retrieved from gencode version v32lift37 were quantified. In order to quantify TEs expression, we used two different software: SQuIRE(47) and TEspeX(45). The two tools employ distinct strategies to assess TEs transcriptional activity. TEspeX quantifies TEs transcription levels by counting the number of reads mapping to a reference transcriptome built merging the RepBase(29) human TEs fasta sequences and the Ensembl(72) transcript sequences containing all the human coding and non-coding annotated transcripts using STAR (v2.6.0c). Only reads flagged as primary (-F 0 × 100 parameter) and not mapping to either coding or non-coding transcripts are counted using Python scripts and Picard FilterSamReads (v2.18.4). The amount of reads mapping to each TE subfamily fasta sequence is therefore obtained. On the other hand, SQuIRE (software for quantifying interspersed repeat expansion) allows for the quantification of TE upon spliced alignment of RNA-seq data to a reference genome. SQuIRE is a set of computational tools which use different strategies to combine counts from uniquely mapping and multi‐mapping reads and then generate counts for individual TE, as well as TEs by class and subfamily. Therefore using SQuIRE it has been possible to analyze transposable element expression on both subfamily and locus-specific (single copies of TEs annotated on RepeatMasker) level. SQuIRE (v0.9.9.92) was used. We performed an analysis which allowed us to estimate the amount of DE genes and of full-length L1HS/L1PA in each individual sample. Firstly, for each dataset, only L1HS/L1PA longer than 5,000 bp and featuring at least an average of 200 reads mapped and only protein-coding genes with at least 15 reads mapped were selected. Raw read counts were normalized with DEseq2(48). Then, for each ASD/KO sample, we calculated a Z-score for each gene and L1 by calculating the difference between each normalized expression value and the average normalized expression of control samples and dividing the result for the standard deviation of control samples. In the case of controls, for each dataset, we iteratively compared each control sample with the distribution of all other controls. Full-length L1HS/L1PA with a |Z-score| > 3 and genes with a |Z-score| > 5 were considered differentially expressed. In the case of the Deneault et al dataset, we also performed a differential expression analysis with DEseq2(48) (release 3.12). DEseq2 analysis was performed with standard parameters on the raw gene counts assessed by HTseq. Only protein-coding genes featuring an average of at least 200 mapped reads among the all the samples were considerate suitable for differential expression analysis. Genes with an FDR < 0.05 and LogFoldChange > 0.5 or < 0.5 were considered as differentially expressed.

ChIP-seq data enrichment profile

We computed enrichment profiles from ChIP-seq derived data using DeepTools(51). Data was retrieved in BigWig format. Enrichment profiles were computed for genomic segments of interest using DeepTools computeMatrix and plotProfile utilities.

Intron retention analysis

In order to estimate the rate of intron retention (IR) in groups of samples showing a strong upregulation of L1s, we used IRFinder(54). We built the reference for human genome GRCh37 (release69) and run IR quantification in fastq-mode, as outlined in the user manual. This mode outputs an end-to-end IR measurement upon aligning fastq reads to the built reference with STAR. Finally, differential intron retention analysis has been performed with the generalized linear model (GLM), as suggested for groups of samples with more than 3 replicates per condition. Introns have been considered differentially retained if featuring an FDR < 0.05.

L1 exonization analysis

We wanted to assess whether upregulated intronic L1 elements were exonized. We used a custom method to estimate the percentage of fragments of the same read pair mapping to two distinct genomic loci. Briefly, for each dataset, all reads mapping within the boundaries of a genomic locus were extracted from the bam file of all samples (only properly paired mapped reads, minimum mapping quality: 30). Then, the ID of these reads was used to extract all fragments belonging to these read pairs from the same bam files; these are the ‘total fragments’. Consequently, the coordinates of the total fragments were intersected with the genomic coordinates of a second genomic locus. Finally, the percentage of the read IDs found on both loci on the number of total fragments was calculated. The higher the percentage, the higher is the number of reads (or read pairs) spanning multiple regions of the same mRNA. We performed this analysis for three kinds of loci ‘couples’: a. 1,000 random expressed exons with their closest exon; b. upregulated FL L1HS/L1PA with their closest exons; c. 1,000 random 6,000 bp sequences randomly selected within the boundaries of expressed genes with their closest exons.

ASD: Autism Spectrum Disorder; TE: Transposable Element; FDR: False Discovery Rate; iPSC: Induced Pluripotent Stem Cell; ACC: Anterior Cingulate Cortex; PFC: Pre-Frontal Cortex; WES: Whole-Exome Sequencing; FL: Full-Length.

Acknowledgements

The authors would like to thank the guest editors and the anonymous reviewers for their constructive comments and all researchers who made the data used in this work publicy available, without it this study would have not been possible.

Availability of data and materials

Most of the analyses have been performed on publicly available data. Namely: Velmeshev et al. (PRJNA434002, DOI:10.1126/science.aav8130); Deneault et al. (PRJNA422099, DOI:10.1016/j.stemcr.2018.10.003); Valle-Garcìa et al. (PRJNA289924, DOI:10.1080/15592294.2016.1169351) and NIH Roadmap epigenomics mapping consortium (GSM772833, GSM772834, GSM773012, GSM773013, GSM773014, GSM773015; DOI:10.1038/nbt1010-1045).

Author contributions

RS, GS, SG designed the study; GS performed all analyses and wrote the manuscript; RS and SG supervised the analyses and reviewed the manuscript; DM reviewed the manuscript; ED, MF reviewed the manuscript and provided blood-derived RNA-seq data used in the study; all the authors approved the final version of the manuscript.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Author details

¹Area of Neuroscience, Scuola Internazionale Superiore di Studi Avanzati (SISSA), Trieste, Italy. ²Central RNA Laboratory, Istituto Italiano di Tecnologia – IIT, Via Enrico Melen 83, Building B 16152, Genova, Italy. ³CMP³vda, via Lavoratori Vittime del Col Du Mont 28, Aosta, Italy. ⁴Department of Cellular, Computational and Integrative Biology (CIBIO), University of Trento, Trento (TN), Italy. ⁵Eurac Research, Institute for Biomedicine, Bolzano (BZ), Italy. ⁶Fondazione The Microsoft Research - University of Trento Centre for Computational and Systems Biology (COSBI), Rovereto (TN), Italy.

Competing interests

The authors declare no competing interests.

Funding

Publication costs are funded by Istituto Italiano di Tecnologia (IIT).

Zeidan J, Fombonne E, Scorah J, Ibrahim A, Durkin MS, Saxena S, et al. Global prevalence of autism: A systematic review update. Autism Res Off J Int Soc Autism Res. 2022 May;15(5):778–90.
Sandin S, Lichtenstein P, Kuja-Halkola R, Hultman C, Larsson H, Reichenberg A. The Heritability of Autism Spectrum Disorder. JAMA. 2017 Sep;318(12):1182–4.
Luo W, Zhang C, Jiang Y-H, Brouwer CR. Systematic reconstruction of autism biology from massive genetic mutation profiles. Sci Adv. 2018 Apr;4(4):e1701799.
Feliciano P, Zhou X, Astrovskaya I, Turner TN, Wang T, Brueggeman L, et al. Exome sequencing of 457 autism families recruited online provides evidence for autism risk genes. NPJ genomic Med. 2019;4:19.
Yu TW, Chahrour MH, Coulter ME, Jiralerspong S, Okamura-Ikeda K, Ataman B, et al. Using whole-exome sequencing to identify inherited causes of autism. Neuron. 2013 Jan;77(2):259–73.
Alonso-Gonzalez A, Rodriguez-Fontenla C, Carracedo A. De novo Mutations (DNMs) in Autism Spectrum Disorder (ASD): Pathway and Network Analysis. Front Genet. 2018;9:406.
Codina-Solà M, Rodríguez-Santiago B, Homs A, Santoyo J, Rigau M, Aznar-Laín G, et al. Integrated analysis of whole-exome sequencing and transcriptome profiling in males with autism spectrum disorders. Mol Autism. 2015;6:21.
Satterstrom FK, Kosmicki JA, Wang J, Breen MS, De Rubeis S, An J-Y, et al. Large-Scale Exome Sequencing Study Implicates Both Developmental and Functional Changes in the Neurobiology of Autism. Cell. 2020 Feb;180(3):568-584.e23.
Fu JM, Satterstrom FK, Peng M, Brand H, Collins RL, Dong S, et al. Rare coding variation provides insight into the genetic architecture and phenotypic context of autism. Nat Genet. 2022 Sep;54(9):1320–31.
Sanders SJ. Next-Generation Sequencing in Autism Spectrum Disorder. Cold Spring Harb Perspect Med. 2019 Aug;9(8).
Banerjee-Basu S, Packer A. SFARI Gene: an evolving database for the autism research community. Vol. 3, Disease models & mechanisms. England; 2010. p. 133–5.
Devlin B, Scherer SW. Genetic architecture in autism spectrum disorder. Curr Opin Genet Dev. 2012 Jun;22(3):229–37.
Lasalle JM. Autism genes keep turning up chromatin. OA autism. 2013 Jun;1(2):14.
Parikshak NN, Swarup V, Belgard TG, Irimia M, Ramaswami G, Gandal MJ, et al. Genome-wide changes in lncRNA, splicing, and regional gene expression patterns in autism. Nature. 2016 Dec;540(7633):423–7.
Parikshak NN, Luo R, Zhang A, Won H, Lowe JK, Chandran V, et al. Integrative functional genomic analyses implicate specific molecular pathways and circuits in autism. Cell. 2013 Nov;155(5):1008–21.
Wong CCY, Smith RG, Hannon E, Ramaswami G, Parikshak NN, Assary E, et al. Genome-wide DNA methylation profiling identifies convergent molecular signatures associated with idiopathic and syndromic autism in post-mortem human brain tissue. Hum Mol Genet. 2019 Jul;28(13):2201–11.
Gandal MJ, Zhang P, Hadjimichael E, Walker RL, Chen C, Liu S, et al. Transcriptome-wide isoform-level dysregulation in ASD, schizophrenia, and bipolar disorder. Science. 2018 Dec;362(6420).
Gokoolparsadh A, Sutton GJ, Charamko A, Green NFO, Pardy CJ, Voineagu I. Searching for convergent pathways in autism spectrum disorders: insights from human brain transcriptome studies. Cell Mol Life Sci. 2016 Dec;73(23):4517–30.
Nardone S, Sams DS, Reuveni E, Getselter D, Oron O, Karpuj M, et al. DNA methylation analysis of the autistic brain reveals multiple dysregulated biological pathways. Transl Psychiatry. 2014 Sep;4(9):e433.
Ladd-Acosta C, Hansen KD, Briem E, Fallin MD, Kaufmann WE, Feinberg AP. Common DNA methylation alterations in multiple brain regions in autism. Mol Psychiatry. 2014 Aug;19(8):862–71.
Corley MJ, Vargas-Maya N, Pang APS, Lum-Jones A, Li D, Khadka V, et al. Epigenetic Delay in the Neurodevelopmental Trajectory of DNA Methylation States in Autism Spectrum Disorders. Front Genet. 2019;10:907.
Castelijns B, Baak ML, Timpanaro IS, Wiggers CRM, Vermunt MW, Shang P, et al. Hominin-specific regulatory elements selectively emerged in oligodendrocytes and are disrupted in autism patients. Nat Commun. 2020 Jan;11(1):301.
McClintock B. The origin and behavior of mutable loci in maize. Proc Natl Acad Sci [Internet]. 1950;36(6):344–55. Available from: http://www.pnas.org/cgi/doi/10.1073/pnas.36.6.344
Kazazian HH. Mobile elements: drivers of genome evolution. Science [Internet]. 2004;303(5664):1626–32. Available from: http://www.sciencemag.org/content/303/5664/1626.full.pdf%5Cnhttp://www.ncbi.nlm.nih.gov/pubmed/15016989
Trizzino M, Park Y, Holsbach-beltrame M, Aracena K. Transposable elements are the primary source of novelty in primate gene regulation. Renome Res. 2017;
Trizzino M, Kapusta A, Brown CD. Transposable elements generate regulatory novelty in a tissue-specific fashion. BMC Genomics. 2018;19(1):1–12.
Kim SY, Pritchard JK. Adaptive evolution of conserved noncoding elements in mammals. PLoS Genet. 2007;3(9):1572–86.
Eric M. Ostertag andHaig H. Kazazian Jr. Biology of mammalian L1 retrotransposons. 2001;
Bao W, Kojima KK, Kohany O. Repbase Update, a database of repetitive elements in eukaryotic genomes. Mob DNA. 2015;6:11.
Beck CR, Collier P, Macfarlane C, Malig M, Kidd JM, Eichler EE, et al. LINE-1 retrotransposition activity in human genomes. Cell. 2010 Jun;141(7):1159–70.
Ozata DM, Gainetdinov I, Zoch A, O’Carroll D, Zamore PD. PIWI-interacting RNAs: small RNAs with big functions. Nat Rev Genet. 2019 Feb;20(2):89–108.
Garcia-Perez JL, Widmann TJ, Adams IR. The impact of transposable elements on mammalian development. Development [Internet]. 2016;143(22):4101–14. Available from: http://dev.biologists.org/lookup/doi/10.1242/dev.132639
Muotri AR, Chu VT, Marchetto MCN, Deng W, Moran J V, Gage FH. Somatic mosaicism in neuronal precursor cells mediated by L1 retrotransposition. Nature. 2005 Jun;435(7044):903–10.
Baillie JK, Barnett MW, Upton KR, Gerhardt DJ, Richmond TA, De Sapio F, et al. Somatic retrotransposition alters the genetic landscape of the human brain. Nature. 2011 Oct;479(7374):534–7.
Terry DM, Devine SE. Aberrantly High Levels of Somatic LINE-1 Expression and Retrotransposition in Human Neurological Disorders. Front Genet. 2019;10:1244.
Coufal NG, Garcia-Perez JL, Peng GE, Yeo GW, Mu Y, Lovci MT, et al. L1 retrotransposition in human neural progenitor cells. Nature. 2009 Aug;460(7259):1127–31.
Zhao B, Wu Q, Ye AY, Guo J, Zheng X, Yang X, et al. Somatic LINE-1 retrotransposition in cortical neurons and non-brain tissues of Rett patients and healthy individuals. PLoS Genet. 2019 Apr;15(4):e1008043.
Faulkner GJ, Kimura Y, Daub CO, Wani S, Plessy C, Irvine KM, et al. The regulated retrotransposon transcriptome of mammalian cells. Nat Genet. 2009 May;41(5):563–71.
Honson DD, Macfarlan TS. A lncRNA-like Role for LINE1s in Development. Dev Cell. 2018 Jul;46(2):132–4.
Wang L, Dou K, Moon S, Tan FJ, Zhang ZZ. Hijacking Oogenesis Enables Massive Propagation of LINE and Retroviral Transposons. Cell. 2018 Aug;174(5):1082-1094.e12.
Jachowicz JW, Bing X, Pontabry J, Bošković A, Rando OJ, Torres-Padilla M-E. LINE-1 activation after fertilization regulates global chromatin accessibility in the early mouse embryo. Nat Genet. 2017 Oct;49(10):1502–10.
Shpyleva S, Melnyk S, Pavliv O, Pogribny I, Jill James S. Overexpression of LINE-1 Retrotransposons in Autism Brain. Mol Neurobiol. 2018 Feb;55(2):1740–9.
Tangsuwansri C, Saeliw T, Thongkorn S, Chonchaiya W, Suphapeetiporn K, Mutirangura A, et al. Investigation of epigenetic regulatory networks associated with autism spectrum disorder (ASD) by integrated global LINE-1 methylation and gene expression profiling analyses. PLoS One. 2018;13(7):e0201071.
Velmeshev D, Schirmer L, Jung D, Haeussler M, Perez Y, Mayer S, et al. Single-cell genomics identifies cell type-specific molecular changes in autism. Science. 2019 May;364(6441):685–9.
Ansaloni F, Gualandi N, Esposito M, Gustincich S, Sanges R. TEspeX: consensus-specific quantification of transposable element expression preventing biases from exonized fragments. Bioinformatics. 2022 Jul;
Deneault E, White SH, Rodrigues DC, Ross PJ, Faheem M, Zaslavsky K, et al. Complete Disruption of Autism-Susceptibility Genes by Gene Editing Predominantly Reduces Functional Connectivity of Isogenic Human Neurons. Stem cell reports. 2018 Nov;11(5):1211–25.
Yang WR, Ardeljan D, Pacyna CN, Payer LM, Burns KH. SQuIRE reveals locus-specific regulation of interspersed repeat expression. Nucleic Acids Res. 2019 Mar;47(5):e27.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.
Sadic D, Schmidt K, Groh S, Kondofersky I, Ellwart J, Fuchs C, et al. Atrx promotes heterochromatin formation at retrotransposons. EMBO Rep. 2015 Jul;16(7):836–50.
Valle-García D, Qadeer ZA, McHugh DS, Ghiraldini FG, Chowdhury AH, Hasson D, et al. ATRX binds to atypical chromatin domains at the 3’ exons of zinc finger genes to preserve H3K9me3 enrichment. Epigenetics. 2016 Jun;11(6):398–414.
Ramírez F, Dündar F, Diehl S, Grüning BA, Manke T. deepTools: a flexible platform for exploring deep-sequencing data. Nucleic Acids Res. 2014 Jul;42(Web Server issue):W187-91.
Van der Auwera GA, Carneiro MO, Hartl C, Poplin R, Del Angel G, Levy-Moonshine A, et al. From FastQ data to high confidence variant calls: the Genome Analysis Toolkit best practices pipeline. Curr Protoc Bioinforma. 2013;43(1110):11.10.1-11.10.33.
Bernstein BE, Stamatoyannopoulos JA, Costello JF, Ren B, Milosavljevic A, Meissner A, et al. The NIH Roadmap Epigenomics Mapping Consortium. Nat Biotechnol. 2010 Oct;28(10):1045–8.
Middleton R, Gao D, Thomas A, Singh B, Au A, Wong JJ-L, et al. IRFinder: assessing the impact of intron retention on mammalian gene expression. Genome Biol. 2017 Mar;18(1):51.
Lu JY, Shao W, Chang L, Yin Y, Li T, Zhang H, et al. Genomic Repeats Categorize Genes with Distinct Functions for Orchestrated Regulation. Cell Rep. 2020 Mar;30(10):3296-3311.e5.
Leblond CS, Le T-L, Malesys S, Cliquet F, Tabet A-C, Delorme R, et al. Operative list of genes associated with autism and neurodevelopmental disorders based on database review. Mol Cell Neurosci. 2021 Jun;113:103623.
Bogu GK, Reverter F, Marti-Renom MA, Snyder MP, Guigó R. Atlas of transcriptionally active transposable elements in human adult tissues. bioRxiv [Internet]. 2019 Jan 1;714212. Available from: http://biorxiv.org/content/early/2019/07/24/714212.abstract
Rangwala SH, Zhang L, Kazazian HH. Many LINE1 elements contribute to the transcriptome of human somatic cells. Genome Biol [Internet]. 2009;10(9):R100. Available from: http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=2768975&tool=pmcentrez&rendertype=abstract
Goodier JL. Restricting retrotransposons: A review. Mob DNA [Internet]. 2016;7(1). Available from: http://dx.doi.org/10.1186/s13100-016-0070-z
Percharde M, Lin C-J, Yin Y, Guan J, Peixoto GA, Bulut-Karslioglu A, et al. A LINE1-Nucleolin Partnership Regulates Early Development and ESC Identity. Cell. 2018 Jul;174(2):391-405.e19.
Rodin RE, Dou Y, Kwon M, Sherman MA, D’Gama AM, Doan RN, et al. The landscape of somatic mutation in cerebral cortex of autistic and neurotypical individuals revealed by ultra-deep whole-genome sequencing. Nat Neurosci. 2021 Feb;24(2):176–85.
Bae T, Fasching L, Wang Y, Shin JH, Suvakov M, Jang Y, et al. Somatic mutations reveal hypermutable brains and are associated with neuropsychiatric disorders. medRxiv [Internet]. 2022; Available from: https://www.medrxiv.org/content/early/2022/02/10/2022.02.09.22269114
Ghosh A, Michalon A, Lindemann L, Fontoura P, Santarelli L. Drug discovery for autism spectrum disorder: challenges and opportunities. Nat Rev Drug Discov. 2013 Oct;12(10):777–90.
Sung M, Chin CH, Lim CG, Liew HSA, Lim CS, Kashala E, et al. What’s in the pipeline? Drugs in development for autism spectrum disorder. Neuropsychiatr Dis Treat. 2014;10:371–81.
Qin L, Ma K, Wang ZJ, Hu Z, Matas E, Wei J, et al. Social deficits in Shank3-deficient mouse models of autism are rescued by histone deacetylase (HDAC) inhibition. Nat Neurosci [Internet]. 2018;21(4):564–75. Available from: http://dx.doi.org/10.1038/s41593-018-0110-8
Marasca F, Sinha S, Vadalà R, Polimeni B, Ranzani V, Paraboschi EM, et al. LINE1 are spliced in non-canonical transcript variants to regulate T cell quiescence and exhaustion. Nat Genet. 2022 Feb;54(2):180–93.
Filosi M, Kam-Thong T, Essioux L, Muglia P, Trabetti E, Spooren W, et al. Transcriptome signatures from discordant sibling pairs reveal changes in peripheral blood immune cell composition in Autism Spectrum Disorder. Transl Psychiatry. 2020 Apr;10(1):106.
Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics. 2009;25(14):1754–60.
Wang K, Li M, Hakonarson H. ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data. Nucleic Acids Res. 2010 Sep;38(16):e164.
Landrum MJ, Lee JM, Riley GR, Jang W, Rubinstein WS, Church DM, et al. ClinVar: public archive of relationships among sequence variation and human phenotype. Nucleic Acids Res. 2014 Jan;42(Database issue):D980-5.
Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013 Jan;29(1):15–21.
Zerbino DR, Achuthan P, Akanni W, Amode MR, Barrell D, Bhai J, et al. Ensembl 2018. Nucleic Acids Res. 2018;46(D1):D754–61.

Table 1 – Differentially expressed genes overlapping upregulated L1s in specific samples.

L1 chr	L1 start	L1 end	L1 ID	Gene name	Gene alteration	Sample ID
chr11	14799221	14805319	L1PA5:L1:LINE	PDE3B	Up-regulated	SRR9292614
chr1	68553278	68559388	L1PA6:L1:LINE	WLS	Up-regulated	SRR9292620
chr11	14799221	14805319	L1PA5:L1:LINE	PDE3B	Up-regulated	SRR9292620
chr11	132190054	132196077	L1PA4:L1:LINE	NTM	Up-regulated	SRR9292620
chr5	139895099	139900527	L1PA3:L1:LINE	ANKHD1-EIF4EBP3	Up-regulated	SRR9292620
chr1	220814954	220821404	L1PA10:L1:LINE	MARK1	Down-regulated	SRR9292620
chr10	49580263	49585413	L1PA7:L1:LINE	MAPK8	Down-regulated	SRR9292620
chr11	132626675	132633167	L1PA7:L1:LINE	OPCML	Down-regulated	SRR9292620
chr18	3695243	3700586	L1PA7:L1:LINE	DLGAP1	Down-regulated	SRR9292620
chr19	40543822	40549961	L1PA4:L1:LINE	ZNF780B	Down-regulated	SRR9292620
chr5	139895099	139900527	L1PA3:L1:LINE	ANKHD1-EIF4EBP3	Up-regulated	SRR9292621

Table 2 – Description of the 5 genes associated to a negative Z-score and overlapping an upregulated L1 element in all the three samples characterized by pervasive L1 upregulation.

Gene name	SFARI score	Gene description (SFARI gene, release 2022 Q3)
CADM2	2	Adhesion molecule that engages in homo- and heterophilic interactions with the other nectin-like family members, leading to cell aggregation. Important for synapse organization, providing regulated trans-synaptic adhesion.
CSMD1	2	Weakly expressed in most tissues, except in brain (expressed at intermediate levels in brain, including cerebellum, substantia nigra, hippocampus, and fetal brain). Variants in this gene have been shown to associate with schizophrenia and bipolar disorder.
DLG2	2	This gene encodes a member of the membrane-associated guanylate kinase (MAGUK) family. The encoded protein forms a heterodimer with a related family member that may interact at postsynaptic sites to form a multimeric scaffold for the clustering of receptors, ion channels, and associated signaling proteins.
DLGAP1	2	The protein encoded by this gene is expressed in the brain, localizes to the postsynaptic density, and interacts with a number of ASD-associated proteins, including DLG1, DLG4, SHANK1, SHANK2 and SHANK3.
GPHN	2	This gene encodes a neuronal assembly protein that anchors inhibitory neurotransmitter receptors to the postsynaptic cytoskeleton via high affinity binding to a receptor subunit domain and tubulin dimers. In nonneuronal tissues, the encoded protein is also required for molybdenum cofactor biosynthesis. Mutations in this gene may be associated with the neurological condition hyperplexia and also lead to molybdenum cofactor deficiency.

No competing interests reported.

AdditionalFile1.xls
Additional file 1 (.xls): Normalized expression values for each TE subfamily in each sample of the Velmeshev et al and Deneault et al datasets, quantified with TEspeX.
AdditionalFile2.xls
Additional file 2 (.xls): List of FL L1s expressed in the Velmeshev et al and Deneault et al datasets
AdditionalFile3.xls
Additional file 3 (.xls): Potentially deleterious mutations identified within SFARI genes in all individuals of the Velmeshev et al dataset.
AdditionalFile4.xls
Additional file 4 (.xls): List of all FL L1 associated to a Z-score >3 in at least one of the samples of the Velmeshev et al - ACC dataset characterized by strong L1 upregulation (SRR9292614, SRR9292620, SRR9292621) and Z-scores of the genes overlapping each upregulated L1s.
Supplementaryfigure1.pdf
Supplementaryfigure2.pdf
Supplementaryfigure3.pdf

Download PDF

Journal Publication

published 28 Jun, 2023

Read the published version in Molecular Autism →

Editorial decision: Major revision
10 Feb, 2023
Reviews received at journal
24 Jan, 2023
Reviews received at journal
21 Dec, 2022
Reviewers agreed at journal
28 Nov, 2022
Reviewers agreed at journal
23 Nov, 2022
Reviewers invited by journal
14 Nov, 2022
Editor assigned by journal
14 Nov, 2022
Submission checks completed at journal
06 Nov, 2022
First submitted to journal
03 Nov, 2022

You are reading this latest preprint version

Exploratory Analysis of L1 retrotransposons expression in Autism

Status:

Journal Publication

Version 1

Abstract

Figures

Introduction

Results

Discussion

Methods

Variant discovery from whole exome sequencing data

TEs and gene quantification and differential expression

L1 exonization analysis

Abbreviations

Declaration

References

Tables

Additional Declarations

Supplementary Files

Status:

Journal Publication

Version 1