Genetic and environmental influences on functional connectivity within and between canonical cortical resting-state networks throughout adolescent development in boys and girls

The human brain is active during rest and hierarchically organized into intrinsic functional networks. These functional networks are largely established early in development, with reports of a shift from a local to more distributed organization during childhood and adolescence. It remains unknown to what extent genetic and environmental influences on functional connectivity change throughout adolescent development. We measured functional connectivity within and between eight cortical networks in a longitudinal resting-state fMRI study of adolescent twins and their older siblings on two occasions (mean ages 13 and 18 years). We modelled the reliability for these inherently noisy and head-motion sensitive measurements by analyzing data from split-half sessions. Functional connectivity between resting-state networks decreased with age whereas functional connectivity within resting-state networks generally increased with age, independent of general cognitive functioning. Sex effects were sparse, with stronger functional connectivity in the default mode network for girls compared to boys, and stronger functional connectivity in the salience network for boys compared to girls. Heritability explained up to 53% of the variation in functional connectivity within and between resting-state networks, and common environment explained up to 33%. Genetic influences on functional connectivity remained stable during adolescent development. In conclusion, longitudinal age-related changes in functional connectivity within and between cortical resting-state networks are subtle but wide-spread throughout adolescence. Genes play a considerable role in explaining individual variation in functional connectivity with mostly stable influences throughout adolescence.


