Distinctive heritability patterns of subcortical-prefrontal cortex resting state connectivity in childhood: A twin study

Connectivity between limbic/subcortical and prefrontal-cortical brain regions develops considerably across childhood, but less is known about the heritability of these networks at this age. We tested the heritability of limbic/subcortical-cortical and limbic/subcortical-subcortical functional brain connectivity in 7- to 9-year-old twins (N = 220), focusing on two key limbic/subcortical structures: the ventral striatum and the amygdala, given their combined influence on changing incentivised behavior during childhood and adolescence. Whole brain analyses with ventral striatum (VS) and amygdala as seeds in genetically independent groups showed replicable functional connectivity patterns. The behavioral genetic analyses revealed that in general VS and amygdala connectivity showed distinct influences of genetics and environment. VS-prefrontal cortex connections were best described by genetic and unique environmental factors (the latter including measurement error), whereas amygdala-prefrontal cortex connectivity was mainly explained by environmental influences. Similarities were also found: connectivity between both the VS and amygdala and ventral anterior cingulate cortex (vACC) showed influences of shared environment, while connectivity with the orbitofrontal cortex (OFC) showed heritability. These findings may inform future interventions that target behavioral control and emotion regulation, by taking into account genetic dispositions as well as shared and unique environmental factors such as child rearing.


