A Whole Transcriptome Analysis in Peripheral Blood Suggests That Energy Metabolism and Inflammation Are Involved in Major Depressive Disorder

Introduction: Previous studies on transcriptional profiles suggested dysregulation of multiple RNA species in major depressive disorder (MDD). However, the interaction between different types of RNA was neglected. Therefore, integration of different RNA species in transcriptome analysis would be helpful for interpreting the functional readout of the transcriptome in MDD.

Methods: A whole transcriptome sequencing were performed on the peripheral blood of 15 patients with MDD and 15 matched healthy controls (HCs). The differential expression of miRNAs, lncRNAs, circRNAs, and mRNAs was examined between MDD and HCs using empirical analysis of digital gene expression data in R (edgeR). Weighted correlation network analysis (WGCNA) was used to identify RNA co-expression modules associated with MDD. A ceRNA network was constructed for interpretation of interactions between different RNA species. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were conducted to explore potential biological mechanisms associated with MDD.

Results: Multiple RNAs and co-expression modules were identified to be significantly dysregulated in MDD compared to HCs. Based on the differential RNAs, a ceRNA network that were dysregulated in MDD were constructed. The pathway networks that related to oxidative phosphorylation and the chemokine signaling were found to be associated with MDD.

Conclusion: Our results suggested that the processes of energy metabolism and inflammation may be involved in the pathophysiology of MDD.

Introduction

Major depressive disorder (MDD) is a complex psychiatric disorder with high morbidity, disability rate, and great clinically heterogeneous features involving disturbances of mood and cognitive function. According to the World Health Organization (1), MDD was predicted to be the disease with the largest burden of disease in the world by 2030 (1). However, there is still a lack of biological indicators for the diagnosis and treatment of MDD. The diagnosis of MDD depends on self-reports from patients and clinical observations (2–4). Moreover, only 35–45% of patients achieved full clinical remission in the first treatment even if they received standard antidepressant treatment (5–8). Therefore, exploring the underlying mechanisms of MDD and identifying effective biomarkers are urgently required to facilitate diagnosis and treatment, which would also enhance our understanding of MDD.

Previous research reported that the prevalence of MDD demonstrated familial aggregation, which was due to the interaction between genetic effects and substantial individual-specific environmental effects, such as stressful life events, poverty and unstable intimacy (9, 10). Transcriptomics is a discipline that studies the transcription and regulation of all genes in a given cell population. The regulation of the transcriptome can effectively explain the dysregulated expression of genes in different physiological or pathological states of different individuals and is helpful for the exploration of biological markers related to physical status or diseases (11, 12). The advances in bioinformatics methods in recent years have greatly promoted the interpretation of the biological meaning of transcriptome profiles (11, 13, 14).

A number of previous studies have focused on transcriptional profiles in MDD and have shown that non-coding RNAs (ncRNAs) have potential value as biomarkers for the diagnosis and therapy of MDD (15). Accumulated evidence indicates that multiple types of RNA, including messengers (mRNAs) (16–18), microRNAs (miRNAs) (19–21), long non-coding RNAs (lncRNAs) (22, 23) and circular RNAs (circRNAs) (22, 23), are dysregulated in brain tissues, fibroblasts and peripheral blood in MDD. A correlation between RNA expression profiles and depressive symptoms was also reported (24, 25). These results suggested that alterations in transcriptome profiles might participate in the pathophysiology of MDD and be amenable for the identification of biomarkers of MDD.

The human transcriptome contains a large number of protein-coding genes as well as non-coding transcripts. Converge evidence suggests that the interaction among various types of RNA molecules plays critical roles in gene expression regulation (26). Recently, the competing endogenous RNA (ceRNA) hypothesis was proposed to interpret interactions among multiple types of RNA (27). Multiple RNA species, including mRNAs, lncRNAs and circRNAs, contain a large number of miRNA binding sites, i.e., miRNA response elements (MREs). In the ceRNA hypothesis, RNA transcripts containing MREs regulate each other by competitively binding miRNAs through shared MREs (28). Therefore, bioinformatics prediction of ceRNA networks in whole transcriptome analysis became helpful and meaningful for better interpreting the functional readout of the transcriptome (29).

In this study, a whole transcriptome sequencing was performed on the peripheral blood of patients with MDD and healthy controls (HCs) to explore the underlying molecular mechanism of MDD. The differentially expressed RNAs, dysregulation of co-expression networks, ceRNA networks and biological pathways between MDD and HCs were identified.

Materials and Methods

Participants

Patients were recruited from West China Hospital of Sichuan University from August 2017 to October 2018. All participants were right-handed Han Chinese and aged from 18 through 60 years. The Wechsler Adult Intelligence Scale Revised in China (WAIS-RC) was used to investigate the intelligence quotient (IQ) of participants (30, 31). To exclude the influence of medication, all patients in this study were required to be antidepressant-naive (never taken any antidepressant medication) or antidepressant-free for at least 3 months before enrollment. The depressive and anxiety symptoms of all patients were assessed using the 17-item Hamilton Depression Scale (HAMD-17) and 14-item Hamilton Rating Scale for Anxiety (HAMA-14), respectively (32, 33). The exclusion criteria for patients were (1) having serious physical diseases, such as neurogenic diseases, endocrine diseases, metabolic disorders, autoimmune diseases, cardiovascular and cerebrovascular diseases, etc.; (2) having any other psychiatric disorders in Axis I or II, such as schizophrenia, bipolar disorder, alcohol or drug abuse and history of loss of consciousness; (3) IQ ≤ 90; and (4) unsuitable or inability to complete the study for other reasons.

