Deciphering the landscape of lncRNA-driven ceRNA network in schizophrenia etiology

Background The unifying hypothesis of competing endogenous RNA (ceRNA) wherein crosstalk between coding (mRNAs) and long non-coding RNAs (lncRNAs) via microRNA (miRNA) response elements, creates a pervasive regulatory network across the transcriptome, has been implicated in complex disorders including schizophrenia. Even with a wide range of high-throughput data, the etiology of schizophrenia remains elusive, necessitating a more holistic understanding of the altered genetic landscape, shifting focus from solely candidate gene studies and protein-coding variants. Objective We developed lncRNA-associated ceRNA networks to elucidate global molecular/regulatory signatures underlying schizophrenia using diverse data in the public domain. Methods Microarray dataset associated with peripheral blood mononuclear cells (PBMCs) of schizophrenia and con-trol patients was used to identify differentially expressed mRNAs. Weighted gene co-expression network analysis (WGCNA) was used to identify highly correlated hubs, and genes from these overlapping Kyoto Encyclopedia of Genes and Genomes (KEGG) and gene ontology (GO) term genesets were considered key mRNA players. StarBase, Human MicroRNA Disease Database, and miRWalk were used to derive mRNA-miRNA and miRNA-lncRNA relationships. Finally, the key mRNAs, interacting lncRNAs and miRNAs were chosen to reconstruct sub-ceRNA networks based on network centrality scores. Results Bioinformatics analysis revealed the involvement of three differentially expressed mRNAs, namely ADRA1A , HAP1 and HOMER3 in the schizophrenia ceRNA networks with lncRNAs NEAT1 , XIST , and KCNQ1OT1 modulating their activity by a suggestive sequestering of miR-3163, miR-214-3p and miR-2467-3p, respectively. Conclusions Furthermore, based on contextual evidence, we propose how ceRNAs could orchestrate crosstalk between neurostructural dynamics and immune/inflammatory processes and enable unifying these disparate models of schizophrenia etiology.


