Integrative analysis of long noncoding RNAs dysregulation and synapse-associated ceRNA regulatory axes in autism

Abstract

Autism spectrum disorder (ASD) is a complex disorder of neurodevelopment, the function of long noncoding RNA (lncRNA) in ASD remains essentially unknown. In the present study, gene networks were used to explore the ASD disease mechanisms integrating multiple data types (for example, RNA expression, whole-exome sequencing signals, weighted gene co-expression network analysis, and protein-protein interaction) and datasets (five human postmortem datasets). A total of 388 lncRNAs and five co-expression modules were found to be altered in ASD. The downregulated co-expression M4 module was significantly correlated with ASD, enriched with autism susceptibility genes and synaptic signaling. Integrating lncRNAs from the M4 module and microRNA (miRNA) dysregulation data from the literature identified competing endogenous RNA (ceRNA) network. We identified the downregulated mRNAs that interact with miRNAs by the miRTarBase, miRDB, and TargetScan databases. Our analysis reveals that MIR600HG was downregulated in multiple brain tissue datasets and was closely associated with 9 autism-susceptible miRNAs in the ceRNA network. MIR600HG and target mRNAs (EPHA4, MOAP1, MAP3K9, STXBP1, PRKCE, and SCAMP5) were downregulated in the peripheral blood by quantitative reverse transcription polymerase chain reaction analysis (false discovery rate <0.05). Subsequently, we assessed the role of lncRNA dysregulation in altered mRNA levels. Experimental verification showed that some synapse-associated mRNAs were downregulated after the MIR600HG knockdown. BrainSpan project showed that the expression patterns of MIR600HG (primate-specific lncRNA) and synapse-associated mRNA were similar in different human brain regions and at different stages of development. A combination of support vector machine and random forest machine learning algorithms retrieved the marker gene for ASD in the ceRNA network, and the area under the curve of the diagnostic nomogram was 0.851. In conclusion, dysregulation of MIR600HG, a novel specific lncRNA associated with ASD, is responsible for the ASD-associated miRNA-mRNA axes, thereby potentially regulating synaptogenesis.

Fig. 1. Flowchart of the overall approach.

ASD autism spectrum disorder, DE differential expression, GO gene ontology, GSEA gene set enrichment analysis, GS gene significance, GSVA gene set variation analysis, HC healthy control, KEGG Kyoto Encyclopedia of Genes and Genomes pathway, LGD likely gene-disruptive mutations, MM module membership, PCA principal component analysis, PPI protein-protein interaction, PSD postsynaptic density, ROC receiver operating characteristic, SVM-REF support vector machine recursive feature elimination, WGCNA weighted gene co-expression network analysis.

Fig. 2. Differential gene expression analysis between ASD and CTL samples and enrichment analysis (Fisher’s exact test).

A Volcano plot for DElncRNAs. B Volcano plot for DEmRNAs. Differentially expressed (DE) mRNAs and lncRNAs are highlighted in red. C PCA analysis. D Heatmap of DEmRNAs and DElncRNAs. E Heatmap showing enrichment of markers for different types of neural cells. F Heatmap showing enrichment of ASD risk genes from SFARI (ASD SFARI), intellectual disability genes (ID all), schizophrenia genes (SCZ), attention deficit hyperactivity disorder genes (ADHD), and Alzheimer’s disease genes (AD). “ASD/ID overlap,” the overlap between the “ASD SFARI” and “ID all” sets. “ASD only” and “ID only,” non-overlapping ASD SFARI and ID genes, respectively. G Heatmap shows the enrichment of genes affected by de novo variants (DNVs), including likely gene-disrupting (LGD), missense, synonymous in ASD (M, male; F, female), and genes encoding proteins in the postsynaptic density (PSD), genes expressed preferentially in human embryonic brains (Embryonic). Fisher’s exact test (two-tailed) with FDR correction was applied for enrichment. In (E–G), enrichment odds ratios (OR) and FDR-corrected p-values are shown for enrichment with FDR < 0.05. PCA principal component analysis, FC fold-change.

Fig. 3. Integrated enrichment analysis.

A Heatmap plot of BP enriched terms for the upregulated DEmRNAs. B Heatmap plot of BP enriched terms for the downregulated DEmRNAs. C, D Enrichment map of KEGG analysis for the up- and downregulated DEmRNAs. E, F Significant KEGG pathway and BP terms by GSEA analysis. G Significant KEGG pathway by GSVA analysis. BP biological process, KEGG Kyoto Encyclopedia of Genes and Genomes pathway, GSEA gene set enrichment analysis, GSVA gene set variation analysis.

Fig. 4. lncRNA and mRNA co-expression modules dysregulated in postmortem ASD cortex.

A Pearson’s correlation analysis between module eigengenes and different covariates (upper part). Correlation coefficients and p-values are shown at p < 0.05. The right side is named according to the BP of each module. The module enrichment analysis (Fisher’s exact test, FDR < 0.05) is shown on the lower part. Enrichment OR and FDR-corrected p-values are shown for enrichment with FDR < 0.05. B–F PPI network construction for five modules (M1, M3, M4, M8, and M10) was correlated with the disease status. ND neurodevelopment, BP biological process, PPI protein-protein interaction.

Fig. 5. ceRNA network, functional enrichment analysis, and validation of the ceRNA network.

A The preliminary lncRNA-miRNA-mRNA (ceRNA) regulatory networks. B Key ceRNA network. C, D KEGG enrichment analysis of preliminary or key ceRNA network. E Expression of 10 hub DEmRNAs and 1 hub DElncRNA in patients with autism (n = 70) and healthy controls (n = 75). Mann–Whitney U test was used for statistical analysis, and FDR was used for multiple testing corrections with each dot representing an individual. Data are presented as the median and interquartile range. * FDR < 0.05, ** FDR < 0.01, *** FDR < 0.001, **** FDR < 0.0001. F Validation of the expression of ceRNA network marker genes in four datasets. * |FC| > 1.3 and p-value < 0.05; x, the chip data cannot detect mRNA and lncRNA in the sample. G Experimental validation of a downregulated lncRNA (MIR600HG) and the corresponding mRNA regulatory axis. Data are presented as individual data points, with bar plots showing the mean and standard deviation (n = 11 for each group). Statistical significance was calculated by Student’s t-test, * P < 0.05, ** P < 0.01, **** P < 0.0001, Mean ± SEM.