All HCs were recruited through advertising and were also interviewed by professional psychiatrists using SCID/NP. The exclusion criteria were as follows: (1) current and past history of depression or any psychiatric disorders in Axis I or II; (2) any psychiatric illness in their first-degree relatives; (3) serious physical diseases; (4) pregnancy or breast feeding; (5) reproductive age without adequate contraception; (6) IQ less than 90; and (4) unsuitable or unable to complete the study for other reasons.

The flowchart of this study is shown in Figure 1.

FIGURE 1

Figure 1. Study flow diagram. Flowchart of the comprehensive analysis process and all methods utilized in present study.

RNA Extraction and Quality Control

Peripheral blood (3 ml) samples from all participants were collected using Tempus Blood RNA Tubes (Applied Biosystems, Foster City, CA, United States), mixed vigorously for 10 s immediately and stored at –80°C until use. Total RNA was isolated using the MagMAX for Stabilized Blood Tubes RNA Isolation Kit (Thermo Fisher Scientific, Waltham, MA, United States) according to the manufacturer’s instructions. The concentration of RNA samples was determined using a NanoDrop ND-2000. The integrity of the RNA samples was assessed using a Bioanalyzer 2100.

RNA Sequencing

One microgram of RNA was used for RNA sequencing. The mRNA, lncRNA, and circRNA in total RNA samples were enriched by removing rRNA using the Ribo-Zero™ rRNA Removal Kit (Illumina). The RNA samples were then fragmented into 200–500 nt fragments using fragmentation buffer. A reverse transcription reaction was performed to generate the first strand cDNA complementary to the fragmented RNA, and then the second strand cDNA was synthesized with incorporation of dUTP using DNA polymerase I. The double stranded cDNA was purified with a QiaQuick PCR extraction kit (Qiagen), followed by end repairing and A-tailing, and Illumina sequencing adapters were ligated to the cDNA fragments. The second strand of cDNA marked with dUTP was digested by uracil-DNA-glycosylase (UDG). The remaining strand was amplified by PCR to generate a cDNA library suitable for sequencing. RNA sequencing was performed using the NovaSeq 6000 Sequencing System (Illumina).

Identification of mRNAs, lncRNAs, and circRNAs

Reads obtained from the sequencing machines were filtered to generate clean reads using Bowtie2 (RRID:SCR_016368). The mapped reads were assembled into transcripts using Cufflinks software (RRID:SCR_014597).

Small RNA Sequencing

Two micrograms of RNA was used for small RNA sequencing. RNA samples were separated by polyacrylamide gel electrophoresis (PAGE), and RNA molecules in a size range of 18 ∼ 30 nt were extracted. Then, 3′ adapters and 5′ adapters were ligated to RNAs. The ligation products were reverse transcribed and amplified through PCR using primers complemented to 3′ and 5′ adapters. The amplification products with sizes of 140–160 bp were extracted to generate a cDNA library and sequenced using Illumina HiSeq™ 2500 (Illumina).

Identification of miRNAs

Reads obtained from the sequencing machines were filtered to generate clean reads according to the following rules: (1) removing reads without 3′ adapters; (2) removing reads containing 5′ adapters; (3) removing low-quality reads containing more than one low-quality (Q value ≤ 30) base or containing unknown nucleotides; (4) removing reads with no insertion between 5′ and 3′ adapters; (5) removing reads containing poly A; and (6) removing reads shorter than 18 nt (not including adapters). The clean reads were aligned with small RNAs in the GenBank database (Release 209.0) and Rfam database (11.0) to identify and remove reads that contained sequences of rRNAs, scRNAs, snoRNAs, snRNAs, and tRNAs. Then, the reads were aligned with the reference genome to remove reads containing fragments of mRNAs. Finally, the reads were aligned with the miRBase database (Release 21) to identify miRNAs.

The lncRNA sequencing, small RNA sequencing and identification of lncRNAs, mRNAs, circRNAs and miRNAs were conducted by Gene Denove Biotechnology Co (Guangzhou, China).

Differential Expression Analysis

The differential expression of RNAs in the comparison between MDD and HCs was conducted using the R package EdgeR (RRID:SCR_012802). First, lowly expressed RNAs with a count per million (CPM) <1 in more than 20% of the total samples were excluded. Then, a trimmed mean of M-values (TMM) normalization was performed. The biological coefficient of variation (BCV) plot was applied to examine the outliers and overall relationships of the samples. The quantile-adjusted conditional maximum likelihood (qCML) method was used for the estimation of dispersions. Differential expression was then determined using the exact test. The exact test employed in EdgeR is a method that was developed for differential expression, which is highly parallel to Fisher’s exact test. FDR < 0.05 and fold change (FC) larger than 1.5 (| log₂FC| > 0.58) were set as thresholds for significantly differentially expressed miRNAs, mRNAs, lncRNAs, and circRNAs.

Co-expression Module Analysis

Weighted gene co-expression network analysis (WGCNA) is a systems biology method for exploring the correlation patterns among genes across samples, which facilitates network-based identification of candidate biomarkers. A co-expression network of genes was constructed using the R package Weighted Gene co-expression Network Analysis (WGCNA, RRID:SCR_003302). The RNAs whose expression profiles were highly correlated across the samples were clustered as co-expression modules. The correlation of repression modules and MDD was then analyzed using a biweight midcorrelation method. The soft thresholding power was chosen based on the criterion of approximate scale-free topology. The minimal module size was set to 50. Highly correlated modules were merged using a cut height of 0.25. An absolute value of r larger than 0.5 was set as the threshold for a significant correlation between the expression modules and MDD.

Construction of ceRNA Networks

