Recent initiatives from the National Institute of Mental Health such as the Research Domain Criteria highlight the importance of identifying dimensions of cognitive and psychosocial functioning that cut across multiple psychiatric categories. These “endophenotypes” are believed to carry pathways of risk for many psychopathological outcomes (Gottesman & Gould, Reference Gottesman and Gould2003). Endophenotypes are measureable psychological or biological phenomena that are intermediate in the causal chain linking genetic liabilities to traditionally defined clinical disorders. These help to explain heterogeneity within disorders, comorbidity across disorders, and the systematization of disease classification (Miller & Rockstroh, Reference Miller and Rockstroh2013). As a result, understanding the genetic contributions to endophenotypes has considerable significance in elucidating the etiology and ultimate classification of psychiatric disorders (Cuthbert, Reference Cuthbert2014).

Executive functioning (EF) is a prototypical endophenotype that has been implicated in several psychiatric and neurodevelopmental conditions, including schizophrenia, bipolar disorder, autism, and attention-deficit/hyperactivity disorder (Luna, Doll, Hegedus, Minshew, & Sweeney, Reference Luna, Doll, Hegedus, Minshew and Sweeney2007; Minzenberg, Laird, Thelen, Cater, & Glahn, Reference Minzenberg, Laird, Thelen, Carter and Glahn2009; Mur, Portella, Martinez-Arán, Pifarré, & Vieta, Reference Mur, Portella, Martínez-Arán, Pifarré and Vieta2007; Willcutt, Doyle, Nigg, Faraone, & Pennington, Reference Willcutt, Doyle, Nigg, Faraone and Pennington2005). EF refers to a clustering of higher order cognitive abilities that support goal-directed action, problem solving, and social functioning. Childhood EF sets the foundation for neuropsychological growth into adulthood (Garon, Bryson, & Smith, Reference Garon, Bryson and Smith2008) and may help to explain continuities in psychiatric difficulties across the life course (Rutter, Kim-Cohen, & Maughan, Reference Rutter, Kim-Cohen and Maughan2006). Thus, identifying genetic liabilities for EF problems in childhood constitutes an important avenue in predicting risk of psychopathology in many domains.

Behavioral genetic studies suggest that EF is heritable (Friedman, Miyake, Robinson, & Hewitt, Reference Friedman, Miyake, Robinson and Hewitt2011; Friedman et al., Reference Friedman, Miyake, Young, DeFries, Corley and Hewitt2008), yet the discrete molecular genetic constituents of EF remain largely unspecified. The preponderance of research points to the monoaminergic genes for dopamine, noradrenaline, and serotonin in the pathogenesis of EF impairment. These genes are believed to exert their effects on EF through prefrontal cortical functioning (Arnsten, Reference Arnsten2011). Nevertheless, the proportion of variance accounted for by individual single nucleotide polymorphisms (SNPs) is small, usually on the order of <1%. This underscores the considerable gap between behavioral and molecular genetic studies known as “missing heritability.” These findings suggest a need to broaden the search for molecular genetic influences on EF.

