Pilot Study of Metabolomic Clusters as State Markers of Major Depression and Outcomes to CBT Treatment

Major depressive disorder (MDD) is a common and disabling syndrome with multiple etiologies that is defined by clinically elicited signs and symptoms. In hopes of developing a list of candidate biological measures that reflect and relate closely to the severity of depressive symptoms, so-called “state-dependent” biomarkers of depression, this pilot study explored the biochemical underpinnings of treatment response to cognitive behavior therapy (CBT) in medication-free MDD outpatients. Plasma samples were collected at baseline and week 12 from a subset of MDD patients (N = 26) who completed a course of CBT treatment as part of the Predictors of Remission in Depression to Individual and Combined Treatments (PReDICT) study. Targeted metabolomic profiling using the AbsoluteIDQ® p180 Kit and LC-MS identified eight “co-expressed” metabolomic modules. Of these eight, three were significantly associated with change in depressive symptoms over the course of the 12-weeks. Metabolites found to be most strongly correlated with change in depressive symptoms were branched chain amino acids, acylcarnitines, methionine sulfoxide, and α-aminoadipic acid (negative correlations with symptom change) as well as several lipids, particularly the phosphatidlylcholines (positive correlation). These results implicate disturbed bioenergetics as an important state marker in the pathobiology of MDD. Exploratory analyses contrasting remitters to CBT versus those who failed the treatment further suggest these metabolites may serve as mediators of recovery during CBT treatment. Larger studies examining metabolomic change patterns in patients treated with pharmacotherapy or psychotherapy will be necessary to elucidate the biological underpinnings of MDD and the -specific biologies of treatment response.

Major depressive disorder (MDD) is a common and disabling syndrome with multiple etiologies that is defined by clinically elicited signs and symptoms. In hopes of developing a list of candidate biological measures that reflect and relate closely to the severity of depressive symptoms, so-called "state-dependent" biomarkers of depression, this pilot study explored the biochemical underpinnings of treatment response to cognitive behavior therapy (CBT) in medication-free MDD outpatients. Plasma samples were collected at baseline and week 12 from a subset of MDD patients (N = 26) who completed a course of CBT treatment as part of the Predictors of Remission in Depression to Individual and Combined Treatments (PReDICT) study. Targeted metabolomic profiling using the AbsoluteIDQ R p180 Kit and LC-MS identified eight "co-expressed" metabolomic modules. Of these eight, three were significantly associated with change in depressive symptoms over the course of the 12-weeks. Metabolites found to be most strongly correlated with change in depressive symptoms were branched chain amino acids, acylcarnitines, methionine sulfoxide, and α-aminoadipic acid (negative correlations with symptom change) as well as several lipids, particularly the phosphatidlylcholines (positive correlation). These results implicate disturbed bioenergetics as an important state marker in the pathobiology of MDD. Exploratory analyses contrasting remitters to CBT versus those who failed the treatment