The software RNAhybrid (RRID:SCR_003252) + svm_light,¹ miRanda (RRID:SCR_017496) and TargetScan (RRID:SCR_010845) were combined for binding prediction between miRNA and lncRNA, miRNA, or mRNA. Differentially expressed lncRNAs and circRNAs were used as candidate ceRNAs. The candidate ceRNA-binding miRNAs were predicted, and the R package ceRNASeek was used to calculate ceRNA interactions (34). Then, miRNA-targeted mRNAs were predicted. Cytoscape (RRID:SCR_003032) was used for the ceRNA network visualization.

Kyoto Encyclopedia of Genes and Genomes Pathway Analyses

The RNA transcripts identified in the differential expression analysis, co-expression module analysis and ceRNA network analysis were subjected to Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. The significantly enriched pathways were subjected to KEGG pathway network analysis to identify relationships between pathways. KEGG pathway enrichment analysis and KEGG pathway network analysis were performed using the online Omicshare Tools of Gene Denove Biotechnology Co.²

Statistical Analysis

Statistical analysis of demographic data was performed using IBM SPSS (version 22.0). The independent-samples T-test was used to compare age, education years and body mass index (BMI) between the MDD and HC groups.

Results

Clinical and Demographic Characteristics

Fifteen patients with MDD (6 males and 9 females) and 15 HCs (6 males and 9 females) were included in this study. There was no significant difference in age (t = 1.43, p = 0.16), education (t = 1.52, p = 0.13), or BMI (t = 0.68, p = 0.50) between MDD and HCs (Table 1). The average age of onset (year), disease duration (month), and HAMD and HAMA scores in MDD were 21.1 ± 5.2, 23.3 ± 25.3, 21.1 ± 4.2, and 14.1 ± 6.1, respectively (Table 1).

TABLE 1

Table 1. Clinical characteristics of patients with MDD and HCs.

Summary of RNA Sequencing and Transcriptome Profiles

RNA sequencing generated 97.71 ± 7.90 G reads for MDD and 89.87 ± 23.72 G reads for HCs on average. After filtering, 96.61 ± 7.78 and 88.92 ± 23.49 G high-quality clean reads remained for MDD and HCs, respectively. There was no significant difference in raw reads and clean reads between MDD and HCs (Table 2).

TABLE 2

Table 2. Summary of RNA sequencing data.

Small RNA sequencing generated 12.77 ± 2.57 G reads for MDD and 11.54 ± 1.75 G reads for HCs on average. After filtering, 12.49 ± 2.52 and 11.27 ± 1.68 G high-quality clean reads remained for MDD and HCs, respectively. There was no significant difference in row reads and clean reads between MDD and HCs (Table 2).

By alignment to the reference genome or miRBase database, 2,660 mature miRNAs (1,181 exist, 1,085 known and 394 novel), 26,788 lncRNAs (24,743 known and 2,045 novel), 53,545 circRNAs (11,734 exist and 41,811 novel) and 100,769 mRNA transcripts were identified, which were used for the subsequent analyses.

Differentially Expressed RNAs in Major Depressive Disorder

After filtering lowly expressed RNAs, 513 miRNAs, 3,369 lncRNAs, 5,240 circRNAs and 23,453 mRNAs remained for differential expression analysis. The overall BCV plot suggested no outlier in the samples (Figures 2A,D). In the analysis of RNA sequencing data, the value of BCV was 0.39, and the common dispersion was 0.14 (Figure 2B). In the analysis of small RNA sequencing data, the value of BCV was 0.25, and the common dispersion was 0.51 (Figure 2E).

FIGURE 2

Figure 2. Differential expression analysis of RNA transcripts. The RNA package EdgeR was employed for differential expression analysis of RNA transcripts between major depressive disorder (MDD) and healthy controls (HCs). (A,D) The plot of biological coefficient of variation (BCV) distance of all the samples for RNA sequencing and small RNA sequencing, respectively; (B,E) the BCV plot of dispersion estimates for RNA sequencing and small RNA sequencing, respectively; (C,F) plots of log-fold change against log-counts per million for RNA sequencing and small RNA sequencing, respectively. The differentially expressed RNAs are highlighted.

By using the exact test, 7 miRNAs (2 up-regulated and 5 down-regulated), 25 lncRNAs (20 up-regulated and 5 down-regulated), 16 circRNAs (11 up-regulated and 5 down-regulated) and 956 mRNA transcripts (439 up-regulated and 517 down-regulated) were identified as significantly dysregulated in MDD compared to HCs (Figures 2C,F). The top 10 differentially expressed miRNAs, lncRNAs, circRNAs, and mRNAs are presented in Table 3.

TABLE 3

Table 3. Top 10 differential miRNAs, circRNAs, lncRNAs, and mRNAs between MDD and HCs.

Identification of Differential Co-expression Modules

Highly correlated RNAs in expression profiles were clustered as expression modules using the R package WGCNA. Based on the criterion of approximate scale-free topology, 4 was selected as the soft thresholding power (Figure 3A). Eighteen expression modules of RNAs were identified (Figure 3B). Then, the correlation between the expression modules and MDD was analyzed. The purple (r = –0.55), turquoise (r = –0.63), blue (r = 0.71), and yellow (r = 0.55) modules were significantly correlated with MDD (Figure 3C).

FIGURE 3

Figure 3. Identification of co-expression modules of RNA transcripts. (A) Analysis of RNA network topology for various soft-thresholding powers; (B) RNA expression module clustering dendrograms based on topological overlap; (C) association of RNA expression modules with major depressive disorder (MDD).

ceRNA Networks

