Identification and characterization of noncoding RNAs-associated competing endogenous RNA networks in major depressive disorder

Abstract

BACKGROUND

Major depressive disorder (MDD) is a common and serious mental illness. Many novel genes in MDD have been characterized by high-throughput methods such as microarrays or sequencing. Recently, noncoding RNAs (ncRNAs) were suggested to be involved in the complicated environmental-genetic regulatory network of MDD occurrence; however, the interplay among RNA species, including protein-coding RNAs and ncRNAs, in MDD remains unclear.

AIM

To investigate the RNA expression datasets downloaded from a public database and construct a network based on differentially expressed long noncoding RNA (lncRNAs), microRNAs (miRNAs), and mRNAs between MDD and controls.

METHODS

Gene expression data were searched in NCBI Gene Expression Omnibus using the search term “major depressive disorder.” Six array datasets from humans were related to the search term: GSE19738, GSE32280, GSE38206, GSE52790, GSE76826, and GSE81152. These datasets were processed for initial assessment and subjected to quality control and differential expression analysis. Differentially expressed lncRNAs, miRNAs, and mRNAs were determined, Gene Ontology and Kyoto Encyclopedia of Genes and Genomes enrichment analyses were performed, and protein-protein interaction network was generated. The results were analyzed for their association with MDD.

RESULTS

After analysis, 3 miRNAs, 12 lncRNAs, and 33 mRNAs were identified in the competing endogenous RNA network. Two of these miRNAs were earlier shown to be involved in psychiatric disorders, and differentially expressed mRNAs were found to be highly enriched in pathways related to neurogenesis and neuroplasticity as per Gene Ontology and Kyoto Encyclopedia of Genes and Genomes enrichment analyses. The expression of hub gene fatty acid 2-hydroxylase was enriched, and the encoded protein was found to be involved in myelin formation, indicating that neurological development and signal transduction are involved in MDD pathogenesis.

CONCLUSION

The present study presents candidate ncRNAs involved in the neurogenesis and neuroplasticity pathways related to MDD.

Key Words: Major depressive disorder; Noncoding RNA; Competing endogenous RNA; Bioinformatics; Data mining

Core Tip: Competing endogenous RNAs (ceRNAs) are novel regulatory molecules involved in a wide range of biological processes and diseases. This study explored the potential ceRNA networks (ceRNETs) involved in the pathogenesis of major depressive disorder (MDD) using bioinformatics data mining. A ceRNET comprising 3 miRNAs, 12 lncRNAs, and 33 mRNAs was constructed based on two public datasets obtained from the Gene Expression Omnibus database. Elucidating the correlation of the ceRNET with MDD opens new avenues to discover specific diagnostic biomarkers for MDD and expands our knowledge about this disease.

INTRODUCTION

Major depressive disorder (MDD), characterized by persistent and intense feelings of sadness for extended periods, is one of the most common causes of morbidity and mortality worldwide[1], yearly affecting approximately 5% of adults worldwide[2]. The exact cause of MDD remains unknown; however, the occurrence of MDD is widely believed to involve crosstalk between genetic, environmental, social, and developmental vulnerabilities and resilience factors[2,3]. Numerous risk factors for MDD have been identified from three major perspectives: medical, social, and substance[4,5]. Available treatment options primarily include pharmacotherapy, psychotherapy, and lifestyle changes. However, prevention and treatment of this disease remain difficult because of varying presentations, unpre-dictable course and prognosis, and variable responses to treatment[6]. Therefore, further studies are needed to determine the biological mechanisms underlying MDD to develop better patient therapies.

