Introduction

For children living in poverty, psychosocial stress is associated with dysregulation of biological stress responses.1 Downregulation of the hypothalamic-pituitary-adrenal (HPA) axis has been described as an adaptation to chronic stress.2,3 In our prior work with preschool-aged children,4 we found that children experiencing greater psychosocial stress had lower morning cortisol levels, consistent with downregulation of the HPA axis.2 Greater psychosocial stress was associated with greater mother-reported child food responsiveness, which in turn was associated with the child being overweight.4 Cortisol is known to have direct effects on metabolism affecting energy balance,5 suggesting the possibility that psychosocial stress may lead directly to obesity in young children through metabolic changes induced either directly by psychosocial stress or via changes in HPA axis functioning.

There is a growing interest in identifying the biological response to stress using metabolomics; the profile of small molecular weight compounds in a biological sample. Metabolite levels in tissues are influenced by dietary intake6 and can both reflect and influence changes in metabolism.7 Multiple studies have assessed how acute stress influences the metabolome in animals8,9 and humans10 by the mobilization of energy substrates via increasing lipolysis to prepare for the typical ā€œfight-or-flightā€ response. Very few studies have assessed how chronic psychosocial stress influences the metabolome.11,12 For instance, adults classified as ā€œstressedā€ had elevations in the metabolites 6-phosphogluconic acid, aminoadipic acid, d-arabitol, cysteine, sorbitol, d-fructose, threonic acid, 2-methylglytaric acid, chenodeoxycholic acid, and l-dihydroorotic acid.12 However, to our knowledge, the association between diurnal cortisol patterns paired with a self-classified stress level and the metabolome is unknown. Furthermore, classifying the metabolome at an early age may identify previously unrecognized biological mechanisms linking psychosocial stress exposure and health outcomes.

Our conceptual model is that psychosocial stress exposure alters the metabolome, which in turn contributes to obesity and cardiometabolic disease risk. Specifically, psychosocial stress exposure may shape the metabolome, via both direct mechanisms (i.e., stress exposure-induced hormonal changes influencing intermediary metabolism) and/or indirect mechanisms (i.e., changes in dietary patterns). In addition, the changes in metabolites associated with stress may shape obesity risk independent of HPA axis functioning, or they may occur within pathways from the HPA axis to obesity. Metabolite profiles associated with psychosocial stress may be related to obesity, which would suggest biological mechanisms linking stress and cardiometabolic health risk. Our objective in this study was to test the hypothesis that early childhood stress exposure, as indexed both by a biological marker (cortisol) and parent report (chaos in the home), is associated with a distinct metabolomic profile.

Methods

Study design

Participants were recruited from community agencies serving low-income families when children were 21 months old (2011ā€“2014). Inclusion criteria were that the biological mother was legal guardian, had an education level < a 4-year college degree, and was at least 18 years old; the family was eligible for Early Head Start, Women Infants and Children (WIC) Program, and/or Medicaid and was English-speaking; and the child was 21 months old, born at a gestational age ā‰„36 weeks, with no food allergies or significant health problems, perinatal or neonatal complications, or developmental delays. At 21 months, mothers completed questionnaires and children provided saliva samples for measurement of cortisol diurnal pattern and participated in anthropometry (initial visit). In 2015ā€“2017, families were invited to participate in a follow-up study where children participated in providing a fasting blood sample and anthropometry (follow-up visit). The research was approved by the University of Michigan Institutional Review Board. Written informed consent was provided by the mother and age-appropriate assent from children; families were compensated for their time. A total of 186 families participated at child age 21 months and of these 40 families enrolled in the follow-up study; 28 children provided blood samples and are included in the current analysis. At the initial visit, mothers reported childrenā€™s sex, race/ethnicity (categorized as non-Hispanic white vs. not), age, maternal education (<high school, high school diploma or generalized equivalency diploma, or some college), and family structure (single parent vs. not). The children who provided blood samples (nā€‰=ā€‰28) did not differ from those who did not (nā€‰=ā€‰12) nor the overall original sample (nā€‰=ā€‰186) in terms of child age, race/ethnicity, or sex. There was no difference in family CHAOS, cortisol intercept, and cortisol slope between the children who provided blood samples (nā€‰=ā€‰28) and the children who did not (nā€‰=ā€‰12).

Anthropometry