The differentially expressed lncRNAs and circRNAs were subjected to ceRNA interaction analysis. Four lncRNAs and three circRNAs were identified to interact and used as ceRNAs for ceRNA network construction. The ceRNA network indicated that 7 ceRNAs may bind to 14 miRNAs and regulate the expression of 33 mRNAs (Figure 4).

FIGURE 4

Figure 4. Construction of ceRNA networks. A ceRNA network was constructed using 7 ceRNAs (3 circRNAs and 4 lncRNAs), 14 miRNAs and 33 mRNA transcripts.

Kyoto Encyclopedia of Genes and Genomes Pathway Analysis

To help interpret differential RNA transcript profiles in the context of biological processes between MDD and HCs, we performed KEGG pathway analysis. The results showed that the differentially expressed RNAs were enriched in multiple pathways, i.e., Parkinson’s disease, ribosome, thermogenesis, oxidative phosphorylation, renal cell carcinoma, hepatocellular carcinoma, non-alcoholic fatty liver disease (NAFLD), chemokine signaling pathway, and Jak-STAT signaling pathway (Figure 5A and Supplementary Table 1). The RNA transcripts in the blue co-expression module were enriched in four pathways, i.e., chemokine signaling pathway, T-cell receptor signaling pathway, leukocyte transendothelial migration, and sulfur metabolism; the RNA transcripts in the yellow co-expression module were enriched in the porphyrin and chlorophyll metabolism pathways. No RNAs in the turquoise and purple modules were found to be significantly enriched in KEGG pathways (Figures 5B–E and Supplementary Table 2). The mRNAs regulated by ceRNA networks were enriched in the glucagon signaling pathway, insulin resistance, the AMPK signaling pathway, and fructose and mannose metabolism (Figure 5F and Supplementary Table 3).

FIGURE 5

Figure 5. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of RNA transcripts. The top 10 pathways in KEGG pathway enrichment analysis.

KEGG pathway network analysis was performed to explore the relationship between the significant pathways identified in KEGG pathway enrichment analysis. The results indicated that two networks were identified. One network contains six pathways, i.e., glucagon signaling pathway, thermogenesis, Parkinson’s disease, oxidative phosphorylation, non-alcoholic fatty liver disease (NAFLD) and hepatocellular carcinoma; the oxidative phosphorylation pathway had the largest degree in this pathway network (Figure 6 left), which was critical for energy metabolism. Another network contains three pathways, i.e., Jak-STAT signaling pathway, chemokine signaling pathway, and leukocyte transendothelial migration, the chemokine signaling pathway had the largest degree in this pathway network (Figure 6 right), which suggested that inflammation might involve in MDD.

FIGURE 6

Figure 6. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway networks of dysregulated pathways in major depressive disorder (MDD). The circles indicate KEGG pathways. The edge between two pathways indicates that there is a relationship between the two pathways. The color of the circles indicates the degree of the pathways, i.e., how many edges linking to this pathways.

Discussion

In this study, a whole transcriptome sequencing was performed on the peripheral blood of MDD and matched HCs. The differentially expressed mRNAs, lncRNAs, miRNAs, and circRNAs, the dysregulated co-expression modules of RNA transcripts and the ceRNA network in MDD compared to HCs were identified. We also explored biological processes associated with dysregulation of transcriptome profiles by employing KEGG pathway analysis.

The ceRNA network indicated that 7 ceRNAs, including 4 lncRNAs and 3 circRNAs, might bind to 14 miRNAs and regulate the expression of 33 mRNAs. The identified lncRNA/circRNA-associated ceRNA network would offer novel insights into the pathophysiology of MDD and thus suggest novel treatments for the disease. In this network, hsa-let-7a-5p, hsa-let-7b-5p, hsa-let-7d-5p, hsa-let-7g-5p, and hsa-let-7i-5p belong to the let-7 family, which is the most highly expressed miRNA in the human brain and exerts a pivotal action on neuronal differentiation and maturation during neurodevelopment and on neurogenesis and neuronal plasticity functions in the adult brain (35, 36). Previous studies showed that hsa-let-7a-5p, hsa-let-7d-5p, and hsa-let-7g-5p were significantly dysregulated in MDD (37, 38) and gradually recovered to normal after 12 weeks of treatment with escitalopram, which suggested that the three miRNAs may be involved in both the pathogenesis of MDD and the effects of antidepressant drugs (37, 39). This is consistent with accumulated evidence that depression is associated with synaptic deficits, while neuroinflammation may influence neuronal activity as well as synaptic plasticity and contribute to neurobiological changes underlying MDD (40, 41). Second, let-7 is a family of miRNAs known to be involved in metabolic control (42). Previous studies showed that hsa-let-7a-5p was dysregulated in individuals with prediabetes and T2D (43), and the levels increased after 12 months of diabetic medication (44). The expression of hsa-let-7b-5p, hsa-let-7i-5p, and hsa-miR-92a-3p was significantly decreased in obese patients after bariatric surgery (45). Plasma hsa-let-7d-5p, hsa-let-7g-5p, and hsa-let-7i-5p were significantly dysregulated in Alzheimer’s disease (AD), and their target genes involved in the disturbance of multiple enzymatic pathways, including lipid metabolism, could play a role in AD etiology (46–48). Let-7b is also involved in inflammation and the immune response by regulating the activation of NF-kB and mediating the down-regulation of IL-6 in macrophages (49, 50). Researchers reported that hsa-let-7b-5p was significantly dysregulated in patients with active Crohn’s disease compared with healthy controls and was correlated with serum C-reactive protein (CRP) levels (51, 52). Finally, other miRNAs were also associated with metabolism, inflammation and immunity, and oxidative phosphorylation. For instance, hsa-miR-25-3p and hsa-miR-92a-3p are related to type 1 diabetes and systemic lupus erythematosus (SLE), respectively (53, 54). While hsa-miR-199a-3p was shown to be able to effectively distinguish schizophrenia or Parkinson’s disease (PD) from healthy controls, its target genes were mostly enriched in the FoxO signaling pathway, which is a potential mediator in the pathogenesis of depression (55–57). These pieces of evidence suggested that RNAs in the network were involved in sugar and lipid metabolism, inflammation and immunity, and oxidative phosphorylation, which might play an important role in the occurrence and development of depression.

