Birth weight is associated with adolescent brain development: A multimodal imaging study in monozygotic twins.

Previous research has shown that the prenatal environment, commonly indexed by birth weight (BW), is a predictor of morphological brain development. We previously showed in monozygotic (MZ) twins associations between BW and brain morphology that were independent of genetics. In the present study, we employed a longitudinal MZ twin design to investigate whether variations in prenatal environment (as indexed by discordance in BW) are associated with resting-state functional connectivity (rs-FC) and with structural connectivity. We focused on the limbic and default mode networks (DMNs), which are key regions for emotion regulation and internally generated thoughts, respectively. One hundred and six healthy adolescent MZ twins (53 pairs; 42% male pairs) followed longitudinally from birth underwent a magnetic resonance imaging session at age 15. Graph theoretical analysis was applied to rs-FC measures. TrackVis was used to determine track count as an indicator of structural connectivity strength. Lower BW twins had less efficient limbic network connectivity as compared to their higher BW co-twin, driven by differences in the efficiency of the right hippocampus and right amygdala. Lower BW male twins had fewer tracks connecting the right hippocampus and right amygdala as compared to their higher BW male co-twin. There were no associations between BW and the DMN. These findings highlight the possible role of unique prenatal environmental influences in the later development of efficient spontaneous limbic network connections within healthy individuals, irrespective of DNA sequence or shared environment.

birth weight

monozygotic

resting‐state functional connectivity

INTRODUCTION

Studies have indicated that the in utero environment is a predictor of mental and physical health later in life (e.g., Newman et al., 2016). Associations have been found, for example, between in utero adversity (e.g., malnutrition, maternal smoking) and risk for attention deficit hyperactivity disorder, anxiety disorders, and conduct disorders (Mathewson et al., 2017; Newman et al., 2016; Van den Bergh et al., 2017). While numerous investigations have emerged indexing the potential link between in utero adversity and brain structure later in life (Adamson, Letourneau, & Lebel, 2018; Miguel, Pereira, Silveira, & Meaney, 2019), less is known about the relationship between in utero environment and intrinsic functional connectivity, which allows for the study of the brain's natural rhythms/connectivity, unrestrained by the administration of a particular task (He, Snyder, Zempel, Smyth, & Raichle, 2008). As there is substantial evidence showing a link between altered resting‐state functional connectivity (rs‐FC) and psychopathology (Broyd et al., 2009; Xia et al., 2018), learning the specific mechanisms of these associations would help to identify children at risk and help develop early preventive interventions.

Of the various general indicators of the prenatal environment, birth weight (BW) is among the best studied, correlating with a variety of in utero conditions (e.g., nutrition, stress), and a predictor of developmental, behavioral, and cognitive outcomes in children, adolescents, and adults (Ment et al., 2009; Shenkin, Starr, & Deary, 2004). For example, studies have shown that children with very low BW (i.e.,<1,500 g) or extremely low BW (i.e., <1,000 g) are more likely to have difficulties in regulating emotions (Clark, Woodward, Horwood, & Moor, 2008), and are more likely to show poorer social adjustment (Ritchie, Bora, & Woodward, 2015), as compared to age‐matched children of normal BW. Furthermore, various studies have now shown that, analogous to the observed associations in very low or extremely low BW children, associations between BW and developmental outcome are also evident among individuals in the normal BW range. For example, adolescents who had lower (yet still within the normal range of) BW showed impaired inhibitory control as compared to adolescents with normal BW, suggesting difficulty regulating one's actions (Schlotz, Godfrey, & Phillips, 2014). Notably, rs‐FC in the default mode network (DMN) and in the limbic network have been linked to the ability to regulate oneself, along with one's emotions, attention, and social interactions (Broyd et al., 2009; Pannekoek et al., 2014; Xia et al., 2018). Other research has shown that rs‐FC measures are relatively stable in children and adolescents, and possible markers for neural network development (Thomason et al., 2011). Together, these findings highlight the potential relevance of DMN and limbic brain networks in investigations of the link between BW and behavioral and cognitive outcomes across development.