Child anthropometry (weight and length at the initial visit; weight, height, and waist circumference at the follow-up visit) was obtained by trained staff. All measurements were taken twice and the average was calculated for use in analyses. At the initial visit, weight-for-length (WFL) was calculated, and at the follow-up visit, body mass index (BMI) was calculated. Weight, length/height, and WFL or BMI z-scores were calculated based on the US Centers for Disease Control (CDC) reference growth curves.13

Psychosocial stress exposure

Psychosocial stress exposure was measured at the initial visit by parent report and using a biological marker. Obtaining biological markers of stress exposure is important because reported stress exposure is not always consistent with the physiological responses to the same stressors.14 Mothers completed the Confusion, Hubbub, and Order Scale (CHAOS),15 which consists of 15 statements (trueā€‰=ā€‰1, falseā€‰=ā€‰0) designed to assess chaos in the home environment as a psychosocial stress exposure index; scores range from 0 to 15 (higher scores reflect more chaos; Cronbachā€™s alphaā€‰=ā€‰0.81).

Chronic psychosocial stress exposure can be measured by assessing diurnal cortisol.16 Cortisol is reliably assessed through saliva.17 Salivary cortisol is highly correlated with free serum cortisol.18 At the initial visit, child saliva samples were collected by trained research staff; details have been described previously.19 The research assistant called the mom and asked where the child was (i.e., Is the child with you at home? Is the child at day care? Is the child still asleep?), and then based on the momā€™s response, went out to the childā€™s location and collected the saliva. Children provided saliva samples three times per day on 3 days within a week (weekdays only). The first sample each day was collected approximately 45ā€‰min after the child woke for the day, with the following samples occurring at 3 to 4-h intervals afterwards. There were 15 saliva samples missing out of an anticipated 252 samples (6%). Following collection, saliva samples were stored in Thermo Scientific Matrix Racks at āˆ’80ā€‰Ā°C until assayed. Cortisol was assayed using an Expanded Range High Sensitivity Salivary Cortisol Enzyme Immunoassay Kit (Catalog No. 1-3002, 96-Well Kit, Salimetrics LLC, PA, USA) with a detection limit of 0.007ā€‰Ī¼g/dL. Inter-assay coefficient of variation was 14% and intra-assay coefficient was 7.5%.

The generation of cortisol intercept and slope values in this sample has been described previously.19 Briefly, individual cortisol values were excluded if they were >3ā€‰SDs from the mean for the sample for that timepoint; to be included in the cortisol analysis, a child needed to have at least five valid cortisol data points over ā‰„2 days. Using methods described previously,4 hierarchical linear modeling generated random intercept and slope parameters to capture individual diurnal curves for each participant using the restricted maximum likelihood method.20 The diurnal curve follows a pattern such that cortisol increases initially after morning awakening, reaches a peak usually within 30ā€‰min, and afterward decays exponentially over the course of the day.21 A log transformation was implemented to capture the log-linear pattern (for time ā‰„60ā€‰min) of the diurnal cortisol to ensure normality for the residuals. Thus, using log transformed cortisol as the outcome and time (since awakening) at which the sampling occurred as the independent variable, the diurnal cortisol pattern is linear on time in a log-scale with the linear trajectory captured by two parameters, intercept and slope.

Untargeted metabolomics and shotgun lipidomics

At the follow-up visit, fasting plasma samples were collected for untargeted metabolomics and shotgun lipidomics. For untargeted metabolomics, samples were analyzed using a 1200 liquid chromatography/6530 qTOF MS instrument (Agilent Technologies, Inc., Santa Clara, CA USA) as previously described.22 The eluent was analyzed using positive and negative mode electrospray ionization. Untargeted shotgun lipidomics was performed as previously detailed23 using reconstituted lipid extract. Samples were ionized in positive and negative ionization model using a Triple TOF 5600 (ABSciex, Concord Canada). For both the untargeted metabolomics and lipidomics pipelines, pooled human plasma samples and pooled experimental samples were randomized and ran for quality control. Data normalization followed a recently described method.24 Detected features were excluded with (1) a relative standard deviation >40% in the pooled samples and/or (2) <70% presence in the samples. Missing data were imputed using the K-nearest neighbor method.25