Introduction
The contributions of limbic brain regions and the prefrontal cortex (PFC) to enhanced coordination in affective/motivational behaviors change considerably from childhood to adulthood (van Duijvenvoorde et al., 2016b). Resting State functional MRI (RS-fMRI) studies on limbic/subcortical-cortical functional brain connectivity in adults have provided insights into the connectivity patterns between different limbic/subcortical (sub) regions and the PFC, with positive connectivity between limbic/subcortical regions and affective PFC regions, and negative connectivity between limbic/subcortical regions and dorsal control regions of the PFC (Choi et al., 2012;Di Martino et al., 2008;Roy et al., 2009). Despite the consistent findings in general connectivity patterns in adults, not much is known about the robustness of these effects in children, and the role of genetic and environmental influences on limbic/subcortical-PFC brain connectivity. To date, the size of environmental and genetic contributions to limbic/subcortical-PFC connectivity has not been examined in children. In this study, we therefore investigated the robustness of findings regarding limbic/subcortical-PFC functional brain connectivity in childhood, and the heritability of these connections in 7-to-9-year-old twins (N ¼ 220). The current paper is the first to investigate childhood RS connectivity in two independent samples and additionally explore genetic and environmental influences on that connectivity, thereby providing important insights in the underlying mechanisms of functional brain connectivity in childhood.
RS-fMRI studies in adults have shown that the striatum is functionally connected to distributed regions throughout the entire brain, including motor, cognitive, and affective systems (Barnes et al., 2010;Choi et al., 2012;Di Martino et al., 2008). Different sub regions within the striatum show distinct functional connectivity patterns (Choi et al., 2012;Di Martino et al., 2008). A pioneering study of Choi et al. (2012) revealed distinct cortical-connectivity for five different sub regions in the striatum. For example, a dorsal sub region of the striatum was mainly connected to a network of the dorsolateral PFC (dlPFC), the dorsal medial PFC (dmPFC), and parietal regions, whereas a more ventral sub region of the striatum was primarily connected to medial/orbitofrontal regions of PFC (Choi et al., 2012;Di Martino et al., 2008). In the current study we focused on the ventral striatum, since this striatal sub region is consistently implicated in affective/motivational behavior (Haber and Knutson, 2010). Adult studies revealed that the ventral striatum is positively connected to limbic-affective regions such as the ventral medial PFC (vmPFC), the ventral anterior cingulate cortex (vACC), the orbitofrontal cortex (OFC), and the insula (Choi et al., 2012;Di Martino et al., 2008). In contrast, negative connectivity has been reported between the ventral striatum and cortical regions related to cognitive control, such as the dlPFC, the dorsal anterior cingulate cortex (dACC), the parietal cortex, and the precuneus (Di Martino et al., 2008). The amygdala also shows negative connectivity with dorsal cortical regions, including the dlPFC, dACC, dmPFC, the parietal cortex, and to the cerebellum (Roy et al., 2009). The positive connectivity patterns from the amygdala are ventrally oriented, including the vmPFC, the rostral ACC, and the OFC, but also more temporally oriented, towards the insula and inferior frontal gyrus (IFG) (Roy et al., 2009;Stein et al., 2007).
The development of limbic/subcortical-prefrontal cortex functional brain connectivity from childhood to adulthood has also been studied with RS-fMRI (e.g., Fareri et al. (2015), Gabard-Durnam et al. (2014), van Duijvenvoorde et al. (2016a), Fareri et al. (2015), Gabard-Durnam et al. (2014), van Duijvenvoorde et al. (2016a)). Developmental studies consistently report an overall shift from local limbic/subcortical-subcortical connectivity in childhood towards more distributed long-range limbic/subcortical-cortical connectivity in adulthood (Fair et al., 2009;Menon, 2013;Rubia, 2013;Vogel et al., 2010). However, this age-related shift from local to distributed connectivity was called into question after several studies had shown that these developmental changes were largely influenced by age-related changes in head-motion (Power et al., 2012;Van Dijk et al., 2010). That is to say, head motion can result in substantial changes in RS-fMRI connectivity (Power et al., 2012;Van Dijk et al., 2010). Specifically, volume-to-volume micro movement (i.e., head motion between two frames) can overestimate short-distance connectivity and underestimate long-distance connectivity (Satterthwaite et al., 2013). Young children usually have more difficulty lying still, resulting in more volume-to-volume micro movement, which may have resulted in an underestimation of subcortical-cortical brain connectivity in childhood. Therefore, there is a need to better understand connectivity patterns in childhood, using large samples and replication designs.
The PFC gradually develops both structurally and functionally until maturation in early adulthood (Lenroot and Giedd, 2006;van Duijvenvoorde et al., 2016a). Both the striatum and the amygdala show plasticity to the environment (for a review, see Tottenham and Galvan (2016)). For example, caregiving adversity during childhood (neglect, institutional care or low parental warmth) has been associated with amygdala hyper reactivity during adolescence (Casement et al., 2014;Garrett et al., 2012;Tottenham et al., 2011). In addition, adults and adolescents with a history of childhood stress show less striatum activity when receiving a monetary reward (Boecker et al., 2014;Goff et al., 2013;Hanson et al., 2016). Given these environmental influences on ventral striatum and amygdala activity, the connectivity between these limbic regions and cortical PFC regions may also be influenced by environmental factors. Alternatively, the high commonality of psychiatric disorders that rely on limbic/subcortical-PFC connections in families may suggest a heritability factor as well (Bouchard and McGue, 2003;Flint and Kendler, 2014). It is important to note that heritability estimates for brain anatomy and connectivity differ across development such that heritability estimates are stronger in adulthood than in childhood (Lenroot et al., 2009;van den Heuvel et al., 2013).
The few studies that examined these contributions in monozygotic (MZ) and dizygotic (DZ) twins in adults reported significant influences of genetics on functional connectivity, with little shared environmental influences (for a review see Richmond et al. (2016)), although some studies reported influences of both genetics and shared environment (Yang et al., 2016). Prior findings are mostly based on adult twin studies, whereas limbic/subcortical-PFC connectivity changes considerably during child and adolescent development. That is to say, functional connectivity from the ventral striatum and the amygdala with (medial) prefrontal regions increases substantially during development (Fareri et al., 2015;Gabard-Durnam et al., 2014;van Duijvenvoorde et al., 2016a). This increase in long range interactions between the ventral striatum, the amygdala, and the PFC may contribute to the improved ability of children to regulate behavior and emotions in the transition to adolescence (Casey, 2015;Ernst, 2014;Somerville et al., 2010). Together, these findings underscore the importance of studying heritability of RS brain connectivity in childhood. Taken together, the aims of the current study were to investigate (1) the robustness of limbic/subcortical-cortical and limbic/subcorticalsubcortical brain connectivity in childhood, and (2) the heritability of these connections in 7-to-9-year-old twins (N ¼ 220). We included 7-to-9-year-old twins since they are old enough to produce relatively good MRI data, while still representing (middle) childhood as a developmental phase. The study pursued two goals: 1) to investigate subcortical-cortical and subcortical-subcortical brain connectivity in childhood using two key limbic structures: the ventral striatum and the amygdala, and 2) to examine the heritability of these connections comparing MZ and DZ twins. We specifically focused on connectivity between limbic/subcortical regions and six PFC regions: the vmPFC, the vACC, the OFC, the dmPFC, the dACC and the dlPFC. These regions have been shown to be functionally connected to both the ventral striatum and the amygdala in adults (Di Martino et al., 2008;Roy et al., 2009) and display developmental changes related to increased cognitive control and emotion regulation (Casey, 2015;Ernst, 2014;Somerville et al., 2010), making them key targets to study in our sample.
The first question, regarding replicability of childhood RS connectivity, was addressed in two independent samples in order to examine connectivity patterns without genetic components. This allowed us to test for replication, thereby contributing to the debate about reproducibility of neuroscientific patterns (Open Science, 2015). Next, we specifically focused on RS-fMRI connectivity from the ventral striatum and amygdala to the six PFC regions and two additional subcortical regions (thalamus and hippocampus); since prior studies have shown that these regions show important developmental effects (Fareri et al., 2015;Gabard-Durnam et al., 2014). Based on prior studies, we expect to find replicable and robust resting state connectivity in childhood (Misic and Sporns, 2016), with distinctive patterns for ventral striatum and amygdala (Choi et al., 2012;Porter et al., 2015;Roy et al., 2009).
To address the second question, concerning the heritability of limbic connectivity, we compared MZ and DZ twin pairs using ACE modeling. This decomposition model provides an estimate of the proportions of the variance in the data that are attributed to heritable, shared environmental, and unshared/unique environmental factors. Previous studies have shown both influences of genetics (Richmond et al., 2016) and environmental contributions (Tottenham and Galvan, 2016), indicating that there could be an interplay between genetics and environment (Yang et al., 2016).

