Vol 74, No 3 (2023)
Original paper
Published online: 2023-05-04

open access

Page views 1511
Article views/downloads 408
Get Citation

Connect on Social Media

Connect on Social Media

Original paper

Endokrynologia Polska

DOI: 10.5603/EP.a2023.0030

ISSN 0423–104X, e-ISSN 2299–8306

Volume/Tom 74; Number/Numer 3/2023

Submitted: 25.01.2023

Accepted: 25.02.2023

Early publication date: 04.05.2023

Leukocyte transcriptome of Cushing’s disease are associated with nerve impairment and psychiatric disorders

Kunyu He1*2Tao Zhou*1Fuyu Wang1Fangye Li1Yanyang Zhang1Xinguang Yu1
1Department of Neurosurgery, The First Medical Centre of Chinese PLA General Hospital, Beijing, China
2Chinese PLA Medical School, Beijing, China
*These authors contributed equally.

Yanyang Zhang, Department of Neurosurgery, The First Medical Centre of Chinese PLA General Hospital, 28 Fuxing Road, Haidian District, Beijing, China, tel: +86 18632218218; e-mail: sjwkzyy@163.com
Xinguang Yu, Department of Neurosurgery, The First Medical Center of Chinese PLA General Hospital, 28 Fuxing Road, Haidian District, Beijing, China, tel: +86 19933309008; e-mail: yuxinguang_301@163.com

This article is available in open access under Creative Common Attribution-Non-Commercial-No Derivatives 4.0 International (CC BY-NC-ND 4.0) license, allowing to download articles and share them with others as long as they credit the authors and the publisher, but without permission to change them in any way or use them commercially

Introduction: The hypothalamus-pituitary-adrenal (HPA) axis and its end product cortisol is a major response mechanism to stress and plays a critical role in many psychiatric disorders. Cushing’s disease (CD) serves as a valuable in vivo “hyperexpression” model to elucidate the effect of cortisol on brain function and mental disorders. Changes in brain macroscale properties measured by magnetic resonance imaging (MRI) have been detailed demonstrated, but the biological and molecular mechanisms underlying these changes remain poorly understood.
Material and methods: Here we included 25 CD patients and matched 18 healthy controls for assessment, and performed transcriptome sequencing of peripheral blood leukocytes. Weighted gene co-expression network analysis (WGCNA) was performed to construct a co-expression network of the relationships between genes and we identified a significant module and hub gene types associated with neuropsychological phenotype and psychiatric disorder identified in enrichment analysis. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis preliminarily explored the biological functions of these modules.
Results: The WGCNA and enrichment analysis indicated that module 3 of blood leukocytes was enriched in broadly expressed genes and was associated with neuropsychological phenotypes and mental diseases enrichment. GO and KEGG enrichment analysis of module 3 identified enrichment in many biological pathways associated with psychiatric disorders.
Conclusion: Leukocyte transcriptome of Cushing’s disease is enriched in broadly expressed genes and is associated with nerve impairment and psychiatric disorders, which may reflect some changes in the affected brain. (Endokrynol Pol 2023; 74 (3): 294–304)
Key words: Cushing syndrome; pituitary ACTH hypersecretion; leukocytes; transcriptome; mental disorders


Repetitive stress exposure plays a critical role in many psychiatric disorders, such as schizophrenia [1], bipolar disorder [2], and anxiety disorders [3]. The hypothalamus–pituitary–adrenal (HPA) axis and its end product cortisol comprise a major response mechanism to stress [4]. It is generally assumed that chronic activation or hyperresponsivity of the HPA axis plays a crucial role in the relationship between stress and the aetiology of multiple psychiatric disorders [1, 5]. Prolonged exposure to excess cortisol causes various deleterious effects on the vulnerable brain rich in cortisol receptors, leading to severe neurocognitive impairment and neuropsychological symptoms [6, 7]. Underlying these symptoms are objective alterations in brain structure and brain connectivity measured by magnetic resonance imaging (MRI) [8–11], such as decreased grey matter volume, altered morphology, and impaired intrinsic quality of resting state connectivity. However, because it is almost impossible to perform biopsies on the living brains of patients, the biological and molecular mechanisms underlying these changes remain poorly understood. Thus, peripheral blood gene expression study becomes a necessary proxy, despite obvious limitations.

Peripheral blood samples are widely used in the study of psychiatric disorders, but more so to search for specific differentially expressed genes for use as biomarkers to aid diagnosis [12, 13]. It has been demonstrated in studies comparing blood-brain transcriptomes, that the transcriptional profiles between them show significant overlap, ranging from 20% [14] to 55% [15]. Because cortisol can affect the whole body via a similar pathway, changes in the blood leukocyte transcriptome may reflect changes in the brain [16–18]. Thus, seeking associations between the blood leukocyte transcriptome and neurological phenotype may help elucidate how HPA affects the brain and leads to mental illness. The omnigenic model [19] predicts that in a complex trait, broadly expressed genes in one or more tissues, like the brain and blood, contain a large amount of the heritability signal and may contribute more to overall risk than tissue-specific genes of small number. Thus, these genes do not necessarily need to be specific to brain tissue, implying the possibility of studying intracerebral mechanisms through peripheral blood. We considered that neural phenotypes relevant to psychiatric disorders could be explained through large-scale coordinated genes in the blood leukocyte.