A separate but related field of inquiry centers on genes of the neurohypophysial hormones, such as oxytocin (OXT). OXT is a known neuromodulator of social behavior and stress regulation, and has been proposed as a key contributor to many of the same psychiatric disorders that are characterized by EF deficits, including autism spectrum disorders, affective disorders, and schizophrenia (Feldman, Monakhov, Pratt, & Ebstein, Reference Feldman, Monakhov, Pratt and Ebstein2016). These effects are frequently attributable to disruptions in social functioning, consistent with the role of OXT in face recognition, altruism, empathy, and the ability to infer others’ mental states (Bartz, Zaki, Bolger, & Ochsner, Reference Bartz, Zaki, Bolger and Ochsner2011; Meyer-Lindenberg, Domes, Kirsch, & Heinichs, Reference Meyer-Lindenberg, Domes, Kirsch and Heinrichs2011). However, these phenomenologically complex behaviors may partly depend on executive functions such as attentional control, set shifting, working memory, and response inhibition. For instance, it has been suggested that empathy and mentalization draw on the ability to inhibit and shift perspectives (Benson, Sabbagh, Carlson, & Zelazo, Reference Benson, Sabbagh, Carlson and Zelazo2013; Vetter, Altgassen, Phillips, Mahy, & Kliegel, Reference Vetter, Altgassen, Phillips, Mahy and Kliegel2013), while altruistic sharing has been shown to engage inhibitory control strategies (Aguilar-Pardo, Martinez-Arias, & Colmenares, Reference Aguilar-Pardo, Martinez-Arias and Colmenares2013). OXT has recently been linked to “nonsocial” aspects of cognition, including spatial and episodic memory, short-term verbal memory, cognitive flexibility, and the executive component of working memory in schizophrenia (Chini, Leonzino, Braida, & Sala, Reference Chini, Leonzino, Braida and Sala2014; Michalopoulou et al., Reference Michalopoulou, Averbeck, Kalpakidou, Evans, Bobin, Kapur and Shergill2015). These effects are further supported by recent evidence suggesting that, within the postpartum period, blockage of the OXT receptor (OXTR) inhibits cognitive flexibility among mothers (Albin-Brooks, Nealer, Sabihi, Haim, & Leuner, Reference Albin-Brooks, Nealer, Sabihi, Haim and Leuner2017). Thus, there is growing evidence that traditionally “social” nonapeptides such as OXT may be related to discrete executive abilities that are impaired in psychiatric conditions and required for adaptive social behavior.

Most genetic research on the behavioral correlates of OXT has focused on variability in the OXTR gene. This gene (located on chromosome 3p25.3) has many SNPs in both the intronic and exonic regions. The most extensively studied SNPs, and those with the most support in both healthy and clinical samples, are rs53576 and rs2254298 (see Feldman et al., Reference Feldman, Monakhov, Pratt and Ebstein2016; Kumsta & Heinrichs, Reference Kumsta and Heinrichs2013). The current study focused on the latter SNP (rs2254298) for two reasons. First, there is strong meta-analytic evidence that this SNP is associated with neurodevelopmental conditions typified by EF impairment, especially after accounting for study bias (LoParo & Waldman, Reference LoParo and Waldman2015); and second, neuroimaging studies specifically point to rs2254298 as relating to structural and functional brain differences in regions implicated in EF. In particular, rs2254298 has been associated with regional volumetric differences in the amygdala and dorsal anterior cingulate cortex (dACC) in both adults and children (Furman, Chen, & Gotlib, Reference Furman, Chen and Gotlib2011; Inoue et al., Reference Inoue, Yamasue, Tochigi, Abe, Liu, Kawamura and Rogers2010). The functionality of rs2254298, including the coupling of dACC to limbic structures, has been replicated in both Caucasian and Japanese samples (Tost et al., Reference Tost, Kolachana, Verchinski, Bilek, Goldman, Mattay and Meyer-Lindenberg2011; Yamasue et al., Reference Yamasue, Suga, Yahata, Inoue, Tochigi, Abe and Takei2011). Wang et al. (Reference Wang, Braskie, Hafzalla, Faskowitz, McMahon, de Zubicaray and Thompson2017) recently demonstrated that, among 10 genotyped OXTR SNPs, only rs2254298 was associated with resting-state functional connectivity between the posterior cingulate cortex and the dACC. Moreover, A-allele carriers demonstrated reduced cortical thickness in the bilateral dACC compared to individuals with a GG genotype. These findings are interesting given that nearly all models of neurocognitive functioning posit a central role of the dACC in the modulation of EF (Bush, Luu, & Posner, Reference Bush, Luu and Posner2000; Niendam et al., Reference Niendam, Laird, Ray, Dean, Glahn and Carter2012). Finally, a recent study found that rs2254298 was associated with parenting behavior through maternal decision making, supporting the notion that the effect of OXTR on social behavior may operate through discrete executive functions in adults (Cost et al., Reference Cost, Unternaehrer, Plamondon, Steiner, Meaney, Atkinson and Fleming2016). Despite these clear connections, and that many studies have examined rs2254298 in clinical conditions with documented EF impairment, no study has examined whether genetic variability in OXTR is associated with EF in children. Thus, the first goal of the present study was to examine whether genetic variability in rs2254298 is associated with individual differences in childhood EF. A sample of preschool-aged children was selected in light of the importance of early-emerging cognitive skills in the developmental onset of neuropsychiatric disorders and social behavior.