Introduction
Schizophrenia presents itself as a severe psychiatric disorder in early to late adolescence with a complex manifestation of positive symptoms (e.g., hallucinations and delusions), negative symptoms (e.g., social withdrawal, blunted affect) and pervasive cognitive deficits associated with functional decline [1,2].Though not as common as other psychiatric disorders with a lifetime prevalence of 0.32% [3] affected worldwide, it executes a high toll on society via multifaceted indices such as economic liabilities, human rights violations and suicides [4,5].Its pathophysiology is perplexing as uncovered through advanced genomic [6], developmental neurobiology [7], and systems biology investigations [8,9] with irregularities in cellular, molecular, neuroanatomical, and neurophysiological domains [10] being implicated.This diverse milieu of factors and their enigmatic interplay contributes to its elusive etiology.
As far as the genetics of schizophrenia are concerned the strong links between schizophrenia and over 100 sus- ceptibility loci, along with identified CNVs and SNVs, show promise.Thousands of common alleles with small effects collectively contribute substantially to schizophrenia risk.These findings may lead to new therapeutic insights.However, we must remember: (i) associations between genetic variants and schizophrenia do not necessarily imply causal pathways; (ii) many associations extend beyond schizophrenia to other mental disorders.The specifics of schizophrenia's origins and genotypeenvironment interactions are largely unknown, warranting caution in assessing the various genetic contributors to its development [11].The evolution of pharmacological interventions to alleviate the patient's condition has been extremely slow since the advent of the first antipsychotics [12].Efforts toward biomarker discovery for early identification of individuals at risk, improving diagnostic accuracy and precision, predicting treatment response, and to obtain new druggable targets are notable [13].However, the prevalent body of work regarding this has mainly focused on the coding part of the genome but the mosaic manifestation of schizophrenia warrants a far more holistic understanding of the underlying pathophysiology.
Toward this end, integrating existing knowledge of schizophrenia etiology with the recent shift in our understanding of the ncRNAs as junk/black matter of the genome to being vital regulatory molecules is of prime importance [14][15][16].By virtue of their ability to silence/alter expressions of multiple targets simultaneously, ncRNAs can affect entire signaling pathways [17].Not surprisingly, understanding the role of ncRNAs in physiological and pathological conditions has significantly enriched our recognition and understanding of alternate/additional molecular pathways involved.More than a hundred miRNAs and though in its infancy ∼ 30 lncRNAs have been found dysregulated across the prefrontal cortex, superior temporal gyrus, parietal cortex, amygdala, serum, and peripheral blood in schizophrenia [18].SNPs in both miRNAs and lncRNAs have also been shown to have significant associations with the schizophrenia phenotype [19].Collectively, substantial evidence has been generated as to the dysregulation of ncRNAs in schizophrenia and that the disruption of the networks they regulate is critical to schizophrenia as they are enriched with target genes pertaining to various neurophysiological processes [20,21].Considering all these, a new paradigm wherein the clinical utility of ncRNAs in the form of a next-generation is now being seriously assessed [22].
Nevertheless, making sense of the huge amount of data on ncRNAs in the public domain for targeted therapeutics is a herculean task [23][24][25][26].Integrating the raw sequencing data from multiple sources and fitting them under logical models is one way of understanding how lncRNAs collectively interact with the coding genome.A possible route is the ceRNA theory, wherein ncRNAs such as lncRNAs that share common MREs with mRNAs can act as siphons and competitively sequester miRNAs, thus forming a complex regulatory network [27].Any differential expression in these RNAs harboring common MREs could thereby lead to imbalances in the regulatory network and disease development [28][29][30].ceRNA networks developed for diseases such as sarcopenia have identified lncRNAs, mRNAs and miRNAs which add to disease risk with very high accuracy [31].Several lines of evidence have already shown a close association of ceRNA networks with several forms of cancer but very little is known for schizophrenia [32].As substantial evidence of an altered immune landscape in schizophrenia has accumulated [33][34][35] and recent research has solidified the probable role of lncRNAs in orchestrating an altered immune landscape in the disease [36], we narrowed down the analysis to focus on identifying immune/ inflammatory signals by using a relevant expression dataset from publicly available database with the propensity to identify these peripheral signals, if any and performed bioinformatics analyses to highlight the lncRNA-associated ceRNA networks.It is hypothesized that such strategies may potentially reveal putative novel pathways and/or give credence to erstwhile proposed theories for disease etiology.These may identify/uncover novel candidate biomarkers and aid in understanding molecular pathways underlying schizophrenia.
The novelty of our study lies in usage of HTS mRNA expression profile followed by hub module identification via WGCNA, PPIN, enrichment, and ceRNA network analyses.This study is the first of its kind where schizophrenia biomarkers are extrapolated via a robust and multi-stage protocol.These predictive biomarkers may be beneficial for patients' early diagnosis and treatment as they play a vital role in schizophrenia pathogenesis.Further validation is required via wet laboratory experimentations in order to prove the efficacy and accuracy of these biomarkers, which comes in as a limitation of the present study.Also, the availability of fewer datasets poses another limitation and a comprehensive meta-analysis study would be much more beneficial in near future.