INTRODUCTION
Major depressive disorder (MDD) is a clinical syndrome that has multiple etiologies and responds to a diverse range of treatments that affect various biological pathways. Nevertheless, it is highly likely that specific biological processes underpin the clinical presentation of the disorder. Identifying these state biological processes could provide a more precise gauge of the pathophysiological processes underpinning the clinical-symptomatic expression of MDD, and that also could reflect treatments' biological effect beyond the information gleaned from the simple assessment of signs and symptoms (Rush and Ibrahim, 2018).
Biological, physiological, neuro-functional, and other measures that are most closely tied to and reflective of the clinical expression of MDD or to the symptomatic expression of other medical syndromes are often referred to "state-dependent" markers (Rush and Ibrahim, 2018). Central venous pressure (CVP), for example, is a state-dependent measure for congestive heart failure (CHF). The greater the CVP, the more severe the symptoms that define CHF such as pedal edema, pulmonary effusion, orthopnea, and dyspnea. On the other hand, "trait-like" markers, are those measures that are persistently abnormal both during and between clinically symptomatic episodes (Rush and Ibrahim, 2018). "Trait-like" markers often reflect the underlying pathobiology of the condition that either sets the stage for the initial clinical expression of the disorder or that reflect the effect/consequence of the clinical episode itself even after the episode ends. The latter are sometimes said to be "scars" or consequences of the clinical episode. Left or right ventricular hypertrophy, for instance, can be consequences of repeated episodes of congestive heart failure (Senni and Redfield, 1997). For MDD, hypercortisolemia is known to be highly state dependent in psychotic or melancholic depressions (Pariante, 2017). On the other hand, some sleep EEG parameters appear to be more trait-like (i.e., persistent even between clinically apparent major depressive episodes) than state-dependent (only apparent during clinically apparent depressive episodes) (Kraemer et al., 1994;Thase et al., 1998). However, neither state nor trait markers for MDD have been found that are as yet of sufficient value to enter clinical practice.
Metabolomics have the potential to define specific biochemical processes that underpin MDD and the effects of treatment on those biological processes. Some work indicates metabolomic profiling may have utility in differentiating MDD from bipolar disorder or healthy controls (Gui et al., 2018;Hashimoto, 2018;Pan et al., 2018;MacDonald et al., 2019). Others have shown pretreatment metabolomic profiles and changes during treatment can predict response to antidepressant medication, or changes in profiles occur (Kaddurah-Daouk et al., 2011;Abo et al., 2012;Kaddurah-Daouk et al., 2013;Zhu et al., 2013;Rotroff et al., 2016;Moaddel et al., 2018;Czysz et al., 2019). To date, however, metabolomic profiling has largely been conducted with depressed patients who have been taking antidepressant medications, which can directly affect metabolomics profiles and thus interfere with the identification of state dependent measures.
Pharmacometabolomic studies from our group have previously reported perturbations in intermediates of TCA cycle, urea cycle, amino acids, and lipids in depressed patients exposed to sertraline (Kaddurah-Daouk et al., 2011Zhu et al., 2013). Another study utilizing intravenous ketamine treatment in depressed patients reported changes in tryptophan metabolism, acylcarnitines, urea cycle intermediates, and lipid metabolism (Rotroff et al., 2016). In a cross-sectional study, the branched chain amino acids (BCAAs) Valine, Leucine, and Isoleucine were significantly lower in MDD patients compared to healthy controls and were negatively correlated to Hamilton Depression Rating Scale scores. In a rat model of depression, biogenic amines like putrescine, spermine, and spermidine were significantly reduced in the hippocampus of stressed animals compared to non-stressed ones, but the biogenic amines were restored by the antidepressant effect of S-adenosyl-L-methionine (Genedani et al., 2001). Plasma lipid and acylcarnitine profiles, which have also been implicated in animal models of depression, suggest inflammatory conditions and incomplete mitochondrial β-oxidations as primary phenomena associated with the pathophysiology of MDD (Chen et al., 2014).
This pilot study utilized a sample from the cognitive behavior therapy (CBT) arm of the Emory Predictors of Remission in Depression to Individual and Combined Treatments (PReDICT) study, a randomized controlled trial of previously untreated patients with MDD (Dunlop et al., 2017a). This sample avoids the likely confounding effects of medications on endogenous metabolomic processes under study; i.e., the sample ensures that when patients improve symptomatically-(or not), there is no confounding effect of concurrent antidepressant medication. Therefore, any biological changes whether found to be state-independent or state dependentwould entirely reflect the depressive symptom severity while being unaffected by the pharmacological effects of medications. To our knowledge this study is the first to explore the underlying biochemical (metabolomic) changes associated with relevant clinical endpoints in treatment-naïve patients treated with CBT.
Analyses were conducted to identify which changes in serum metabolites and pathways were most closely related to changes in depressive symptoms between baseline to week 12 in the acute treatment of MDD with CBT. Herein, we report observed metabolite alterations within a biomarker panel targeting 186 plasma metabolites from 5 distinct metabolite classes (specifically, amino acids, biogenic amines, acylcarnitines, glycerophospholipids, and sphingolipids) available in the Biocrates AbsoluteIDQ R p180 Kit 1 .

