Full Text

Introduction

Schizophrenia (SZ) is a neurodevelopmental disorder listed among the top 15 most burdening disabilities worldwide [1]. Despite decades of research, the etiology of the disease remains elusive due to its complexity, heterogeneity, and polygenicity. Genomic abnormalities may fractionally explain the substantial heritability of SZ but show limitations as diagnostic markers and etiology indicators due to a low fraction of explained disorder risk probability, thus suggesting a substantial role of epigenetic and environmental factors [2]. Gene expression analysis can bridge the gap between genomic and environmental risks, making it a promising approach to studying the pathophysiology of the disease.

Multiple regions in the brains of SZ patients display structural and functional abnormalities in neuroimaging studies, yet current molecular analysis remains restricted to a few selected brain areas. There is widespread cortical thinning in SZ individuals, with significant volumetric decreases in the frontal, temporal, and parietal lobes [3, 4]. Several subregions of these lobes, namely the dorsolateral prefrontal cortex and the superior temporal gyrus, are commonly selected for the transcriptomic profiling of psychiatric diseases [5, 6]. Similarly, a large-scale imaging study of subcortical structures revealed smaller hippocampus, amygdala, and thalamus [7]. Consequently, multiple gene expression studies investigated particular locations within these regions, particularly those with functional relevance to cognitive and emotional functions, but the findings were surprisingly inconsistent [8]. Furthermore, most of the studies focused on a single region, limiting our capability to unfold the molecular networks underlying such a multiplex disorder like SZ. Furthermore, SZ-associated transcriptome alterations have been found in other “not-as-popular” parts of the brain, such as the parietal lobe [9] and the cerebellum [10, 11]. Yet, these regions are virtually neglected in psychiatric research, leaving their mechanistic involvement in SZ pathology unknown.

Long non-coding RNAs (lncRNAs) play critical roles in gene expression regulation in the brain, with disruption of this regulation implicated in various mental disorders including SZ [12, 13]. Even though lncRNAs are usually synthesized by RNA polymerase II, similar to the protein-coding transcripts, they vastly exceed the mRNAs in terms of diversity, especially in nervous tissue [14]. Nonetheless, most of the expression analyses of post-mortem brains focused on protein-coding RNAs leaving the non-coding component of the transcriptome unexplored. Within the brain, many lncRNAs are specifically expressed in particular regions and at defined developmental stages [14, 15]. Therefore, alterations of lncRNA expression patterns could be affiliated with the disrupted developmental programming postulated by the neurodevelopmental hypothesis of SZ. A growing number of lncRNAs is documented to be regulators at multiple levels of gene expression affecting biological processes encompassing neuronal differentiation and immune response [15, 16]. Due to the widespread comorbidities and substantial overlap of behavioral symptoms among psychiatric illnesses, many candidate lncRNA regulators could be connected to more than one disease [16].

The largest class of lncRNA is long intergenic non-coding RNA (lincRNAs) which, in addition to the length and non-translated requirement of lncRNA, do not overlap with protein-coding sequences. Compared to mRNAs, lincRNAs are less conserved, less efficiently spliced, and more tissue-specific despite sharing similar biogenesis pathways [17]. Aside from the overlapping features, a few other aspects are supporting the distinction of lincRNAs from the other intragenic lncRNAs [18]. Herein, we identified lincRNAs associated with SZ and annotated their biological functions by comparing the transcriptomes of 35 anatomical regions corresponding to 10 anatomical sections in postmortem brains of four healthy and four SZ-diagnosed individuals (Table 1, Table S1).

Results

Region-dependent lincRNA expression differences cluster in white matter regions of the brain.