The second goal of this study was to determine whether OXTR interacted with birth weight variability on EF. Birth weight is a proxy for intrauterine growth related to the quality of the prenatal environment. Variance in neurocognitive functioning and psychopathology in adulthood involves the programming of biobehavioral systems that can often be traced to fetal development (Schlotz & Phillips, Reference Schlotz and Phillips2009). Studies suggest that normative birth weight variation is associated with individual differences in EF in the preschool period (Wade, Browne, Madigan, Plamondon, & Jenkins, Reference Wade, Browne, Madigan, Plamondon and Jenkins2014). Perhaps more interesting, it has been shown that birth weight is linearly associated with increased regional surface area and volume in the anterior cingulate cortex of children and adolescents, which in turn is associated with improved cognitive control (Walhovd et al., Reference Walhovd, Fjell, Brown, Kuperman, Chung, Hagler and Akshoomoff2013). Thus, variation in prenatal neurodevelopment appears to be linked to higher order executive functions across childhood (Raznahan, Greenstein, Lee, Clasen, & Giedd, Reference Raznahan, Greenstein, Lee, Clasen and Giedd2012). It has also been suggested that cortical differences associated with birth weight may interact with genetic risks for psychiatric disorders (Haukvik et al., Reference Haukvik, Rimol, Roddey, Hartberg, Lange, Vaskinn and Agartz2013).

A gene–environment interaction between OXTR and birth weight is indicated by several strands of evidence. First, evolutionary models suggest that the rs2254298 A (minor) allele may be a susceptibility allele that, if combined with adversity during childhood or the prenatal period, increases the risk of developing psychopathologies (Brüne, Reference Brüne2012). Second, studies on the epigenetic modulation of OXTR have demonstrated that the degree of DNA methylation is associated with neural activity in the dACC during a social attention task (Jack, Connelly, & Morris, Reference Jack, Connelly and Morris2012). This provides molecular evidence that OXTR modification may be associated with brain functioning in regions that support EF. In addition, recent studies have documented global DNA methylation patterns associated with perinatal outcomes such as birth weight (Michels, Harris, & Barault, Reference Michels, Harris and Barault2011). Consistent with both the evolutionary and molecular accounts, it has recently been shown that prenatal adversity is specifically linked to DNA methylation of OXTR in cord blood at birth (Unternaehrer et al., Reference Unternaehrer, Bolten, Nast, Staehli, Meyer, Dempster and Meinlschmidt2016). Given that birth weight is a marker of fetal development that is influenced by myriad maternal adversities, these results raise the possibility that intrauterine conditions related to fetal growth may contribute to an epigenetic mechanism of fetal programming that impacts later cognitive and mental health outcomes of children. Thus, in the current study we tested the hypothesis that birth weight variability interacts with OXTR in the prediction of children's EF. We hypothesized that children would show more problems with EF as a function of lower birth weight in combination with more copies of the previously demonstrated risk allele (A) of OXTR marker rs2254298.

Methods

Study sample