For the untargeted metabolomics dataset, annotated metabolites were classified by comparing their retention time, ionization mass, and MS/MS spectra to external standards ran on the same instrument. Features were classified as ā€œredundantā€ and removed using ā€œBinnerā€; a new method to visualize closely eluting features allowing for the detection of adducts. Unannotated features are listed as ā€œexact mass @ retention time.ā€ For the shotgun lipidomics dataset, lipids were classified using LIPIDBLAST,26 a computer-generated tandem mass spectral library containing 119,200 compounds from 26 lipid classes. Raw feature data files with mass and retention time were searched against LIPIDBLAST using Multiquant 1.1.0.26 (ABsciex, Concord, Canada).27 All lipids mentioned have the nomenclature ā€œX:Yā€; where X is the length of the carbon change and Y is the number of double bonds. The final untargeted metabolomics dataset contained 513 annotated metabolites and 2022 unannotated features. The final lipidomics dataset contained 309 lipids from lipid classes including ceramides (CER), cholesterol esters (CE), diacylglycerols (DG), free fatty acids (FA), lysophosphatidylcholine (LysoPC), lysophosphatidylethanolamine (LysoPE), monoacylglycerol (MG), phosphatidic acid (PA), phosphatidylcholine (PC), phosphatidylethanolamine (PE), phosphatidylglycerol (PG), phosphatidylinositol (PI), phosphatidylserine (PS), sphingomyelin (SM), and triacylglycerol (TG).

Statistical analysis

TableĀ 1 summarizes the participantā€™s age, race/ethnicity, anthropometry, and chronic psychosocial stress exposure during the initial and follow-up visits. Linear models were created to test the association between chronic psychosocial stress exposure (CHAOS score, cortisol slope, and cortisol intercept) with the metabolome and lipidome, adjusting for sex and age at the follow-up visit. Individual linear models were performed for each metabolite identified, creating 2844 models. To limit Type 1 errors, associations were considered significant using a false discovery rate (FDR)ā€‰<ā€‰0.100;28 a commonly used method of adjusting for multiple comparisons in untargeted metabolomics studies. In certain lipid groups of biological interest, beta coefficients from the linear models were explored vs. chain length and number of double bonds in the lipid. The connectivity between the significantly associated annotated and unannotated metabolites was assessed using Pearsonā€™s correlations and clustered using hierarchical clustering. Using the National Institute of Standards and Technology (NIST) reference manual,29 groups of positively correlated unannotated metabolites were classified. To address our conceptual model that psychosocial stress exposure alters the metabolome, which in turn contributes to obesity and cardiometabolic disease, linear models were created to determine if metabolites associated with stress are associated with anthropometry at the follow-up visit, adjusting for sex and age at follow-up (unadjusted p value Ī±ā€‰=ā€‰0.05). Statistical analyses were performed using SAS 9.4 (Cary, North Carolina). Figures were created using GraphPad Prism version 7.4 (La Jolla, California). Heatmaps were created using an in-house software, CoolMap.

Table 1 ABC participant characteristics.
Full size table

Results

TableĀ 1 summarizes the participant characteristics. The sample consisted of 16 boys and 12 girls. A portion of the participants were from a single-family household (18%). The majority of participants were white, non-Hispanic (54%) and black, non-Hispanic (29%). The mean WFLz score was 1.17 at the initial visit and mean BMIz was 0.67 at the follow-up visit.

As an exploratory analysis, the linear relationship of chronic psychosocial stress exposure, measured by the CHAOS score and diurnal cortisol patterns, with the metabolome was quantified, adjusting for sex and age at the follow-up visit (Supplemental TableĀ S1 (online)). CHAOS score is associated with 11 metabolite features, cortisol intercept is associated with 8 metabolite features, and cortisol slope is associated with 20 metabolite features (FDRā€‰<ā€‰0.1). Of the significant metabolites, four unannotated features are associated with both cortisol intercept and slope (Fig.Ā 1a). FigureĀ 1b reports the direction of association between each chronic psychosocial stress exposure and the 35 significant features. The associations of the metabolome with cortisol intercept and slope tend to be in the same direction. Three annotated metabolites have a significant association with chronic psychosocial stress; cortisol slope is positively associated with dicarboxylic FA (FA-DiC) 21:0 (Ī²ā€‰=ā€‰26.3435, FDRā€‰=ā€‰0.0375) and inversely associated with FA 20:3 (Ī²ā€‰=ā€‰āˆ’29.5488, FDRā€‰=ā€‰0.0038) and FA 20:4 (Ī²ā€‰=ā€‰āˆ’23.8913, FDRā€‰=ā€‰0.0434), suggesting a relationship between diurnal cortisol patterns and fatty acid metabolism.

Fig. 1: Correlations between metabolite features associated with stress exposure.
figure 1