Previously, we published transcriptome data covering 33 regions in the healthy human brain[19]. Here, we analyzed gene expression in the same 33 regions in four individuals diagnosed with schizophrenia, and two additional brain regions in the diagnosed and control groups (Table 1). Based on these data, we detected the expression of 861 annotated human lincRNAs passing the intensity threshold. Visualization of expression variation using the principal component analysis (PCA) showed substantial overlap of healthy control (HC) and schizophrenia (SZ) samples (Figure 1A), while the segregation of the samples with regard to anatomical regions was more evident (Figure 1B). In accordance with this observation, we observed three clusters of brain regions produced by unsupervised hierarchical clustering based on lincRNA expression levels, aligning with the anatomical subdivision of the brain (Figure 1C, S1). The first cluster contains mainly the neocortical areas, the second – all the connective nerve tracts and the cerebellum, and the third – the diencephalon and most of the basal ganglia. Substructures of the limbic system belong to both clusters I and III, with the regions spatially related within the clusters.

Out of 861 detected lincRNAs, the expression of 135 depended significantly on both conditions and brain regions (two-way ANOVA, BH-adjusted p < 0.05 for the interaction term, excluding genes with unequal variance). Among these lincRNAs, we identified 37 differentially expressed lincRNAs (DELs) showing significant differences between HC and SZ individuals in at least one region (Figure 2A, Table S2). Of them, four DELs were dysregulated in two regions, while the rest were dysregulated in one. Further, 31 of the 37 DELs showed a two-fold or greater expression level difference between the two conditions (Table S2). Notably, these significant differences were not distributed uniformly within the brain but were associated with 10 out of 35 examined regions. Most of the associations were found in the regions containing connective axonal tracks: cerebellar white matter contained two down- with 14 up-regulated lincRNAs and three regions of the cerebral white matter contained 13 down- with seven up-regulated lincRNAs (Figure 2B).

The three DELs showing the most significant difference between SZ and HC expression in our study included: MEG3, RP11-247L20.4, and LINC01252. These lincRNAs were all dysregulated in the white matter of the cerebellum, and RP11-247L20.4 was also significantly downregulated in the cerebellar gray matter. MEG3 was previously reported to be differentially expressed in the hippocampus [20], superior temporal gyrus [21], and amygdala [22, 23] of SZ patients. Similarly, LINC01252 was reported to be upregulated in SZ in the amygdala[24]. Two other DELs that appeared in related literature included MEG9 dysregulation in the amygdala [22] and MALAT1 – in the dorsolateral prefrontal cortex [25]. We illustrated the expression profiles of these five genes in Figure 2C. Notably, MALAT1 was also the gene with the biggest difference amplitude among DELs, showing a 16-fold increase in the cerebellum of SZ patients compared to HC individuals. The direction of effect reported in the published studies coincided with the one found in our analysis, with the sole exception of the MEG3 expression difference in the hippocampus, where the difference was not statistically significant in our study. While most SZ gene expression studies did not focus on lincRNA expression, three studies contained lincRNA datasets. Re-analysis of these data revealed a positive and significant correlation of SZ-associated fold changes between our and published lincRNA data in the amygdala (two studies, Spearman correlation test, ρ > 0.45, p < 0.03) [22, 24] but not in the dorsolateral prefrontal cortex (one study, Spearman correlation test, ρ = 0.15, p = 0.39) [25]. The absence of statistically significant agreement in the prefrontal cortex could be due to insufficient power of the comparison and, more importantly, lack of substantiation of lincRNA expression alterations in this region. Our general analysis, as well as expression profiles of the five selected DELs, show the concentration of large-amplitude SZ-associated expression differences in white matter regions and cerebellar gray matter, with some significant differences also found in the amygdala, but none in the prefrontal cortex).

Functional annotation of DELs links them to neuroplasticity and neurotransmission.

