Abstract

   Schizophrenia (SCZ) is a serious psychiatric condition with a 1%
   lifetime risk. SCZ is one of the top ten global causes of disabilities.
   Despite numerous attempts to understand the function of genetic factors
   in SCZ development, genetic components in SCZ pathophysiology remain
   unknown. The competing endogenous RNA (ceRNA) network has been
   demonstrated to be involved in the development of many kinds of
   diseases. The ceRNA hypothesis states that cross-talks between coding
   and non-coding RNAs, including long non-coding RNAs (lncRNAs), via
   miRNA complementary sequences known as miRNA response elements, creates
   a large regulatory network across the transcriptome. In the present
   study, we developed a lncRNA-related ceRNA network to elucidate
   molecular regulatory mechanisms involved in SCZ. Microarray datasets
   associated with brain regions ([38]GSE53987) and lymphoblasts (LBs)
   derived from peripheral blood (sample set B from [39]GSE73129) of SCZ
   patients and control subjects containing information about both mRNAs
   and lncRNAs were downloaded from the Gene Expression Omnibus database.
   The [40]GSE53987 comprised 48 brain samples taken from SCZ patients (15
   HPC: hippocampus, 15 BA46: Brodmann area 46, 18 STR: striatum) and 55
   brain samples taken from control subjects (18 HPC, 19 BA46, 18 STR).
   The sample set B of [41]GSE73129 comprised 30 LB samples (15 patients
   with SCZ and 15 controls). Differentially expressed mRNAs (DEmRNAs) and
   lncRNAs (DElncRNAs) were identified using the limma package of the R
   software. Using DIANA-LncBase, Human MicroRNA Disease Database (HMDD),
   and miRTarBase, the lncRNA- associated ceRNA network was generated.
   Pathway enrichment of DEmRNAs was performed using the Enrichr tool. We
   developed a protein–protein interaction network of DEmRNAs and
   identified the top five hub genes by the use of STRING and Cytoscape,
   respectively. Eventually, the hub genes, DElncRNAs, and predictive
   miRNAs were chosen to reconstruct the subceRNA networks. Our
   bioinformatics analysis showed that twelve key DEmRNAs, including BDNF,
   VEGFA, FGF2, FOS, CD44, SOX2, NRAS, SPARC, ZFP36, FGG, ELAVL1, and
   STARD13, participate in the ceRNA network in SCZ. We also identified
   DLX6-AS1, NEAT1, MINCR, LINC01094, DLGAP1-AS1, BABAM2-AS1, PAX8-AS1,
   ZFHX4-AS1, XIST, and MALAT1 as key DElncRNAs regulating the genes
   mentioned above. Furthermore, expression of 15 DEmRNAs (e.g., ADM and
   HLA-DRB1) and one DElncRNA (XIST) were changed in both the brain and
   LB, suggesting that they could be regarded as candidates for future
   biomarker studies. The study indicated that ceRNAs could be research
   candidates for investigating SCZ molecular pathways.

   Subject terms: Computational biology and bioinformatics, Genetics,
   Psychology

Introduction

   Schizophrenia (SCZ) is a severe psychotic disorder with a lifetime risk
   of almost 1%. SCZ is considered one of the ten most common causes of
   disability throughout the world^[42]1. The symptoms of SCZ are usually
   classified into three categories: positive symptoms such as
   hallucinations, delusions, and thinking problems; negative symptoms
   such as diminished emotions, reduced motivation, decreased speaking,
   and social deterioration; cognitive symptoms such as attention,
   emphasis, processing speed, and memory dysfunctions^[43]1,[44]2. The
   etiology of SCZ is demonstrated to be multifactorial, in which both
   environmental and genetic factors have been identified to participate
   in the manifestation of symptoms^[45]3. Dysregulated gene expression
   and protein production are considered to be associated with the
   pathophysiology of SCZ and have spatial variation across brain regions
   and temporal variation during disease progression^[46]4–[47]6. A recent
   postmortem transcriptional profiling study shows a considerable burden
   of differentially expressed genes (DEGs) across hippocampus (HPC),
   Brodmann area 46 (BA46), and striatum (STR) regions and reveals that
   multiple signaling and inflammatory pathways are dysregulated in all
   three regions in SCZ^[48]7. On the other hand, a growing body of
   evidence has indicated that the brain is linked to the periphery via
   the circulatory system containing hormones and other secreted
   regulatory molecules produced in the diffuse neuroendocrine system
   which impact the gene expression pattern of the peripheral
   markers^[49]8,[50]9. These findings are confirmatory evidence that SCZ
   is a systemic disorder and support the notion that investigation of the
   SCZ can be through studying the biomarkers in the peripheral samples
   such as whole blood, peripheral blood mononuclear cells (PBMCs),
   lymphoblasts (LBs), and olfactory epithelium^[51]10–[52]13. Numerous
   studies demonstrated alterations in non-coding RNA (ncRNA) in people
   with SCZ^[53]14. These findings provided information concerning the
   underlying mechanisms of dysregulation in gene expression and protein
   production. The ncRNA entails several RNA molecule families with a
   genetic distinction based on mature sequence size^[54]14. A growing
   body of evidence suggests that microRNAs (miRNAs) (20–22 nucleotides)
   are aberrantly expressed in the brains of SCZ patients^[55]15,[56]16.
   Furthermore, new research indicates that the aberrant expression of
   other ncRNAs, including long non-coding RNAs (lncRNAs) with more than
   200 nucleotides, could be implicated in SCZ^[57]17,[58]18.

   Remarkable data have been provided about the molecular mechanisms and
   biological procedures of SCZ through advanced genetic, genomic,
   developmental neurobiology, and systems biology investigations.
   High-throughput discovery of genes that contribute to the molecular
   pathology of SCZ is now possible using microarray technologies and
   next-generation sequencing^[59]7,[60]13,[61]19,[62]20. The role of
   genetic in SCZ pathophysiology is still vague, and newly discovered
   mechanisms responsible for regulating gene transcription could be
   beneficial in resolving how these alterations cause the development and
   progression of SCZ. As a result, potencies for more efficient treatment
   and diagnosis are found. Pier Paolo Pandolfi and colleagues recommended
   a novel regulatory mechanism, known as competing endogenous RNA
   (ceRNA), in 2011^[63]21. It was suggested that cross-talk between both
   coding and non-coding RNAs such as lncRNAs, circular RNAs (circRNAs),
   and pseudogenes via miRNA response elements (MREs), as miRNA
   complementary sequences, leads to the formation of a large-scale
   regulatory network in diverse parts of the transcriptome^[64]22.
   According to the ceRNA hypothesis, RNAs transcripts can regulate each
   other by competing for shared miRNAs. In addition, the expression
   levels of these two RNA transcripts will have a positive relationship
   with each other^[65]21. Over the last few years, numerous studies have
   validated the ceRNA theory. It is well known that disruption of the
   balance of ceRNA cross-talk is associated with various developmental
   processes and pathological conditions, including tumorigenesis,
   neurodegenerative diseases (e.g., Alzheimer’s disease (AD))^[66]23, and
   mental disorders (e.g., SCZ and depression)^[67]24,[68]25.

   Given that SCZ is one of the ten most prevalent causes of disability
   globally and since the nature of the genetic component in SCZ
   pathophysiology is not fully understood, we performed a bioinformatics
   analysis to elucidate the lncRNA-associated ceRNA network to clarify
   molecular regulatory mechanisms involved in SCZ.