Clinical
The PReDICT study was a randomized clinical trial that enrolled 344 adults ages 18-65 years with a primary psychiatric diagnosis of MDD without psychotic features. The design and clinical results of the study have been published (Dunlop et al., 2012(Dunlop et al., , 2017a. The Structured Interview for DSM-IV (First et al., 1995) assessed MDD diagnosis, which was confirmed with a psychiatrist's interview. Patients meeting all eligibility criteria were randomized in a 1:1:1 manner to receive either CBT (delivered in up to 16 one-hour individual sessions), escitalopram, or duloxetine for 12 weeks. One-hundred-fifteen patients were assigned to CBT, of whom 26 had serum samples available for metabolomic analyses at baseline and week 12; these 26 patients are the subjects of the current analysis.
Key inclusion criteria for the trial included no lifetime history of having received treatment for depression (either ≥4 weeks of antidepressant medication at a minimally effective dose or ≥4 sessions of an evidence-based psychotherapy), and fluency in either English or Spanish. At screening, patients had to score ≥18 on the HAM-D 17 (Hamilton, 1967) and at the baseline randomization visit had to score ≥15. Key exclusion criteria included: a lifetime history of bipolar disorder, psychotic disorder, or dementia; a current significant medical condition that could affect study participation or data interpretation; a diagnosis of obsessive-compulsive disorder, an eating disorder, substance dependence, or dissociative disorder in the 12 months before screening; or substance abuse within the 3 months prior to baseline. The only other psychotropic agents permitted during the trial were sedatives (eszopiclone, zolpidem, zaleplon, melatonin, or diphenhydramine) up to three times per week. The Emory Institutional Review Board and the Grady Hospital Research Oversight Committee approved the study protocol, and all patients provided written informed consent prior to beginning study procedures.
The therapy was delivered in accordance with Beck's protocolbased CBT (Beck et al., 1979) and therapists' fidelity to the protocol was assessed by independent raters at the Beck Institute using the Cognitive Therapy Scale (Young and Beck, 1980). Raters blinded to treatment assignment assessed depression severity using the HAM-D 17 at baseline, weeks 1-6, 8, 10, and 12. For the individual patient outcomes, the protocol defined remitters as patients who achieved HAM-D 17 score ≤7 at both week 10 and 12 (Dunlop et al., 2012). Consistent with the prior analyses of this dataset (Dunlop et al., 2017a), outcomes for nonremitters were using percent change in HAM-D 17 score from 1 https://www.biocrates.com/products/research-products/absoluteidq-p180-kit baseline to week 12, as follows: Non-remitting responder, ≥50% reduction, but not meeting remitter criteria; Partial responder: 30-49% reduction; Treatment failure: <30% reduction.

Metabolomic Profiling Using Absolute IDQ p180 Kit
Metabolites were measured with a targeted metabolomics approach using the AbsoluteIDQ R p180 Kit (BIOCRATES Life Science AG, Innsbruck, Austria), with an ultra-performance liquid chromatography (UPLC)/MS/MS system [Acquity UPLC (Waters), TQ-S triple quadrupole MS/MS (Waters)]. This procedure provides measurements of up to 186 endogenous metabolites in quantitative mode (amino acids and biogenic amines) and semi-quantitative mode (acylcarnitines, sphingomyelins, phosphatidylcholines and lysophosphatidylcholines across multiple classes). The AbsoluteIDQ R p180 kit has been fully validated according to European Medicine Agency Guidelines on bioanalytical method validation. Additionally, the kit plates include an automated technical validation to assure the validity of the run and provide verification of the actual performance of the applied quantitative procedure including instrumental analysis. The technical validation of each analyzed kit plate was performed using MetIDQ R software based on results obtained and defined acceptance criteria for blank, zero samples, calibration standards and curves, low/medium/high-level QC samples, and measured signal intensity of internal standards over the plate. De-identified samples were analyzed following the manufacturer's protocol, with metabolomics labs blinded to the clinical data.

Preprocessing of P180 Profiles
The raw metabolomic profiles included 182 metabolite measurements of serum samples. Each assay plate included a set of duplicates obtained by combining approximately 10 µl from the first 76 samples in the study (QC pool duplicates) to allow for appropriate inter-plate abundance scaling based specifically on this cohort of samples (n = 24 across all plates). Metabolites with >40% of measurements below the lower limit of detection (LOD) were excluded from the analysis (n = 160 metabolites passed QC filters). To adjust for the batch effects, a correction factor for each metabolite in a specific plate was obtained by dividing the metabolite's QC global average by QC average within the plate. Missing values were imputed using each metabolite's LOD/2 value followed by log 2 transformation to obtain a normal distribution of metabolite levels. The presence of multivariate outlier samples was checked by evaluating the squared Mahalanobis distance of samples. Samples were flagged as "outliers" when their Mahalanobis distances exceeded the critical value corresponding to a Bonferroni-corrected threshold (0.05/n, n: number of samples) of the Chi-square distribution with m degrees of freedom (m = 160: number of metabolites).

Data Analysis
For statistical analysis we adopted a two-pronged approach. Initially, a multivariate "co-expression network" analysis was employed with CBT treated patients to detect clusters (or modules) of metabolites demonstrating similar patterns of perturbations that correlated with changes in depressive symptom scores based on the HAM-D 17 . Univariate analyses were also performed to detect whether the metabolites within or outside the clusters were individually and significantly correlated to the depressive symptom outcome. The traditional univariate analysis method to metabolomic profiling focuses on the individual metabolites; thus, the interactions among metabolites are largely ignored, even though it is appropriate to assume that metabolites play their roles not in isolation but via interactions with each other. Consequently, metabolite "co-expression" analysis is a powerful, multivariate approach to identify groups of perturbed metabolites belonging to same class or pathways. This approach has the additional benefit of alleviating the multiple testing problem (DiLeo et al., 2011). The workflow of data analysis is presented in Figure 1.

Univariate analysis
To define the association between changes in metabolite levels from baseline to week 12 of CBT treatment and the changes in depressive symptom of total HAM-D 17 scores over that time, linear mixed effects models were fitted to each metabolite change, adjusting for age and gender and with subjects as a random variable. All p-values were checked for false discovery rates by Benjamini-Hochberg method (Hochberg and Benjamini, 1990). Correlation between metabolite changes and depressive symptoms changes were also assessed by Pearson's correlation coefficients.

Co-expression network analyses
Changes in modules of "co-expressed" metabolites were identified using the R package WGCNA (weighted gene co-expression network analysis) (Langfelder and Horvath, 2008). Signed and weighted Pearson's correlation networks were constructed with the subject-wise changes of baseline to week 12 metabolite concentrations (in the logarithmic scale). First a weighted adjacency matrix was created based on pairwise Pearson's correlation coefficients between the metabolites. A scale-free topology criterion was used to choose the soft threshold of beta = 18 for the correlations as per the WGCNA protocol. The obtained adjacency matrix was used to calculate the topological overlap measure (TOM) for each pair of metabolite log2 fold changes comparing their adjacencies with all of the other metabolite log2 fold changes. Densely interconnected groups (or modules) of metabolites were identified by hierarchical clustering using 1-TOM as a distance measure through the use of the dynamic hybrid tree cut algorithm with a deep split of 2 and a minimum cluster size of 3. Each module is summarized by the module eigenvector, which is the first principal component of the metabolite changes across all the subjects. Similar clusters were subsequently merged if the correlation coefficient between the clusters' eigenvectors exceeded 0.75. The association between the resultant modules and the changes in HAM-D 17 scores was measured by the pairwise Pearson correlation coefficients and presented in a heatmap. In all analyses for this small scale pilot study an uncorrected p-value threshold of 0.10 was used as the significance cutoff.

Patient Characteristics
Plasma metabolite data were available at baseline and week 12 from 26 patients. The mean number of therapy sessions attended was 14.0 ± 1.5. Table 1 summarizes the characteristics of the study sample. Four patients (15.4%) were in a chronic depressive episode. Twelve (46.2%) of the sample achieved remission and 7 (26.9%) were classified as treatment failures.

Detection of Metabolite "Co-expression" Modules
To investigate the functional response of the MDD metabolome during receipt of CBT, we adopted a multivariate approach. Using WGCNA methodology we focused on identifying modules (or clusters) of metabolites that showed a similar pattern of change from baseline to week 12. Thus, each module represented metabolite changes (week12/baseline ratios) in the logarithmic scale. Eight such metabolite modules were identified in which the member metabolites showed statistically significant strong correlations (mean R 2 ranged between 0.74 and 0.94, all p < 0.05) amongst each other in their perturbation patterns, and each module was assigned a unique color. Black, blue, brown, greenyellow, midnight-blue, purple, royal-blue and yellow were the metabolite modules representing 8, 15, 12, 96, 5, 6, 3 and 9 metabolites, respectively. The grey module represented 6 metabolites that could not be assigned to any module. The detected modules were represented by metabolites belonging primarily to the same metabolite class; this may indicate that these metabolites have a functional relationship to each other. Additionally, for each module we identified a "hub" metabolite (also known as a "driver" metabolite) that had the maximum number of connections in the module. The hub metabolites are important, and they merit further investigation because they may influence the function of other metabolites, or even may be significant contributors to the trait of interest. The eight modules, their hub metabolites, and their major metabolite classes are presented in Table 2. A list of the network metrics, the intra-modular correlations between member metabolites, and the module membership for each metabolite is presented in Supplementary Files (Supplementary Tables S1-S3).

Metabolite Modules That Were Associated With Changes in Depressive Symptoms (HAM-D 17 Scores)
Next, we evaluated the association between the identified metabolite modules and changes in HAM-D 17 scores from baseline to week 12. Three metabolite modules were found to be significantly associated (R 2 > 0.3, at p < 0.1) with changes in the symptom severity scores: (a) the purple module containing the short chain acylcarnitines (C3, C4, and C0), α-aminoadipic acid, and the two amino acids, Glutamate and Proline; (b) the yellow module containing the BCAAs, Isoleucine and Valine, the BCAAderived C5-carnitine (Isovalerylcarnitine), the neurotransmitterrelated amino acids Tryptophan, Tyrosine, Phenylalanine, Methionine, Methionine-sulfoxide, and the biogenic amine Sarcosine; and (c) the green-yellow module containing 96 lipid molecules including the phosphatidylcholines, sphingomyelins, and acylcarnitines. A heatmap showing the correlations between each of the metabolite modules and the changes in HAM-D 17 scores is presented in Figure 2. The correlation coefficients were moderate, ranging between 0.31 and 0.36 for each of the three modules. We examined the correlation of each of the metabolite members of the purple, yellow and green-yellow modules to HAM-D 17 scores. Figure 3A presents a composite plot of the correlations between the changes of each member-metabolite of the purple module to each other and also to HAM-D 17 changes. Each of the metabolite changes was negatively associated with changes in HAM-D 17 , but they were positively correlated to each other. α-aminoadipic acid showed the strongest correlation to HAM-D 17 changes (R 2 = -0.52, p < 7E-06). The BCAAderived C3-(propionyl) and C4-(butyryl) carnitines were strongly correlated to each other (R 2 = 0.93, p < 1.3E-11) as well as to α-aminoadipic acid. C3-carnitine was the hub metabolite in this module. Figure 4A presents a similar composite figure for the yellow module. All members were negatively correlated to HAM-D 17 changes, with the BCAA valine and the known oxidative stress biomarker methionine sulfoxide perturbations showing the strongest correlations with HAM-D 17 changes (R 2 < -0.40, p < 0.05). All the yellow module metabolites showed strong positive correlations with each other (mean R 2 = 0.75, range 0.44-0.93, all p < 0.02) in their perturbation patterns. Notable were the strong positive correlations between the branched chain amino acids and methionine, sarcosine, and the neurotransmitter related amino acids (phenylalanine, tryptophan and tyrosine) indicating that they may be functionally related. Valine was the hub metabolite for this module.
The members of the 96 lipids-containing green-yellow module also showed strong correlations amongst each other and also to HAM-D 17 changes. Unlike the amino acids and the short chain acylcarnitines from the purple and yellow modules, these lipid molecules' changes were positively associated with HAM-D 17 changes and also to each other (Supplementary Figures S1A,B).
In univariate linear mixed model analysis, the log2 fold changes of the α-aminoadipic acid, the branched chain amino acids, isoleucine and valine, methionine-sulfoxide and several of the phosphatidylcholines, containing long chain fatty acids, from the green-yellow module, were found to be significantly associated (p < 0.10) with the changes in HAM-D 17 scores FIGURE 2 | A heatmap showing correlations between the metabolite modules and the changes in symptom severity  scores. Each module represents a number of metabolites that showed strongly correlated perturbation patterns, from baseline to week 12, in response to CBT. The metabolite members of the three modules, purple, yellow and green-yellow, that showed significant association (R 2 > 0.3, at p < 0.1) with changes in HAM-D 17 scores are presented in adjacent boxes.
after adjusting for age and sex. These findings highlight the potential involvements of mitochondrial energy metabolism and lipids in the response of depressed patients receiving CBT. Table 3 depicts the list of metabolites that were significant in univariate models and their associations with depressive symptom scores (HAM-D 17 ).
To maximize our ability to detect the effects of treatment outcomes, we plotted the mean trajectories of the metabolites among the CBT remitters (N = 12) and the treatment failures (N = 7), leaving out the patients with intermediate outcomes, consistent with the approach used in other biomarker studies (Dunlop et al., 2017b;Vadodaria et al., 2019). There were interesting trends among the metabolites in the purple, yellow, and green-yellow modules. In the green-yellow module, consisting of lipids, ∼75% of the phosphatidylcholines were higher at baseline in the remitters compared to the treatment failures (Supplementary Figure S2), with some of them being statistically significant at p < 0.10 (PC aa-C30:0, C34:1, C36:2, C36:3, PC ae-C36:0, C38:2). In the purple ( Figure 3B) and yellow ( Figure 4B) modules, consisting mostly of amino acids including BCAAs and the short-chain acylcarnitines, the remitters mostly had lower baseline levels but showed an upward trend in their metabolite trajectories from baseline to week 12 while the treatment failures all trended to decrease to lower levels posttherapy. The limited sample size, however, precluded detection of statistical significance to these differences in trajectories.

DISCUSSION
Using liquid-chromatography coupled to mass-spectrometry analyses, we examined the biochemical changes that occurred in the plasma of depressed outpatients completing a course of CBT. Changes in several metabolite modules, containing primarily short-chain acylcarnitines and α-aminoadipic acid (purple module) as well as branched-chain and neurotransmitterrelated amino acids (yellow module) and lipids (green-yellow module), were significantly associated with changes in depressive symptom severity over the 12-weeks of CBT treatment. The metabolites within each module were highly correlated, and therefore it is likely that the similarity in their perturbations stemmed from their functional relatedness or being members of the same affected pathways.
Changes in individuals' metabotypes during periods of intense stress, and their return to the original homeostatic levels upon stress resolution (Ghini et al., 2015), support the possibility of identifying metabolomic state markers in MDD. Our analyses found specific amino acids, acylcarnitines, phosphatidylcholines, and sphingomyelins were associated with the depressed state and with changes after CBT treatment. Most notably, concentrations of the BCAAs isoleucine and valine, along with methionine sulfoxide and α-aminoadipic acid, showed strong inverse correlations with change in depression severity. Conversely, many lipid metabolites were directly correlated with changes in depression severity.
Alterations in BCAAs have been linked to altered mitochondrial energy metabolism and have been previously implicated in MDD (Baranyi et al., 2016) and in response to treatment (Kaddurah-Daouk et al., 2011. We have shown that MDD patients have alterations in the phenylalanine,  According to common lipid nomenclature, a LPC stands for lysophosphatidylcholine with a single fatty acid chain, b PC aa stands for Phosphatidylcholine diacyl (two fatty acid chains) and c PC ae stands for Phosphatidylcholine acyl-alkyl. The lipid species are described as CX:Y where X is the length of the carbon chain C, Y is the number of double bonds; "a" means the acyl chain is attached via an ester bond to the backbone while "e" means the attachment is via an ether bond.
tyrosine, and tryptophan pathways, which are involved in the biosynthesis of the monoamine neurotransmitters (Lucca et al., 1992;Maes et al., 1997;Bhattacharyya et al., 2019). Our metabolomic results overlap with the state metabolic markers identified in obesity, type 2 diabetes, and overall worsening metabolic health (Schooneman et al., 2013;Libert et al., 2018). Taken together, the patterns of change observed in our sample implicate bioenergetics as a focus of the pathobiology of the depressed state. The lipid perturbations, especially those of the phosphatidylcholines (consisting of either diacyl or alkylacyl moieties) showed positive correlations to the changes in HAM-D 17 scores. Phosphatidylcholines are a large class of lipid molecules commonly known as the glycerophospholipids. They have important functions in membrane stability, permeability, and signaling. Phosphatidylcholines have been implicated in MDD in several studies. Recently Knowles et al. (2017) suggested that a subclass of the phosphatidylcholines (the ether-phosphatidylcholines) might have a shared genetic etiology with MDD, and thus might be candidates for improved diagnosis and treatment of depression. Our previous work with ketamine have implicated these lipids along with acylcarnitines in the mechanism of response to ketamine (Rotroff et al., 2016).
These results, suggesting that the metabolites may serve as state markers of depression, received tentative support from our exploratory contrast of the differential trajectories of the changes in metabolites between the patients with the clearest treatment outcomes: remitters versus treatment failures. The BCAAs, their catabolic byproducts, the short chain acylcarnitines, the lysine metabolite α-aminoadipic acid, and the aromatic amino acids (phenylalanine, tyrosine and tryptophan) were present at comparatively higher levels at baseline in the treatment failures compared to the remitters. These findings suggest that metabolic wellbeing may be an important factor contributing to CBT response. Interestingly, there was a general downward trend in the trajectories of these metabolites over the course of treatment among the CBT treatment failures, whereas in the remitters they all exhibited stable or upward trajectories. It is possible that the differences in metabolic trajectories observed in remitters versus treatment failures to CBT indicate a state of "metabolic resilience" (Ghini et al., 2015) in the remitters. The two outcome groups also differed in their metabolomic profiles at baseline, with approximately 75% of the phosphatidylcholines being elevated in the remitters compared to the treatment failures, suggesting that levels of these lipid components may serve as moderators of outcome to CBT. It is also possible that these baseline differences reflect true subtype differences among the MDD patients that might have a role in treatment selection. Larger scale, longitudinal studies will be necessary to test these hypotheses.
The primary limitation to this study is the relatively small number of subjects analyzed. Consequently, testing for the statistical significance of the effects observed was limited, particularly for the categorical comparisons of differences by treatment outcome group; we could not analyze non-remitting responders and partial responders as additional groups. The study also lacked a healthy control group that could have permitted quantification of how far the observed metabolite concentrations were outside the "normal" ranges. Six of the subjects took a non-benzodiazepine sedative on an as-needed basis at some point during their CBT treatment. Although unlikely, we cannot rule out that this infrequent use could have impacted the metabolomic profiles. Despite these limitations, the study is important as an original exploration of the metabolomic changes in depression with a proven and welldelivered psychotherapy treatment, in the absence of the powerful metabolomic effects of psychopharmacotherapy.
In summary, this pilot evaluation assessed ∼180 metabolites from the Biocrates Absolute p180 kit that clustered into 8 "co-expression modules" based on their propensities to change over 12 weeks of treatment with CBT. The results were largely confirmed by additional univariate analyses of the individual metabolites in the co-expression analyses. Specifically, BCAAs, methionine sulfoxide, α-aminoadipic acid, and multiple phosphatidylcholines were all altered in association with changes in the HAM-D 17 scores at a level that equaled or exceeded a correlation coefficient of 0.4. Hence, these metabolites may represent markers for the depressed state, and perhaps may act as moderators or mediators for improvement from depression.These results will be useful in comparing and contrasting metabolomic changes that occur during and after treatment with antidepressant medication, and perhaps serving as a biomarker to inform treatment selection for MDD patients. In addition, metabolomic profiles in medication-free patients may prove to have utility as a biomarker for impending depressive relapse during long-term follow-up studies. Metabolomic profiles may also find a clinical application as surrogate measures for more expensive and cumbersome biomarkers, such as neuroimaging, if links between these biological systems can be identified (Dunlop and Mayberg, 2014). The results of this study provide an empirically testable set of hypotheses for MDD, namely the utility of these metabolites as potential state markers, in order to understand the mechanistic underpinnings of MDD and their change associated with symptomatic improvement.

DATA AVAILABILITY
The datasets generated for this study are available on request to the corresponding author.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Emory University Institutional Review Board. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
SB, BD, and AR conceived of the study and wrote the manuscript. SB and SM performed the data analysis. SB, BD, WC, AR, and RK-D assisted in biochemical and clinical interpretations of the results. All authors edited the manuscript.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnins. 2019.00926/full#supplementary-material FIGURE S1 | (A) A heatmap showing correlations between the metabolite modules and the changes in symptom severity (HAM-D17) scores, highlighting the greenyellow module. The member metabolites of this large module is presented in the adjacent box. (B) A heatmap showing the pairwise Pearson correlations between the member metabolites of the greenyellow module. The empty cells indicate that the p-values of the correlations for that metabolite-pair was not statistically significant (p > 0.1).
FIGURE S2 | The trajectories of each member metabolite of the greenyellow module, consisting of mostly lipids, or lipid-related metabolites, with its mean (± SEM) at baseline and week 12 are presented for remitters and treatment-failures.