Brain differential gene expression and blood cross-validation of a molecular signature of patients with major depressive disorder.

INTRODUCTION: The agreement between clinicians diagnosing major depressive disorder (MDD) is poor. The objective of this study was to identify a reproducible and robust gene expression marker capable of differentiating MDD from healthy control (HC) subjects.
MATERIALS AND METHODS: Brain and blood gene expression datasets were searched, which included subjects with MDD and HC. The largest database including different areas of brain samples (GSE80655) was used to identify an initial gene expression marker. Tests of robustness and reproducibility were then implemented in 13 brain and 7 blood independent datasets. Correlations between expression in brain and blood samples were also examined. Finally, an enrichment analysis to explore the marker biological meaning was completed.
RESULTS: Twenty-eight genes were differentially expressed in GSE80655, of which 23 were critical to differentiate MDD from HC. The accuracy obtained using the 23 genes was 0.77 and 0.8, before and after the forward selection model, respectively. The gene marker's robustness and reproducibility were between the range of 0.46 and 0.63 in the other brain datasets and between 0.45 and 0.78 for the blood datasets. Brain and blood expression tended to correlate in some samples. Thirteen of the 23 genes were related to stress and immune response.
CONCLUSION: A 23 gene expression marker was able to distinguish subjects with MDD from HC, with adequate reproducibility and low robustness in the independent databases investigated. This gene set was similarly expressed in the brain and blood and involved genes related to stress and immune response.

Introduction

Major depressive disorder (MDD) is highly prevalent (‘NIMH » Major Depression’, n.d.) and a leading cause of disability, affecting about 4.4% of the world population in the year 2017 (Friedrich, 2017). However, the diagnostic reliability of MDD has been questioned. According to the field trials of the DSM-5 Mood Disorders Working Group, the agreement between clinicians diagnosing MDD only achieved an intraclass kappa of 0.28 [95% confidence interval (CI): 0.20–0.35] (Regier ). This limited agreement has been also observed in other studies (Enns ; Canuto ; Behera ).

As in other psychiatric disorders, the development of biomarkers has been pursued to improve the diagnostic reliability of MDD. Hence, many brain and peripheral features have been examined such as growth factors, cytokines, inflammatory markers, oxidative stress markers, endocrine markers, energy balance hormones, structural and functional imaging markers (Hacimusalar and Eşel, 2018). Though several biomarker studies for MDD have been performed (Spijker ; Savitz ; Sullivan ; Chang ; Liu ; Wingo ; Pantazatos ; Hacimusalar and Eşel, 2018), the results thus far have shown limitations in reproducibility and robustness (Savitz ; Sullivan ; Chang ; Pantazatos ). Robustness refers to the persistence of an effect in settings different from and outside of an experimental framework. Results’ reproducibility refers to obtaining the same results from the conduct of an independent study that uses the same procedures (Goodman ). There are no currently available biomarkers for MDD or any psychiatric disorder to assist the clinician diagnostically (Hacimusalar and Eşel, 2018). Because the heritability of MDD is 37% (95% CI: 31–42%) (Sullivan ), genomic studies have been at the forefront of biomarker development for depression (Gadad ). Many psychiatric genomic studies have used peripheral blood mononuclear cell (PBMC) as a practical source of genetic material. Because of the immune disturbances reported in MDD (Bahn ; Li ; Tomita ; Franz ; Strawbridge ; Gao ; Bhak ), PBMCs offer other potential advantages as a tissue source for biomarkers. However, the impermeability of the blood–brain barrier prevents the PBMCs access to the brain in normal conditions (Erdő ), questioning their relevance to disorders of the brain, like MDD. Still, it has been reported that gene expression from blood and brain are somewhat correlated in several medical conditions (Canuto ; Arion ; Ramaker ), including MDD (Tolentino and Schmidt, 2018). Furthermore, concordance between gene expression variation in postmortem brain and in PBMCs would validate any potential peripheral biomarkers (Deep-Soboslay ; Harrison, 2011; McCullumsmith and Meador-Woodruff, 2011). However, for obvious practical reasons, no biomarker studies of gene expression have jointly examined postmortem brain and PBMC in MDD. Still, the complementary examination of gene expression in brain and PBMCs from different cohorts of subjects with MDD could advance the development of genomic biomarkers.