In addition to behavioral and cognitive outcomes, several studies conducted in children, adolescents, and adults have found that BW is linked to brain development. For example, individuals with very low BW have been found to have smaller brains overall and alterations in cortical morphology, as compared to individuals with normal BW (Martinussen et al., 2005; Nagy et al., 2009; Shang et al., 2019). Notably, even when BW is within the normal range, structural associations with BW have been found throughout the brain, particularly in midline structures such as the anterior cingulate cortex, paracentral lobule, and precuneus, along with the orbitofrontal, inferior parietal, and temporal cortices (Haukvik et al., 2014; Shang et al., 2019; Walhovd et al., 2012). Additionally, a recent study found that adolescents born very preterm showed altered rs‐FC between the amygdala and various regions (e.g., hippocampus, prefrontal cortex), relative to adolescents born at term with normal weight (Johns, Lacadie, Ment, & Scheinost, 2019). Overall, these studies support the relevance of in utero environment for long‐term typical and atypical structural and functional brain development.

As brain structure and function are driven by both genetic and environmental factors, a methodological challenge when studying associations with prenatal environment is to tease apart the relative contributions of genetics and environment (Chiarella, Tremblay, Szyf, Provencal, & Booij, 2015). By using a monozygotic (MZ; i.e., identical) twin design, one can dissociate the effect of the individual environment from shared environmental and genetic contributions (Chiarella et al., 2015; Craig, Calais‐Ferreira, Umstad, & Buchwald, 2020; Vitaro, Brendgen, & Arseneault, 2009).

MZ twins are the ideal genetic controls for studying environmental contributions to phenotypes, because MZ twin pairs share essentially the same DNA sequence (Chiarella et al., 2015; Craig et al., 2020). In MZ twins, environmental differences in utero within a twin pair could lead to discordance in BW and differential later‐life health conditions, even if they share the same genetic setup. Although specific mechanisms are not known, differences in BW may occur due to differences in vascular architecture or to differences in placental mass (Hack et al., 2008; Nikkels, Hack, & van Gemert, 2008). Interestingly, studies of MZ twin pairs have shown that the lower BW twin has lower intelligence quotient scores, relative to the higher BW twin (Edmonds et al., 2010; Newcombe, Milne, Caspi, Poulton, & Moffitt, 2007). Other studies have found that greater BW discordance was associated with greater within‐pair differences in child problematic behavior (Lim et al., 2018; Tore et al., 2018; van Os et al., 2001). We and others have previously found in MZ twin pairs that overall brain volume and cortical surface area were positively associated with BW discordance in frontal and temporal brain regions in adolescents and adults (Casey et al., 2017; Levesque et al., 2015; Raznahan, Greenstein, Lee, Clasen, & Giedd, 2012). Although within‐pair differences are relatively small, the preceding results highlight the role of the unique prenatal environment in later brain morphology development. Given the importance of the limbic network and the DMN for emotional and cognitive functioning, along with the previously reported associations in singleton designs between BW and emotional, cognitive and brain development, investigating the association between BW discordance and rs‐FC in the limbic network and DMN in a MZ twin design could provide important insights into the mechanisms of how unique environmental factors in utero affect long‐term outcomes.

In the present study, we investigated in MZ adolescent twins recruited from a longitudinal cohort and followed since birth whether differences in the prenatal environment were associated with within‐pair differences in rs‐FC within the limbic network and DMN. Graph theory has been used to characterize brain networks alterations in different mental health disorders (Sato et al., 2018), suggesting that graph theory metrics are sensitive to BW variation. Graph theory measures are based on functional connectivity and allow for the characterization of functional brain networks using graphs with neural elements represented by nodes and functional connections between nodes (i.e., synapses, axonal projections) represented by edges or links (Bullmore & Sporns, 2009). Graph theory uses small‐world properties of networks to extract key features of networks at local and global scales (Watts & Strogatz, 1998). This specifically allowed us to study rs‐FC at the network level and link it to variation in BW, in two well‐characterized neural networks that are related to emotional and cognitive regulation (Sheline et al., 2009; Spittle et al., 2009): the limbic network and the DMN. Following this, we subsequently examined differences in underlying structural network connectivity as a function of BW variation. We took this approach based on prior work suggesting that while the structure and function of the “resting brain” often converges, diverging accounts can indicate certain functional connections may be mediated by indirect, third‐region, structural connections (Damoiseaux & Greicius, 2009). Thus, supplementing the resting‐state findings with underlying structural connectivity information could highlight any potential functional connections that are robust only due to additional connections. Finally, we also conducted exploratory analyses investigating potential sex differences based on previous research in singletons showing sex‐specific associations between in utero adversity and brain structure (Kersbergen et al., 2016; Reiss et al., 2004) and on research showing sex differences in rs‐FC in typically developing adolescents (Alarcón, Cservenka, Rudolph, Fair, & Nagel, 2015; Kilpatrick, Zald, Pardo, & Cahill, 2006).