Multiparous women giving birth to infants in the cities of Toronto and Hamilton, Ontario, Canada, between 2006 and 2008, who had been contacted by the Healthy Babies Healthy Children public health program (run by Toronto and Hamilton Public Health Units), were considered for participation. Inclusion criteria for the intensive sample of the Kids, Families, Places study were as follows: (a) an English-speaking mother, (b) a newborn weighing >1500 g, (c) one or more children <4 years old in the home, and (d) agreement to the collection of observational and biological data. Thirty-four percent of mothers whose information was passed by the Healthy Babies Healthy Children program consented to participate in the study. Reasons for nonenlistment included inability to contact families, ineligibility once contacted, and refusals. The current sample was nested within a larger longitudinal project, the goals of which were to examine contextual and biological influences on children's social–emotional and cognitive development using within-family methodology, thus requiring a minimum of two children per family. However, as we were interested in EF in the preschool period, only target newborns were included. The University of Toronto Research Ethics Board approved all procedures, including informed consent.

Data for the current study were primarily drawn from Wave 4 of the Kids, Families, Places longitudinal study, at which point children were about 4.5 years old (N = 323; M _age = 4.79 years, SD = 0.28). Attrition up to Wave 4 related to younger maternal age at first pregnancy, t (494) = –5.10, p < .001, lower socioeconomic status, t (498) = –5.07, p < .001, and lower maternal education, t (498) = –2.99, p < .005. Of the 323 children participating at Wave 4, outcome data on EF were unavailable for 12 of them due to child noncompliance or tester administration error. One additional child was excluded because he was 5.87 SD above the birth weight mean. The final sample therefore consisted of 310 children (51% female). Demographics of the study sample are provided in Table 1. The sample consisted of typically developing children covering a wide range of birth weight (1900–5000 g), consistent with recent findings that variability in this range is associated with individual differences in neurocognitive development.

Table 1.

Descriptive statistics for study variables

^a This is a standardized score with a mean of zero.

^b Reported on a scale from 1 (no income) to 16 ($105,000 or more).

Procedure

A home visit of approximately 2 hr was conducted in which parents (usually the mother) filled out questionnaires pertaining to family life, sociodemographic information, and child health and behavior. Children also underwent a battery of observational and standardized tasks assessing multiple domains of functioning, including EF (described below).

Measures

EF

This was assessed using two previously developed and widely used tasks that were appropriate to the age of the child: Bear/Dragon (Reed, Pien, & Rothbart, Reference Reed, Pien and Rothbart1984) and the Dimensional Change Card Sort (DCCS; Zelazo, Reference Zelazo2006). For the Bear/Dragon task, children were instructed to do what they were told by the nice bear (e.g., “touch your nose”), but not to do what they were told by the mean dragon. Children were scored on their total number of correct responses (0–10) on five dragon and five bear trials. For the DCCS, children were required to sort a series of bivalent test cards, first according to one dimension (e.g., color), and then according to the other (e.g., shape). Children who passed the postswitch phase of the standard version of the DCCS proceeded immediately to the border version, which used the same target cards as the standard version. The border version consisted of 12 trials. Children were required to sort cards based on “border” criteria (“If there's a border, play the color game. If there's no border, play the shape game”). Previous studies have shown that the Bear/Dragon and DCCS load onto the same latent factor measuring set shifting/cognitive flexibility, working memory, and inhibitory control (Bernier, Carlson, Deschênes, & Matte-Gagné, Reference Bernier, Carlson, Deschênes and Matte-Gagné2012). In the present study, Bear/Dragon and DCCS were significantly correlated, r (299) = .34 p < .001. Thus, the two tasks were z scored and averaged to create a composite, with higher scores reflecting better EF ability.

Birth weight

At Wave 1 (child was a newborn), mothers reported on their child's birth weight in kilograms or grams. Following the removal of one baby who was an outlier, described above, this variable was normally distributed: skewness = 0.022 (SE = 0.14), kurtosis = 0.60 (SE = 0.28). The range on the variable was 1.90–4.99 kg (M = 3.42, SD = 0.48). We residualized birth weight for gestational age to give a purer metric of fetal growth (Phua, Rifkin-Graboi, Saw, Meaney, & Qiu, Reference Phua, Rifkin-Graboi, Saw, Meaney and Qiu2012).Footnote ¹