The increasing popularity of high-throughput microarray technologies has facilitated the identification of genome variations in MDD, improving the understanding of its pathogenesis and deve-lopment, and revealing promising biomarkers. Gao et al[7] identified 241 differentially expressed genes (DEGs) in the hippocampus (hip), 218 DEGs in the prefrontal cortex (pfc), and 327 DEGs in the striatum (str) from patients with MDD. These DEGs were enriched in glycan biosynthesis, RIG-I-like receptor signaling, and pyrimidine metabolism pathways, which significantly contribute to MDD pathogenesis. Additionally, the DEGs AR, PTK2, IRAK1, IL12A, GALNT12, GALNT2, CD19, and PTDSS2 were identified as novel therapeutic targets[7,8]. Moreover, Segman et al[9] found 73 DEGs between patients with postpartum major depression and controls using peripheral blood cells; the immune response, transcriptional effects on cell proliferation, and DNA replication and repair were significantly enriched in these patients. Although novel genes have been discovered and well-studied, genetic predisposition to MDD only partly explains the occurrence of this illness. With the development of epigenetic theory and technology, non-coding RNAs (ncRNAs) have been identified to help explain the interaction between genetics and the environment. Li et al[10] found an inverse relationship between serum brain-derived neurotrophic factor (BDNF) levels and miR-132/miR-182 levels in depression, supporting the notion that miR-182 is a putative regulatory microRNA (miRNA) of BDNF. Zhou et al[11] found that miRNAs act as key regulators of synaptic plasticity in MDD pathogenesis. Liu et al[12] detected four long non-coding ribonucleic acids (lncRNAs) related to the expression of mRNAs associated with MDD. Although ncRNAs have been proposed to play important roles in MDD, the regulatory network among RNA species remains unclear.

The competing endogenous RNA (ceRNA) network proposed by Salmena et al[13] has been verified in an increasing number of studies. This complicated post-transcriptional regulatory network plays a critical role in the progression and pathogenesis of various illnesses, including cancer. However, to date, no ceRNA network has been proposed for MDD. Here, we investigated RNA expression datasets downloaded from the Gene Expression Omnibus (GEO) database and constructed a network based on differentially expressed lncRNAs, miRNAs, and mRNAs between MDD and controls. This network was used to mine for functional ncRNAs that may contribute to the epigenetic mechanism underlying MDD pathogenesis and may be useful as potential therapeutic targets and biomarkers for MDD.

MATERIALS AND METHODS

Data resources and preprocessing

Gene expression data were searched in NCBI GEO using the search term “major depressive disorder”. Six array datasets from humans were related to the search term: GSE19738 (expression profiling by array, n = 42; MDD, 21; control, 21), GSE32280 (expression profiling by array, n = 24; MDD, 8; subsyndromal symptomatic depression, SSD, 8; control, 8), GSE38206 (expression profiling by array, n = 18; MDD, 9; control, 9), GSE52790 (expression profiling by array and ncRNA profiling by array, n = 22; MDD, 10; control, 12), GSE76826 (expression profiling by array, n = 22; MDD, 10; control, 12), and GSE81152 (ncRNA profiling by array, n = 50; MDD, 31; control, 19). These datasets were processed for initial assessment; if they met the requirements, they were subjected to quality control and differential expression analysis.

Identification of differentially expressed lncRNAs, miRNAs, and mRNAs

