Epigenetic Gene-Regulatory Loci in Alu Elements Associated with Autism Susceptibility in the Prefrontal Cortex of ASD

Alu elements are transposable elements that can influence gene regulation through several mechanisms; nevertheless, it remains unclear whether dysregulation of Alu elements contributes to the neuropathology of autism spectrum disorder (ASD). In this study, we characterized transposable element expression profiles and their sequence characteristics in the prefrontal cortex tissues of ASD and unaffected individuals using RNA-sequencing data. Our results showed that most of the differentially expressed transposable elements belong to the Alu family, with 659 loci of Alu elements corresponding to 456 differentially expressed genes in the prefrontal cortex of ASD individuals. We predicted cis- and trans-regulation of Alu elements to host/distant genes by conducting correlation analyses. The expression level of Alu elements correlated significantly with 133 host genes (cis-regulation, adjusted p < 0.05) associated with ASD as well as the cell survival and cell death of neuronal cells. Transcription factor binding sites in the promoter regions of differentially expressed Alu elements are conserved and associated with autism candidate genes, including RORA. COBRA analyses of postmortem brain tissues showed significant hypomethylation in global methylation analyses of Alu elements in ASD subphenotypes as well as DNA methylation of Alu elements located near the RNF-135 gene (p < 0.05). In addition, we found that neuronal cell density, which was significantly increased (p = 0.042), correlated with the expression of genes associated with Alu elements in the prefrontal cortex of ASD. Finally, we determined a relationship between these findings and the ASD severity (i.e., ADI-R scores) of individuals with ASD. Our findings provide a better understanding of the impact of Alu elements on gene regulation and molecular neuropathology in the brain tissues of ASD individuals, which deserves further investigation.


Introduction
Autism spectrum disorder (ASD) is a neurodevelopmental disorder characterized by behavioral deficits based on the Diagnostic and Statistical Manual of Mental Disorders, Fifth Edition (DSM-5) criteria, including deficits in social interactions and communication, restricted interests, and repetitive behaviors [1]. According to the Centers for Disease Control and Prevention (CDC), the prevalence of ASD increased dramatically in 2020, reaching approximately 1 in 36 children in the United States [2]. ASD is a complex condition, even though its root causes are still unknown. Despite substantial genetic evidence over the past two decades, only 10-20% of ASD cases involve an identified genetic cause, and more than half of ASD cases are idiopathic [3]. A combination of both genetic and environmental factors influences the heterogeneity of ASD's clinical phenotypes and severity [4][5][6][7][8][9][10][11], suggesting that epigenetic mechanisms may be the link between these two factors. DNA methylation is a key epigenetic process that contributes to the etiology of neurodevelopmental disorders, including ASD [12].
Several studies have demonstrated that alterations in DNA methylation at transcription start sites (TSSs) in ASD brains compared to control brains result in the altered expression of many genes, including OXTR, SNRPN, MAGEL2, RELN, and GAD1 [12]. However, epigenetic changes in a single gene may not occur frequently, and genome-wide DNA methylation profiling in ASD is considered an alternative tool used for identifying differentially methylated loci contributing to ASD susceptibility. Genome-wide transcriptome and DNA methylome analyses in postmortem ASD brain tissue have revealed that differentially methylated and expressed genes are involved in synaptic development and immune function [13][14][15][16]. Nonetheless, most DNA methylation and transcriptome analyses of postmortem ASD brain tissues have focused on protein-coding regions rather than noncoding regions, including the transposable elements that account for most of the CpG sites in the human genome.
Transposable elements (TEs) are mobile genetic elements that copy and paste themselves into a new genomic location. Alu elements and long interspersed nuclear element-1 (LINE-1) comprise the majority of TEs that remain active in the human genome, comprising more than 30% of the genome at a copy number of over 1 million elements [17]. By introducing alternative promoters or enhancers, novel splicing sites, and epigenetic alterations, these elements affect the expression of host or nearby protein-coding genes [17]. Alu elements, which are 300 bp in length, are the most common short interspersed nuclear element (SINE). Alu has a dimeric structure that is derived from the 7SL RNA gene, with each monomer linked by an A-rich region [18]. The left monomer of Alu contains the A-and B-box promoters for RNA polymerase III promoters. The right monomer of Alu elements has a poly(A) tail that is important for their retrotransposition mechanism [19]. Alu is a nonautonomous retrotransposon. The retrotransposition mechanism of Alu is similar to that of LINE-1 retrotransposition but differs in that Alu is unable to produce ORF1p and ORF2p proteins. Alu RNA requires the ORFp protein produced from LINE-1 for transport to the nucleus and target-site-primed reverse transcription [17][18][19]. At present, understanding of the precise role of LINE-1 and Alu elements in ASD development remains limited. In previous studies, we reported alterations in the global DNA methylation of LINE-1 and Alu elements in lymphoblastoid cell lines (LCLs) of an ASD subgroup based on clinical phenotype (ADI-R scores). Alterations in LINE-1 methylation were specifically observed in the ASD subtype with severe language impairment, whereas an alteration in Alu methylation was specific to ASD with mild symptoms. Moreover, several genes with LINE-1 and Alu insertions were differentially expressed in the peripheral blood and LCLs of ASD. These genes are known to be associated with neurodevelopmental disorders [20,21]. Locus-and family-specific Alu and LINE-1 methylation and the expression level of their host genes were found to be altered in the blood of ASD individuals identified using the Infinium 450 K platform [22]. Global LINE-1 methylation was significantly decreased in the blood samples from ASD individuals with mental regression [23]. In addition, whole-genome sequencing analysis of ASD brain tissues shows a higher number of novel insertions of Figure 1. The conceptual framework of this study. This figure shows the structure and transcriptional regulation or retrotransposition of Alu elements through epigenetic mechanisms. Alu elements can act as enhancers or alternative promoters for host/target genes, resulting in transcriptome alteration because they contain several transcription factor binding sites [17,26,27]. It is unclear whether Alu elements are associated with the transcriptome, neuropathology, and clinical phenotype (ADI-R scores) of the brain tissues of ASD individuals. This figure was created with BioRender.com (accessed on 19 March 2023).