As previously stated, there is a low concordance in the diagnosis of MDD between clinicians. A biomarker could improve this concordance and the accuracy in the diagnosis of MDD. Most of the previous gene expression studies in MDD lack tests of reproducibility and robustness, and are limited to the population where the biomarker was identified. In the present study, with access to multiple datasets, we performed extensive robustness and reproducibility tests. In an effort to identify a reproducible and robust gene expression biomarker capable of differentiating MDD from healthy control (HC) subjects, we searched all the relevant online databases available up to October 2019. We started by examining protein-coding genes expressed in brain that could differentiate MDD from HC. We then tested the expression of this gene set in PBMCs in several other datasets. Subsequently, we examined the correlation of brain and blood samples of subjects with MDD and HC participants using the gene expression marker. Finally, we performed an enrichment analysis to determine the biological implications of the set of genes.

Methods

Identification of brain and blood datasets on Gene Expression Omnibus DataSets of NCBI

The workflow of the present study is summarized in Fig. 1. The search of databases was completed employing the terms used on Gene Expression Omnibus (GEO) DataSets (https://www.ncbi.nlm.nih.gov/gds) of NCBI: ‘MDD’ AND ‘Gene expression’, ‘MDD’ AND ‘RNA expression’, ‘Major Depressive Disorder’ AND ‘Gene expression’, ‘Major Depressive Disorder’ AND ‘RNA expression’, ‘Depression’ AND ‘Gene expression’, and ‘Depression’ AND ‘RNA expression’. The robustness and reproducibility tests were completed in datasets of brain and blood, also identified with the previously mentioned terms.

Fig. 1

Summarization of material and methods.

Analysis and filtering of variables

After the selection of appropriate datasets was completed, a visual inspection of the data was performed, to confidently ascertain datasets with no missing gene expression data, meaning no empty cells in the data matrix. The platforms used to measure gene expression were microarrays or RNASeq (Table 1). The microarray datasets were normalized using quantile normalization, to control as much as possible the technical noise of the microarrays. Moreover, all the datasets included in this study were scale rank normalized with a distribution from 0 to 1, a process denominated uniformization, to address the platform differences between RNASeq and microarrays. Dataset GSE80655 was normalized with the DESeq2 R package (Love ), with the function varianceStabilizationTransformation.

Table 1

Classification accuracy in datasets

Dataset ID	Platform/tissue	Total no. of probes/total no. of genes of the marker present	Subjects per class	Accuracy
Brain
GSE80655	RNASeq/DLPFC, nAcc, anCg	23/23	HC 70, MDD 69	0.8
GSE87610	Micro/DLPFC	48/22	HC 15, MDD 19	0.63
GSE101521	RNASeq/DLPFC	23/23	HC 28, MDD 31	0.62
GSE54564	Micro/Amyg	31/20	HC 21, MDD 21	0.59
GSE44593	Micro/Amyg	27/19	HC 14, MDD 14	0.58
GSE54570	Micro/DLPFC	19/15	HC 13, MDD 13	0.58
GSE54562 + 54563	Micro/anCg	31/20	HC 35, MDD 35	0.55
GSE54575	Micro/orbitoPFC	19/15	HC 12, MDD 12	0.51
GSE53987	Micro/Hip, Stra, PFC	31/18	HC 55, MDD 50	0.49
GSE54565 + 54572	Micro/anCg	31/20	HC 28, MDD 28	0.47
GSE54567 + 54568	Micro/DLPFC	30/20	HC 29, MDD 29	0.46
Blood PBMC
GSE38206	Micro/PBMC	38/19	HC 9, MDD 9	0.78
GSE76826	Micro/PBMC	48/21	HC 12, MDD 10	0.64
GSE39653	Micro/PBMC	29/19	HC 24, MDD 21	0.62
GSE98793	Micro/PBMC	31/18	HC 64, MDD 128	0.59
GSE19738	Micro/PBMC	25/19	HC 34, MDD 33	0.51
GSE67663	Micro/PBMC	12/9	HC 72, MDD 112	0.51
GSE52790	Micro/PBMC	20/19	HC 12, MDD 10	0.45

Amyg, amygdala; anCg, anterior cingulate; DLPFC, dorsolateral prefrontal cortex; FSM, forward selection model; HC, healthy control; Hip, hippocampus; MDD, patient with major depressive disorder; Micro, microarrays; nAcc, nucleus accumbens; OrbitoPFC, orbito prefrontal cortex; PBMC, peripheral blood mononuclear cell; PFC, prefrontal cortex; RNASeq, RNA expression measured by sequencing of RNA; Stra, striatum.

Due to its larger sample, greater number of brain structures contained in the dataset, and use of RNASeq data, database GSE80655 was selected to first examine a molecular signature based on gene expression. This dataset was filtered for differential gene expression between patients with MDD and HCs using the function DESeq of the DESeq2 R Package. This procedure was completed to include in the analysis only genes with different expressions in both classes, that is, a gene with high expression in HC and low expression in MDD or vice versa. This means, the classification between MDD and HC used the combined differential weighted gene expression level as a continuous variable; no cutoff of gene expression group difference was used. The final combination of genes was selected with a classification algorithm designed to improve accuracy, which is explained in more detail in the next section. The protein-coding genes with an adjusted P value in the differential gene expression less than 0.05, were selected to identify the gene expression marker. The P value of the differential gene expression was calculated with a false discovery rate method (Benjamini and Hochberg, 1995) only in protein-coding genes, to avoid type I errors in null hypothesis testing when conducting multiple comparisons. Finally, the decision of working only with protein-coding genes was determined by the lack of additional datasets with noncoding RNA to test the performance of the biomarker.

GALGO R package for marker identification and feature selection

To identify the gene expression marker in the dataset GSE80655, the Genetic ALGOrithms (GALGO) R package was used to classify the patients with MDD and HC, with the classification approach of ‘nearest centroid’ (Trevino and Falciani, 2006). GALGO is a genetic algorithm based on evolutionary-like principles. It combines genes, in this case, their expression, to find the best combination to classify any given classes (HC and MDD), with the highest possible accuracy. A chromosome in GALGO is defined as a combination of genes (biomarker size); a generation is the process of internal combinations of chromosomes, which is controlled by the user with a maximum number of solutions and a goal of fitness (accuracy). The user can control the desired fitness, and the number of solutions if that fitness is not reached during the specified generation. If the accuracy (fitness) is not reached, and the max solutions target is reached, GALGO moves on to the next generation, and the process continues until it reaches the generations defined by the user, in our case 300.

An internal validation of two-thirds of the sample for training and one-third for testing was done to avoid overfitting. The accuracy of the gene expression marker was used to measure the statistical performance to classify between the two groups. The accuracy is the average between sensitivity and specificity, which is defined as Accuracy = (TP + TN)/(TP + TN + FP + FN), where TP is True Positive, TN is True Negative, FP is False Positive and FN is False Negative.

The configuration of GALGO used chromosome size 10 (10-gene marker size), max solutions 300 and goal fitness 0.9. Also, to estimate the final error, the bootstrap method was used with 0.632 weight for test and 0.368 for train (Bradley and Tibshirani, 1993). The training error was a 10-fold cross-validation and the first-level splitting selected was 0.5 and 0.5 in only the 20 first splits. We used a forward selection model (FSM) algorithm only in dataset GSE80655, because this was potentially the most informative one because of its larger sample, greater number of brain structures included and use of RNASeq data. FSM is a method of choosing the expression of the genes based on their contribution in improving the classification accuracy. The process starts with the gene that has the highest accuracy, meaning a univariate classification between MDD and HC. Then, each next gene that significantly improves accuracy is added until there is no further improvement.

Finally, a heatmap of the gene expression signature was generated to visually display the difference in the gene expression between MDD and HC subjects, using gplots R Package, with the data z-scaled per row.

Robustness and reproducibility capacity tests

The robustness and reproducibility tests were implemented in brain and blood datasets. Reproducibility (also referred to as replicability) refers to obtaining the same results from the conduct of an independent study whose procedures are as closely matched to the original experiment as possible. Robustness refers to the stability of experimental conclusions to variations in either baseline assumptions or experimental procedures. It is somewhat related to the concept of generalizability (also known as transportability), which refers to the persistence of an effect in settings different from and outside of an experimental framework (Goodman ).

Due to the lack of a quantitative definition of reproducibility and robustness in the literature, we defined some parameters. Reproducibility was considered high, moderate or low, taking into account the SD calculated from the accuracy result in all datasets. Hence, high reproducibility is an accuracy within 1 SD, moderate within 2 SDs and low within 3 SDs or higher, in comparison with the accuracy of the normative dataset GSE80655 (which was 0.8).

We defined robustness as high when more than 90% of the independent datasets achieve reproducibility within 1 SD in comparison with the accuracy of GSE80655. For moderate and low robustness, we considered that more than 90% of the independent datasets achieved reproducibility within 2 and 3 SDs, respectively, in comparison with the accuracy of GSE80655.

Correlation analysis between samples in brain and blood datasets

To examine a relationship between expression in brain and blood tissue samples using the shared genes of each dataset of the identified marker, an analysis of correlation was calculated in R with the function Pearson’s, creating a correlation matrix. The correlation matrix is presented as a heatmap with the function heatmap.2 of gplots R package.

Enrichment analysis of the identified marker

Finally, to investigate the biological implications and significance of the identified gene expression marker, a Gene Ontology (GO) analysis was performed using GO terms on the STRING online database (https://string-db.org).

STRING determines protein–protein interaction by genomic context (neighborhood, fusion and co-occurrence), genomic experiments (co-expression and experiments), and prior consolidated knowledge on protein–protein associations (knowledge and text-mining). In the case of association by neighborhood, two proteins are given an association score when their encoding genes are in close proximity to each other on the chromosome. For the fusion channel, all genomes are scanned for open reading frames that appear to be the result of gene fusion events. In co-occurrence, pairs of genes are identified whose occurrence patterns throughout evolution show similarities. For co-expression, pairs of genes that show consistent similarities between their expression profiles are assigned association scores. The experiments collect protein–protein interaction evidence from experiments and assays in the laboratory. Knowledge channel parses association evidence from curated pathway databases (i.e. Kyoto Encyclopedia of Genes and Genomes and Reactome), where it has been collected and consolidated manually by expert curators. Finally, the text-mining channel is based on PubMed abstracts and other sources. Pairs of proteins mentioned together in the same sentence, the same paragraph or merely the same publication are assigned a benchmarked association score (Szklarczyk ).

Results

For the current in-silico analysis, 21 datasets were identified on the GEO Datasets of NCBI (Table 1), 14 from brain and 7 from blood tissues. The clinical data can be reviewed in Supplementary Table 1, Supplemental digital content 1, http://links.lww.com/PG/A260.

The datasets GSE54562 and GSE54563 were merged and analyzed as a one dataset. Merging was also done for the datasets GSE54567 with GSE54568, and GSE54565 with GSE54572, resulting in three datasets instead of six. This was done because the datasets shared the same microarrays platform to measure gene expression, as well as the same brain structures (Table 1). This increased the sample size in the three final datasets and decreased the internal variance.

The gene expression marker was identified in the dataset GSE80655, because it had more diverse brain areas, and the RNA levels were measured with RNASeq, which is more accurate than microarrays. From the 57 905 total probes of the dataset, 56 658 probes corresponded genes with protein product, 958 to pseudogenes, 1 to an RNA long non-coding and 288 with unknown biological implications. Twenty-eight genes were differentially expressed between MDD and HC.

After the filtering process was completed, the 28 genes were analyzed with GALGO R package. As a result, 23 genes were selected as the most important to differentiate between patients with MDD and HC, with an accuracy of 0.8 (Table 1).

The results of the accuracy of the tests are presented in Table 1 for each dataset. Figure 2 illustrates the heatmap of the gene expression differences across the samples in the MDD and HC classes.

Fig. 2

Heatmap with a hierarchical clustering of the gene expression of the 23 gene marker in the dataset GSE80655. The classes are MDD (subjects with major depressive disorder) and HC (healthy control), and are grouped by diagnosis in the x axis. The genes, grouped hierarchically in the y axis, have a yellow color if the expression is high in that subject, whereas the genes with low expression have a blue color.

The dataset GSE38206 had gene expression of subjects both during baseline depression and 8 weeks later after remission. In the present analysis, only the baseline sample was included. GSE87610 included postmortem data in schizophrenia, bipolar disorder, MDD and HC. Gene expression was from pyramidal neurons in layers 3 and 5 of the dorsolateral prefrontal cortex (DLPFC). We averaged expression from both layers because they were highly correlated in both MDD and HC. GSE76826 examined gene expression in late-onset depression (LOD) at baseline, the same subjects in remission and HCs. We included only the baseline data of LOD as well as HCs. The dataset GSE19738 originally compared baseline and activation of the PBMCs to lipopolysaccharide in subjects with and without depression. We used gene expression of PBMC before activation. Finally, in the dataset GSE67663, the subjects had comorbidity of depression and posttraumatic stress disorder.

The single SD of the accuracy for all datasets was 0.096; for 2 SDs it was 0.19, and for 3 SDs it was 0.29. Hence, datasets with an accuracy between 0.8 and 0.7 had high reproducibility; those with an accuracy of 0.61–0.69, moderate reproducibility; and those with less than 0.6 were considered to have low reproducibility. The only dataset with high reproducibility was GSE38206. The datasets GSE87610, GSE101521, GSE76826 and GSE39653 had moderate reproducibility, and the rest (12 datasets) had low reproducibility.

Based on our definition of robustness, our marker had a low robustness performance. Only 5 of the 17 datasets had a high or moderate reproducibility (30%). However, datasets GSE54564, GSE44593, GSE54570 and GSE98793 had accuracies between 0.58 and 0.59, very close to the 0.61 threshold for moderate reproducibility. Hence, the inclusion of these four datasets improved the overall robustness to 53% with moderate or high reproducibility.

Several complementary approaches were used to examine the correlation of gene expression in brain and blood tissues. First, we explored the impact of the platform used (microarray or RNASeq). The datasets GSE80655 (DLPFC measured with RNASeq), GSE101521 (DLPFC measured by RNASeq), and GSE87610 (DLPFC measured with microarrays) were selected because they shared more genes from the 23 gene marker. Results are depicted in Fig. 3 showing a high correlation between the genes [correlation mean = 0.7 (range = 0.16–1)].

Fig. 3

Heatmap of the correlation between the subjects MDD and HC in the datasets GSE80655, GSE101521 and GSE87610 of only the DLPFC of the mentioned datasets. This correlation was calculated using the gene expression of 20 of the 23 genes of original the marker, due to these are the shared genes in these three datasets. The purpose of correlating the data of RNASeq (GSE80655 and GSE101521) and microarrays (GSE87610) in the same brain area (DLPFC) was to test if there is a correlation dependent of the used platform to measure the gene expression. The correlations are hierarchically grouped in x and y axes. A positive linear correlation is represented as red color whereas a negative linear correlation is blue; the color white is used to represent the lack of linear correlation. DLPFC, dorsolateral prefrontal cortex; HC, healthy control; MDD, major depressive disorder.

Second, a correlation between brain samples (datasets GSE80655, GSE87610, GSE101521 and GSE53987), and selected PBMC samples (datasets GSE38206, GSE76826, GSE39653, GSE52790 and GSE98793) was done, including 14 shared genes in the original 23 gene marker. The results are depicted in Fig. 4 showing a higher linear correlation with Pearson coefficients between samples in different tissues; specifically, the brain datasets GSE87610 and GSE53987, with the blood datasets GSE39653 and GSE38206 with an average correlation of 0.38 (range = −0.3 to 0.8). The selection of these datasets was based on (a) higher accuracy in comparison with the other datasets (greater than or equal to 0.61; hence, at least moderate reproducibility); (b) higher number of genes present in the dataset corresponding to the marker (greater than or equal to 18 genes); or (c) higher number of samples in the datasets (n ≥ 100; this size facilitates the testing of correlations among samples independent of accuracy). More specifically, GSE80655, GSE87610, GSE38206, GSE76826 and GSE39653 were selected due to their higher accuracy. (i.e. greater than or equal to 0.61). GSE101521 was selected because the 23 genes of the marker were present in the dataset and also measured the gene expression with RNASeq. GSE52790 was included based on the number of genes of the biomarker present in the dataset (19/23). Finally, GSE53987 and GSE98793 were selected because of their particularly large samples, 105 and 182, respectively.

Fig. 4

Heatmap of the correlation using 14 of the 23 genes from the original marker, between the subjects MDD and HC, between the datasets of brain tissue GSE80655, GSE87610, GSE101521 and GSE53987 and the PBMC datasets GSE38206, GSE76826, GSE39653, GSE52790 and GSE98793. The correlations are hierarchically grouped in x and y axes. A positive linear correlation is represented as red color whereas a negative linear correlation is blue; the color white is used to represent the lack of linear correlation. HC, healthy control; MDD, major depressive disorder; PBMC, peripheral blood mononuclear cell.

Finally, there was no statistical difference between brain and blood samples in the accuracy of the gene expression set we found using a t test (t = −0.36, P = 0.77).

The results from the GO analysis are presented in Table 2.

Table 2

Gene Ontology analysis of the gene expression marker

Term description	Observed gene count	Background gene count	False discovery rate	Matching proteins in your network (labels) gene ID
Response to stress	13	3267	0.0013	CCL4, CCL4L1, CPS1, DMBT1, GIG25, HAMP, HSPA6, LCN2, OLR1, PRF1, SERPING1, SGK1, TATDN2
Defense response	9	1234	0.0013	CCL4, CCL4L1, DMBT1, GIG25, HAMP, LCN2, OLR1, PRF1, SERPING1
Immune response	10	1560	0.0013	CCL4, CCL4L1, DMBT1, GIG25, HAMP, HSPA6, LCN2, OLR1, PRF1, SERPING1
Immune effector process	7	927	0.0072	DMBT1, GIG25, HSPA6, LCN2, OLR1, PRF1, SERPING1

Discussion

We examined all the available data up to October 2019 in the GEO Datasets https://www.ncbi.nlm.nih.gov/geo/, from subjects with MDD and HC and from brain and blood tissues. We identified a set of 23 protein-coding genes, capable of distinguishing MDD from HC with an average accuracy greater than 0.50, in 7 of 11 brain and in 6 of 7 blood datasets. Furthermore, the expression of 14 of the 23 genes was correlated in some brain and blood tissue samples. In the enrichment analysis, the identified genes were associated with immune response, a pathophysiological mechanism known to be involved in the development of depression. It is worth highlighting that the majority of the databases employed in this study were not originally used for the distinction of MDD from HCs.

Several studies have used gene expression sets to distinguish MDD from HC groups. Yi ) examined PBMCs with microarrays in MDD (n = 8), HC (n = 8), and subsyndromal symptomatic depression (n = 8), and reported a model of 48 genes with an accuracy of 100%, using a leave-one-out cross-validation (LOOCV) approach. However, LOOCV is known to be inconsistent with the linear model assumptions used. Hence, this study was limited by a very small sample and lack of validation in an independent dataset.

Wang ) originally used dataset GSE98793 (included in our analyses) which they split into two subgroups of HC (n = 32 each) and two of MDD (n = 64 each). They identified a group of 20 hub genes (genes with high protein–protein interactions with other genes) in PBMCs capable of classifying MDD from HC, with an accuracy (from receiving operating curve analyses) between 0.833 and 0.876 with validation in the second subgroups.

