Longitudinal Trajectories of Hair Cortisol: Hypothalamic-Pituitary-Adrenal Axis Dysfunction in Early Childhood

The objective of this study was to examine if longitudinal trajectories of hair cortisol concentrations (HCC) measured at two or three yearly time points can identify 1-3 year old children at risk for altered hypothalamic-pituitary-adrenal (HPA)-axis function due to early life stress (ELS). HCC was measured (N = 575) in 265 children using a validated enzyme-linked immunosorbent assay. Hair was sampled in Clinic Visits (CV) centered at years 1, 2, and 3 (n = 45); 1 and 2 (n = 98); 1 and 3 (n = 27); 2 and 3 (n = 95). Log-transformed HCC values were partitioned using latent class mixed models (LCMM) to minimize the Bayesian Information Criterion. Multivariable linear mixed effects models for ln-HCC as a function of fixed effects for age in months and random effects for participants (to account for repeated measures) were generated to identify the factors associated with class membership. Children in Class 1 (n = 69; 9% Black) evidenced declining ln-HCC across early childhood, whereas Class 2 members (n = 196; 43% Black) showed mixed trajectories. LCMM with only Class 2 members revealed Class 2A (n = 17, 82% Black) with sustained high ln-HCC and Class 2B (n = 179, 40% Blacks) with mixed ln-HCC profiles. Another LCMM limited to only Class 2B members revealed Class 2B1 (n = 65, 57% Black) with declining ln-HCC values (at higher ranges than Class 1), and Class 2B2 (n = 113, 30% Black) with sustained high ln-HCC values. Class 1 may represent hair cortisol trajectories associated with adaptive HPA-axis profiles, whereas 2A, 2B1, and 2B2 may represent allostatic load with dysregulated profiles of HPA-axis function in response to varying exposures to ELS. Sequential longitudinal hair cortisol measurements revealed the allostatic load associated with ELS and the potential for developing maladaptive or dysregulated HPA-axis function in early childhood.

The objective of this study was to examine if longitudinal trajectories of hair cortisol concentrations (HCC) measured at two or three yearly time points can identify 1-3 year old children at risk for altered hypothalamic-pituitary-adrenal (HPA)-axis function due to early life stress (ELS). HCC was measured (N = 575) in 265 children using a validated enzyme-linked immunosorbent assay. Hair was sampled in Clinic Visits (CV) centered at years 1, 2, and 3 (n = 45); 1 and 2 (n = 98); 1 and 3 (n = 27); 2 and 3 (n = 95). Logtransformed HCC values were partitioned using latent class mixed models (LCMM) to minimize the Bayesian Information Criterion. Multivariable linear mixed effects models for ln-HCC as a function of fixed effects for age in months and random effects for participants (to account for repeated measures) were generated to identify the factors associated with class membership. Children in Class 1 (n = 69; 9% Black) evidenced declining ln-HCC across early childhood, whereas Class 2 members (n = 196; 43% Black) showed mixed trajectories. LCMM with only Class 2 members revealed Class 2A (n = 17, 82% Black) with sustained high ln-HCC and Class 2B (n = 179, 40% Blacks) with mixed ln-HCC profiles. Another LCMM limited to only Class 2B members revealed Class 2B1 (n = 65, 57% Black) with declining ln-HCC values (at higher ranges than Class 1), and Class 2B2 (n = 113, 30% Black) with sustained high ln-HCC values. Class 1 may represent hair cortisol trajectories associated with adaptive HPA-axis profiles, whereas 2A, 2B1, and 2B2 may represent allostatic load with dysregulated profiles of HPA-axis function in response to varying exposures to ELS. Sequential longitudinal hair cortisol measurements revealed the allostatic load associated with ELS and the potential for developing maladaptive or dysregulated HPA-axis function in early childhood. Keywords