Transposable Element Profiling in the Prefrontal Cortex of ASD
To characterize the TE expression profile in the prefrontal cortex of ASD individuals, we performed a transcriptomic study using postmortem brain tissues from ASD and unaffected individuals, as described in the Methods. We selected one (Liu's study) from a total of four studies (Table S1) because it provides RNA-seq data on SRA deposits and has the largest sample size, including 21 males with ASD and 21 sex-and age-matched unaffected controls (female data were excluded in this study to remove sex bias). A total of 76,877 detectable TEs (excluding TEs with low expression) were used for differential expression analysis. This TE expression profiling revealed 864 differentially expressed TEs (DETEs) at the statistical threshold (adjusted p < 0.01, log2FC > 1) (Figure 2a, Table S2), including 218 downregulated and 646 upregulated DETEs. The distribution of upregulated DETEs per family showed that the majority of upregulated TEs belong to the Alu (48% or 308 elements) and LINE (24% or 154 elements) families (Figure 2b). Approximately 7% (or 47) and 21% (or 137) of the DETEs are members of the mammalian interspersed repeat (MIR) and other TE families, respectively. Moreover, the downregulated DETEs were distributed across other TE families, including LINE-1 (33%), MIRs (27%), Figure 1. The conceptual framework of this study. This figure shows the structure and transcriptional regulation or retrotransposition of Alu elements through epigenetic mechanisms. Alu elements can act as enhancers or alternative promoters for host/target genes, resulting in transcriptome alteration because they contain several transcription factor binding sites [17,26,27]. It is unclear whether Alu elements are associated with the transcriptome, neuropathology, and clinical phenotype (ADI-R scores) of the brain tissues of ASD individuals. This figure was created with BioRender.com (accessed on 19 March 2023).
In this study, we aimed to investigate relationships between TEs, the transcriptome, neuropathology, and clinical phenotype (ADI-R scores) in the brain tissues of individuals with ASD. We first collected RNA-seq datasets from postmortem brain tissues of individuals with ASD and those who were not affected from the NCBI Gene Expression Omnibus (GEO) DataSets and the NCBI Sequence Read Archive (SRA) repository to explore the expression profile of TEs that suggest TE dysregulation in ASD. The datasets were then reanalyzed for gene and TE profiling and the bedtools package was used to identify host/target genes among the intersected genomic loci between differentially expressed TEs and differentially expressed genes. To determine whether TEs can influence host/target gene expression, we performed correlation analysis between the expression value of each TE and genes by classifying TEs into those residing within exons, introns, and those found upstream from the TSS of genes. Next, to understand whether TEs are associated with ASD, we used Ingenuity Pathway Analysis (IPA) software to predict the biological functions and gene regulatory networks of host/target genes associated with TEs. Moreover, we analyzed in more detail the conserved regions and upstream regulators of dysregulated TEs in ASD. Due to the abundance of TEs in the genome, we collected DNA sequences of differentially expressed TEs (DETEs) in ASD brain tissues and compared them with sequences of nondifferentially expressed TEs using multiple sequence alignment analysis. The conserved regions were then employed to predict transcription factor binding sites (TFBSs) using the JASPAR CORE database. In addition, we conducted a pilot study using postmortem human brain tissues of ASD and unaffected individuals obtained from the Harvard Brain Tissue Resource Center and confirmed the expression levels of host/target genes selected from RNA-seq analysis. To determine DNA methylation of Alu elements, we applied combined bisulfite restriction analysis (COBRA) for global Alu methylation and locus-specific methylation of Alu located near host genes. Finally, Nissl staining was performed to evaluate the neuropathology in paraffin-embedded postmortem brain tissues. To determine the relationship between TEs, neuropathology, and ASD severity (ADI-R scores), we performed correlation analysis between the results. We expect that the results of our study will provide a better understanding of the relationship between TEs and the molecular and neuropathology of ASD.

Transposable Element Profiling in the Prefrontal Cortex of ASD
To characterize the TE expression profile in the prefrontal cortex of ASD individuals, we performed a transcriptomic study using postmortem brain tissues from ASD and unaffected individuals, as described in the Methods. We selected one (Liu's study) from a total of four studies (Table S1) because it provides RNA-seq data on SRA deposits and has the largest sample size, including 21 males with ASD and 21 sex-and age-matched unaffected controls (female data were excluded in this study to remove sex bias). A total of 76,877 detectable TEs (excluding TEs with low expression) were used for differential expression analysis. This TE expression profiling revealed 864 differentially expressed TEs (DETEs) at the statistical threshold (adjusted p < 0.01, log 2 FC > 1) (Figure 2a, Table S2), including 218 downregulated and 646 upregulated DETEs. The distribution of upregulated DETEs per family showed that the majority of upregulated TEs belong to the Alu (48% or 308 elements) and LINE (24% or 154 elements) families ( Figure 2b). Approximately 7% (or 47) and 21% (or 137) of the DETEs are members of the mammalian interspersed repeat (MIR) and other TE families, respectively. Moreover, the downregulated DETEs were distributed across other TE families, including LINE-1 (33%), MIRs (27%), and others (33%).
TEs are abundant in the human genome; thus, to ensure that the large number of DETEs of each family was not caused by this, we compared the number of DETEs of each family with the total number of elements in the family (Figure 2c). Interestingly, when compared with other TEs (LINE and MIR families), young families (AluY and AluS) exhibited more DE-Alu elements than the total number elements in the family. For the LINE family, most DE-LINEs belonged to the inactive family (L2 and L3), with a larger number of total elements than DE-Alu elements. These DETEs indicate epigenetic alteration at TE loci in the prefrontal cortex of ASD individuals. Moreover, only a small number of DE-LINE-1 (L1) elements, which is an active LINE family, were differentially expressed. However, the most active L1 subfamily (L1H) was not significantly expressed in the prefrontal cortex of ASD individuals in our analysis. In this study, we focused on the Alu family, which is a major DETE according to the number of DE-Alu elements and the total number of elements in the family. loci in the prefrontal cortex of ASD individuals. Moreover, only a small number of DE-LINE-1 (L1) elements, which is an active LINE family, were differentially expressed. However, the most active L1 subfamily (L1H) was not significantly expressed in the prefrontal cortex of ASD individuals in our analysis. In this study, we focused on the Alu family, which is a major DETE according to the number of DE-Alu elements and the total number of elements in the family.

Correlations between Alu Expression and Differentially Expressed Genes in the Prefrontal Cortex of ASD Individuals in Cis-and Trans-Regulatory Manners
Because TEs can provide enhancer elements, act as alternative promoters, and epigenetically regulate nearby and distant target genes, we performed correlation analysis between the expression of DETEs and differentially expressed genes (DEGs). We hypothesized that a dysregulation of TEs may reflect a change in an epigenetic mark that induces the transcriptome profile in the brains of ASD individuals. In this study, we performed Pearson's correlation analysis for cis-and trans-regulation. We also reanalyzed the same RNA-seq data for differential gene expression analysis. To determine whether DETEs regulate the expression level of genes located near DETEs, we overlapped the genomic locations of DETEs and RefSeq genes to identify targeted genes using bedtools prior to correlation analysis (Figure 3a,b). First, DETE loci were classified based on genomic features, including regions upstream of the TSS (50 kb, 10 kb, and 1 kb), intronic regions, and exonic regions. We found 18,045 DETE positions located within or near 1587 locations corresponding to 1493 DEGs in the prefrontal cortex of ASD individuals compared to unaffected controls (Table S3). Most DETEs are located in the 50 kb upstream region (8199 loci) and intronic region (7232 loci) of 507 and 744 DEGs, respectively (Table 1). Correlation results showed that 534 DETE loci correlated significantly with DEGs in the prefrontal cortex of ASD individuals (Benjamini-Hochberg (BH)-adjusted p < 0.05). As expected, most of the correlating DETEs are located in the 50 kb and intronic regions of the DEGs. When the correlation was performed separately for ASD and unaffected individuals, we observed a different number of correlating DEGs. In ASD, a high number of DEGs correlated with DETEs, especially in the Alu family. When focusing on Alu elements, we found a greater difference in the number of significantly correlating DEGs between ASD and unaffected controls for 659 loci corresponding to 456 DEGs (Table 1, Table S4). In ASD, 161 Alu elements correlated with DEGs; in unaffected individuals, 86 Alu elements correlated with DEGs. A correlation matrix was generated that showed the difference in the top 30 significant loci between ASD and unaffected individuals (Figure 3c-e). These findings indicate that DEGs in the prefrontal cortex of ASD individuals may be associated with dysregulation of the Alu family.
For trans-regulation analysis, we performed correlation analysis between all probabilities between 659 DE-Alu and 1493 DEGs. A~500 k correlation probability between DE-Alu was tested, and most of the DEGs showed a significant correlation with these elements (Table 2) based on Pearson's correlation analysis with multiple hypothesis correction (FDR = 0.05). We found 40,599 significant probabilities that Alu influences 1465 DEGs in the prefrontal cortex of ASD individuals and 7141 probabilities for 1239 DEGs in unaffected individuals. These findings suggested that Alu elements contribute to the trans-regulation of several genes in the genome. This was expected, because TEs contribute to normal neurodevelopment and it is possible that these elements control other genes in the genome in addition to the DEGs.
To fully comprehend how TEs contribute to ASD, it may not be adequate to focus only on the number of DEGs involved in trans-regulation; the abundance of TE components in ASD and unaffected individuals is also interesting. Alu elements were found in higher abundance in ASD individuals than in unaffected controls, even though the results showed a similar number of DEGs in trans-regulation (316 and 196 elements, respectively). We thus overlapped significant probability loci between ASD and unaffected individuals ( Figure S1), with 191 Alu elements found in both ASD and unaffected individuals (only 5 Alu elements from among the 196 found were exclusively in unaffected controls). These elements might be normal Alu elements that contribute to normal processes in the genome. Additionally, we discovered 125 Alu elements to be substantially associated in individuals with ASD but not in unaffected controls. These findings suggest that trans-regulatory elements from Alu are increased in the prefrontal cortex of ASD, leading to the dysregulation of the related genes.    To determine whether DEGs in the prefrontal cortex of ASD individuals are associated with ASD and neuropathology, we obtained a list of significantly correlating DEGs for biological function and pathway analyses using IPA software (161 Alu-correlating loci corresponding to 133 DEGs in Table 1). The results showed that Alu-correlating DEGs were enriched for diseases and functions associated with ASD (16 genes, p = 1.19 × 10 −3 ) and comorbid disorders, including pervasive developmental disorder, mental retardation, movement disorders, and epilepsy ( Figure 4a, Table S5). We also found that Alu-correlated DEGs are related to apoptosis (30 genes, p = 5.75 × 10 −4 ), and the gene regulatory network predicted by IPA also highlighted the survival of neural cells and the cell death of cerebral cortex cells (Figure 4c). Furthermore, IPA results showed the association of Alu-correlated DEGs with canonical pathways such as the synaptogenesis signaling pathway, neuroinflammation signaling pathway, cell cycle, and glutamate receptor signaling (p < 0.05) (Figure 4b, Table S6). These findings indicate that in the prefrontal cortex of ASD individuals, a dysregulation of Alu retrotransposons affects the expression of host/nearby protein-coding genes associated with the neuropathology implicated in ASD. among autism candidate genes (p = 3.67 × 10 −6 ); the list of overlapping genes associated with Alu elements is shown in Figure S2b. This finding suggests that DE-Alu influences autism gene expression in the prefrontal cortex, which may lead to risk of ASD.

Consensus Sequence and Transcription Factor Binding Sites at DE-Alu Elements
Although Alu elements are abundant in the genome, only a small portion of Alu elements are active and escape from restriction by an epigenetic mechanism. To determine In addition to the predicted biological functions and pathways, we intersected the list of significantly correlating DEGs with autism candidate genes from the SFARI database using hypergeometric distribution analysis ( Figure S2a). A Venn diagram revealed that 35 genes from among 309 DEGs correlated with DETEs were significantly enriched among autism candidate genes (p = 3.67 × 10 −6 ); the list of overlapping genes associated with Alu elements is shown in Figure S2b. This finding suggests that DE-Alu influences autism gene expression in the prefrontal cortex, which may lead to risk of ASD.

Consensus Sequence and Transcription Factor Binding Sites at DE-Alu Elements
Although Alu elements are abundant in the genome, only a small portion of Alu elements are active and escape from restriction by an epigenetic mechanism. To determine the characteristics of the DE-Alu sequence, we conducted multiple sequence alignment analysis and identified the consensus sequence from each DE-Alu compared with nondifferentially expressed elements (non-DE) by random selection. To draw the closest comparison between DE-Alu and non-DE-Alu, we randomly chose the non-DE-Alu sequence equally to the DE-Alu sequence (3881). This type of selection with similar numbers (~5000) has been performed in a previous study [28]. We performed an analysis of the AluS and AluY families, which are the major DETEs in the RNA-seq results. The analysis focused on the left arm monomer of Alu elements, which is the conserved region and contains promoters. A comparison between the consensus sequence and the conservation score indicated by the sequence logo and bar chart, respectively, showed a highly conserved region in the DE-AluS family compared with non-DE-AluS. The red lines in Figure 5a,b indicate higher conservation. AluY showed a slightly conserved region between DE-AluY and non-DE-AluY ( Figure S3). This difference may be the reason for the abnormal expression of DETEs, as controlled by upstream regulators that ultimately affect downstream genes.  We next predicted transcription factor binding sites at DE-AluS and non-DE-AluS consensus sequences using the JASPAR CORE database. We selected only the conserved region for scanning using 951 matrix profiles corresponding to 681 TFs (filtered by Homo sapiens) and found 191 matrix profiles (120 TFs) with high relative scores (relative score thresholds = 0.8) in the DE-AluS and 243 matrix profiles (99 TFs) in non-DE-AluS. When overlapping TFs between DE-AluS and non-DE-AluS, we detected 68 unique TFs binding to the DE-AluS consensus sequence (Table 3, Figure 5c). We further identified unique TFs associated with ASD by overlapping with ASD candidate genes from the SFARI database ( Figure 5d). In total, 12 TFBSs within the DE-AluS consensus sequence corresponding to nine ASD-associated TFs (Table S7), including NR2F1, VDR, NR4A2, NR1D1, RORA, SATB1, ARNT2, NR4A2, AR, TCF7L2, and NR4A2, were found. These findings provide a better understanding of the upstream regulators of DE-Alu elements associated with ASD. Table 3. List of unique ASD-associated TFs binding at the DE-AluS consensus sequence predicted by JASPAR CORE database. "+" represents the sense strand and "-" represents the antisense strand.

DNA Methylation of Alu Is Altered in the Prefrontal Cortex of ASD Individuals
We conducted a pilot study using human postmortem brain tissues of ASD and unaffected controls that were obtained from the Harvard Brain Tissue Resource Center, consisting of 13 samples of prefrontal cortex tissues from ASD individuals (n = 7) and unaffected controls (n = 6). There were no significant differences in age, sex, or postmortem interval (PMI) between the ASD group and the unaffected control group. Our postmortem ASD brain tissue was matched with the ASD sample used for RNA-seq reanalysis using brain IDs, and we applied the COBRA method to determine the global DNA methylation of the AluSx subfamily, which is the most abundant Alu element in the genome. The COBRA method was designed to measure DNA methylation levels ( m Cs) and patterns ( m C m C, m C u C, u C m C, u C u C) at AluSx subfamily promoters according to our previous study [20]. The percentage of AluSx methylation levels and patterns from the COBRA analysis are shown in Table 4. The percentages of partially methylated patterns (% u C m C) were significantly increased in ASD compared with sex-and age-matched controls (∆% u C m C = 0.285, p = 0.037) ( Table 4). Due to the heterogeneity of ASD, we classified ASD individuals based on the ADI-R score for verbal-ASD (V-ASD) and nonverbal ASD (NV-ASD), which were obtained from multiple ASD studies (Table S1). Interestingly, we found differences in AluSx methylation patterns in the V-ASD individuals compared with the sex-and age-matched controls, including total methylation (∆% m C = −1.159, p = 0.019), hypermethylated patterns (∆% m C m C = −2.650, p = 0.004), and partially methylated patterns (∆% u C m C = 0.462, p = 0.023; ∆% m C u C = 1.225, p = 0.031); however, we did not observe AluSx methylation changes in the NV-ASD individuals compared with the sex-and age-matched controls. Because changes in the DNA methylation level of Alu elements were observed, we also determined the expression levels and relative copy numbers of AluSx elements in postmortem brain tissues. Table 4. COBRA-derived percentages of AluSx methylation levels and patterns in the prefrontal cortex of ASD individuals and sex-and age-matched unaffected controls. There was no statistically significant difference between ASD individuals and matched controls with regard to AluSx expression level or relative copy number (Figure 6a,c) or among ASD subtypes based on ADI-R scores (Figure 6b,d); however, correlation analysis revealed an association between DNA methylation and AluSx activity in the prefrontal cortex of ASD individuals and unaffected controls (Figure 6e-g). The correlation between partially methylated patterns (% u C m C and % m C u C) that contained one unmethylated CpG showed a moderate positive correlation with AluSx copy number (R = 494 and R = 456, Figure 6e). The correlation between DNA methylation and expression of AluSx showed a weak correlation in the ASD and unaffected individuals (R = 0.177). Additionally, correlation analysis involving the ASD subtype revealed that significant AluSx methylation levels and patterns showed moderate to high correlation levels with AluSx expression and relative copy number in the V-ASD subtype (Pearson's correlation; |R| ranges = 0.417-0.833) (Figure 6f,g). These findings indicate that DNA methylation of AluSx, which is altered in the prefrontal cortex of ASD individuals, contributes to Alu expression and copy number in prefrontal cortex tissues. is altered in the prefrontal cortex of ASD individuals, contributes to Alu expression and copy number in prefrontal cortex tissues.

Correlation Analysis between ADI-R Score and DNA Methylation and Expression of Alu Elements in ASD Brain Tissues
To determine whether Alu elements are related to ASD symptoms, we correlated the DNA methylation, expression, and relative copy number of Alu elements with the ADI-R scores of ASD individuals. First, the correlation between ADI-R and Alu elements in the prefrontal cortex revealed the Alu methylation level and pattern (% m C and % m C m C) to be related to worse ADI-R scores, especially in the repetitive behavior domain (ADI-R C) and the communication domain for nonverbal individuals (ADI-R B NV). Partially methylated patterns (% m C u C and % u C m C) correlated positively with ADI-R (ADI-R C, ADI-R B NV, and ADI-R B V). Moreover, the expression and relative copy number of Alu elements correlated highly with the communication domain for both NV and V (ADI-R B NV and ADI-R B V) (Figure 7). These findings indicate that Alu elements may be related to ASD symptoms based on ADI-R scores and deserve further study using animal models and behavioral examinations.

Confirmation of Expression Levels of DEGs Associated with Alu Elements in Postmortem Brain Tissues
To confirm the expression of significantly correlated DEGs, we performed qPC analyses using our postmortem brain tissues of ASD (n = 7) and unaffected (n = 6) ind viduals. We selected DEGs from the most significantly correlating DEGs or IPA resul associated with the neuropathology of the cortex (i.e., ATP8A, RNF135, GRIN2A, GRIN2 HMOX1, DOCK1, KCNQ5, and PARP9) ( Table 5). The qPCR results showed significant

Confirmation of Expression Levels of DEGs Associated with Alu Elements in Postmortem Brain Tissues
To confirm the expression of significantly correlated DEGs, we performed qPCR analyses using our postmortem brain tissues of ASD (n = 7) and unaffected (n = 6) individuals. We selected DEGs from the most significantly correlating DEGs or IPA results associated with the neuropathology of the cortex (i.e., ATP8A, RNF135, GRIN2A, GRIN2B, HMOX1, DOCK1, KCNQ5, and PARP9) ( Table 5). The qPCR results showed significantly decreased PARP9 gene expression in ASD compared with unaffected controls (Figure 8a) (log 2 FC = −0.87, p = 0.038). Moreover, we classified ASD individuals based on ADI-R B-for nonverbal (NV) and verbal (V) individuals as NV-ASD and V-ASD individuals. Expression of the RNF135 and HMOX1 genes was found to be significantly increased in NV-ASD individuals compared to sex-and age-matched controls (Figure 8b) (RNF135: log 2 FC = 3.75, p < 0.01, and HMOX1: log 2 FC = 2.60, p = 0.045). There were no significant genes observed in V-ASD compared with sex-and age-matched controls (Figure 8c). These results show the reproducibility of DEGs associated with Alu elements in both RNA-seq data and our postmortem brain tissues. show the reproducibility of DEGs associated with Alu elements in both RNA-seq data and our postmortem brain tissues.

RNF135 and HMOX1 Genes Are Related to Neuron Number in the Prefrontal Cortex
According to biological function and gene regulatory network analyses, we found that DEGs correlating with Alu elements are associated with the survival of neural cells and the cell death of cerebral cortex cells. We hypothesized that Alu elements may be associated with genes involved in the neuropathology found in the postmortem brain tissues of ASD individuals. Total neuron density was quantified in prefrontal cortex tissues (ASD; n = 3, unaffected; n = 5) stained with the Nissl staining method (Figure 9a). According to this analysis, neuron density in the prefrontal cortex between ASD and unaffected individuals was significantly increased in the former (p = 0.042) (Figure 9b). We next performed correlation analysis between neuron density and the expression level of DEGs associated with Alu elements using Pearson's correlation (Figure 9c), with a moderate to high correlation (|R| = 0.599-0.931) with neuron density in the prefrontal cortex, except for DOCK1 (R = 0.213). In particular, the RNF135 and HMOX1 genes, which were significantly differentially expressed, as confirmed by qPCR, showed a higher correlation with neuron density than the other genes. It is important to note that there was a limitation regarding the number of paraffin-embedded tissues available, especially for matching with frozen tissues for qPCR analyses. Thus, we did not analyze ASD individuals by classifying them into NV-ASD and V-ASD groups. Further study using a larger sample size is needed.

Locus-Specific DNA Methylation of AluYk3 Located in the Upstream Region of the RNF-135 Gene
To investigate the susceptibility loci of Alu elements in DEGs in brain tissues, we applied the COBRA method to assess locus specificity at three CpG sites of AluYk3 elements located upstream of the RNF-135 gene ( Figure 10a). As in the previous analyses, we compared methylation levels and patterns between all ASD, NV-ASD, and V-ASD individuals versus sex-and age-matched controls. The results showed no significant differences in AluYk3 methylation between all ASD individuals versus unaffected individuals (Figure 10b-d, Table S8). Interestingly, we found one differentially methylated CpG position at AluYk3 elements, whereby the CpG second position was hypomethylated in the

Locus-Specific DNA Methylation of AluYk3 Located in the Upstream Region of the RNF-135 Gene
To investigate the susceptibility loci of Alu elements in DEGs in brain tissues, we applied the COBRA method to assess locus specificity at three CpG sites of AluYk3 elements located upstream of the RNF-135 gene ( Figure 10a). As in the previous analyses, we compared methylation levels and patterns between all ASD, NV-ASD, and V-ASD individuals versus sex-and age-matched controls. The results showed no significant differences in AluYk3 methylation between all ASD individuals versus unaffected individuals (Figure 10b-d, Table S8). Interestingly, we found one differentially methylated CpG position at AluYk3 elements, whereby the CpG second position was hypomethylated in the NV-ASD group (∆%methylated CpG no.2 = −6.79, p = 0.033) but not in all ASD or V-ASD groups compared with sex-and age-matched controls (Figure 10e-g). We performed correlation analysis between AluYk3 methylation and RNF135 expression levels and found that the %methylated CpG second position at AluYk3 showed a negative correlation with RNF135 (R-value = −0.575) (Figure 10h,i). These findings suggest that AluYk3 located 30 kb upstream of the RNF135 gene might contribute to the regulation of RNF135 gene expression levels in the prefrontal cortex of ASD individuals via DNA methylation. Nevertheless, we did not detect any significant difference in other selected Alu elements located in other DEGs, such as the HMOX1 gene (Table S8). Notably, it is possible that other CpGs or other epigenetic mechanisms of these Alu elements may have contributed to host DEG expression. Due to the limitation of the COBRA method, only 1-3 CpG sites at the promoter region of Alu elements were detected for these analyses. We also performed correlation analysis between AluYk3 methylation and the ADI-R scores of ASD individuals.

Discussion
At present, it is not known exactly what causes abnormal gene expression levels in postmortem brain tissues of ASD individuals. Most analyses of DNA methylation and transcriptome studies in postmortem brain tissues of ASD individuals have focused on protein-coding regions rather than noncoding regions, with the latter contributing most of the CpG sites in the human genome, including TEs. TEs can be classified into two types based on their structure and the requirement of reverse transcription for transposition: DNA transposons and retro (RNA) transposons [17]. There are two main types of re-

Discussion
At present, it is not known exactly what causes abnormal gene expression levels in postmortem brain tissues of ASD individuals. Most analyses of DNA methylation and transcriptome studies in postmortem brain tissues of ASD individuals have focused on protein-coding regions rather than noncoding regions, with the latter contributing most of the CpG sites in the human genome, including TEs. TEs can be classified into two types based on their structure and the requirement of reverse transcription for transposition: DNA transposons and retro (RNA) transposons [17]. There are two main types of retrotransposons: long terminal repeats (LTRs: ERV) and non-long terminal repeats (non-LTRs), and most retrotransposons in the human genome are the latter [17,29]. The most common non-LTR retrotransposons are long interspersed nuclear element 1 (LINE-1 or L1) and short interspersed nuclear element (SINE: Alu and MIR) transposons, which remain active and comprise approximately 30% of the human genome at a copy number of over 1 million elements [17,29]. Dysregulation of TEs can affect developmental processes via their gene product through new insertions that cause genetic changes. Moreover, retrotransposons introduce a promoter or premature terminator (poly-A site) into host genes, causing transcriptional disruption. Retrotransposons can also act as enhancers of host/nearby genes because they contain several transcription factor binding sites [17].
Several pieces of evidence suggest that Alu elements and other TEs influence mammalian genome evolution by introducing novel gene formation, transcriptome diversity, and transcriptional regulatory elements [18,30]. Alu elements are beneficial for the development of the central nervous system, including brain formation and function [30,31]. Additionally, previous studies found that Alu controls several genes that are important for neuronal function, such as the serotonin transporter gene [32][33][34]. However, a loss of control of Alu element activity by epigenetics, due to environmental and other nongenetic factors, contributes to many human diseases, including neurological diseases [30]. This evidence supports the idea that epigenetic dysregulation in ASD brain tissues may modulate the activities of Alu elements and alter the transcriptome.
In this study, we characterized the TE expression profile of the prefrontal cortex of ASD individuals. We found differentially expressed TEs from the retrotransposon family, including the Alu, LINE, and MIR families (Figure 2b). Due to the high abundance of TEs in the human genome, we also considered the total number of elements in each family together with the number of differentially expressed elements (Figure 2c). Although we observed many dysregulated elements belonging to the LINE-2 (L2), LINE-3 (L3), and MIR families, it appears that the large number may be due to the total number of elements. Another interesting TE family is the dysregulated LINE-1 family (L1), which is an active family of LINEs. When we considered subfamilies of L1, dysregulated L1 was classified as belonging to the L1M subfamily, which is the oldest and least active family of L1 subfamilies [35]. The most active and youngest subfamily of L1 (L1H) was not detected in our analysis. This is consistent with a study by Shpyleva, who performed qPCR analysis of the L1H subfamily in ASD using the cerebellum and prefrontal cortex, with L1H being significantly overexpressed only in the former [25]. Indeed, we focused on the dysregulated elements of the Alu family because this TE family has a large number of dysregulated elements compared with the total number of Alu and other TE families. Moreover, the highest number of dysregulated Alu belonging to the AluSx and AluY subfamilies constitutes a middle and young subfamily of Alu elements [35].
We further investigated how this Alu dysregulation affects ASD neuropathology. At present, our understanding of the precise role of Alu elements in ASD development remains limited. Their possible function may be as a result of LINE-1 and Alu elements acting as enhancers or promoters for host genes [17]. Alu retrotransposons are mostly distributed in gene-rich regions close to the TSS of genes [36,37] and contain many transcription factor (TF) binding sites [38]. According to the study of Su and colleagues in 2014, Alu retrotransposons are significantly conserved at the upstream region of genes and enriched for enhancer-like epigenetic regulation [26]. Furthermore, the epigenetic state of Alu retrotransposons located upstream of genes might regulate the expression and chromatin structure of genes associated with the cell cycle [39]. In our study, we expected that dysregulation of the Alu family would reflect epigenetic alterations at Alu elements, which may impact host/nearby gene expression. We thus performed correlation analysis between the expression of Alu and their host genes. From 1493 DEGs in the prefrontal cortex of ASD individuals, as many as 161 loci corresponding to 133 DEGs correlated significantly with Alu elements located within/nearby host DEGs. The biological function of these genes regulated by Alu is associated with ASD and comorbid disorders. In addition to canonical pathway analyses, these correlating DEGs were significantly enriched in the synaptogenesis signaling pathway, neuroinflammation signaling pathway, cell cycle, and glutamate receptor signaling. Accumulating evidence suggests that a deficiency of genes in pathways involving synaptic protein synthesis and degradation, postsynaptic scaffold architecture, and neurotransmitter receptors contribute to synapse deficits in ASD individuals [40]. Moreover, several studies have demonstrated that children with ASD experience immune system changes at both systemic and cellular levels. In fact, individuals with ASD exhibit signs of active inflammation and altered gene expression in immune signaling and functional pathways [41]. Accordingly, we investigated how these Alu elements are dysregulated and affect gene activity in the brain tissues of ASD individuals.
In this study, we found that DE-AluS elements have a conserved region within the left arm monomer (promoter) when compared with non-DE-AluS. This suggests that the conserved region may contribute to a dysregulation of DE-AluS in ASD. We also evaluated transcription factor binding sites at the conserved region of DE-AluS in the prefrontal cortex of ASD individuals. Interestingly, we found that nine TFs binding to this region comprise autism candidate genes in the SFARI database. An interesting TF is retinoic acid-related orphan receptor alpha (RORA), a nuclear hormone receptor which is downregulated in both brain tissues and lymphoblastoid cell lines from individuals with ASD versus controls [42]. Interestingly, RORA not only transcriptionally regulates the aromatase enzyme, which converts testosterone to estrogen, but is itself under both negative and positive feedback regulation by male and female sex hormones, respectively [42]. These observations led Sarachana et al., (2011) to suggest that a decrease in RORA and its transcriptional target aromatase may be related to the reported elevated testosterone in amniotic fluid associated with ASD traits, which has been proposed as a factor influencing the male sex bias in ASD [43]. Sex differences in TEs in the human genome may be related to the biology of sex differences, such as higher TE density on the Y chromosome, differences in DNA methylation regulators (DNA methyltransferase enzymes and methylcytosine dioxygenases; TET), and histone modification marks [44]. Another interesting TF is vitamin D receptor (VDR), a nuclear receptor/steroid hormone receptor superfamily. VDR regulates several targeted genes involved in cellular proliferation and differentiation [45]. Genetic polymorphisms in VDR are associated with risk of ASD [46]. There is evidence that approximately 50% of children with ASD have vitamin D deficiency and that 30% have insufficient vitamin D, with serum vitamin D correlating negatively with ASD severity [47]. In the brain, vitamin D has several benefits, including calcium signaling and the synthesis of neurotrophic factors and neurotransmitters [48]. Our findings indicate that upstream regulators of TEs may be key factors for transcriptome diversity in ASD brain tissues, and further study is warranted. Lowe et al., (2007) reported that retrotransposons are most often located in genes involved in development and transcriptional regulation, which may reflect that retrotransposons are involved in a common mechanism for development [49]. In addition to epigenetic alterations of retrotransposons at the integration site, including DNA methylation and chromatin remodeling, retrotransposon expression can be restricted and lead to the expression of host/neighboring genes [50]. Retrotransposon activity is suppressed through epigenetic mechanisms that silence their expression, and Alu hypomethylation has been reported in response to environmental exposures [51][52][53]. Current evidence links aberrant LINE-1 and Alu DNA methylation to human brain disorders, especially in patients with schizophrenia and ASD [54]. In the ASD brain, whole-genome sequencing analysis has revealed a higher number of novel insertions of LINE-1 and Alu elements than in normal brain tissues [24].
We have previously reported alterations in the global DNA methylation of LINE-1 and Alu retrotransposons in lymphoblastoid cell lines (LCLs) of ASD subtypes, and such alterations were specific to the ASD case subgroup based on clinical phenotype [20,21]. An alteration in LINE-1 methylation in LCLs was specific for the ASD subtype with severe language impairment, whereas an alteration in Alu methylation was specific for ASD with mild symptoms. Moreover, several genes with LINE-1 and Alu insertions were differentially expressed in the peripheral blood and LCLs of ASD. These genes are known to be associated with neurodevelopmental disorders. In the present study, we used COBRA to determine global AluSx methylation levels and patterns in the prefrontal cortex of ASD individuals. Although there was no significant difference in the overall methylation level (% m C) of the Alu promoter, we found that the percentage of partially methylated patterns containing one unmethylated CpG in the promoter was significantly increased (% u C m C) in all ASD individuals versus unaffected individuals. Interestingly, these methylation patterns exhibited moderate correlations with AluSx expression and copy number in the prefrontal cortex of ASD individuals.
When we classified ASD individuals based on the ADI-R score for nonverbal and verbal ASD, we found differences in the Alu methylation in the brain of individuals with each ASD subphenotype. For global Alu methylation, Alu methylation was changed specifically in the prefrontal cortex of V-ASD individuals. Correlation analysis between Alu methylation, expression, and relative copy number showed a high correlation in V-ASD; however, the Alu expression level and relative copy number did not differ between ASD individuals and unaffected controls, even when ASD individuals were subtyped. This is unnecessary for disrupting host gene expression because Alu elements or other TEs can provide regulatory elements for cis-and trans-regulation to interfere with Pol-II transcription, whereas Alu elements are transcribed by Pol-III [55]. Moreover, correlation analysis with the ADI-R scores of ASD individuals showed Alu elements to be associated with ASD symptoms. Based on this finding, the role of Alu elements in autism-related behaviors deserves further study.
In this research, we sought to identify autism-susceptibility loci of TEs in the brains of ASD individuals. We selected host DEGs and investigated the methylation levels of Alu elements located within or near host DEGs. First, we confirmed the expression level of host DEGs using qPCR and found that HMOX1, PARP9, and RNF135 were differentially expressed in the prefrontal cortex of ASD individuals. Unfortunately, when we implemented the COBRA method to assess the locus-specific DNA methylation of Alu elements, we found only one differentially hypomethylated locus at AluYk3 elements. AluYk3 elements were found to be located~30 kb upstream of the RNF135 gene. RNF135 was overexpressed in both RNA-seq data and qPCR analysis; interestingly, AluYk3 methylation correlated significantly with RNF135 expression. This result suggests that AluYk3 may provide cisregulatory elements for the RNF135 gene located nearby. RNF135, or Ring Finger Protein 135 or RING-Type E3 Ubiquitin Transferase RNF135 gene, is an autism candidate gene in the SFARI database involved in protein-protein and protein-DNA interactions. RNF135 has been associated with ASD [56,57] and neurofibromatosis [58] and reported to be involved in intronic Alu-mediated recombination in NF1 patients [58]. Regardless, the role of RNF135 in neurodevelopment remains unclear, though there is evidence showing that RNF135 can promote the proliferation of human glioblastoma cells in vivo and in vitro [59]. We also observed that this region of pathogenic copy number variation contains several genes (NF1, RNF135, ADAP2, and SUZ12) and shows many Alu-mediated recombination events [58,60].
Finally, we sought to determine whether Alu elements are associated with neuropathology in the ASD brain. In our analysis, neuron density in the prefrontal cortex of ASD individuals was increased compared with unaffected individuals. Our correlation analysis with qPCR results showed that neuron number in the prefrontal cortex correlated highly with Alu-correlating genes. These findings suggest that Alu elements may contribute to neuropathology in the prefrontal cortex of ASD individuals by disrupting the gene expression involved in ASD and neuropathology, such as cell death and survival. The prefrontal cortex area is responsible for memory, verbal attention, and language processing, which are strongly linked to ASD [61][62][63]. Several studies have shown changes in cytoarchitecture in the brains of ASD individuals, such as neuron and glial cell number or density. Studies by Falcone et al. and Courchesne et al. also demonstrated increased neuron density in the prefrontal cortex of ASD individuals compared with unaffected controls [64,65]; however, there is inconsistency among studies with respect to the number of neurons in the cerebral cortex of ASD individuals [66,67]. It is very important to note that our finding is from a preliminary study highlighting the impact of TEs, especially Alu elements, on the neuropathology in the prefrontal cortex of ASD individuals.
Because many comparable sequences have low mappability, TE profiling is exceedingly difficult. The techniques used to prepare RNA libraries, requirements for certain parameters to address their mappability, tools used to detect their transcripts, and genomic locations can all have an impact on these limitations. Thus, some active young TEs may not have been detected in our analysis due to low mappability and low abundance. Moreover, our analysis did not detect the potential expression of TEs not in the hg38 reference genome. Other limitations include use of postmortem brain tissues, the small sample size because of the limited availability of brain tissues, no pilot study for calculating a suitable sample size, the characterization of methylated DNA-binding proteins, and functional consequences. There was also a limitation regarding the COBRA method used for DNA methylation analysis due to restriction sites of the enzyme, and the CpG sites determined in our study may not reflect the entire methylation status of the promoter of Alu elements.
Based on our findings, we propose a molecular mechanism related to neuropathology in ASD brain tissues ( Figure 12). Environmental risk factors for ASD induce a dysregulation of epigenetic mechanisms in the brain, leading to changes in the activities of Alu elements, as well as other TEs. The dysregulation of TEs contributes to the transcriptome through two possible mechanisms: (i) altering the transcriptome profile of the brain via cis-and trans-regulation of targeted genes and (ii) directly changing the transcriptome profile via functional products, including transcripts (that directly interfere with RNA) and DNA (which alters DNA sequence). Both mechanisms may be related to autism-associated TFs that lead to a disruption of the interactome, which is related to the neuropathology of ASD brain tissues. Moreover, our analysis did not detect the potential expression of TEs not in the hg38 reference genome. Other limitations include use of postmortem brain tissues, the small sample size because of the limited availability of brain tissues, no pilot study for calculating a suitable sample size, the characterization of methylated DNA-binding proteins, and functional consequences. There was also a limitation regarding the COBRA method used for DNA methylation analysis due to restriction sites of the enzyme, and the CpG sites determined in our study may not reflect the entire methylation status of the promoter of Alu elements. Based on our findings, we propose a molecular mechanism related to neuropathology in ASD brain tissues ( Figure 12). Environmental risk factors for ASD induce a dysregulation of epigenetic mechanisms in the brain, leading to changes in the activities of Alu elements, as well as other TEs. The dysregulation of TEs contributes to the transcriptome through two possible mechanisms: (i) altering the transcriptome profile of the brain via cis-and trans-regulation of targeted genes and (ii) directly changing the transcriptome profile via functional products, including transcripts (that directly interfere with RNA) and DNA (which alters DNA sequence). Both mechanisms may be related to autism-associated TFs that lead to a disruption of the interactome, which is related to the neuropathology of ASD brain tissues. Figure 12. A schematic diagram illustrating a proposed mechanism by which Alu elements contribute to neuropathology in ASD brain tissues, as created with BioRender.com (accessed on 4 January 2023). Figure 13 describes the overall strategy and experimental workflow of this study. To investigate relationships between TEs, the transcriptome, neuropathology, and clinical phenotype (ADI-R scores) in the brain tissues of individuals with ASD, we obtained RNAseq data from GEO DataSets for TE profiling and identified host/target genes that may be influenced by TEs. Then, we examined the sequence characteristics and TF binding sites Figure 12. A schematic diagram illustrating a proposed mechanism by which Alu elements contribute to neuropathology in ASD brain tissues, as created with BioRender.com (accessed on 4 January 2023). Figure 13 describes the overall strategy and experimental workflow of this study. To investigate relationships between TEs, the transcriptome, neuropathology, and clinical phenotype (ADI-R scores) in the brain tissues of individuals with ASD, we obtained RNAseq data from GEO DataSets for TE profiling and identified host/target genes that may be influenced by TEs. Then, we examined the sequence characteristics and TF binding sites of differentially expressed TEs and predicted the biological functions of host and target genes. Finally, postmortem brain tissues obtained from the Harvard Brain Tissue Resource were used to investigate DNA methylation of Alu elements and neuropathology. Correlation analysis was conducted to investigate the relationship between the results. of differentially expressed TEs and predicted the biological functions of host and target genes. Finally, postmortem brain tissues obtained from the Harvard Brain Tissue Resource were used to investigate DNA methylation of Alu elements and neuropathology. Correlation analysis was conducted to investigate the relationship between the results.

Data Collection
To investigate TE expression profiling in the prefrontal cortex of ASD individuals, we obtained publicly available data from the NCBI Gene Expression Omnibus database (GEO data: https://www.ncbi.nlm.nih.gov/gds) accessed on 29 January 2020 using the following criteria: (i) ASD studies that used postmortem brain tissues from the prefrontal cortex of ASD and unaffected individuals; (ii) studies that used RNA sequencing and provided raw data; (iii) sample sizes greater than or equal to 40; and (iv) RNA-seq data of male individuals. A total of four studies were obtained [16,[68][69][70]. RNA-seq raw data were obtained from the NCBI Sequence Read Archive repository (SRA: https://www.ncbi.nlm.nih.gov/sra) accessed on 29 August 2020.

RNA-Seq Data Reanalysis for Transposable Element Profiling in the Prefrontal Cortex of ASD Individuals
All RNA-seq data (FASTQ files) for ASD and unaffected individuals in the SRA repository were downloaded to the Galaxy program (https://usegalaxy.org) accessed on 12 August 2020, and analyses were performed using Galaxy program tools [71]. First, the FASTQ files were cleaned, and quality control was performed using the fastp package (fast

Data Collection
To investigate TE expression profiling in the prefrontal cortex of ASD individuals, we obtained publicly available data from the NCBI Gene Expression Omnibus database (GEO data: https://www.ncbi.nlm.nih.gov/gds) accessed on 29 January 2020 using the following criteria: (i) ASD studies that used postmortem brain tissues from the prefrontal cortex of ASD and unaffected individuals; (ii) studies that used RNA sequencing and provided raw data; (iii) sample sizes greater than or equal to 40; and (iv) RNA-seq data of male individuals. A total of four studies were obtained [16,[68][69][70]. RNA-seq raw data were obtained from the NCBI Sequence Read Archive repository (SRA: https://www.ncbi.nlm.nih.gov/sra) accessed on 29 August 2020.

RNA-Seq Data Reanalysis for Transposable Element Profiling in the Prefrontal Cortex of ASD Individuals
All RNA-seq data (FASTQ files) for ASD and unaffected individuals in the SRA repository were downloaded to the Galaxy program (https://usegalaxy.org) accessed on 12 August 2020, and analyses were performed using Galaxy program tools [71]. First, the FASTQ files were cleaned, and quality control was performed using the fastp package (fast all-in-one preprocessing for FASTQ files). For TE expression profiling, clean reads were mapped to the human reference genome (hg38) using RNA STAR aligner [72] with a parameter for multimapped reads according to Teissandier's study [73]. To quantify TE transcripts, the mapped reads of individuals were counted using the Subread package FeatureCounts [74] with multimode (-M, -f). RepeatMasker annotation for hg38 was downloaded from the UCSC Table Browser [75]. To quantify protein-coding genes, the unique mode (-U) and the primary alignment mode of FeatureCounts were used for quantification with FeatureCounts built-in annotations for hg38. We excluded a few individuals from the differential expression analysis because their unique and multimapped ratios were markedly different.
Differential expression analyses of TEs and protein-coding genes were performed separately using the R package DESeq2 for ASD versus unaffected groups [76]. To normalize sequencing depth and RNA composition, read counts were normalized by the geometric mean method. Moreover, we normalized the data for other variable factors (i.e., library preparation and other unwanted technical effects) using the remove unwanted variation (RUV) tool [77]. Genes and TEs with a P and FDR of less than 0.01 were considered significant.

Identification of Genomic Loci of Differentially Expressed TEs Located Within/near Protein-Coding Genes
To determine whether TEs are associated with the expression level of protein-coding genes, we identified the genomic locations of each differentially expressed TE (DETE) located within or near a protein-coding gene region using the bedtools package [78]. First, the genomic loci of each TE were retrieved from RepeatMasker annotation for hg38 (http://repeatmasker.org) accessed on 1 November 2020, and NCBI RefSeq genes were downloaded from the UCSC Table Browser [79]. Then, the genomic loci of DETEs were intersected with the NCBI RefSeq genes (exonic regions) using the bedtools Intersect intervals function [78]. These processes were repeated with each intronic region and regions upstream of the TSS (50 kb, 10 kb, and 1 kb).

Cis-and Trans-Regulation of DETEs and the Expression Level of Nearby and Distant Genes
To determine whether dysregulation of TEs can affect host gene expression in cis-and trans-regulation, we performed Pearson's correlation analyses between the normalized expression value (rlog transformation) of individual DETEs and their host genes. Correlations between expression of Alu elements and host genes with an adjusted P less than 0.05 (FDR = 0.05) were considered statistically significant. In this study, we also propose trans-regulation by Alu elements. We performed correlation analysis between all probabilities between 864 DETEs and 1493 DEGs, and an adjusted p < 0.05 (FDR = 0.05) was considered statistically significant.

Identification of Biological Functions, Pathways, and Autism Candidate Genes of Differentially Expressed Genes Associated with TEs
To predict the biological function and pathway of differentially expressed genes correlating significantly with TEs, the list of genes with log2-fold change and adjusted P were selected. Biological functions, pathways, and gene regulatory networks were predicted by IPA software (QIAGEN Inc., Hilden, Germany, https://www.qiagenbioinformatics. com/products/ingenuitypathway-analysis) accessed on 18 June 2021, and p < 0.05 was considered statistically significant. To identify autism candidate genes regulated by TEs in the prefrontal cortex of ASD individuals, genes were obtained from the SFARI database. The list of differentially expressed genes associated with TEs overlapped with 1003 autism candidate genes from the SFARI database [80] when using a Venn diagram drawing tool (Venny 2.1 software) [81], and gene set enrichment analysis was performed using the hypergeometric distribution calculator (http://keisan.casio.com/exec/system/1180573201) accessed on 22 February 2021.

Multiple Sequence Alignment Analysis and Transcription Factor Binding Site Prediction of Alu Elements
To determine the sequence characteristics of Alu elements differentially expressed in the postmortem brain tissues of ASD individuals, we extracted the sequence (FASTA) of each DE-Alu element from the human genome hg38 using the bedtools GetFastaBed function. The DE-Alu element sequences were analyzed using the multiple sequence alignment package (msa, R package) [82] by comparing the consensus sequence of DE-Alu elements with non-DE-Alu elements randomly chosen equally to the number of DETEs. The sequence logo and conserved scores were created using the ggseqlogo package [83].
For TF binding site prediction at the conserved region of DE-Alu elements, sequences were scanned using the matrix profile of a known TF binding site (Homo sapiens) in the JASPAR CORE database [84]. We selected all matrix profiles of Homo sapiens for scanning the conserved region identified by msa analysis and performed scanning with a relative score threshold >0.8.

DNA Extraction from Postmortem Brain Tissues
In this study, prefrontal cortices from ASD individuals (n = 7) and unaffected individuals (n = 6) were kindly provided by Dr. Valerie W. Hu (George Washington University), as originally obtained from the Harvard Brain Tissue Resource Center. First, approximately 100 mg of tissue was washed with PBS on ice three times. Then, lysis solution and 0.5 mm glass beads were added to the microcentrifuge tube. Homogenization was performed by using Disruptor Genie (Scientific Industries, Bohemia, NY, USA) at high speed at 4 • C for 2 min; this step was repeated if homogenization was incomplete. Acid phenol/chloroform was added to the homogenate, which was centrifuged at 10,000× g for 5 min. Genomic DNA was isolated using GENEzol reagent and precipitated by adding 100% ethanol to the interphase and the lower phenol-chloroform phase, followed by centrifugation at 2000× g for 5 min. The DNA pellet was resuspended in 0.1 M sodium citrate in 10% ethanol and incubated for 30 min before centrifuging and discarding the supernatant; this step was repeated twice. The DNA pellet was washed with 75% ethanol and air-dried for 5-10 min; the DNA pellet was resuspended and stored in 8 mM NaOH and 1 mM EDTA. The quality and quantity of purified DNA was measured using a NanoDrop One (Thermo Fisher Scientific, Inc., Rockingham County, NH, USA).

Combined Bisulfite Restriction Analysis (COBRA) for Alu Methylation Level and Patterns
The first step of the COBRA technique is bisulfite conversion, in which 1 µg of gDNA from human brain tissues was treated with sodium bisulfite using an EZ DNA Methylation-Gold™ Kit (Zymo, Irving, CA, USA) according to the manufacturer's protocol. Then, the bisulfite-treated DNA was amplified by 35 cycles of PCR (Hot Start Taq DNA polymerase, QIAGEN, USA) using specific primers for the global AluSx family based on the AluSx consensus sequence. The PCR amplicons were digested with a restriction enzyme and incubated at 65 • C overnight, and the digested products were electrophoresed through an 8% nondenaturing polyacrylamide gel. The band intensity was measured by using GelAnalyzer 19.1 (http://www.gelanalyzer.com; accessed on 24 August 2020). The percentages of Alu methylation levels and patterns were calculated using the following formulas according to previous studies: percentage methylated loci (% m C) = 100 × (E + B)/(2A + E + B + C + D), percentage of the hypermethylated pattern (% m C m C) = 100 × F/(A + C + D + F), percentage of the partially methylated pattern (% u C m C) = 100 × C/(A + C + D + F), percentage of the partially methylated pattern (% m C u C) = 100 × D/(A + C + D + F), and percentage of the partially hypomethylated pattern (% u C u C) = 100 × A/(A + C + D + F) [20,21]. All samples were analyzed in duplicate, and the SK-N-FI cell line was used for interassay normalization.
For locus-specific DNA methylation of DE-Alu elements located within or near host genes, candidate DE-Alu elements were selected using the following criteria: (i) Alu elements must correlate with their host/distant gene expression. (ii) Alu elements must have CpG sites at the promoter region that can be detected by the COBRA method. (iii) Host/distant genes of Alu elements must be the candidate differentially expressed genes found in IPA and interactome analyses. (iv) Host/distant genes of Alu elements must be associated with the neuropathology of ASD postmortem brain tissues. The genomic locations of selected Alu elements were obtained from RepeatMasker annotation and used to retrieve the sequence of each Alu element from the UCSC table browser. To obtain the sequence for designing specific primers, we also selected the upstream region of Alu elements (100 bp upstream region plus + Alu sequence). The COBRA primer was designed using Bisulfite Primer Seeker (https://www.zymoresearch.com/pages/bisulfite-primer-seeker) accessed on 14 February 2022, with specific forward primers spanning the upstream region (host DEG sequence) and reverse primers within the Alu sequence. This strategy increases the specificity of the primer because Alu elements have a high copy number. We also verified the primers by performing gradient PCR and gel electrophoresis to select the optimal temperature. Finally, COBRA was performed and the %methylation level assessed as described in global methylation analysis.

Quantitative RT-PCR Analysis
Quantitative RT-PCR was performed to determine expression levels and the relative copy number of AluSx and to confirm the expression levels of selected differentially expressed genes from RNA-seq. Briefly, total RNA was isolated from postmortem brain tissues (ASD; n = 7, unaffected controls; n = 6) using a mirVana miRNA isolation kit (Ambion, Foster City, CA, USA) according to the manufacturer's protocol. A total of 1 µg RNA was treated with DNaseI enzyme (Thermo Fisher Scientific, Waltham, MA, USA) before conversion to cDNA using a RevertAid RT reverse transcription kit (Thermo Fisher Scientific, USA). The qPCR assay was performed using RealMOD™ Green W 2 2× qPCR mix (iNtRON Biotechnology, Gyeonggi-do, Republic of Korea). Amplification was performed using a CFX Connect Real-Time PCR System (Bio-Rad, Hercules, CA, USA), consisting of an initial step at 95 • C for 15 min followed by 40 cycles of 45 s at 95 • C for denaturation and 30 s at 60 • C for annealing/extension. Product formation was confirmed by melting curve analysis (55 to 94 • C). The AluSx transcript-specific primers used were as follows: forward 5 -GTGGCTCACGCCTGTAATC-3 and reverse 5 -GTAGAGACGGGGTTTCACCA-3 . The number of AluSx transcripts was normalized to 18S RNA and calculated using the 2 −∆∆Ct method.
For AluSx copy number quantification, we performed direct qPCR to amplify AluSx copies from genomic DNA as described in previous studies [85,86]. Briefly, 250 pg of genomic DNA was used directly for the qPCR assay using AccuPower ® 2× GreenStar™ qPCR MasterMix (Bioneer, Daejeon, Republic of Korea) according to the manufacturer's instructions. For each sample, reactions were prepared in triplicate, in which the AluSx copy number was normalized by SATA elements and calculated using the 2 −∆∆Ct method.

Nissl Staining and Quantification of Neuron Number
To quantify neuronal cell number in the prefrontal cortex of ASD and unaffected individuals, paraffin-embedded tissues, also from the Harvard Brain Tissue Resource Center, were deparaffinized in xylene and stained with cresyl violet as follows: The tissue sections were deparaffinized in xylene for 10 min (2 steps) and then rehydrated in 100% ethanol solutions for 2 min (2 change steps) followed by decreasing concentrations of ethanol solution (1 step in 90% ethanol, 1 step in 80% ethanol, and 1 step in 70% ethanol; 1 min each step) and then placed in distilled water for 2 min. The deparaffinized tissue sections were stained with 0.5% cresyl violet (Sigma-Aldrich, St. Louis, MO, USA) for 2 min. Finally, the sections were rinsed in running tap water for 2 min followed by differentiation and dehydration in increasing concentrations of ethanol solution (1 step in 70% ethanol, 1 step in 80% ethanol, 1 step in 90% ethanol, and 2 steps in 100% ethanol; 5 min each step). The sections were captured by light microscopy (DM1000; Leica Microsystems, Wetzlar, Germany). We quantified neuronal cell numbers in three selected regions of each prefrontal cortex that covered the complete thickness of the cortical gray matter. We identified neurons based on their unique morphological characteristics by Nissl staining [87].

Statistical Analysis
Statistical analyses were performed using IBM SPSS version 22. For comparison of means, a two-tailed Student's t test was used. A p value less than 0.05 was considered statistically significant. Pearson's correlation analysis was employed for cis-and transregulation, DNA methylation, gene expression, neuron density, and ADI-R scores. For multiple test comparisons, the Benjamini-Hochberg correction method was applied at an FDR < 0.05. The statistical threshold for DETEs and DEGs was FDR = 0.01.

Conclusions
TE expression profiling revealed the Alu family to be dysregulated in the prefrontal cortex of ASD individuals. Moreover, a dysregulation of individual Alu elements correlated with host/nearby gene expression levels in cis-and trans-regulation. The DEGs associated with Alu elements are related to ASD and the neuropathology of ASD, and some of these genes are autism candidates in the SFARI database. In addition, COBRA analysis showed a change in the percentage of methylated patterns of two CpG sites in the Alu promoter in a subpopulation of individuals with ASD. Our findings provide a better understanding and highlight the role of Alu elements in gene dysregulation in the prefrontal cortex of ASD individuals. Further investigation of the relationship between Alu elements, gene regulation, and neuropathology in larger ASD sample sizes is warranted.

Informed Consent Statement: Not applicable.
Data Availability Statement: The datasets used as the main analysis in this study were obtained from GEO DataSets (GSE51264 and GSE59288).