Perinatal depression (PND) is a common, morbid condition affecting >400,000 mother-child dyads in the United States each year [1]. In addition to the direct effect of PND on mothers’ health and well-being, it presents challenges for early child development, including reduced maternal sensitive caregiving [2]; impaired parent-child attachment relationships; and compromised child emotion regulation abilities [3,4], socioemotional development [5,6], and executive functioning (EF) [7,8]. These effects may occur even when mothers receive medication or psychotherapy for depressed mood [9], and they may extend well into adolescence and early adulthood [10,11]. Research has implicated the dysregulation of the oxytocin system [12] and the hypothalamic-pituitary-adrenal (HPA) axis [13] in the context of maternal depression as potential contributors to individual and parent-child relational dysfunction that increases the risk of nonoptimal child outcomes. The Mood, Mother and Child (MMC) study (funded by the National Institutes of Health; grant R01HD093901; grant review provided in Multimedia Appendix 1) addresses specific aims regarding the complex roles oxytocin functioning plays within the mother-child dyad, including the role of maternal oxytocin and HPA axis functioning as predictors of postpartum depression and anxiety symptoms (PDAS) and parenting behaviors and the role of child oxytocin functioning as both an outcome, mediator, and moderator of exposure to maternal depression and maladaptive parenting behaviors. This paper provides a brief synopsis of these goals and a detailed description of the MMC study protocol.

Both the severity and chronicity of maternal depression are associated with child developmental outcomes [14,15], and this association is (in part) mediated by dysregulated dyadic interactions between mother and child [16]. Symptomatic mothers are more likely to display withdrawn or intrusive behaviors when interacting with their infants in the first postpartum year [17,18], and these associations persist through early childhood, with the lowest levels of sensitive caregiving observed among mothers with chronic depression [19,20]. Building on previous research demonstrating the developmental risk associated with exposure to maternal depression, the MMC study will examine PDAS and associated parenting behaviors across early childhood as predictors of a variety of biopsychosocial outcomes in middle childhood, including EF, general cognitive development, and socioemotional adjustment. However, far less is known about the biopsychosocial processes underlying the manifestation of depressive symptoms in mothers and the transmission of risk in offspring. To address this gap, the MMC study will examine the roles of maternal and child oxytocin functioning, heart rate variability (HRV), and HPA axis functioning as potential sources of psychobiological risk and resilience.

A subset of systems contributing to the biopsychosocial process underlying experiences of PND may include oxytocin and its effects on HPA functioning. Research suggests that PND is associated with dysregulation of the HPA axis [21] and that this association may (in part) be driven by dysregulated oxytocin functioning [22]. Previous studies on this topic—including research by the authors of this paper—have been mixed, including research demonstrating (1) no association between oxytocin during breastfeeding and PND [23], (2) an inverse association between oxytocin during breastfeeding and maternal symptoms of depression and anxiety at 2 months post partum [12], and (3) that oxytocin levels differentially predicted HPA axis reactivity in women with or without mood symptoms [13]. There are also candidate gene studies reporting correlations between oxytocin-related genetic variants and maternal depression [24] as well as oxytocin-related variants that predict developmental risk factors associated with maternal depression, such as stress reactivity [25,26], empathetic functioning [27], and parenting behaviors [28-30]. In addition, there is mixed evidence that exogenous administration of oxytocin can alter parenting behaviors [31-33] (for a systematic review, see the study by Szymanska et al [34]), although we are unaware of studies examining oxytocin administration and depressive symptoms. Furthermore, research on oxytocin and the psychobiology of PND has largely been conducted on mothers of infants and toddlers, with little empiricism on these topics in subsequent years of development. The MMC study will examine oxytocin functioning, oxytocin-related genetic variants, and HPA axis functioning as predictors of PDAS and parenting behaviors. This will be done using both a passive longitudinal design and a placebo-blinded trial examining the effects of exogenous administration of oxytocin on HPA axis functioning for mothers with and without depression diagnoses.

There is mixed evidence from candidate gene studies of oxytocin-related genetic variants functioning as moderators of the associations between adverse developmental exposures and subsequent child outcomes. For example, oxytocin receptor (OXTR) rs53576 G allele carriers who had experienced childhood maltreatment reported more depressive symptomatology compared with AA homozygotes, suggesting that the G allele confers increased vulnerability to stressors [35]; however, contradictory findings have reported that A allele carriage increased intergenerational transmission of maternal depression to child depressive symptomatology [36]. Meanwhile, a study of OXTR rs2254298 suggested that carriers of the A allele have increased developmental risk of anxiety and depression in the context of exposure to early-life adversity [37,38]. Finally, several studies indicate that these associations may vary by child sex [39,40] in the form of both sex-specific diathesis-stress and differential susceptibility models of gene × environment interaction. The MMC study will examine oxytocin-related candidate genes (as well as child sex) as potential moderators of child exposure to PDAS and maternal parenting behaviors with respect to a variety of outcomes in middle childhood, including EF and socioemotional adjustment.

