Sixty-nine overweight or obese youth between the ages of 9 to 17 years with currently untreated depressive symptoms were recruited for this study. All participants provided written assent and at least one parent or legal guardian provided written informed consent prior to all study procedures. This study was approved by Stanford University's Institutional Review Board. Participants were recruited from pediatric mood and weight control programs and community advertisements.

Youth were included if their body mass index (BMI) was at the 85th percentile or higher for their age and sex based on the Center for Disease Control and Prevention BMI calculator for children and teens (https://nccd.cdc.gov/dnpabmi/calculator.aspx). During the screening visit, height (with accuracy of 0.1 cm) and weight (with accuracy of 0.1 kg) were measured with the Seca 284, an electronic measuring station, after the removal of shoes and jackets. Two measures of height were obtained and averaged. We assessed Tanner stage using the youth self-reported Pubertal Development Scale (40) in conjunction with a clinician's physical examination of secondary sex characteristics to confirm the self-reported rating. Participants were also evaluated on levels of depression severity using the Children's Depression Rating Scale-Revised (CDRS-R) (41) administered separately to parents and youth by a study clinician or trained coordinator. Participants were included if their raw CDRS-R summary scores were greater than 35, signifying at least moderate levels of depression severity at the time of enrollment. All youth in this sample were treatment-seeking for active unremitted symptoms, but unmedicated, and were generally early in the course of their depressive illness.

Individuals were assessed for current and lifetime psychiatric disorders with the Kiddie Schedule for Affective Disorders and Schizophrenia-Present and Lifetime version (KSADS-PL) (42). Youth were excluded if they were already being treated for a mood disorder when evaluated at the screening visit. Youth were also excluded if they had type 1 or type 2 diabetes, were taking medication that affected their mood, weight, or metabolism at the time of screening, had a contraindication for an MRI (e.g., metal in their body or anterior-posterior diameter > 46 cm), or if their Full-4 IQ score on the Wechsler Abbreviated Scale of Intelligence (WASI) (43) was <70.

The Childhood Trauma Questionnaire (CTQ) was used as our primary clinical assessment of early childhood abuse. The CTQ is a 28-item scale of childhood trauma self-reported by youth participants (44). The CTQ separates trauma into five categorical scores with five questions per category: emotional abuse, physical abuse, sexual abuse, emotional neglect, and physical neglect. We secondarily used the Multidimensional Peer Victimization Scale (PVS) to understand other environmental forms of abuse, such as bullying, that are commonly experienced by youth outside of the home (45). This scale is a 16-item self-report measure with 4 subscales, to assess physical victimization, verbal victimization, social manipulation, and property attacks. We used a totaled PVS score in our linear models to adjust for the possible confounding effect of peer victimization on the relation between abuse and neural and endocrine outcomes, as childhood bullying has been shown to predict systemic inflammation and being overweight in early adulthood (46, 47). We used the Three-Factor Eating Questionnaire (TFEQ) (48) in youth to contexualize neural findings with behavioral measures of food-seeking behavior, including the constructs of cognitive restraint, uncontrolled eating or eating disinhibition, and emotional eating. Finally, we assessed socioeconomic status by the Hollingshead Four-Factor Index of Social Status and reported participants' raw scores (49).

After the screening visit, eligible youth were assessed for insulin resistance. Serum markers of insulin resistance were assessed using a 2-h oral glucose tolerance test (OGTT). After a 10-h fasting period and an initial fasting blood draw, participants consumed 75 g of oral glucose and had their blood drawn every 30 min for 2 h to measure insulin and glucose (37). Insulin values were measured by immunoassay.

Insulin response to OGTT was determined by calculating the area under the serum insulin curves, plotted at each time point during the OGTT. Graphing plasma glucose levels vs. time and finding the area under the curve similarly calculated glucose response to OGTT. Insulin resistance was determined by fasting insulin measure and by HOMA-IR, which was calculated using the equation fasting insulin (mU/mL) × fasting glucose (mg/dL)/405 (50).

Using the cumulative CTQ categories of emotional abuse, physical abuse, and sexual abuse, we grouped participants into “high abuse” and “low abuse” groups. Here, we focused on abuse rather than emotional or physical neglect, due to previous findings of childhood abuse as risk factors for metabolic outcomes such as prediabetes and diabetes (51, 52), and due to the demonstrated role of emotional abuse in predicting psychopathology (16). We classified participants as highly abused if they met the threshold for “low to moderate” abuse, as defined by the CTQ manual, for at least one of the three aforementioned categories. We used the following “high abuse” cutoff scores: an emotional abuse categorical score > 9; a physical abuse categorical score > 8; and a sexual abuse categorical score > 6. We also totaled the three abuse scores for a total CTQ abuse score that we used as a continuous variable in linear modeling. In order to be included in the “low abuse” group, youth could not endorse abuse above any of the cut-off scores for the three CTQ categories of abuse. Using such criteria, we were able to define with confidence a “low abuse” group of youth who did not report any significant early life stress from childhood abuse.

We used neuroimaging to investigate neural connectivity group differences at resting state, assessing for neural markers of approach motivation. Using functional MRI data of a subset of our participants (n = 47), we seeded three different regions that have been implicated in reward and depression: amygdala, insula, and nucleus accumbens (NAcc). After participants were familiarized with the scanning environment in an MRI simulator, whole-brain images were acquired on a 3T GE Signa Excite (General Electric Co., Milwaukee, WI) scanner equipped with an 8-channel head coil. Functional images were collected at rest using a spiral pulse sequence with the following parameters: repetition time (TR) = 2,000 ms, echo time (TE) = 30 ms, flip angle (FA) = 80°, field of view (FOV) = 22 cm, number of slices = 30 slices in the axial plane, and slice thickness = 4 mm with a gap of 1 mm. The first four volumes of each resting state scan were discarded at the scanner to allow for stabilization of longitudinal magnetization. High-order shimming was used before acquisition of resting state data to improve field inhomogeneity. High-resolution structural images were also collected to assist in registration of functional data to standard space using a fast spoiled gradient recalled (3D FSPGR) pulse sequence with the following parameters: TR = 8.5 ms, TE = 3.32 ms, TI = 400 ms, flip angle = 15°, field of view(x) = 25.6 cm, matrix of 256 × 256, number of slices = 186 slices in the axial plane, and a slice thickness of 1 mm.

Pre-processing of resting-state data was carried out using FEAT Version 6.00 within FSL (FMRIB's Software Library; www.fmrib.ox.ac.uk/fsl). The 210-volume functional dataset for each participant was realigned to compensate for small head movements using MCFLIRT (53), skull-stripped using the Brain Extraction Tool (BET) (54); spatially smoothed using a Gaussian kernel of 5 mm FWHM, intensity normalized by a single multiplicative factor, and band-pass filtered to correct for baseline drift and high frequency noise (high-pass temporal filter: Gaussian-weighted least-squares straight line fitting, with sigma = 50.0 s; low-pass temporal filter: Gaussian with sigma = 2.8 s). Functional images were registered to corresponding high-resolution T1-weighted structural images and then normalized to Montreal Neurological Institute (MNI) space using a 12-parameter transformation. Masks of white matter and cerebrospinal fluid (CSF) generated from anatomical images were applied to the functional data to extract white matter and CSF time-series. These time-series were used together with 6 motion parameters as nuisance regressors in a voxel-wise regression of the fMRI data. Data scrubbing was also performed following the method of Power et al. (55), excluding any volume in which either the value for DVARS (the root mean squared change in BOLD signal from the prior volume) or the value for framewise displacement exceeded the upper boxplot threshold (the 75th percentile plus 1.5 times the interquartile range), along with the previous volume and the 2 following volumes. Forty-Seven out of sixty-nine participants in the resting state analysis had < 33% of the volumes requiring removal, enabling inclusion in this analysis. Subjects included in the RSFC analysis did not differ from subjects excluded due to motion (P > 0.05 on all demographic, clinical, and OGTT-related variables). Similarly, high and low abuse groups did not differ in number of censored volumes [t₍₄₇₎ = 0.403, P = 0.689]. Residuals of the voxel-wise regression were used in subsequent seed-based connectivity analyses.

A seed-based intrinsic connectivity approach was used to examine resting state functional connectivity (RSFC) with the bilateral amygdala, NAcc, and insula in three separate whole brain analyses. Seed regions were anatomically defined using probabilistic maps from the Harvard-Oxford Subcortical Structural Atlas (https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/Atlases), incorporating voxels that had 25% or greater probability of being labeled as the amygdala (left: 3,456 mm³, right: 3,936 mm³), the NAcc (left: 1,088 mm³, right: 1,016 mm³) or the insula (left: 9,928 mm³, right: 10,080 mm³). Amygdala, NAcc, and insula ROIs were registered from the MNI template to each participant's preprocessed fMRI data, and the mean time-series of voxels in the right and left hemispheres for each of these regions were extracted for each participant for use as primary regressors in a GLM analysis of all other voxel time-series, resulting in individual whole-brain amygdala, NAcc, and insula RSFC maps. Talairach Daemon Labels (http://www.talairach.org) were used to find corresponding Brodmann areas.

Group differences in bilateral amygdala, bilateral NAcc, and bilateral insula RSFC were examined in separate voxel-wise t-tests, covarying age and sex. Subjects included in the RSFC analysis showed no difference in Full-4 IQ scores by abuse group (P > 0.05), so we did not include IQ as a covariate in the imaging subanalyses. Resulting statistical maps were thresholded with a height threshold of z > 2.3 and an extent threshold of P < 0.01667 (Bonferroni-corrected to control for multiple comparisons with three seeds), using Gaussian random field theory to correct for multiple comparisons. Parameter estimates (proportional to fMRI signal change) of BOLD signal response were extracted separately for each cluster and for each participant using featquery (fsl.fmrib.ox.ac.uk/fsl/fsl4.0/feat5/featquery.html) and analyzed in separate GLMs in SPSS (v.22; www.ibm.com/analytics/us/en/technology/spss/) that modeled the parameter estimate as the dependent variables and age, sex, and group as independent variables.

All statistical analyses were carried out using R version 3.5. t-Tests were used to test the association between abuse group and demographic, clinical, metabolic, and behavioral characteristics. Linear regression models were performed to test for association between abuse and OGTT-derived measures of insulin and glucose. Specifically, totaled CTQ abuse score was included in linear models as a predictor of area under the OGTT insulin curve, insulin at OGTT time points of 90 and 120 min, and glucose at OGTT 120 min. To account for various possible confounders and mediators, we included age, sex, BMI percentile, IQ, pubertal development (assessed by Tanner stage), depression severity, and total Peer Victimization Scale score as covariates in subsequent linear regression analyses.

Statistical assumptions for linear regression models were tested using regression diagnostics, including tests of normality of residuals, heteroscedasticity, linearity, and collinearity. While performing regression diagnostics, we found that insulin measures, including area under the OGTT insulin curve, were not normally distributed across all study participants. We performed log transformations and checked with the Shapiro-Wilk test that log-transformed area under the insulin curve met the normality assumption for linear regression. We also identified one subject as a statistical outlier in the association between totaled CTQ abuse score and glucose at OGTT endpoint. By reviewing Cook's distance and leverage statistics, we determined the subject's T+120 glucose level to exert undue influence on the regression, and we removed the subject from this and subsequent linear regression analyses with T+120 glucose.

Following our linear regression analysis, to determine if food-seeking behavior mediated the differences in OGTT insulin response between abuse groups, we used t-tests to investigate group differences in eating behavior.

Finally, with the connectivities that we identified with RSFC analysis as significantly different between abuse groups, we performed an interaction analysis between totaled abuse score and parameter estimates of BOLD signal response. These included (1) the interaction between abuse and amygdala-precuneus connectivity on fasting glucose, (2) the interaction between abuse and insula-precuneus connectivity on area under insulin curve, (3) the interaction between NAcc-prefrontal connectivity on fasting glucose, and (4) the interaction between abuse and NAcc-paracingulate gyrus connectivity on T+120 glucose.

Demographic and clinical characteristics of the 69 participants are summarized in Table 1. Fifty-nine percent of the youth in this sample were female, with diverse ethnic backgrounds (55% Caucasian), with above average IQ (102.46 ± 14.14), with an average BMI in the obese range (29.95 ± 6.17), and with moderate to severe depression severity (CDRS-R Raw Score = 54.19 ± 11.28). By grouping the participants by high or low levels of abuse (see Methods), we noted significant differences in age, IQ, and Tanner stage, so we adjusted for these variables in all subsequent linear modeling. We also examined how OGTT insulin and glucose measures varied by participants' level of abuse (Table 2). Participants in the high abuse group had significantly lower insulin levels at the 2-h time point (P = 0.015) and almost significantly lower areas under the insulin curve (P = 0.083). Glucose levels at the 2-h time point were also significantly different between the two abuse groups (P = 0.049), but the difference between abuse groups in areas under the glucose curve was not significant (P = 0.280). We found no mediating effect of food-seeking behavior in our analysis, as the results showed no significant group differences in self-reported eating behaviors or significant associations with OGTT measures.

In order to evaluate differences in insulin and glucose response to OGTT between the two abuse groups, we visualized each abuse group's mean serum insulin and mean plasma glucose at time point during the oral glucose tolerance test, creating metabolic response curves. A more elevated insulin response to the oral glucose challenge was observed in participants who reported a low level of abuse (Figure 1A). Insulin and glucose shapes varied between groups; mean insulin curve in the high abuse group showed earlier peak and decline, whereas the curve for the low abuse group was elevated for a sustained period. Mean glucose curves showed more similarity between the two groups (Figure 1B). Peaks for both curves were observed at the 30-min time point, followed by declines in plasma glucose.

Linear regression modeling showed that insulin and glucose response was associated with level of abuse (Table 3). Totaled CTQ abuse scores significantly correlated with log-transformed area under the insulin curve, insulin at 90 min, and glucose at OGTT endpoint. We found that age, sex, race, IQ, and depression severity did not confound or mediate the relationship between abuse and OGTT measures. On the other hand, BMI percentile and Peer Victimization Scale (PVS) total score were identified as possible negative confounders of the association between abuse and log-transformed insulin curve area (Δß = +12.5% adjusting for BMI; +29% adjusting for PVS), as well as between abuse and log-transformed T+90 insulin (Δß = +68.5% adjusting for PVS). In other words, adjusting for BMI percentile and peer victimization accentuated the significant effect of abuse on insulin response. Whereas we found significant differences in T+120 insulin between abuse groups when we treated abuse as a categorical variable, our linear modeling with totaled CTQ abuse as a continuous variable did not show correlations with insulin at OGTT endpoint.

Group differences in bilateral amygdala, bilateral NAcc, and bilateral insula RSFC were examined in separate voxel-wise t-tests, co-varying age, and sex. After correcting for multiple comparisons, we found a significant main effect of group between the bilateral amygdala and precuneus (k = 644 voxels, peak x/y/z MNI coordinate = 14/−58/24, z = 3.56, P = 0.000594, Brodmann areas 7, 23, 31) (Figure 2A). We also found a significant main effect of group between the bilateral insula and precuneus (k = 566 voxels, peak x/y/z MNI coordinate = 10/−66/34, z = 4.05, P = 0.002, Brodmann areas 7, 31) (Figure 2B). For both of these main effects, the high abuse group showed reduced negative connectivity to the precuneus compared to the low abuse group.

We also found significant group differences in connectivity between bilateral NAcc and two regions: a prefrontal cluster which included the left dorsolateral prefrontal cortex (DLPFC), precentral gyrus, and inferior frontal gyrus (IFG) (k = 597 voxels, peak x/y/z MNI coordinates = −54/4/30, z = 3.46, P = 0.0013, Brodmann areas 6, 9, 45, 46) and a cluster encompassing the bilateral paracingulate gyri in the medial prefrontal cortex (mPFC) (k = 386 voxels, peak x/y/z MNI coordinates = −12/48/−4, z = 3.89, P = 0.023, Brodmann areas 10, 32) (Figures 2C,D). In the high- compared to low- abuse group, there was increased negative connectivity between NAcc and the prefrontal cluster, and reduced positive connectivity between NAcc and the paracingulate gyrus cluster.

We extracted parameter estimates of BOLD signal response for these four connectivity clusters showing abuse group differences. Taking total abuse score as a continuous variable, we explored a dose relation between level of abuse and connectivity values that was significant for insula-precuneus (P = 0.0032, r = 0.421), NAcc-prefrontal (P = 0.0035, r = −0.418), and NAcc-paracingulate gyrus (P = 0.0046, r = −0.406) connectivities, and near significant for amygdala-precuneus connectivity (P = 0.064, r = 0.272).

With extracted neural connectivity values from the four connectivity group differences shown above, we explored associations between connectivity estimates and OGTT measures of insulin and glucose. We made scatterplots with least-squares lines for each abuse group and found that the two groups had visibly different correlations between the insulin or glucose response and brain connectivity (Figure 3). The high abuse group showed correlations between OGTT measures and amygdala-precuneus, NAcc-prefrontal cortex, and NAcc-paracingulate gyrus connectivity, while the low abuse group did not (Figures 3A,C,D). On the other hand, the low abuse group showed a negative correlation between area under the insulin curve and insula-precuneus connectivity, while the high abuse group did not (Figure 3B). Adjusting for age, sex, BMI, and peer victimization, we found a significant negative correlation across all participants between log-transformed area under the insulin curve and bilateral insula-precuneus connectivity (P = 0.041), as well as between log-transformed T+120 insulin and bilateral amygdala-precuneus connectivity (P = 0.019).

Mediation analysis using connectivity values and abuse scores as covariates in linear models of OGTT measures, both study-wide and by separated abuse group, did not show significant mediation by connectivity of the effect of abuse on insulin or glucose. However, when we performed moderation analysis to see if abuse moderated the association between connectivity estimates and OGTT measures, we found that each of the four identified connectivity differences interacted with abuse in predicting OGTT measures, in line with our scatterplot findings. Cumulative CTQ abuse score as a continuous variable moderated the associations between: amygdala-precuneus connectivity and fasting glucose (interaction term P = 0.043), insula-precuneus connectivity and log-transformed area under the insulin curve (interaction term P = 0.037), and NAcc-prefrontal connectivity and fasting glucose (interaction term P = 0.0005) (Tables 4–6). We also found moderation by cumulative CTQ abuse score of the association between NAcc-paracingulate gyrus connectivity and T+120 glucose, trending toward significance (interaction term P = 0.114) (Table 7).

MS receives research support from Stanford's Maternal Child Health Research Institute, National Insitute of Mental Health, National Institute of Aging, Neuronetics, Johnson and Johnson, and the Brain and Behavior Foundation. She is on the advisory board for Sunovion and Google X. NR receives research support from the National Institute of Aging. The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.