INTRODUCTION
The stress response forms a highly conserved regulatory system in all eukaryotic multicellular organisms, evolutionarily designed to cope with a broad range of stimuli that may threaten, or be perceived as threatening to their survival, growth, and/or dynamic equilibrium ("homeostasis"). The stress response includes the neuroendocrine, neuroimmune, and other systems and a central modulator of these systems is the hypothalamicpituitary-adrenal (HPA)-axis (1). Allostatic load is defined as the cumulative burden of exposures to repeated or chronic stressors, associated with deleterious consequences for growth and survival as well as longterm changes in the stress response system (2). Chronic or repetitive exposures to early life stress (ELS) alter HPA-axis regulation, associated with the early onset of chronic non-communicable diseases (NCDs) and psychopathology (3,4). The American Academy of Pediatrics (AAP) made an appeal for distinguishing adaptive, maladaptive, and toxic stress in children and declared that "Identifying children at high risk for toxic stress is the first step in providing targeted support for their parents and other caregivers" (5).
Cortisol is the hormonal end-product of the HPA-axis that regulates the endocrine and metabolic stress responses, brain and other organ development, and immune functions. Transitory changes in salivary or serum cortisol may reflect exposures to acute stress (6) but cannot capture cumulative exposures to chronic or repetitive stress. Cortisol binds to growing hair; thus, hair cortisol concentrations (HCC) can provide a summative measure of responses to stressful experiences over the previous 3-6 months (7,8). A recent review proposed that HCC may represent a measure of recent/acute stress, but this metaanalysis included studies from 16 animal species, ignored major differences between animal fur and human scalp hair, and analyzed studies with cross-sectional rather than longitudinal data (9). HCC in humans are positively correlated with repeated measures of salivary cortisol (8).
Adverse childhood experiences (ACEs) alter subsequent stress responses associated with HPA-axis dysregulation, particularly among Black adults (10)(11)(12). A dysregulated HPA-axis often shows persistent hyper-or hypo-responsiveness, associated with a reduced cortisol awakening response, flattened diurnal cortisol slope, higher evening cortisol levels, and dampened responses to corticotrophin (ACTH) or corticotrophin-releasing hormone (CRH) (13). Preclinical data illustrate progressive phases of altered HPA-axis function ultimately leading to HPAaxis dysregulation. However, a substantial gap exists between preclinical and clinical data because longitudinal human studies have not mapped these phases of HPA-axis function, particularly among young children. Heightened neuroplasticity in children aged 0-3 years makes them vulnerable to long-term changes in the HPA-axis regulation following exposure(s) to ELS, linked epidemiologically with poor physical and mental health across the lifespan (1,5,(13)(14)(15).
Our overarching goal is to measure the long-term effects of ELS on the HPA axis. We have previously published novel methods for measuring HCC (16) and socioeconomic adversity (17), and also identified the demographic and psychosocial factors related to HCC in a cross-sectional population of preschool children (18,19). Black children showed higher HCC than White/other children because of exposure to psychosocial stressors and daily microaggressions; these differences persisted even after adjusting for socioeconomic factors. HCC values measured at earlier ages were significant predictors of HCC at later ages (18), therefore we sought to examine longitudinal HCC trajectories using latent class mixed models (LCMM). Longitudinal analyses may be more likely to identify individual children with altered states of HPA axis regulation during early childhood, and therefore, may be more vulnerable to the longterm effects of ELS on the subsequent cognitive, behavioral, educational, psychosocial, physical health and psychiatric outcomes. If these children can be identified, they would be eligible for potential therapeutic interventions to prevent or ameliorate these longterm effects.
Compared to cluster analyses, structural equation models, or generalized estimating equation models, LCMM is more suitable for investigating developmental trajectories because it is robust to missing values and categorical variables, and is specifically designed to partition groups of individuals with similar profiles in their longitudinal data (20,21). We hypothesized that latent patterns in HCC values measured longitudinally during early childhood will identify the groups of children at higher risk for altered HPA-axis function. By comparing these groups of children, we sought to identify the maternal and psychosocial factors that may be related to these unique longitudinal HCC profiles. We speculate whether longitudinal measurements of HCC values may identify the children at risk for HPAaxis dysregulation and its longterm consequences, thereby identifying participants for future therapeutic interventions to overcome the pernicious and persistent effects of ELS across the lifespan.

CANDLE Cohort
Following approval from the Institutional Review Board (IRB) at the University of Tennessee Health Sciences Center, participants or legal guardians gave informed consent for participation in the conditions affecting neurocognitive development and learning in Early childhood (CANDLE) study as well as hair sampling from the CANDLE children during yearly clinic visits. Pregnant women (N = 1,503) were enrolled during their second trimester between 2006 and 2011 for a longitudinal birth cohort study. CANDLE study participants were described previously (17)(18)(19)22) with detailed descriptions of the instruments and survey items for characterizing children and their mothers (17-19, 22, 23). Characteristics of the children and mothers who gave consent for hair sampling are listed in Supplementary Table A; the variables measured from clinic visits (CV) centered at ages 1, 2, 3 years are listed in Supplementary Table B. Parents declined hair sampling for children on some annual visits due to cosmetic reasons, e.g., preference for buzz cuts for boys or elaborate hair designs for girls.