Interpretation of biological readouts of transcriptome profiles is one of the challenges in transcriptome analyses. In this study, KEGG pathway network analysis was applied for intergradation of the results of KEGG pathway enrichment analysis, which was helpful for exploration of the potential biological mechanism underlying MDD. Two pathway networks were identified, and the oxidative phosphorylation pathway and chemokine signaling pathway have the most connection to other pathways in the two networks, which may suggest the critical roles of the two pathways in the pathophysiology of MDD. The oxidative phosphorylation pathway is widely involved in energy metabolism processes. It is known that the brain requires a continuous supply of energy in the form of ATP, most of which is predominantly generated within mitochondria by oxidative phosphorylation of glucose via the tricarboxylic acid (TCA) cycle (58). Previous studies also indicated that alterations in the oxidative phosphorylation system might represent the pathophysiological basis of MDD, and mitochondrial energy metabolism could be considered one of the major mechanisms underlying chronic corticosterone (CORT)-mediated depression (59, 60). And abnormal energy metabolism is considered to be one of the important mechanisms for the occurrence and development of depression (61). In addition, some researchers suggested that MDD and AD might share a common pathological basis, and MDD was considered to be a precursor of AD (62). Evidence has been accumulating that impaired brain energetics are involved in the cause and progression of neurodegenerative disorders of aging, especially AD (63, 64). The chemokine signaling pathway is a critical pathway that participates in inflammation. Accumulating evidence has confirmed the association between chronic inflammation and MDD, which is manifested by increased levels of multiple proinflammatory cytokines, such as interleukin (IL)-2, IL-6, CCL4, and CCL8 (65). Despite the considerable heterogeneity of experimental samples and methodologies, previous studies have demonstrated the important roles of chemokines in the pathophysiology of MDD (66). Chemokines could be considered potential peripheral markers and treatment targets for MDD. Therefore, the roles of oxidative phosphorylation-related energy metabolism and chemokine signaling pathway-related inflammation in MDD deserve refining study in the future.

Our research has some limitations. First, the sample size was relatively small due to the high costs of sequencing, which may result in insufficient sample representativeness. Second, we could not exclude the long-term effect of antidepressants on the results even though patients had not taken any antidepressant medication for at least 3 months before enrollment. Therefore, large sample studies are warranted in the future.

Conclusion

In conclusion, our results suggested that RNAs in peripheral blood mainly affected metabolism, inflammation and immunity, and oxidative phosphorylation, thereby being involved in the pathophysiology of MDD. The results will be helpful for better understanding the pathophysiology of MDD and provide novel clues for exploring molecular biomarkers.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics Statement

The studies involving human participants were reviewed and approved by the Ethics Committee of West China Hospital of Sichuan University. The patients/participants provided their written informed consent to participate in this study.

Author Contributions

YW, JW, TC, TL, and XM contributed to the conception and design of the review. YW, XY, MW, YKD, and YD contributed to the data collection. YW, JW, LZ, and RN performed the statistical analysis. YW and JW wrote the first draft of the manuscript. XM critically revised the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This study was supported by the Key Research and Development Program of Science and Technology Department of Sichuan Province (Nos. 22ZDYF1531 and 22ZDYF1696), the Program of Chengdu Science and Technology (No. 2021-YF05-00272-SN), the National Natural Science Foundation of China (No. 82001432), the China Postdoctoral Science Foundation (Nos. 2020TQ0213 and 2020M683319), the Open Project Program of the National Laboratory of Pattern Recognition (No. 202000034), and the West China Hospital Postdoctoral Science Foundation (No. 2020HXBH104).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We would like to express our sincere thanks to all the hospital staff and the patients for their contributions to the study.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpsyt.2022.907034/full#supplementary-material

Footnotes

1. ^ https://www.cs.cornell.edu/people/tj/svm_light/old/svm_light_v6.01.html
2. ^ https://www.omicshare.com/tools

References

1. World Health Organization. Depression and Other Common Mental Disorders: Global Health Estimates. Geneva: World Health Organization (2017).

Google Scholar

2. Belmaker RH, Agam G. Major depressive disorder. N Engl J Med. (2008) 358:55–68.

Google Scholar