Cushing’s disease (CD) is a common endocrine disorder characterized by the unique excess of endogenous cortisol [20], naturally serving as a valuable in vivo “hyperexpression” model to study the prolonged effect of cortisol on the brain [21]. In the current study, we included 25 CD patients and matched 18 healthy controls for neuroendocrine and neuropsychological assessment, and performed transcriptome sequencing of peripheral blood leukocytes. Weighted gene co-expression network analysis (WGCNA) is a powerful tool for constructing free-scale gene co-expression networks to examine the relationship between gene sets and clinical features [22]. Unlike differentially expressed gene (DEG) analysis, WGCNA can incorporate as many genes as possible for analysis and fits well with the idea of the “Omnigenic Model”. Thus, to better study mechanisms of brain changes in patients, we performed WGCNA to construct a co-expression network, and we identified a significant module and hub gene types associated with neuropsychological phenotype and psychiatric disorder identified in enrichment analysis. Gene ontology (GO) and Kyoto Encyclopaedia of Genes and Genomes (KEGG) enrichment analysis preliminarily explored these modules’ biological functions and molecular mechanisms.

Material and methods


This study was approved by the Ethics Committee of the Chinese PLA General Hospital, Beijing, China (ethics approval number: S2021-677-01) and strictly followed the principles in the Declaration of Helsinki. Patients were included in the study after signing informed consent. Forty-three subjects were recruited, including 25 CD patients with adrenocorticotropic hormone (ACTH) pituitary tumours (3 males and 22 females; mean age 41.88 years; age range 25–59 years) and 18 matched healthy control (HC) individuals (2 males and 16 females; mean age 41.44 years; age range 25–57 years) in terms of sex, age, and education. All participants were checked to ensure no present or past mental disorder or psychiatric drug exposure history. Diagnosis of CD patients was determined by both experienced endocrinologists and neurosurgeons according to the latest clinical guideline after excluding exogenous cortisol exposure [23]. Detailed diagnostic criteria could be found in the supplementary material. Because the status differences between CD patients and healthy control were too significant, data collection and analysis were not blinded.

Neuroendocrine and neuropsychological assessment

Functional status of HPA in CD patients was quantified by biochemical assessments including 24-h urinary free cortisol (UFC), serum cortisol (nmol/L), and plasma ACTH (pmol/L) at 0:00, 8:00, and 16:00. As comparisons, in healthy control individuals, cortisol and ACTH at 8:00, and 24-h UFC were measured. Specifically, according to the manufacturer’s instructions, ACTH was detected with an IMMULITE 2000 Analyzer (Siemens Healthcare Diagnostics Inc., LA, United States), and cortisol was measured using an ADVIA Centaur Analyzer (Siemens Healthcare Diagnostics, Tarrytown, NY, United States). All the assays and blood sample collection for RNA-seq were completed in a standardized examination workflow before surgery.

All participants CD patients and healthy controls received the Self-Rating Depression Scale (SDS, scores < 53: “normal”, 53–62: “mild depression”, 63–72 “moderate depression”, > 73 “severe depression”) and the Self-Rating Anxiety Scale (SAS, scores < 50: “normal”, 50–59: “mild anxiety”, 60–69 “moderate anxiety”, > 69 “severe anxiety”) to reflect their mental status. The references to the above questionnaire and the following R packages can be found in the supplementary material.

Blood sample collection, sequencing, and data processing

Four to six millilitres of blood were collected into EDTA-coated tubes from the median cubital vein of CD patients and HC group when they had no fever, cold, flu, infections, or use of medications for illnesses 72 h before the blood draw. Leukocytes were captured by centrifugation, permeabilized by RNAlater, stored at –80°C (store duration up to 24 months), and finally transcriptome sequenced simultaneously to avoid batch effects. Detailed methods for sequencing can be found in the supplementary material. Duplicated genes, sex chromosome, and genes of undefined function and labelled as novel transcripts were removed from the raw count matrix. A total of 12,107 genes with count-per-million (CPM) expression greater than 0.5 in more than 90% of samples were retained for subsequent analysis.

Weighted gene coexpression-network analysis

Using the WGCNA R package (v1.71), we acquired 12 modules of tightly co-expressed genes from 12,107 genes (Supplementary Files Tab. S1). Detailed methods for WGCNA can be found in the supplementary material and codes. Further WGCNA analysis was performed in each group to check the preservation of detected modules across groups at the soft power threshold of 26 (Supplementary Files Fig. S2). Then the preservation was assessed through module preservation and quality statistics, which were computed using the module Preservation () function in the WGCNA package [24]. The adjacency matrix of the network of all subjects was taken as the reference, and the dataset of either CD or HC group was selected as test data with 1000 permutations. The stability of the modules was tested through the statistics median rank and Zsummary. We also tested the preservation of hub leukocyte module of all subjects in the post-mortem brain of interested mental illness using the same method.

Analysis of module eigengenes (MEs) group differences were conducted with t-test, and modules showed group differences were defined as differentially expressed (DE) modules. Module-trait associations were estimated by the correlations between the modules and clinical traits. Pearson correlation coefficients were calculated between the MEs of each module and each clinical feature, allowing the identification of modules that were significantly associated with the feature. DEG analysis was also performed using the DESeq2 R package (v1.34.0) (Supplementary File Tab. S14). We randomly chose several DEGs and performed polymerase chain reaction (PCR) to validate the transcriptome sequencing results. Detailed PCR methods can be found in the supplementary material. The PCR results are consistent with the transcriptome sequencing results (Supplementary File Fig. S3).