Hair Collection, Processing, and HCC Quantification
Hair samples were collected from the posterior vertex by trimming hair close to the scalp surface. Hair samples (10 mg) were cut to powder consistency and proteins extracted by sequential methanol (52 • C for ∼15 h) acetone (10 min at room temperature) extractions, repeated once. Supernatants were air evaporated at 4 • C; the protein residue reconstituted in 70 mcl of phosphate buffered saline (PBS, pH 7.8) per 10 mg hair, kept at 4 • C for immediate use or at −20 • C for long-term storage. The Pierce Micro BCATM Protein Reagent Kit was used to determine total protein content and the Alpco Cortisol (Saliva) ELISA assay was used to quantify HCC. An Epoch BioTek plate reader with associated software Gen 5 Biotech (Version 1.11, Winooski, VT, USA) determines total protein and cortisol levels (ng/mg hair) (16). HCC assays had intra-assay and inter-assay coefficients of variation at <4 and <8%, respectively.

Statistical Methods
The Stanford University IRB approved these data analyses. Data analyses were conducted on CANDLE subjects providing hair samples at 2 or 3 clinic visits (CV) with age at CV1 = 11-18 months; CV2 = 23-30 months; and CV3 = 35-42 months. Longitudinal HCC values were examined visually via scatterplot and log-transformed HCC values (ln-HCC) were partitioned into different trajectories using LCMM (20,21). The model was framed as a linear mixed effects model for ln-HCC as a function of fixed effects for age in months and age in months squared, coupled with random effects for CANDLE participants (to account for repeated measures from the same subject) and latent classes identified by age in months. The number of latent classes were prespecified to 2, 3, 4, or 5, and the pre-specified value that minimized the Bayesian Information Criterion (BIC) was selected (Supplementary Table C). Vrieze states that "BIC has α = log(N), where N is the number of observations contributing to the sum in the likelihood equation. The number of model parameters κ is taken to be a direct indicator of model complexity; the more parameters the more flexible and complex the model" (24).
The same process of fitting LCMM models within subgroups was repeated until consistent ln-HCC patterns remained within each identified group. This approach is commonly recommended in various statistical texts (25)(26)(27) and in a reference text: "Applied Latent Class Analysis" edited by Hagenaars and McCutcheon (28). These authors state that LCMM analyses may be performed sequentially or simultaneously within the subclasses identified from an initial model. Sequential LCMM modeling procedures within stratified subpopulations are also reported in several previous studies (29)(30)(31).
Associations between candidate variables and the LCMMidentified classes were examined via univariate and multivariate tests at each clinic visit (Supplementary Tables D-F). Associations between each candidate variable and group membership were tested using Fisher's exact test for categorical variables and Wilcoxon-rank sum for continuous variables.
Candidate variables with p ≤ 0.1 were entered into logistic regression models using backwards selection to eliminate the variable with the largest p-value at each step until all the remaining variables were p < 0.1. Variable selection was implemented in the R package "rms" (Regression Modeling Strategies; R package version 5.1-3.1.). All statistical analyses were conducted using R version 3.6.1 (32). Study integrity was assured by completing the STROBE checklist for observational studies (Supplementary Table G

RESULTS
Visual inspection of the ln-HCC values plotted by age in months (Supplementary Figure 1) illustrated the complexity of longitudinal ln-HCC trajectories for individual children with varying exposures to contextual stressors over time.

DISCUSSION
Karlen et al. found that sequential HCC measurements in children aged 1, 3, 5, and 8 years had significant linear correlations across ages, but they did not examine longitudinal trajectories in individual children or relate these effects to ELS (33). Our cross-sectional analyses also revealed that higher ln-HCC in the preschool children at earlier ages significantly predicted their ln-HCC values at subsequent ages (18). Consequently, we conducted longitudinal analyses to reveal latent patterns in their ln-HCC trajectories and to examine the maternal and childrelated factors that may determine membership in these latent classes. We present the first-ever longitudinal analyses of hair cortisol from a geographically defined and racially diverse urban/suburban population of children aged 1-3 years, identifying latent classes based on their ln-HCC trajectories. We previously reported cross-sectional analysis of ln-HCC in this age group, showing that Black children were exposed to greater adversity (17) leading to chronic stress (18,19,22). Here, we extend our studies to explore latent patterns in ln-HCC trajectories from the children with longitudinal data. Membership in latent classes based on ln-HCC trajectories defined groups of children with progressive profiles of altered HPA-axis function, presumably resulting from increasing allostatic load as suggested by maternal and psychosocial factors as well as their physical and behavioral outcomes (14,34).
The typically developing human HPA-axis matures by 4 years, but it remains malleable before age 4. Therefore, it is important to determine the trajectories indicative of HPA-axis function in earlier developmental windows, when the set-points for longterm HPA-axis regulation are vulnerable to contextual stressors in early life (15,34,35), predisposing children to progress from adaptive, to maladaptive, and to toxic stress responses (5). Adaptive stress responses promote mastery and help build a child's resilience, maladaptive stress responses may interrupt these developmental processes, whereas toxic stress responses disable resilience and problem-solving skills, promoting passive submission or self-injurious/risk-taking behaviors (1,4,36).
Children assigned to LCMM Class 1 had decreasing ln-HCC across ages 1-3 years (Figure 1), following the overall patterns of HPA-axis function expected in normal healthy children (18).
Karlen et al. confirmed a pattern of decreasing HCC during early childhood (33), which was not seen in cross-sectional studies with sparse sampling from the preschool age groups (37). Given their favorable demographic and psychosocial factors, we propose that Class 1 identifies those children with an adaptive profile of HPA-axis function, easily capable of maintaining homeostasis after exposure to contextual stressors. Class 1 children had fewer social-emotional problems in infancy and manifested fewer internalizing or externalizing behaviors than those in Class 2. Other studies also found impaired social-emotional development associated with ELS (38,39), further linked with the poorer cognitive and emotional outcomes in older children and adolescents, and psychopathology in adults (4,40).
Subsequent partitioning of Class 2 identified the children with a hyperresponsive HPA-axis in Class 2A compared to 2B, showing high ln-HCC at age 1 year and further increasing values across early childhood rather than the expected decline in typically developing children (Figure 2). Most of these children were Black, had been exposed to maternal psychosocial risk factors, possibly nutrition-related intrauterine growth restriction at birth (41), and were more likely to develop social-emotional problems, internalizing and externalizing behaviors in infancy, as well as problems with aggression and attention in later childhood (38,39). HPA-axis hyperreactivity follows failure of the cortisolmediated negative feedback loop that regulates the secretion of CRH and ACTH (42). These children may have higher risks for developing long-term HPA-axis dysregulation, epidemiologically linked with early obesity, metabolic syndrome, adult-type NCDs (1,11,43) and psychopathology (4,40,44). Another study measuring salivary cortisol levels revealed that "Race differences in cortisol were not present at 4 months but emerged at 12 months of age, with Black infants having higher cortisol. . . suggesting that racial discrimination is already "under the skin" by 1 year of age" (45). We previously reported preliminary evidence for HPA-axis dysregulation in 1-year-old Black infants (19), and now present further evidence for HPA-axis dysfunction from longitudinal HCC trajectories in Class 2A mostly including 1-3 year-old Black children (82.4%), as compared to Class 1 membership mostly including White/other children (91.3%).
Class 2B1 members also showed the expected decline across early childhood, but their HCC values were higher than in Class 1 (Figure 3 vs. Figure 1) (18,33). This group contained similar numbers of White/other and Black children, perhaps in the compensatory phase of HPA-axis dysregulation with attempted adaptation to contextual stressors (46). Children in Class 2B1 were exposed to maternal and psychosocial risk factors, evidenced intra-uterine growth restriction at birth and, consistent with previous research, showed more internalizing behaviors at age 2 and anxiety at age 3 (40,47,48). Class 2B2, including 70% White/other children, showed a flat trajectory across early childhood-somewhat like Class 2A which mainly included Black children. The mothers of Class 2B2 children uniquely showed a hyperthymic temperament, which may lead to overly intrusive or overprotective parenting, associated with HPA-axis dysregulation in another study (49).
The interplay of health and psychosocial determinants occurs in a bidirectional fashion; therefore, correlates of HCC in toddlers may become prognosticators of later health outcomes. Given the patterns suggestive of HPA dysfunction in Classes 2A, 2B1, and 2B2, we speculate that children experiencing repeated contextual stressors may develop maladaptive HPA-axis responses that were associated with depression or posttraumatic stress disorder (PTSD) (48,(50)(51)(52)(53). Lower HCC values were also noted in individuals with PTSD as compared to normal or acute stress groups (54). These findings appear to align with the theoretical phases of HPA-axis dysregulation (Figure 4), in which cumulative allostatic loads related to ELS may progressively impair HPA-axis function, as illustrated in Class 2 and its subclasses. The current criteria for HPA-axis dysregulation are based on plasma and salivary cortisol levels, which are measures of acute stress which can be altered by diurnal cycles, high state reactivity, food intake, pulsatile secretion patterns, superimposed on the robust changes across age, sex, reproductive cycles, and pubertal stages (12,13,55,56). Because prolonged or chronic stressors are more likely precursors of HPA-axis dysfunction than acute stressors, putative criteria for identifying HPA-axis dysregulation should also be developed from HCC, which serves as a summative measure of chronic stress (57). Based on our data, we present a starting point for the HCC ranges in 1-3 year-old children that could be used for developing such criteria (Figure 4).
The results of these analyses must be interpreted in the light of four limitations. First, as in other cohort studies of early stress (19,58,59), the generalizability of these results is limited by the differences in demographic characteristics between children with HCC measurements in the CANDLE study and those of other children (17). Second, HCC data from all clinic visits were only present in 45 children, with other children contributing measurements at any combination of two timepoints. This highlights the challenges of implementing longitudinal studies in a cohort of inner-city children and obtaining consent for hair sampling at each of their clinic visits. Though missing data may affect the confidence intervals of the HCC trajectories identified, our use of LCMM modeling was specifically suited for handling missing data in our subjects (21,25,26,28). Third, because these are exploratory analyses to identify the characteristics of LCMM class membership, we did not adjust the p-values for multiple testing. Therefore, the explanatory factors identified in our univariate or logistic regression analyses must be confirmed in larger cohort studies. Fourth, backwards selection is one of many variable selection techniques that can be used, each with its own limitations. Whereas elastic net may be superior for minimizing prediction error, our choice of backwards selection was designed to identify those variables that retain significant associations with class membership. Our sensitivity analyses included variable selection using LASSO and elastic net, and found that similar variables were selected consistently across all methods. However, our use of logistic regression modeling with backwards selection allows easier interpretability of the preliminary variables that characterize latent class membership.
We present an objective method for longitudinally measuring chronic stress in a cohort of children aged 1-3 years, identify factors associated with adaptive and altered forms of HPA function, and suggest an approach to interpreting HCC values in this cohort. These results expand on our previous publications on HCC and social adversity in preschool children, clearly demarcating racial and socioeconomic disparities in urban/suburban populations (17,18). Longitudinal measures of HCC in young children may provide an objective tool for classifying the early developmental trajectories in immature stress pathways and could also be used assess the efficacy of therapeutic interventions. Taken together, our linked approaches for measuring early life adversity and using HCC trajectories to identify associated risk factors in early childhood may ultimately enable monitoring the effects of policy reform developed to address the widening health inequities between racial, socioeconomic, immigrant, and other demographic subgroups.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: Data for the CANDLE

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Institutional Review Board (IRB) at the University of Tennessee Health Sciences Center. Written informed consent to participate in this study was provided by the participants' legal guardian/next of kin.

AUTHOR CONTRIBUTIONS
CR contributed to the study conceptualization and methodology, data curation, analyses and interpretation, drafting the initial manuscript, making critical revisions, and final approval of the submitted manuscript. JR contributed to statistical design and methodology, data curation, analyses and interpretation, drafting the initial manuscript, making critical revisions, and final approval of the submitted manuscript. J-MR and VC contributed to data curation, analyses and interpretation, making critical revisions, and final approval of the submitted manuscript. MR contributed to data curation, analyses and interpretation, drafting the manuscript, making critical revisions, and final approval of the submitted manuscript. KA obtained grant funding and contributed to study conceptualization, statistical design and methodology, data curation, analyses and interpretation, supervision of the team, drafting of the initial manuscript, making critical revisions, and final approval of the submitted manuscript. All authors contributed to the article and approved the submitted version.

FUNDING
The current studies received funding from the Maternal and Child Health Research Institute at Stanford University (P.I. KA) and the Eunice Kennedy Shriver National Institute for Child Health and Human Development (R01 HD099296, P.I. KA). Study sponsors had no role in the design and conduct of the study; the collection, management, analysis, or interpretation of the data; the preparation, review, approval, or decision to publish this manuscript.

ACKNOWLEDGMENTS
Grant Somes, PhD [University of Tennessee Health Science Center (UTHSC)] conceived of the CANDLE study and obtained initial funding, Frances A. Tylavsky, PhD was the CANDLE Study PI during the hair sampling, and W. Alex Mason, PhD is the current CANDLE Study PI. We gratefully acknowledge the participants, particularly the mothers who consented to participate in the CANDLE study; the cognitive examiners from the UTHSC Boling Center for Developmental Disabilities