Miyata ) studied late-onset MDD (LOD; in subjects aged greater than or equal to 50 years; n = 18) and HC (n = 12). They first examined a 3066 probes microarray in PBMCs from a mice model of depression. Fourteen genes were differentially expressed and concordant between the animal model and LOD subjects. Finally, the ability to classify LOD and HC was tested using reverse transcription polymerase chain reaction (RT-PCR) in independent LOD subjects (n = 23). The CIDEC gene (cell death-inducing DFFA-like effector C) reached the highest accuracy (0.9185). However, their sample was small.

Spijker ) studied whole blood with microarrays from a primary cohort of MDD (n = 21) and HC (n = 21). They reported a set of seven genes differentially expressed. These findings were validated in an independent cohort of MDD (n = 13) and HC (n = 14) using RT-PCR. The reported specificity was 71.4% and sensitivity was 76.9%; accuracy was not reported. As was already mentioned, small studies tend to inflate effects (Button ).

Bhak ) used methylome acquired with methyl-seq data and transcriptome obtained with whole-transcriptome sequencing from PBMCs of subjects with a history of suicidal attempt (SA, n = 56), MDD (n = 39) and HC (n = 87). They first reported differentially methylated sites (DMSs) for each clinical group versus HC. The MDD had 80 DMS and the SA group had 95 compared to HCs. Gene expression models failed to differentiate the clinical groups from HCs. However, gene expression combined with DMS achieved accuracies of 87.3% for MDD versus HC and 86.7% for SA versus HC. This underlines the potential advantage of linking gene expression with epigenetic markers.