3. Altamura AC, Buoli M, Albano A, Dell’Osso B. Age at onset and latency to treatment (duration of untreated illness) in patients with mood and anxiety disorders: a naturalistic study. Int Clin Psychopharmacol. (2010) 25:172–9. doi: 10.1097/YIC.0b013e3283384c74

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Huerta-Ramírez R, Bertsch J, Cabello M, Roca M, Haro JM, Ayuso-Mateos JL. Diagnosis delay in first episodes of major depression: a study of primary care patients in Spain. J Affect Disord. (2013) 150:1247–50. doi: 10.1016/j.jad.2013.06.009

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Kemp AH, Gordon E, Rush AJ, Williams LM. Improving the prediction of treatment response in depression: integration of clinical, cognitive, psychophysiological, neuroimaging, and genetic measures. CNS Spectr. (2008) 13:1066–1086quiz1087–8. doi: 10.1017/s1092852900017120

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Rush AJ, Trivedi MH, Wisniewski SR, Stewart JW, Nierenberg AA, Thase ME, et al. Bupropion-SR, sertraline, or venlafaxine-XR after failure of SSRIs for depression. N Engl J Med. (2006) 354:1231–42. doi: 10.1056/NEJMoa052963

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Thase ME, Nierenberg AA, Vrijland P, van Oers HJ, Schutte AJ, Simmons JH. Remission with mirtazapine and selective serotonin reuptake inhibitors: a meta-analysis of individual patient data from 15 controlled trials of acute phase treatment of major depression. Int Clin Psychopharmacol. (2010) 25:189–98. doi: 10.1097/YIC.0b013e328330adb2

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Simon GE, Perlis RH. Personalized medicine for depression: can we match patients with treatments? Am J Psychiatry. (2010) 167:1445–55. doi: 10.1176/appi.ajp.2010.09111680

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Kendler KS, Ohlsson H, Sundquist K, Sundquist J. Sources of parent-offspring resemblance for major depression in a National Swedish extended adoption study. JAMA Psychiatry. (2018) 2:194–200. doi: 10.1001/jamapsychiatry.2017.3828

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Sullivan PF, Neale MC, Kendler KS. Genetic epidemiology of major depression: review and meta-analysis. Am J Psychiatry. (2000) 157:1552–62. doi: 10.1176/appi.ajp.157.10.1552

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Velculescu VE, Zhang L, Zhou W, Vogelstein J, Basrai MA, Bassett DE, et al. Characterization of the yeast transcriptome. Cell. (1997) 88:243–51. doi: 10.1016/s0092-8674(00)81845-0

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Zhang L, Chen X, Yin J. Prediction of potential miRNA-disease associations through a novel unsupervised deep learning framework with variational autoencoder. Cells. (2019) 8:1040. doi: 10.3390/cells8091040

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Li X, Zuo H, Zhang L, Sun Q, Xin Y, Zhang L. HMMRValidating expression and its prognostic significance in lung adenocarcinoma based on data mining and bioinformatics methods. Front Oncol. (2021) 11:720302. doi: 10.3389/fonc.2021.720302

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Qi YX, Liu YB, Rong WH. RNA-Seq and its applications: a new technology for transcriptomics. Yi Chuan. (2011) 33:1191–202. doi: 10.3724/sp.j.1005.2011.01191

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Shi Y, Wang Q, Song R, Kong Y, Zhang Z. Non-coding RNAs in depression: promising diagnostic and therapeutic biomarkers. EBioMedicine. (2021) 71:103569. doi: 10.1016/j.ebiom.2021.103569

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Liu W, Yan H, Zhou D, Cai X, Zhang Y, Li S, et al. The depression GWAS risk allele predicts smaller cerebellar gray matter volume and reduced SIRT1 mRNA expression in Chinese population. Transl Psychiatry. (2019) 9:333. doi: 10.1038/s41398-019-0675-3

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Chen L, Li XS, Zheng GE, Xie GJ, Cheng Y. Peripheral blood BDNF-TrkB signaling in first-episode, drug-free patients with major depressive disorder and schizophrenia. Neurosci Lett. (2020) 714:134618. doi: 10.1016/j.neulet.2019.134618

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Farmer R, Burbano SD, Patel NS, Sarmiento A, Smith AJ, Kelly MP. Phosphodiesterases PDE2A and PDE10A both change mRNA expression in the human brain with age, but only PDE2A changes in a region-specific manner with psychiatric disease. Cell Signal. (2020) 70:109592. doi: 10.1016/j.cellsig.2020.109592

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Lopez JP, Lim R, Cruceanu C, Crapper L, Fasano C, Labonte B, et al. miR-1202 is a primate-specific and brain-enriched microRNA involved in major depression and antidepressant treatment. Nat Med. (2014) 20:764–8. doi: 10.1038/nm.3582

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Garbett KA, Vereczkei A, Kálmán S, Brown JA, Taylor WD, Faludi G, et al. Coordinated messenger RNA/MicroRNA changes in fibroblasts of patients with major depression. Biol Psychiatry. (2015) 77:256–65. doi: 10.1016/j.biopsych.2014.05.015

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Baudry A, Mouillet-Richard S, Schneider B, Launay JM, Kellermann O. MiR-16 targets the serotonin transporter: a new facet for adaptive responses to antidepressants. Science. (2010) 329:1537–41. doi: 10.1126/science.1193692

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Cui X, Sun X, Niu W, Kong L, He M, Zhong A, et al. Long non-coding RNA: potential diagnostic and therapeutic biomarker for major depressive disorder. Med Sci Monit. (2016) 22:5240–8. doi: 10.12659/msm.899372

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Ni X, Liao Y, Li L, Zhang X, Wu Z. Therapeutic role of long non-coding RNA TCONS_00019174 in depressive disorders is dependent on Wnt/β-catenin signaling pathway. J Integr Neurosci. (2018) 17:203–15. doi: 10.31083/JIN-170052

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Zhang Y, Du L, Bai Y, Han B, He C, Gong L, et al. CircDYM ameliorates depressive-like behavior by targeting miR-9 to regulate microglial activation via HSP90 ubiquitination. Mol Psychiatry. (2020) 25:1175–90. doi: 10.1038/s41380-018-0285-0

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Cui X, Niu W, Kong L, He M, Jiang K, Chen S, et al. hsa_circRNA_103636: potential novel diagnostic and therapeutic biomarker in major depressive disorder. Biomark Med. (2016) 10:943–52. doi: 10.2217/bmm-2016-0130

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Guil S, Esteller M. RNA-RNA interactions in gene regulation: the coding and noncoding players. Trends Biochem Sci. (2015) 40:248–56. doi: 10.1016/j.tibs.2015.03.001

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Salmena L, Poliseno L, Tay Y, Kats L, Pandolfi PP. A ceRNA hypothesis: the rosetta stone of a hidden RNA language? Cell. (2011) 146:353–8. doi: 10.1016/j.cell.2011.07.014

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Tay Y, Rinn J, Pandolfi PP. The multilayered complexity of ceRNA crosstalk and competition. Nature. (2014) 505:344–52. doi: 10.1038/nature12986

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Zuo H, Li X, Zheng X, Sun Q, Yang Q, Xin Y. A Novel circRNA-miRNA-mRNA hub regulatory network in lung adenocarcinoma. Front Genet. (2021) 12:673501. doi: 10.3389/fgene.2021.673501

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Yaoxian G. Wechsler Adult Intelligence Scale-Revised in China Version. Changsha: Hunan Medical College (1992).