Covariates

Several covariates were controlled in the adjusted analysis including family income/assets, which was measured as a composite that included annual family income, assessed on a scale from 1 (no income) to 16 ($105,000 or more); number of rooms in the family's residence; and whether the family owns/co-owns their house/apartment and/or a vehicle (yes/no variables). These four variables were standardized, and averaged to create a composite for income/assets. Higher values represented higher SES; maternal education, assessed as the total number of years of formal schooling, not including kindergarten; child age (in years and months); and gender (0 = male, 1 = female).

Genotyping

DNA was extracted from cheek swab samples from both children and their biological parents. The rs2254298 SNP was genotyped with the ABI 7900-HT Sequence Detection System^® (Applied Biosystems) using the TaqMan 5′ nuclease assay for allelic discrimination. The polymerase chain reactions (10 μl volume) contained 30 ng of genomic DNA, 10 μmol of TaqMan^® Universal PCR Master Mix (Applied Biosystems), and 0.25 μl of allelic discrimination mix (Applied Biosystems) containing 36 μM of each primer and 8 μM of each probe. The thermal cycling conditions were 50 °C for 2 min, 95 °C for 10 min, followed by 40 cycles of 94 °C for 15 s and the annealing temperature of 60 °C for 1 min. Each 96 well plate contained two negative controls. Plates were then read on the ABI 7900HT Sequence Detection System using the allelic discrimination end-point analysis mode of Sequence Detection System software package version 2.0 (Applied Biosystems).

Regarding SNP characteristics of rs2254298, the major allele (G) had a frequency of .88. This marker is located in intron 3 of the OXTR gene, with a chromosome position of 8,802,227. This SNP conformed to Hardy–Weinberg equilibrium (p = .12) in this sample.

Statistical analysis

The analysis was carried out using the Pedigree-Based Association Test (PBAT) with Golden Helix SVS. PBAT is robust to population admixture and stratification, meaning population subgroups and ethnic differences in allelic distributions are controlled. Specifically, the alleles transmitted to probands are compared to the distribution of alleles expected among offspring, conditional on parental genotypes and based on Mendelian transmission. In other words, parents are used as pseudocontrols. The result is an unbiased test statistic that allows inferences for both association and linkage. We tested an additive model in which genotypes were coded to reflect the number of copies of each allele. The quantitative trait, EF, was offset by the phenotypic mean to generate a phenotypic residual for each person (i.e., mean centered). This increases the statistical power of the PBAT statistic (Lange, DeMeo, & Laird, Reference Lange, DeMeo and Laird2002). To further improve power, we controlled for nongenomic factors (covariates) related to attrition and known to be associated with EF (Laird & Lange, Reference Laird and Lange2008).

Computation of the main genetic and interaction effect followed the procedure outlined by Vansteelandt et al. (Reference Vansteelandt, DeMeo, Lasky-Su, Smoller, Murphy, McQueen and Silverman2008). The approach uses a semiparametric model: there is a nonparametric component for the main environmental effect that is robust to misspecification (Moerkerke, Vansteelandt, & Lange, Reference Moerkerke, Vansteelandt and Lange2010), and a parametric component for the main genetic and the gene–environment interaction effect. We fit the model:

$$E\left( {{\rm Y}_{ij} \vert X_{ij}\comma \, Z_{ij}\comma \, S_i} \right){\rm} = {\rm \mu} (Z_{ij}\comma \, S_i ) + {\rm \beta} _{{\rm SNP}} X_{ij} + {\rm \beta} _{{\rm inter}} X_{ij} Z_{ij}\comma $$

where i indexes the family, j the offspring, Y is the trait (i.e., EF), X is the additive coding of the genotype (0, 1, or 2 copies), Z is the environmental exposure (i.e., birth weight), S encodes the dependency on the parental genotypes, and μ(.) is a nonparametric function that encodes the dependence of the trait on the environmental exposure and parental genotype. A positive value for the interaction indicates that there is an increase in EF with more copies of that allele combined with higher levels of the exposure variable (birth weight).