Watanabe ) through PCR array plates in a primary dataset assessed 40 candidate genes in PBMCs from subjects with MDD (n = 25) and HC (n = 25). A group of seven genes achieved a sensitivity of 80% and specificity of 92% in the primary dataset. In an independent dataset test, including subjects with MDD (n = 20) and HC (n = 18), the specificity and sensitivity were 85 and 89%, respectively. Accuracy was not reported.

This literature has not succeeded in the identification of a reliable gene expression biomarker for MDD in the blood (Spijker ; Yi ; Watanabe ; Miyata ; Bhak ; Wang ). The studies have some limitations such as using LOOCV as an internal test (Yi ; Bhak ) which is inconsistent when used in linear models (Shao, 1993); subsplitting the data instead of using an independent dataset for validation (Wang ); and use of small datasets (Spijker ; Watanabe ; Miyata ; Bhak ) which are biased in favor of inflated effects (Button ). We found no reports using brain tissue for the diagnostic classification of MDD.

Our accuracy findings are similar to the studies described above (0.8 in the primary brain dataset, with a range between 0.46 and 0.8 in the other brain datasets and between 0.45 and 0.78 for the blood datasets). However, our approach differed from others in that the expression set was restricted to protein-coding genes, first identified in brain and then replicated in both brain and blood. This is important because expression in brain supports the neurobiological validity of the marker, whereas expression in blood supports its potential to be eventually developed as a diagnostic aid. Conceptually, the fact that the GO analyses link the identified marker to immune, stress and defense responses is consistent with the broad literature supporting the involvement of the immune system in the pathophysiology of depression (Stein ; Maes ; Raison , 2009; Quan and Banks, 2007; Dantzer ; Dowlati ; Leonard, 2010; Hughes ; Miller and Raison, 2016). Pragmatically, because ours is a compact marker, the implementation of RT-PCR in blood seems feasible. Hence, we believe a further study of this marker in a larger dataset with PCR is warranted.