While human protein-coding genes tend to have substantial functional annotation, this is not the case for the vast majority of lincRNA. Nonetheless, the co-expression of mRNA and lincRNA transcripts could indicate the functional roles of non-coding counterparts. To perform such an annotation, we tested the correlations of DEL expression difference profiles recorded across the 35 brain regions (for example, see profiles of five selected DELs in Figure 2C) with the expression profiles of mRNA extracted from the same data. The mRNAs strongly correlated with a DEL expression difference profile (Pearson’s r ≥ 0.85) were defined as targets of the respective DEL. Whilst most DELs had few or no targets, three DELs correlated with outstanding numbers of mRNAs, which altogether constituted 218 out of 231 identified DEL-mRNA correlations. These three DELs included RP11-74E22.8 lincRNA upregulated in the cerebellar white matter, and LINC01963 and RP11-416I2.1 both downregulated in the internal capsule, also a white matter region. Furthermore, the profiles of RP11-74E22.8 and LINC01963 were strongly positively correlated (Pearson r = 0.76) and they shared 17 common mRNA targets out of 71 and 105 targets, respectively.

            The mRNA targets of the three DELs were significantly associated with the neuronal activity terms listed in the Gene Ontology (GO) database (hypergeometric test, BH-adjusted p-values < 0.05). Specifically, mRNA targets of RP11-416I2.1 were associated with voltage-gated channels and neuroplasticity, while targets of the other two DELs were linked to synaptic transmission and signaling terms (Figure 3A, S2A, S2B). Analysis of mRNA target enrichment using another functional annotation database, the Kyoto Encyclopedia of Genes and Genomes (KEGG), linked four pathways, including “synaptic vesicle cycle” and “metabolism of alanine, aspartate, and glutamate” with targets of LINC01963 (Figure S2B). Functional annotation using Reactome Knowledgebase [26] yielded similar results: the targets of RP11-416I2.1 were enriched in potassium channels and G-protein coupled receptors, while the other two target groups shared mutual entries related to the neurotransmitter release cycle (Figure S2D).

            We further investigated the association of DEL target mRNAs with eight main brain cell types by testing them using a customized list of markers genes extracted from publications. The analysis revealed an evident and significant association between the targets of all three target-rich DELs and general neuronal markers, as well as markers of excitatory neurons (Figure 3B). This result aligns with the functional annotation outcome dominated by terms related to neuronal functionality. In addition, the targets of RP11-74E22.8 overlapped significantly with inhibitory neuron markers, which might be related to their enrichment in GABA, dopaminergic, and norepinephrine pathways. These results suggested that the three DELs could modulate a network of genes expressed in neurons and involved in synaptic signal transduction.

Discussion

Our analysis of lincRNA expression in the SZ patients’ brains revealed few alterations in the cerebral cortex and basal ganglia regions commonly thought to be associated with the disorder. Instead, our analysis shows substantial lincRNA dysregulation in the cortical white matter regions and the cerebellum. Our study deferrers from most of the previous SZ brain expression analyses in two substantial aspects. First, by measuring gene expression in multiple regions of the same brain we base our analysis on the expression profiles of transcripts within the brain, thus minimizing the interindividual variation. This approach allows us to focus on expression differences particular to specific brain regions, including the ones neglected by previous studies. Interindividual variation poses a serious problem in human studies due to uncontrollable and diverse genetic and environmental factor effects, resulting in the loss of biologically meaningful differences with marginal to modest effect sizes [27]. Most existing gene expression studies of the SZ brain focused on either a single or few regions of the cerebrum [8]. However, sporadic omics screening of the neglected brain regions, such as the transcriptome and proteome assessment of the cerebellum [10, 11], have reported meaningful expression alterations. In this study, focusing on the expression alteration patterns recorded across the 35 brain regions, we show that the white matter and the cerebellum might need to receive more attention in future SZ studies.

The main component of the white matter is myelinated axons extending from the neuronal cell bodies. Thus, it might seem unusual that observed lincRNA expression alterations were not accompanied by changes in the corresponding gray matter. One hypothesis is that changes in gene expression of these regions arose mainly from local glial cells. This notion, however, is not likely for the cerebellum gray matter given that the non-neuronal cells account for less than one-fifth of the total cell population in this structure [28]. Thus, gene dysregulations in glial cells have to be immense in order to explain the observed differences. Another possible explanation of observed differences is the redistribution of transcripts leading to the accumulation of the DELs in axons, possibly due to molecular transport impairment. This explanation aligns with the white matter pathology of SZ [29] and could be linked to a transcriptomic study of the cerebellum reporting dysregulation of genes involved in Golgi function and presynaptic vesicular transport in SZ [11]. Alternatively, the disconnected patterns of white matter regions and the cerebellum from the rest of the brain could be a consequence of cell-type-specific expression. There are many nuclei in the cerebellar and cerebral cortices, including all of the cerebellar granular layer neurons, that do not project into the other brain regions. Thus, the transcriptomes of such non-projecting neurons might not be reflected in the white matter transcriptome. Future experiments should investigate the cellular locations of these molecules to reveal the underlying mechanism and clarify these speculations.