METHODS

Participants

One hundred and six 15‐year‐old MZ twins (53 pairs; 31 female pairs), mean age ± SD: 15.7 ± 0.3 years, range 15.3–16.7 years; mean BW ± SD: 2,550 ± 54 g, range 1,000–3,730 g; mean discordance ± SD: 300 ± 250 g, range 10–1,070 g) were recruited from the Québec Newborn Twin Study (QNTS; Boivin et al., 2019). The majority of the twin pairs were Caucasian (94.3%), and average maternal education ± SD at the time of birth of the twins was 12.1 ± 2.6 years. The QNTS recruited participants born in Québec, Canada between April 1, 1995 and December 31, 1998 through the Québec Ministry of Health and Social Services registry, and have followed them longitudinally since recruitment (Boivin et al., 2019). Zygocity was determined by analyzing eight or nine polymorphic genetic markers, and monozygocity was confirmed for those twins who were concordant for every genetic marker (Boivin et al., 2019). Good current health was reported for all twins who underwent brain imaging, with no history of physical or psychiatric illness and no use of psychotropic medications. All twins underwent a behavioral assessment to rule out lifetime disorders, using the Kiddie‐Schedule for Affective Disorders and Schizophrenia (Kaufman et al., 1997). The parents and twins provided written informed consent and assent, respectively. The study protocol was approved by the appropriate ethics committees (i.e., Sainte‐Justine Hospital Research Centre and the Montréal Neurological Institute [MNI]) and was in compliance with the Code of Ethics of the World Medical Association (Declaration of Helsinki).

Neuroimaging procedure

Acquisition