Introduction
The human brain is active during rest (Biswal et al., 1995(Biswal et al., , 1997. Data-driven approaches have been applied to resting-state functional MRI scans to obtain spatial patterns of temporally coherent signals that divide the brain into distinct intrinsic functional networks (DeLuca et al., 2005;Fox et al., 2005;Power et al., 2011;van den Heuvel and Hulshoff Pol, 2010;Yeo et al., 2011). The hierarchical organization of these functional networks is already present around birth. Primary functional networks, such as the sensorimotor, visual, and auditory networks are the first to develop in utero (Gilmore et al., 2018;Keunen et al., 2017;Thomason et al., 2015). After birth, the default mode network (DMN), dorsal mode network (DAN), and salience network (SN) mature into "adult-like" networks by the age of two years (Gao et al., 2011;Gilmore et al., 2018;Keunen et al., 2017). The executive control network (ECN) matures later on in life, in line with the protracted development of executive functions during childhood and adolescence (Gilmore et al., 2018;Zhang et al., 2017). These functional networks can be reliably identified in children and adolescents aged 9-15 years for both short-term (i.e. consecutive scan sessions) and long-term repeated measures at 2.5 years interval (Thomason et al., 2011). Thus, by adolescence, these spatially distributed and functionally linked brain regions that share information already closely resemble their adult state.
Cross-sectional studies have provided indications that functional connectivity of the human brain is undergoing subtle alterations during childhood and adolescence (for reviews see Cao et al., 2016;Ernst et al., 2015;Grayson and Fair, 2017;Stevens, 2016). Based on these cross-sectional studies, it is generally believed that the functional brain shifts from a local to a more distributed organization (Cao et al., 2016;Ernst et al., 2015;Fair et al., 2009;Menon, 2013). This is supported by decreases in functional connectivity separating functionally distinct regions (i.e. segregation), and increases in functional connectivity improving communication between functionally related regions (i.e. integration) (Cao et al., 2014;Dosenbach et al., 2010;Fair et al., 2009Fair et al., , 2008Gu et al., 2015;Kelly et al., 2009;Marek et al., 2015;Sato et al., 2014;Supekar et al., 2009;Uddin et al., 2011;Wig, 2017). The processes of segregation and integration are reflected in graph theoretical metrics by a decrease in local clustering coefficient, an increase in modularity, and an increase in global efficiency, and are furthermore accompanied by the emerging of hubs of increasing importance (i.e. consolidation of the network into rich-club networks) that shift from primary to higher order cortical regions (Cao et al., 2016(Cao et al., , 2014Grayson et al., 2014;Hwang et al., 2013;Sato et al., 2014Sato et al., , 2015Supekar et al., 2009;Wu et al., 2013;Zuo et al., 2011). However, results are inconsistent regarding the direction of change and affected regions (Stevens, 2016). In part, this may be due to the limited ability of cross-sectional studies to control for inter-individual variation (i.e. the "cohort effect") and are thereby restricted in their interpretation of "true" development (i.e. within subject developmental trajectories). In contrast, longitudinal studies acquire repeated measures of the same individuals and can utilize these measures as control to measure development changes over time within the individual (Crone and Elzinga, 2015;Mills and Tamnes, 2014;Telzer et al., 2018). Longitudinal studies on resting-state or task-regressed functional connectivity in typically developing children and adolescents (aged 9-15 years) reveal high levels of consistency and stability of functional connectivity estimates within and between several cortical resting-state networks over a 2-3 years interval (Thomason et al., 2011). There are reports of longitudinal age-related increases in functional connectivity (or integration) within several networks (Bernard et al., 2016;Long et al., 2017;Sherman et al., 2014;Sylvester et al., 2017;Wendelken et al., 2017Wendelken et al., , 2016. However, age-related decrease in functional connectivity (or segregation) within (Sylvester et al., 2017;Wendelken et al., 2016) and between networks have also been reported (Sherman et al., 2014), as well as mixed results reported for cortical-subcortical connectivity (Jalbrzikowski et al., 2017;Peters et al., 2017;Strikwerda-Brown et al., 2015). Thus, although functional brain connectivity during childhood and adolescence is largely stable and adult-like, there are several indications from longitudinal studies of reorganization of functional cortical networks during childhood and adolescent development (Table 1).
Genes partly control individual differences in brain functioning, at least in adults (Blokland et al., 2012;Douet et al., 2014;Jansen et al., 2015;Thompson et al., 2010;Richmond et al., 2016;Thompson et al., 2013). Heritability estimates for functional connectivity within the default mode network range from 10% to 80%, depending on the population and methodology used, and typically identify connections involving the posterior cingulate cortex (PCC) and bilateral parietal cortices (LLP and RLP) as strongest heritable connections (Adhikari et al., 2018a;Fu et al., 2015;Ge et al., 2017;Glahn et al., 2010;Korgaonkar et al., 2014;Meda et al., 2014;Sudre et al., 2017;Yang et al., 2016, Table 2). Findings in children and adolescents are still sparse. We were among the first to report that genes explain up to 40% of individual difference in brain activity during resting-state at the age of 12 years (van den Heuvel et al., 2013). These findings were confirmed and extended for cortical-subcortical connections in younger children, aged 7-9 years, with heritability ranging from 32% to 67% (Achterberg et al., 2018). And in 16-year-old adolescents reporting peaks of local clusters with heritability ranging between 55% and 83%but note that several cortical resting-state networks revealed overall low heritability estimates <10% (Fu et al., 2015). In the only longitudinal twin study on functional connectivity to date, in infants from birth to 2 years, age-dependent genetic effects on functional connectivity within cortical networks were found Table 1 Studies on longitudinal development of cortical resting-state or task-regressed functional connectivity in typically developing children and adolescents; ordered by midrange age at baseline of each cohort.  (Gao et al., 2014). It is unknown whether these age-dependent dynamic influences of genes on functional connectivity extend into childhood and adolescence in the absence of any longitudinal twin studies during this developmental period.
Here we report on genetic and environmental influences on functional connectivity within and between eight canonical cortical restingstate networks in a longitudinal adolescent cohort of twins and their older siblings measured on two occasions (mean ages of twins 12 and 17 years; mean ages of siblings 15 and 20 years; mean ages combined 13 and 18 years). We utilize a model that accounts for measurement imprecision by analyzing data from split-half sessions to obtain a reliable component of functional connectivity. This is the first longitudinal study on cortical resting-state networks to estimate the importance of genetic factors for functional connectivity within and between cortical resting-state networks during adolescence. The longitudinal data allowed us to investigate possible dynamic influences of genetic factors throughout adolescence (Teeuw et al., 2018;van Soelen et al., 2012b). We investigated the effects of sex and age on functional connectivity while controlling for measurement imprecision and residual head motion. Finally, we investigated the relation between intelligence and functional development of resting-state networks.

Participants
This project is part of the longitudinal BrainSCALE study on development of brain and cognition in twins and their older sibling (van Soelen et al., 2012a), a collaborative project between the Netherlands Twin Register (NTR; Boomsma et al., 2006;van Beijsterveldt et al., 2013) at the Vrije Universiteit (VU) Amsterdam and University Medical Center Utrecht (UMCU). The BrainSCALE cohort is a representative sample of typically-developing children from the Dutch population. A total of 112 families with twins and an older sibling participated in the study. The twins and siblings were assessed with a battery of cognitive and behavior tests and extensive neuroimaging protocol at baseline assessment when the twins were 9 years old (Peper et al., 2009). Two follow-up assessments were conducted when the twins were 12 and 17 years old (Koenis et al., 2017;Teeuw et al., 2018;van Soelen et al., 2012b. Here, we report results of the analysis of resting-state functional MRI scans that were acquired during the second and third assessment of the BrainSCALE study, when the twins and siblings were on average 13 and 18 years of age, hereafter referred to as time point 1 (TP1) and time point 2 (TP2). Intelligence was assessed using an abbreviated version of the Weschler Intelligence Scale for Children -Third edition (WISC-III; Wechsler, 1991) IQ test at age 13 years, and an abbreviated version of Weschler Adult Intelligence Scale -Third edition (WAIS-III; Wechsler, 1997) IQ test at age 18 years. The use of subtasks of the WISC-III as proxy for full WISC-III IQ test has previously been established as a valid construct (Koenis et al., 2015).
The BrainSCALE study was approved by the Central Committee on Research Involving Human Subjects of The Netherlands (CCMO), and studies were performed in accordance with the Declaration of Helsinki. Children and their parents signed informed consent forms. Parents were financially compensated for travel expenses, and children received a present or gift voucher at the end of the testing days. In addition, a summary of cognition scores and a printed image of their T1 brain MRI scan, when available, were provided afterwards.

MRI acquisition
Whole-brain magnetic resonance imaging (MRI) scans were acquired on two identical 1.5 T Philips Achieva scanners (Philips, Best, Netherlands) at the University Medical Center Utrecht (UMCU). Threedimensional T1-weighted structural scans (Spoiled Gradient Echo; TE ¼ 4.6 ms; TR ¼ 30 ms; flip angle ¼ 30 ; 160 to 180 contiguous coronal slices of 1.2 mm; in-plane resolution of 1.0 Â 1.0 mm 2 ; acquisition matrix of 256 Â 256 voxels; field-of-view of 256 mm with 70% scan percentage (Peper et al., 2009;Teeuw et al., 2018) and resting-state functional MRI scans (PRESTO-SENSE; TE ¼ 31.1 ms; TR ¼ 21.1 ms; flip angle ¼ 9 ; 4.0 mm isotropic voxels; 900 vol over 540 s; effective TR ¼ 600 ms) of the whole head were acquired (van den Heuvel et al., 2013). The same scanners and scan sequence parameters were used to acquire MRI scans at both ages to limit possible effects of differences in scan acquisition. Subjects were instructed to lie still with theirs eyes closed and keep their mind free from thoughts while preventing from falling asleep during acquisition of the resting-state functional MRI scans. Invited participants were excluded from the scanning protocol when contraindications for MRI were present at the time of examination. In particular, the presence of dental braces incompatible with the magnetic field of the MRI scanner resulted in a decline of participants for the neuroimaging assessment, specifically at age 13 years.

Processing of resting-state functional MRI scans
Processing of the MRI scans was performed using the CONN toolbox version 18a (Whitfield-Gabrieli and Nieto-Castanon, 2012; https://we b.conn-toolbox.org/) and SPM toolbox version 12 (http://www.fil .ion.ucl.ac.uk/spm/) in MATLAB version 2015b (The MathWorks Inc., Massachusetts, United States). The CONN toolbox is an open-source toolbox for processing and analysis of resting-state functional MRI scans. The toolbox is based on the aCompCor method for artefact correction that performs linear regression of undesired confounders, such as head motion and signal from white matter and cerebrospinal fluids, to recover the neuronal BOLD signal of interest (Behzadi et al., 2007). This artifact correction method has shown to reduce motion-related artifacts in resting-state fMRI in children (Muschelli et al., 2014).
To obtain the signal from white matter and cerebrospinal fluid (CSF), brain tissue from the structural T1-weighted MRI scans was segmented into CSF, gray matter (GM), and white matter (WM) tissue maps using a partial volume segmentation algorithm that incorporates a non-uniform partial volume distribution . The structural T1-weighted MRI scans were registered to MNI-152 space using non-linear transformation. The non-linear transformation was then applied to the tissue maps to warp them to MNI-152 space and resampled to 3.0 mm isotropic voxels. The white matter and CSF tissue maps were threshold at 50% (i.e. selecting only voxels with >50% of tissue proportion attributed to white matter or CSF) and binarized to create masks. The white matter tissue masks were eroded by two voxels to reduce the number of voxels at the white-gray matter tissue interface. No erosion was performed for the CSF tissue masks due to the occasional small volume of the lateral ventricles in children at age 13 years (Giedd et al., 1996;Lenroot and Giedd, 2006;Sowell et al., 2002). Instead, CSF tissue masks were constrained to contain only voxels inside the lateral ventricles.
The volumes within the resting-state functional MRI scans were first realigned to the mean image of the volumes using a rigid-body realignment procedure without reslicing the data. The rigid-body transformation parameters were used to retrospectively estimate head movement during scan acquisition using framewise displacement (Power et al., 2012). Mapping between resting-state functional MRI scans and structural MRI scans was determined by linear registration of the mean of the realigned resting-state functional MRI scan to the structural T1-weighted MRI scan. By concatenating all transformations (realignment, functional-to-T1 and T1-to-MNI), the mapping between individual functional space and MNI-152 space was obtained. The resulting transformation was used to warp the resting-state functional MRI scans into MNI space and resampled to 3.0 mm isotropic voxels. Global signal fluctuations time series were extracted from the warped functional MRI scans using the DVARS method (Power et al., 2012).
Correction of undesired confounders was performed using linear regression of the top ten principal components from the BOLD signal of white matter and (ventricular) cerebrospinal fluids (Behzadi et al., 2007;Chai et al., 2012), 24 head motion parameters (Friston et al., 1996;Yan et al., 2013), and scrubbing of subject-dependent number of high motion frames (Power et al., 2012). Linear regression was performed on the individual voxels of the brain after linear and quadratic detrending of the BOLD time series to reduce effects of scanner drift. The six rigid-body transformation parameters (R) derived during realignment of resting-state volumes, its first-order temporal derivative (R 0 ), and the squared product of all terms (R 2 and R 02 ) were included as regressors to control for head motion during scan acquisition (Friston et al., 1996;Yan et al., 2013). In addition, scrubbing of frames with high motion (FD > 0.30 mm) or unusually large whole-brain BOLD signal changes (DVARS Z-score > 3.0) was performed by including a regressor for each of the flagged frames, the frame immediately preceding the flagged frame, and the two frames following the flagged frame (Power et al., 2012); see supplementary information for more details of head motion and scrubbing. The average number of flagged frames is 79 (9%) out of 900 frames at age 13 years, and 57 (6%) out of 900 frames at age 18 years. By including frames surrounding the flagged frames the average number of scrubbed frames is 214 (24%) out of 900 at age 13 years, and 154 (17%) out of 900 frames at age 18 years. The residuals of the linear regression provided the voxel-wise denoised time series with a duration of 900 frames regardless of the number of frames scrubbing used in the regression. Temporal bandpass filtering was applied at the frequency range of 0.008-0.080 Hz after linear regression was performed that contained on average 39% of the total power spectral density after denoising (Biswal et al., 1995;Ciric et al., 2017;Waheed et al., 2016).
All resting-state functional MRI scans were processed independently from each other, including scans from subjects with longitudinal data.

Functional connectivity estimates
Functional connectivity matrices were obtained for the resting-state networks atlas provided by the CONN toolbox (Whitfield-Gabrieli and Nieto-Castanon, 2012; https://web.conn-toolbox.org/). The atlas is based on ICA analysis of resting-state scans of 497 unrelated young adults from the Human Connectome Project and provides regions of interest (ROI) for 7 canonical cortical resting-state networks and the cerebellum: the core Default Mode network (4 components), Sensorimotor (3), Visual (4), Salience (7), Dorsal Attention (4), Frontoparietal (4), Language (4), and Cerebellar (2); see supplementary information for details on the CONN atlas (Supplementary Fig. S1; Supplementary Table S1), and for comparison, a group-ICA decomposition was performed on the Brain-SCALE resting-state functional MRI scans used in this analysis (Supplementary Fig. S2).
The CONN toolbox atlas is based on a data-driven decomposition of resting-state functional data from the Human Connectome Project that consists of most large-scale canonical networks covering a large area of the cortical surface. Further decomposition of each network into separate regions that include homologous contra-lateral regions allows for studying global patterns of the developing functional brain within and between resting-state networks. Its limited number of regions makes it suitable for the computational complexity of twin modelling, and the moderate size of the regions provide the benefit of increased signal-tonoise ratio (SNR) through spatial averaging of the inherently noisy BOLD signal from neighboring voxels.
Full-score and half-score measures of functional connectivity were obtained using full Pearson correlation between spatially averaged denoised time series of two regions of interest (ROI). Full-score measures were obtained by considering the entire denoised time series of the 900volume resting-state functional MRI scan. Half-score measures were obtained by splitting the denoised time series into two independent halves of equal length; the first half-score measure (H1) corresponding to the first 450 volumes of the denoised time series, and the second half-score measure (H2) consisting of the remaining 450 volumes. All functional connectivity correlations were transformed using Fisher's r-to-Z transformation prior to any statistical analysis. Mean functional connectivity for subsets of connections (e.g. mean functional connectivity between all resting-state networks) was calculated as the average of the r-to-Ztransformed correlations across the subset of connections.

Quality control
Incomplete resting-state fMRI scans and scans with discernable defects related to scanner artefacts or receiver coil malfunction were discarded beforehand, resulting in the exclusion of 18 of the 380 (5%) scans from 17 subjects; 11 of the 152 (7%) scans at age 13 years, and 7 of 228 (3%) scans at age 18 years (Supplementary Table S2). In addition, individual full-score and half-score measures were excluded when the corresponding mean FD of the time series exceeded 0.30 mm/volume head motion or the number of scrubbed frames exceeded more than half the number of frames in the time series (i.e. more than 450 scrubbed frames for full-score measures, and more than 225 scrubbed frames for halfscore measures). After filtering, 97 full-score measures remain at TP1 (age 13 years), 108 half-score measures at TP1 H1, 88 half-score measures at TP1 H2, 202 full-score measures at TP2 (age 18 years), 203 halfscore measures at TP2 H1, and 200 half-score measures at TP2 H2 (Supplementary Table S2; Supplementary Fig. S4). See supplementary information for more details on quality control procedure and analysis of head motion.

Genetic modelling
Genetic modelling of twin and sibling data can provide information on the extent that variation of a trait in the population is explained by genetic factors (Boomsma et al., 2002;. Monozygotic (MZ) twins share 100% of their genetic material and dizygotic twins and full siblings share on average 50% of their segregating genes. Inclusion of these relatives into an extended twin design enables decomposition of the phenotypic variance (V P ) of a trait into three variance components: additive genetic (V A ), common environmental (V C ), and unique environmental (V E ) components of variance. Additive genetic influences represent effects of multiple alleles at different loci across the genome that act upon the phenotypic trait. Common environment represents influences that are shared between twins and siblings from the same family and causes them to be more alike than children grown up in different families. Unique environmental influences are not shared between family members and may include measurement error (Boomsma et al., 2002;Falconer and Mackay, 1996). If monozygotic (MZ) twins resemble each other more than dizygotic (DZ) twins and siblings for a trait, then the hypothesis that the trait is influenced by genetic factors is supported. If both MZ and DZ twins are more alike in resemblance than expected based on genetics alone, common environmental influences are likely to play a role.

Structural equation modelling
Within genetic structural equation modelling (SEM), a phenotype can be modelled as influenced by latent additive genetic factors, and common and unique environmental factors. These factors are modelled as unobserved or latent variables with unit variance where path coefficients going from latent trait to phenotype and symbolized by a, c, and e quantify their respective influence on the phenotype. The model is made identifiable by putting constraints on the correlation ρ A between the latent variable A of family members; ρ A ¼ 1:0 for monozygotic twins, and ρ A ¼ 0:5 for dizygotic twins and twin-sibling pairs. The correlation ρ C between latent variable C of family members is constrained to ρ C ¼ 1:0 for all twins and siblings from the same family. The latent variable E is uncorrelated between individuals. The sum of the squared path coefficients a 2 , c 2 , and e 2 , representing the variance components A, C, and E, is equal to the phenotypic variance (V) of a trait; i.e.
Heritability (h 2 ) of the trait is estimated as the proportion of phenotypic variance (V) that is due to additive genetic variance (A); i.e.
Structural equation models were defined using OpenMx version 2.8.3 (Neale et al., 2015; https://openmx.ssri .psu.edu/), a package for structural equation modelling in R version 3.4.2 (R Core Team, 2017; https://www.r-project.org/). Model fitting was performed using full-information maximum likelihood (FIML) to take advantage of all available information in case of missing data.

Modelling of twin and sibling data
Longitudinal data from extended twin designs can be modelled by Cholesky decomposition (Supplementary Fig. S3) to estimate the genetic and environmental influences on repeated observations (Neale and Cardon, 1992). A longitudinal Cholesky decomposition, with multiple measurements of the same trait acquired at different ages within the same individuals, allows for estimating the dynamics of genetic and environmental influences traits over age. This model provides estimates of heritability at the individual ages and of the genetic and environmental correlations that explain the sources of stability across ages. Genetic correlations represent the extent to which the same genes influence a trait at multiple ages. A longitudinal Cholesky decomposition can also provide indication of fluctuating influences of the same genes over time, or the presence of novel genetic influences (i.e. innovation) unique to a specific age (Teeuw et al., 2018;van Soelen et al., 2012b). Here, we modelled the split-half phenotypic information as a common factor at mean ages 13 and 18 years (van Baal et al., 1998;van Beijsterveldt et al., 2001) such that a latent phenotypic factor represented the reliable component of the two half-score measures at each age. Residual variance of a measurement (E s ) that is not attributed to a latent phenotypic factor is considered to be measurement error due to imprecisions of the measurement instrument. Estimation of genetic and environmental components was carried out for two submodels: the model with two latent phenotypic factors loading on the half-scores at each age, and a second model with a single latent phenotypic factor loading on all four half-scores (Fig. 1). Heritability of a latent phenotypic factor F j (h 2 j ) is estimated as the proportion of additive genetic variance of the latent phenotypic factor (A j ) over the variance of the latent phenotypic factor The heritability estimate of a latent phenotypic factor can be projected back to obtain heritability estimates for the individual half-score measures. According to path tracing rules, heritability of an individual measurement M k (h 2 k ) due to heritability of the latent phenotypic factor is the sum of the multiplication of the path coefficients along all the paths that visit an additive genetic variance component; e.g. the heritability of the first half-score at age 13 years (M 1 ) is h 2 1 ¼ f 11 Á a 11 Á a 11 Á f 11 , or simplified, h 2 1 ¼ f 2 11 Á a 2 11 , and the heritability of the first half-score at age 18 years (M 3 ) in the two-factor model is h 2 3 ¼ ðf 32 Á a 21 Á a 21 Á f 32 Þ þ ðf 32 Á a 22 Á a 22 Á f 32 Þ, or simplified, h 2 3 ¼ f 2 32 Áða 2 21 þa 2 22 Þ (see Fig. 1). The genetic correlation (r a ) between the latent phenotypic factors in the two-factor model is defined as the additive genetic covariance between the two factors over the square root product of the additive genetic variances of the two factors: r a ¼ covðA1;A2Þ ffiffiffiffiffiffiffiffiffiffi A1 Á A2 p .

Associations between functional connectivity and intelligence
The phenotypic associations between full-score functional connectivity (FC) estimates and intelligence (IQ) at each age, and the phenotypic associations between change in full-score functional connectivity with age (ΔFC) and change in IQ scores (ΔIQ) was assessed using a bivariate Cholesky model. Change scores were computed as the difference between the two ages. The bivariate Cholesky model included sex and age as fixed effects on each of the measures to account for possible mean differences between sexes and ages.

Statistical testing
Significance of parameters was determined using the log-likelihood ratio test by comparing the likelihood of the model with additional constraints on the parameters to the likelihood of the less constrained model. For bounded variance components (e.g. heritability estimates), the difference in À2 times the log likelihood (-2LL) between models with a single constraint follows a 50:50 mixture of χ 2 distributions with zero and one degree of freedom, and a 50:45:5 mixtures of χ 2 distributions with zero, one and two degrees of freedom for models with two constraints, etcetera; all effectively allow p-values to be cut in half (Dominicus et al., 2006).
Correction for multiple comparison was performed using FDR (Genovese et al., 2002) per condition (e.g. testing for effects of sex on functional connectivity) for all between and within resting-state network connections including mean functional connectivity estimates for a total of 92 tests per condition.

Model selection
The log-likelihood ratio test was used to determine the most parsimonious model amongst the models with different configuration of variance components on the latent phenotypic factors (i.e. ACE, AE, CE or E). The optimal number of latent phenotypic factors was determined using the log-likelihood ratio test on the models with ACE variance components.
In advance of the results, note that in 24 of the 92 cases (26%) a model with two factors was optimal (Supplementary Table S4; Supplementary Table S5) and that about half of these two-factor models do not have statistically significant heritability or common environmental influences at both ages (i.e. either E-AE, E-CE, or E-E configuration; Supplementary Table S6), and show mostly stable genetic influences primarily driven by the increased power at age 18 years. Although the presence of two factors might also be due to changes in unique environmental influence, we cannot distinguish between true unique environmental influences and age-specific measurement error in the common pathway reliability model with two factors. We therefore present our main analysis using the results of the common pathway reliability model with a single factor first, followed by the results from two-factor model in a post-hoc fashion.

Post-hoc analyses
We performed a post-hoc analysis to validate the main findings from the CONN atlas when global signal regression (GSR) is applied during the preprocessing stage. We performed a second post-hoc analysis to validate the main findings using Yeo resting-state networks atlas (Yeo et al., 2011). This atlas has a slightly different parcellation of the cortical surface into networks, which includes an extended default mode network, i.e. the parahippocampal and temporal regions in addition to the regions of the core default mode network.

Demographics
Data from 240 participants with either one or two resting-state functional MRI scans that passed quality control were included in the analysis, providing a total of 315 scans (Table 3). The twins were on average 12.2 AE 0.23 and 17.2 AE 0.17 (mean AE SD) years old at time point 1 (TP1) and time point 2 (TP2), with their older siblings on average 2.7 AE 1.2 (mean AE SD) years older.

Stability and reliability of functional connectivity
Group-level mean full-score functional connectivity matrices between and within resting-state networks reveal minor changes in functional connectivity estimates between the two timepoints (Fig. 2). Mean functional connectivity at group-level ranges from þ0.25 to þ0.71 at age 13 years and from þ0.17 to þ0.73 at age 18 years for connections between resting-state networks, and ranges from þ0.32 to þ0.76 at age 13 years Fig. 1. A common pathway reliability model with two (left) and one (right) latent phenotypic factor. Measurements are represented by rectangles for age 13 years (TP1) and age 18 years (TP2) for the first (H1) and second (H2) half-score measures. Latent variables are represented by circles. The variance of latent common factor F j represent the reliable component of the measurements and explains part of the variance of the measurement variables proportional to the square of the path coefficients on paths f kj . Latent factors representing additive genetics (A i ), common environment (C i ), and unique environment (E i ) together explain the variance of the common factor F j proportional to the square of their respective paths a ji , c ji , and e ji . Measurement-specific variances (i.e. residual variances not attributed to the common factor) are explained by the latent variables E Sl explaining the residual variances of measurements proportional to the square of the loadings on paths e Skl . Family members are linked through bidirectional paths on their latent variance components with values constrained to 1.0 for the additive genetic factor(s) (A i ) between monozygotic twins, 0.5 for the additive genetic factor(s) (A i ) between dizygotic twins and siblings, 1.0 for common environmental factors (C i ) for all pairs within one family, and 0.0 for unique environment (E i ). and from þ0.33 to þ0.84 at age 18 years for connections within restingstate networks ( Fig. 2 Table S9). The strongest increases in functional connectivity occur mostly within the sensorimotor network (β age mean FC within SMN ¼ þ0.0256; p < 0.001; FDR-corrected p < 0.001; Fig. 3; Supplementary Table S9) and the visual network (β age mean FC within VN ¼ þ0.0190; p < 0.001; FDRcorrected p < 0.001; Fig. 3; Supplementary Table S9). In addition, most contralateral connections between homotopic regions of the hemispheres are amongst the strongest to increase with age ( Fig. 3; Supplementary   Table S9). Several ipsilateral connections within the frontoparietal, salience, and language networks (e.g. the connection between lateral prefrontal cortex (LPFC) and posterior parietal cortex (PPC) show increase in functional connectivity ( Fig. 3; Supplementary Table S9). Including mean framewise displacement as an additional covariate in the functional connectivity analysis does not alter the results. Moreover, the same pattern of age-related effects is found in a post-hoc analysis when global signal regression was included during the preprocessing stage, and when repeating the analysis with Yeo's resting-state networks atlas.

Sex and functional connectivity
Sex effects on functional connectivity are sparse, and after multiple comparison correction found only for connections within the default mode network (DMN) and salience network (SN) (Fig. 3; Supplementary  Table S7; Supplementary Table S9). Increased functional connectivity for girls compared to boys is found within the default mode network (DMN) (β sex ¼ À0.0748; p < 0.001; FDR-corrected p ¼ 0.008; Fig. 3; Supplementary Table S9), and is mostly due to the connections between the medial prefrontal cortex (MPFC) and the PCC (β sex ¼ À0.0998; p < 0.001; FDR-corrected p ¼ 0.006; Fig. 3; Supplementary Table S9) and between the left lateral parietal cortex (LLP) and the PCC (β sex ¼ À0.0934; p ¼ 0.001; FDR-corrected p ¼ 0.025; Fig. 3; Supplementary Table S9). An opposing sex effect, where functional connectivity for boys is greater than for girls, is found for the connection between the left and right anterior insula (aIns) within the salience network (SN) (β sex ¼ 0.1018; p ¼ 0.001; FDR-corrected p ¼ 0.025; Fig. 3; Supplementary Table S9).

IQ and functional connectivity
None of the associations between functional connectivity and IQ test scores survived correction for multiple comparison (FDR-corrected p ! 0.3049 [n.s.]; Supplementary Table S11; Supplementary Table S12).
Within the core default mode network, connections involving the precuneus (PCC) show additive genetic influences (MPFC-PCC IQ scores and mean FD for full-scores measures of participants with longitudinal data were averaged across both measurements in the longitudinal data column. Flagged frames are the number of frames within the resting-state fMRI scan that exceeded the threshold for head motion of FD > 0.30. Scrubbed frames are the number of frames for which regressors were included in the preprocessing stage, and is derived from expansion of the flagged frames by also included one frame prior and the two frames following flagged frames (Power et al., 2012 Table S10). Mean functional connectivity averaged over all six connections within the default mode networks show significant influences of common environment (c 2 ¼ 37%; p ¼ 0.003; FDR-corrected p ¼ 0.031; Fig. 4; Supplementary Table S10).
By definition of the common factor model, heritability at the individual half-score measures depends on the proportion of variance explained by the common factor and the heritability of the reliable factor. Averaged standardized factor loading across the four half-scores of the individual connections ranged from 18% to 46% (mean standardized factor loadings 33%), with higher loadings on measurements at age 18 years (mean standardized factor loadings: 41%) than at age 13 years (mean standardized factor loadings: 24%) (Supplementary Fig. S7; Supplementary Fig. S8). Heritability estimates at individual half-score measurements ranged from 5% to 53% and from 5% to 33% for common environment estimates (Supplementary Fig. S7; Supplementary Fig. S8; Supplementary Table S8; Supplementary Table S10).
We find a similar pattern of additive genetic and common environmental influences on the reliable components of functional connectivity in a post-hoc analysis when using Yeo's resting-state networks atlas. However, when applying global signal regression during the preprocessing stage, common environment estimates on the reliable component of functional connectivity, which no longer contains the global signal, is drastically reduced, and the reliable component of functional connectivity for connections previously influenced by common environment is now primarily influenced by additive genetics instead.

Dynamic genetic and environmental influences on functional connectivity throughout adolescence
A two-factor common pathway reliability model was a better fit to the data than a single factor for 24 of the 92 connections (26%) (Supplementary Table S4; Supplementary Table S5). However, only half of these connections have statistically significant heritability or common environment estimates at both ages to warrant longitudinal investigation into dynamics of genetic and common environmental influences. For seven connections there is indication of possible dynamics in additive genetic or common environmental influences with age (Supplementary  Table S6). Two connections within the default mode network (DMN) show common environmental influences on changes in age-related functional connectivity, between the medial prefrontal cortex (MPFC) and right lateral parietal (RLP) with distinct genetic influences at each age due to innovation (c 2  and occipital regions within the visual network (VN) shows additive genetic influences on changes in functional connectivity due to increasing heritability originating from the same set of genes with age (a 2 ðΔFCÞ ¼ 29% [2%; 96%]; p ¼ 0.015; fluctuating influences p ¼ 0.028; Supplementary Table S6). The connection between occipital and left lateral regions of the visual network (VN) shows common environmental influences on change in functional connectivity with innovation in environment factors over time (c 2 ðΔFCÞ ¼ 41% [17%; 93%]; p ¼ 0.012; innovation p ¼ 0.019; Supplementary Table S6). Three connections within the frontoparietal network (FPN) show with additive genetic influences on changes in functional connectivity due to fluctuating influences from the same genes: between the left lateral prefrontal cortex (LPFC) and left posterior parietal cortex (PPC) (a 2 ðΔFCÞ ¼ 58% [17%; 91%]; p ¼ 0.002; fluctuating influences p ¼ 0.003; Supplementary  Table S6). The remaining connections with genetic of common environmental influences at both ages do not reveal any significant dynamics in heritability or common environment (Supplementary Table S6). However, these results should be interpreted with caution due to limited power to detect significant genetic or environmental estimates at age 13 years, in part due to the reduced sample size at age 13 years.

Discussion
With this longitudinal resting-state fMRI study, we measured the heritability of functional connectivity throughout adolescent development for the first time. Approximately half of the functional connections within and between canonical cortical resting-state networks are influenced by either additive genetic (h 2 up to 53%) or common environmental influences (c 2 up to 33%) during adolescence. During adolescence, functional connectivity between resting-state networks decreases with age, whereas functional connectivity within cortical restingstate networks increases with age, except for several connections within the salience network that decrease with age. There is limited evidence for dynamics in genetic or common environmental influences, suggesting mostly stable influences across adolescence. Girls had significantly stronger functional connectivity than boys within the default mode network between the precuneus and medial prefrontal cortex and between the precuneus and left lateral parietal cortex. Boys had significantly stronger functional connectivity than girls within the salience network for the connection between the bilateral insula. Associations between functional connectivity and intelligence did not survive multiple comparison correction. Head motion is heritable across the ages and shows a small but statistically significant decline with age (Supplementary Table S3; Supplementary Fig. S5). The aCompCor method used by CONN is effective at removing head motion effect cross-sectionally, however, longitudinal changes in functional connectivity estimates between the canonical resting-state networks remain associated with the longitudinal changes in degree of head motion of individuals (Supplementary Table S13; Supplementary Table S14). The results remained consistent after including mean framewise displacement as an additional covariate in the functional connectivity analysis and after including global signal regression during the preprocessing stage.
We find significant heritability on functional connectivity in adolescence, h 2 ranging from 6% to 53% for 23 out of 55 (42%) connections within resting-state networks, and common environment estimates c 2 ranging from 5% to 33% for 8 out of 55 (15%) connections. Previous studies found heritability estimates ranging from 10% to 40%, and occasionally up to 60% or 80% (Adhikari et al., 2018a;Colclough et al., 2017;Fu et al., 2015;Ge et al., 2017;Glahn et al., 2010;Korgaonkar et al., 2014;Meda et al., 2014;Sudre et al., 2017;Yang et al., 2016, Table 2), thus overall these findings are within the same range across the ages. The notable exception is the default mode network. In our cohort the default mode network is partially influenced by common environmental instead of additive genetics (mean FC within DMN c 2 ¼ 37%). Previous studies have established the default mode network being influenced by additive genetics (Adhikari et al., 2018a;Fu et al., 2015;Ge et al., 2017;Glahn et al., 2010;Korgaonkar et al., 2014;Meda et al., 2014;Sudre et al., 2017;Xu et al., 2016;Yang et al., 2016). The discrepancy with previous studies may be due to increased sensitivity in finding common environmental effects in the current study because of three reasons. One, the extended twin design (i.e. including twin pairs and one of their singleton siblings) used in this study provides additional power to detect significant common environment estimates  allowing detection of statistically significant common environmental estimate at individual measures as low as 5% when separating measurement error from the common factor. However, most likely functional connections in the brain are influenced by both additive genetics and common environment. Two, since previous studies were conducted mostly in adults, it may be possible that heritability estimates on functional connectivity increase with age, as is the case with heritability of cognitive performance (Briley and Tucker-Drob, 2013). Indeed, the cohorts closest to our age range find similarly low estimates for heritability of functional connectivity within the default mode network (Fu et al., 2015;Yang et al., 2016). Finally, atlas choice and preprocessing strategies, including global mean signal regression, all varied between these studies, which could have introduced differences in heritability estimates. Although we were able to replicate our results using Yeo's resting-state atlas (Yeo et al., 2011), including global signal regression during the preprocessing stage substantially decreased the common environment estimates in favor of additive genetics for some connections.
Few other studies have investigated heritability of functional connectivity with resting-state networks beyond the default mode network (Adhikari et al., 2018a;Ge et al., 2017;Sudre et al., 2017;Yang et al., 2016). Our heritability estimates for functional connectivity within the frontoparietal network of h 2 ffi 14% at around age 13 years and h 2 ffi 40% at around age 18 years are slightly lower than the estimates of h 2 ¼ 32%-58% found across lifespan in families with ADHD family members (Sudre et al., 2017), and substantially less than the estimate of h 2 ¼ 65% reported in a sample of healthy young adults (Yang et al., 2016), but more in line with results from the Genetics of Brain Structure (GOBS) and the Human Connectome Project (HCP) studies (Adhikari et al., 2018a). The sensorimotor network in our cohort is influenced by a mixture of additive genetics h 2 ¼ 18%-20% and common environment c 2 ¼ 5%-16%, with influences of common environment previously reported in young adults (Yang et al., 2016), and low to no heritability in the GOBS and HCP cohort, although they did not test for common environmental influences (Adhikari et al., 2018a). This is in stark contrast to the study performed on the Brain Genomics Superstruct Project (GSP) cohort and alternative analysis of the HCP cohort where they found high heritability estimates of h 2 ¼ 60%-70% (Ge et al., 2017).
We find significant heritability of functional connectivity between resting-state networks in adolescence, with 8 out of 28 (29%) connections influenced by additive genetics with h 2 ¼ 5%-50%, and 6 out of 28 (21%) connections influenced by common environment with c 2 ¼ 6%-25%. In particular, connections between the frontoparietal, dorsal attention, and salience networks, all involved in higher order cognitive control, were influenced primarily by additive genetics. Common environment played a considerable role for most sensory networks, including the language network and cerebellum. So far, the only other study that investigated heritability of functional connectivity between resting-state networks was performed in young adults aged 18-29 years (Yang et al., 2016), where 8 out of 21 (38%) connections between resting-state networks were influenced by additive genetics with h 2 ¼ 26%-42%, and 11 out of 21 (52%) were influenced by common environment with c 2 ¼ 18%-47%, showing overall similar connectivity profiles. Interestingly, synchronous resting-state activity in the brain has been associated with gene expression levels, where distal functionally connected regions show similar gene expression profiles , and the strength in functional connectivity was influenced by polymorphisms of a set of genes enriched for ion channels in healthy adolescents (Richiardi et al., 2015). Here, our study adds that functional connectivity within and between cortical resting-state networks is strongly influenced by genes and common environment throughout adolescent development. Including global signal regression during the preprocessing stage resulted in decreased estimates for common environment and increased heritability estimates on the reliable component of functional connectivity for some of the connections. This insight will have potential consequences for genetic studies that aim to find genetic variants implicated in functional connectivity.
Varying estimates of heritability or common environment between half-score measures of different sessions (i.e. between the two ages) may be an indication of dynamics in genes or common environment that was tested with a two-factor common pathway model; e.g. the change in heritability estimates for functional connectivity within the frontoparietal network increased significantly from h 2 ffi 14% at age 13 years to h 2 ffi 40% at age 18 years. However, the lower heritability estimates at age 13 years are likely due to the reduced sample size as a result of motion scrubbing and exclusion due to presence of dental braces incompatible with high magnetic fields or increased residual noise rather than represent "true" changes in additive genetic or common environmental variances. This effect is reflected in the two-factor common pathway model as "increasing" influences of the same additive genetic or common environmental factor over time, and is consequently found in the single factor common pathway model as varying estimates of heritability due to differences in factor loadings on the individual half-score measures. Therefore, the results from the two-factor model on dynamics of genetic and environmental influences are suggestive at best and generally demonstrate stable additive genetic or common environmental influences from a single source (Supplementary Table S6). Varying estimates for half-score of a single session (i.e. within the same age) are very unlikely to represent short-term fluctuating genetic or environmental influences, but can most likely be attributed to fluctuating noise (e.g. slight increase in head motion or restlessness during second half of scan).
The longitudinal age effects that we found are subtle but wide-spread throughout the brain despite most resting-state networks already appearing "adult-like" by age 2 years (Gao et al., 2015;Gilmore et al., 2018). We found age-related decreases in functional connectivity for about half of the connections between cortical resting-state networks, which likely reflect segregation between functionally distinct modules of the brain. A previous longitudinal study during early adolescence reported segregation between the frontoparietal (FPN) and default mode network (DMN) (Sherman et al., 2014). Although our results are not statistically significant, it suggests a possible decrease between the FPN and DMN (β age ¼ À0.0076; p ¼ 0.061 [n.s.]; FDR-corrected p ¼ 0.089 [n.s.]) that is consistent with prior reports. Previously, a decrease in functional connectivity between the dorsolateral prefrontal cortex and posterior cerebellum, but not anterior cerebellum, was reported (Bernard et al., 2016). We do not find age-related changes in functional connectivity between the frontoparietal and cerebellar networks. However, we do not distinguish between sub-regions of the network for between network connectivity. Nearly two-thirds of the connections within resting-state networks show age-related increase in functional connectivity that likely reflect integration within functional modules of the brain. Previously longitudinal studies have reported on integration within the default mode (DMN), frontoparietal (FPN), and language network (LN) during childhood and adolescence (Long et al., 2017;Sherman et al., 2014;Wendelken et al., 2017Wendelken et al., , 2016Xiao et al., 2016). However, age-related changes in functional connectivity within the DMN, FPN, and salience network (SN) are not always consistently found during early adolescence (Sylvester et al., 2017). Similar to previous reports, we find integration within the frontoparietal network (FPN) for the ipsilateral connections between the frontal and posterior regions of the FPN (Wendelken et al., 2017(Wendelken et al., , 2016, and integration within the language network between the inferior frontal gyrus (IFG) and posterior superior temporal gyrus (pSTG) (Xiao et al., 2016). Better integration between frontal and parietal regions has been proposed to support better cognitive performance in the Parieto-Frontal Integration Theory (P-FIT; Deary et al., 2010;Jung and Haier, 2007). A notable exception to integration within resting-state networks is the age-related decrease in functional connectivity within the salience network for connections involving the anterior cingulate cortex (ACC), left rostral prefrontal cortex (RPFC) and left supramarginal gyrus (SMG). The only other longitudinal study that investigated the salience network in children and adolescents was between 8 and 13 years, preceding our age range, reporting absence of significant age-related effects (Sylvester et al., 2017). The anterior cingulate cortex plays an important role in motor control and cognition, in particular reward-based decision making and response inhibition (Bush et al., 2002;Stevens et al., 2011). A decreasing connectivity between the anterior cingulate cortex and insula could possibly reflect a decoupling between the integration of external sensory information and internal emotional and bodily state signals (Uddin et al., 2017) or indicate segregation of bottom-up stimuli processing and top-down cognitive control processing in the salience network that may coincide with improved self-control during adolescence (Casey, 2013). The decreasing connectivity strength between networks and increasing connectivity strength within networks that we find in this longitudinal study is the opposite pattern of what is typically reported in cross-sectional literature, where there is a shift from a local oriented (i.e. stronger functional connectivity between proximal regions) to a more distributed organization (i.e. stronger functional connectivity between distal regions) during childhood and adolescence (Cao et al., 2016;Ernst et al., 2015;Grayson and Fair, 2017;Stevens, 2016). This discrepancy could be due to residual head motion effects on changes in functional connectivity with age that are not always properly accounted for in studies predating 2012 (Power et al., 2012), although despite our stringent control for head motion there are still residual effects present. Secondly, cross-sectional studies do not always show consensus on the direction of change and affected regions (Stevens, 2016), which may be due to the cohort effect. Few longitudinal studies have been conducted to date (see Table 1), with even fewer conducting brain-wide analysis. Several longitudinal studies show increasing functional connectivity within functional networks or decreasing functional connectivity between functional networks with age consistent with our results (Bernard et al., 2016;Long et al., 2017;Sherman et al., 2014;Wendelken et al., 2017Wendelken et al., , 2016. We find significant sex effects within the default mode (girls showing stronger functional connectivity) and salience network (boys showing stronger functional connectivity). Sex effects in functional connectivity analyses are typically discarded as covariate of no interest, despite extensive support for sex effects in behavior (Gur et al., 2012;Gur and Gur, 2016), brain gray matter (Giedd et al., 2012;Ruigrok et al., 2014) and white matter structure (Herting et al., 2012;Ingalhalikar et al., 2014), and function (Sacher et al., 2013;Stevens and Hamann, 2012). A few studies have reported sex effects for functional connectivity within the default mode network (stronger functional connectivity in females compared to males) and salience network (stronger functional connectivity in males compared to females) in adults (Biswal et al., 2010) and across the lifespan . These previous reports are consistent with our findings, and corroborate that sex differences in brain functioning are already present during adolescence (Gur and Gur, 2016), although sex effects are not always found in these networks during development (Sol e-Padull es et al., 2016;Sylvester et al., 2017). The default mode network plays an important role in auto-biography memory and emotion regulation (Mak et al., 2017;Raichle, 2015;Zhou et al., 2018). The increased functional connectivity within the default mode network for girls may therefore explain their better performance at memory and emotive tasks (Gur et al., 2012;Gur and Gur, 2016). In contrast, the salience network plays an important role in overt attention/stimuli processing, integration of multimodal sensory information, and switching between passive and active tasks (Marek et al., 2015;Seeley et al., 2007;Uddin, 2014;Zhou et al., 2018). The increased functional connectivity within the salience network for boys may therefore explain their better performance at visuospatial and motor tasks (Gur et al., 2012;Gur and Gur, 2016).
We do not find any association between functional connectivity and IQ scores that remain significant after multiple comparison correction. One of the leading theories on cognitive development postulates that better integration between frontal and parietal areas supports better cognitive performance (Deary et al., 2010;Jung and Haier, 2007). Previous studies on resting-state functional connectivity have provided support for this theory in both children and adults (Dubois et al., 2018;Langeslag et al., 2013;Santarnecchi et al., 2017;Song et al., 2008;Vakhtin et al., 2014). However, despite better integration in the frontoparietal network (FPN) with age, we found no association with IQ test scores that survived multiple comparison correction (Supplementary Table S11; Supplementary Table S12).
Head motion is a major point of concern in resting-state fMRI measurements when studying functional connectivity during childhood and adolescent development (Power et al., 2012;van Dijk et al., 2012;Satterthwaite et al., 2013). Heritability of head motion during childhood, adolescence, and young adulthood was previously established (Achterberg et al., 2018;Couvy-Duchesne et al., 2014, 2016Engelhardt et al., 2017;Zhou et al., 2016), and is consistent with our findings during adolescence (h 2 ¼ 86% at age 13 years; p < 0.0001; and h 2 ¼ 43% at age 18 years; p ¼ 0.0003). Preprocessing the resting-state data with linear regression of white matter and cerebrospinal fluid signal components and derivatives of re-alignment parameters using CONN toolbox appears to be mostly effective to control for head motion at individual ages (i.e. cross-sectionally; Supplementary Table S13; Supplementary Table S14). However, in a longitudinal setting, change in the degree of head motion was still associated with changes in functional connectivity for more than half of the generally long-distance connections between resting-state networks, all with positive association (i.e. reduction in head motion results in smaller change in functional connectivity; Supplementary  Table S13; Supplementary Table S14), with only the association between head motion and the connection between the sensorimotor and salience network surviving multiple comparison correction. Including global signal regression during the preprocessing stage -despite its controversy to introduce artificial negative correlations considered to be most effective at reducing head motion effects -did not substantially change the results in addition to the aCompCor method, as previously found (Ciric et al., 2017). This suggests that despite stringent control of head motion during preprocessing of the resting-state fMRI scans, not all variance due to head motion is accounted for, possible due to complex non-linear interaction of head motion with the BOLD signal (Satterthwaite et al., 2013). However, most connections within resting-state network, including long-distance connections between frontal and parietal regions, and cross-hemisphere connections, showed few associations with head motion. Moreover, age-related changes in functional connectivity were unaffected whether or not head motion was included as covariate. All scans, including longitudinal scans, were processed independently due to lack of publicly available tools that support longitudinal pipelines such as used for structural imaging (Reuter et al., 2012), although longitudinal methods for preprocessing of resting-state fMRI are becoming available (Hart et al., 2018).
A number of other limitations should be taken into consideration when interpreting these results. First, resting-state fMRI scans were acquired at 1.5T MRI scanner, using an at the time state-of-the-art fast repetition time (TR) T2-weighted PRESTO-SENSE acquisition protocol. The field strength was intentionally not upgraded to higher field strengths to minimize effects of scanner differences for longitudinal data acquisition. Although the fast TR aims to minimize head motion between acquired volumes, the 3D acquisition makes it more sensitive to abrupt motion resulting in ghosting and blurring of the BOLD signal that may require a different approach to denoising than traditional 2D single-shot EPI acquisition protocol (van Gelderen et al., 2012). In addition, effort was taken to reduce in-scanner head motion using mock-scanner session before actual acquisition to acclimate the children during previous visitation at age 9 years when the children also experienced an MRI scan but no resting-state fMRI scan was acquired (Durston et al., 2009;van Soelen et al., 2012a). Secondly, although sample size is large for neuroimaging standard (i.e. N ¼ 108 at age 13 years and N ¼ 207 at age 18 years), it is modest in size for twin studies. In particular, the reduced sample size at age 13 years due to exclusion of high motion subjects has limited the number of subjects with longitudinal scans and will have an impact on the power of this study to detect age-related effects. Combining the data across the ages using common factor analysis with a single factor had the benefit to mitigate this reduction in power for detecting statistically significant heritability or common environment estimates under the assumption that these estimates remain stable across ages, which is the case for most connections. Thirdly, resting-state fMRI is inherently noisy with no optimal strategy to remove non-neuronal signal (e.g. due to head motion or vascular system) (Ciric et al., 2017;Parkes et al., 2018). Alternative approaches or more advanced methods, such as improving temporal signal-to-noise ratio using machine learning (Adhikari et al., 2018a) or the use of a longitudinal preprocessing pipeline (Hart et al., 2018) may proof beneficial for future studies. Finally, choice of atlas may influence the results as both the CONN and Yeo atlas consist of predefined regions that do not allow for individual variation in localization of brain function and the moderate size of the regions prevent analysis of possible localized effects with behavioral measures such as cognitive performance on IQ tests. In addition, both atlases do not include subcortical structures for which both longitudinal development and heritability has been reported. Unfortunately, we could not use more fine-grained atlases due to computational complexity of twin modelling and the limited spatial resolution of our scans. However, the comparable results for both the CONN toolbox atlas and Yeo's resting-state network atlas suggests at least reasonable robustness of findings independent of atlas choice.
In conclusion, there are wide-spread influences of additive genetics and common environment on the functioning of cortical resting-state networks during adolescent development that generally remain stable over time. Wide-spread subtle age-related changes in functional connectivity occur in the presence of sizable individual variation, with presence of sex effects to be taken into consideration in developmental studies on functional connectivity.

Data and source code availability
Requests for access to the data and code used in this analysis should be directed to the corresponding author. Definitions of the structural equation models in OpenMx are provided as supplementary files.