The second particular aspect of our study is its focus on non-coding RNA expression. Unlike mRNA, lincRNA expression displays more pronounced tissue and brain region specificity, as well as greater response amplitudes, making them better perspective markers of disorder-related alterations[30]. On the other hand, the evident drawback of the lincRNA research is an almost complete lack of functional annotation, hindering the results’ interpretation. In our study, however, we were able to largely overcome this limitation by taking advantage of lincRNA-mRNA co-expression analysis relying on the transcript profiles measured across the 35 brain regions. The reliance on these profiles, instead of variation-prone inter-individual comparisons, allowed us to unambiguously assign co-expressed lincRNA targets to neurons and neuron-specific functions, such as synaptic signal transduction. It is also noteworthy that a substantial fraction of differentially expressed lincRNAs identified in our study (DELs) overlapped with lincRNAs previously reported by the few corresponding analyses. Specifically, the DEL with the lowest p-value in our study, MEG3, was also differentially expressed in various brain regions [20–23] as well as peripheral blood tissue [31] of SZ patients. This gene is proposed to regulate long-term potentiation via modulating a specific signaling cascade [32] and displayed a significant downregulation in SZ patients treated with risperidone, an antipsychotic medication [31]. The DEL showing the greatest fold change, MALAT1, was reported in another transcriptomic study of the dorsolateral prefrontal cortex in SZ [25]. These primitive findings are supporting primitive evidence for more confirmative research focusing on disease-associated non-coding RNA to be conducted.

The main limitation of our work is the low count of investigated brains per group. Although our study included 35 regions from each individual and the number of samples per group was balanced, we had only four biological replicates in each group. This number is certainly low for such a heterogeneous disease as SZ. The clinical presentation of SZ is extremely diverse, so whether the subjective, behavior-based diagnosis of the disease agrees with the pattern of molecular alterations remains a topic of controversy [33]. However, this limitation is universal and represents a problem largely unresolved in most post-mortem brain studies [2, 8]. Unlike other studies, our analysis minimizes inter-individual variation by using the average expression level of each transcript within a given brain as an internal control. As a result, despite the limited sample size, we identified numerous differentially expressed lincRNAs whose involvements in SZ could be supported by previous SZ studies, as well as their evident neuronal-specific functionality suggesting biologically meaningful signals.

Conclusion

The decades of research showed that complex systemic brain disorders, such as SZ, require systematic analysis of associated alterations. Our analysis of long non-coding RNA expression patterns across 35 diverse brain regions supports this notion, revealing the clustering of SZ-associated expression alterations in brain structures commonly neglected by transcriptome studies: white matter and cerebellar brain regions. Due to their high expression specificity, long non-coding RNAs represent prospective markers of functional alterations associated with SZ. Indeed, despite their prevalence in white matter, identified lincRNA expression alterations showed stark association with neurons and neuron-specific functions, such as synaptic transmission.

Material and Methods

Brain samples

Our study used samples dissected from eight frozen human brains, four healthy ones (HC) obtained from the Chinese Brain Bank Center and four schizophrenia ones (SZ) obtained from the Moscow Psychiatric Hospital No. 1 (Table S1). Diagnosis of SZ was made by experienced clinicians. Healthy controls had no history of psychiatric diseases or brain-related disputations. From each brain, we dissected 35 regions, listed in Table 1, without thawing. The dissected specimens were then preserved at -80^oC until RNA extraction.