Enrichment analyses

We included a series of gene lists obtained from other research, and we analysed the enrichment between genes from hub modules and these gene lists. For these enrichment analyses, a custom R code written by M.V.L. (https://github.com/mvlombardo/utils/blob/master/genelistOverlap.R), which computes hypergeometric p-values and enrichment odds ratios (ORs), was used. The background pool for these enrichment tests was always set to 12,107. Only enrichments passing a false discovery rate (FDR) < 0.05 were interpreted.

We first used 3 gene classes reported by Boyle et al. [19] in GTEx data to annotate DE modules. These classes were broadly expressed genes, whole-blood-specific genes, and lymphocyte-specific genes. The HPA axis is an important factor in the aetiology of multiple psychiatric disorders. In our research, CD patients showed higher SDS and SAS scores than the HC group. We then performed enrichment analysis of some stress- and cortisol-related diseases. To test for enrichment between hub modules and gene sets of schizophrenia (SCZ) and bipolar disorder (BD), we examined differentially expressed genes in the post-mortem brain [25]. The list of differential genes for psychiatric disorders vs. controls was obtained from the authors’ article.

To annotate the biological functions of the hub modules, we performed enrichment analysis using the GO and KEGG database. Items under FDR q < 0.05 were considered for inclusion in the interpretation and analysis.

Statistical analysis

The data were analysed with SPSS version 26.0 software (SPSS Inc., Chicago, IL, United States). Continuous variables of demographic and clinical indicators were compared between CD patients and healthy control individuals using 2-sample T-tests. Categorical variables were compared using the Pearson c2 test. The Mann-Whitney U test was used if the data did not pass the normality test (Shapiro-Wilk test, p < 0.05). The R package WGCNA (v1.71) and stats (4.1.3) were used for statistical analysis in the transcriptome. When performing multiple comparisons, FDR (Benjamini & Hochberg method, BH) was used to adjust the p-value. Unless explicitly noted, for all analyses, statistical tests were 2-tailed, with a 0.05 alpha level to define statistical significance.


Group differences in demographics and clinical characteristics

Twenty-five CD patients and 18 healthy control participants were finally included in our research. The demographics and clinical characteristics of CD patients and healthy control participants are shown in Table 1. There was no significant difference in sex composition, age, education level, and RNA integrity number. Compared to the HC group, the CD patients showed significant anxiety, and depression based on their higher scores of SDS and SAS in neuropsychological assessment. Serum cortisol, plasma ACTH, and 24-h UFC levels were elevated in patients with CD, and normal cortisol circadian rhythm was disturbed, corresponding with the endocrinological features of Cushing’s syndrome. Serum cortisol and plasma ACTH at 8:00 in CD patients were significantly higher than in the HC group.

Table 1. Demographics and clinical characteristics of the Cushing’s Disease (CD) and healthy controls (HC) groups

CD patients (n = 25)

Healthy controls (n = 18)



Age [years]

41.88 ± 11.05

41.44 ± 9.44



Sex [male/female]





Education [years]

12.68 ± 4.78

10.89 ±4.09



Duration of illness [months]

60.56 ± 70.08




Neuropsychological tests


47.84 ± 11.75

27.78 ± 3.70




49.16 ± 13.66

28.56 ± 4.00



Serum cortisol [nmol/L]


570.80 ± 221.84





726.90 ± 247.59

373.46 ± 107.07




705.48 ± 248.02




ACTH [pmol/L]


13.91 ± 6.77





19.30 ± 12.83

5.70 ± 3.30




21.10 ± 12.77




24-h UFC [nmol/24 h]*

2140.88 ± 1038.66

222.96 ± 100.94



Transcriptome sequencing


8.50 ± 0.55

8.65 ± 0.85



Identification of differentially expressed (DE) modules and modules correlated with neuropsychological traits

Next, using the WGCNA R package, we acquired 12 modules of tightly co-expressed genes from 12,107 genes (Supplementary File Tab. S1). The co-expression modules were summarized by the module eigengene (ME). Module 1-4, 8-9, and 12 show significant ME differences between the CD and HC groups (FDR q < 0.05) and are defined as DE modules (Fig. 1C, 2; Supplementary File Tab. S1). The following analysis considers only the DE module. Further WGCNA analysis was performed in each group to check the preservation of detected modules across groups at the soft power threshold of 26. These analyses indicated that all detected modules had a high level of preservation (Zsummary > 10) (Fig. 3AB; Supplementary File Tab. S4).

Figure 1. Weighted gene co-expression network analysis (WGCNA) results. A. Soft power of all subjects. Left: Network topology analysis for various soft-thresholding powers (weighting coefficient, b). The X-axis represents different soft-thresholding powers, while the y-axis represents the correlation coefficient between log (k) and log [P(k)]. The red line indicates a correlation coefficient of 0.8. Right: The average network connectivity under different weighting coefficients. B. Topological superposition matrix (TOM) dendrogram. Clustering dendrograms of all genes with dissimilarity based on topological overlap. Altogether, 12 coexpression modules were constructed and displayed in different colours. C. Plots of module eigengenes for differential expressed modules [false discovery rate (FDR) q < 0.05]. The box in the boxplots indicates the interquartile range (IQR; Q1 indicates the 25th, while Q3 indicates the 75th percentile), and the whiskers indicate Q1-(1.5*IQR) or Q3+(1.5*IQR). The line within the box represents the median. The y-axis indicates the module eigengene, while the X-axis shows each group and module
Figure 2. Summary of module enrichment. Table indicating differentially expressed (DE) modules, neuropsychological correlation modules, and which gene co-expression modules (rows) were enriched in a variety of different gene lists (columns). The first column is module names. The second column shows the module colour. The third column shows the number of genes in each module. In the first 3 columns, yellow blocks mean passed false discovery rate (FDR) < 0.05 and represents DE modules. In the rest of the columns, modules with enrichments or module-trait analysis passing FDR q < 0.05 for multiple-comparison correction are indicated in yellow. The red box represents module 3, which is not only related to mental status but is enriched in broadly expressed genes

Then we calculate the Pearson correlation coefficient of all module MEs and neuropsychological traits (SDS, SAS) to determine which modules are correlated with these clinical traits (Fig. 2, 4A). The results showed that DE modules 1–4 and 9 were significantly associated with SDS (OR = –0.65 to 0.63, FDR q = 3.26e-05 to 0.04) and SAS (OR = –0.58 to 0.58, FDR q = 2.4e-03 to 0.19) (Supplementary File Tab. S1). These modules might represent the corresponding clinical traits of CD patients, and highly co-expressed genes inside the module have potential biological significance.

Figure 3. Module preservation analysis. Zsummary > 10 (the red line) indicates strong preservation, while 2 < Zsummary < 10 indicates moderate preservation. AB. Zsummary preservation statistics between all subjects and either group CD(A) or group HC(B); C–E. Zsummary preservation statistics between Cushing’s disease (CD) leukocytes and either schizophrenia (SCZ) hippocampus (C), SCZ dorsolateral prefrontal cortex (D), or SCZ striatum (E)

Neuropsychological trait-associated modules are enriched in broadly expressed genes

We next looked for which class of genes contributes greatly to modules correlated with neuropsychological traits. This would allow us to explain why the blood leukocyte transcriptome correlates with neuropsychological changes. Based on the idea of the omnigenic model [19], genes that are broadly expressed in the whole body (including the brain) could also be expressed and measured in blood leukocytes and, therefore, may be highly relevant for modules correlated with neuropsychological traits.

In the blood of patients with Cushing’s disease the neutrophil proportion is usually elevated and the proportion of lymphocytes is reduced. To identify blood-relevant modules and thus exclude the effect of cell proportion, we included whole-blood-specific genes and lymphocyte-specific genes for tissue enrichment (Fig. 2). Among DE modules, module 2 enriches both broadly expressed genes (OR = 2.47, FDR q = 6.63e-35) and blood-specific genes (OR = 7.19, FDR q = 4.01e-33), while modules 1, 4, and 8 were not enriched for any type of genes, module 9 was enriched for blood-specific genes (OR = 5.71, FDR q = 1.63e-8), and module 12 was enriched for lymphocyte-associated genes (OR = 5.02, FDR q = 4.41e-5). The ME difference in blood and lymphocyte-enriched modules was in accordance with the cell proportion changes in Cushing’s disease (Fig. 1C). Although module 2 is also enriched with broadly expressed genes, the influence of blood-specific genes cannot be ruled out. Only module 3 enriched broadly expressed genes alone (OR = 4.85, FDR q = 9.7e-101), which was perhaps the module that really correlated with neuropsychological traits. Detailed enrichment analysis results can be found in Supplementary File Table S1.

Module 3 is enriched in SCZ genes in post-mortem hippocampus and striatum

In our research, CD patients showed higher SDS and SAS scores than the HC group. To run enrichment between module 3 and gene sets of some common psychiatric disorders that may be related to HPA axis dysfunction, SCZ, and BD, we examined genes differentially expressed in the post-mortem brain, which contains several brain regions, including the dorsolateral prefrontal cortex (DLPFC), hippocampus (HIP), and striatum (STR). ME showed that genes in module 3 showed lower expression compared to the control group; therefore, we only included down-regulated genes in the post-mortem brain. We found that module 3 was enriched in the SCZ hippocampus (OR = 3.25, FDR q = 1.39e-30) and striatum (OR = 2.02, FDR q = 2.36E-03) (Fig. 4B; Supplementary File Tab. S2 and S3).

Figure 4. Module-trait and disease enrichment. A. Correlation between modules and clinical traits. Each column corresponds to a trait, and each row corresponds to a module. The number in the rectangle is the correlation coefficient, and the number in brackets is the corresponding false discovery rate (FDR) value. The table is colour-coded by correlation based on the colour legend. SDS_score Self-Rating Depression Scale; SAS_score Self-Rating Anxiety Scale; B. Enrichments between module 3 and dorsolateral prefrontal cortex (DLPFC), hippocampus (HIP), and striatum (STR) of schizophrenia (SCZ), bipolar disorder (BD)

We also tested the preservation of module 3 in the SCZ post-mortem brain. Module 3 shows strong preservation (Zsummary > 10) in all three parts of the SCZ brain (Fig. 3C–E; Supplementary File Tab. S5), though not enriched in DLPFC. Interestingly, the other module (module 7) that showed strong preservation was also enriched in broadly expressed genes alone.

This enrichment between broadly expressed module of Cushing’s disease and SCZ suggests a role of elevated HPA activity in patients with SCZ [26], emphasizing the adverse consequences of chronic stress and its mechanism.

GO and KEGG enrichment analysis

It is meaningful to figure out what biological processes and molecular mechanisms were included in module 3. Broadly expressed genes serve as a bridge that allows us to gain a preliminary understanding of the mechanisms of brain alterations in disease states through peripheral blood. Then we performed GO and KEGG enrichment analysis. Interestingly, in CD patients, the downregulated pathways centred at ribosomal, proteasome and mitochondrial functions, as well as protein synthesis process (Fig. 5 and 6). All enrichment results are shown in Supplementary File Tab. S6–9. In conclusion, the GO and KEGG enrichment analysis gave us a preliminary indication of the key module’s biological function and molecular mechanism.

Figure 5. Kyoto Encyclopaedia of Genes and Genomes (KEGG) enrichment analysis barplot (up) and dotplot (down) show the 10 most enriched pathways
Figure 6. Gene ontology (GO) enrichment analysis. Barplot (left) and dotplot (right) show the 10 most enriched pathways


In this study, we compared blood leukocyte-based gene expression profiles between patients with CD and matched healthy controls and tried to seek biological and molecular mechanisms underlying the brain changes in patients. We performed WGCNA to identify DE gene co-expression modules between groups and modules correlated with neuropsychological traits. Among these modules, tissue class enrichment found module 3, a broadly expressed gene enrichment module, which is consistent with the “omnigenic model” idea. This could be a link between peripheral blood and the brain and might be related to neuropsychological traits. A series of disease enrichment analyses showed that module 3 was associated with SCZ hippocampus and striatum, suggesting the role of HPA axis disturbances in some psychiatric disorders. GO and KEGG enrichment analysis was performed to explore these genes’ biological functions and molecular mechanisms.

The omnigenic model considers that most heritability signals for complex traits or diseases are widely distributed across large genome regions. Although the numerous widespread “peripheral” genes only have a small effect, taken together they have larger effect than a smaller set of “core” genes and probably interact within gene regulatory networks with “core” genes. Such associations may be detectable in many tissue types other than the brain, such as blood leukocytes. The omnigenic model suggests that what is associated with a complex trait are likely to be broadly expressed genes that account for a large percentage of whole genes. Here, we found that module 3 was enriched in broadly expressed genes only, and it can be enriched in the SCZ hippocampus and striatum, showing a strong module preservation in the whole brain of SCZ. These results indicate that the patient transcriptome of blood leukocytes contains information about mental illness and highlight the value of broadly expressed genes in the search for the mechanism of brain changes.

To further explore the biological functions represented by the broadly expressed genes in this study, we performed GO and KEGG enrichment analysis of module 3. Module 3 comprises both the DE module and broadly expressed gene-related module, as well as the neuropsychological correlation module and disease enrichment-associated module. Module 3 is enriched for many pathways that are differentially expressed in psychiatric disorders. Proteasome dysfunction may be associated with defective clearance of misfolded cytoplasmic, nuclear, and ER proteins, leading to intracellular protein aggregation, cytotoxicity, and cell death [27, 28]. Under chronic, repetitive, or severe stress, the proteasome could lose its effectiveness and responsiveness [29, 30]. Functional decline of proteasome has been detected both in the peripheral blood and cerebral cortex in psychiatric disorders SCZ, BD, and anxiety [31–33]. Failure to efficiently remove potentially toxic proteins may lead to increased aggregated insoluble proteins mediated by oxidative stress in neurodegenerative diseases [34], with impairment of the proteasome as the central mechanism of neurodegeneration. Numerous studies have shown that the proteasome proteolytic pathway is also involved in multiple aspects of protein synthesis [35]. Proteasome inhibition is sufficient to impair neuron protein synthesis [36]. This is consistent with our enrichment of other pathways, with GO and KEGG showing a wide range of decreased protein synthesis functions in CD patients. This is also consistent with the reduced grey matter volume seen in CD patients [37, 38]. Mitochondria play a crucial role in many fundamental processes leading to cell death or in synaptic plasticity, such as intracellular calcium buffering, energy production, apoptosis, neurotransmitter transmission, and ROS production. Mitochondria are abundantly present in central nervous cells and not only provide essential energy to neurons but also influence synaptic plasticity. Injury to these organelles may impair cognition through impaired neurotransmission and altered calcium ion (Ca2+) homeostasis [39]. These processes potentially alter neuronal activity, thereby generating aberrant neuronal circuits and plasticity that ultimately affect behavioural outcomes [40]. Moreover, they greatly contribute to the process of neural apoptosis, which is the basis of many mental disorders, such as schizophrenia [41], Alzheimer’s dementia, BD [42], and the most severe forms of major depression [43]. Several important pathways related to mitochondrial function, such as oxidative phosphorylation, mitochondrial electron transport, and mitochondrial membrane, were shown to be significantly downregulated in CD patients. Impairment of mitochondrial oxidative phosphorylation causes cells to switch to anaerobic respiration and produce lactate [42]. In summary, module 3 was enriched not only for pathways related to neurological impairment but also for many pathways that are widely dysregulated in psychiatric disorders. This exemplifies the value of CD as a natural hypercortisolism model in the study of mental illness, and how cortisol affects the brain and subsequently contributes to several psychiatric disorders.

These findings and methods could be used in other research. Blood sampling with high-throughput techniques to quantify the leukocyte transcriptome is feasible for patients with different types and severity of psychiatric disorders. It is also essential to associate in vivo examination of the molecular mechanisms with higher-level macroscale neural systems and heterogeneity in clinical phenotypes for furthering progress toward precision medicine [44]. Considering the inability to directly and noninvasively detect gene expression in brain tissue of living patients, the current approach provides a new in vivo window for understanding molecular mechanisms.

Several limitations in our study should be noted. First, the sample size was relatively small, which may have reduced the statistical power of comparing gene expression and WGCNA between the CD patients and healthy control groups. More subjects should be included in future studies. For enrichment analysis, we included a series of gene lists obtained from research studying psychiatric disease. This contains only a limited number of studies, and exhaustive enrichment analysis for more types of disease should be systematically performed in the future. Furthermore, the blood transcriptome can only indirectly make limited inferences about brain changes, and more in-depth studies of the Cushing’s disease brain are needed if samples are available.


We compared blood leukocyte-based gene expression profiles between patients with CD and matched healthy controls. The co-expression network and enrichment analysis indicated that module 3 of blood leukocytes was enriched in broadly expressed genes and was associated with neuropsychological phenotypes and mental disease enrichment. GO and KEGG enrichment analysis of module 3 identified enrichment in many biological pathways associated with psychiatric disorders, initially elucidating the molecular mechanism of cortisol affecting the brain and causing some mental diseases.

Conflict of interests

The authors declare that they have no competing interests.

Author contributions

K.Y.H., Y.Y.Z., and X.G.Y. contributed to the conception and design of the study. T.Z. collected the sample and organized the database. K.Y.H. performed the R and statistical analysis. K.Y.H. wrote the first draft of the manuscript. F.Y.W., F.Y.L., and Y.Y.Z. wrote sections of the manuscript. All authors contributed to manuscript revision, and read and approved the submitted version.


This study was supported by the National Natural Science Foundation of China (No. 82001798 and No. 81871087). The funder played no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.


We thank Wuhan Metware Biotechnology Co., Ltd for assisting in sequencing.

Data and code availability

The datasets presented in this article can be found at https://ngdc.cncb.ac.cn/search/?dbId=hra&q=HRA003787+. Any request to access the datasets for non-commercial use will be permitted. However, we can only give permission to researchers in mainland China now because the data is being policy reviewed according to the law before publicly providing data. All codes to generate the results in this article can be found at https://github.com/mingchengzu/4CD; part of the code was modified from Dr. Lombardo MV’s great work in autism spectrum disorder analysis [45].


  1. Holtzman CW, Trotman HD, Goulding SM, et al. Stress and neurodevelopmental processes in the emergence of psychosis. Neuroscience. 2013; 249: 172–191, doi: 10.1016/j.neuroscience.2012.12.017, indexed in Pubmed: 23298853.
  2. Agnew-Blais J, Danese A. Childhood maltreatment and unfavourable clinical outcomes in bipolar disorder: a systematic review and meta-analysis. Lancet Psychiatry. 2016; 3(4): 342–349, doi: 10.1016/S2215-0366(15)00544-1, indexed in Pubmed: 26873185.
  3. Moreno-Peral P, Conejo-Cerón S, Motrico E, et al. Risk factors for the onset of panic and generalised anxiety disorders in the general adult population: a systematic review of cohort studies. J Affect Disord. 2014; 168: 337–348, doi: 10.1016/j.jad.2014.06.021, indexed in Pubmed: 25089514.
  4. Zorn JV, Schür RR, Boks MP, et al. Cortisol stress reactivity across psychiatric disorders: A systematic review and meta-analysis. Psychoneuroendocrinology. 2017; 77: 25–36, doi: 10.1016/j.psyneuen.2016.11.036, indexed in Pubmed: 28012291.
  5. Heim C, Newport DJ, Mletzko T, et al. The link between childhood trauma and depression: insights from HPA axis studies in humans. Psychoneuroendocrinology. 2008; 33(6): 693–710, doi: 10.1016/j.psyneuen.2008.03.008, indexed in Pubmed: 18602762.
  6. León-Carrión J, Atutxa AM, Mangas MA, et al. A clinical profile of memory impairment in humans due to endogenous glucocorticoid excess. Clin Endocrinol (Oxf). 2009; 70(2): 192–200, doi: 10.1111/j.1365-2265.2008.03355.x, indexed in Pubmed: 18702680.
  7. Piasecka M, Papakokkinou E, Valassi E, et al. Psychiatric and neurocognitive consequences of endogenous hypercortisolism. J Intern Med. 2020; 288(2): 168–182, doi: 10.1111/joim.13056, indexed in Pubmed: 32181937.
  8. Valli I, Crossley NA, Day F, et al. HPA-axis function and grey matter volume reductions: imaging the diathesis-stress model in individuals at ultra-high risk of psychosis. Transl Psychiatry. 2016; 6(5): e797, doi: 10.1038/tp.2016.68, indexed in Pubmed: 27138796.
  9. Thai M, Schreiner MW, Mueller BA, et al. Coordination between frontolimbic resting state connectivity and hypothalamic-pituitary-adrenal axis functioning in adolescents with and without depression. Psychoneuroendocrinology. 2021; 125: 105123, doi: 10.1016/j.psyneuen.2020.105123, indexed in Pubmed: 33465581.
  10. Madsen KS, Jernigan TL, Iversen P, et al. Hypothalamic-pituitary-adrenal axis tonus is associated with hippocampal microstructural asymmetry. Neuroimage. 2012; 63(1): 95–103, doi: 10.1016/j.neuroimage.2012.06.071, indexed in Pubmed: 22776454.
  11. Belleau EL, Treadway MT, Pizzagalli DA. The Impact of Stress and Major Depressive Disorder on Hippocampal and Medial Prefrontal Cortex Morphology. Biol Psychiatry. 2019; 85(6): 443–453, doi: 10.1016/j.biopsych.2018.09.031, indexed in Pubmed: 30470559.
  12. Mamdani F, Weber MD, Bunney B, et al. Identification of potential blood biomarkers associated with suicide in major depressive disorder. Transl Psychiatry. 2022; 12(1): 159, doi: 10.1038/s41398-022-01918-w, indexed in Pubmed: 35422091.
  13. van de Leemput J, Glatt SJ, Tsuang MT. The potential of genetic and gene expression analysis in the diagnosis of neuropsychiatric disorders. Expert Rev Mol Diagn. 2016; 16(6): 677–695, doi: 10.1586/14737159.2016.1171714, indexed in Pubmed: 27017833.
  14. Rollins B, Martin MV, Morgan L, et al. Analysis of whole genome biomarker expression in blood and brain. Am J Med Genet B Neuropsychiatr Genet. 2010; 153B(4): 919–936, doi: 10.1002/ajmg.b.31062, indexed in Pubmed: 20127885.
  15. Witt SH, Sommer WH, Hansson AC, et al. Comparison of gene expression profiles in the blood, hippocampus and prefrontal cortex of rats. In Silico Pharmacol. 2013; 1: 15, doi: 10.1186/2193-9616-1-15, indexed in Pubmed: 25505659.
  16. Girgenti MJ, Duman RS. Transcriptome Alterations in Posttraumatic Stress Disorder. Biol Psychiatry. 2018; 83(10): 840–848, doi: 10.1016/j.biopsych.2017.09.023, indexed in Pubmed: 29128043.
  17. Cattane N, Minelli A, Milanesi E, et al. Altered gene expression in schizophrenia: findings from transcriptional signatures in fibroblasts and blood. PLoS One. 2015; 10(2): e0116686, doi: 10.1371/journal.pone.0116686, indexed in Pubmed: 25658856.
  18. Guidotti A, Auta J, Davis JM, et al. Toward the identification of peripheral epigenetic biomarkers of schizophrenia. J Neurogenet. 2014; 28(1-2): 41–52, doi: 10.3109/01677063.2014.892485, indexed in Pubmed: 24702539.
  19. Boyle EA, Li YI, Pritchard JK. An Expanded View of Complex Traits: From Polygenic to Omnigenic. Cell. 2017; 169(7): 1177–1186, doi: 10.1016/j.cell.2017.05.038, indexed in Pubmed: 28622505.
  20. Lacroix A, Feelders RA, Stratakis CA, et al. Cushing’s syndrome. Lancet. 2015; 386(9996): 913–927, doi: 10.1016/S0140-6736(14)61375-1, indexed in Pubmed: 26004339.
  21. van der Werff SJA, Andela CD, Nienke Pannekoek J, et al. Widespread reductions of white matter integrity in patients with long-term remission of Cushing’s disease. Neuroimage Clin. 2014; 4: 659–667, doi: 10.1016/j.nicl.2014.01.017, indexed in Pubmed: 24936417.
  22. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008; 9: 559, doi: 10.1186/1471-2105-9-559, indexed in Pubmed: 19114008.
  23. Nieman LK, Biller BMK, Findling JW, et al. Endocrine Society. Treatment of Cushing’s Syndrome: An Endocrine Society Clinical Practice Guideline. J Clin Endocrinol Metab. 2015; 100(8): 2807–2831, doi: 10.1210/jc.2015-1818, indexed in Pubmed: 26222757.
  24. Langfelder P, Luo R, Oldham MC, et al. Is my network module preserved and reproducible? PLoS Comput Biol. 2011; 7(1): e1001057, doi: 10.1371/journal.pcbi.1001057, indexed in Pubmed: 21283776.
  25. Lanz TA, Reinhart V, Sheehan MJ, et al. Postmortem transcriptional profiling reveals widespread increase in inflammation in schizophrenia: a comparison of prefrontal cortex, striatum, and hippocampus among matched tetrads of controls with subjects diagnosed with schizophrenia, bipolar or major depressive disorder. Transl Psychiatry. 2019; 9(1): 151, doi: 10.1038/s41398-019-0492-8, indexed in Pubmed: 31123247.
  26. Jacobson L. Hypothalamic-pituitary-adrenocortical axis: neuropsychiatric aspects. Compr Physiol. 2014; 4(2): 715–738, doi: 10.1002/cphy.c130036, indexed in Pubmed: 24715565.
  27. Kravtsova-Ivantsiv Y, Cohen S, Ciechanover A. Modification by single ubiquitin moieties rather than polyubiquitination is sufficient for proteasomal processing of the p105 NF-κB precursor. Adv Exp Med Biol. 2011; 691: 95–106, doi: 10.1007/978-1-4419-6612-4_10, indexed in Pubmed: 21153313.
  28. Demuro A, Mina E, Kayed R, et al. Calcium dysregulation and membrane disruption as a ubiquitous neurotoxic mechanism of soluble amyloid oligomers. J Biol Chem. 2005; 280(17): 17294–17300, doi: 10.1074/jbc.M500997200, indexed in Pubmed: 15722360.
  29. Raynes R, Juarez C, Pomatto LCD, et al. Aging and SKN-1-dependent Loss of 20S Proteasome Adaptation to Oxidative Stress in C. elegans. J Gerontol A Biol Sci Med Sci. 2017; 72(2): 143–151, doi: 10.1093/gerona/glw093, indexed in Pubmed: 27341854.
  30. Pickering AM, Vojtovich L, Tower J, et al. Oxidative stress adaptation with acute, chronic, and repeated stress. Free Radic Biol Med. 2013; 55: 109–118, doi: 10.1016/j.freeradbiomed.2012.11.001, indexed in Pubmed: 23142766.
  31. Bousman CA, Chana G, Glatt SJ, et al. Preliminary evidence of ubiquitin proteasome system dysregulation in schizophrenia and bipolar disorder: convergent pathway analysis findings from two independent samples. Am J Med Genet B Neuropsychiatr Genet. 2010; 153B(2): 494–502, doi: 10.1002/ajmg.b.31006, indexed in Pubmed: 19582768.
  32. Fedoce Ad, Ferreira F, Bota RG, et al. The role of oxidative stress in anxiety disorder: cause or consequence? Free Radic Res. 2018; 52(7): 737–750, doi: 10.1080/10715762.2018.1475733, indexed in Pubmed: 29742940.
  33. Arion D, Corradi JP, Tang S, et al. Distinctive transcriptome alterations of prefrontal pyramidal neurons in schizophrenia and schizoaffective disorder. Mol Psychiatry. 2015; 20(11): 1397–1405, doi: 10.1038/mp.2014.171, indexed in Pubmed: 25560755.
  34. Grune T, Jung T, Merker K, et al. Decreased proteolysis caused by protein aggregates, inclusion bodies, plaques, lipofuscin, ceroid, and ‘aggresomes’ during oxidative stress, aging, and disease. Int J Biochem Cell Biol. 2004; 36(12): 2519–2530, doi: 10.1016/j.biocel.2004.04.020, indexed in Pubmed: 15325589.
  35. Ding Q, Cecarini V, Keller JN. Interplay between protein synthesis and degradation in the CNS: physiological and pathological implications. Trends Neurosci. 2007; 30(1): 31–36, doi: 10.1016/j.tins.2006.11.003, indexed in Pubmed: 17126920.
  36. Ding Q, Dimayuga E, Markesbery WR, et al. Proteasome inhibition induces reversible impairments in protein synthesis. FASEB J. 2006; 20(8): 1055–1063, doi: 10.1096/fj.05-5495com, indexed in Pubmed: 16770004.
  37. Andela CD, van der Werff SJA, Pannekoek JN, et al. Smaller grey matter volumes in the anterior cingulate cortex and greater cerebellar volumes in patients with long-term remission of Cushing’s disease: a case-control study. Eur J Endocrinol. 2013; 169(6): 811–819, doi: 10.1530/EJE-13-0471, indexed in Pubmed: 24031092.
  38. Santos A, Granell E, Gómez-Ansón B, et al. Depression and Anxiety Scores Are Associated with Amygdala Volume in Cushing’s Syndrome: Preliminary Study. Biomed Res Int. 2017; 2017: 2061935, doi: 10.1155/2017/2061935, indexed in Pubmed: 28607927.
  39. Petschner P, Gonda X, Baksa D, et al. Genes Linking Mitochondrial Function, Cognitive Impairment and Depression are Associated with Endophenotypes Serving Precision Medicine. Neuroscience. 2018; 370: 207–217, doi: 10.1016/j.neuroscience.2017.09.049, indexed in Pubmed: 28987512.
  40. Ben-Shachar D. The interplay between mitochondrial complex I, dopamine and Sp1 in schizophrenia. J Neural Transm (Vienna). 2009; 116(11): 1383–1396, doi: 10.1007/s00702-009-0319-5, indexed in Pubmed: 19784753.
  41. Gonçalves VF, Giamberardino SN, Crowley JJ, et al. Examining the role of common and rare mitochondrial variants in schizophrenia. PLoS One. 2018; 13(1): e0191153, doi: 10.1371/journal.pone.0191153, indexed in Pubmed: 29370225.
  42. Jeong H, Dimick MK, Sultan A, et al. Peripheral biomarkers of mitochondrial dysfunction in adolescents with bipolar disorder. J Psychiatr Res. 2020; 123: 187–193, doi: 10.1016/j.jpsychires.2020.02.009, indexed in Pubmed: 32078836.
  43. Marazziti D, Baroni S, Picchetti M, et al. Psychiatric disorders and mitochondrial dysfunctions. Eur Rev Med Pharmacol Sci. 2012; 16(2): 270–275, indexed in Pubmed: 22428481.
  44. Courchesne E, Pramparo T, Gazestani VH, et al. The ASD Living Biology: from cell proliferation to clinical phenotype. Mol Psychiatry. 2019; 24(1): 88–107, doi: 10.1038/s41380-018-0056-y, indexed in Pubmed: 29934544.
  45. Lombardo MV, Eyler L, Pramparo T, et al. Atypical genomic cortical patterning in autism with poor early language outcome. Sci Adv. 2021; 7(36): eabh1663, doi: 10.1126/sciadv.abh1663, indexed in Pubmed: 34516910.