In addition to the genetic-based moderation of exposure to PDAS and parenting behavior, there is evidence suggesting that epigenetic alterations in children’s oxytocin-related gene systems may mediate the associations between these exposures and child outcomes. Exposure to maternal depression is associated with differences in glucocorticoid receptor (NR3C1) methylation, which are thought to mediate associations between early-life adversity and subsequent psychopathology [41-45]. Furthermore, a recent genome-wide methylation study of children aged 1 year of mothers with or without major depressive disorder found differences across the epigenome [46]. It is also possible that the association between PDAS, maternal sensitivity, and child changes in epigenetic form may vary by child sex with respect to oxytocin functioning. Although speculative, research suggests that circulating oxytocin levels are comparable across boys and girls during the first 5 years of life [47,48]; however, oxytocin responsiveness to stress is significantly greater for girls who have a history of parental maltreatment as compared with girls without maltreatment or boys in either condition [49]. This alteration in functioning in girls may be a product of epigenetic changes in response to an extreme and maladaptive environment that accelerates reproductive development and individual fitness in the context of a dangerous environment; as such, if female individuals are genetically more sensitive to such environments, it is possible that markers of epigenetic change will be stronger for girls than for boys even at the age of 5 years. The MMC study will examine oxytocin-related changes in the epigenome from early infancy to middle childhood as potential mediators of associations between exposure to PDAS and parenting behavior and developmental outcomes, as well as the potential moderation of this pathway by child sex.

MMC study measures are presented in the following sequence. First, measures pertinent to the RCT examining oxytocin and HPA axis functioning are presented. Next, measures pertinent to testing the hypotheses in aims 1, 2, and 3 are presented in order of appearance in the passive longitudinal design (Figure 1). This begins with maternal oxytocin functioning and oxytocin genetics as primary predictors (maternal oxytocin psychobiology and oxytocin genetics) followed by maternal biopsychosocial outcomes (PDAS and maternal parenting behaviors) and then child biopsychosocial outcomes (child oxytocin epigenetics, stress reactivity, and cognitive and socioemotional outcomes). Summaries of potential moderators and control variables are also provided.

We plan to integrate PDAS measures obtained from repeated assessments of the Edinburgh Postnatal Depression Scale, Beck Depression Inventory–II, and State-Trait Anxiety Inventory by constructing a hierarchical confirmatory factor analysis (HCFA) model with a general PDAS factor giving rise to anxiety and depression subdomains (Figure 2). We will incorporate orthogonal method factors (not shown in the figure) to account for the scale of origin. The HCFA model will be estimated using 1 randomly drawn observation from each mother’s first year of postpartum data. Once we have obtained a good-fitting HCFA model, we will use the parameter estimates from the calibration model to generate time-varying maternal PDAS factor score estimates—1 per month during the first year post partum [84].

Once time-varying PDAS scores have been obtained, we will determine the optimal method for characterizing maternal PDAS. We will consider three approaches: (1) latent mean, (2) growth parameters, and (3) growth mixture classes. We describe each of these approaches in the following paragraphs.

First, we will generate spaghetti plots to examine the shape and heterogeneity of PDAS scores over time. We will use linear mixed models to construct univariate growth models for PDAS to determine empirically whether scores follow a linear, quadratic, or cubic average trajectory and whether there is significant individual heterogeneity in change over time. The following set of equations presents an example of a quadratic growth trajectory for PDAS, where i=1...N indexes individuals and t=0...12 indexes time in months, starting immediately post partum. Time is centered at 6 months post partum in the equation. The b parameters are the fixed effects that characterize the PDAS trajectory for the average participant. The u parameters are individually varying random effects that define each participant’s unique PDAS trajectory. We assume that the u terms follow a multivariate normal distribution centered at 0 and that the residual error is normally distributed at each time point. The residual error is assumed to be uncorrelated with the random effects.

PDAS_ti = b₀ + b₁ × (t − 6) + b₂ × (t − 6)² + u_0i + u_1i × (t − 6) + u_2i × (t − 6)² + r_ti

(1)

If we do not find a significant fixed effect of time (ie, if b₁ and b₂ are not significantly different from 0), we will output empirical Bayes estimates of the random intercept and use these scores to represent each participant’s average PDAS score . We refer to this first approach as the latent mean approach.

If individual trajectories follow the same general form (eg, all linear but with continuously varying slope values), we will output empirical Bayes estimates of mother-specific intercepts (int_i = b₀ + u_0i, representing the average PDAS at 6 mo post partum) and slopes (slope_i = b₁ + u_1i, representing the direction and rate of change over time). In this “growth parameter” approach, we will represent PDAS trajectories using a combination of 2 or more growth parameters.