Below, we report results using the family-based association test (FBAT) and FBAT-I p value, heritability coefficient (h ²), and standard effect size estimates where appropriate. The FBAT p value is the primary test statistic for examining the main genetic effect in family-based designs. It is a generalization of the transmission disequilibrium test, appropriate for quantitative traits (Lange et al., Reference Lange, DeMeo and Laird2002). The FBAT-I p value tests the null hypothesis of no gene–environment interaction. The heritability coefficient (h ²) reflects the proportion of phenotypic variance explained by the analyzed SNP or interaction, and is thus a measure of effect size.

Results

Table 2 reports bivariate correlations among study variables. For the outcome phenotype, EF, significant associations with age, maternal education, and family income/assets were observed, while there were no gender differences on EF. There were no gender differences in the distribution of the G allele, χ² (df = 1) = 0.022, p = .88, or the A allele, χ² (df = 1) = 0.072, p = .79. There were also no gender differences in overall genotype, χ² (df = 2) = 0.080, p = .96. All covariates were retained in the adjusted models outlined below.

Table 2.

Bivariate correlations between study variables

†p < .10. *p < .05. **p < .01. ***p < .001.

Computation of the interaction effect using PBAT assumes genotype–environment independence, which is generally a reasonable assumption in family-based tests (Umbach & Weinberg, Reference Umbach and Weinberg2000). To confirm this assumption, we performed a preliminary association analysis between the selected OXTR SNP, rs2254298, and the environmental exposure, birth weight. There was no significant relation between rs2254298 and birth weight in either the unadjusted model (FBAT p = .54, h ² = 0.003) or the adjusted model that controlled for covariates (p = .53, h ² < 0.001), supporting model assumptions.

Regarding the main genetic effect, there was a significant association between the rs2254298 marker and EF in the unadjusted analysis (FBAT p = .049) and in the adjusted analysis controlling for covariates (FBAT p = .023). With regard to direction of effects, more copies of the major allele (G) were associated with higher EF, while more copies of the minor allele (A) were associated with lower EF. We then tested the FBAT association test of the interaction effect (FBAT-I), which revealed a significant interaction between rs2254298 and birth weight in both the unadjusted (FBAT-I p = .011) and adjusted analyses (FBAT-I p = .003).

Next, we examined the main genetic and interaction effects separately, each controlling for the other. In the unadjusted model, the main genetic effect (known as QBAT) was not significant for the G allele, β (SE) = 0.20 (0.13), p = .10, h ² = 0.014, or the A allele, β (SE) = –0.16 (0.14), p = .27, h ² = 0.010. Note that the p-values are not identical for the two alleles because of the semiparametric model used. Effects were likewise nonsignificant in the adjusted model for the G allele, β (SE) = 0.18 (0.13), p = .16, h ² = 0.014, and the A allele, β (SE) = –0.13 (0.14), p = .36, h ² = 0.010. In contrast, strong gene–environment interaction effects were observed between rs2254298 and birth weight. In the unadjusted model, the interaction (known as QBAT-E) was significant for the G allele, β (SE) = 0.26 (0.09), p = .0025, h ² = 0.27, and for the A allele, β (SE) = –0.46 (0.17), p = .007, h ² = 0.052. Results were robust in the adjusted model, with a significant interaction for the G allele, β (SE) = 0.31 (0.11), p = .0045, h ² = 0.52, and for the A allele, β (SE) = –0.62 (0.22), p = .0051, h ² = 0.12. The direction of effects was consistent with study hypotheses, namely, that EF increased as a function of more copies of the G (major) allele in combination with higher birth weight, and decreased with more copies of the A (minor) allele and lower birth weight. There were no significant Gene × Gender × Birth Weight interactions.