Google Scholar

31. Liang S, Brown MRG, Deng W, Wang Q, Ma X, Li M, et al. Convergence and divergence of neurocognitive patterns in schizophrenia and depression. Schizophr Res. (2018) 192:327–34. doi: 10.1016/j.schres.2017.06.004

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Hamilton M. A rating scale for depression. J Neurol Neurosurg Psychiatry. (1960) 23:56–62.

Google Scholar

33. Hamilton M. The assessment of anxiety states by rating. Br J Med Psychol. (1959) 32:50–5. doi: 10.1111/j.2044-8341.1959.tb00467.x

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Zhang MY, Jin XY, Lin JY, Tian Y, Wang Q, Li XH, et al. CeRNASeek: an R package for identification and analysis of ceRNA regulation. Brief Bioinfor. (2021) 22:bbaa048. doi: 10.1093/bib/bbaa048

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Anacker AM, Beery AK. Life in groups: the roles of oxytocin in mammalian sociality. Front Behav Neurosci. (2013) 7:185. doi: 10.3389/fnbeh.2013.00185

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Shao NY, Hu HY, Yan Z, Xu Y, Hu H, Menzel C, et al. Comprehensive survey of human brain microRNA by deep sequencing. BMC Genomics. (2010) 11:409. doi: 10.1186/1471-2164-11-409

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Maffioletti E, Cattaneo A, Rosso G, Maina G, Maj C, Gennarelli M, et al. Peripheral whole blood microRNA alterations in major depression and bipolar disorder. J Affect Disord. (2016) 200:250–8. doi: 10.1016/j.jad.2016.04.021

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Mendes-Silva AP, Pereira KS, Tolentino-Araujo GT, Nicolau Ede S, Silva-Ferreira CM, Teixeira AL, et al. Shared biologic pathways between alzheimer disease and major depression: a systematic review of microrna expression studies. Am J Geriatr Psychiatry. (2016) 24:903–12. doi: 10.1016/j.jagp.2016.07.017

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Bocchio-Chiavetto L, Maffioletti E, Bettinsoli P, Giovannini C, Bignotti S, Tardito D, et al. Blood microRNA changes in depressed patients during antidepressant treatment. Eur Neuropsychopharmacol. (2013) 23:602–11. doi: 10.1016/j.euroneuro.2012.06.013

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Duman RS, Aghajanian GK, Sanacora G, Krystal JH. Synaptic plasticity and depression: new insights from stress and rapid-acting antidepressants. Nat Med. (2016) 22:238–49. doi: 10.1038/nm.4050

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Wohleb ES, Franklin T, Iwata M, Duman RS. Integrating neuroimmune systems in the neurobiology of depression. Nat Rev Neurosci. (2016) 17:497–511. doi: 10.1038/nrn.2016.69

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Jiang S, Yan W, Wang SE, Baltimore D. Let-7 suppresses B-cell activation through restricting the availability of necessary nutrients. Cell Metab. (2018) 27:393–403. doi: 10.1016/j.cmet.2017.12.007

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Mononen N, Lyytikäinen L, Seppälä I, Mishra PP, Juonala M, Waldenberger M, et al. Whole blood microRNA levels associate with glycemic status and correlate with target mRNAs in pathways important to type 2 diabetes. Sci Rep. (2019) 9:8887. doi: 10.1038/s41598-019-43793-4

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Santovito D, De Nardis V, Marcantonio P, Mandolini C, Paganelli C, Vitale E, et al. Plasma exosome microRNA profiling unravels a new potential modulator of adiponectin pathway in diabetes: effect of glycemic control. J Clin Endocrinol Metab. (2014) 99:E1681–5. doi: 10.1210/jc.2013-3843

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Langi G, Szczerbinski L, Kretowski A. Meta-analysis of differential miRNA expression after bariatric surgery. J Clin Med. (2019) 8:1220. doi: 10.3390/jcm8081220

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Kumar P, Dezso Z, MacKenzie C, Oestreicher J, Agoulnik S, Byrne M, et al. Circulating miRNA biomarkers for Alzheimer’s disease. PLoS One. (2013) 8:e69807. doi: 10.1371/journal.pone.0069807

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Chen J, Qi Y, Liu C, Lu J, Shi J, Shi Y. MicroRNA expression data analysis to identify key miRNAs associated with Alzheimer’s disease. J Gene Med. (2018) 20:e3014. doi: 10.1002/jgm.3014

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Gámez-Valero A, Campdelacreu J, Vilas D, Ispierto L, Reñé R, Álvarez R, et al. Exploratory study on microRNA profiles from plasma-derived extracellular vesicles in Alzheimer’s disease and dementia with Lewy bodies. Transl Neurodegener. (2019) 8:31. doi: 10.1186/s40035-019-0169-5

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Li D, Jia H, Zhang H, Lv M, Liu J, Zhang Y, et al. TLR4 signaling induces the release of microparticles by tumor cells that regulate inflammatory cytokine IL-6 of macrophages via microRNA let-7b. Oncoimmunology. (2012) 1:687–93. doi: 10.4161/onci.19854

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Teng GG, Wang WH, Dai Y, Wang SJ, Chu YX, Li J. Let-7b is involved in the inflammation and immune responses associated with Helicobacter pylori infection by targeting Toll-like receptor 4. PLoS One. (2013) 8:e56709. doi: 10.1371/journal.pone.0056709

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Guo Z, Gong J, Li Y, Gu L, Cao L, Wang Z, et al. Mucosal MicroRNAs expression profiles before and after exclusive enteral nutrition therapy in adult patients with Crohn’s Disease. Nutrients. (2016) 8:519. doi: 10.3390/nu8080519

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Guo Z, Wu R, Gong J, Zhu W, Li Y, Wang Z, et al. Altered microRNA expression in inflamed and non-inflamed terminal ileal mucosa of adult patients with active Crohn’s disease. J Gastroenterol Hepatol. (2015) 30:109–16. doi: 10.1111/jgh.12644

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Nielsen LB, Wang C, Sørensen K, Bang-Berthelsen CH, Hansen L, Andersen ML, et al. Circulating levels of microRNA from children with newly diagnosed type 1 diabetes and healthy controls: evidence that miR-25 associates to residual beta-cell function and glycaemic control during disease progression. Exp Diabetes Res. (2012) 2012:896362. doi: 10.1155/2012/896362

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Kim BS, Jung JY, Jeon JY, Kim HA, Suh CH. Circulating hsa-miR-30e-5p, hsa-miR-92a-3p, and hsa-miR-223-3p may be novel biomarkers in systemic lupus erythematosus. HLA. (2016) 88:187–93. doi: 10.1111/tan.12874

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Pala E, Denkçeken T. Evaluation of miRNA expression profiles in schizophrenia using principal-component analysis-based unsupervised feature extraction method. J Comput Biol. (2020) 27:1253–63. doi: 10.1089/cmb.2019.0412