a Venn diagram displays small overlap between metabolite features associated with stress exposure (FDRā€‰<ā€‰0.10). b Features significantly associated with cortisol intercept, cortisol slope, and CHAOS score were correlated using Pearsonā€™s correlations and clustered using hierarchical clustering. Beta coefficients describing the relationship between stress exposure and the metabolome with blue denoting a positive association and yellow denoting a negative association, adjusting for sex and age at the follow-up visit. Significance denoted by ā€œ*ā€ (FDRā€‰<ā€‰0.1).

Full size image

Pearsonā€™s correlations were found between the annotated and unannotated metabolites significantly associated with psychosocial stress and correlations were clustered using hierarchical clustering to highlight the relationships between (Fig.Ā 1b, Supplemental TableĀ S2 (online)). Positive correlations were observed between the annotated metabolites FA 20:3 and FA 20:4 with 15 unannotated features, mainly from positive ionization mode. Using the NIST reference manual,29 we annotated three of the 15 highly correlated features (TableĀ 2). The three features were classified as very long-chain (VLC) hydroxyl- or epoxy-fatty acids with four, five, or six double bonds. One of the features is inversely associated with CHAOS score; 19,20-dihydroxy-FA 22:5 (Ī²ā€‰=ā€‰āˆ’0.1750, FDRā€‰=ā€‰0.0342). Two features are inversely associated with cortisol slope; 5(6)-Epoxy-FA 20:4, methyl ester (Ī²ā€‰=ā€‰āˆ’26.6333, FDRā€‰=ā€‰0.0076) and hydroxyl-FA 22:6 (Ī²ā€‰=ā€‰āˆ’26.0580, FDRā€‰=ā€‰0.0375). This suggests that these three metabolite features are intermediates in the VLC-PUFA metabolism, emphasizing the relationship between self-reported stress and blunted cortisol slope with PUFA levels. All potential annotations followed standard protocol for quality assurance in untargeted metabolomics studies.30

Table 2 Classification of unannotated features associated with early life stress exposure. Unannotated features were classified and scored by NIST software using mass and retention time from untargeted metabolomics analysis. "NIST identification" reports the potential identification of the feature with the "score" classifying the confidence of identification.
Full size table

Based on the above findings, we sought to identify patterns in the association between stress exposure with free fatty acids and triglycerides with fatty acid side chains of different degrees of saturation to further elucidate the relationship between stress and the lipidome. Beta coefficients describing the linear relationship between stress exposure and the metabolome (Supplementary TableĀ S1 (online)) were plotted against the number of double bonds within FA and TGs (Supplementary Fig.Ā S1 (online)). Unique trends were observed dependent on the number of double bonds and the chain length, in particular within TGs. For instance, we observed that long-chain PUFA within triglycerides tend to be inversely associated with CHAOS score, with several lipids reaching statistical significance (Supplementary Fig.Ā S1d (online)), represented by the wave-like fluctuations.

Lastly, we determined if the 35 metabolites significantly associated with stress exposure are moderately associated with measures of anthropometry concurrent with plasma draw. Using linear regression and adjusting for age and sex, two features were associated with BMIz and five features were associated with waist circumference (pā€‰<ā€‰0.05) (Supplementary TableĀ S3 (online)). The metabolite 19,20-dihydroxy-FA 22:5, which was inversely associated with CHAOS score, was positively associated with both BMIz (Ī²ā€‰=ā€‰0.3629, p valueā€‰=ā€‰0.0447) and waist circumference (Ī²ā€‰=ā€‰28.3160, p valueā€‰=ā€‰0.0040).

Discussion

The current exploratory analysis is the first examination of the association between chronic psychosocial stress exposure in early childhood and their subsequent metabolome in prepubertal children. Inverse associations were observed between both maternal-reported CHAOS score and diurnal cortisol patterns and VLC polyunsaturated fatty acids, including FA 20:3, FA 20:4, and hydroxy- and epoxy-PUFA metabolites. These associations suggest that a blunted cortisol slope, indicative of lower HPA activity across the day, is related to lower PUFA levels. Similar trends were observed with a lower cortisol intercept, which reflects morning levels. Further, results indicate that greater maternal-reported chronic psychosocial stress exposure is associated with lower child PUFA levels.

