Gastrointestinal and mental health problems are highly comorbid. For example, anxiety is five times higher in individuals with irritable bowel syndrome (IBS) than in those with no IBS symptoms, and among anxiety sufferers, rates of IBS (and other functional gastrointestinal disorders) are increased fourfold (Lee et al., Reference Lee, Wu, Ma, Tsang, Guo and Sung2009; Mak et al., Reference Mak, Wu, Chan, Chan, Sung and Lee2012), possibly reflecting a mechanistic association. Several studies in adults have shown support for a bidirectional relationship between mental illness and gastrointestinal disorders such as IBS. For example, in a population of healthcare seekers, anxiety and depression diagnoses preceded IBS diagnoses by approximately 3.5 years (Jones et al., Reference Jones, Tack, Van Oudenhove, Walker, Holtmann, Koloski and Talley2017), and experimental studies have shown that psychosocial stress can increase IBS symptoms (Williams, Villar, Peterson, & Burks, Reference Williams, Villar, Peterson and Burks1988), supporting the pathway of psychological to gastrointestinal (GI) dysfunction. The opposite pathway of influence (GI to psychological dysfunction) also has strong support in the literature. For example, in a prospective cohort study, IBS was shown to precede mood disturbances in over 60% of cases, whereas mood disturbance preceeded IBS in only 30% of cases (Koloski, Jones, & Talley, Reference Koloski, Jones and Talley2016). Such bidirectional associations may be indicative of a causal relationship between each symptom group, and/or may be reflective of a shared etiology.
Research has consistently shown that exposure to early life adversity is a potent risk factor for both GI and mental illnesses. Rats exposed to parental-deprivation stress exhibit increased IBS symptoms (i.e., visceral hypersensitivity) as well as elevated levels of anxiety at postweaning age, and into adulthood, compared to nondeprived peers (Yi et al., Reference Yi, Zhang, Sun, Zhou, Chen, Xuan and Xu2017). In humans, early adversity (particularly in the late postnatal period) is associated with up to a threefold increase in risk of irritable bowel symptoms (Bradford et al., Reference Bradford, Shih, Videlock, Presson, Naliboff, Mayer and Chang2012; Chitkara, van Tilburg, Blois-Martin, & Whitehead, Reference Chitkara, van Tilburg, Blois-Martin and Whitehead2008; Klooker et al., Reference Klooker, Braak, Painter, de Rooij, van Elburg, van den Wijngaard and Boeckxstaens2009; Park et al., Reference Park, Videlock, Shih, Presson, Mayer and Chang2016), and contributes to over a third of lifetime mental illness diagnoses (Green et al., Reference Green, McLaughlin, Berglund, Gruber, Sampson, Zaslavsky and Kessler2010; McLaughlin et al., Reference McLaughlin, Green, Gruber, Sampson, Zaslavsky and Kessler2010; Repetti, Taylor, & Seeman, Reference Repetti, Taylor and Seeman2002), as well as subthreshold and clinical levels of anxiety in children (Goff et al., Reference Goff, Gee, Telzer, Humphreys, Gabard-Durnam, Flannery and Tottenham2013; Wiik et al., Reference Wiik, Loman, Van Ryzin, Armstrong, Essex, Pollak and Gunnar2011). Such adversities also impact neurobiology implicated in affective responding and irritable bowel symptoms. In the rodent, early caregiving adversities alter prefrontal, amygdala, and hippocampal development (Roceri et al., Reference Roceri, Cirulli, Pessina, Peretto, Racagni and Riva2004), change cognitive and emotional behaviors that are supported by those neural structures (Callaghan & Richardson, Reference Callaghan and Richardson2012, Reference Callaghan and Richardson2014; Cowan, Callaghan, & Richardson, Reference Cowan, Callaghan and Richardson2013), as well as lead to increased visceral pain responses and elevated activation of the prefrontal cortex and amygdala when experiencing such pain (Felice et al., Reference Felice, Gibney, Gosselin, Dinan, O'Mahony and Cryan2014). Similarly, in humans, caregiver neglect, or its more potent form, early institutional care (which is characterized by parental deprivation, and is a significant risk factor for anxiety symptoms), is followed by a heightened risk for atypical prefrontal cortex and amygdala functional and structural development, prefrontal cortex connectivity, amygdala reactivity, and changes in default resting state networks (Bick et al., Reference Bick, Zhu, Stamoulis, Fox, Zeanah and Nelson2015; Callaghan et al., Reference Callaghan, Gee, Gabard-Durnam, Telzer, Humphreys, Goff and Tottenhamin press; Gee, Gabard-Durnam, et al., Reference Gee, Gabard-Durnam, Flannery, Goff, Humphreys, Telzer and Tottenham2013; Hanson et al., Reference Hanson, Adluru, Chung, Alexander, Davidson and Pollak2013; Hodel et al., Reference Hodel, Hunt, Cowell, Van Den Heuvel, Gunnar and Thomas2015; Maheu et al., Reference Maheu, Dozier, Guyer, Mandell, Peloso, Poeth and Ernst2010; Silvers et al., Reference Silvers, Lumian, Gabard-Durnam, Gee, Goff, Fareri and Tottenham2016; Tottenham et al., Reference Tottenham, Hare, Quinn, McCarry, Nurse, Gilhooly and Eigsti2010; Van der Werff et al., Reference Van der Werff, Pannekoek, Veer, van Tol, Aleman, Veltman and Van der Wee2013). Those neural targets overlap with the hubs of dysregulated neural activity seen in adult IBS sufferers (Mertz et al., Reference Mertz, Morgan, Tanner, Pickens, Price, Shyr and Kessler2000).
While clear links have been drawn between gastrointestinal and psychiatric symptoms across the life span (Garber, Zeman, & Walker, Reference Garber, Zeman and Walker1990; Rajindrajith et al., Reference Rajindrajith, Devanarayana, Lakmini, Subasinghe, de Silva and Benninga2014; Shelby et al., Reference Shelby, Shirkey, Sherman, Beck, Haman, Shears and Walker2013), the role of early adversity in increasing vulnerability to functional gastrointestinal complaints is only established in adults, with examinations in childhood almost nonexistent (but see Rajindrajith et al., Reference Rajindrajith, Devanarayana, Lakmini, Subasinghe, de Silva and Benninga2014). Hence it is unclear when in life such associations emerge. This is a critical gap in the literature, as early detection of such associations could lead to early intervention. In a previous study, approximately 80% of children and adolescents (5–18 years) who were presenting to a primary care physician with recurrent abdominal distress (which includes gastrointestinal disturbances) had comorbid psychopathology (with anxiety being the most common diagnosis; Campo et al., Reference Campo, Bridge, Ehmann, Altman, Lucas, Birmaher and Brent2004). Considering that GI distress is such a common presentation for youth in primary care settings (McFerron & Waseem, Reference McFerron and Waseem2012), identifying how it relates to early adversity exposure, and whether it is associated with concurrent and future anxiety could help to identify a high-risk target group for early interventions. As the median age of anxiety onset is 11 years (Kessler et al., Reference Kessler, Berglund, Demler, Jin, Merikangas and Walters2005), and adversity effects on anxiety and related neurocircuitry have been detected in children under 10 years (e.g., Gee et al., 2013), it is likely that associations between adversity, anxiety, and GI complaints will be detected during childhood, but this has not been empirically tested. We sought to establish such associations in the current study.
Building a mechanistic understanding of the pathways through which adversity might affect both GI and anxiety symptoms is an important step toward improving outcomes following early adversity. In that regard, one important biological factor to consider is the gastrointestinal microbiome. The gastrointestinal microbiome, which is composed of a constellation of bacteria, archaea, fungi, and viruses, is essential for virtually all aspects of gastrointestinal function, from digestion and motility to inflammation and disease (Ursell, Metcalf, Parfrey, & Knight, Reference Ursell, Metcalf, Parfrey and Knight2012). The GI microbiome is increasingly recognized as important for brain function and mental health, and vice versa. For example, GI bacteria produce neurochemicals that are essential for brain function and emotional behaviors, such as ϒ-aminobutyric acid, short chain fatty acids, and, 5-hydroxytryptophan (the biological precursor to serotonin), which then reach the central nervous system (CNS) through humoral and vagal nerve pathways (bidirectional communication channels between the brain and gut; Cryan & Dinan, Reference Cryan and Dinan2012). In addition, the microbiome influences immune and inflammatory pathways (Belkaid & Hand, Reference Belkaid and Hand2014) that are directly associated with anxiety and depression (Vogelzangs, Beekman, de Jonge, & Penninx, Reference Vogelzangs, Beekman, de Jonge and Penninx2013; Vogelzangs, de Jonge, Smit, Bahn, & Penninx, Reference Vogelzangs, de Jonge, Smit, Bahn and Penninx2016), and manipulations of such bacteria have been shown to influence anxiety levels in adult humans (Collins, Kassam, & Bercik, Reference Collins, Kassam and Bercik2013; Lyte, Li, Opitz, Gaykema, & Goehler, Reference Lyte, Li, Opitz, Gaykema and Goehler2006; Messaoudi et al., Reference Messaoudi, Lalonde, Violle, Javelot, Desor, Nejdi and Cazaubiel2011), as well as fear behaviors in developing rodents (Callaghan, Cowan, & Richardson, Reference Callaghan, Cowan and Richardson2016; Cowan, Callaghan, & Richardson, Reference Cowan, Callaghan and Richardson2016).
When considered in the context of development, the role of the gastrointestinal microbiome in mental health is particularly pertinent, as there may be a sensitive period for microbiome-dependent maturation of the CNS. Rodent models using animals that were born and raised germ free have identified a causal role of GI bacteria in CNS maturation. In particular, disrupted development of several neural networks with known roles in emotional functioning, including the amygdala, prefrontal cortex, and hippocampus, have been observed in germ-free animals (Hoban et al., Reference Hoban, Stilling, Ryan, Shanahan, Dinan, Claesson and Cryan2016, Reference Hoban, Stilling, Moloney, Moloney, Shanahan, Dinan and Clarke2017, Reference Hoban, Stilling, Moloney, Shanahan, Dinan, Clarke and Cryan2018; Ogbonnaya et al., Reference Ogbonnaya, Clarke, Shanahan, Dinan, Cryan and O'Leary2015), aberrations that can be lifelong if microbial reconstitution does not occur before adolescence.
Beyond germ-free conditions, the early psychosocial environment also shapes microbiome maturation during sensitive developmental windows. For the microbiome, this period of maximally enhanced sensitivity to environments appears to extend from birth to around the fourth year of life, with both the overall numbers of species represented in the gut and their taxonomic relation to adult bacteria increasing dramatically across this time frame (Yatsunenko et al., Reference Yatsunenko, Rey, Manary, Trehan, Dominguez-Bello, Contreras and Anokhin2012). Thus, it seems that, like the brain, the young microbiome is particularly vulnerable to environmental perturbation during this early stage of developmental plasticity. Rodent and monkey models have shown that early caregiving adversity (e.g., maternal separation) significantly alters microbial communities across the life span (Bailey & Coe, Reference Bailey and Coe1999; Gareau, Jury, MacQueen, Sherman, & Perdue, Reference Gareau, Jury, MacQueen, Sherman and Perdue2007). Moreover, interventions that directly alter microbial communities (e.g., probiotics), if performed early enough, can reverse the effects of early adversity on affective behaviors (threat learning and extinction), anxiety phenotypes, and the brain (hippocampus and prefrontal cortex; Ait-Belgnaoui et al., Reference Ait-Belgnaoui, Colom, Braniste, Ramalho, Marrot, Cartier and Tompkins2014; Bravo et al., Reference Bravo, Forsythe, Chew, Escaravage, Savignac, Dinan and Cryan2011; Callaghan et al., Reference Callaghan, Cowan and Richardson2016; Cowan et al., Reference Cowan, Callaghan and Richardson2016; Waworuntu et al., Reference Waworuntu, Hain, Chang, Thiede, Hanania and Berg2014) in rodent models, suggesting the microbiome might be a treatment target for anxiety in adversity-exposed populations.
The most likely pathway linking microbiome and mental health is via an influence over CNS functioning. However, there is currently no evidence of brain–gut bacteria associations in middle childhood, or following adversity. In addition, it is not clear whether adversity is associated with GI symptoms in middle childhood, which would be important to establish if we are to identify children who would most likely benefit from gastrointestinal manipulations. Here we present the results of two studies, which together provide evidence for the utility of examining the brain–gut–microbiome axis in affective functioning following early adversity. In Study 1, we test whether adversity is associated with GI complaints in a large sample of youth aged 3–18 years. Then, we establish whether GI complaints, which are a common reason for presentation at primary care clinics, are associated with concurrent anxiety, and are predictive of future anxiety, validating clinical utility of examining GI complaints for detecting anxiety risk. In Study 2, which occurred in a smaller convenience sample taken from the larger Study 1 population, we establish proof-of-concept for microbiome–brain associations in middle childhood, following adversity, paving the way for future investigations to explore this mechanistic pathway across development using currently available techniques (e.g., functional neuroimaging).
Study 1: Associations of Adversity With Gastrointestinal Complaints and Anxiety
Method
Participants
Participants were children and adolescents who had either been exposed to early adverse (EA) caregiving experiences (i.e., institutional care, N = 103, or foster care, N = 12, total EA N = 115, followed by international adoption into families), or who had been raised with their biological families without any report of adverse caregiving (comparison group; COMP, N = 229). Institutional care is a significant deviation from species-expected caregiving experiences (e.g., high infant:caregiver ratios; Gunnar, Bruce, & Grotevant, Reference Gunnar, Bruce and Grotevant2000; van Ijzendoorn et al., Reference van IJzendoorn, Palacios, Sonuga-Barke, Gunnar, Vorria, McCall and Juffer2011) representing a potent stressor for the infant. Foster care, although often representing a vast improvement from institutional care settings (depending on the quality of caregiving), nonetheless involves a significant disruption in caregiving and is also associated with elevated risk for emotional difficulties in children (Brand & Brinich, Reference Brand and Brinich1999; Van Den Dries, Juffer, van IJzendoorn, Bakermans-Kranenburg, & Alink, Reference Van Den Dries, Juffer, van IJzendoorn, Bakermans-Kranenburg and Alink2012; Zeanah et al., Reference Zeanah, Egger, Smyke, Nelson, Fox, Marshall and Guthrie2009). Thus, we analyzed the data provided by those children as one group of EA youth, although the majority of the EA sample had experienced institutional care. Upon enrollment, the parents of comparison youth reported no child/adolescent psychiatric diagnoses, and as a group they scored in the average range on the Child Behavior Checklist (CBCL; mean = 44.87, SD = 10.44; Achenbach, Reference Achenbach1991). Inclusion criteria for all youth (N = 344) were that parents completed the CBCL and Revised Children's Anxiety and Depression Scale—Parent version (RCADS-P) for a baseline assessment of GI symptoms (Time 1).
The sample used in Study 1 was gathered from two projects that were being completed in the laboratory at the same time. As the measures used in this study were collected across both projects, we included participants from both projects in the current study. In one of those projects, N = 234 participants (ages 4–18 years old) were enrolled in a longitudinal imaging study, aiming for a longitudinal sample target of N =130 participants to receive two visits, and N = 65 to receive three visits (depending on entry date; i.e., not all participants were actually intended for follow-up). Participants were overenrolled at Time 1 in order to reach follow-up recruitment goals given the anticipated difficulties in scanning this age range (e.g., braces, motion artifact, refusal, and scheduling). Thus, this sample provided N = 234 participants with Time 1 data, out of which N = 117 participants also provided anxiety symptom data at either two (N = 114 [COMPS N = 69, EA N = 45]; i.e., 88% of N =130 targeted for follow-up assessment) or three (N = 60 [COMPS N = 33, EA N = 27]; i.e., 92% of N = 65 targeted for follow-up assessment) time points. For the current analysis, we counted the first time point at which participants had provided reports of gastrointestinal symptoms (via RCADS-P and CBCL items) as their Time 1 assessment, independent of the visit in which those data had been collected. If participants provided more than one time point of anxiety symptom data, those were labelled Time 2 (Time 1 + 2 years) for the second time point, and Time 3 (Time 1 + 4 years) for the third time point. In addition to those N = 234 youth, another N = 110 preschoolers (aged 3–6 years old) were recruited for a separate single time point study that examined emotional learning (these data were included in the Time 1 assessment for the current study). Hence, for the current study, the Time 1 assessment included a total of 344 participants between the ages of 3 and 18 years old (COMP N = 229, 120F/109M; EA N = 115, 75F/40M; see Table 1 for demographics). Across both of the past projects, recruitment of participants occurred through local birth records, local classifieds, international adoption agencies, family networks, posted flyers, and friend referral.
Table 1. Demographic information stratified by caregiving history for Study 1 (rows 2 and 4) and Study 2 (rows 3 and 5)
Note: Country of origin is listed for the early adversity (EA) sample, whereas race is listed for the US-born comparison sample. Sex distribution; mean IQ (SD); mean age in years (and range) at the time of the scan, at the time of the gastrointestinal/anxiety assessment, at MRI scan, and microbiome sample, and the interval between the scan and sample; as well as mean (SD) anxiety score on the Screen for Child Anxiety Related Disorders—Parent version (SCARED-P) at baseline, are listed separately for each group. Statistical significance of group differences between caregiving conditions for Study 1 and 2 are presented in the bottom two rows. Mean age in years (and range) at the time of adoption and years with the adoptive family by the time of stool sample are listed for the EA group only.
The Wechsler Abbreviated Scale of Intelligence provided estimated IQ for participants over 6 years of age (data available on N = 232). Mean levels of measured intelligence were average for both groups (COMP: mean ± SD = 111.48 ± 16.45; IA: mean ± SD = 103.35 ± 16.93) but were significantly higher in the comparison youth than in children who were internationally adopted (controlling for participant sex), F (1, 229) = 17.54, p < .0001. Scores on the Wechsler Abbreviated Scale of Intelligence were not associated with gastrointestinal symptoms, controlling for caregiving history and sex; GI Factor 1: β = .00, t (3) = –0.12, p = .902; GI Factor 2: β = .00, t (3) = 0.95, p = .342, or with anxiety measured at the Time 1 assessment, β = .00, t (3) = –1.95, p = .053. Modal household income per year for both caregiving groups was between $100,000 and $150,000, but was statistically higher in adopting families than in the comparison sample, F (1, 331) = 19.69, p < .0001 (though both were well above the national mean $38,321; US Census Data, for 1 or more child families in 2015). Income was not associated with gastrointestinal symptoms, Factor 1: β = –.02, t (4) = –0.91, p = .361, Factor 2: β = .00, t (4) = 0.14, p = .890, or anxiety at the baseline assessment, β = .01, t (4) = –1.87, p = .63 (controlling for age, sex, and caregiving group). Median parental education was a 4-year college degree in comparisons and a master's degree in the EA group, which also differed across groups, F (1, 337) = 15.63, p < .0001. Differences in household education were also not associated with gastrointestinal or anxiety symptoms: gastrointestinal distress Factor 1: β = –.04, t (4) = –1.14, p = .256; digestive issues Factor 2: β = .03, t (4) = –1.04, p = .298; anxiety at baseline assessment: β = .01, t (4) = 1.84, p = .070 (controlling for participant age, sex, and caregiving group).
The institutional review board at the University of California, Los Angeles, approved the protocol. Parents provided written consent, children 7+ years old provided written assent, and children under 7 provided verbal assent.
Child gastrointestinal symptoms (parent report)
To assess gastrointestinal symptomatology, we examined gastrointestinal distress items included within the CBCL for age 1.5–5 and 4–18 years (Achenbach, Reference Achenbach1991) and the RCADS-P (Chorpita, Moffitt, & Gray, Reference Chorpita, Moffitt and Gray2005). The CBCL uses a 3-point scale (0 = not true, 1 = somewhat true, 2 = very true). Only items asssessing gastrointestinal symptoms that were present in both CBCL 1.5–5 years and 4–18 years were included in the analysis, such that the same GI symptoms could be assessed across the entire age range. The CBCL items assessed incidence of nausea (Question 45 [1.5–5 years], or Question 56c [4+ years]: “child suffers nausea with no known medical cause”), stomach aches/cramps (Question 78 [1.5–5 years] or Question 56f [4+ years]: “child suffers stomach aches with no known medical cause”), expulsion (Question 93 [1.5–5 years] or Question 56g [4+ years]: “child suffers vomiting with no known medical cause”), and constipation (Question 12 [1.5–5 years] or Question 49 [4+ years]: “child suffers constipation with no known medical cause”). The RCADS-P uses a 4-point scale (0 = never, 1 = sometimes, 2 = often, 3 = always). The RCADS-P item assessed gastrointestinal somatic symptoms (Question 3: “When my child has a problem he/she gets a funny feeling in their stomach”). The five items from the CBCL and one item from the RCADS-P were subjected to a factor analysis to summarize their underlying covariation.
Child anxiety assessment (parent report)
Anxiety-related behavior was assessed through parent report on the Screen for Child Anxiety Related Disorders—Parent version (SCARED-P; Birmaher et al., Reference Birmaher, Khetarpal, Brent, Cully, Balach, Kaufman and McKenzie Neer1997), which provides a continuous measure of anxiety. Because we were interested in anxiety behavior counts (not diagnoses) across this wide age range, we used the SCARED-P, which was originally validated for clinical cutoffs on a sample of children aged between 8 and 18 years. To rule out the possibility that including SCARED-P data from children younger than 8 years old altered the nature of the association between gastrointestinal symptoms and SCARED-P, our analyses included a confirmatory analysis in children aged 8 and over; the results of the two analyses did not differ. Parents reported on the frequency of each behavior for their child on a 3-point scale (0 = never, 1 = sometimes, 2 = often), and these responses were summed to create a behavior count and then averaged to create a mean score. One question on the SCARED-P (Item 11: “my child gets stomach aches at school”) was removed from the total summed subscale to avoid overlap with the gastrointestinal measure. Parents completed the measure at up to three time points, spaced approximately 2 years apart, providing a longitudinal assessment of child anxiety. Thus, anxiety could be assessed across a 5-year time frame.
Longitudinal attrition analysis
Of the N = 130 participants with planned follow-up contrasts (i.e., a subset of the large parent grant study [N = 234] who were targeted for follow-up), 117 were actually followed-up (either one or two additional time points). Thus, of the 117 subjects for which no follow-up data was obtained, only N = 13 subjects represent true attrition (i.e., they were intended for follow-up but did not return). To ensure that true attrition (i.e., within the youth intended for follow-up) was not being driven by systematic differences in anxiety at baseline between youth who returned for intended follow-up and those who did not return for intended follow-up, we performed an analysis of variance (ANOVA) with baseline anxiety as the dependent variable and attrition status (lost to follow-up vs. successfully followed) as the predictor (controlling for caregiving group, sex, and age at Time 1). Youth who were lost to follow-up did not differ in Time 1 anxiety scores from youth who were successfully followed, F (1, 181) = 1.22, p = .270. To ensure that there were no systematic differences in Time 1 anxiety between youth who were followed-up (regardless of whether they were intended to be followed) we performed another ANOVA with Time 1 anxiety as the dependent variable, and follow-up status (follow-up vs. no follow-up data) as the predictor (controlling for caregiving group, sex, and age at Time 1). Youth who provided follow-up data did not differ from youth who did not provide follow-up data on Time 1 anxiety levels, F (1, 271) = 0.52, p = .470. (see Figure 1 for a graphical depiction of the time line of Study 1 and Study 2).
Figure 1. Study timeline for Study 1 (full sample) and Study 2 (microbiome–fMRI proof-of-concept subsample). Associations between gastrointestinal symptoms at Time 1 and anxiety symptoms at Times 1, 2, and 3 (Years 1, 3, and 5) of the study were made. In addition, associations between fMRI at Time 1 and microbiome composition at Time 3 are explored. Times 2 and 3 are presented on increasingly translucent backgrounds to highlight the decreasing number of participants who were intended for follow-up at each time point.
Statistical tests
All analyses were conducted in SPSS. All statistical tests were two-sided and the alpha value was set as α = 0.05. We used separate linear regressions to assess associations between gastrointestinal distress/digestive issues and anxiety measured at Visit 1. We then tested a cross-sectional mediation model using the PROCESS macro within SPSS (Hayes, Reference Hayes2012; Model 4) with caregiving group as the independent variable, anxiety at Visit 1 assessment on the SCARED-P as the outcome variable, and gastrointestinal symptoms as the mediator (controlling for participant age and sex at the level of the mediator and dependent variable). Bias corrected 95% confidence intervals (CI) using 5,000 bootstrapped samples are reported. Cross-sectional, rather than longitudinal, mediation was used as we had the highest number of participants at the baseline (i.e., Time 1) assessment, to address our claim of determining whether significant variance in the association between adversity and anxiety was explained by gastrointestinal symptoms. That is, we were not interested in the predictive nature of the mediation model.
To assess prospective associations between gastrointestinal distress and anxiety across time, we performed a linear mixed model in SPSS with maximum likelihood estimation to accommodate the nested structure of the data (individual change in anxiety from Time 1 to Time 3). This method captures individual variance while allowing for missing data points, thus dealing with attrition (note: as explained above, due to study design, the fewer data points available at Time 2 and Time 3 than at Time 1 was expected, and there were no systematic differences in anxiety between followed and nonfollowed participants). Caregiving condition (comparison vs. early adversity), time (baseline, Time 2, and Time 3), and gastrointestinal distress scores were entered as fixed effects predictors of SCARED-P anxiety scores, with random slope and intercept between individuals. Age and sex (male vs. female) were entered as covariates in the model.
Results
Factor analysis of gastrointestinal symptoms
Principal components factor analysis was used to identify and compute composite scores for factors underlying the items that were assessing gastrointestinal distress extracted from the CBCL and RCADS-P (nausea, stomach ache/cramp, funny feeling in stomach, vomiting, and constipation). Four of the five items examined correlated at least .3 with at least one other item, indicating that factor analysis was a reasonable approach to condense these data. Moreover, the Kaiser–Meyer–Olkin measure of sampling adequacy was .66, which is above the recommended cutoff value of .6, and Bartlett's test of sphericity was also significant, χ2 (10) = 130.84, p < .0001. The diagonals in the anti-image correlation matrix were all over .5, and the extracted communalities were all over .3 also suggesting shared variance among the items that could be reasonably condensed through factor analysis.
Initial eigenvalues indicated that two factors explained 36.07% and 20.33% of the variance, respectively, whereas Factors 3–5 had eigenvalues <1 and explained 17.10%, 14.58%, and 11.92% of the variance, respectively. The two-factor solution was examined using varimax rotation of the factor-loading matrix and was found to be a good fit to the data with all primary loadings exceeding .5, and no cross-loading exceeding .4 (factor loading matrix is presented in Table 2). The first factor was labeled “gastrointestinal distress” (GI distress) due to its being composed of items measuring nausea, somatic complaints, and stomach aches and cramps. The second factor was labeled “digestive issues” due to its being composed of items that reflected digestive issues: vomiting and constipation.
Table 2. Factor loadings on the varimax rotated component matrix for gastrointestinal symptoms assessed on the Child Behavior Checklist (CBCL) and the Revised Children's Anxiety and Depression Scale (RCADS) in N = 344 youth
Note: The two-factor solution of “gastrointestinal distress” and “digestive issues” was the best fit of the data.
Association of caregiving history with gastrointestinal distress and digestive issues
Composite scores on the two factors were entered as outcome variables in separate linear regressions with caregiving group as the independent variable of interest (COMP vs. EA), covarying for the effects of sex and age. EA exposure was significantly associated with higher GI distress scores (i.e., Factor 1), β = .25, t (340) = 2.20, p = .028, 95% CI [0.03, 0.47] (Figure 2a), as was increasing child age, β = .04, t (340) = 3.22, p = .001, 95% CI [0.02, 0.07] (online-only Supplemental Figure S.1), but sex and GI distress were not associated, β = .17, t (340) = 1.59, p = .113, 95% CI [–0.04, 0.38]. EA exposure, was also associated with higher scores on the digestive issues factor (i.e., Factor 2), β = .41, t (340) = 3.57, p = .0004, 95% CI [0.18, 0.63] (Figure 2b), but age and sex were not, age: β = .00, t (340) = –0.07, p = .943, 95% CI [–0.02, 0.03], sex: β = .11 t (340) = 0.98, p = .330, 95% CI [–0.11, 0.32]. Within the EA group, using child sex and age as covariates, age of adoption was not associated with GI Distress, β = .00 t (109) = 0.95, p = .342, 95% CI [–0.01, 0.01], or with digestive issues, β = .00, t (109) = –0.55, p = .582, 95% CI [–0.01, 0.01]. See Supplemental Figure S.1 for the association between caregiving history and GI distress broken down by child age range: preschool, middle childhood, and adolescence.
Figure 2. Association of caregiving history group membership (early adversity exposed [EA] or comparison [COMP] group) with (a) Factor 1, gastrointestinal distress, and (b) Factor 2, digestive issues. Bars represent the group averages. An * denotes a statistically significant result. P values and 95% confidence intervals (CI) of the beta estimate from the regression analysis are reported above each graph.
Association of GI distress and digestive issues with youth anxiety
Separate linear regressions were used to test the association of scores on each factor (GI distress or digestive issues) with anxiety (measured with the SCARED-P), controlling for caregiving group, age, and sex.
Factor 1: GI distress
GI distress was significantly associated with Time 1 anxiety scores on the SCARED-P, β = .08, t (325) = 6.04, p < .0001, 95% CI [0.06, 0.11] (Figure 3a), as was sex, β = .07, t (325) = 2.78, p = .006, 95% CI [0.02, 0.13], and caregiving group, β = .18, t (325) = 6.25, p < .0001, 95% CI [0.12, 0.23], but not age, β = .00, t (325) = 0.67, p =.506, 95% CI [–0.01, 0.01]. The association of GI distress and anxiety was confirmed when including only children over 8 years of age (N = 168); β = .09, t (164) = 4.91, p < .0001, 95% CI [0.05, 0.12].
Figure 3. Scatter plots illustrating the association between mean parent-rated child anxiety scores on the Screen for Child Anxiety Related Disorders—Parent version (SCARED-P) at Visit 1 with (a) Factor 1, gastrointestinal (GI) distress, and (b) Factor 2, digestive issues. Dots are colored by caregiving group membership, where children who had experienced early adversity (EA) are represented in open circles, and children from the comparison group (COMP) are presented in gray transparent circles (transparency is utilized so that overlapping values can be viewed). (c) The cross-sectional mediation model where children who were exposed to early adversity (coded positively) exhibit elevated anxiety through a pathway of increased GI distress.
Factor 2: Digestive issues
Digestive issues were not significantly associated with variance in SCARED-P anxiety scores, β = .03, t (325) = 1.83, p = .068, 95% CI [0.00, 0.05] (Figure 3b), though digestive Issues were generally low incidence overall. Similarly, age was not associated with SCARED-P anxiety scores, β = .01, t (325) = 1.76, p = .080, 95% CI [0.00, 0.01], but sex, β = .09, t (325) = 3.10, p = .002, 95% CI [0.03, 0.14], and caregiving group, β = .19, t (325) = 6.35, p < .0001, 95% CI [0.13, 0.25] were. As with Factor 1, the nonsignificant association between digestive issues and anxiety was confirmed when including only children over 8 years of age (N = 168); β = .02, t (164) = 1.16, p = .250, 95% CI [–0.02, 0.06].
Elevated GI distress explains variance in the association between early adversity and concurrent anxiety
To determine whether GI distress was a mediator of the relationship between caregiving history and anxiety, we tested a full cross-sectional mediation model with caregiving group as the independent variable, Time 1 anxiety scores on the SCARED-P as the dependent measure, scores on the GI distress factor at Time 1 as the mediator, controlling for participant age and sex. Cross-sectional mediation was performed as Time 1 had the highest number of data points, and to explain variance in, rather than “predict,” the anxiety outcome. The mediation was significant; the path between caregiving group and scores on the GI Distress measure (path a) was significant, β = .28, SE = .12, t (326) = 2.42, p = .016, 95% CI [0.05, 0.51], as were the paths between GI distress and anxiety (path b), β = .08, SE = .01, t (325) = 6.04, p < .0001, 95% CI [0.05, 0.11], and between EA caregiving history and anxiety (path c’), β = .18, SE = .03, t (325) = 6.25, p < .0001, 95% CI [0.12, 0.23]. The bootstrapped confidence interval for the mediation effect (pathway a*b; i.e., the indirect association of EA group with increased anxiety via elevated GI distress), did not contain zero; β = .02, SE = .01, bootstrapped CI [0.003, 0.05] (Figure 3c).
Association of anxiety change across time with GI distress
To assess whether GI distress was associated with change in anxiety scores across time, we performed a linear mixed-model regression. There were significant main effects of GI distress, F (29, 275.27) = 4.14, p < .0001, ηp2 = .32, caregiving group, F (1, 353.40) = 24.50, p < .0001, ηp2 = .06, sex, F (1, 308.34) = 7.26, p = .007, ηp2 = .02, and time, F (2, 111.63) = 14.66, p < .0001, ηp2 = .21, such that individuals with higher GI distress at Time 1, those from the EA group, and females, had higher anxiety symptoms, and anxiety symptoms decreased across time. There was also a significant GI Distress × Time interaction, F (31, 113.00) = 5.80, p < .0001, ηp2 = .61, and a Caregiving Group × GI Distress interaction, F (8, 313.60) = 3.39, p = .001, ηp2 = .19 (Figure 4a). There was no main effect of age at Time 1, F (1, 330.28) = 0.05, p = .829, ηp2 < .0001, nor were the interactions between Caregiving Group × Time, F (2,104.89) = 2.99, p = .054, ηp2 = .05, and GI Distress × Caregiving Group × Time, F(6, 129.15) = 1.47, p = .193, ηp2 = .01, significant (but see Figure 4b for associations broken down by caregiving condition). Post hoc linear regressions of the association between anxiety and GI distress (controlling for age, sex, and caregiving group) revealed a positive association between GI distress and anxiety at baseline (i.e., Time 1), β = .08, t (325) = 6.04, p < .0001, CI [0.06, 0.11], and Time 2, β = .05, t (109) = 2.57, p = .012, CI [0.01, 0.09], but by Time 3 there was no longer an association between anxiety and GI distress, β = .00, t (55) = 0.114, p = .909, CI [–0.04, 0.05], such that GI distress in early life was associated with heightened anxiety symptoms concurrently and when measured 2 years later. See the online-only supplemental data for a confirmatory analysis using a repeated-measures ANOVA with listwise deletion, which shows that imputed data were not driving the results reported here.
Figure 4. (a) Association between Factor 1, gastrointestinal (GI) distress scores (median split into high GI distress and low GI distress scores), and predicted values of parent-rated child anxiety on the Screen for Child Anxiety Related Disorders—Parent version (SCARED-P) at Visits 1, 2, and 3. (b) The same association broken down by caregiving group with children from the comparison (COMP) group represented in the solid lines, and children in the early adversity (EA) exposed group in dotted lines, with low GI distress in black and high GI distress in red.
To ensure that the association of gastrointestinal symptoms at Time 1 with anxiety symptom change across Time (Time 1–Time 2) was not confounded by continued gastrointestinal symptoms, we tested another model that covaried for gastrointestinal symptoms at later time points (Time 2 and 3). Because the same factor structure may not emerge at all time points, we instead created a standardized sum score of the items from the Factor 1 structure that emerged at Time 1 (baseline) and used those as covariates at Time 2 and Time 3. Due to the study design, which resulted in a smaller sample size at Time 2 and Time 3 than at Time 1, entering covariates for Time 2 and Time 3 substantially reduced the sample size (N = 62). Nonetheless, the same basic associations persisted. There was a significant effect of time, F (2, 67.73) = 14.67, p < .0001, ηp2 = .30, baseline GI distress, F (13, 69.46) = 3.24, p = .001, ηp2 = .38, and GI distress at Time 2, F (1, 57.92) = 13.26, p = .001, ηp2 = .19, such that anxiety symptoms decreased across time in everyone, but remained higher in individuals with elevated GI distress at baseline and Time 2. Most important, the significant Time × Baseline GI Distress interaction reported above was also seen in these models, F (26, 67.82) = 5.29, p < .0001, ηp2 = .67. There were no main effects of caregiving group, F (1, 69.68) = 0.55, p = .462, ηp2 = .007, sex, F (1, 58.21) = 1.47, p = .230, ηp2 = .02, and age, F (1, 58.39) = 1.00, p = .321, ηp2 = .02, nor were the interactions Caregiving Group × Baseline GI Distress, F (2, 69.13) =1.00, p = .370, ηp2 = .03, Caregiving Group × Time, F (2, 67.97) = 2.50, p = .090, ηp2 = .07, or the three-way interaction Caregiving Group × Time × Baseline GI Distress, F (4, 68.35) = 2.12, p = .088, ηp2 = .11, significant.
Study 2: Microbiome-Brain “Proof-of-Concept” Associations
Method
Participants
We tested a “proof-of-concept” association between early caregiving adversity, the gastrointestinal (GI) microbiome, and brain reactivity to threat stimuli (fear faces) in a smaller “convenience sample” that was a subset of the Study 1 population (N = 16; mean age = 8 years, range = 5–11 years). The children used in Study 2 were selected if they had provided usable task-based magnetic resonance imaging (MRI) data in childhood (ages 5–11 years), and agreed to donate a stool sample. This age range was chosen as we have previously documented it as a sensitive period for prefrontal cortex development in youth (Gee, Humphreys, et al., Reference Gee, Humphreys, Flannery, Goff, Telzer, Shapiro and Tottenham2013). The functional MRI (fMRI) data from this subsample has been published previously to address a different question of interest (i.e., not related to the microbiome; Gee, Gabard-Durnam, et al., Reference Gee, Gabard-Durnam, Flannery, Goff, Humphreys, Telzer and Tottenham2013; Gee, Humphreys, et al., Reference Gee, Humphreys, Flannery, Goff, Telzer, Shapiro and Tottenham2013). As this was a convenience sample, the microbiome data were collected after children had come in for their MRI scan (average 3.95 years, SD = 0.55 years, after scanning/anxiety assessments). The GI microbiome was analyzed through taxonomic identification of bacteria in stool samples using 16S rRNA as a marker gene for bacterial identification. No exclusions were made for any of the analyses. Adversity-exposed children in Study 2 were all exposed to prior institutionalization (i.e., no child experienced foster care). Sample demographics (split into Study 1 and 2) are presented in Table 1.
Child diet assessment (parent report)
Parents completed a daily diary of their children's food intake in the 2 days preceding stool sample collection. Using online nutrition programs, the food diary was coded into amount of protein, carbohydrate, and fat (which can influence bacterial composition) for each meal and was averaged across 2 days. These variables were then included as covariates in offline analyses of potential confounding factors. In addition, parents answered questions that assessed the presence and frequency of children's gastrointestinal symptoms/pain (e.g., diarrhea, bloating, and pain in the upper or lower abdomen). Frequency of gastrointestinal symptoms was coded on a 4-point scale corresponding to the judgments never, infrequent, frequent, and always. Frequency of pain was coded on a 5-point scale corresponding to the judgments never, 1–3 times per month, once per week, several times a week, and every day. No children were on antibiotic or antifungal medication, and only 1 child had a strictly vegetarian diet. As these data are only available in N = 16 subjects from Study 2, results of this gastrointestinal symptom questionnaire are presented in Table 3, but were not included in the factor analysis in Study 1.
Table 3. Mean (SD) responses for individual questions on the parent proxy-reported gastrointestinal symptoms questionnaire stratified by caregiving history in the subsample (Study 2) only
Note: For both groups n = 7. T tests were used to assess group differences between children from the comparison (COMP) and early adversity (EA) exposed groups.
Stool sample collection
Parents were asked to help their child collect a small stool sample. Parents were given an OMNIgene.GUT (DNA Genotek) home collection kit, which included a spatula, collection tube with stabilizing liquid (ethyl alcohol <24%), and instructions on how to collect, homogenize, and store the sample. This method of sample stabilization has been shown to successfully preserve bacteria concentrations in stool for up to 8 weeks at room temperature (Song et al., Reference Song, Amir, Metcalf, Amato, Xu, Humphrey and Knight2016). Each child provided one sample, collected at home, and then sent directly to the sequencing lab: Research and Testing Laboratory (RTL Genomics, Texas).
16s rRNA sequencing and preprocessing
Bacterial DNA was extracted from stool samples provided by 16 children using RNA Power Soil kits (Mo Bio Laboratories). Illumina MiSeq, paired-end sequencing across 250 base pairs was used. Samples were amplified for sequencing in a two-step process using primers spanning the V1-V2 hypervariable region of the 16S rRNA gene. The forward primer was constructed with (5’−3’), the Illumina i5 sequencing primer (TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG), and the 28F primer (GAGTTTGATCNTGGCTCAG). The reverse primer was constructed with (5’−3’), the Illumina i7 sequencing primer (GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG), and the 388R primer (TGCTGCCTCCCGTAGGAGT). Amplifications were performed in 25 μl reactions with Qiagen HotStar Taq master mix (Qiagen Inc, Valencia, California), 1 μl of 5uM primer, and 1μl of template. Reactions were performed on ABI Veriti thermocyclers (Applied Biosytems, Carlsbad, California) under the following thermal profile: 95 ○C for 5 min, then 25 cycles of 94 ○C for 30 s, 54 ○C for 40 s, 72 ○C for 1 min, followed by one cycle of 72 ○C for 10 min and 4 ○C hold. Products from the first stage amplification were added to a second polymerase chain reaction (PCR) based on qualitatively determined concentrations. Primers for the second PCR were designed based on the Illumina Nextera PCR primers as follows: Forward, AATGATACGGCGACCACCGAGATCTACAC[i5index]TCGTCGGCAGCGTC, and Reverse, CAAGCAGAAGACGGCATACGAGAT[i7index]GTCTCGTGGGCTCGG. The second stage amplification was run the same as the first stage except for 10 cycles. Amplification products were visualized with eGels (Life Technologies, Grand Island, New York). Products were then pooled equimolar, and each pool was size selected in two rounds using Agencourt AMPure XP (BeckmanCoulter, Indianapolis, Indiana) in a 0.7 ratio for both rounds. Size selected pools were then quantified using the Quibit 2.0 fluorometer (Life Technologies) and loaded on an Illumina MiSeq (Illumina, Inc. San Diego, California) 2 × 300 flow cell at 10 pM.
Illumina reads were quality controlled by removal of short sequences, singleton sequences, and noisy reads. Denoising was performed on each region by trimming back to the last base where the total average of the sequence was greater than 25. Prefix dereplication was then performed using the USEARCH algorithm (Edgar, Reference Edgar2010), clustering at the 4% divergence level. Selection of operational taxonomic units (OTUs) was performed using the UPARSE OTU selection algorithm (Edgar, Reference Edgar2013). OTU selection was performed de novo, and OTU sequences were then classified against an in-house database of high-quality sequences derived from the National Center for Biotechnology Information and maintained at Research & Testing Laboratory Genomics. Chimera detection was then performed on the selected OTUs using the UCHIME chimera detection software executed in de novo mode (Edgar, Haas, Clemente, Quince, & Knight, Reference Edgar, Haas, Clemente, Quince and Knight2011), and chimeric sequences removed. The remaining sequences were then corrected base by base to aid in the removal of noise. Denoised and chimera checked reads were then condensed into FASTA formatted files for taxonomic identification and community analyses. To assign taxonomic information the sequences were run through the USEARCH global alignment program. All preprocessing steps were performed blind to the caregiving history and characteristics of the children.
Microbial richness and diversity analyses
Figure 5 provides a description of different metrics for measuring community composition of the microbiome that are used in this paper. We visualized genera-level taxonomic community composition through the VAMPS (Visualization and Analysis of Microbial Population Structure) community visualization tool (Huse et al., Reference Huse, Mark Welch, Voorhis, Shipunova, Morrison, Eren and Sogin2014). Ribosomal Database Project (RDP) taxonomy was assigned, and within-subject differences in the number of bacterial groups (in this case, genus) present in the sample (i.e., alpha diversity) were calculated using two estimators: observed richness (which provides a count of the number of different genera represented in the sequenced sample), and Shannon's index (which provides a metric of how rich and evenly distributed the groups are). Statistical differences in alpha diversity between the groups were analyzed using separate linear regressions (controlling for age and sex). Variances between the caregiving groups were similar. Between-subject differences in bacteria present in the sample (i.e., beta diversity) was quantified using an approach that takes into account the phylogenetic relations between genera, Unifrac distance, calculated within the phyloseq package in R (McMurdie & Holmes, Reference McMurdie and Holmes2013). Paired samples that contain bacteria that are closely related on the phylogenetic tree (i.e., share many tree branches) score low on beta diversity, whereas paired samples that are not related (or are distantly related) on the phylogenetic tree will receive high beta diversity scores, resulting in a distance matrix of between-samples beta diversity scores for every sample pair in the study. A nonparametric permutational multivariate analysis of variance was used to compute an R 2 and significance value for the amount of variance in the distance matrices that was explained by the grouping variable of interest (caregiving group: COMP or EA). Two of such analyses were conducted: weighted (which weights the dominant bacteria) and unweighted (which allows the rare bacteria to be influential in the beta diversity score).
Figure 5. Explanatory figure of microbial community assessment metrics for alpha and beta diversity. The within-subject metric of alpha diversity called “observed richness” quantifies the total count of bacterial groups at a given taxonomic level (e.g., species, genera, or phyla) in a sample from an individual person in the study. Higher counts of bacterial groups present in the sample indicate higher observed richness. Note that only the bacterial groups, not the numbers of members within each bacterial group, are quantified. The second within-subject metric of alpha diversity measured in this study is Shannon's diversity, which provides a measure of how rich and evenly distributed the bacterial groups are. Samples scoring highest on Shannon's diversity will have a high number of bacterial groups represented (richness), and the members of each bacterial group will be approximately equal (i.e., evenly distributed). Samples with very low counts of bacterial groups (richness) and where a small subset of bacterial groups dominates the community (uneven distribution of members) will result in the lowest Shannon's diversity score. The between-subject metric of beta diversity measured in this study was Unifrac distance. This metric calculates the relative closeness of bacterial groups found in paired samples from people on a phylogenetic tree. Samples that are closely related on the tree will share many branches and will receive a low Unifrac distance score (unshared/all branches). Samples from two individuals where the bacterial groups are far away from each other on the phylogenetic tree will share few branches and will receive a high Unifrac distance score. A distance matrix between each pair of samples in the study is constructed.
Bacterial biomarker analysis
Differences in specific bacteria between caregiving groups were identified with a linear discriminant analysis of effect sizes (LEfSE) performed through the Galaxy platform (Segata et al., Reference Segata, Izard, Waldron, Gevers, Miropolsky, Garrett and Huttenhower2011). In short, LEfSE uses the nonparametric sum rank Kruskal–Wallis test to detect bacteria that differ between the classes of interest (in this case, care history), followed by linear discriminant analysis (Fisher, Reference Fisher1936) to detect the effect size of each of these differences. Due to the relatively small sample size, we carefully checked the resulting potential biomarkers for outlying values (>3 SD from the mean) and excluded any “biomarker” that had at least one value exceeding the threshold.
MRI task paradigm
During the fMRI scan, participants completed two runs of an emotional faces task. The task consisted of a mixed design with one blocked variable (emotional valence: fear and happy) and one event-related variable (emotion and neutral). During each run, participants viewed faces presented one at a time that were either emotional or neutral. The order of runs (i.e., happy/neutral, fear/neutral) was counterbalanced across participants, and the stimulus order within each run was randomized and fixed. To ensure that participants were paying attention, they were instructed to press a button with their index finger when they saw a neutral face. The faces were selected from the Karolinska Directed Emotional Faces database (Lundqvist, Flykt, & Ohman, Reference Lundqvist, Flykt and Ohman1998). The faces were presented in color at a visual angle of approximately 15 degrees. The probability of an emotional face was 50% on any given trial. Stimuli were jittered (variable intertrial interval ranging from 3000 to 9000 ms) and randomized based on a genetic algorithm in order to allow for unique estimates of the hemodynamic response for each trial type. Each run contained 48 trials (24 neutral faces, 24 fearful or happy faces). Each face remained on the screen for 500 ms. Because fear faces have been shown to activate prefrontal cortex circuits in developing youth (Gee, Gabard-Durnam, et al., Reference Gee, Gabard-Durnam, Flannery, Goff, Humphreys, Telzer and Tottenham2013; Gee, Humphreys, et al., Reference Gee, Humphreys, Flannery, Goff, Telzer, Shapiro and Tottenham2013), and are known to engage threat circuitry relevant to anxiety (Baird et al., Reference Baird, Gruber, Fein, Maas, Steingard, Renshaw and Yurgelun-Todd1999), analyses focused on fear faces (i.e., we examined blood oxygen level dependent responses to fear faces, relative to implicit baseline). As participants were instructed to press a button when they saw a neutral face, effectively rendering the neutral face salient, we decided not to examine blood oxygen level dependent responses to neutral faces. Rather than trying to isolate any one neural process in our examination of the fear versus implicit baseline contrast, we are using fear faces as a probe of affective circuitry function.
fMRI data acquisition
Scanning was performed in a Siemens Trio 3.0 tesla MRI scanner. A standard radiofrequency 12-channel head coil was used. A whole-brain high resolution T1 weighted anatomical scan (MPRAGE; 256 × 256 in-plane resolution; 256 mm field of view; 192 × 1 mm sagittal slices) was acquired for each participant in order to register and localize their functional data to standard space. The emotional faces task was presented on a computer screen within MR-compatible goggles. T2* weighted echoplanar images (interleaved) were collected at an oblique angle of ~30 degrees (130 volumes/run; resonance time 2000 ms; echo time, 30 ms; flip angle, 90 degrees; matrix size, 64 × 64; field of view, 192 mm; 34 slices; 4-mm slice thickness; skip 0 mm; 24 observations per event type).
fMRI data analysis
Preprocessing
Each participant's fMRI scan data was motion corrected using the Analysis of Functional Neuroimages (AFNI; Arumugam et al., Reference Arumugam, Raes, Pelletier, Le Paslier, Yamada, Mende and Bork2011) 3dvolreg; the first volume was the reference image. Images were then slice time corrected in AFNI, and aligned to each subjects MPRAGE image using FLIRT in FMRIB Software Library (FSL; Jenkinson, Beckmann, Behrens, Woolrich, & Smith, Reference Jenkinson, Beckmann, Behrens, Woolrich and Smith2012), before being scaled to the mean of each voxel. MPRAGE structural images were skull stripped in FSL using BET, and FAST was used to segment gray matter, white matter, and ventricles (cerebrospinal fluid; CSF). These anatomical parcellations were then used to extract signals from the white matter and CSF in the fMRI scan. MPRAGE scans were warped to standard template space from the Montreal Neurological Institute (Montreal, Canada) using a nonlinear registration procedure (FLIRT and FNIRT) with boundary-based registration cost algorithm in FSL, and resulting warp transformation saved for application to the functional data. Preprocessed functional data were spatially aligned and normalized to each participant's warped MPRAGE using the nonlinear registration (FNIRT). Functional data were spatially smoothed with a 6-mm full width at half max kernel. Montreal Neurological Institute transformed images had a resampled resolution of 3 mm3. Time series were normalized to percent signal change to allow for comparison across individuals and runs. The functional runs were concatenated before creating two individual-level models for each participant to model functional reactivity and connectivity.
Motion correction
Systematic procedures were implemented to reduce motion due to the age of the children. Children participated in a mock scanning session in the lab before the actual MRI. This mock session allowed children to become acclimated to the scanning environment, noises, and procedures, to get real-time feedback on their ability to lie still, and to practice that skill. During data collection, an air vacuum pillow (Siemens Comfort Pack) was used to pad and secure the child's head in a comfortable, steady position. Additional padding was placed around the child's head as required. Children also received reminders to reduce motion and feedback on their level of motion throughout the scanning session.
After the scanning session, multiple steps were taken to correct for motion. Following recommended standards (Power, Barnes, Snyder, Schlaggar, & Petersen, Reference Power, Barnes, Snyder, Schlaggar and Petersen2012) all analyzed data were free of motion >0.9 mm in any direction. Volumes where the motion derivatives had a Euclidean norm >0.9 mm were excluded through censoring that volume and the subsequent volume. All participants had <40% of the total volumes censored (mean percentage of censored volumes, 8.12%, mode, 0%). Demeaned motion in each of six directions, and their temporal derivatives, were also included as regressors in the trial level regression (12 motion regressors total). At the group-level analysis, each subject's average motion (mean displacement derivative) was also entered as a regressor. We also verified that our clusters of bacteria-related reactivity did not overlap with clusters that were related to average motion in the whole brain maps.
Whole-brain reactivity
To examine whether the microbiome was associated with reactivity to fear faces, a whole-brain reactivity analysis was performed. A Generalized linear model analysis was performed in AFNI for each participant with regressors for task, CSF signal, white matter signal, 6 motion regressors, plus their 6 backward-temporal derivatives. The four task regressors modeled which face comprised each trial (i.e., fearful, happy, neutral faces [in the fearful run], and neutral faces [in the happy run]). Analyses focused on the contrast of fearful faces against implicit baseline (fixation). The generalized linear model analyses fit the percentage of signal change in the brain to each regressor, and linear and quadratic trends were modeled for the time series of each voxel to control for correlated drift.
Following individual-level analyses, the regression coefficients were then subjected to random-effects, group level analysis using the 3dttest command in AFNI. Two group-level analyses were run with participant level of bacterial genera (Bacteroides or Lachnospiraceae) entered as the regressor of interest in each, whereas participant caregiving condition, age at scan, time interval between scan and stool sample, and average motion displacement were entered as regressors of no interest in each. Monte Carlo simulations were conducted in AFNI's 3dClustSim program, taking into account the average smoothness of the data set, with voxel thresholds of p = .01 and α = 0.05. Output cluster sizes were used to correct for multiple comparisons (k > 721 voxels). Based on the animal literature, which has highlighted the amygdala, hippocampus, and prefrontal cortex as central nodes involved in microbiome brain interactions across development (Hoban et al., Reference Hoban, Stilling, Ryan, Shanahan, Dinan, Claesson and Cryan2016, Reference Hoban, Stilling, Moloney, Moloney, Shanahan, Dinan and Clarke2017, Reference Hoban, Stilling, Moloney, Shanahan, Dinan, Clarke and Cryan2018; Ogbonnaya et al., Reference Ogbonnaya, Clarke, Shanahan, Dinan, Cryan and O'Leary2015), and the human literature implicating those same brain regions in anxiety and threat responsivity (reviewed in Callaghan & Tottenham, Reference Callaghan and Tottenham2015), we identified bilateral amygdala, hippocampus, and medial prefrontal cortex a priori as regions of interest. We decided to perform a restricted search within those three brain regions if no microbiome associations with reactivity were observed when cluster correction occurred at the whole-brain level (this restricted search was used when examining associations with Bacteroides only). Using the Harvard Oxford Cortical and subcortical atlas in FSL, we made masks of the frontal medial cortex (20% threshold), bilateral amygdala (50% threshold), and bilateral hippocampus (50% threshold). We then binarized and combined the three masks to set a threshold for multiple comparison correction within the restricted search regions (5735 voxels in combined mask; p < .01, α = 0.05, k > 100 voxels). Across whole-brain and restricted searches, a relatively lenient voxel threshold of p < .01 (rather than p < .001) was chosen in this small proof-of-concept sample as the primary aim was to demonstrate feasibility of the approach and drive future hypotheses.
Statistical tests
We assessed group differences in the community composition of microbes (i.e., alpha diversity) through ANOVA with the bacterial community measure (e.g., Shannon's diversity) as the outcome and caregiving history as the factor of interest, controlling for age at the time of stool sample, sex, and gastrointestinal distress and digestive issues factor scores.
Group differences in the community composition of bacteria between COMP and EA youth
Across all samples, bacteria from the genus Bacteroides were the most abundantly expressed, and the level of Bacteroides did not differ between groups of individuals with different early caregiving histories (controlling for age, sex, and gastrointestinal symptom factor scores), F (1, 10) = 1.49, p = .251, ηp2 = .13 (Figure 6a shows the relative abundance of the top 30 genera split by caregiving group). When considering the community structure of the microbiome, individuals that were EA had lower counts of bacteria: alpha diversity, observed richness, F (1, 10) = 7.18, p = .023, ηp2 = .42, but did not differ in the relative abundance of those bacteria: alpha diversity, Shannon index, F (1, 10) = 2.37, p = .155, ηp2 = .02. Group differences also emerged when comparing beta diversity. When rare bacteria were allowed to be influential (unweighted Unifrac distances), caregiving group significantly accounted for variation in bacterial distance matrices between individuals, adonis test, F (1, 14) = 1.72, p = .034, r 2 = .11, ηp2 = .11 (see Figure 6b for the unweighted phylogenetic tree), but that association was reduced to trend level when dominant bacteria were weighted: weighted Unifrac distances, adonis test, F (1, 14) = 2.80, p = .060, ηp2 = .17.
Figure 6. (a) Relative abundance of the 30 most dominant bacterial genera in early adversity (EA, n = 8) exposed and comparison (COMP, n = 8) children. (b) Heatmap of phylogenetic relatedness of samples based on unweighted Unifrac distances. The y-axis represents group membership with EA youth in pink and COMP youth in blue (i.e., each row represents an individual subject in the analysis), sorted by unweighted Unifrac distance. The bottom x-axis represents the individual bacterial groups (genera), and the top x-axis illustrates their sorting by Euclidean distance. The colors in the map represent the relative abundance of each bacterial group with warmer colors indicating higher relative abundance. (c) Caregiving group means of bacteria from order Clostridiales family Lachnospiraceae genus unknown and order Clostridiales family unknown genus unknown (bar) with individual subject values overlaid in gray transparent circles (transparency and random jittering of dots along the x-axis is used so that overlapping values can be viewed).
The LEfSE analysis indicated 10 bacterial genera as potential biomarkers for caregiving group membership, 8 of which were affected by outlying values for a single subject (>3 SD from the mean value) and were excluded from further analyses. The two biomarkers that were not affected by outlying values were from the order Clostridiales one of which was an unknown genus in the family Lachnospiraceae, and the other of which was an unknown genus and family. Both of these biomarkers were higher in children from the COMP than EA groups (see Figure 6c).
Brain–bacteria associations: General approach using genus Bacteroides
To establish associations between gastrointestinal bacteria and functional reactivity within brain networks implicated in anxiety, a whole-brain analysis was conducted on data from a task in which participants were looking at fear faces (relative to an implicit baseline). Individual subject Bacteroides levels were used as the regressor of interest (controlling for participant caregiving group, age at scan, the interval between scan and stool sample, and participant average motion after volume censoring). When correcting at the whole-brain level, Bacteroides was not associated with any reactivity clusters. Hence, we performed the restricted search within the bilateral amygdala and hippocampus, and frontal medial cortex, performing cluster correction within these regions only (resulting in a more liberal threshold). This approach revealed a small cluster of Bacteroides-associated reactivity in the medial prefrontal cortex (p < .01, restricted search correction; Figure 7a).
Figure 7. Whole brain statistical maps of thresholded reactivity to fear faces over implicit baseline in whole brain regression with (a) Bacteroides or (b) Lachnospiraceae as the regressor of interest. Warm colors indicate a positive association between Bacteroides or Lachnospiraceae and activity to fear faces, whereas cold colors indicate a negative association. (a) Reactivity in the medial prefrontal cortex (mPFC). (b) Reactivity in the left lateral prefrontal cortex (lPFC), mPFC, precuneus, and posterior cingulate cortex (PCC)/cerebellum. (c) Individual subject statistical maps masked with the group cluster map, showing unthresholded reactivity in left lPFC and mPFC, with comparison individuals featured in the top row and early adversity exposed individuals featured in the bottom row. Individual subject level of Lachnospiraceae is written underneath each subject's statistical map with high Lachnospiraceae levels to the left and low Lachnospiraceae levels to the right.
Brain–bacteria associations: Biomarker approach using genus Lachnospiraceae
To establish associations between gastrointestinal bacteria that acted as a biomarker for caregiving history, we performed another whole-brain analysis (fear > implicit baseline), but this time using Lachnospiraceae as the regressor of interest. Lachnospiraceae was selected for this analysis, as it was the biomarker with the greatest amount of taxonomic information available. Individual levels of Lachnospiraceae exhibited positive associations with the left lateral prefrontal cortex (PFC), medial PFC (mPFC), precuneus/cerebellum, and negative correlation with the post central gyrus; all p <.01, whole brain corrected; Figure 7b–c; see Table 4 for peak and center of mass coordinates for the whole-brain thresholded maps).
Table 4. Voxel number, center of mass, and peak coordinates (in the X, Y, Z directions) for clusters passing significance threshold for the group analyses associating brain reactivity to fear faces with individual Bacteroides and Lachnospiraceae levels
Note: mPFC, medial prefrontal cortex. PFC, prefrontal cortex. N = 16.
Assessment of potential confounding variables
Diet
A bivariate correlation was performed between current diet variables (average protein, carbohydrate, and fat intake) and Bacteroides, Lachnospiraceae, and extracted reactivity estimates (i.e., β values) from the clusters associated with bacteria (Bacteroides–mPFC; Lachnospiraceae–mPFC, left lateral PFC, posterior cingulate cortex, and precuneus; see Table 5 for correlation matrix). Average carbohydrate consumption was positively correlated with the level of Bacteroides (r = .552, p = .01). Entering diet variables into the regression between bacteria and brain did not change any reported associations (see Table 6 for statistics).
Table 5. Correlation matrix of diet variables (proportion of protein, carbohydrate, and fat in the diet) with bacterial genera (Bacteroides and Lachnospiraceae), and each of the significant clusters of activity in the brain, N = 16
*Correlation is significant at the .05 level (2 sided). **Correlation is significant at the .01 level (2 sided).
Table 6. Statistical table of regressions between bacteria (Bacteroides or Lachnospiraceae) and significant clusters of reactivity in the brain when controlling for potential confounding variables: dietary protein, carbohydrate, fat, child country of origin, child IQ, and child sex
Participant demographics
Entering participant demographic variables (country of origin, IQ, and sex) into the regression did not change the association of Bacteroides with mPFC reactivity (see Table 6 for statistics). Similarly, those demographic variables did not change the relationship between Lachnospiraceae and mPFC, left lateral PFC, precuneus, and posterior cingulate cortex (see Table 6 for statistics).
Discussion
Several studies have noted the strong association between gastrointestinal and mental health (Lee et al., Reference Lee, Wu, Ma, Tsang, Guo and Sung2009; Mak et al., Reference Mak, Wu, Chan, Chan, Sung and Lee2012). Moreover, the association of gastrointestinal and mental health problems with experiences of early adversity is well established (at least in adults; Bradford et al., Reference Bradford, Shih, Videlock, Presson, Naliboff, Mayer and Chang2012; Chitkara et al., Reference Chitkara, van Tilburg, Blois-Martin and Whitehead2008; Park et al., Reference Park, Videlock, Shih, Presson, Mayer and Chang2016). Stress programming effects at the level of the gastrointestinal microbiome and brain have been posited as potential biological mechanisms through which such associations emerge, with rodent models demonstrating the strongest evidence to this effect (Callaghan, Reference Callaghan2017; Foster & Neufeld, Reference Foster and Neufeld2013; O'Mahony, Hyland, Dinan, & Cryan, Reference O'Mahony, Hyland, Dinan and Cryan2011). However, while it is probable that early adversity effects on gastrointestinal and mental health manifest early in life, much of the research on this topic has been performed in adults, relying on retrospective reports of childhood adversity experiences, which lack the accuracy of prospective designs and have been shown to often reflect adult mental health state (Reuben et al., Reference Reuben, Moffitt, Caspi, Belsky, Harrington, Schroeder and Danese2016; Susser & Widom, Reference Susser and Widom2012). In addition, such studies lack the ability to isolate the effects of stress on gastrointestinal health to the early life period.
In the current study, we examined associations between adversity exposure in infancy and gastrointestinal symptoms across childhood and adolescence in a unique population where the end date of adversity is known: youth who received foster or institutional care before international adoption. Similar to past research in rodents (Yi et al., Reference Yi, Zhang, Sun, Zhou, Chen, Xuan and Xu2017), we observed that early adversity was associated with gastrointestinal symptoms in youth aged 3–18 years (with the largest effects observed in late childhood). In addition to those findings, we observed that GI distress at the Time 1 assessment was associated with concurrent anxiety, as well as with future anxiety (measured over a 5-year time frame), over and above the persistence of GI distress and baseline anxiety scores. Within the cross-sectional data (where the largest number of data points were observed) we also saw that GI distress mediated the association between early adversity and elevated anxiety symptoms, demonstrating that the indirect pathway from caregiving adversity to elevated anxiety through concurrently increased GI distress explains a significant amount of the variance in the adversity–anxiety association. Moreover, we obtained proof of concept that the gastrointestinal microbiome was altered by early experiences of adverse caregiving, and that such microbial variation was associated with brain reactivity within emotion networks in the brain: the prefrontal cortex, posterior cingulate cortex, and precuneus.
Research in adults has shown that early adversity is associated with a constellation of functional gastrointestinal complaints from visceral pain and feelings of sickness to constipation, diarrhea, and bloating (Wu, Reference Wu2012). In line with that past research, in Study 1 we observed adversity associations with a range of GI symptoms, including those that fell into the “gastrointestinal distress” factor (aches and pains, funny feelings in the stomach, and nausea), and those that fell into the “digestive distress” factor (i.e., vomiting and constipation). In Study 2 we also saw that early adverse caregiving history was associated with the frequency of diarrhea. Nonetheless, in Study 1 only a subset of these GI symptoms (“aches and pains in the stomach,” “funny feelings in the stomach,” and “nausea”) were associated with anxiety (note: items assessing abdominal sensations were excluded from the anxiety assessment to avoid inflated correlations), both concurrently and in the future, suggesting that there is specificity in the GI distress–anxiety association. We suggest that this specific association reflects the bidirectional nature of anxiety and GI distress correlations already reported in the literature. That is, GI distress is likely a correlate, a symptom, and in some cases a cause of anxiety. While the current study did not aim to differentiate between these possibilities, the longitudinal association between GI distress and anxiety supports the use of GI distress as a predictor of future psychopathology, particularly in pediatric primary care settings in which GI distress is a common reason for presentation. The fact that GI distress predicted future anxiety symptoms over and above concurrent anxiety, as well as early adverse caregiving, furthers the proposed clinical utility of measuring GI symptoms, as it suggests that intervening at the level of the GI system could potentially affect future anxiety symptoms, regardless of their etiology. Moreover, the finding that GI distress was a mediator of the adversity–anxiety relationship (cross-sectionally) suggests that addressing GI symptoms may be especially important in treating elevated anxiety within adversity-exposed populations. Although EA youth as a group had higher anxiety levels than comparisons, these differences were dampened in those experiencing lower GI distress, providing a diagnostic marker that may be more readily available to clinicians as well as health workers outside of the psychology/psychiatry field.
While there are several possible mechanisms through which caregiving adversity, GI distress, and anxiety may be associated, one candidate with strong clinical promise is the GI microbiome. In Study 2, we obtained proof-of-concept that adversity is associated with altered microbial patterns in developing youth, and that bacteria–brain associations are observable within task-based functional reactivity (emotional faces probe task) particularly in brain regions well known to be implicated in emotional functioning (prefrontal cortex and posterior cingulate cortex; Bush, Luu, & Posner, Reference Bush, Luu and Posner2000; Hiser & Koenigs, Reference Hiser and Koenigs2018). Moreover, by using adversity-associated (Lachnospiraceae) and adversity-independent (Bacteroides) bacteria, we observed microbiome–brain associations that may be a biomarker or fingerprint of adversity exposure, as well as those that may be somewhat experience independent (i.e., that transcend adversity effects). Rodent models have consistently demonstrated a causal link between early microbial communities and typical brain development, anxiety, and mood (e.g., Clarke et al., Reference Clarke, Grenham, Scully, Fitzgerald, Moloney, Shanahan and Cryan2013; Collins, Kassam, & Bercik, Reference Collins, Kassam and Bercik2013; Hoban et al., Reference Hoban, Stilling, Ryan, Shanahan, Dinan, Claesson and Cryan2016; Kelly et al., Reference Kelly, Borre, O'Brien, Patterson, El Aidy, Deane and Moloney2016), and both rodent and monkey models have shown that caregiving stress alters such microbial communities (Bailey & Coe, Reference Bailey and Coe1999; O'Mahony et al., Reference O'Mahony, Hyland, Dinan and Cryan2011). In humans, although previous work has demonstrated correlations between the microbiome and neural phenotypes in adult females and children without a history of adversity (brain structure; Carlson et al., Reference Carlson, Xia, Azcarate-Peril, Goldman, Ahn, Styner and Knickmeyer2018), here we observed brain–microbiome relationships following adverse caregiving in brain regions with known importance for anxiety and mental health and that have been shown in rodent models to exhibit developmental programming in response to the microbiome (e.g., prefrontal cortex). Considering that microbial populations can be easily manipulated (e.g., via probiotics and diet), the association between adversity, microbiome, brain, and anxiety seen here in developing populations has potential clinical value, with significant intervention implications that warrant further investigation.
The fact that early experiences of caregiving adversity increase the risk for both GI and mental health problems, as well as being associated with changes to the GI microbiome, suggests that the timing of this exposure (infancy to early childhood) coincides with the sensitive period for the maturation of each of these systems. While such “sensitive periods” have causal evidence in animal models (e.g., Bailey & Coe, Reference Bailey and Coe1999; Heijtz et al., Reference Heijtz, Wang, Anuar, Qian, Björkholm, Samuelsson and Pettersson2011; Neufeld, Kang, Bienenstock, & Foster, Reference Neufeld, Kang, Bienenstock and Foster2011), there are few human populations in which the timing of adversity can be isolated to the early life period, making investigation of “sensitive periods” exceedingly difficult. Children who were internationally adopted following early foster or institutional care represent a unique group in which the end date of adversity is known (i.e., the date of adoption), providing a rare opportunity to identify whether such “sensitive periods” operate in human GI and psychosocial health domains. In our study, almost all of the EA children experienced institutional or foster care prior to 13 months of age, but the age at which they were adopted from those care settings was not associated with GI distress or digestive issues. These data suggest that time-limited exposure to caregiving adversities during the infancy–early childhood stages of life, but perhaps not the chronicity of such exposures within this period, may be important for development of the brain–gut–microbiome axis in humans. However, considering that 90% of our EA population was adopted by age 5 years (with median age of adoption at 13 months old), it is also possible that longer durations of foster/institutional care do contribute to enhanced risk for GI health issues that were not captured here. Ascertaining such developmental timing effects will assist in the prevention and targeted treatment of adversity-associated pathology. Moreover, assessing the influence of adversity timing on the composition of GI bacteria will be critically important to examine in future studies with larger sample sizes.
There are some limitations to the current studies. In Study 1, we did not have a direct, validated measure of functional gastrointestinal disorders through which GI symptoms could be assessed. Nonetheless, the questions in the CBCL and RCADS-P do assess functional GI symptoms and were factorable, suggesting that they were assessing an underlying construct. Another limitation was the study design, in which only some children were scheduled for follow-up visits, resulting in a smaller sample size at Time 2 and Time 3 assessments, approximately 2 and 4 years after the Time 1 assessment, respectively. Thus, additional longitudinal studies should confirm the associations observed here. For Study 2, we are limited by the small sample size of the population, and delay between fMRI and microbiome data collection, which occurred as a result of using a convenience sample of children who had already contributed a usable fMRI scan during childhood, and who subsequently agreed to donate a fecal sample. Given some estimates indicating high microbiome stability after approximately 4 years of age (Yatsunenko et al., Reference Yatsunenko, Rey, Manary, Trehan, Dominguez-Bello, Contreras and Anokhin2012), it is possible that the fecal sample at time of collection was reflective of earlier bacterial communities. However, it is also possible that the developmental stability of the microbiome may change following stress exposure. These possibilities will need to be examined in future cohorts. In addition, adversity exposure in the current sample is confounded with country of origin, which could influence the microbiome, but participant country of origin was not associated with community composition of the gut, or with microbiome–brain reactivity patterns. Moreover, bacteria were associated with brain reactivity patterns in regions implicated in emotion regulation (e.g., PFC), which are not expected to vary by country of origin. As half of the sample used to examine brain–microbiome associations were previously institutionalized, we have no preadoption information on the incidence of several early life factors known to influence the microbiome, such as mode of birth (vaginal vs. cesarean), pre- and postnatal antibiotic usage, and breastmilk versus formula feeding. However, it is unlikely that any one infancy environment factor alone accounts for the association between caregiving group and microbiome composition reported here. In contrast, some postadoption factors, such as age at adoption, and diet around the time of stool sample, are known and (in future studies with larger sample sizes that have the power to examine individual differences) should be explored as predictors of microbiome composition, and in the case of diet, as a potential ameliorative factor in adversity-associated changes in bacterial composition. Another limitation relates to our amplification of the 16S region of ribosomal RNA to identify bacteria. As other components of the microbiome, including fungi, viruses, and parasites, do not possess the 16S gene, the conclusions of this study are limited to bacteria. Even with these limitations, the current data strongly support further investigation of GI distress in youth as a predictor of future psychopathology, as well as a more general investigation of the gastrointestinal system in foster care, previously institutionalized, and otherwise early adversity-exposed populations of developing youth. In addition, these data support future investigations of the microbiome–brain functional relationship as a pathway through which adversity may get under the skin to affect the development of emotion neurobiology.
Acknowledgments
The authors thank Dima Amso and Andrew Lynn for help with setting up the fMRI processing pipeline within AFNI and FSL. We thank RTL Genomics for 16s rRNA amplification and sequencing.
Financial support
This research was supported by the National Institute of Mental Health Grant R01MH091864 (to N.T.), Dana Foundation (to N.T.), National Institute of Mental Health Grant 1K99MH113821-01 (to B.C.), National Health and Medical Research Council Early Career Fellowship 1091571 (to B.C.), and American Australian Association Fellowship (to B.C.).