Participants
Participants were part of the Leiden Consortium on Individual Development (L-CID) twin study. Families with a same-sex twin pair born between 2006 and 2009, living within two hours travel time from Leiden, were recruited through the Dutch municipal registry and received an invitation by mail to participate. 256 families with a twin pair (512 children) were included in the L-CID study, of which 443 children underwent the RS scan (Table S1). The Dutch Central Committee on Human Research (CCMO) approved the study and its procedures (NL50277.058.14). Written informed consent was obtained from both parents. Families received financial compensation (€80.00) for their participation in the L-CID study. All participants were fluent in Dutch, had normal or corrected-to-normal vision, and were screened for MRI contra indications. All anatomical MRI scans were reviewed and cleared by a radiologist from the radiology department of the Leiden University Medical Center (LUMC). Three anomalous findings were reported and these participants were excluded. Participants' intelligence (IQ) was estimated with a verbal intelligence subtest (Similarities) and a performance intelligence subtest (Block Design) of the Wechsler Intelligence Scale for Children, third edition (WISC-III, Wechsler (1991)).
Since head motion can result in substantial changes in RS-fMRI connectivity (Power et al., 2012;Van Dijk et al., 2010), we investigated micro-movement using the motion outlier tool in FSL version 5.0.9 (FMRIB's Software Library, Smith et al. (2004)). Volumes with more than 0.5 mm framewise displacement (FD) were flagged as outliers. In line with recent studies (Couvy-Duchesne et al., 2014;Engelhardt et al., 2017), our twin analyses indicated that motion (amount of FD) was heritable. That is to say, there was a stronger correlation within MZ than DZ twins (r mz ¼ .44, p < .001; r dz ¼ .25, p ¼ .02). Behavioral genetic modeling of the amount of motion in the initial sample pointed towards genetic influences (A ¼ 38%, 95 confidence interval (CI): 26-56%, see Table S2). Children with more than 20% of their volumes flagged were excluded from further analyses (Power et al., 2012). In total, 209 participants (47.5%) were excluded based on excessive head motion. An additional 11 participants were excluded due to registration problems. The final sample consisted of 220 children (41% boys, mean age 8.00 AE 0.67, age range 7.02-9.08), of which 64 complete twin pairs (128 children, 58% MZ). There was no association between age and motion in the final sample (r ¼ .06, p ¼ .35). Moreover, there were no significant influences of heritability for head motion in the final sample (A ¼ 0%, 95% CI: 0-35%, see Table S2), implying that only more extreme motion is heritable, and this is not true of more subtle motion. For an overview of sample selection and dropout, see Table S1.
For the first set of analyses (examining replicability of childhood RS connectivity) we divided the sample into two subsamples of genetically independent individuals. Of the 64 complete twin pairs, we randomly chose either the youngest or oldest child within a twin pair. The other half of the twin pair was left out of the replication analyses. The replication sample therefore consisted of 156 (220-64) genetically independent children who were divided over two samples of N ¼ 78. Table 1 provides an overview of demographic characteristics, estimated IQ and motion in samples I and II. There were no significant differences in demographic characteristics between the samples (Table 1). Moreover, the distribution of gender did not significantly differ from chance (Sample I -45% boys, t (77 For the second set of analyses (testing heritability of childhood RS connectivity), we estimated the contributions of genetic and environmental factors to subcortical-cortical and subcortical-subcortical functional brain connectivity using behavioral genetic modeling on seed-ROI connections. The complete twin pairs were therefore divided in monozygotic (N ¼ 37) and dizygotic (N ¼ 27) twin pairs. Table 2 provides an overview of demographic characteristics, estimated IQ and motion in MZ and DZ twins. There were no significant differences in demographic characteristics between the samples (Table 2). For the twin samples, the distribution of gender significantly differed from chance, with the inclusion of fewer boys than girls in both samples (MZ -35% boys, t

Data acquisition
MRI scans were acquired with a standard 32 channel whole-head coil on a Philips Ingenia 3.0 Tesla MR system. Resting state data was acquired at the end of a fixed imaging protocol. Children were instructed to lie still with their eyes closed for 5 min. They were explicitly told not to fall asleep. To prevent head motion, foam inserts surrounded the children's    heads. A total of 142 T2 -weighted whole-brain echo planar images (EPIs) were acquired, including 2 dummy volumes preceding the scan to allow for equilibration of T1 saturation effects (scan duration 316.8 s; repetition time (TR) ¼ 2.2 s; echo time (TE) ¼ 30 ms; flip angle ¼ 80 ; field of view (FOV, in mm) ¼ 220.000 (rl) x 220.00 (ap) x 111.65 (fh); 37 slices). In addition, a high-resolution EPI scan was obtained for registration purposes (scan duration 46.2 s; TR ¼ 2.2 s; TE ¼ 30 ms, flip angle ¼ 80 , FOV ¼ 220.000 (rl) x 220.00 (ap) x 168.00 (fh), 84 slices), as well as a T1weighted anatomical scan (scan duration 296.6 s; TR ¼ 9.72 s; TE ¼ 4.59 ms, flip angle ¼ 8 , FOV ¼ 177.333 (rl) x 224.000 (ap) x 168.000 (fh), 140 slices). Since motion causes substantial artifacts within structural scans, we visually inspected the quality of the T1-weighted anatomical scan directly after acquisition. If the scan was affected by motion (blurry T1 image), we repeated the T1 scan. This was the case for 3% of the included participants.

Data preprocessing
The preprocessing of resting-state fMRI data was carried out using FMRIB's Expert Analysis Tool (FEAT; version 6.00) as implemented in FSL version 5.09 (Smith et al., 2004). The following preprocessing steps were used: motion correction (MCFLIRT; Jenkinson et al. (2002)), slice time correction, removal of non-brain tissue using the Brain Extraction Tool (BET; Smith (2002)), spatial smoothing using a Gaussian kernel of 6 mm full width at half maximum, and high-pass temporal filtering (Gaussian weighted least-squares straight line fitting, with sigma ¼ 100 s, 0.01 Hz cut-off). To register fMRI scans to standard space, each subject's functional scan was registered to the corresponding high resolution EPI scan, by using FMRIB's Linear Image Registration Tool (FLIRT, Jenkinson et al. (2002)). Next, an integrated version of boundary based registration (BBR; Greve and Fischl (2009)) was performed to improve the accuracy of the registration from high resolution EPI to subjects' structural space. Lastly, FMRIB's Nonlinear Imaging Registration Tool (FNIRT) with a 10 mm warp resolution was used to further refine registration from subjects' structural space to standard MNI-152 space (Jenkinson et al., 2002;Jenkinson and Smith, 2001). To ensure accurate alignment, we visually inspected the summery of the registration for all participants. Examples of correct and incorrect registration can be found in the supplementary materials ( Figure S1). In total, 11 participants were excluded due to registration problems (Table S1).

First-level seed based analysis
To investigate limbic/subcortical-cortical and limbic/subcorticalsubcortical functional brain connectivity we used two subcortical seeds: the ventral striatum (VS) and the amygdala (AMY). The VS seed was based on the "limbic striatum" of the Oxford-GSK-Imanova structural connectivity striatal atlas (Tziortzi et al., 2014). The AMY seed was based on the Harvard-Oxford subcortical structural atlas. Seeds were anatomical, bilateral and thresholded at !75% probability, resulting in a VS seed of 197 voxels and an AMY seed of 254 voxels (Fig. 1). To extract subject specific time series, seeds were first registered to subject space by using FLIRT (Jenkinson et al., 2002). The subject-specific seeds were then used to extract time series from preprocessed RS data.
First-level general linear models (GLM) were performed separately on time-series from each seed. The following nuisance signals were included: global signal, white matter (WM), cerebral spine fluid (CSF), 6 motion parameters and FD outliers. The global signal was included to reduce the influence of artifacts caused by physiological processes (i.e., cardiac and respiratory fluctuations) and scanner drifts (Birn et al., 2006;Fox and Raichle, 2007). In order to extract the time series for WM and CSF, we used subject specific WM and CSF masked, which were generated with FMRIB's Automated Segmentation Tool (FAST, Zhang et al. (2001)). Additionally, each frame with an FD outlier, (FD > 0.5 mm) was represented by a single regressor in the first-level GLM (see also Chai et al. (2014), Chai et al. (2014)). With this approach the amount of regressors is different between participants (ranging from 0 to 28). To account for this difference in first-level GLMs, the number of FD outliers (and thus the number of extra regressors) was added to the higher level statistical analyses as an additional covariate.

Higher-level seed based analysis
For both seeds, two higher-level group analyses were carried out using FMRIB's Local Analysis of Mixed Effects (FLAME) stage 1; one for sample I and one for sample II. Higher-level analyses were performed using FLAME stage 1 with automatic outlier detection and included the number of extra regressors induced by the FD outlier modeling as covariate of no interest. Corrections for multiple comparisons were thresholded with Gaussian Random Field Theory cluster-wise correction with a minimal Z > 3.09 (corresponding to p < .001) and cluster significance of p < .05. Next, we inspected the overlap between whole brain connectivity from sample I and sample II using conjunction analyses. Conjunction analyses were performed using the easythresh_conj script in FSL (Nichols et al., 2005), using the same threshold described for the previous analyzes (Z > 3.09, p < 0.05) in order to identify regions commonly connected in both samples.

Region of interest analysis
To further investigate limbic/subcortical-cortical and limbic/ subcortical-subcortical brain connectivity we examined the zstats in predefined ROIs. Since studies have shown that different regions of the PFC have distinct functions, we investigated six specific subdivisions of the PFC (Fig 4a): the ventral and dorsal medial prefrontal cortex (vmPFC, dmPFC), the orbitofrontal cortex (OFC), the dorsal lateral prefrontal cortex (dlPFC), and the ventral and dorsal anterior cingulate cortex (vACC, dACC). All ROIs were bilateral. Regions were based on the Harvard-Oxford cortical structural atlas and were thresholded on !25% probability, resulting in the following sizes of anatomical ROIs: vmPFC 1189 voxels; dmPFC 5378 voxels; OFC 3502 voxels; dlPFC 5741 voxels; vACC 1313 voxels; and dACC1925 voxels. The following regions were used: Frontal Medial Cortex for vmPFC, Superior Frontal Gyrus for dmPFC, Frontal Orbital Cortex for OFC, Middle Frontal Gyrus for dlPFC, and the Cingulate Cortex anterior division for the ACC. The ACC was divided in a dorsal and ventral division with a cutoff at y ¼ 30.
Since both the VS and AMY also have shown to be connected the hippocampus (HPC) and the thalamus (TH) (Fareri et al., 2015;Gabard-Durnam et al., 2014;Roy et al., 2009), we included exploratory analyses of limbic/subcortical-subcortical connectivity, with additional subcortical ROIs of the TH and HPC (Fig 4b). Regions were based on the Harvard-Oxford subcortical structural atlas and were thresholded on !75% probability, resulting in a bilateral, anatomical TH ROI of 1646 voxels and a HPC ROI of 494 voxels. We used a stricter probability for the subcortical regions in order to prevent subcortical regions would overlap. In addition, we investigated functional connectivity between the VS and AMY. Zstats were extracted from subjects' specific first level for each seed with the different ROIs as a mask using Featquery (as implemented in FSL v5.09). This way we extracted subject-specific connectivity estimates for 12 different subcortical-PFC connections and 5 different subcortical-subcortical connections.
To explore possible outliers, we calculated z-values of the subject specific zstats at the group level. When outliers were detected (Z-value <-3.29 or >3.29), scores were winsorized (Tabachnick and Fidell, 2013). One sample t-tests were used to investigate whether connectivity between a seed and a ROI was significantly different from zero (separately for both samples). Independent sample t-tests were used to test whether there were differences in connectivity between sample I and II. Paired sample t-tests were used to test whether there were differences in connectivity between ROIs and the VS and AMY seeds.

Genetic modeling
Within the final sample (N ¼ 220), there were 64 complete twin pairs (37 MZ and 27 DZ, Table 2). Zygosity was determined by DNA analyses. DNA was tested with buccal cell samples collected via a mouth swab (Whatman Sterile Omni Swab). Buccal samples were collected directly after the MRI session, thereby ensuring that the children had not eaten for at least one hour prior to DNA collection.
Similarities among twin pairs can be due to shared genetic factors (A) and shared environmental factors (C), while dissimilarities are ascribed to unique environmental influences and measurement error (E), see Fig  S2. Behavioral genetic modeling with the OpenMX package (Neale et al., 2016) in R (R Core Team, 2015) provides estimates of these A, C, and E components. Since several heritable psychiatric disorders are associated with limbic/subcortical-PFC connections (Bouchard and McGue, 2003;Flint and Kendler, 2014), VS and AMY connectivity might also be heritable. However, these regions have also shown plasticity to the environment (Tottenham and Galvan (2016)), which could indicate influences of (shared or unique) environment. Therefore, we calculated the ACE models for each of the 17 seed-ROI connections and report the point estimates and 95% confidence intervals of A, C and E. High estimates of A indicate that genetics play an important role, whilst C estimates indicate influences of the shared environment. If the E estimate is the highest, variance in connectivity is mostly accounted for by unique environmental factors and measurement error. Comparisons of the ACE models with more parsimonious models (AE model, CE model, and E model) are described in the Supplementary Materials.

Whole brain analyses
First, we performed whole brain analyses for the subcortical seeds (VS and AMY) in sample I and II. Next we investigated the overlap between the two samples by using conjunction analyses.

Ventral striatum
Whole brain functional connectivity with the VS as seed for sample I is displayed in Fig. 2a (left top panel) and Table S3. Whole brain results for sample II are displayed in Fig. 2a (right top panel) and Table S4. To formally assess which connectivity patterns replicated across samples, conjunction analyses were performed. As visualized in Fig. 2a, whole brain VS connectivity in the two samples showed pronounced consistent positive connectivity with vACC, vmPFC, thalamus, insula, inferior temporal gyrus, parietal operculum cortex, putamen, pallidum, caudate, nucleus accumbens, amygdala, and the OFC (Table 3). Negative connectivity was consistent over two samples between VS and dACC, dlPFC, paracingulate gyrus, para-hippocampus, and hippocampus (Table 3).

Amygdala
Whole brain functional connectivity with the AMY as seed for sample I is displayed in Fig. 2b (left top panel) and Table S3. Whole brain results for sample II are displayed in Fig. 2b (right top panel) and Table S4. As visualized in Fig. 2b, whole brain AMY connectivity patterns showed overlap across the two samples, showing pronounced positive connectivity with the thalamus, pallidum, putamen, caudate, hippocampus, para-hippocampus, brainstem, frontal pole, insula, inferior frontal gyrus (IFG), fusiform cortex, and superior temporal gyrus (STG) (Table 3). Moreover, we found consistent negative connectivity between AMY and dmPFC, dlPFC, paracingulate gyrus, precuneus cortex, parietal cortex, posterior cingulate cortex, and lateral occipital cortex (Table 3).

Post-Hoc examination of subcortical-cortical connectivity
We investigated limbic/subcortical-cortical brain connectivity in more detail by visualizing connectivity patterns between subcortical seeds (VS and AMY) and prefrontal cortical ROIs of the vmPFC, dmPFC, vACC, dACC, OFC, and dlPFC. Connectivity patterns replicated across sample I and II, with the exception of VS-dmPFC and AMY-vACC connectivity (Fig. 3a, Table S5). Overall, subcortical regions exhibited positive connectivity with ventral cortical regions (vmPFC, vACC, OFC) and negative connectivity with dorsal cortical regions (dmPFC, dACC, dlPFC), see Fig. 3a. Paired sample t-tests were used to investigate differences in VS-PFC and AMY-PFC connectivity. For the vmPFC and vACC, positive connectivity with the VS was significantly stronger than connectivity with AMY (Table 4). Note that connectivity between AMY and the vmPFC and vACC was not significantly different from zero in one of the samples (Table S6). There were no differences between the VS and the AMY in connectivity with the OFC. The VS and AMY showed pronounced negative connectivity with dorsal cortical regions (Fig. 3a). For the dlPFC and dmPFC, negative connectivity with the AMY was significantly stronger than connectivity with the VS (Fig. 3a, Table 4). Note that connectivity between VS and the dmPFC was not significantly different from zero in one of the samples (Table S6). Connectivity between dACC and AMY was stronger than connectivity between dACC and VS in sample II, but not in sample I (Table 4). There were no significant gender or agerelated differences in subcortical-cortical connectivity (sample I and II combined).

Post-Hoc examination of subcortical-subcortical connectivity
To investigate limbic/subcortical-subcortical brain connectivity in more detail, we used two additional ROIs of the HPC, TH. Moreover, we investigated connectivity between the VS and the AMY. Connectivity patterns replicated across sample I and II (Fig. 3b, Table S6). The overall pattern showed pronounced positive connectivity between subcortical regions, see Fig. 3b. Interestingly, the HPC ROI showed strong positive connectivity with AMY (Fig. 3b, Table 4). More stringent thresholded (smaller) HPC ROIs resulted in similarly strong positive connectivity patterns (see supplementary materials, Fig S3), indicating that this strong connectivity was not inflated by cross-boundary blurring. VS-Hippocampus showed negative connectivity (Fig. 3b, Table 4), however, note that VS-HPC connectivity was not significantly different from zero in Sample II (Table S6). VS-TH connectivity was significantly stronger than AMY-TH connectivity, which was negative, and not significantly different from zero in sample II (Table S6). The connectivity estimate between the VS and AMY was small and not significantly different from zero in both samples (Fig. 3b, Table S6). There were no significant gender differences in limbic/subcortical-subcortical connectivity (sample I and II combined). We found weak negative correlations between age and VS- Table 3 MNI coordinates and local maxima for whole brain connectivity clusters from conjunction analyses (Sample I and Sample II) with Z > 3.09, p < .05 cluster correction. Anatomical regions were derived from the Harvard-Oxford atlas in FSL.   Heritability of subcortical-cortical connectivity An overview of ACE models for limbic/subcortical-cortical brain connectivity between seed (VS and AMY) and cortical ROIs (vmPFC, vACC, OFC, dmPFC, dACC, dlPFC) is provided in Table 5. Comparisons of the full ACE model with more parsimonious AE, CE and E models are displayed in Table S7 (VS) and Table S8 (Amygdala). Note that the estimates of the different components add up to 1 (100%). The overall pattern showed that the variance in VS-PFC connectivity was best accounted for by genetic and unique environmental factors (including measurement error). That is to say, the A estimate was moderately high for connectivity between VS and vmPFC (A ¼ 67%, Table 5. In addition to genetic influences, VS-vACC connectivity also showed influences of shared environment (A ¼ 12%, C ¼ 17%, E ¼ 71%). Variance in AMY-dorsalPFC connectivity was less influenced by genetics, with small contributions of the A component for connectivity between AMY and dmPFC ( . AMY-vACC connectivity showed moderately high estimates of the shared environment (C ¼ 35%, E ¼ 65%), with no influence of genetics (A ¼ 0%). AMY-vmPFC connectivity showed moderate influences of genetics (A ¼ 23%, C ¼ 0%, E ¼ 77%), and AMY-OFC connectivity showed high heritability Table 5.

Heritability of subcortical-subcortical connectivity
An overview of ACE models for limbic/subcortical-cortical brain connectivity between seed (VS and AMY) and the subcortical ROIs (HPC, TH, AMY) is provided in Table 6. Comparisons of the full ACE model with more parsimonious AE, CE and E models are displayed in Table S9 Note that the estimates of the different components add up to 1 (100%). Subcortical-subcortical connectivity was moderately influenced by genetics, with A estimates ranging from 32 to 42% (VS-HPC A ¼ 37%, E ¼ 63%; VS-AMY A ¼ 42%, E ¼ 58%; AMY-HPC A ¼ 32%, E ¼ 68%; AMY-TH A ¼ 35%, E ¼ 65%), and no influence of the shared environment (C ¼ 0%), with the exception of VS-TH connectivity, which was mostly influenced by environmental factors (A ¼ 4%, C ¼ 15%, E ¼ 81%), see Table 6.

Discussion
We investigated genetic and environmental influences on limbic/ subcortical-cortical and limbic/subcortical-subcortical RS-fMRI in a relatively large sample of 7-to-9-year-old MZ and DZ twins. As a complement to prior studies of genetic and environmental influences in adults (for example, Yang et al. (2016)), here we assessed twin concordance in children during a time of rapid development of these connections.

Replicability of childhood resting state connectivity
First we addressed childhood resting state brain connectivity, by studying patterns of connectivity from the ventral striatum and the amygdala, in two genetically independent samples. Reassuringly, and consistent with adult research (Misic and Sporns, 2016;Power et al., 2010;Thomason et al., 2011), we observed strongly replicable brain connectivity patterns over two samples of 7-to-9-year-old children, both in the whole brain seed based analyses and in the post-hoc ROI analyses. The general patterns showed positive connectivity between amygdala and ventral striatum and orbitofrontal cortex; and negative connectivity between these limbic/subcortical regions and dorsal medial and lateral regions. Previous studies showed that orbitofrontal cortex is more strongly involved in affective processes, whereas dorsal medial and lateral prefrontal cortex is more strongly associated with behavioral control, and the current findings fit with the hypothesized top-down control of dorsal lateral prefrontal cortex over the limbic subcortical brain regions (Casey, 2015;Ernst, 2014;Somerville et al., 2010).
In line with adult striatal-cortico connectivity patterns we found positive connectivity between the ventral striatum and vACC, vmPFC, and OFC (Di Martino et al., 2008), suggesting that these connections are already in place during middle childhood. The post-hoc ROI analyses indicated negative connectivity between the VS and the dACC, dlPFC and dmPFC, but these were less pronounced in the whole brain analyses. The difference between the current results and the connectivity patterns in adults could be due to developmental processes, since dorsal medial and lateral PFC regions continue to develop throughout adolescence (Casey, 2015;Ernst, 2014). Moreover, these differences in results might derive from the differences in limbic/subcortical seed regions. To date there is no consensus about the different sub regions of the striatum and different studies have used different approaches. Prior studies have suggested a more detailed subdivision of the striatum with, for example, additional distinctions within the ventral striatum (Choi et al., 2012;Di Martino et al., 2008). For the current paper we specifically chose only the ventral striatum, since this striatal sub region is specifically associated with developmental differences in affective/motivational behaviors. Future research could shed light on developmental differences in connectivity from different sub regions within the striatum, by directly comparing children and adults, using the same methodology in both samples (as was previously done for the VS by Fareri et al. (2015)).
Regarding amygdala-cortico connectivity, our developmental results were generally in line with the findings in adults. That is, we found positive connectivity with the OFC, the insula and the IFG, and negative connectivity with the dlPFC, dACC, dmPFC and parietal cortex (Roy et al., 2009;Stein et al., 2007). This is also in line with previous findings spanning ages from childhood to adulthood, showing that amygdala connectivity over development was largely stable . We did, however, find differences in amygdala-cerebellum connectivity compared to results in adults (Roy et al., 2009). Our whole brain analyses revealed a band of positive connectivity from the amygdala through the brainstem to the dorsal cerebellum, whereas adult results showed negative connectivity between the amygdala and the dorsal cerebellum (Roy et al., 2009). Interestingly, a recent study on amygdala functional connectivity in 4-to-7-year-old children also showed positive connectivity between amygdala and the cerebellum (Park et al., 2018). We submit that this is a developmental effect, reflecting positive connectivity to the dorsal cerebellum in childhood that becomes negative over development. Indeed age dependent changes in amygdala connectivity have been documented, with increasingly negative connectivity between the amygdala and cerebellum with increasing age . Notably, a recent cross-sectional longitudinal study of Jalbrzikowski et al. (2017) reported strong amygdala-mPFC connectivity in childhood, which declined to zero by adulthood (age range 10-19). However, we did not find strong amygdala-vmPFC connectivity in neither of the samples. This could be due to differences in age ranges, differences in the amygdala and vmPFC sub regions that were examined, as well as methodological differences in RS-fMRI analyses. In the current paper, we chose to use the whole amygdala as seed, to strike a balance between completeness and the number of connections and additional genetic analyses. However, it should be noted that the amygdala is not a single unit, but consists of several nuclei (Ball et al., 2007;Roy et al., 2009). Some studies have shown distinct connectivity patterns from different amygdala sub nuclei in adults (Roy et al., 2009), and over development .
In sum, our results showed robust and replicable whole brain connectivity in children, for the amygdala as well as the ventral striatum. In addition to previous studies that have shown that limbic/subcorticalcortical connectivity increases during adolescence (Fair et al., 2009;Gabard-Durnam et al., 2014;Menon, 2013;Rubia, 2013;Vogel et al., 2010); the findings from this study show that the vast architecture of this connectivity is already present before adolescence.

Heritability of childhood resting state connectivity
The second aim of this study was to examine the heritability of childhood resting state connections, specifically focusing on connections between the ventral striatum and amygdala with prefrontal cortex and other subcortical regions. Variance in the majority of connections from the ventral striatum to the prefrontal cortex was best described by genetics, with moderately strong heritability factors (up to 67%). Weaker ventral striatum-prefrontal cortex connections have been linked to psychiatric disorders such as depression (Russo and Nestler, 2013) and substance abuse (Deadwyler et al., 2004), which are thought to have a genetic component (Bouchard and McGue, 2003;Flint and Kendler, 2014). The association between genotypic characteristics and psychiatric disorders might be mediated by genetically based connectivity in the brain (Hyman, 2000). Interestingly, connectivity from the ventral striatum to the vACC and thalamus was mostly influenced by shared and unique environmental factors, which is in line with previous findings that reported environmental plasticity of the striatum (Tottenham and Galvan, 2016). These results suggest that long-range cortical-striatal connectivity is more strongly influenced by genetic profiles, while short range thalamic and vACC connectivity is more influenced by environmental factors.
With the exception of ventral striatum-thalamic connectivity, limbic/ subcortical-subcortical connectivity was notably influenced by genetics, with heritability estimates ranging from 32 to 42%. For instance, we found heritability for amygdala-hippocampus connectivity (A ¼ 32%), indicating that this emotional memory network (Phelps, 2004) is influenced by genetic factors. Interestingly, a broad literature has shown that these two regions independently are affected by environmental influences such as stress and early adversity (Barch et al., 2016;Lupien et al., 2009;Tottenham and Sheridan, 2009). This raises new questions with respect to how the amygdala-hippocampus circuitry is shaped and develops during child development. Moreover, while ventral striatum-prefrontal cortex connective showed large genetic influences, amygdala-prefrontal cortex connectivity showed mostly effects of the environment, with high estimates of the E component (up to 92%). There were two exceptions to this general pattern. First, in line with the ventral striatum, amygdala-vACC connectivity showed influences of the shared environment. The vACC has been shown to signal for socially salient cues such as peer feedback, both in adults as well as in children (Achterberg et al., 2016Somerville et al., 2006). Connectivity between the vACC and limbic/subcortical regions might also be susceptible to social context and social environmental factors, as these connections are significantly influenced by environment . Secondly, 54% of the variance in amygdala-OFC connectivity was explained by genetic influences. Interestingly, Whittle et al. (2014) have reported longitudinal effects of positive parenting on structural development of the amygdala and OFC. Our study is the first to show that variance in amygdala-OFC functional connectivity in childhood is explained by genetic factors. This finding has important implications for intervention research: Certain genetic profiles might be more susceptible to environmental influences than others, as is proposed by the differential susceptibility theory (Bakermans-Kranenburg and van Ijzendoorn, 2007;Ellis et al., 2011). A next step could be to examine whether children with specific genetic profiles are more susceptible to both the adverse effects of unsupportive environments and the beneficial effects of supportive rearing (see the study protocol of Euser et al. (2016)). Important aspects to take into account in those studies are the developmental differences in heritability estimates for brain anatomy and connectivity (Lenroot et al., 2009;van den Heuvel et al., 2013). That is, previous studies have found lower heritability estimates in children than in adults (van den Heuvel et al., 2013). However, the literature on heritability of functional brain connectivity is still relatively sparse, and most studies have examined whole brain RS and/or used different RS methods (Colclough et al., 2017;Ge et al., 2017;Glahn et al., 2010;Richmond et al., 2016;Yang et al., 2016), making comparisons between studies difficult. Studying differences in heritability estimates between children and adults, nevertheless, is an important issue for future studies, providing important insights in the developmental phase during which connections might be most sensitive to environmental influences.
Overall, the patterns of genetic and environmental influences for ventral striatum and amygdala were distinct: Long-range PFC connectivity with the ventral striatum was genetically influenced, whereas longrange amygdala connectivity was mostly environmentally influenced. These results may be the starting point for a better understanding of how brain development is both biologically based and environmentally driven.

Methodological considerations
Some methodological considerations should be noted. First, due to excessive motion, we had to exclude almost half of our initial sample. Nevertheless, due to our large sample size we could still perform analyses on a relatively large group of children, thereby increasing the statistical power of our analyses. It should be noted that the current standard of remaining motion in (adult) RS studies is even stricter, often using a cutoff of 0.3 mm FD. However, in terms of motion, the current results are based on a very clean dataset compared to earlier developmental studies. After exclusion of participants with excessive motion the gender distribution was significantly different from chance in the MZ and DZ twin samples, with more girls than boys included. Although there were no significant differences in gender between the MZ and DZ samples, and therefore this gender distribution is unlikely to have influenced our results, future studies on heritability of brain measures in childhood should opt to oversample young boys, since our results show the highest attrition rate in boys. Secondly, even after controlling for motion and including additional regressors with CSF and WM signals, our whole brain analyses show minimal but potentially artefactual correlations with non grey matter tissue. Future studies could include additional analytic steps to further minimize these effects, for example by controlling for cortical signal bleeding, i.e., regressing out signal from surrounding voxels Choi et al., 2012).
Third, we included the global signal as nuisance signals to reduce artifacts of cardiac and respiratory fluctuations and scanner drifts (Birn et al., 2006;Fox and Raichle, 2007), however, inclusion of global signal regression can introduce negative correlations between regions (Murphy et al., 2009) and therefore the interpretation of these negative connectivities should be done with caution.
Fourth, some of our genetic analyses of neural responses resulted in high estimates for the E component (up to 92%), reflecting influences from the unique environment and measurement error. The statistical power of genetic studies is influenced by, amongst others, the sample size (Verhulst, 2017;Visscher, 2004). Although our sample size can be considered relatively large for a developmental RS-fMRI study, it is modest for behavioral genetic modeling. Our sample size may have been insufficient to detect significant contributions of A (genetics) and C (shared environment), resulting in inflated estimates of the E component. Future studies should try to discriminate between the influence of unique environment and measurement error, for example by accounting for intra-subject fluctuations using repeated measures, as has recently been described by Ge et al. (2017).
Lastly, the current study made use of post hoc ROI analyses to further investigate limbic/subcortical-cortical connectivity, based on structural brain atlases. Although recent studies have provided functional atlases of the brain (Choi et al., 2012;Yeo et al., 2011), these are based on adults. To our best knowledge, there are no functional atlases based on developmental samples, and the vast majority of developmental studies have used anatomical regions to mask and/or extract functional connectivity (Fareri et al., 2015;Gabard-Durnam et al., 2014;van Duijvenvoorde et al., 2016a). By using these structural ROIs our results can be compared or combined with previously published studies. Nevertheless, we acknowledge that the functional architecture of the brain does not follow structural subdivisions, and this may be considered as a limitation of the current design.

Conclusions
Taken together, this study was the first to investigate twin effects in subcortical-subcortical and subcortical-cortical RS-fMRI in children, providing important insights in genetic and environmental influences on childhood brain connectivity. The behavioral genetic analyses showed moderate to substantial heritability of striatum-prefrontal cortex brain connectivity, and environmental influences on amygdala-orbitofrontal cortex connectivity, with implications for our understanding of the etiology of disorders that are associated with disrupted connectivity, such as drug abuse and depression. Prior studies have mainly estimated heritability for brain connectivity in adults (Yang et al., 2016), whereas child development provides unique possibilities for understanding the role of shared environment (Polderman et al., 2015). Examining how limbic/subcortical brain regions are functionally connected to the prefrontal cortex and whether a positive childrearing environment can foster these connections are important issues to address in future research. The current findings provide the first step in laying the groundwork for understanding genetic and environmental influences in shaping brain connectivity and may be the starting point for a better understanding of how brain development is both biologically based and environmentally driven.