Methods

   In the present study, we employed a system biology approach to mine
   data of the two microarray datasets of multiple regions of the brain
   ([69]GSE53987) and the LBs derived from peripheral blood (sample set B
   from [70]GSE73129) from subjects with SCZ and matched controls. We
   intended to recognize differentially expressed mRNAs (DEmRNAs) and
   lncRNAs (DElncRNAs) and make a lncRNA-associated ceRNA network. All
   methods were carried out in accordance with relevant guidelines and
   regulations.

Gene expression profile data collection

   The gene expression profiles mentioned above were obtained from the
   NCBI Gene Expression Omnibus database (GEO,
   [71]https://www.ncbi.nlm.nih.gov/geo/). A chip-based platform
   [72]GPL570 (HG-U133_Plus_2) Affymetrix Human Genome U133 Plus 2.0 Array
   was applied for both datasets. The [73]GSE53987 included 48 brain
   samples from SCZ patients (15 HPC, 15 BA46, 18 STR) and 55 brain
   samples from control subjects (18 HPC, 19 BA46, 18 STR)^[74]7. The
   sample set B from [75]GSE73129 contained 30 LB samples, of which 15
   were SCZ patients, and 15 were control^[76]13.

Data preprocessing and DEmRNAs and DElncRNAs identification

   We used Robust Multichip Average (RMA) for background correction, and
   quantile normalization of all the raw data files^[77]26. Brain regions
   and LB samples were normalized separately. An interquartile range (IQR)
   filter (IQR across the samples on the log base two scale greater than
   median IQR) was used to reduce the number of evaluated genes followed
   by an intensity filter (a minimum of > 100 expression signals in a
   minimum of 25% of the arrays) aiming to remove the non-significant
   probe sets which are not expressed or not changing^[78]27. AgiMicroRna
   Bioconductor package was used for assessing the quality control^[79]28.
   We employed principal component analysis (PCA) to conduct a dimensional
   reduction analysis^[80]29, aiming to find similarities between each
   sample group using R software’s ggplot2 package^[81]30. Besides, a
   heatmap was drawn to show the correlation between samples using the
   Pheatmap package of R^[82]31. Differential expression gene analysis
   (DEGA) was done between SCZ and normal samples using the linear models
   for microarray data (limma) R package^[83]32 in Bioconductor
   ([84]https://www.bioconductor.org/)^[85]33. We used the
   removeBatchEffect() function from the limma package to remove batch
   effects. Age, gender, and race were used as covariates in [86]GSE73129
   study and race and post-mortem interval in [87]GSE53987 study to modify
   the data for confounding variables. We did a linear model on the main
   factor “group” (disease vs. control) using limma’s lmFit() function and
   included covariates in the linear model as well. The limma’s eBayes ()
   function was then used to calculate differentially expressed data from
   the linear fit model. We utilized the previously used approach to
   identify lncRNA probes^[88]34. We downloaded the complete list of
   lncRNA genes with approved HUGO Gene Nomenclature Committee (HGNC)
   symbols from ([89]https://www.genenames.org/)^[90]35. Then, we compared
   the lncRNA gene list with our dataset gene symbols and chose the
   overlapped genes. Student t-test was used and we set the aberrantly
   expressed RNAs cut-off as: (1) a false discovery rate (adjusted
   P-value) < 0.001, and (2) |log2 fold change (log2FC)|> 0.69. Note that
   we discarded probes with no gene symbols.

lncRNA-associated ceRNA network construction

   DIANA-LncBase v3^[91]36 was used to identify the experimentally
   validated interactions between the miRNAs and lncRNAs. Homo Sapiens
   “Species”, Brain or Peripheral Blood “Tissue”, and high “miRNA
   Confidence Levels” were selected as criteria for the DIANA-LncBase
   query. SCZ-related miRNAs were collected from the Human microRNA
   Disease Database (HMDD) v3.2 database^[92]37. Besides, we obtained the
   interactions between miRNAs that were collected using the HMDD database
   and target mRNAs supported by strong experimental evidence from
   miRTarBase^[93]38. These mRNAs were then compared with the previously
   obtained mRNAs. Then, we used the duplicated mRNAs to construct the
   lncRNA-miRNA-mRNA network. LncRNAs, targeted mRNAs, and the interacted
   miRNAs were removed from the ceRNA network in the opposite expression
   pattern between the targeted mRNAs and lncRNAs. Cytoscape software
   (version 3.8.0)^[94]39 was used to generate ceRNA regulatory network.

Pathway enrichment analysis of DEmRNAs

   The Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis
   was performed using the Enrichr tool^[95]40,[96]41 for pathway
   enrichment to analyze the DEmRNAs that were in the ceRNA network.

Protein–protein interaction (PPI) network construction and hub genes
identification

   A PPI network was constructed using online STRING database
   ([97]https://string-db.org/)^[98]42 for predicting the interactive
   relationships of common DEmRNAs encoding proteins. A minimum
   interaction score of above 0.7 (high confidence) was needed for PPI
   network construction. Cytoscape software was employed for PPI network
   visualization and hub genes analysis. For PPI network simplification,
   we removed the non-interacting genes. Then, we considered the top five
   genes with the highest degree of connection to the others as the hub
   genes according to CytoHubba (version 0.1) analysis^[99]43 from
   Cytoscape.

Reconstruction of the subceRNA network

   We retrieved the corresponding lncRNA and miRNA from the original
   network of ceRNA. We used them in the reconstruction of the subceRNA
   network to confirm the ceRNA network of hub genes.

Results

DEmRNAs and DElncRNAs identification

   Background correction, normalization, gene filtration, and batch
   correction were performed prior to DEGA. We applied AgiMicroRna
   Bioconductor package for quality control assessment. Box plots of gene
   expression data were depicted to evaluate data distribution following
   normalization (Supplementary Fig. [100]S1). In the box plots, separate
   arrays had similar medians of expression level, confirming that
   correction was carried out properly. We checked sample correlation at
   the probe level. The samples were divided into two main clusters, one
   of which had SCZ, and the other consisted of the control samples.
   Moreover, a PCA plot was used to demonstrate the spatial distribution
   of samples (Fig. [101]1). PCA presents information about the structure
   of the analyzed data. It helps find similarities between samples. One
   of the STR samples in the SCZ group was removed due to being spatially
   far from other SCZ samples. The number of DEmRNAs and DElncRNAs are
   listed in Supplementary Table [102]S1. Our analysis with the limma
   package in R identified 249 DEmRNAs (228 upregulated and 21
   downregulated) (e.g., BDNF: brain-derived neurotrophic factor, VEGFA:
   vascular endothelial factor A, FGF2: fibroblast growth factor 2, FOS,
   and CD44), eight DElncRNAs (six upregulated and two downregulated)
   (e.g., DLX6-AS1: DLX6 antisense RNA 1, NEAT1: nuclear-enriched abundant
   transcript 1, MINCR: MYC-induced long non-coding RNA, LINC01094: long
   intergenic non-protein coding RNA 1094, DLGAP1-AS1: DLGAP1 antisense
   RNA 1, BABAM2-AS1: BABAM2 antisense RNA 1, and PAX8-AS1: PAX8 antisense
   RNA 1) in HPC region, 286 DEmRNAs (271 upregulated and 15
   downregulated) (e.g., VEGFA, FGF2, CD44, SOX2: SRY-box transcription
   factor 2, and NRAS), five DElncRNAs (five upregulated) (including,
   NEAT1, MINCR, LINC01094, PAX8-AS1, and ZFHX4-AS1: ZFHX4 antisense RNA
   1) in BA46 region, 271 DEmRNAs (250 upregulated and 21 downregulated)
   (e.g., VEGFA, SPARC: secreted protein acidic and rich in cysteine,
   ZFP36, FGG: fibrinogen gamma chain, and ELAVL1: ELAV-like protein 1),
   seven DElncRNAs (six upregulated and one downregulated) (e.g., NEAT1,
   MINCR, LINC01094, DLGAP1-AS1, PAX8-AS1, and XIST: X-inactive specific
   transcript) in STR region, 159 DEmRNAs (39 upregulated and 120
   downregulated) (e.g., STARD13: StAR-related lipid transfer protein 13),
   and two lncRNAs (two upregulated) (e.g., MALAT1: metastasis associated
   lung adenocarcinoma transcript 1) in LB. Among the genes mentioned
   above, only BDNF and DLX6-AS1 genes had decreased expression and others
   showed increased expression. The complete lists of DEGs are shown in
   Supplementary Tables [103]S2–[104]S5. Moreover, the Venn diagram
   indicated 15 DEmRNAs (ADM, CLEC2B, F13A1, HLA-DQB1, HLA-DRB1, HLA-DRB3,
   HLA-DRB4, HLA-DRB5, CSF2RB, RHPN2, DDX3Y, EIF1AY, KDM5D, RPS4Y1, USP9Y)
   and one DElncRNAs (XIST) were dysregulated in both brain and LB. The
   details of these overlapped genes are shown in Fig. [105]2 and Table
   [106]1.

Figure 1.

   [107]Figure 1
   [108]Open in a new tab

   Hierarchically clustered correlation heatmap and principal component
   analysis (PCA) for the brain dataset ([109]GSE53987) are shown in (a)
   and (b), and for the lymphoblast dataset ([110]GSE73129) are depicted
   in (c) and (d), respectively. Each column in the heatmap represents one
   sample and shows the correlation to all samples (including itself),
   with red for correlation = 1 and blue for the lowest observed
   correlation. All SCZ and CTL samples correlated well with each other,
   except one SCZ sample from the [111]GSE53987 study marked by an
   asterisk (*). In PCA, all samples are segregated by condition group (on
   PC1). [112]GSM1305052 (a patient striatum sample) was removed from
   further analysis in order to its wrong spatial enrichment. This figure
   was made using R version 4.0.3 ([113]https://www.r-project.org/). CTL,
   control; SCZ, schizophrenia.

Figure 2.

   [114]Figure 2
   [115]Open in a new tab

   Analysis of differentially expressed genes (DEGs) in Schizophrenia
   (SCZ) brain and lymphoblast samples by Venn. Co-expression of (a)
   DEmRNAs and (b) DElncRNAs in different brain regions and lymphoblast of
   SCZ patients. This figure was made using R version 4.0.3
   ([116]https://www.r-project.org/). BA46, Brodmann area 46; HPC,
   hippocampus; LB, lymphoblast; STR, striatum.

Table 1.

   List of DEmRNAs and DElncRNAs that were dysregulated in both brain and
   blood.
   Samples DEmRNAs DElncRNA
   HPC and LB ADM, CLEC2B, F13A1, HLA-DQB1, HLA-DRB1, HLA-DRB3, HLA-DRB4,
   HLA-DRB5 –
   BA46 and LB ADM, CSF2RB, F13A1, HLA-DQB1, HLA-DRB1, HLA-DRB3, HLA-DRB4,
   HLA-DRB5, RHPN2 –
   STR and LB ADM, CLEC2B, CSF2RB, DDX3Y, EIF1AY, F13A1, HLA-DQB1,
   HLA-DRB1, HLA-DRB3, HLA-DRB4, HLA-DRB5, KDM5D, RPS4Y1, USP9Y XIST
   HPC, BA46, and LB ADM, F13A1, HLA-DQB1, HLA-DRB1, HLA-DRB3, HLA-DRB4,
   HLA-DRB5 –
   HPC, STR, and LB ADM, CLEC2B, F13A1, HLA-DQB1, HLA-DRB1, HLA-DRB3,
   HLA-DRB4, HLA-DRB5 –
   BA46, STR, and LB ADM, CSF2RB, F13A1, HLA-DQB1, HLA-DRB1, HLA-DRB3,
   HLA-DRB4, HLA-DRB5 –
   HPC, BA46, STR, and LB ADM, F13A1, HLA-DQB1, HLA-DRB1, HLA-DRB3,
   HLA-DRB4, HLA-DRB5 –
   [117]Open in a new tab

   BA46 Brodmann area 46, DElncRNA differentially expressed lncRNA,
   DEmRNAs differentially expressed mRNAs, HPC hippocampus, LB
   lymphoblast, STR striatum.

lncRNA-associated ceRNA network construction

   The DIANA-LncBase v3 online tool was used to anticipate DElncRNA-miRNA
   interaction pairings. DELncRNAs, targeted DEmRNAs, and also the
   interacted miRNAs were deleted from the ceRNA network in the opposite
   expression pattern present between DElncRNAs and the targeted DEmRNAs.
   Following that, miRTarBase was utilized to identify interactions
   between target mRNAs and SCZ-associated miRNAs (discovered by the
   HMDD). We developed a regulatory ceRNA network based on
   DElncRNA-miRNA-DEmRNA interactions to explore the potential mechanism
   behind SCZ pathophysiology (Fig. [118]3 and Table [119]2). CeRNA
   network contained 50 nodes (seven lncRNAs, 17 miRNAs, and 26 mRNAs) and
   95 edges in HPC, 54 nodes (five lncRNAs, 19 miRNAs, and 30 mRNAs) and
   101 edges in BA46, 54 nodes (seven lncRNAs, 19 miRNAs, and 28 mRNAs)
   and 97 edges in STR, and three nodes (one lncRNA, one miRNA, and one
   mRNA) and two edges in LB. There was no common ceRNA axis between LB
   and brain regions. We identified one ceRNA axis in LB
   (MALAT1/hsa-miR-9-5p/STARD13). Among the list, 90 different ceRNA axes
   were common between all three areas of the brain. These ceRNA axes are
   bolded in Table [120]1.

Figure 3.

   [121]Figure 3
   [122]Open in a new tab

   The long non-coding RNA-associated competing endogenous RNA (ceRNA)
   network. ceRNA axes in (a) hippocampus, (b) Brodmann area 46, (c)
   striatum regions and (d) lymphoblasts. The red and blue nodes represent
   the upregulation and downregulation, respectively. Gray edges represent
   interactions between RNAs. This figure was generated using Cytoscape
   version 3.8.0 ([123]https://cytoscape.org/)^[124]39. LncRNAs, miRNAs,
   and mRNAs are represented by hexagon, round rectangle, and ellipse,
   respectively.

Table 2.

   Details of competing endogenous RNA (ceRNA) axes.
   DElncRNA(s) Shared miRNA(s) DEmRNA(s) Sample Expression of DElncRNA(s)
   and DEmRNA(s)
   DLX6-AS1 hsa-miR-22-3p ACVR1C, BDNF HPC Down
   DLX6-AS1

   hsa-miR-124-3p

   hsa-miR-132-3p

   hsa-miR-15a-5p
   BDNF HPC Down
   DLX6-AS1 hsa-miR-708-5p NNAT HPC Down
   BABAM2-AS1, NEAT1, LINC01094 hsa-miR-30b-5p BCL6 HPC Up
   NEAT1 hsa-miR-9-5p BCL6, TGFBI, TGFBR2 HPC Up
   LINC01094, NEAT1 hsa-miR-132-3p CDKN1A HPC Up
   NEAT1 hsa-miR-22-3p CDKN1A HPC Up
   BABAM2-AS1, LINC01094, NEAT1 hsa-miR-17-5p CDKN1A, TGFBR2, VEGFA HPC Up
   NEAT1, LINC01094 hsa-miR-29a-3p FGG, TET1, VEGFA, ZFP36, TNFAIP3, MCL1,
   CD93 HPC Up
   LINC01094, NEAT1 hsa-miR-7-5p FOS HPC Up
   NEAT1, LINC01094 hsa-miR-29c-3p MCL1, FGG, VEGFA HPC Up
   LINC01094, MINCR, NEAT1 hsa-miR-23a-3p MT2A, HMGB2 HPC Up
   NEAT1 hsa-miR-126-3p NFKBIA, PIK3R2, ADM, VEGFA HPC Up
   PAX8-AS1, LINC01094, NEAT1 hsa-miR-29b-3p TET1, FOS, ANGPTL4, AQP4,
   FGG, VEGFA, TNFAIP3, MCL1 HPC Up
   BABAM2-AS1, PAX8-AS1, LINC01094, NEAT1 hsa-miR-20b-5p VEGFA, CDKN1A HPC
   Up
   LINC01094, NEAT1 hsa-miR-195-5p VEGFA, FGF2 HPC Up
   NEAT1, DLGAP1-AS1 hsa-miR-15a-5p VEGFA, IL10RA, HSPA1B HPC Up
   PAX8-AS1, LINC01094, NEAT1, MINCR hsa-miR-34a-5p VEGFA, NAMPT, FOS,
   CD44 HPC Up
   NEAT1 hsa-miR-126-5p VEGFA HPC Up
   NEAT1, LINC01094, DLGAP1-AS1 hsa-miR-15b-5p VEGFA HPC Up
   NEAT1 hsa-miR-126-3p ADM, SOX2, VEGFA, MERTK, NFKBIA BA46 Up
   NEAT1, LINC01094 hsa-miR-30b-5p BCL6, CAT BA46 Up
   NEAT1 hsa-miR-22-3p CDKN1A, BMPR1B BA46 Up
   LINC01094, NEAT1 hsa-miR-132-3p CDKN1A BA46 Up
   NEAT1, LINC01094 hsa-miR-17-5p CDKN1A, VEGFA BA46 Up
   PAX8-AS1, LINC01094, NEAT1 hsa-miR-20b-5p CDKN1A, VEGFA BA46 Up
   NEAT1 hsa-miR-138-5p FERMT2, VIM, YAP1 BA46 Up
   ZFHX4-AS1, LINC01094, NEAT1 hsa-miR-7-5p FOS BA46 Up
   NEAT1, MINCR, LINC01094, ZFHX4-AS1 hsa-miR-23a-3p HMGB2, MT2A BA46 Up
   NEAT1, LINC01094 hsa-miR-29a-3p MCL1, ZFP36, TET1, FGG, VEGFA, TNFAIP3,
   CD93 BA46 Up
   NEAT1, PAX8-AS1, MINCR, LINC01094 hsa-miR-98-5p NRAS BA46 Up
   MINCR, PAX8-AS1, LINC01094, ZFHX4-AS1, NEAT1 hsa-miR-34a-5p SOX2,
   VEGFA, NAMPT, CD44, FOS BA46 Up
   NEAT1, PAX8-AS1, LINC01094, MINCR hsa-let-7 g-5p TGFBR1 BA46 Up
   LINC01094, NEAT1 hsa-miR-15b-5p VEGFA BA46 Up
   LINC01094, PAX8-AS1, NEAT1 hsa-miR-29b-3p VEGFA, MCL1, ELAVL1, TNFAIP3,
   FOS, FGG, AQP4, ANGPTL4, TET1 BA46 Up
   NEAT1, LINC01094 hsa-miR-29c-3p VEGFA, MCL1, FGG BA46 Up
   NEAT1 hsa-miR-126-5p VEGFA BA46 Up
   ZFHX4-AS1, NEAT1 hsa-miR-9-5p VIM, BCL6, ELAVL1, TGFBI BA46 Up
   NEAT1, MINCR hsa-miR-124-3p VIM BA46 Up
   NEAT1, LINC01094 hsa-miR-195-5p YAP1, VEGFA, FGF2 BA46 Up
   NEAT1 hsa-miR-15a-5p YAP1, VEGFA, IL10RA, UCP2, HSPA1B BA46 Up
   DLX6-AS1 hsa-miR-34a-5p EPHA5 STR Down
   NEAT1 hsa-miR-126-3p ADM, VEGFA, NFKBIA, PIK3R2 STR Up
   LINC01094, NEAT1 hsa-miR-30b-5p BCL6 STR Up
   NEAT1 hsa-miR-22-3p CDKN1A STR Up
   NEAT1, LINC01094 hsa-miR-132-3p CDKN1A STR Up
   MINCR, LINC01094, NEAT1 hsa-miR-23a-3p FAS, HMGB2 STR Up
   NEAT1, LINC01094 hsa-miR-7-5p FOS STR Up
   LINC01094, MINCR, PAX8-AS1, NEAT1 hsa-miR-34a-5p FOS, VEGFA, NAMPT,
   CD44 STR Up
   LINC01094, NEAT1 hsa-miR-29a-3p MCL1, VEGFA, TNFAIP3, CD93, SPARC,
   ZFP36, FGG STR Up
   MINCR, PAX8-AS1, LINC01094, NEAT1 hsa-let-7 g-5p TGFBR1 STR Up
   DLGAP1-AS1, NEAT1 hsa-miR-15a-5p UCP2, VEGFA, IL10RA, HSPA1B, YAP1 STR
   Up
   NEAT1, LINC01094 hsa-miR-17-5p VEGFA, CDKN1A STR Up
   PAX8-AS1, NEAT1, LINC01094 hsa-miR-20b-5p VEGFA, CDKN1A STR Up
   DLGAP1-AS1, LINC01094, NEAT1 hsa-miR-15b-5p VEGFA STR Up
   LINC01094, PAX8-AS1, NEAT1 hsa-miR-29b-3p VEGFA, MCL1, ELAVL1, SPARC,
   ANGPTL4, AQP4, FGG, FOS, TNFAIP3 STR Up
   LINC01094, NEAT1, XIST hsa-miR-29c-3p VEGFA, MCL1, FGG, SPARC STR Up
   NEAT1 hsa-miR-126-5p VEGFA STR Up
   LINC01094, NEAT1 hsa-miR-195-5p VEGFA, YAP1, FGF2 STR Up
   NEAT1 hsa-miR-9-5p VIM, ELAVL1, BCL6, TGFBI STR Up
   NEAT1, MINCR hsa-miR-124-3p VIM STR Up
   NEAT1 hsa-miR-138-5p VIM, YAP1 STR Up
   MALAT1 hsa-miR-9-5p STARD13 LB Up
   [125]Open in a new tab

   Bold values are common ceRNA axes between all three areas of the brain.

   BA46 Brodmann area 46, DElncRNAs differentially expressed lncRNAs,
   DEmRNAs differentially expressed mRNAs, down downregulation, HPC
   hippocampus, LB lymphoblast, STR striatum, up upregulation.

Pathway enrichment analysis of DEmRNAs

   A KEGG enrichment analysis of all expressed DEmRNAs that were in the
   ceRNA network was conducted. The top five enriched KEGG pathways were
   “Kaposi sarcoma-associated herpesvirus infection”, “Hepatitis B”, “MAPK
   signaling pathway”, “Relaxin signaling pathway”, and “FoxO signaling
   pathway” (Fig. [126]4).

Figure 4.

   [127]Figure 4
   [128]Open in a new tab

   Overall results of pathway enrichment analysis. The bar chart shows the
   top five enriched pathways, along with their corresponding P-values.
   Colored bars correspond to terms with significant P-values (< 0.05). An
   asterisk (*) next to a P-value indicates the term also has a
   significant adjusted P-value (< 0.05). This chart was generated with
   Appyter ([129]https://appyters.maayanlab.cloud/#/)^[130]44.

PPI network construction and hub genes identification

   The PPI network of three brain regions is shown in Fig. [131]5a–c. We
   assessed their degree of connection and recognized 11 hub genes in
   three brain regions (BDNF, VEGFA, FGF2, FOS, CD44, SOX2, NRAS, SPARC,
   ZFP36, FGG, and ELAVL1), which were all upregulated except BDNF. The
   strongest connections with others were observed for the hub genes
   marked with red nodes. In addition, orange and yellow nodes represent
   the hub genes with moderate and weak connections, respectively
   (Fig. [132]5d–f). VEGFA had the strongest connections with others in
   the HPC region, FGF2 had moderate connections, while the hub genes
   CD44, BDNF, and FOS had weak connections. The VEGFA hub gene had the
   strongest connections with others in the BA46 region, whereas the hub
   genes CD44 and FGF2 had moderate connections, and the hub genes NRAS
   and SOX2 had weak connections. The VEGFA gene exhibited the strongest
   connections with others in the STR region, whereas the hub genes
   ELAVL1, SPARC, FGG, and ZFP36 had weak connections. Based on these
   findings, VEGFA had the strongest connections in all three regions of
   the brain. Two other genes (FGF2 and CD44) were also identified as hub
   genes in both HPC and BA46 regions.

Figure 5.

   [133]Figure 5
   [134]Open in a new tab

   Protein–protein interaction (PPI) network of differentially expressed
   mRNAs (DEmRNAs) and hub genes in three brain regions. Blue nodes show
   the interaction among DEmRNAs in the PPI network in (a) hippocampus
   (HPC), (b) Brodmann area 46 (BA46), and (c) striatum (STR) regions. Hub
   genes in (d) HPC, (e) BA46, and (f) STR regions identified from the PPI
   network. Red, orange, and yellow nodes represent the hub genes with
   strong, moderate, and weak connections, respectively. This figure was
   generated using Cytoscape version 3.8.0
   ([135]https://cytoscape.org/)^[136]39.

Reconstruction of the subceRNA network

   Based on the hub genes, the lncRNA-miRNA-hub genes’ subnetwork of three
   brain regions was reconstructed. In total, nine lncRNAs (DLX6-AS1,
   NEAT1, MINCR, LINC01094, DLGAP1-AS1, BABAM2-AS1, PAX8-AS1, ZFHX4-AS1,
   and XIST), 17 miRNAs (hsa-miR-124-3p, hsa-miR-132-3p, hsa-miR-15a-5p,
   hsa-miR-22-3p, hsa-miR-126-3p, hsa-miR-15b-5p, hsa-miR-17-5p,
   hsa-miR-195-5p, hsa-miR-20b-5p, hsa-miR-29a-3p, hsa-miR-29b-3p,
   hsa-miR-29c-3p, hsa-miR-34a-5p, hsa-miR-7-5p, hsa-miR-126-5p,
   hsa-miR-98-5p, and hsa-miR-9-5p) and 11 hub genes (BDNF, VEGFA, FGF2,
   FOS, CD44, SOX2, NRAS, SPARC, ZFP36, FGG, and ELAVL1) were included
   (Fig. [137]6). We used “key lncRNAs”, “key miRNAs”, and “key mRNAs”
   keywords to refer to these transcripts, respectively. CeRNA axes in
   which the lncRNA NEAT1 regulates the VEGFA via sponging different
   miRNAs (e.g., miR-29 family and miR-20b), as well as
   MALAT1/hsa-miR-9-5p/STARD13 (the only axis identified in the LB), could
   be good candidates for future experimental studies in different brain
   regions and peripheral blood of SCZ patients, respectively.

Figure 6.

   [138]Figure 6
   [139]Open in a new tab

   The lncRNA-miRNA-hub gene subceRNA network. ceRNA axes in (a)
   hippocampus, (b) Brodmann area 46, and (c) striatum regions. The red
   and blue nodes represent the upregulation and downregulation,
   respectively. Gray edges represent interactions between RNAs. LncRNAs,
   miRNAs, and mRNAs are represented by hexagon, round rectangle, and
   ellipse, respectively. This figure was generated using Cytoscape
   version 3.8.0 ([140]https://cytoscape.org/)^[141]39. CeRNA, competing
   endogenous RNA.

Discussion

   The expression patterns of ceRNA differ depending on tissue-related,
   cellular, and subcellular circumstances. A network might consist of
   different types of ceRNAs, such as lncRNAs, circRNAs, pseudogenes, as
   well as mRNAs^[142]45. LncRNAs, one of the major forms of RNAs, are
   identified in the ceRNA machinery and perform a significant role in
   physiological and pathological cellular mechanisms. It has been
   revealed that they are implicated in a variety of mental disorders
   (e.g., AD, SCZ, and autism spectrum disorder (ASD))^[143]46. It is now
   well accepted that lncRNA expression varies by tissue, cellular types,
   and developmental levels. This tissue specificity, besides subcellular
   distributions, clearly indicates that lncRNA expression is tightly
   regulated^[144]47. According to the as-mentioned theoretical concepts,
   the ceRNA regulation network related to lncRNAs can play a crucial role
   in SCZ pathogenesis. In this study, we used a public database to
   download the expression profiles of four different types of SCZ tissues
   to evaluate the DElncRNAs and DEmRNAs in SCZ and normal tissues and
   then create a lncRNA-miRNA-mRNA regulatory network. We identified
   possible lncRNA-miRNA-mRNA axes in the pathogenesis of SCZ consisting
   of nine key lncRNAs, 17 key miRNAs, 11 key mRNAs in the brain, and one
   key lncRNA, one key miRNA, one key mRNA in LB.

   LncRNA expression has been hypothesized to be dysregulated in SCZ and
   might be involved in the pathogenesis^[145]48. Among the key lncRNAs
   identified by our analysis, only the association of NEAT1, MALAT1, and
   DLX6-AS1 with SCZ has been the specific focus of previous studies.
   NEAT1 is a member of the lncRNA family highly enriched in the mammalian
   brain^[146]49. It is a crucial scaffolding unit as well as a structural
   determinant of paraspeckles^[147]50, which is a subnuclear structure
   implicated in several cellular processes such as splicing and
   transcriptional modulation through chromatin structure
   modification^[148]49. This type of lncRNA has already been found to be
   upregulated in other types of mental illnesses, namely Huntington’s
   disease^[149]51 and ASD^[150]52. Studies related to the comparison of
   NEAT1 expression in the patients with SCZ against the control group
   have shown contradictory findings. Some studies identified
   downregulation of NEAT1 in peripheral blood^[151]53 and multiple
   subareas of the cerebral cortex such as BA46 and HPC^[152]54 of SCZ
   patients compared with healthy controls. However, others revealed no
   significant difference in peripheral blood^[153]18, or upregulation in
   BA46, HPC, and STR regions^[154]7, which is similar to our results.
   Further studies are needed to explain these contradictory findings.
   MALAT1 is predominantly expressed in nerve cells and is
   transcriptionally elevated in nuclear speckles. It has been reported
   that MALAT1 employs splicing proteins of the Ser/Arg (SR) family in
   transcription sites in order to regulate synaptogenesis-related gene
   expression in hippocampal neurons cultivated. Furthermore, in this
   system, knocking down MALAT1 results in a reduction in synaptic
   density, while overexpression promotes the density^[155]55. A study
   recently has shown no substantial variation in the expression level of
   MALAT1 in SCZ patients’ peripheral blood compared with healthy
   control^[156]56. This finding is inconsistent with our results which
   might be attributed to the treatment of the patients with clozapine.
   The upregulation of DLX6-AS1, a regulator of GABAergic interneuron
   development^[157]57, in cerebral organoids derived from induced
   pluripotent stem cells showed that it might play a key role in ASD and
   SCZ pathogenesis^[158]58. In addition, we found XIST as another key
   lncRNA in the ceRNA network, which is associated with mental disorders.
   LncRNA XIST regulates the inactivation of the X chromosome, which is a
   critical developmental mechanism in the brain. XIST is highly expressed
   in lymphoblastoid cells from female patients diagnosed with bipolar
   disorder or recurrent severe depression as well as postmortem human
   brains of female psychiatric patients. Our results are in line with
   these findings. To our knowledge, the association of the other six key
   lncRNAs (MINCR, LINC01094, DLGAP1-AS1, BABAM2-AS1, PAX8-AS1, and
   ZFHX4-AS1) has not been the specific focus of investigations in the
   context of mental disorders. MINCR is a Myc-induced lncRNA related to
   some cancers via acting as a ceRNA^[159]59. Similarly, recent studies
   report that LINC01094^[160]60, DLGAP1-AS1^[161]61, and PAX8-AS1^[162]62
   can act as a ceRNA in some cancers. Besides, it has been shown that
   BABAM2-AS1 (also known as BRE-AS1)^[163]63 and ZFHX4-AS1^[164]64
   regulate cancer cells through multiple mechanisms. The association
   between some of the lncRNAs with SCZ has been reported in this study
   for the first time; therefore, further studies are needed to validate
   our findings.

   MiRNAs control target gene expression by attaching to the
   non-transcript site of the targeted gene, which then influences
   biological and signals transduction pathways in a cell, with potential
   effect on the onset and progression of SCZ^[165]65. In our work, we
   found that key lncRNAs may regulate hub genes via sponging 17
   SCZ-associated key miRNAs (hsa-miR-124-3p, hsa-miR-132-3p,
   hsa-miR-15a-5p, hsa-miR-22-3p, hsa-miR-126-3p, hsa-miR-15b-5p,
   hsa-miR-17-5p, hsa-miR-195-5p, hsa-miR-20b-5p, hsa-miR-29a-3p,
   hsa-miR-29b-3p, hsa-miR-29c-3p, hsa-miR-34a-5p, hsa-miR-7-5p,
   hsa-miR-126-5p, hsa-miR-98-5p, and hsa-miR-9-5p). Our prediction is
   that most of these miRNAs show decreased expression, except for miRNAs
   that target hub genes that have reduced expression (e.g., BDNF).
   According to the HMDD database, some of these miRNAs, such as
   miR-132^[166]66, miR-195^[167]67, miR-20b, miR-29 family^[168]68, are
   downregulated, while some others, such as miR-15 family,
   miR-126^[169]15, miR-17^[170]69 are upregulated in some brain regions
   of SCZ patients. Although most of these findings can confirm our
   results, the predicted ceRNA axes must still be validated by molecular
   techniques such as luciferase reporter systems, co-immunoprecipitation
   assays, and PCR.

   In the present study, KEGG enrichment analysis was performed as well.
   The results showed that most of the DEmRNAs involve in viral infection
   and three different signaling pathways. A previous study found a higher
   prevalence of human herpesvirus 8 (HHV8) infection in the patients with
   SCZ against the control group^[171]70. HHV8 has been associated with
   the development of Kaposi sarcoma^[172]71. Previously, a meta-analysis
   demonstrated that people with SCZ are more likely to have infectious
   hepatitis, which could be attributed to behavioral, immunological
   parameters, or both^[173]72. It is also shown that expression levels of
   mitogen-activated protein kinase (MAPK)^[174]73 and forkhead box O
   (FOXO) pathway-related genes^[175]74 are altered in SCZ. Association
   between relaxin polymorphisms and SCZ also suggests that the relaxin
   pathway plays a role in the development of SCZ^[176]75. Dysfunction of
   these signaling pathways might lead to functional pathology in SCZ.

   We found that 15 mRNAs (ADM, CLEC2B, F13A1, HLA-DQB1, HLA-DRB1,
   HLA-DRB3, HLA-DRB4, HLA-DRB5, CSF2RB, RHPN2, DDX3Y, EIF1AY, KDM5D,
   RPS4Y1, USP9Y) and one lncRNA (XIST) were differentially expressed in
   the LB and brain of patients with SCZ using a Venn diagram, suggesting
   that these RNAs could be regarded as a biomarker for SCZ. Among them,
   ADM^[177]76 and HLA-DRB1^[178]77 have been suggested as biomarkers in
   previous studies.

   Twelve key mRNAs (BDNF, VEGFA, FGF2, FOS, CD44, SOX2, NRAS, SPARC,
   ZFP36, FGG, ELAVL1, and STARD13) were reported in this study. BDNF is
   capable of regulating neuroplasticity, inhibiting the cascade of
   apoptosis, and raising the expression of many cellular proteins
   essential for neurogenesis, neuronal proliferation, as well as
   survival^[179]78. A myriad of studies have shown that low BDNF levels
   are associated with poor cognitive skills in SCZ patients^[180]79,
   which is consistent with our result. The VEGFA molecule is necessary
   for angiogenesis, brain plasticity, and neurogenesis, and it differs
   from other neurotrophic factors due to its significant angiogenic role.
   VEGFA can enhance the proliferation of neural stem cells, endothelial
   cells, mature astrocytes, and neurons. VEGF signaling is changed in
   SCZ, according to the results from transcriptomic, epigenomic, and
   genetic variation studies in hippocampus and prefrontal cortex
   samples^[181]80. In addition, FGF2 is involved in angiogenesis. FGF2
   was later confirmed to be a neurotrophic factor, significant in neural
   development and central nervous system maintenance. Serum levels of
   FGF2 protein were found to be substantially higher in patients with
   first-episode, drug-free SCZ^[182]81. As one of the immediate early
   genes most studied in the brain, FOS is believed to play a significant
   role in SCZ pathophysiology. In line with our findings, FOS is
   upregulated in some areas of an SCZ rat model brain^[183]82.
   Contrarily, some research found that FOS is highly expressed in
   non-neural peripheral specimens in addition to the reduction of the FOS
   in the brain of SCZ patients^[184]83. More research is required to
   understand these conflicting findings. Also, contradictory results have
   been obtained from the study of CD44^[185]84–[186]86 and SPARC^[187]87
   dysregulation components of the brain extracellular matrix, in some
   brain regions. Another key gene in our study was SOX2. The association
   of this gene, which is a promoter of neurogenesis, with SCZ has been
   investigated in some animal studies^[188]88,[189]89. The NRAS-encoded
   protein, a Ras superfamily intrinsic GTPase, is commonly associated
   with cancer progression and advancement. It is also essential in
   neurodevelopmental disorders due to its involvement in signal
   transduction from activated receptors toward the MAPK cascade. A
   previous study reported higher expression of this gene in the frontal
   cortex of SCZ patients^[190]90. This data is in agreement with our
   results. The genetic network of hemostatic procedure demonstrated that
   FGG, as a coagulation-related mediator gene, may have a role in
   SCZ^[191]91. ELAVL1 plays a crucial role in mRNA stabilization
   (mediated by AU-rich elements) and translation as well as neuronal mRNA
   post-transcriptional regulation. The capability of ELAVL1 in
   stabilizing mRNA targets may be influenced by posttranslational
   modifications. It is postulated that an inefficient regulation of MeCP2
   transcription in the SZ cohort is caused by alterations in one of its
   interactant targets, ELAVL1^[192]92. To the best of our knowledge, the
   role of two other key genes (i.e., ZFP36 and STARD13) has not been
   investigated in SCZ thus far. Similar to ELAVL1, ZFP36 is
   disease-relevant RNA-binding proteins (RBPs) that interact with AU-rich
   sequences, but unlike ELAVL1, ZFP36 mediates the decay of target
   mRNAs^[193]93. In addition, its association with mental conditions
   (e.g., AD^[194]94 and suicidal behavior^[195]95) has been reported
   previously. Moreover, STARD13, a GTPase-activating protein (GAP),
   inhibits RhoA and Cdc42 specifically in glioblastoma cells and other
   types of tumor models^[196]96.

   There are some limitations in the present study. First, variations in
   methodology and the preparation of samples, the platform, patient
   characteristics, and data processing may all have an effect on gene
   expression profiles^[197]97–[198]99. Secondly, it is challenging to
   explain dynamic variations in disease progression and the emergence of
   complications using postmortem brains^[199]13. Finally, a validation
   analysis is needed to substantiate the reliability of our findings.

Conclusion

   In conclusion, this study identified a lncRNA-associated ceRNA network
   potentially relevant for SCZ. This network includes key mRNAs
   (including BDNF, VEGFA, FGF2, FOS, CD44, SOX2, NRAS, SPARC, ZFP36, FGG,
   ELAVL1, and STARD13), and lncRNAs (DLX6-AS1, NEAT1, MINCR, LINC01094,
   DLGAP1-AS1, BABAM2-AS1, PAX8-AS1, ZFHX4-AS1, and XIST). While the
   functions of this network need to be investigated, the present study
   provides potential research targets to study molecular mechanisms which
   could be relevant for the pathogenesis of SCZ.

Supplementary Information

   [200]Supplementary Figure 1.^ (1.4MB, docx)
   [201]Supplementary Table 1.^ (14.6KB, docx)
   [202]Supplementary Table 2.^ (19.3KB, txt)
   [203]Supplementary Table 3.^ (20.5KB, txt)
   [204]Supplementary Table 4.^ (19.4KB, txt)
   [205]Supplementary Table 5.^ (12.1KB, txt)

Author contributions

   M.T., H.S. and M.R. wrote the manuscript and revised it. M.M.M.,
   B.M.H., Ma.M.M., N.K.A. and M.R.A. performed the bioinformatic analysis
   and collected the information. The authors contributed equally and are
   fully aware of submission.

Competing interests

   The authors declare no competing interests.

Footnotes

   Publisher's note

   Springer Nature remains neutral with regard to jurisdictional claims in
   published maps and institutional affiliations.

Contributor Information

   Mohammad Taheri, Email: Mohammad_823@yahoo.com.

   Maryam Rezazadeh, Email: Rezazadehm@tbzmed.ac.ir.

Supplementary Information

   The online version contains supplementary material available at
   10.1038/s41598-021-03993-3.

References