PubMed Abstract | CrossRef Full Text | Google Scholar

56. He S, Huang L, Shao C, Nie T, Xia L, Cui B, et al. Several miRNAs derived from serum extracellular vesicles are potential biomarkers for early diagnosis and progression of Parkinson’s disease. Transl Neurodegener. (2021) 10:25. doi: 10.1186/s40035-021-00249-y

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Polter A, Yang S, Zmijewska AA, van Groen T, Paik JH, Depinho RA, et al. Forkhead box, class O transcription factors in brain: regulation and behavioral manifestation. Biol Psychiatry. (2009) 65:150–9. doi: 10.1016/j.biopsych.2008.08.005

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Bordone MP, Salman MM, Titus HE, Amini E, Andersen JV, Chakraborti B, et al. The energetic brain – a review from students to students. J Neurochem. (2019) 151:139–65.

Google Scholar

59. Kuffner K, Triebelhorn J, Meindl K, Benner C, Manook A, Sudria-Lopez D, et al. Major depressive disorder is associated with impaired mitochondrial function in skin fibroblasts. Cells. (2020) 9:884. doi: 10.3390/cells9040884

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Gong Q, Yan XJ, Lei F, Wang ML, He LL, Luo YY, et al. Proteomic profiling of the neurons in mice with depressive-like behavior induced by corticosterone and the regulation of berberine: pivotal sites of oxidative phosphorylation. Mol Brain. (2019) 12:118. doi: 10.1186/s13041-019-0518-4

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Gu X, Ke S, Wang Q, Zhuang T, Xia C, Xu Y, et al. Energy metabolism in major depressive disorder: recent advances from omics technologies and imaging. Biomed Pharmacother. (2021) 141:111869. doi: 10.1016/j.biopha.2021.111869

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Gatchel JR, Rabin JS, Buckley RF, Locascio JJ, Quiroz YT, Yang HS, et al. Longitudinal association of depression symptoms with cognition and cortical amyloid among community-dwelling older adults. JAMA Netw Open. (2019) 2:e198964. doi: 10.1001/jamanetworkopen.2019.8964

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Camandola S, Mattson MP. Brain metabolism in health, aging, and neurodegeneration. EMBO J. (2017) 36:1474–92. doi: 10.15252/embj.201695810

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Johnson ECB, Dammer EB, Duong DM, Ping L, Zhou M, Yin L, et al. Large-scale proteomic analysis of Alzheimer’s disease brain and cerebrospinal fluid reveals early changes in energy metabolism associated with microglia and astrocyte activation. Nat Med. (2020) 26:769–80. doi: 10.1038/s41591-020-0815-6

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Syed SA, Beurel E, Loewenstein DA, Lowell JA, Craighead WE, Dunlop BW, et al. Defective Inflammatory pathways in never-treated depressed patients are associated with poor treatment response. Neuron. (2018) 99:914–24. doi: 10.1016/j.neuron.2018.08.001

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Leighton SP, Nerurkar L, Krishnadas R, Johnman C, Graham GJ, Cavanagh J. Chemokines in depression in health and in inflammatory illness: a systematic review and meta-analysis. Mol Psychiatry. (2018) 23:48–58. doi: 10.1038/mp.2017.205

PubMed Abstract | CrossRef Full Text | Google Scholar