Selection of schizophrenia-associated mRNA expression profile and DEA
We accessed NCBI-GEO [37] (https:// www.ncbi.nlm.nih.gov/ geo/) to retrieve schizophrenia-associated mRNA expression profiles utilizing "schizophrenia" as a suitable keyword.All results were further shortlisted as per the following inclusion criteria: (i) the dataset(s) must be "ncRNA expression profiling by array" type with all its samples belonging to "Homo sapiens"; (ii) the dataset(s) must include raw as well as processed data; (iii) the dataset(s) must be submitted to GEO from last 10 years till present (i.e., 2012 to 2024); (iv) the dataset(s) must contain at least 25 samples; (v) dataset(s) must com- prise schizophrenia and healthy control patient samples; (vi) patient samples within dataset(s) must be acquired from PBMCs.Studies devoid of any review articles, abstracts, case reports, non-human samples and cellline-based experimental study designs were excluded.The series matrix file of the selected dataset was downloaded followed by further QCs (i.e., normalization and log 2 transformation) .The ARSyNseq function available within the NOISeq package [38,39] was utilized (with unknown batch settings) to acquire batch-corrected expression values in R. Probe IDs were mapped to their corresponding HGNC symbols via ArrayExpress online portal [40] (https:// www.ebi.ac.uk/ biost udies/ array expre ss) with respect to the sequencing platform of the selected dataset.An average was taken for the expression value of genes mapping to more than one probe ID to avoid redundancy.An unpaired two-sample t test was utilized for computing p-values and log 2 (fold change) values of all genes between schizophrenia and normal patients via limma package in R [41].All p-values were also corrected via BH method in R. The genes were considered as differentially expressed corresponding to a BH − p − value < 0.01 and log 2 (fold change) > 0.5 [42].DEGs having log 2 fold change > 0.5 and log 2 fold change < −0.5 were designated as upregu- lated and downregulated, respectively.PCA method was utilized to assess the sample aggregation degree.It is an unsupervised method that can be used to understand the difference between two or more sample groups [43,44].Unsupervised PCA/dimensionality reduction was performed via R software based on the DEGs expression with respect to samples [45,46].

WGCN construction and hub module selection
WGCNA implemented in R [45,47] assisted in constructing a WGCN from schizophrenia-associated DEGs followed by the detection of representative modules.Primarily, goodSamplesGenes function was used to check the data for any missing values and LV genes.Next, the pickSoftThreshold function assisted in choosing β based on SFT.The power adjacency function assisted in transforming similarity matrix into a weighted adjacency matrix.Thereafter, we computed TOM and dissTOM followed by clustering dendrogram (hierarchical) construction via hclust function.DTC algorithm was applied to identify gene modules from branches of the tree [48].ME and MEdiss were computed followed by joining module(s) with similar high-expression profiles based on ME dendrogram.The module with significantly highest correlation between MM and k.in was finalized, and the genes with MM > 0.9 were regarded as hub genes for further analysis.

Schizophrenia-associated 3-node ceRNA network construction and topological analysis
The miRNAs interacting with the functionally enriched schizophrenia-associated hub mRNAs were acquired from miRWalk 2.0 [58] (http:// mirwa lk.umm.uni-heide lberg.de/) and ENCORI databases [59] (https:// rnasy su.com/ encori/), respectively.Combined miRWalk and ENCORI make up for each others drawbacks thereby leading to a robust prediction scheme.ENCORI is the second most cited database for such predictions, and miRWalk harbors experimental data from luciferase assays, microarrays, NGS, pSILAC, and western blot experiments.Furthermore, where miRWalk made up for ENCORI's last update of 2021 by being updated twice a year and containing predictions from MirTarBase, ENCORI was crucial in understanding the lncRNA interactions which miRWalk lacked [60].miRNAs binding on 3′UTR region, having an interaction score > 0.95 , and binding gap = 1 were considered as significant and retrieved from miRWalk.The list of miRNAs were validated from available literature and those associated with schizophrenia were finally retained.lncRNAs interacting with schizophrenia-associated miRNAs were cataloged using ENCORI.All data were extracted using a low stringency of >= 1 in both CLIP data and degradome data.Finally, with multiple lines of evidence indicating that TFs may also be involved in ceRNA crosstalk, TFs interacting with the hub mRNAs as their targets were acquired using TF2DNA [61] (https:// www.fiser lab.org/ tf2dna_ db/) at a significant p − value < 0.0001.

Schizophrenia-associated mRNA expression profile selection and DEA
Based on the abovementioned inclusion and exclusion criteria, a schizophrenia-associated mRNA expression profile possessing accession number GSE54913 was selected.It comprised 12 healthy controls and 18 schizo- phrenia PBMC samples.The QC analysis revealed that the expression values were already normalized and log 2 -transformed in the preprocessed series matrix file.Post batch correction, all the probe IDs were mapped to their corresponding HGNC symbols (protein-coding) with respect to Arraystar Human LncRNA microarray V2.0 platform.Finally, the expression values corresponding to duplicate gene symbols were averaged leading to 12971 unique genes.Corresponding to a BH − p − value < 0.01 and log 2 (fold change) > 0.5 , 2169 DEGs were finally screened.The distribution of upregulated ( n = 950 ) and downregulated ( n = 1219 ) schizophrenia-associated DEGs along with nonsignificant genes ( n = 10802 ) were visualized using a volcano plot (Fig. 1A).The log 2 (fold change) value of total genes varied with respect to the −log 10 ( BH − p − value ).

WGCN construction and hub module selection
The expression data of 2169 schizophrenia-associated DEGs and sample information were given as an input to WGCNA.WGCN was established at β = 12 (cor- responding to R 2 = 0.8 ) with no sample outliers and LV genes.Figure S1A-D shows plots for β in consideration with SFT criteria.The clustering dendrogram (hierarchical) and DTC algorithm gave fourteen color-coded modules (i.e., black, blue, brown, cyan, green, greenyellow, magenta, pink, purple, red, salmon, tan, turquoise, yellow) ranging in sizes from 34 to 522 (Fig. 2A).The modules with highly co-expressed gene patterns were merged together by cutting the dendrogram at a height of 0.2 (which corresponds to a correlation of 0.8 ) as shown in Fig. S2.After merging, fourteen modules were clubbed into four color-coded modules (i.e., black, brown, salmon, magenta) ranging in sizes from 36 to 1152 .The clustering dendrogram (hierarchical) based on dissTOM and ME with original ( n = 14 ) and merged ( n = 4 ) color modules is shown in Fig. 2B.WGCN is shown as a heatmap plot depicting TOM among merged color modules in Fig. S3.Based on the most significant correlation between MM and k.in (Table S1), the brown module ( MM vs k.in = 0.85, p − value = 6.0 × 10 −155 ) was picked as the hub module.Scatterplot of MM versus k.in across brown module is shown in Fig. 2C.Heatmap of brown module DEGs along with their corresponding ME levels is shown in Fig. 2D.Finally, a total of 119 DEGs were identified from the brown hub module corresponding to MM > 0.9.

PPIN, pathway, and GO term enrichment analyses
A total of 118 out of 119 DEGs from hub module suc- cessfully mapped to their corresponding protein names via STRING database.Corresponding to an interaction score > 0.4 , the PPIN comprised 35 nodes and 27 edges as shown in Fig. S4. Figure S5 shows cen- trality distributions like betweenness, closeness, ND, TC, NC, and ASPL of PPIN.We entered all 119 DEGs from hub module into Enrichr.A total of 11 , 16 , 37 , 12 KEGG pathways, GO-CC, GO-MF, GO-BP terms were obtained (Tables S2-S5) corresponding to p − value < 0.05 .Venn plot as shown in Fig. 3A illustrates three overlapping hub genes (i.e., ADRA1A, HAP1, HOMER3) between significant pathway, GO-BP, GO-MF, and GO-CC genesets.The box-and-whisker plots show the relative expression distribution of ADRA1A, HAP1, and HOMER3 across control and schizophrenia patient samples (Fig. 3B-D).As observed, mRNA expression levels of all three hub genes were significant across control and schizophrenia samples.

Schizophrenia-associated 3-node ceRNA network construction and topological analysis
The schizophrenia-associated 3-node ceRNA network comprised 767 nodes and 3030 edges as shown in Fig. 4. The breakup of nodes and edge pairs is summarized in Table S6.Tables S7−S8 shows top 3 lncRNAs and miR- NAs within ceRNA network ranked based on betweenness, degree, and closeness centralities.Within this network, degree of lncRNAs, miRNAs and mRNAs ranged from 1 to 117 , 2 to 50 , and 24 to 90,respectively.Average degrees of lncRNAs, miRNAs, and mRNAs were 4.64 , 21.18 , 47.66 , respectively.As observed from these centralities, ADRA1A hub gene was repressed and regulated by maximum miRNAs while lncRNA NEAT1 interacted with the highest number of miRNAs.Numerous  14) and merged ( 4 ) module colors.The merged module sizes were as follows: brown ( 551 ), black ( 430 ), salmon ( 36 ), magenta ( 1152 ).C Scatterplot showing significant ( p − value < 0.05 ) correlation of MM with k.in across brown module genes.D Expression heatmap of brown module genes where the rows and columns correspond to genes and samples.The red and green color bands in the heatmap signify higher and lower expression level across brown module genes.Also, the corresponding ME expression levels (along y-axis) across all samples (along x-axis) are displayed at the bottom panel of module heatmap in the form of barplot miRNAs, lncRNAs, and mRNAs participating in higherorder subnetwork motifs were observed and the top three higher-order subnetwork motifs based on highest centrality scores of betweenness, degree, and closeness have been reported.The first-ranked subnetwork motif comprised one miRNA (miR-3163), one lncRNA (NEAT1), and one hub gene (ADRA1A).The second-ranked subnetwork motif comprised one miRNA (miR-214-3p), one lncRNA (XIST), and one hub gene (HOMER3).And, the third-ranked subnetwork motif comprised one miRNA (miR-2467-3p), one lncRNA (KCNQ1OT1), and one hub gene (ADRA1A) (Fig. 5A-C).Figure S6 shows centrality distributions like betweenness, closeness, ND, TC, NC, and ASPL of 3-node ceRNA network.

Discussion
Rare/ultra-rare protein-coding variants de novo or captured through WES of familial forms of schizophrenia have provided insights into a few genes from dopaminergic and neurodevelopmental pathways in schizophrenia [69][70][71][72] but heritability and etiology remain unexplained.Regulatory variants such as ncRNAs are emerging to be essential players in our understanding of the biology/ etiology of common conditions such as schizophrenia [73,74].As one of the most common types of ncRNAs, lncRNAs are believed to play a pivotal role in the ceRNA machinery and elicit a significant effect in both physiological and pathological mechanisms.Multiple lines of evidence have implicated them in various psychiatric disorders [75][76][77].Their differential expression in tissue, cell types, and developmental levels indicates that lncRNA expression is tightly regulated [78][79][80].These give further credence to the idea that lncRNA-associated ceRNA networks may play a crucial role in schizophrenia etiology.However, the dire lack of studies on lncRNAs in the public domain has made it difficult for bioinformatic analyses to annotate their role in disease biology sufficiently.This is exacerbated further by the dearth of HTSeq studies in schizophrenia, including lncRNAs and the lack of postmortem data.To account for these, we employed an approach to shortlist coding genes and their interacting miRNAs and then extract the lncRNA-miRNA interactions reported in the public domain, thus making the mRNA-lncRNA-disease interaction hypothesis-free.
Furthermore, recent studies have indicated the link between brain and periphery via the circulatory system, which contains secreted regulatory molecules and hormones produced in the diffused NES that impact the peripheral markers' gene expression pattern [81][82][83][84].These findings confirm that schizophrenia is a systemic disorder and support the notion that biomarkers in peripheral samples such as WB, PBMCs, lymphoblasts and olfactory epithelium may be insightful.Another line of evidence that dictated the choice of blood expression profiles for the analysis was based on the current evidence wherein immune/inflammatory processes are located in the disorder [33][34][35] and the strong connections established between altered immunity and lncR-NAs [36,85].Understanding the networks at play in the peripheral system might help generate a holistic view of the underlying connection.
In the enrichment analysis using the highest-order WGCNA module, three genes were found overlapping in all the significant pathway and GO term libraries tested.
ADRA1A was enriched for terms such as neuroactive ligand-receptor interaction, cytosolic Ca 2+ ion transport, Fig. 4 Schizophrenia-associated 3-node ceRNA network comprising 767 nodes and 3030 edges.Magenta-colored diamond nodes represent the lncRNAs, red circular nodes represents the miRNAs, and green-colored octagonal nodes represents the hub genes Fig. 5 A Top higher-order subnetwork motif based on betweenness, degree, and closeness comprising one miRNA (miR-3163), one lncRNA (NEAT1), and one hub gene (ADRA1A).B The second higher-order subnetwork motif comprising one miRNA (miR-214-3p), one lncRNA (XIST), and one hub gene (HOMER3).C Third higher-order subnetwork motif comprising one miRNA (miR-2467-3p), one lncRNA (KCNQ1OT1), and one hub gene (ADRA1A).Magenta-colored diamond nodes represent the lncRNAs, red circular nodes represents the miRNAs, and green-colored octagonal nodes represents the hub genes positive regulation of GABAergic synaptic transmission.HOMER3 was enriched for glutamatergic synapse and G protein-coupled glutamate receptor binding and was a cellular component of dendrites.HAP1 was associated with pathways involved in neurodegeneration, neurogenesis, neurotrophin binding and similar to HOMER3 was found located in the dendrites.All of the pathways have either been directly reported in schizophrenia etiology before or are of substantial importance in the processes involved in schizophrenia pathophysiology [86][87][88][89][90][91].Thereby, any perturbation in their expression could potentially disturb these pathways and initiate or maintain the disease phenotype.With this hypothesis, we assessed the lncRNAs and miRNAs that could putatively dysregulate these key mRNAs leading to disease manifestation.
Among the three lncRNAs, NEAT1 has recently been reported to be upregulated [ log 2 (fold change) > 2 ] in Brodmann area 46, hippocampus and striatum.NEAT1 is highly enriched in the mammalian brain and is an indispensable structural component of paraspeckles which are membrane-less cellular bodies involved in several cellular processes such as splicing and transcriptional modulation through chromatin structure modifications with emerging evidence suggesting their altered abundance with several innate immune activating responses stimuli such as sequestering to IFNGR1 [92].NEAT1 itself has been cited as lncRNA-type immunoregulator (i) affecting monocyte-macrophage functions and T cell differentiation [93], (ii) assembly of inflammasomes by recruitment, maturation, and stabilization of CASP1 in activated macrophages [94], (iii) elicits pro-proliferative and anti-apoptotic roles and migration, invasion, and inflammatory cytokines secretion [95], (iv) exhibits innate immunity responses against viral infections [96].Furthermore, multiple studies have shown miRNA sequestering tendencies of NEAT1, thereby attenuating target gene activity [97][98][99][100][101].Based on our bioinformatic analyses, we propose that NEAT1 could be sequestering miR-3163 because of its higher number of transcripts in the disease state, thereby elevating the repression of ADRA1A which was downregulated in cases in the DEA (Fig. 3B).
The second lncRNA, XIST has also been previously associated with multiple mental disorders [102][103][104].XIST is involved in the inactivation of the X chromosome, which has a long-standing reputation for harboring genes important for brain development and function [105].Outside of its silencing roles, XIST i) stimulates proliferation and differentiation of naive CD4 + T cells [106], (ii) is delocalized in B cells of female-biased autoimmunity [107], (iii) in-part promotes CD11c + atypi- cal B cell formation [107], and (iv) has been shown to perturb PDL1 levels by probable competitive binding of miR-34a-5p [108].Even though the expression profile of HOMER3 in schizophrenia is unknown, HOMER1 (member of the three-member HOMER family) has been shown to be up and downregulated in schizophrenia depending on the tissue type with variants in both found to be associated with schizophrenia [109].Overexpression of XIST has been reported in bipolar disorder and major depressive disorder (phenotypes closely associated with schizophrenia) as well but is highly tissue-specific [102].Therefore, in a similar fashion as NEAT1, XIST could be competitively sequestering to miR-214-3p, a miRNA already known to target the Qki [110], thereby leading to the altered HOMER3 levels.
The third lncRNA, KCNQ1OT1, targeting KCNQ1, though actively involved in epigenetic phenomenon via chromatin modifications, HMT G9a, and PRC2 [111], has no direct association with schizophrenia yet.However, it does seem to sponge miR-15a, leading to immune evasion and malignant progression of prostate cancer via upregulating PDL1, an essential immune checkpoint [112].This might explain its expression correlation with CD4 + , CD8 + , and cytotoxic T cell levels and several other immune cell subsets in another ceRNA reported in colorectal cancer [85].Furthermore, it might be indirectly associated with increased sudden cardiac arrest in schizophrenia patients [113].The KCNQ1 protein forms functional potassium channels [114].Multiple lines of evidence, structural variants and mice knockouts, have shown KCNQ1 to be associated with LQT1, a condition synonymous with increased adverse cardiac events [115][116][117].It is established that all atypical antipsychotics affect the cardiac potassium pump and that about ∼ 6 − 10% of schizophrenia patients show a longer QT interval under treatment [118,119].We propose that the expression of KCNQ1 as dictated by altered KCNQ1OT1 levels could be a putative cause of these adverse drug reactions.
Further investigations into this aspect could potentially lead to pre-emptive treatment strategies.Though the levels of KCNQ1OT1 in schizophrenia are not known, we can extrapolate from available knowledge that rs8234 [120] leads to lower expression of KCNQ1 and is also associated with reduced processing speed, reduced white matter FA and higher risk for schizophrenia [121], thereby implying that lower levels of KCNQ1 are associated with these impaired phenotypes.Elevated KCNQ1OT1 levels could also propagate this scenario.Considering these derived associations, studies establishing KCNQ1OT1 levels in schizophrenia could be informative.
We could thereby imply that elevated KCNQ1OT1 transcripts could be competitively sequestering to miR-2467-3p and inhibiting the expression of ADRA1A, thereby leading to its downregulated state in the DEA.
Interestingly, ADRA1A is also associated with several cardiac conditions [122][123][124].This gives a glimpse into the intricate mechanism in which ceRNAs could be acting in the pathophysiology of schizophrenia.
In conclusion, this study identified lncRNAs NEAT1, XIST and KCNQ1OT1-associated ceRNA networks which could be potentially relevant to schizophrenia by interacting with schizophrenia-relevant genes, ADRA1A and HOMER3.Furthermore, the affinity of the mRNAs to neurodevelopmental processes and that of the lncR-NAs to immune/inflammatory processes might indicate a mechanism to unite the two most significant models proposed in schizophrenia etiology.Of note, though the current analyses is based on data specific to schizophrenia, neuroinflammation and its effect on neurodynamics is a well-established phenomenon in a variety of psychiatric illnesses such as depression, bipolar depression and obsessive-compulsive disorder [125].ceRNAs established through this study and new ones discovered by using similar methods have the potential of uncovering further such pathways.Further refinements in such prediction strategies have the potential of unveiling additional interactions in schizophrenia biology, which, eventually, systems biology approaches coupled with artificial intelligence and machine learning technologies can integrate into a holistic picture.However, it is also important to note that the prediction strategy deployed in this study does not take into account the miRNA and potential ceRNA expression levels.This is important as it is well-established that both miRNAs and ceRNAs have temporal, spatial, and disease-specific expression patterns.Furthermore, studies have shown that ceRNAs and miRNAs with concentrations within a particular range are capable of eliciting such crosstalks.Even though we have provided evidence to give strong credence to the highlighted ceRNA axes, these must still be validated by qRT-PCR, luciferase reporter systems and co-IP assays.Furthermore, we have discussed in favor of the standalone components of the ceRNA networks.We believe that additional investigations into their roles in the diseased state would be valuable in assessing their role as important biomarkers for schizophrenia.Further wet lab experimentations would be an asset in proving the efficacy and accuracy of the predicted biomarkers.Also, design of lead compounds as potential drugs post successful clinical trials could be helpful for the treatment of schizophrenia in near future.