The differentially expressed mRNAs, lncRNAs, or miRNAs were determined based on the threshold criteria of |log2 (fold-change)| ≥ 1 and an adjusted P value < 0.05 using the R limma package (version 3.40.6). The differentially expressed lncRNAs (DE_lncRNAs) obtained from GSE76826 were entered into the miRcode database (http://www.mircode.org/) that provides an integrated, searchable map of putative target sites of miRNAs, to predict the interactions of DE_lncRNAs and DE_miRNAs. Only DE_lncRNAs for which miRNAs were successfully predicted were selected for ceRNA network construction. The miRNAs obtained were overlapped with differentially expressed miRNAs obtained from GSE81152, and the overlapping miRNAs were identified as a differentially expressed miRNA (DE_miRNA). The name of the DE_miRNA was entered into the miRDB (http://mirdb.org/) and TargetScan (http://www.targetscan.org/vert_72/) databases, both of which are online databases used to predict miRNA targets. The mRNAs obtained above were overlapped with differentially expressed RNAs identified from GSE76826. The DE_lncRNAs, DE_miRNAs, and DE_mRNAs obtained in the manner described above were used for ceRNA network construction.

DE_lncRNA/DE_miRNA/DE_mRNA network construction

Cytoscape 3.4.0 was used to visualize the DE_lncRNA/DE_miRNA/DE_mRNA network construction.

Gene Ontology function and Kyoto Encyclopedia of Genes and Genome pathway enrichment analysis

To help identify the potential biological functions of these genes in the MDD ceRNA network, Gene Ontology (GO) functional and Kyoto Encyclopedia of Genes and Genome (KEGG) pathway enrichment analyses were performed using the clusterProfiler package in R. GO terms have three categories; biological process, cellular compartment, and molecular function. P ≤ 0.05 indicated a significant difference in GO terms and KEGG pathways.

Construction of the protein interaction network

The protein interaction (PPI) network was constructed to evaluate the interactions among genes. The DE_mRNAs were selected to construct a PPI network based on Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) ver.11.0 (https://string-db.org/). In the PPI network, a node represents a gene; the undirected link between two nodes is an edge, representing the interaction between two genes, and the degree of a node corresponds to the number of interactions of a gene with other genes in the network.

RESULTS

Quality control and differential expression analysis

The workflow of the analysis is shown in Figure 1. Because the GSE19738 dataset lacks documents to annotate gene symbols and the GSE32280 dataset uses peripheral blood mononuclear cells rather than whole blood, both these datasets were excluded from analysis at the first stage. Next, sample-level quality control was performed using principal component analysis to ensure that the cases and controls were well-separated and to identify sample outliers. After removing the outliers, there were 10 controls and 8 MDD cases in GSE52790, 12 controls and 9 MDD cases in GSE76826, and 10 controls and 23 MDD cases in GSE81152. The four remaining datasets, GSE32280, GSE52790, GSE76826, and GSE81152, were utilized to identify differentially expressed mRNAs, lncRNAs, and miRNAs using the limma R package with the threshold set to an adjusted P value < 0.05 and |log2 (fold-change)| ≥ 1. A hierarchical cluster heatmap showing the expression patterns of the DEGs between MDD and controls is shown in Figure 2. As no DEGs were found in GSE32280, this dataset was excluded from subsequent analysis. A total of 125, 373, and 54 DEGs were identified in the GSE52790, GSE76826, and GSE81152 datasets, respectively. There were 112 upregulated and 13 downregulated genes among the 125 DEGs in GSE52790, 109 upregulated genes (including 89 mRNAs and 10 lncRNAs) and 264 downregulated genes (including 198 mRNAs and 66 LncRNAs) in the 373 DEGs of GSE76826, and 51 upregulated and 3 downregulated genes in the 54 DEGs of GSE81152.

Figure 1 Flowchart of the study. A: Step 1, dataset screening; B: Step 2, competing endogenous RNA network construction. DEGs: Differentially Expressed Genes; mRNA: Messenger RNA; miRNA: MicroRNA; lncRNA: Long non-coding RNA; ceRNA: Competing endogenous RNAs; miRcode: http://www.mircode.org/; miRDB: http://mirdb.org/; TargetScan: http://www.targetscan.org/vert_72/; DE_mRNA: Differentially expressed mRNA; DE_lncRNA: Differentially expressed lncRNA; DE_miRNA: Differentially expressed miRNA.

Figure 2 Difference and cluster analyses. A: Volcanic maps of differentially expressed genes in GSE52790; B: Cluster analysis heatmap of the expression patterns of the differentially expressed genes in GSE52790; C: Volcanic maps of differentially expressed genes in GSE76826; D: Cluster analysis heatmap of the expression patterns of the differentially expressed genes in GSE76826; E: Volcanic maps of differentially expressed genes in GSE81152; F: Cluster analysis heatmap of the expression patterns of the differentially expressed genes in GSE81152. Red indicates upregulated genes; green indicates downregulated genes.

Construction of a ceRNA network for MDD

To better understand the effect of lncRNAs on mRNAs mediated by combinations with miRNAs in MDD, we constructed a ceRNA network based on the GSE76826 and GSE81152 data. Eighty-six differentially expressed lncRNAs from GSE76826 were retrieved and entered into the miRcode database to predict the potential interacting miRNAs. Eighty-four miRNAs (645 LncRNA-miRNA pairs) were identified in the miRcode database. Next, we overlapped these 84 miRNAs to the differentially expressed miRNAs identified in GSE81152, generating 3 shared miRNAs (the corresponding lncRNAs were reduced to 12) (Table 1). Thereafter, we entered these 3 miRNAs into the miRDB and TargetScan databases to identify their potential target mRNAs and combined the results from these two databases to achieve high true-positive coverage. A total of 1,868 mRNAs was generated from the two databases; when the 287 differentially expressed mRNAs identified from GSE76826 were overlapped, 33 shared mRNAs were identified (Table 2). Finally, we assembled all identified ceRNA pairs to generate an mRNA-miRNA-lncRNA network in MDD (Figure 3).

Figure 3 Major depressive disorder-associated enriched competing endogenous RNA network. Red and green indicate upregulated downregulated genes, respectively. The circles represent differentially expressed RNA (DE_RNA), triangles represent DE long noncoding RNA (DE_lncRNA), and rounded squares represent DE_microRNA.

Table 1 The twenty differentially expressed long noncoding RNAs and differentially expressed_microRNA pairs identified in this study.

lncRNA	miRNA
STAG3L4	hsa-miR-4465
MTERFD2	hsa-miR-4465
TCL6	hsa-miR-4465
MDS2	hsa-miR-4465
LINC00202	hsa-miR-4465
LINC00402	hsa-miR-4465
TCL6	hsa-miR-137
ANKRD36BP2	hsa-miR-137
LINC00202	hsa-miR-137
LINC00402	hsa-miR-137
MIR600HG	hsa-miR-137
MTERFD2	hsa-miR-125b-5p
HCG4	hsa-miR-125b-5p
TCL6	hsa-miR-125b-5p
UBE2Q2P2	hsa-miR-125b-5p
TREML2P1	hsa-miR-125b-5p
KIAA0125	hsa-miR-125b-5p
ANKRD36BP2	hsa-miR-125b-5p
LINC00202	hsa-miR-125b-5p
MIR600HG	hsa-miR-125b-5p

lncRNA: Long noncoding RNA; miRNA: microRNA.

Table 2 The forty-three differentially expressed_mRNAs and differentially expressed microRNA pairs identified in this study.

miRNA	mRNA
hsa-miR-4465	FAM98B
hsa-miR-4465	ADM
hsa-miR-4465	INHBB
hsa-miR-4465	FA2H
hsa-miR-4465	HOXA9
hsa-miR-4465	GRB10
hsa-miR-4465	UGT8
hsa-miR-4465	ACVR1C
hsa-miR-4465	COL19A1
hsa-miR-4465	TTC28
hsa-miR-4465	HAS3
hsa-miR-4465	SLC22A23
hsa-miR-4465	ADAM23
hsa-miR-4465	NACC2
hsa-miR-4465	SH3PXD2A
hsa-miR-4465	CELSR1
hsa-miR-4465	PLEKHG1
hsa-miR-4465	ZSWIM5
hsa-miR-137	AKAP2
hsa-miR-137	HLF
hsa-miR-137	FGF9
hsa-miR-137	GLIS2
hsa-miR-137	SIPA1L2
hsa-miR-137	WNT7A
hsa-miR-137	NRG1
hsa-miR-137	ADAM23
hsa-miR-137	LRRC4
hsa-miR-137	SLC22A23
hsa-miR-137	COL19A1
hsa-miR-137	TTC28
hsa-miR-137	APLN
hsa-miR-137	NACC2
hsa-miR-137	BACH2
hsa-miR-137	ACVR1C
hsa-miR-125b-5p	NIPAL4
hsa-miR-125b-5p	NR6A1
hsa-miR-125b-5p	ZSWIM5
hsa-miR-125b-5p	GALNT14
hsa-miR-125b-5p	GRB10
hsa-miR-125b-5p	ELOVL4
hsa-miR-125b-5p	ACVR1C
hsa-miR-125b-5p	RALGPS2
hsa-miR-125b-5p	GLIS2

FA2H: Fatty acid 2-hydroxylase.

Gene set enrichment analysis

Enrichment analysis indicated that the 33 differentially expressed mRNAs were significantly associated with the GO terms for 887 biological process, 105 molecular function, and 55 cell component terms. For the KEGG pathways, the only significant pathway was the mammalian target of rapamycin (mTOR) signaling pathway. Table 3 and Figure 4 list the top 10 most significant biological processes, molecular functions, and cellular components from the GO terms.

Figure 4  Gene Ontology functional enrichment analysis.

Table 3 The significant Gene Ontology terms identified in this study.

Category	Term description	Count of genes	P value
Biology process	Activin receptor signaling pathway	3	5.54E-05
Biology process	Peripheral nervous system development	3	0.000292
Biology process	G protein-coupled receptor internalization	2	0.0003
Biology process	Desensitization of G protein-coupled receptor signaling pathway	2	0.000487
Biology process	Negative adaptation of signaling pathway	2	0.000487
Biology process	Positive regulation of vascular endothelial growth factor receptor signaling pathway	2	0.000487
Biology process	Adaptation of signaling pathway	2	0.000597
Biology process	Uterus development	2	0.000656
Biology process	Sphingolipid biosynthetic process	3	0.000679
Biology process	Reproductive structure development	5	0.00078
Cellular component	Glutamatergic synapse	4	0.002395
Cellular component	Schaffer collateral-CA1 synapse	2	0.00786
Cellular component	Cul2-RING ubiquitin ligase complex	1	0.024149
Cellular component	Integral component of synaptic membrane	2	0.025933
Cellular component	Intrinsic component of synaptic membrane	2	0.029786
Cellular component	Extracellular matrix	3	0.044750
Cellular component	Podosome	1	0.047732
Cellular component	Excitatory synapse	1	0.075303
Cellular component	Integral component of postsynaptic density membrane	1	0.079823
Cellular component	Intrinsic component of postsynaptic density membrane	1	0.084321
Molecular function	Receptor ligand activity	6	0.000191
Molecular function	Hormone activity	3	0.001352
Molecular function	UDP-glycosyltransferase activity	3	0.002314
Molecular function	Guanyl-nucleotide exchange factor activity	4	0.002772
Molecular function	Growth factor activity	3	0.003326
Molecular function	Transferase activity, transferring hexosyl groups	3	0.006074
Molecular function	Phosphatidylinositol-4,5-bisphosphate 3-kinase activity	2	0.006746
Molecular function	Cytokine activity	3	0.007189
Molecular function	Phosphatidylinositol bisphosphate kinase activity	2	0.007334
Molecular function	Transcription factor activity, RNA polymerase II proximal promoter sequence-specific DNA binding	4	0.008356

PPI network construction

A total of 33 DEG mRNAs were entered into the STRING database to build a PPI network. Gene nodes with the highest W values in the network were ACVR1C, ELOVL4, fatty acid 2-hydroxylase (FA2H), INHBB, NIPAL4, and UGT8 (Figure 5).

Figure 5 Protein interaction network construction. FA2H: Fatty acid 2-hydroxylase.

DISCUSSION

Considerable progress has been made in understanding the pathophysiology of MDD; however, no single model or mechanism can fully explain all aspects of the disease[6]. MDD is believed to involve the following: Reduced levels of monoamines, changes in the hypothalamic-pituitary-adrenal axis, inflammation, limited neurogenesis, changes in brain structure and function, heredity factors, environmental milieu, and epigenetics[6]. Epigenetics, a field focusing on gene-environment interactions, has gained attention over the past decades. The recent application of microarray-based genome-wide expression analysis has enabled the identification of MDD-associated genes, including novel ncRNAs. Although numerous studies have demonstrated that ncRNAs play important roles in MDD pathogenesis, the mechanisms of how these genes regulate MDD have not been well-characterized[14,15]. Thus, we performed ceRNA network analysis by data mining to determine the potential regulatory mechanisms of MDD.

Multiple GO terms were identified in enrichment analysis, reflecting the complexity of the disorder. We observed that protein-coding genes in the network were highly correlated with cellular growth and developmental differentiation, including the activin receptor signaling pathway, peripheral nervous system development, G protein-coupled receptor internalization, receptor-ligand activity, hormone activity, growth factor activity, and cellular growth. The term “signaling by activins” in our results is consistent with previous studies showing that activin can modulate depression and anxiety-related behavior[16-21]. Activins are members of the transforming growth factor-β superfamily that are expressed in various tissues, including neuronal cells, and are involved in proliferation, differentiation, metabolism, homeostasis, apoptosis, immune response, and tissue repair. Postnatal neurogenesis, which leads to the production of new neurons in the adult brain, can influence the replacement of damaged neurons, stress responses, memory formation, and depression, among others[22]. Previous studies have reported that the activin activity in the adult forebrain influences locomotor activity, anxiety-related behavior, and hippocampal neurogenesis, which is also associated with age-related cognitive decline[23]. Additionally, consistent with our results, utilizing a genome-wide association study, Smeeth et al[23] suggested a complex relationship among reproductive hormones, hippocampal neurogenesis, and depression. Based on transcriptome array results analysis of Chinese patients with MDD, Liu et al[12] reported the enrichment of differentially expressed coding genes in the signal transduction pathway and basic metabolic process associated with neurodevelopmental disease[12,13]. Our KEGG results also showed that the mTOR signaling pathway is a part of the regulatory network of MDD. This agrees with previous studies which showed that the mTOR signaling pathway, which senses and integrates diverse extracellular stimuli to promote cellular growth or limit catabolic processes, contributes to normal neuronal growth and is associated with neurogenesis[24]. Notably, the mTOR signaling pathway has become an important target for treating depression using drugs such as ketamine[25,26].

The use of ncRNAs as biomarkers of psychological disorders is gaining momentum. In the ceRNA network developed in this study, there were 12 lncRNAs, 3 miRNAs, and 33 mRNAs specific to MDD. Consistent with our results, hsa-mir-137, as a brain-enriched miRNA, was earlier confirmed as a gene related to MDD, bipolar disorder, schizophrenia, and Parkinson’s disease susceptibility[27-30], whose underlying mechanisms regulate synaptic plasticity. Moreover, Zhao et al[31] detected significantly lower hsa-miR-137 Levels in the brain and peripheral blood in post-stroke depression rats, and Kim et al[32] found that it can be used as a diagnostic marker for methamphetamine withdrawal syndrome. Similarly, hsa-miR-125b-5p in the peripheral blood is highly expressed in Alzheimer’s disease compared with that in controls.

In addition to genetic components, different epigenetic mechanisms play an important role in the occurrence and development of MDD. A previous study revealed that the non-protein-coding RNA repressor of NFAT, human accelerated region 1, transcribed antisense of rheelin, and B-secretase-1 are associated with cognitive function, potentially contributing to MDD, and the lncRNA BDNF-AS was reportedly related to synaptic plasticity and potentially associated with MDD[33-36]. In a recent study, Abedpoor and colleagues reported that the lncRNA network could play an indispensable role in regulating depression-like behaviors in mice[37]. Furthermore, two studies reported lower RMRP expression in a mouse model of depression and MDD relative to normal subject samples[14]. In another study, Issler et al[38] observed that LINC00473 was downregulated in the prefrontal cortex of depressed females. Collectively, these results suggest that lncRNAs are potential diagnostic and prognostic biomarkers. Our results revealed 12 different lncRNAs involved in the MDD ceRNA network, suggesting the existence of a regulatory network in MDD.

In the present study, FA2H was highlighted in the PPI network. FA2H encoded by FA2H is highly expressed in the human brain and is involved in the formation of 2-hydroxy galactosylceramides and 2-hydroxy sulfatides in myelin[39,40]. The FA2H gene appears to be related to lipid metabolism and is one of the 10 candidate genes identified in an inherited neurologic disorder known as neurodegeneration with brain iron accumulation and is also involved in hereditary spastic paraplegia SPG35 and leukodystrophy[41,42]. Mutations in FA2H have been identified in autism spectrum disorders[43]. The function of FA2H does not appear to be limited to the brain since elevated FA2H expression has also been found in hepatocellular carcinoma and lung adenocarcinoma involving UGT8[44,45]. However, there have been no reports of FA2H in MDD and thus, this is the first study to suggest that FA2H is related to MDD. Considering its role in myelin formation, this gene may contribute to MDD patho-genesis via impacting neurological development and modulating signal transduction.

Elongation of very long-chain fatty acids-4 (ELOVL4) was another protein of interest identified in the PPI. This protein is essential for the synthesis of very long-chain polyunsaturated (VLC-PUFA) and saturated fatty acids (VLC-SFA) of chain lengths greater than 26 carbons[46]. The VLC-PUFAs play an important role in maintaining neural tissue homeostasis. Studies have suggested that dysregulated ELOVL4 expression may be involved in the lipid alterations observed in neuroblastoma, and mutations in the ELOVL4 gene could cause several distinct neurodegenerative diseases, including stargardt-like macular dystrophy, spinocerebellar ataxia 34, and a neuro-ichthyotic syndrome with severe seizures and spasticity[47,48]. However, like FA2H, none of the studies so far have linked ELOVL4 and MDD, and therefore, future validation benchmark experiments are needed.

There were some limitations to our study. The potentially important genes associated with MDD lack gene expression and functional validation, and hence, further in vivo and in vitro studies are required. Furthermore, the two datasets were not from the same platform, which may have resulted in a few deviations. However, the present study provides novel insight into MDD pathogenesis.

CONCLUSION

To date, MDD diagnosis mainly depends on the patient’s subjective expression, which may be misinterpreted by clinicians. Therefore, identification of a reliable biomarker that can be used to diagnose MDD and guide its treatment is imperative. In conclusion, the present study revealed that hsa-miR-4465, hsa-miR-137, and hsa-miR-125b-5p may be useful genetic biomarkers for MDD as they are potentially involved in the neurogenesis and neuroplasticity pathways in MDD. These results require further validation in future studies.

ARTICLE HIGHLIGHTS

Research background

Major depressive disorder (MDD) is a common and serious mental illness. Many novel genes in MDD have been characterized by high-throughput methods such as microarrays or sequencing. Recently, noncoding RNAs (ncRNAs) were suggested to be involved in the complicated environmental-genetic regulatory network of MDD occurrence.

Research motivation

The interplay among RNA species, including protein-coding RNAs and ncRNAs, in MDD remains unclear.

Research objectives

To investigate the RNA expression datasets and construct a network based on differentially expressed long noncoding RNA (lncRNAs), miRNAs, and mRNAs between MDD and controls through data mining method.

Research methods

Gene expression data were searched and downloaded from NCBI Gene Expression Omnibus database. Six array datasets from humans were related to the search term: GSE19738, GSE32280, GSE38206, GSE52790, GSE76826, and GSE81152. These datasets were processed for initial assessment and subjected to quality control and differential expression analysis. Differentially expressed lncRNAs, miRNAs, and mRNAs were determined, Gene Ontology and Kyoto Encyclopedia of Genes and Genomes enrichment analyses were performed, and protein-protein interaction network was generated. The results were analyzed for their association with MDD.

Research results

After analysis, 3 miRNAs, 12 lncRNAs, and 33 mRNAs were identified in the ceRNA network. Two of these miRNAs were earlier shown to be involved in psychiatric disorders, and differentially expressed mRNAs were found to be highly enriched in pathways related to neurogenesis and neuroplasticity as per Gene Ontology and Kyoto Encyclopedia of Genes and Genomes enrichment analyses. The expression of hub gene fatty acid 2-hydroxylase was enriched, and the encoded protein was found to be involved in myelin formation, indicating that neurological development and signal transduction are involved in MDD pathogenesis.

Research conclusions

The present study presents candidate ncRNAs involved in the neurogenesis and neuroplasticity pathways related to MDD.

Research perspectives

Bioinformatic data mining method can be an cost-effective way to explore potential biomarkers for complicated diseases. However, benchmark experiment is needed in the future to validate the results.