Polyunsaturated fatty acids are nutritionally essential fatty acids that are either ingested or formed from elongation and desaturation from linoleic acid (LA), the precursor for omega-6 fatty acids, and alpha-linolenic acid (ALA), the precursor for omega-3 fatty acids.31 LA is the most abundant PUFA within the Western diet, rich in soybean oil, corn oil, and sunflower oil, which are used in the production of most processed foods.32 Arachidonic acid (AA) is synthesized from LA, which is highly abundant in cell membranes and is a precursor for the production of eicosanoids, a class of signaling molecules.31 ALA is the precursor for eicosapentaenoic acid (EPA) and docosahexaenoic acid (DHA), which is highly abundant in vegetable sources, such as nuts and flax seed, and marine fish.32 EPA and DHA are involved in neural activity, such as cell membrane fluidity and ion transport.33 Omega-3 fatty acids have a large presence in the brain as 30% of the phosphoglycerides within the gray matter of the brain contain DHA, suggesting they are essential for neuronal functioning.34 Our results indicate that a blunted cortisol slope and a higher CHAOS score are associated with lower levels of plasma omega-3 fatty acids (TableĀ 2). These results may suggest lower levels of omega-3 PUFA within neuronal tissue, which could feed back to alter HPA regulation.

Of interest is the implication of PUFA deficiencies on neuronal functioning during critical periods of growth. Brain development varies by region, but overall, begins in early gestation and continues to mature through adulthood,35 highlighting the importance of maintaining fatty acid levels throughout growth. It is known that PUFA play a critical role in neurodevelopment during sensitive periods such as pregnancy and early childhood.36 The developmental origins of health and disease hypothesis suggests that insults during critical periods of development can permanently influence the development of the organism, via epigenetic processes, with long-term health implications. In early childhood, inadequate levels of omega-3 PUFA may alter cell membrane formation and cell signaling as DHA accumulates in the brain during the first 18 months of life.37 The negative impact of deficiencies in omega-3 FAs during development has been studied in animals, finding influences in behavior, learning, and cognition.35,38 The associations observed in this analysis may suggest the influence of deficiencies of PUFA on brain functioning.

Studies assessing the relationship between the HPA axis and PUFAs are focused in adults with psychiatric disorders, including major depressive disease.39,40 Higher evening cortisol levels (indicative of a blunted cortisol slope) are inversely associated with DHA, the FA unsaturation index, the FA chain length index, and the FA peroxidation index,40 suggesting a relationship between HPA activity and FA metabolism. The mechanisms of association remain poorly understood.

Chronic psychosocial stress exposure may contribute to dietary intake. Children tend to prefer more palatable and energy-dense snack foods, which are less healthy and higher in saturated and lower in PUFA.41 Our group has previously described associations of psychosocial stress exposure with greater consumption of processed foods, known as comfort foods.4 Eating these comfort foods may reduce feelings of stress via changes in the HPA axis.5 Potentially, children with more chronic psychosocial stress exposure may be fed foods characterized by lower levels of PUFA.

Cortisol levels have been reported to influence FA metabolism by increasing both the uptake and turnover of fatty acids in adipose tissue (reviewed in ref.ā€‰42). Supplementation of omega-3 PUFA could address potential impacts of psychosocial stress exposure on PUFA metabolism. Fish oil supplements in preterm infants increase levels of PUFAs in red blood cells.43 Future directions could consider altering dietary omega-3 intake to gauge if supplementation is a novel intervention for potential effects of psychosocial stress on metabolism.

Our study adds to the limited literature evaluating associations of chronic psychosocial stress exposure with the metabolome. Strengths of this study include the high-risk sample of young children recruited from the community with followed longitudinally, as well as the availability of both biological and maternal-reported measures of stress exposure. A limitation of this study is the small sample size and we emphasize the need for replication of these metabolomics results in a large study cohort. Our participantsā€™ race/ethnicities were mainly white, non-Hispanic (54%) and black, non-Hispanic (29%), potentially making our results not generalizable to other cohorts of children. Future studies should be focused on recruiting a larger sample population with a variety of race/ethnicities to assess if the relationship between psychosocial stress and the metabolome is exacerbated by race/ethnicity. A second limitation of this study is not having a psychosocial stress measurement at the time of blood draw for metabolomics to assess their concurrent relationship. Incorporating measurements of stress across multiple timepoints may decipher if early exposure to psychosocial stress programs the metabolome during childhood.

In conclusion, this exploratory analysis found some support for the conceptual model that psychosocial stress exposure alters the metabolome. Future work should engage larger sample sizes with variability in adiposity and psychosocial stress exposure, and examine possible mechanisms of association. In the ABC cohort, future work will involve profiling the untargeted metabolome and lipidome on a larger sample size to assess for reproducibility of our results in a subsample.