Differential co-expression and regulation analyses reveal different mechanisms underlying major depressive disorder and subsyndromal symptomatic depression

Recent depression research has revealed a growing awareness of how to best classify depression into depressive subtypes. Appropriately subtyping depression can lead to identification of subtypes that are more responsive to current pharmacological treatment and aid in separating out depressed patients in which current antidepressants are not particularly effective. Differential co-expression analysis (DCEA) and differential regulation analysis (DRA) were applied to compare the transcriptomic profiles of peripheral blood lymphocytes from patients with two depressive subtypes: major depressive disorder (MDD) and subsyndromal symptomatic depression (SSD). Six differentially regulated genes (DRGs) (FOSL1, SRF, JUN, TFAP4, SOX9, and HLF) and 16 transcription factor-to-target differentially co-expressed gene links or pairs (TF2target DCLs) appear to be the key differential factors in MDD; in contrast, one DRG (PATZ1) and eight TF2target DCLs appear to be the key differential factors in SSD. There was no overlap between the MDD target genes and SSD target genes. Venlafaxine (Efexor™, Effexor™) appears to have a significant effect on the gene expression profile of MDD patients but no significant effect on the gene expression profile of SSD patients. DCEA and DRA revealed no apparent similarities between the differential regulatory processes underlying MDD and SSD. This bioinformatic analysis may provide novel insights that can support future antidepressant R&D efforts.

Unfortunately, MDD is often resistant to standard antidepressant medication, and a large percentage of patients respond just as well to placebo [3,4]. As it is unlikely that a polymorphic syndrome like depression reflects a single disease process, recent depression research has revealed a growing awareness of how to best classify depression into depressive subtypes [5]. However, the broad DSM-based diagnosis of MDD does not encourage a search for depressive subtypes that may require their own specific treatments, and most antidepressant trials are commercially-sponsored multicenter studies that aggregate many possible subtypes under an overarching 'depression' umbrella for the sake of powering without regard to depressive subtypes [5]. In contrast, the clinical reality is that physicians regularly subtype depression when describing patients to colleagues [5]. Appropriately subtyping depression should lead to identification of subtypes that are more responsive to current pharmacological treatment and aid in separating out depressed patients in which current antidepressants are not highly effective.
One reason for this failure at depression subtyping is our insufficient understanding of the various neurobiological processes underlying depression [6]. An improved understanding of the pathophysiological mechanisms underlying different depressive subtypes should aid antidepressant pharmacological development. For instance, previous studies by our group have applied proteomic, transcriptomic, and metabololomic approaches to better characterize MDD and rodent models of depression [7][8][9][10][11]. Other advanced technologies such as DNA microarrays and next-generation sequencing can allow for rapidly acquiring detailed biochemical information about DNA polymorphisms and transcriptome profiles in different depressive subtypes [12]. These studies collectively show profound biochemical changes in depressed patients, and animal depression models partly corroborate or complement the alterations identified in depressed patients. Transcriptional changes in the course of MDD and the effects of antidepressant drugs on multiple disease-related transcriptional pathways have also been surveyed. The combination of both can unravel additional mechanisms of disease etiology and can provide a 'bottom-up' approach for the discovery of novel antidepressant drugs.
One depressive subtypesubsyndromal symptomatic depression (SSD)has been identified as a transitory depressive state in the depression spectrum. In this study, we hypothesized that the transcriptomic profiles of leukocytes derived from drug-free first-episode SSD patients and MDD patients were significantly different in these two depressive subtypes. Additionally, we hypothesized that these transcriptomic profiles were significantly different before and after treatment with the popular antidepressant venlafaxine (Efexor™, Effexor™).
Applying differential co-expression analysis (DCEA) and differential regulation analysis (DRA), the transcriptomic profiles of leukocytes derived from drug-free first-episode SSD patients, MDD patients, and healthy controls were compared using whole-genome cDNA microarrays in order to discover venlafaxine's mechanism (s) of action in MDD and SSD patients.