RNA library preparation and sequencing data assessment

From each sample, RNA was extracted from an approximately 30-milligram block following the protocol of RNA extraction with QIAzol Lysis Reagent. After RNA integrity and concentration were measured, sequencing libraries were prepared following the manufacturer’s recommendations with poly-A selection. The libraries were sequenced using Illumina Hiseq 4000 platform.

FastQC [34] was used to assess the quality of raw reads. We used Trimmomatic [35] to remove low-quality reads and all adapters identified previously or provided by Trimmomatic. We mapped the reads to the human reference genome GRCh38 with HISAT2 [36]. The gene count matrix was retrieved using Stringtie [37]. Gene annotations were obtained from Ensembl v91.

Differential gene expression analysis

Genes identified as long intergenic non-coding RNAs (lincRNAs) and protein-coding (mRNAs) according to Ensembl annotation were chosen for downstream analysis. For both the control and disease cohorts, only genes with no more than two zero coverage values out of four replicates in all regions were retained (two-zero threshold). The count data were transformed to the logarithmic scale using the package DESeq2 [38]. We adjusted for the sample quality by regressing the expression levels using the RIN (RNA Integrity Number) values. Finally, we used donor-centered normalization by subtracting the means of log-transformed expression values within each brain from the regional expression values of the respective individual. Based on the regional means of pooled HC and SZ groups, we identified three clusters of brain regions using the hierarchical clustering method, in which distances were calculated as one minus Pearson correlation coefficients and clusters were defined by Ward’s linkage function.

For the lincRNAs dataset, Levene’s test was used to exclude genes with high heteroscedasticity (the threshold for exclusion was p < 0.05). Two-way ANOVA including the effect of region and conditions (diagnosis) was used to identify transcripts with significant differences for the interaction terms (p < 0.05 after Benjamini-Hochberg correction). These transcripts were chosen for the post hoc Tukey test. Differentially expressed lincRNAs (DELs) were determined as those with at least one significant difference between the two conditions of the same region in the Tukey test. The fold changes were measured by subtracting the mean transformed expression of the control groups from that of the SZ groups.

Enrichment analysis

We calculated the Pearson correlation between the normalized expression values of DELs and all mRNAs passing the two-zero threshold described above. Protein coding genes with a correlation coefficient r ≥ 0.85 were defined as targets of the respective DELs. Groups with at least 10 targets were used for downstream analyses. For each DEL, we analyzed the enrichment of its targets using all coding genes passing the detection threshold as the background set. Gene Ontology (GO) terms, Kyoto Encyclopedia of Genes and Genomes (KEGG), and Reactome pathways analyses for all groups were implemented using the clusterProfiler R package [39]. Terms and pathways with adjusted p-value < 0.05 were considered enriched.

We customized a set of marker genes for eight main brain cell types based on reported gene sets. For excitatory and inhibitory neurons, we chose intersections of the corresponding sets from the studies in [19, 40]. Similarly, other neuron markers, microglia, astrocytes, and oligodendrocytes represented the respective intersections of lists reported in [19, 41, 42]. Markers for oligodendrocyte progenitor cells were reported in [19] and the ones for endothelial cells were extracted from [40]. We tested the overrepresentation of these markers in the list of DEL targets using Fisher’s exact test followed by Benjamini-Hochberg correction.

About the authors

Tuan Nguyen

Email: tuan.nguyen@skoltech.ru

Olga Igorevna Efimova

Email: o.efimova@skoltech.ru

Artem Valerievich Tokarchuk

Email: a.tokarchuk@skoltech.ru

Yana Alexandrovna Zorkina

Email: zorkina.ya@serbsky.ru

Anna Yurievna Morozova

Email: hakurate77@gmail.com

Denis Sergeevich Andreyuk

Author for correspondence.
Email: denis.s.andreyuk@yandex.ru

George Petrovich Kostyuk

Email: editorinchief@consortium-psy.com

Philipp Efi movich Khaitovich

Email: p.khaitovich@skoltech.ru

References

Supplementary files