Introduction
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
Subjects
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.
Results
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.
|
|
CD patients (n = 25) |
Healthy controls (n = 18) |
t/Z/c2 |
p-value |
|
Age [years] |
41.88 ± 11.05 |
41.44 ± 9.44 |
0.135 |
0.89 |
|
Sex [male/female] |
3/22 |
2/16 |
N.A. |
1.000 |
|
Education [years] |
12.68 ± 4.78 |
10.89 ±4.09 |
–1.413 |
0.158 |
|
Duration of illness [months] |
60.56 ± 70.08 |
N.A. |
N.A. |
N.A. |
|
Neuropsychological tests |
||||
|
SDS |
47.84 ± 11.75 |
27.78 ± 3.70 |
8.001 |
0.000 |
|
SAS |
49.16 ± 13.66 |
28.56 ± 4.00 |
7.129 |
0.000 |
|
Serum cortisol [nmol/L] |
||||
|
00:00 |
570.80 ± 221.84 |
N.A. |
N.A. |
N.A. |
|
08:00 |
726.90 ± 247.59 |
373.46 ± 107.07 |
6.359 |
0.000 |
|
16:00 |
705.48 ± 248.02 |
N.A. |
N.A. |
N.A. |
|
ACTH [pmol/L] |
||||
|
0:00 |
13.91 ± 6.77 |
N.A. |
N.A. |
N.A. |
|
08:00* |
19.30 ± 12.83 |
5.70 ± 3.30 |
–5.133 |
0.000 |
|
16:00 |
21.10 ± 12.77 |
N.A. |
N.A. |
N.A. |
|
24-h UFC [nmol/24 h]* |
2140.88 ± 1038.66 |
222.96 ± 100.94 |
–5.541 |
0.000 |
|
Transcriptome sequencing |
||||
|
RIN |
8.50 ± 0.55 |
8.65 ± 0.85 |
–0.686 |
0.497 |
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).
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.
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).
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.
Discussion
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.
Conclusion
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.
Funding
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.
Acknowledgments
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].