Finally, if we find visual and empirical evidence of substantial heterogeneity in the shape of PDAS change over time, we will use growth mixture models to estimate a small number of prototypical PDAS trajectory shapes and then use posterior probability scores to assign participants probabilistically to a trajectory class. This third approach will allow for a small number of different trajectories (eg, 1 stable trajectory, 1 increasing linearly, and 1 decreasing to a plateau). Each participant will receive a posterior probability of membership to each trajectory class. Once PDAS trajectories have been operationalized as (1) latent means, (2) growth parameters, or (3) growth mixture classes, we will use the resulting measures to address the study hypotheses.

First, we will use a general linear model to test the hypothesis that participants who are randomized to receive nasal oxytocin experience lower HPA reactivity during the TSST to evaluate this proposed causal pathway (hypothesis 1a). Next, we will also evaluate whether observed oxytocin functioning or genetic variation in oxytocin-associated pathways correlates with PDAS and will then use structural equation modeling (SEM) to assess whether mothers’ HPA functioning mediates the association between oxytocin variation and PDAS (hypotheses 1a and 1b). The SEM framework will enable us to model uncertainty in the growth parameters or latent class membership while simultaneously estimating the 2 mediation pathways (oxytocin to maternal HPA functioning [a path] and maternal HPA functioning to PDAS [b path]). We will use bootstrapped SEs to conduct significance tests for the indirect effects.

Finally, multigroup SEM will be used to examine the potential of moderated mediation, with maternal genetic variation in oxytocin-related pathways moderating the associations between maternal oxytocin and HPA functioning and PDAS (hypothesis 1b). We will also test whether maternal genetic variation moderates the associations between PDAS and maternal HPA functioning, oxytocin functioning, and parenting behaviors (hypothesis 1c). Significant moderation effects will be probed following procedures outlined by Preacher et al [85] to obtain estimated effect conditions on genetic variation within the range of our data.

Our next step of analysis will be to evaluate the extent to which PDAS predicts more adverse child biopsychosocial outcomes during early elementary years (hypothesis 2a), whether high stress reactivity or maladaptive parenting behaviors mediate this association (hypothesis 2b), whether the direct effects of PDAS on child outcomes remain significant after accounting for concurrent maternal depression and anxiety severity (hypothesis 2b), and whether sex moderates the associations between PDAS and child outcomes such that the effects of PDAS on child externalizing are stronger for male individuals and the effects of PDAS on child internalizing are stronger for female individuals. We will continue within an SEM framework for aim 2 analyses as this approach allows for the simultaneous modeling of mediation, moderation, and a variety of response distributions.

In addition to demographic control variables, we will include child temperament (assessed in infancy and concurrently at the age of 5 y) to control for genetic confounds between mothers and children as well as passive gene-environment correlations within the family system.

For aim 3, we will conduct SEM analyses in two ways: (1) including estimated risk scores and (2) stratified by child OXT and OXTR genotypes to determine if child genotype moderates the associations between PDAS and developmental outcomes (hypothesis 3a) and between maternal sensitivity and developmental outcomes (hypothesis 3b). Next, we will test whether the pathways tested in aims 1 and 2 vary by child sex by introducing 2-way interactions between PDAS and child sex and parental sensitivity and child sex and probing significant effects. We will test whether the methylation status of oxytocin-related genetic regions mediates associations among PDAS exposure, sensitivity, attachment quality, and developmental outcomes (hypothesis 3d). Finally, we will conduct exploratory analyses using multigroup models to examine whether the effects of PDAS and maternal sensitivity on epigenetic changes in child oxytocin functioning are stronger for girls than for boys (hypothesis 3e). We will use bootstrapped SEs to conduct significance tests for the indirect effects. Moderation effects will be probed following the procedures outlined by Preacher et al [85].

The proposed study design incorporates extensive repeated assessments. Growth models and SEMs are estimated using full information maximum likelihood, which permits the assumption that missing data are missing at random conditional on observed measures. We will rely on this assumption when data are missing exclusively on endogenous variables. We will use multiple imputation for all analyses involving missing covariates to avoid listwise deletion of cases with missing exogenous variables.

As described previously, attrition and suspension of data collection because of COVID-19 resulted in a lower enrollment in the MMC study than we had hoped; however, we anticipate enrolling 110 dyads from the original sample. With this sample size and a type-I error rate of Cronbach α=.05, we will have 80% power to detect a correlation of r≥0.26, a mean difference of Cohen d≥0.5, and an odds ratio of ≥2.02. Thus, broadly, we expect to have the power to detect moderate, clinically meaningful effects; however, we will have limited power to detect smaller effect sizes that tend to occur in mediation and moderation analyses.