Fig. 1 A
Fig. 1 A Volcano plot showing the distribution of 2169 schizophrenia-associated significant DEGs (upregulated: 950 + downregulated: 1219 ) and nonsignificant ( 10802 ) genes across GSE54913 dataset.B Annotation heatmap showing the expression distribution of top 10 down and top 10 upregulated schizophrenia-associated DEGs.Cluster dendrograms representing Euclidean distance-based hierarchical clustering for both columns and rows are presented along the top and left sides of the plot.Sample type annotation bar is presented at the top of heatmap.C PCA plot showing the expression variability of 2169 schizophrenia-associated DEGs across GSE54913 dataset.The relative expression level of all DEGs dimensionally reduced in compliance with sample type leading to distinct cluster formations were signified by solid circular points in the plot.Green and red-colored points within ellipses signify control and schizophrenia samples.The % of total variation accounted for by the 1st ( 58.8% ) and 2nd ( 7.9% ) PCs are shown on the x-and y-axes, respectively.D Scree plot displaying the % of explained variances captured by their corresponding PCs

Fig. 3 A
Fig. 3 A Overlapping hub genes between significant GO-BP, GO-MF, GO-CC, and pathway genesets.The red, blue, yellow, and green-colored areas signify KEGG, GO-BP, GO-MF, GO-CC genesets, respectively.Box-and-whisker plots showing expression intensity distribution of B HAP1, C HOMER3, D ADRA1A across control and schizophrenia patient samples.The top and bottom of the boxes signify 75th and 25th percentile of distribution.Horizontal lines within the boxes represent the median values while minimum and maximum values label the axes endpoints.P-values shown at the top of boxplots represent significance levels between sample groups for each hub gene