Participants
As part of a larger study on the relationship between stress and hormone function in MDD, adults with Major Depressive Disorder and healthy controls without current psychiatric disorders 21 and 45 years of age and with a BMI between 19 and 45 kg/m2 were recruited from the greater Boston area through online advertisements. From the larger sample of participants who completed all study components, RNA samples for the current analyses were collected on a subset of 21 individuals with Hyperphagic and Hypophagic MDD (8 Hyperphagic MDD [5 female]; 13 Hypophagic MDD [7 female]). Exclusion criteria included the history of substance abuse; current comorbid psychiatric disorders; current psychotropic medications use; history of neurological disease; endocrine disorders; diabetes; cardiovascular disease; current treatment with weight loss medications; glucocorticoids; steroids; current suicidal ideation; traumatic brain injury; for females, pregnancy or breastfeeding, and past amenorrhea greater than three months. All participants with MDD met diagnostic criteria as assessed by the Structured Clinical Interview for Diagnoses (SCID; DSM-IV-TR)11 administered by a trained clinician interviewer with over 20 years of experience. Participants with MDD were assigned into Hyperphagic MDD (e.g., increased appetite/weight gain) and Hypophagic MDD (e.g., decreased appetite/weight loss) groups based on their responses to appetite- and weight-related assessment in the mood disorders module of the SCID. Specifically, MDD group was assigned based on responses to criterion Question A3 under the “A. Mood Episodes” module: “(3) significant weight loss when not dieting, or weight gain (e.g., a change of more than 5% body weight in a month) or decrease or increase in appetite nearly every day for the past 2 weeks”. Participants who endorsed either weight loss OR decrease in appetite nearly every day for 2 weeks (A4) were categorized as Hypophagic MDD. Participants who endorsed either weight gain OR increase in appetite nearly every day for 2 weeks over the past month (A5) were categorized as Hyperphagic MDD. All participants provided written informed consent. All study procedures were approved by the Mass General Brigham Institutional Review Board.
Procedures
Following a screening session to establish MDD status and obtain height, weight, anthropometric measurements, and hematocrit level (< 34% in females and < 37% in males), eligible participants underwent two main study visits (1 week apart) that were identical except one included a stress-inducing version (Stress) and the other included a non-stressful control version (No Stress) of the Maastricht Acute Stress Task (MAST)12. Visit order (Stress, No Stress) was randomized and counterbalanced across groups. Visits occurred between 0800 hrs and 1300 hrs following a 12-hour overnight fast. For females, all visits occurred during the follicular phase of the menstrual cycle. Upon arrival, participants completed the Beck Depression Inventory (BDI-II) to assess current depressive symptoms13; female participants underwent a pregnancy test. Next, a nurse placed an intravenous catheter (IV) with a saline lock for blood sampling. After the IV was placed, participants consumed a breakfast meal (within a 15-minute period to standardize rate of intake) designed to comprise 30% of their recommended daily caloric intake (see detailed description below). The T0 blood draw was completed immediately following the breakfast session. At this time, participants completed an assessment of their mood using an electronic visual analog scale (VAS) system14.
Participants then completed either the Stress or No-Stress version of the MAST. Participants were seated in front of a computer, and a water bath was positioned on the side of the arm that did not have the IV. Participants were then introduced to a female experimenter posing as a doctor who provided detailed instructions for the water and math task. This was followed by 10 minutes during which the participants alternated between five hand immersion trials (60 sec to 90 sec) and four arithmetic trials (30 sec to 90 sec). During the task, trial switching was prompted by a timed slideshow running on the computer (Fig. 1).
During the Stress visit, hand immersion trials required the participant to place their hand and wrist in cold water that was held between 0°C and 2°C. Arithmetic trials required counting backward as quickly and accurately as possible from 2,043 in intervals of 17. If a mistake was made, the experimenter instructed the participant to start again from 2,043. Participants were told their performance was being videotaped by a camera that was mounted to the computer screen. In reality, the camera was not recording. Following the final hand immersion trial, participants were told that their performance was poor and that they would need to repeat the task later during the visit; this manipulation was used to induce sustained levels of stress throughout the visit.
During the No-Stress visit, the water temperature was lukewarm (25°C – 37°C), and the math trials required counting up consecutively from 1 to 25 in steps of 1 at their own pace and starting over once they reached 25. There was no mention of videotaping, and the experimenter gave no feedback on their performance. The study staff member playing the role as doctor was consistent for each participant across Stress and No-Stress visits.
Approximately 85 minutes after the MAST (105 minutes since T0), a second sample of blood was drawn while participants completed a second mood assessment via the VAS system. At the end of each visit, participants were fully debriefed. During the No-Stress visit, participants completed a questionnaire packet including the trait portion of the State-Trait Anxiety Inventory (STAI)15, Perceived Stress Scale (PSS)16, and the Dutch Eating Behavior Questionnaire (DEBQ)17.
Anthropometry
At the screening visit, height (cm) and weight (kg) were measured using a stadiometer. Those values were used to calculate body mass index (BMI [kg/m2]).
Breakfast meal
At the start of each visit, participants ate a breakfast meal that was standardized to meet 30% of their daily caloric intake. Breakfast was prepared by Center for Clinical Investigation (CCI) dietary staff at Brigham and Women’s Hospital and varied according to each participant’s basal metabolic rate and physical activity level, measured by the Harris-Benedict equation with 18% calories from protein, 23% calories from fat, and 59% calories from carbohydrates18. Participants were encouraged to consume as much of the meal that they could in 15 minutes. After the 15-minute period, CCI dietary staff weighed the remaining food items to determine nutrient intake for total kcals and kcals per macronutrient. Planned versus consumed (PVC) kcals for total breakfast meal and then each individual macronutrient (protein, fat, and carbohydrate) was calculated by dividing the intake values by the planned standardized amounts.
Blood sampling, hormone measurement, and RNA processing
Blood was sampled at three time points at each visit. Serum cortisol, acylated ghrelin, and RNA were collected at two time points: T0: After the breakfast meal/before the MAST and T105: 105 minutes after T0. Serum cortisol, alone, was collected at T20 as a manipulation check for the effect of the MAST.
For cortisol, after transferring to a tube containing pefabloc, serum was frozen until assayed. Serum samples for cortisol and ghrelin were stored at − 80°C in plastic tubes containing a 10-mg/ml solution of PMSF (phenylmethanesulfonylfluoride) in methanol. Cortisol samples were assayed by LabCorp (Raritan, NJ) using electrochemiluminescence immunoassay (ECLIA) on a Roche Cobas analyzer; intra-assay CV 1.0-1.7%; inter-assay CV 1.4–2.2%.
Acylated ghrelin was collected in Ethylenediaminetetraacetic acid (EDTA) tubes, aliquoted with HCL, and then centrifuged. Ghrelin samples were assayed by the Brigham Research Assay Core (BRAC) using an enzyme immunoassay (ELISA; Millipore, St. Charles, MO; intra-assay CV 0.8–7.5%; inter-assay CV 3.9 = 12.9%).
Total RNA was collected in Tempus Blood RNA tubes (Thermo Fisher Scientific, Waltham, MA) incubated at room temperature for 2 hours and stored at − 80°C according to the instructions of the manufacturer until further processing. Total RNA was extracted using the Norgen Total RNA Purification Kit (Cat. 17200, Norgen Biotek Corp, Thorold, ON, Canada) and quantified using a Qubit 2.0 Fluorometer (Thermo Fisher Scientific Inc., Waltham, MA). All samples were run on an Agilent 2100 Bioanalyzer (Agilent, Technologies, Santa Clara, CA) to determine RNA Integrity Number (RIN) and all samples passed QC with RIN values > 6 (mean = 8.81, SD = 0.46).
Behavioral measures
Beck Depression Inventory (BDI-II)– The BDI-II is a 21-item self-report questionnaire used to assess presence and severity of depressive symptoms13. Each item presents a characteristic (e.g., ‘Self-Dislike’) with four statements below it (labeled 0, 1, 2, and 3), with 0 representing minimal change in/severity of, and 3 representing maximal change in/severity of that characteristic. Participants were asked to indicate to what extent they experienced that characteristic in the past two weeks by choosing the statement they agreed with. Scores were calculated by summing the answers to each question. Higher scores indicate more severe depressive symptoms.
Mood Assessment – Ratings of sadness and tension were measured during the T0 and T105 blood draws using an electronic visual analog scale (VAS) system14. Participants were asked to rate how sad and how tense they felt by moving a slider on a linear scale that went from 0 to 100, with 0 being the least sad/tense they have ever felt and 100 being the saddest/tensest they have ever felt.
Dutch Eating Behavior Questionnaire (DEBQ) – During the No-Stress visit, participants completed the DEBQ, a 33-item self-report questionnaire that assesses three eating behavior domains: Restrained Eating (10 items), Emotional Eating (13 items), and External Eating (10 items)17. All items are rated on a 5-point scale with responses that range from 1 (‘Never’) to 5 (‘Very Often’). Scores for each subscale were calculated by adding up the item responses within a subscale and then dividing by the number of questions in that subscale resulting in a final score per subscale. Higher scores indicate greater tendency to display the subscale behavior.
State-Trait Anxiety Inventory-Trait (STAI-T) – Trait-based characteristics relating to general anxiety level were assessed using the STAI-T15. The STAI is comprised of 20 item that ask participants to rate level of proneness to anxiety (e.g., general states of calmness, security, or inadequacy) through Likert scale responses that range from 1 (‘Almost Never’) to 4 (‘Almost Always’). Scoring was reversed for the anxiety-absent items, which comprised 9 of the 20 questions. Total scores were calculated by adding the responses after reverse-scoring. The range of scores is from 20–80, with higher scores indicating greater general anxiety.
Perceived Stress Scale (PSS) – Participants completed the PSS as part of the post-visit questionnaire packet. The PSS is a 10-item questionnaire that assesses how often respondents have felt a particular way in relation to stressful situations within the past month (e.g., ‘In the last month, how often have you been upset because of something that happened unexpectedly?’)16. Each item is rated on a 5-point Likert-scale, with 0 being never and 4 being very often. Six of the items indicating last of perceived stress were reverse-scored. Total scores were calculated by summing the answers to all the scored and reverse-scored items. Higher scored correspond with higher perceived stress.
Statistical Methods
Demographic data, questionnaires, and hormone levels were analyzed using SPSS version 28 (IBM Corp., 2021). The p-value threshold for significance was 0.05. Based on small sample sizes and non-normal distributions, non-parametric statistical tests were employed, and serum cortisol and acylated ghrelin were log2 transformed. Demographics, baseline characteristics, BDI-II, STAI-T, DEBQ, PSS, and breakfast meal values (PVC total, PVC protein, PVC fat, and PVC carbohydrate) were assessed using Fisher’s Exact Tests, χ2, and Mann-Whitney U tests. Within group comparisons in hormone levels and VAS ratings of sadness and tension were assessed using Wilcoxon Signed Rank tests. Between-subjects effects in breakfast intake, cortisol and acylated ghrelin samples, and VAS ratings of sadness and tension were analyzed using Mann-Whitney U tests.
Bulk RNA-seq Analysis
Total RNA was sent to Genewiz (Azenta Life Sciences, Burlington, MA) for bulk-RNAseq. All samples were treated with TURBO DNase (Thermo Fisher Scientific, Waltham, MA) to remove residual DNA contaminants. Ribosomal RNA and globin RNA was depleted using QIAseq FastSelect − rRNA HMR and − Globin kit (Qiagen, Germantown, MD). RNA sequencing libraries were constructed with the NEBNext Ultra II RNA Library Preparation Kit for Illumina. The sequencing libraries were multiplexed and clustered on seven lanes of a flowcell. After clustering, the flowcell was loaded on the Illumina HiSeq 4000 instrument according to manufacturer’s instructions. The samples were sequenced using a 2x150 Pair-End (PE) configuration. Image analysis and base calling were conducted by the HiSeq Control Software (HCS). Raw sequence data (.bcl files) generated from Illumina HiSeq was converted into FASTQ files and de-multiplexed using Illumina's bcl2fastq 2.17 software. One mismatch was allowed for index sequence identification. The obtained RNAseq data was processed through the bulk RNA-seq pipeline in bcbio (v1.2.8)19. Reads were aligned using STAR against transcriptome references of the human genome (GRCh38 Ensembl release 94)20 and quantified with Salmon21. Two samples were excluded in further analysis due to poor sequencing quality.
Differential Gene Expression Analysis
Differential gene expression analysis was performed using limma (v3.46.0)22. Counts were filtered using limma::filterByExpr with default settings. Counts from both time points were normalized with limma::voom and repeated measures from same individual were accounted for using limma::duplicateCorrelation. To determine which covariates to include in the final analysis, association tests were performed between all known variables (group [Hyperphagic MDD, Hypophagic MDD], visit [Stress, No-Stress], sequencing quality metrics [percent_gc, r_rna, r_rna_rate, intergenic_rate, intronic_rate, exonic_rate, duplication_rate_of_mapped, average_insert_size, total_reads, mapped_reads, mapped_paired_reads, duplicates_pct], individual, sex, age, BMI, visit, breakfast intake [bfast_pvc_TOT, bfast_pvcPROT, bfast_pvcFAT, bfast_pvcCARB], RIN, RNA_extraction_batch) and the 12 calculated principal components (PCs) for each timepoint separately. Some sequencing metrics, RIN, RNA extraction batch, and breakfast intake tested, did not show a significant association with the top PCs (Supplemental Fig. 1). However, sequencing quality related metrics (intergenic_rate, intronic_rate, exonic_rate, duplication_rate_of_mapped, average_insert_size, and duplicates_pct) showed a significant association with the variation in our data as captured in the top PC and therefore exonic rate was included in the downstream models (pBonferroni < 0.05). Other sequencing related metrics were not included due to their correlation with exonic rate to prevent collinearity in the model.
At baseline T0, the gene expression differences between Hyperphagic MDD and Hypophagic MDD groups were compared using a linear mixed model controlling for age, BMI, sex, calories of breakfast intake, sequencing quality, and repeated measures. At T105, the gene expression differences between Hyperphagic MDD and Hypophagic MDD groups under stress were compared using a linear mixed model with group and visit interaction terms controlling for age, BMI, sex, calories of breakfast intake, sequencing quality, and repeated measures. Nominal p-values from both models were adjusted for inflation with bacon (v1.18.0)23 and differentially expressed genes were considered significant at pFDR<0.05.
Normalized counts of all significant genes identified at both time points were tested for correlation with corresponding cortisol and ghrelin levels, breakfast PVC values and mood ratings in linear regression models controlling for sex and repeated measures. Normalized counts at T0 were then averaged between visits to test for association with corresponding DEBQ subscale scores in linear regression models controlling for sex.
Enrichment Analysis
Differentially expressed genes (DEGs) identified at nominal significance level were used in the enrichment analyses. Overrepresentation tests were used to determine the enrichment of Gene Ontology (GO) terms24,25. All genes were ranked based on log fold change and used as input for Gene Set Enrichment Analysis (GSEA) to assess the level of enrichment of 161 inflammatory pathways identified by de Kluiver et al10 that are available from the Human MSigDB Gene Set26,27. The pathway “Reactome RIGI MDA5 mediated induction of IFN alpha beta pathway” was not found in the Human MSigDB Gene Sets database and excluded in the GSEA analysis.
Data and Code Availability
All RNA-seq statistical analyses were performed in R software version 4.0.2. Bulk RNA-seq data can be accessed through GEO Accession GSE231347. Code used in the analysis has been deposited on Github: https://github.com/klengellab/HypoHyperMDD .