Participants were scanned on a 3 T TIM Trio Siemens MRI Scanner with a 32‐channel head coil. The scan included a magnetization‐prepared rapid acquisition gradient‐echo (MPRAGE) 9‐min anatomical sequence (176 sagittal slices; 1 mm thickness, TR = 2,300 ms, TE = 2.98 ms, TI = 900 ms, flip‐angle = 9°, field of view [FOV] = 240 × 256 mm, matrix size 256 × 256). Resting‐state functional data were acquired with BOLD Mosaic resting echo‐planar imaging (EPI; 165 volumes, TR = 2,000 ms, TE = 30 ms, FA = 90°, 6‐min 8 s' acquisition time, slice thickness = 3.5 mm, descending acquisition, FOV = 224 mm, matrix size 64 × 64). Diffusion EPI was acquired with a b‐value of 1,000 s/mm2 in 64 directions. (63 axial slices; 2 mm thickness, TR = 8,500 ms, TE = 89 ms, FOV = 256 × 256 mm, matrix size 128 × 128).

The twins also performed an implicit emotional processing task, where they were asked to identify the sex of different faces that varied on emotionality (15 min; not discussed here). Portions of the structural and task‐based functional connectivity data have been published (Casey et al., 2017; Ismaylova, Levesque, et al., 2018; Levesque et al., 2015).

Processing

rs‐FC analyses were performed using MATLAB R2016A, SPM 12, and the CONN Functional Connectivity toolbox v17 (Whitfield‐Gabrieli & Nieto‐Castanon, 2012; http://www.nitrc.org/projects/conn). Preprocessing was done in CONN, specifically volume realignment, slice‐timing correction, coregistration of functional and structural scans followed by normalization to the MNI standard space, segmentation into three tissue compartments (gray matter, white matter, and cerebrospinal fluid [CSF]), outlier identification, and smoothing with an 8 × 8 × 8 mm FWHM Gaussian filter.

The default CONN denoising steps were applied, including the CompCor strategy which first extracts the five main components from the white matter and CSF signal and regresses them out from the signal. Then, it performs scrubbing to remove the signal from invalid scans due to too much movement (scan‐to‐scan subject motion threshold >1 or scan‐to‐scan subject rotation threshold >0.02) or due to outlier signals, and were used as additional regressors (noise), and therefore excluded from the rest of the signal. Motion regression (12 regressors), and band‐pass filtering between 0.008 and 0.09 Hz (Behzadi, Restom, Liau, & Liu, 2007) were also applied at the denoising step. In our sample, an average of 4.6 volumes were scrubbed (range 0–36, i.e., 0–20%), with no difference between low and high BW twins. No participant was excluded based on denoising. Mean BOLD time series were extracted from each region of interest (ROI) and correlated with the BOLD time‐series signal of every other ROI in the network to create an ROI‐to‐ROI correlation map showing connectivity between each region within the network.

Graph theoretical analyses were applied to quantify rs‐FC patterns and brain networks. The two measures we opted to investigate were global efficiency and average path length, as they provide an index of overall network efficiency (Achard & Bullmore, 2007; Bullmore & Sporns, 2009). Greater efficiency is indexed as higher global efficiency values and lower average path lengths (Bullmore & Sporns, 2009), meaning that information flow is maximized between regions within a network. Graph theory values were extracted by selecting the specific brain regions corresponding to two networks: the limbic network and the DMN. All brain regions (henceforth referred to as nodes when discussing graph theoretical analyses) were chosen based on the FSL version of the Harvard‐Oxford atlas as implemented in CONN (Desikan et al., 2006). rs‐FC was computed as the average time series in each regional mask (see Figure 1). The DMN included nine nodes: the precuneus, posterior cingulate gyrus, medial frontal cortex, inferior parietal lobule (both the left and right anterior and posterior supramarginal gyrus), and both the left and right angular gyrus. The limbic network included nine nodes: the left and right insula, anterior cingulate cortex, left and right amygdala, left and right hippocampus, and left and right orbitofrontal cortex (OFC). In addition to being core brain regions pertaining to the limbic network and to the DMN, in our previous studies, we noted variations in rs‐FC in those networks in relation to mood in healthy adults (Ismaylova, di Sante, et al., 2018 ). Furthermore, in previous work in the same adolescent cohort, we observed an association between BW and morphological development (Casey et al., 2017; Levesque et al., 2015). To only include highly efficient networks (i.e., well‐connected nodes in a small‐world regime), connectivity matrices were thresholded at a wiring cost >0.15 at the individual level (Achard & Bullmore, 2007). We exported the pertinent values from CONN to use for statistical analysis.

FIGURE 1

The nine nodes for the default mode network (DMN) (a) and limbic network (b). aSMG l, left anterior supramarginal gyrus; pSMG l, left posterior SMG; AG l, left angular gyrus; PC, posterior cingulate gyrus; MedFC, medial prefrontal cortex; ASMG r, right anterior SMG; pSMG r, right anterior SMG; AG r, right angular gyrus; AC, anterior cingulate; IC l, left insular cortex; hippocampus l, left hippocampus; amygdala l, left amygdala; F orb l, left orbitofrontal cortex; IC r, right insular cortex; hippocampus r, right hippocampus; amygdala r, right amygdala; F orb r, right orbitofrontal cortex

Structural network connectivity was examined by extracting diffusion information for those network(s) for which rs‐FC was/were found to be significantly associated with BW variation. Track count indexes structural connectivity in the brain (Gao et al., 2013; Hagmann et al., 2007). FSL tools (http://www.fmrib.ox.ac.uk/fsl; Jenkinson, Beckmann, Behrens, Woolrich, & Smith, 2012) were used for motion and eddy current correction of the Diffusion Weighted Images (DWI), as well as image coregistration. The Harvard Oxford label set was coregistered to each DWI using the same ROIs listed above. Diffusion Toolkit and TrackVis (Wang, Benner, Sorensen, & Wedeen, 2007) were used to generate fiber tracks. Fractional anisotropy maps were first reconstructed (dti_recon) and then used for fiber tracking (dti_tracker) using the FACT algorithm, and a 35° curvature limit. The resulting tracks were then spline filtered (spline_filter). TrackVis was then used to generate maps of the tracks passing through each ROI and between each pair of ROIs in the DMN and limbic networks.

Statistics

Global efficiency and path length were analyzed using linear mixed models, with sex (male/female) and BW (higher/lower) as fixed factors and family as the random factor. Only the main effects of sex and BW were added to the model (not the interaction), and the intercept was included. For the random effect, the covariance type was set to scaled identity, and the intercept was included (as per Casey et al., 2017). We studied BW difference using two analytical methods. First, we split each twin pair by BW, assigning the higher BW twin to one group, and the lower BW twin to another (Casey et al., 2017), and then compared the rs‐FC of the two groups. Second, for those significant network(s), we followed up by calculating a within‐twin‐pair BW difference score (i.e., higher BW–lower BW) to provide a continuous, within‐pair rs‐FC discordance for global efficiency and average path length. Those values were analyzed using the R package “mztwinreg” (Córdova‐Palomera, 2015; https://CRAN.R‐project.org/web/packages/mztwinreg/mztwinreg.pdf). Specifically, the env_diff function, based on the work of Carlin, Gurrin, Sterne, Morley, and Dwyer (2005), was used to fit linear regression models without intercept that allow for testing whether intrapair differences in a predictor variable are related to intrapair differences in an outcome variable. The DTI track counts were analyzed using the same linear mixed effects model employed in the rs‐FC analyses, through MATLAB using the SurfStat toolbox (http://www.math.mcgill.ca/keith/surfstat/).

RESULTS

Default mode network

Neither global efficiency (Figure 2a top; t(103) = 1.03, p = .31) nor average path length (Figure 2b bottom; t(103) < 1, p = .40) showed any relationship with BW. There were no main effects of sex (ps > .7).

FIGURE 2

Violin plots depicting global efficiency (top) and average path length (bottom) for both the (a) default mode network (DMN) and (b) limbic network (limbic), as a function of birth weight (BW). †p = .05, black dots indicate the mean, and black vertical lines the error bars (SD). Values in the plot represent the means

Limbic network

Higher BW twins showed a shorter average path length (t(52) = −2.01, p = .05; Figure 2b bottom), indicating greater efficiency, relative to their lower BW co‐twin. Higher BW twins also showed a nonsignificant tendency toward less global efficiency (Figure 2a top; t(52) = −1.96, p = .06).

To investigate which nodes of the limbic network showed differences in connectivity as a function of BW, we ran linear mixed effects models for all limbic nodes, for both global efficiency and average path length, using the same model as above (see Table 1 for the results, follow‐up analyses performed at the p < .05 level, uncorrected). To quantify rs‐FC values, particular regions of the brain need to be active above a certain baseline threshold. For a subset of our participants, however, activity within a particular node did not reach threshold, thus resulting in missing path length values for some participants for some nodes (average missing data for DMN: 17%; limbic: 16%). We took the following conservative approach to account for the missing values: if one value was missing for a twin pair, we replaced the value with the value of their co‐twin. If both twins were missing a value, we replaced both values with the average of all present values. This strategy thus biased the data against finding a difference. The rs‐FC data revealed a shorter average path length for the higher BW twin relative to the lower BW twin for the right amygdala and right hippocampus (Table 1; Figure 3). In addition, we ran the same analyses specifically for our average path length values in the right amygdala and right hippocampus, removing those individuals without a value; our findings were confirmed (right amygdala, t(45.3) = −2.01, p = .05; right hippocampus, t(49.8) = −2.18, p = .03).

TABLE 1

Results from the linear mixed effects analyses comparing higher versus lower BW twins for the limbic network nodes, for the two graph theoretical measures

	Global efficiency		Average path length
	t	p	t	p
ACC	−0.54	.6	0.42	.7
Left IC	−0.46	.7	—	—
Right IC	−0.81	.4	—	—
Left OFC	−0.83	.4	−1.16	.3
Right OFC	−1.2	.24	−0.63	.5
Left hippo	−0.15	.9	—	—
Right hippo	−1.29	.2	−2.4	.02 a
Left amygdala	0.39	.7	0.46	.7
Right amygdala	−1.87	.07	−2.16	.04 a

Note: Degrees of freedom (df) = 52 for all analyses.

Abbreviations: ACC, anterior cingulate cortex; BW, birth weight; hippo, hippocampus; IC, insular cortex; OFC, orbitofrontal cortex.

Denotes significance at an alpha of .05 (p‐values are uncorrected for multiple comparisons). Empty cells = values with too little variability for reliable analysis (due to a high number of individuals with activity below baseline in the particular nodes).

FIGURE 3

Violin plots depicting average path length for the (a) right (r) hippocampus and (b) r amygdala as a function of birth weight (BW). * p < .05 uncorrected, black dots indicate the mean, and black vertical lines the error bars (SD). Values in the plot represent the means

When BW discordance (i.e., difference in BW for each twin pair) was regressed against discordance in rs‐FC (i.e., difference in rs‐FC for each twin pair), BW discordance was associated with right amygdala path length discordance (effect estimate [β] = −.33, t(51) = −2.13, p = .038), but not with the right hippocampus (p > .17).

Male twin‐pair analyses conceptually mirrored that of the full sample (global efficiency: t(42) = −1.6, p = .12, average path length: t(42) = −1.92, p = .06), with a less efficient limbic network for the low BW males compared to their higher BW co‐twin. For female pairs, there were no significant associations between BW and global efficiency (t(30) = −1.1, p = .27) and average path length (t(30) = −0.9, p = .37). Follow‐up analyses in the male pairs revealed that lower BW male twins had a less efficient right hippocampus compared to their higher BW co‐twins (t(42) = −2.44, p = .02) but no difference in efficiency for the right amygdala (t(42) = −1.54, p = .13). Associations were not significant for female pairs (right amygdala: t(30) = −1.66, p = .11; right hippocampus: t(30) = −1.12, p = .27).

When BW discordance for male twins was regressed against rs‐FC discordance, we found an association between discordance in BW and (a) left hippocampus efficiency (β = −.18, t(21) = −2.48, p = .02), (b) right hippocampus path length (β = −1.14, t(21) = −2.65, p = .015), and (c) left insular cortex path length (β = −.95, t(21) = −3.14, p = .005). Trends were also observed in male pairs for the association between BW discordance and right amygdala path length discordance (β = −.77, t(21) = −2.02, p = .056). In female pairs, associations between BW discordance and the preceding brain regions were not significant (left hippocampus efficiency: β = .01, t(21) = 0.20, p = .84; right hippocampus path length: β = −.03, t(21) = −0.15, p = .88; left insular cortex path length: β = .28, t(21) = 0.94, p = .35; right amygdala path length: β = −.21, t(21) = −1.38, p = .18, nor were there any significant associations in any of the other selected brain regions (p > .15).

SUMMARY

We found that rs‐FC varied as a function of prenatal environment in the limbic network but not the DMN. Two regions appeared to preferentially drive this finding, namely the right amygdala and right hippocampus, particularly in male pairs. To follow‐up, we subsequently examined underlying structural connectivity as a function of BW variation within this limbic network using DTI tractography measures.

DTI Tractography

DTI tractography measures were missing for one (female) twin pair. Thus, the DTI analyses were conducted on 104 participants (i.e., 52 twin pairs). Track count values were extracted for the connection between the right hippocampus and right amygdala (see Figure 4), along with the number of tracks connected to each of these regions separately (i.e., three measures).

FIGURE 4

A representative DTI image of tracts between the right amygdala and right hippocampus for one participant

Within the full sample, we did not find any difference in track counts as a function of BW (ps > .14). Track count differed as a function of sex (see Figure 5); males had greater track counts compared to females for the right hippocampus: t(50) = 2.22, p = .03; however, the right amygdala did not reach statistical significance: t(50) = 1.78, p = .08.

FIGURE 5

Violin plots depicting DTI tract count for the (a) path between the right amygdala and right hippocampus (path rHippo–rAmygd), (b) right hippocampus (rHippocampus), and (c) right amygdala (rAmygdala), split by sex and BW. * p < .05 uncorrected, † p = .05 uncorrected, ‡ p < .1 uncorrected. Black dots indicate the mean, and black vertical lines the error bars (SD). Values in the plot represent the means

Follow‐up tests in males indicated that BW was related to track count between the right hippocampus and right amygdala (t(21) = 2.08, p = .05). Specifically, higher BW male twins had more tracks connecting these nodes than their lower BW siblings. Such association was not found for the female twin pairs (t[29] = 0.043, p = .97).

GENERAL DISCUSSION

We studied whether variations in BW (indexing unique prenatal environment) within MZ adolescent twins were linked to variations in rs‐FC in the DMN and/or the limbic network. We found that the lower BW twins showed less efficient limbic network connectivity as compared to their higher BW co‐twin. This result was mainly driven by the right amygdala and right hippocampus, and primarily present in male twin pairs. Further targeted analyses using white matter connectivity measures and focusing on the structural connectivity between the right amygdala and right hippocampus provided some converging evidence, in that higher BW male twin pairs had larger track counts between the right amygdala and right hippocampus as compared to the lower BW male twins. Taken together, our data extend the current understanding of how unique prenatal environmental exposures may affect later limbic brain development, irrespective of genetics.

Most of the previous work investigating the link between BW and brain development utilized cortical morphometry (Casey et al., 2017; Haukvik et al., 2014; Levesque et al., 2015; Martinussen et al., 2005; Nagy et al., 2009; Raznahan et al., 2012; Walhovd et al., 2012). These studies found that BW is associated with brain structure in both extremely low BW individuals and those in the typical BW range. In MZ twin samples, although within‐pair differences are generally small, lower BW is linked with smaller brain volume (Levesque et al., 2015), decreased frontal‐limbic cortical surface area (Casey et al., 2017; Raznahan et al., 2012), increased cortical thickness (Casey et al., 2017), and decreased gray matter in the right superior frontal cortex (Levesque et al., 2015). The current study extends these findings to functional neural activity at rest, indicating that the unique prenatal environment relates to limbic network neural connectivity, with higher BW linked to more efficient limbic network communication.

We found that the association between BW and limbic network efficiency was primarily driven by two critical nodes of the limbic network: the right amygdala and right hippocampus. This may not be surprising, given that both brain regions show high functional and anatomical connectivity to prefrontal cortical areas (amygdala: Bzdock, Laird, Zilles, Fox, & Eickhoff, 2013; hippocampus: Libby, Ekstrom, Ragland, & Ranganath, 2012). While some research suggests that more than 30% of variability in connectivity between the amygdala and hippocampus can be explained by genetics, the two nodes have been shown to be independently influenced by environmental factors (Achterberg et al., 2018). The environmental influence on limbic network functional connectivity dovetails with our work showing that the prenatal environment is associated with hippocampus and amygdala connectivity. Indeed, exposure to prenatal stress could have deleterious effects on hippocampal structure and function (Frodl & O'Keane, 2013). At the neurochemical level, we have previously shown in a sample of (non‐twin) adults that BW is associated with lower serotonin—an important modulator of early brain development—in the hippocampus (Booij et al., 2012). In addition, the amygdala has been extensively linked to (neural vulnerability to) psychopathology (He, Xu, Zhang, & Zuo, 2016; Zhang et al., 2018). For example, greater global amygdala intrinsic connectivity was associated with lower trait anxiety (He et al., 2016). Altered rs‐FC of the amygdala has also been linked to adolescent depression (Cullen et al., 2014) and externalizing behaviors (Saxbe et al., 2018). Taken together, these findings suggest that the hippocampus and amygdala may be important mediators between early environment and future potential emotional problems.

While our data highlighted the central role both the right amygdala and hippocampus play in the relationship between BW and rs‐FC in the full sample, the data pattern was more consistent for the amygdala. Specifically, and unlike the right hippocampus, analyses with the right amygdala, as a function of BW, were consistent across both functional analyses, namely the one employing high versus low BW group comparison and the one using the within‐pair BW difference score. This finding may suggest that right amygdala rs‐FC may be a sensitive measure of environment‐related alterations in the developing brain. Indeed, cumulative evidence of rs‐fMRI studies of neural correlates of early‐life environment in singleton and twin populations indicated alterations in (right) amygdala‐based rs‐FC with other regions of the limbic neurocircuitry including insula, hippocampus, as well as orbitofrontal and anterior cingulate cortices (e.g., Grewen, Salzwedel, & Gao, 2015; Johns et al., 2019; Qiu et al., 2015; Salzwedel et al., 2019; Scheinost et al., 2016). Moreover, and unlike the hippocampus, since the right amygdala rs‐FC was linked to BW in both analyses, and in particular the within‐pair discordance analysis, we might argue that connectivity showed a greater influence of individual unique environment than genetic factors. Further research on heritability estimates of limbic rs‐FC is necessary to confirm this hypothesis. Lastly, prior research suggests that the impact of ones early‐life environment on the hippocampus may be minimized until adulthood (Lupien, McEwen, Gunnar, & Heim, 2009). Examining the current population of adolescents during adulthood would thus provide further data to test this hypothesis.

In contrast, no associations were found with the DMN. This may suggest that the DMN is primarily driven by genetic or shared environmental factors. The notion of high genetic influences in the DMN is in line with a study in MZ and dizygotic twins showing relatively large genetic contributions to global DMN network cost‐efficiency, with less clear contributions of genetics to regional activity (Fornito et al., 2011). In contrast, Fu et al. (2015) found primarily environmental influences on widespread neural connectivity, with the highest genetic influence found in sensory networks, and less genetic influences in cognition‐related networks such as the DMN. Taken together, more work is needed to better understand the specific genetic versus environmental contribution of DMN rs‐FC.

Previous research has indicated the importance of developmental timing in determining genetic or environmental contributions on the brain and on emotional behaviors (Lacourse et al., 2014; Schumann et al., 2017). At the behavioral level, longitudinal research in twins has shown that that the relative contributions of genetics and environment to aggression and negative affect (behaviors known to be modulated by limbic circuits) are age dependent (Lacourse et al., 2014; Schumann et al., 2017). Longitudinal studies in large twin samples, with repeated brain measures throughout development, may help to further clarify these age‐specific genetic and environmental contributions to brain structure and function. Such studies may also allow the investigation of the moderating role of additional intervening factors (e.g., social‐economic status) on within‐pair associations, over time.

Interestingly, the strongest associations between BW and rs‐FC/structural connectivity were found in males. Sex differences in brain development have been well documented in the literature (Asato, Terwilliger, Woo, & Luna, 2010; Kilpatrick et al., 2006). We found a more robust association between BW and limbic connectivity in males than in females, which is in line with some behavioral studies; for example, the association between BW and various types of externalizing behaviors are more pronounced in adolescent boys than in adolescent girls (Momany et al., 2017). As suggested by Momany et al. (2017), males may be more sensitive to prenatal adversities than females, potentially due in part to higher prenatal levels of testosterone, which, in turn, may contribute to sex dimorphism in brain development and subsequent behavioral outcomes.

Strengths of the study include the use of a well‐characterized, longitudinal sample of adolescents, and the complementary nature of two different analytical strategies for the primary analyses. Further, all twins were 15 years of age and were carefully screened for the absence of psychiatric history. The study was designed to investigate associations between prenatal environment and brain structure/function in a community sample whose BW fell within the normal range. While this then makes it difficult to generalize our findings to clinical populations, it also indicates that variation in BW within the general population may also contribute to variations in neural activity. To note, our analyses were uncorrected for multiple comparisons. While this may seem to be a pitfall of the current study, when studying healthy MZ twins with small BW differences, one should not expect to find very large differences in brain structure or function. To mitigate this limitation, we used a conservative approach to account for missing data, biasing the data away from rejecting the null hypothesis, and followed up with analyses in which we removed those twins with missing values, again confirming our findings. Furthermore, our two resting‐state networks were selected a priori, based on previous work. That being said, while the current study conducted in a healthy MZ twin sample contributes to the understanding of the role of in utero environment on functional default mode and limbic networks in typically developing adolescents, replication in a larger sample is required. A larger sample would also allow one to reliably investigate whether discordance in adolescent brain function predicts behavioral outcomes in adulthood.

Taken together, our findings suggest that subtle variations in BW are associated with changes in rs‐FC within the limbic network in healthy adolescent twins. By extending this method to clinical samples who have difficulty with emotion processing, we may glean further important insights into the functioning of the limbic network and into the clinical significance of our results.

CONFLICT OF INTEREST

The authors declare no conflicts of interest.