Expression dataset
We reutilized the GSE32280 dataset [13], which was published on GEO (http://www.ncbi.nlm.nih.gov/geo/query/ acc.cgi?acc=gse32280), to unveil the dysregulated mechanisms in MDD and SSD and explore the differential processes of venlafaxine action on MDD and SSD. The GSE32280 data consists of whole blood gene expression profiles of five SSD patients and eight MDD patients before venlafaxine treatment (hereinafter 'pre-treatment') and after venlafaxine treatment (hereinafter 'post-treatment') as well as eight healthy controls (Table 1). Raw data were normalized by the RMA method and log2 transformed. After filtering out ambiguous genes, 20283 gene expression values were obtained for all samples.

Differential Co-Expression Analysis (DCEA)
In transcriptomics, differential co-expression analysis (DCEA) has emerged as a unique complement to traditional differential expression analysis (DEA) [14]. Rather than merely calculating expression level changes in individual genes as DEA does, DCEA analyzes differences in gene interconnection by calculating the expression correlation changes of gene pairs between two conditions [14][15][16]. The rationale underlying DCEA is that changes in gene co-expression patterns between contrasting phenotypes (in this case, SSD versus MDD versus healthy control) provide insight into the disrupted regulatory relationships or affected regulatory subnetworks specific to the phenotype (s) of interest.
We utilized the DCGL software package to conduct DCEA, which is a software program released as an R package including two gene filtering functions, three link filtering functions, and five DCEA functions [17,18] In the DCGL software package, differential co-expression profile (DCp) and differential co-expression enrichment (DCe) are designed based on the exact co-expression changes of gene pairs, and thus can differentiate significant co-expression changes from relatively trivial ones and identify co-expression reversal between positive and negative.
For example, for the MDD vs. healthy control comparison, the DCGL package constructed two co-expression networks by calculating gene pair-wise correlation coefficients based on MDD samples and healthy samples. Then, a cutoff for correlation coefficients was applied to filter out some gene pairs and retain the most relevant gene pairs to form two co-expression networks. These two coexpression networks possessed the same topological structure with different correlation coefficients for each gene pair. In these co-expression networks, genes that connected with other genes were deemed to be neighboring genes. Then, DCp was employed to calculate the change between the two co-expression networks that was defined by different co-expression (dC) (see Eq. 1 below). In order to evaluate the statistical significance of dC, DCp implemented a permutation test, in which we shuffled 20283 genes and samples randomly 500 times. Then, we estimated 500 pseudo dCs for every gene using the permuted data. Finally, the p-value for each gene could be estimated, and the FDR value could be obtained accordingly. Differentially co-expressed genes (DCGs) were identified by p-value and/or FDR threshold. DCe used the limit fold change (LFC) model [19] to sort out different coexpression gene pairs or differentially co-expressed links (DCLs) and employed a binomial probability model to identify DCGs based on enrichment of DCLs (see Eq. 2 below).
where x i1 , x i2 … x in (y i1 , y i2 … y in ) denote correlation coefficients between gene i with its n neighbors in condition x (y) and n denoting the number of neighbors. In DCe, N represents the number of links in each coexpression network, K represents the number of DCLs determined by the LFC model. Eq. 2 calculates the enrichment of DCL for gene i with n i links of which k i are DCLs.

Differential regulation analysis (DRA)
Among the many growing directions of DCEA, there is differential regulation analysis (DRA), which integrates the transcription factor-to-target (TF2target) information to probe upstream regulatory events that account for the observed co-expression changes [20,21]. DRA can unveil the central regulatory network which straightforwardly reflected the differential regulation mechanisms of two contrastive conditions via differential co-expression analysis results. DCGL package was also used in the differential regulation analysis step. The DCGL package contains a library of transcription factors known to regulate the relationships between human target genes termed the transcription factor-to-target library (TF2targetlibrary), which includes 214607 binary tuples involving 215 human transcription factors and 16863 targets [18]. Using this TF2target library, the DRsort function of DCGL was used to scrutinize the DCGs and DCLs against the transcription factor-to-target information to highlight a subset of the genes and links that are potentially highly related to the putative differential regulation mechanisms. If a DCG coincides with a transcription factor (TF), it is regarded as a differentially regulated gene (DRG) based on the likelihood that a differential co-expression of this type of DCG could be attributed to disrupted regulatory relationships between the transcription factor and its targets. Hereinafter, 'TF2target DCLs' refers to differentially co-expressed gene links or pairs (DCLs) that coincide with transcription factor-to-target relations. Our rationale here is that the disruption of regulatory relations can affect the co-expression links between a regulator and its targets [18].

Application of DCEA and DRA
To explore the underlying regulatory mechanisms in MDD, we conducted DCEA and DRA to compare the two different MDD states (pre-treatment and post-treatment), pre-treatment MDD versus healthy control, and posttreatment MDD versus healthy control. To explore the effects of venlafaxine therapy on MDD patients, we conducted DCEA and DRA to compare post-treatment MDD versus healthy controls and pre-treatment MDD versus post-treatment MDD. An identical analysis was performed on SSD samples as well. First, we identified DCGs thorough the DCp function, and DCLs and another list of DCGs through the DCe function. Then, in order to avoid lack of individual method, we sorted out true DCGs via overlapping DCp's DCGs and DCe's DCGs. DCLs were defined as true DCLs if one of the genes among the DCLs was identified in overlapping DCGs. Finally, the TF2target library was employed to highlight DCGs and DCLs. If a DCG coincided with a TF, it was listed as a DRG. If a DCL coincided with a TF-regulated target gene in the TF2target library, it was listed as a TF2target DCL or DRL.

Differential regulatory mechanisms in MDD
In order to reveal the underlying regulatory mechanisms of non-medicated MDD, samples of pre-treatment MDD versus healthy controls were compared. First, we identified 13 putative DRGs (termed DRG MDD-pre that included FOSL1, FOXL1, MEF2A, HNF1A, IRF1, JUN, SOX9, SRF, TFAP4, TFCP2, TLX2, HLF, and ZNF423) and 49 putative TF2target DCLs (termed DRL MDD-pre that are displayed in Table 2). Furthermore, by comparing post-treatment MDD and healthy control samples, 10 putative DRGs (termed DRG MDD-pos that included BPTF, EP300, FOXI1, IL10, NFKB1, NFYC, NR3C1, TP53, USF2, and FOXL1) and 294 putative TF2target DCLs (termed DRL MDD-post that are shown in Additional file 1: Table S1) were identified. In order to exclude differences within individuals from different samples, the overlapping elements between pre-treatment MDD versus healthy controls and posttreatment MDD versus healthy controls were deemed to indicate differential regulatory information among individuals. Therefore, the expression results that excluded the overlapping differences in co-expression of pre-treatment MDD versus healthy controls were deemed to be the DRGs and TF2target DCLs responsible for differentiating MDD and healthy controls (MDD panel, Figure 1). Therefore, DRG MDD = DRG MDD − pre − (DRG MDD − pre ∩ DRG MDD − post ) were deemed to be the true MDDrelated DRGs (termed DRG MDD that included FOSL1, MEF2A, HNF1A, IRF1, JUN, SOX9, SRF, TFAP4, TFCP2, TLX2, HLF, and ZNF423). These 12 true DRG MDD are depicted in Figure 1 by the lime green area of pre-treatment/control that excludes the overlapping brown area of pre-treatment/control with posttreatment/control. Based on the foregoing analysis, these 12 DRG MDD may be involved in the pathological processes underlying MDD. According to the same logical analysis applied to DRG, we identified 48 true MDD-related TF2target DCLs (i.e., Denotes DCG in each TF2target DCL. Some rows were discarded (strikethrough) since they are individual basis. Rows in bold represent the 16 key-DRL MDD . The p. DCG column denotes the p-value of DCG for each TF2target DCL; the "Cor.1" and "Cor.2" columns represent the correlation coefficients of the two respective conditions. The "Type" column reflects the mode of regulation with "1" denoting a positive/negative regulation in disease samples switched to a negative/positive regulation in control samples and "0" denoting no change in regulation.  (Table 2 except striked-through rows). By overlapping the 12 DRG MDD with the 48 DRL MDD , six key DRGs (termed key-DRG MDD that included FOSL1, SRF, JUN, TFAP4, SOX9, and HLF) and 16 TF2target DCLs (termed key-DRL MDD ) appear to be the key differential factors in MDD ( Figure 2). GO analysis revealed that the six key-DRG MDD are involved in RNA metabolism (SRF, JUN, HLF, and FOSL1) and developmental processes (SRF and JUN) (Additional file 1: Table S3, Figure S1). Moreover, four key-DRG MDD (SRF, JUN, HLF, and FOSL1) localize to the nucleus and are enriched for transcription regulatory activity (Additional file 1: Table S3, Figure S1).

Potential mechanism of venlafaxine in MDD
In order to discover the pharmacological mechanism (s) of venlafaxine in MDD, samples from the same MDD patients pre-treatment and post-treatment were compared. Consequently, a total of nine DRGs and 532 TF2target DCLs were discovered. The 532 TF2target DCLs contained 65 transcription factors and 466 target genes (Additional file 1: Table S2). Nine differential regulatory transcription factors, which participated in TF2target DCLs, appear to be the key differential regulatory relation pairs reflecting the pharmacological action of venlafaxine in MDD ( Figure 3).

Differential regulatory mechanism in SSD
Similar to the analysis conducted in MDD, the DCGL software package identified three putative DRGs (termed DRG SSD-pre that included FOSB, PATZ1, and TFAP4) and 10 putative TF2target DCLs (termed DRL SSD-pre that are displayed in Table 3) between pre-treatment SSD and healthy control samples. A comparison of posttreatment SSD with healthy control samples yielded five DRGs (termed DRG SSD-post that included EGR1, IRF1, MYOD1, SOX5, and TFAP4) and 15 TF2target DCLs (termed DRL SSD-post ). After exclusion on an individual basis, i.e., and DRL SSD = DRL SSD − pre − (DRL SSD − pre ∩ DRL SSD − post ), two DRGs (DRG SSD = PATZ1 and FOSB) and nine DRL SSD appear to be associated with SSD. The nine DRL SSD are indicated by the lime green area that excludes the overlapping brown area (SSD panel, Figure 1) and listed in Table 3 (except striked-through rows). Applying a similar analysis to that used in MDD by overlapping the 2 DRG SSD with the 9 DRL SSD , we determined one key-DRG SSD (PATZ1) and 8 key-DRL SSD in SSD (Figure 4).

Potential mechanism of venlafaxine in SSD
No DRG or TF2target DCLs were identified between SSD pre-and post-treatment. Therefore, the expression profiles between pre-and post-treatment SSD samples were not significantly different; thus, venlaxafine appears to have no significant effect on the gene expression profile of SSD patients.

Discussion
With regard to the differential regulatory mechanism (s) underlying MDD, six key-DRG MDD (FOSL1, SRF, JUN, TFAP4, SOX9, and HLF) and 16 key-DRL MDD (TF2target DCLs) appear to be the key differential factors in MDD (Table 2). FOSL1 has been previously associated with addiction, depression, and anxiety [22]. SRF and JUN (c-JUN) are functionally involved in the MAPK signaling pathway (Additional file 1: Figure S2) [13]. GO analysis revealed that the six key-DRG SSD are involved in RNA metabolism (SRF, JUN, HLF, and FOSL1) and developmental processes (SRF and JUN) (Additional file 1: Table S3, Figure S1). Moreover, four key-DRG MDD (SRF, JUN, HLF, and FOSL1) localize to the nucleus and are enriched for transcription regulatory activity (Additional file 1: Table S3, Figure S1). With regard to the differential regulatory mechanism (s) underlying SSD, one key-DRG SSD (PATZ1) and eight key-DRL SSD (TF2target DCLs) appear to be the key differential factors in SSD (Table 3). Unfortunately, after an extensive literature search, we found no previous studies showing any association between these genes and SSD, which may be attributed to the paucity of research on SSD. These initial findings can open an entirely new field of investigation on the underlying regulatory mechanism (s) of SSD.
As for a potential mechanism (s) underlying venlafaxine in MDD and SSD, nine DRGs and 532 TF2target DCLs were responsible for distinguishing pre-treatment and post-treatment MDD samples (Additional file 1: Table S2). However, no DRGs or TF2target DCLs were found between pre-treatment and post-treatment SSD samples, suggesting that venlaxafine appears to have a significant effect on the gene expression profile of MDD patients but no significant effect on the gene expression profile of SSD patients.

Different mechanisms underlying MDD and SSD
DCEA and DRA revealed no similarities between the differential regulatory processes underlying MDD and Figure 2 Key TF2target DCLs in MDD. Every line with an arrow represents a differential regulatory relationship in MDD. Squares represent transcription factors (yellow = differential, blue = non-differential), circles represent target genes (yellow = differential, blue = non-differential), and lines represent the action of transcription factors on target genes (green = positive regulation in disease samples transformed to negative regulation in control samples or vice-versa, red = unchanged regulation). Denotes DCG in each TF2target DCL. Some rows were discarded (strikethrough) since they are individual basis. Rows in bold represent the 8 key-DRL SSD . The p. DCG column denotes the p-value of DCG for each TF2target DCL; the "Cor.1" and "Cor.2" columns represent the correlation coefficients of the two respective conditions. The "Type" column reflects the mode of regulation with "1" denoting a positive/negative regulation in disease samples switched to a negative/positive regulation in control samples and "0" denoting no change in regulation.

Figure 3
Key TF2target DCLs of Venlafaxine in MDD. Nine differentially regulated genes (yellow squares) display the action of venlafaxine in MDD. Squares represent transcription factors (yellow = differential, blue = non-differential), circles represent target genes (yellow = differential, blue = non-differential), and lines represent the action of transcription factors on target genes (green = positive regulation in disease samples transformed to negative regulation in control samples or vice-versa, red = unchanged regulation).
SSD. Specifically, six DRGs (FOSL1, SRF, JUN, TFAP4, SOX9, and HLF) and 16 TF2target DCLs appear to be the key differential factors in MDD; in contrast, one DRG (PATZ1) and eight TF2target DCLs appear to be the key differential factors in SSD (Tables 2, 3). Moreover, there was no overlap in MDD target genes and SSD target genes. While there were 62 downstream target genes for MDD or SSD, there was no overlap in MDD target genes and SSD target genes. In addition, via the Wilcoxon test, we found that the distance between TFs and target genes in SSD were significantly larger than those in MDD (p = 0.02). Therefore, although MDD appears to include more key-DRG MDD and key-DRL MDD as compared to the number of key-DRG SSD and key-DRL SSD in SSD, the extent of differential gene regulation in MDD was lower than that in SSD. Based on this evidence, we hypothesize that MDD involves preferential influence to a wider set of genes, while SSD involves preferential influence of gene regulatory elements.
With regard to venlafaxine's mode of action in MDD, nine transcription factors (ARID5B, ATF6, BPTF, GATA3, HAND1, IL10, NFE2L1, NFYC, and RFX1) were found to act on 256 target genes. In contrast, no differences between pre-treatment and post-treatment SSD were identified, suggesting that venlafaxine has no significant effect on the gene expression profile of SSD patients. Therefore, as the aforementioned evidence reveals that MDD and SSD display completely different underlying mechanisms, we presume that venlaxafine has no significant effect on SSD patients due to the different underlying mechanisms of the two depressive subtypes. In future studies, conducting a similar analysis involving additional antidepressants would provide additional evidence to support this view and may provide insight into other antidepressant therapies or potential drug targets that may be more efficacious for SSD patients.

Numerous hypotheses for MDD
The pathoetiology underlying MDD remains largely unknown [23]. MDD can spontaneously develop but often follows a traumatic emotional experience or can be a symptom of other diseases, most often neurological (e.g., stroke, multiple sclerosis, and Parkinson's disease) or endocrine (e.g., Cushing's disease and hypothyroidism). MDD can also be triggered by pharmacological agents or drug abuse [24]. Different mechanisms have been proposed to explain the pathophysiological basis of MDD from a neurobiological point of view [25,26] These hypotheses include monoaminergic deficiency, hypothalamic-pituitary-adrenal (HPA) axis dysregulation [27][28][29], neurogenetic and neurotrophic-growth factor impairment, metabolic disturbances, circadian rhythm desynchronization, and inappropriate stimulation of the immune system [25,26,30]. This multiplicity of hypotheses can be explained by the fact that several of them are intertwined rather than being mutually exclusive [27][28][29][30]. Moreover, it is probable that different endophenotypes correlate with different neurobiological adaptations. From this vantage point, a lack of criteria in selecting appropriate patients may contribute to the difficulties experienced in the quest for new therapies that rely on non-monoamine mechanisms of action.

Lighting the future for antidepressant R&D
The efficacy of available antidepressant therapies has been demonstrated with respect to placebo, especially in cases of severe depression [31]. Nevertheless, a relatively large number of patients still fail to respond to conventional treatment [32]. Moreover, a number of symptoms may not be adequately resolved in patients that experience an overall positive therapeutic response. In addition, although the safety profile of newer antidepressant agents is quite superior to older drugs, the incidence of side effects, such as sexual dysfunction, may cause therapy discontinuation, especially in younger patients. Therefore, considerable efforts have been dedicated to developing therapeutic agents that address novel neurobiological targets with the hope of overcoming the aforementioned issues [33,34].
Unfortunately, approaches aimed at identifying therapies based on different mechanisms have not been successful Figure 4 Key TF2target DCLs in SSD. Eight TF2target DCLs showing the differential regulatory mechanisms in SSD. Squares represent transcription factors (yellow = differential, blue = non-differential), circles represent target genes (yellow = differential, blue = non-differential), and lines represent the action of transcription factors on target genes (green = positive regulation in disease samples transformed to negative regulation in control samples or vice-versa, red = unchanged regulation). [35], and pharmaceutical companies are disengaging from this disease area, which has been perceived as highly risky. One reason is our insufficient understanding of the neurobiological basis of MDD [26,30]. As discussed earlier, MDD is likely a constellation of merging disease states that can be split into endophenotypes [27]. Several lines of evidence support the notion that the occurrence of environmental challenges, often in the form of stressful experiences, needs to be associated with a pre-existing genetic predisposition to bring about the disease [36]. The bioinformatic analysis in this report provides reliable evidence that can support future antidepressant R&D efforts.

Study limitations
First, the sample size of the current study was limited; thus, the low power of this study limits the applicability of our findings. Second, this study only employed two methods -DCEA and DRAto analyze the expression profiles of the samples. Although this did provide novel evidence that can aid in pharmacological target development, no biological experiments were performed to validate the findings. Therefore, further biological experiments with larger sample sizes are required to validate the current findings. Third, it would have been more insightful to measure clinical response to venlaxafine in terms of improvements in depression rating scales in lieu of assessing differential expression preand post-treatment without regard to clinical response. Unfortunately, there were no clinical response data reported in the previously published studies. Therefore, for our future research on this topic, we will collect patient samples from our hospital and record the clinical responses in terms of improvements in depression rating scales.