Some technical aspects negatively impacted the performance of the identified marker. First, in the brain tissue datasets, different regions were included. In brain, there tends to be high variability of gene expression in different areas (Naumova ). However, the identified marker was capable of distinguishing MDD from HC in most of the brain datasets, independently of the included brain areas (Table 1). Moreover, postmortem interval affects the RNA levels and integrity (Franz ) which could be accountable for the lower accuracy in GSE10152, even though, this dataset used RNASeq samples from the DLPFC area, like in GSE80655. Still, the gene set was able to discriminate the groups with an accuracy greater than 0.6. Additionally, inconsistency with the previous literature may, in part, be due to greater tolerance to degradation in the gene expression we identified in GSE80655, compared to some of the genes previously reported. Also, the inherent variability of microarray experiments is large (Jaksik ). As initially detected in GSE80655, we expected greater gene expression differences between MDD and HC in the replication datasets. However, although the microarray probes measured the same genes, the variability of the probes was not consistent across the datasets. We mitigated this with the selection of the probes that reflect the original assumption, with an expression difference between MDD and HC. Furthermore, the microarray experiments used different number of probes (up to 48 for GSE87610 and as low as 12 for GSE67663). This may account for the higher accuracy in GSE98793 compared to GSE67663, although the sample sizes were almost the same. In general, datasets with more probes, such as GSE38206 and GSE76826, had a higher accuracy, but this could be artifactual. Of note, variability in sample sizes was significant (e.g. 192 in GSE98793 and 18 in GSE38206). Smaller datasets achieved higher accuracy perhaps due to higher internal variability, biasing toward inflated effects (Button ). Finally, the association of gene expression between brain and blood samples was implemented in tissues originating from different subjects. Ideally, samples would come from the same subject collected at the same time (Sullivan ; Leday ). This would reduce variation in clinical characteristics that could impact expression (Ciobanu ; Yang ).

The current analysis showed that (a) there is an association between the gene expression of the identified marker from subjects with MDD and HC in the blood and brain tissues; (b) a 23 gene expression marker was able to distinguish subjects with MDD from HC with an accuracy higher than 50% in most of the databases investigated and (c) this marker includes genes related to stress and immune response, consistent with a current model of depression. Future studies in larger samples will be necessary, including brain and blood samples from the same subjects to further develop a marker with minimal sample bias, which is reproducible and universal. Though reliability between clinicians diagnosing MDD according to DSM-5 criteria was poor (Regier ), it improved with the Structured Clinical Interview for DSM-5 – Clinician Version (SCID-5-CV; kappa greater than 0.70 for MDD) (Osório ). Hence, a future study would examine the accuracy of our proposed genetic biomarker compared to the SCID-5-CV as a gold standard.

Acknowledgements

Conflicts of interest

Juan Bustillo receives royalties from UpToDate. For the remaining author, there are no conflicts of interest.