Structural neurodevelopment at the individual level - a life-course investigation using ABCD, IMAGEN and UK Biobank data

Developmental trajectories of whole-brain GMV during adolescence define three clusters

We began by estimating the longitudinal trajectories of GMV in 44 brain regions of interest (ROIs) (34 cortical and 10 subcortical ROIs) that spanned the whole brain of each adolescent across baseline (at age 14y) and two follow-up scans (at age 19y and 23y) in the IMAGEN study, adjusting for intracranial volume (ICV), sex, handedness and site (Methods). Individuals showed strong heterogeneity and clustering patterns in terms of baseline total GMV and GMV developmental trajectories (Supplementary Fig. 1). Next, we reasoned that neurobiologically meaningful clusters could be explained by the developmental patterns in a subset of ROIs. Therefore, we conducted dimension reduction via principal component analysis (PCA) and selected the first 15 principal components (PCs), which explained 80% of the total variation in whole-brain GMV trajectories, in the clustering analysis (Supplementary Table 1). The first and second PCs defined two combinations of GMV trajectories over the entire brain that were significantly associated with baseline total GMV (Supplementary Fig. 2a). However, they exhibited different association patterns with items of the Cambridge Gambling Task, where PC1 was significantly associated with delay aversion (r = 0.07, P_adj = 0.30) and risk adjustment (r = - 0.08, P_adj = 0.020), and PC2 was significantly associated with deliberation time (r = 0.1, P_adj = 0.003), overall betting (proportion bet) (r = 0.07, P_adj = 0.014) and risk-taking (r = 0.08, P_adj = 0.008) (Supplementary Fig. 2b). These PCs were then used in the multivariate clustering to identify groups of adolescents with distinct neurodevelopmental patterns.

Among 1,543 adolescents with at least two sMRI scans, our analyses identified three clusters of structural neurodevelopment (𝑃_permutation < 0.001) (Supplementary Fig. 3). Group 1 consisted of 711 (46.1%) adolescents, had high baseline total GMV and continuously decreasing GMV at follow-ups, which was consistent with the population GMV developmental trend²⁸. Group 2 included 765 (49.6%) adolescents and compared to Group 1, they had lower baseline total GMV, lower peak GMV, and slower rate of GMV decrease. In addition, adolescents in Group 2 are more likely to be older (Diff = 0.11y, P < 0.001) and be males (Diff = 9.5%, P < 0.001), have lower parental education (P = 0.020 for maternal education; P=0.003 for paternal education) and lower WISCIV full score at age 14 (Diff = -1.76, P = 0.011). The remaining 67 (4.3%) belonged to Group 3, among whom we observed the lowest baseline total GMV and surprisingly increasing GMV at follow-ups, which was opposite to the population developmental trend (Fig. 1b, Supplementary Fig. 4-6). Compared to Group 1, Adolescents in Group 3 are more likely to have lower parental education (P = 0.015 for maternal education; P = 0.010 for paternal education) and lower WISCIV full score at age 14 (Diff = -9.22, P < 0.001). The full demographic and baseline characteristics for each group were provided in Supplementary Table 2. Since we aim to investigate group-specific brain developmental patterns in late childhood, we further estimated the age-specific GMV growth rate in each ROI from age 5y to 25y in each group (Methods) using population neurodevelopmental curve as a reference. Consistently we observed continuously decreasing GMV in Group 1 and Group 2 (with slower rate of GMV decrease in Group 2), and increasing GMV in Group 3 for most ROIs (Fig. 1c), indicating delayed neurodevelopment and brain maturation in Group 3 compared to the other groups.

To understand the neurobiological basis of group heterogeneity, we next tested for differences in whole-brain GMV development among these groups. We observed common delayed GMV development in ROIs spanning the inferior temporal, middle temporal, lateral orbitofrontal, precentral and superior frontal areas in Group 3 (relative to Groups 1/2) (Fig. 1d top and Supplementary Table 3). Group 2 showed lower peak GMV and slower rate of GMV decrease in ROIs spanning superior frontal, caudal middle frontal, rostral middle frontal, precentral and inferior parietal areas (relative to Group 1) (Fig. 1d bottom and Supplementary Table 3). These are all among the last areas in the brain to mature and had been implicated to play a key role in executive functions. This led us to ask whether variations in structural neurodevelopment could predict the developmental trajectories of neurocognition and risk of neuropsychiatric disorders in these groups.

Structural neurodevelopment predicts neurocognition and risk factors for neuropsychiatric disorders

To investigate the association between neurodevelopment and executive functions, we tested for differences of neurocognitive performance among these groups at baseline and at the last follow-up. Full results of these comparisons were provided in Supplementary Table 4. We found that compared to Group 1, Group 3 with delayed neurodevelopment showed worse neurocognitive performance (Spatial Working Memory, Cambridge Gambling Task (CGT) and Stop Signal Task (SST)) at baseline, but most of these items improved over time with brain maturation and became statistically equivalent at the last follow-up (Fig. 2a and Supplementary Fig. 7a). This can be predicted by the structural architecture of GMV development in Group 3, where increasing GMV in the top discriminating ROIs showed positive correlation with improvements of neurocognition (Supplementary Fig. 8 and Supplementary Table 5). In contrast, Group 2 with slower rate of GMV decrease showed worsened neurocognitive performance (CGT and SST) at the last follow-up relative to baseline (Fig. 2a and Supplementary Fig. 7b), which could be predicted by the negative correlations between the GMV developmental trajectories in the top discriminating ROIs and neurocognition (Supplementary Fig. 8 and Supplementary Table 6).

• Download figure
• Open in new tab

Fig. 2 Structural neurodevelopment predicts neurocognition and risk factors for neuropsychiatric disorders.

(a) Comparison of neurocognitive performances between Group 3 and Group 1 (orange), and between Group 2 and Group 1 (green) at baseline (BL) and the last follow-up (FU) (Methods). Total number of neurocognitive tests in cantab where Group 3 performed worse than Group 1 decreased from 7/12 at baseline to 1/12 at follow-up 3, while the number of tests where Group 2 performed worse than Group 1 increased from 0 to 3/12. Full results with item-specific comparisons among these groups are provided in Supplementary Fig. 7. CANTAB, Cambridge Neuropsychological Test Automated Battery. (b) Longitudinal trajectories of Depression (Left) and ADHD symptoms (Right) among adolescents in the three groups. Group-specific means at each visit were plotted and * indicated significant differences relative to Group 1 adjusting for sex, handedness, stie and ICV. Baseline mental health score was also adjusted for comparison at the last follow-up. Statistical tests were conducted at baseline (14y) and the last follow-up. Depression, Group 2 vs Group 1 at 14y (d=-0.05, P=0.256), Group 2 vs Group 1 at 23y (d=0.13, *P=0.023), Group 3 vs Group 1 at 14 y (d=-0.01, P=0.566), Group 3 vs Group 1 at 23y (d=0.70, **P=0.001); Parent rated ADHD, Group 2 vs Group 1 at 14y (d=0.04, P=0.220), Group 2 vs Group 1 at 16y (d=-0.03, P=0.574), Group 3 vs Group 1 at 14 y (d=0.34, **P=0.004), Group 3 vs Group 1 at 16y (d=0.01, P=0.954); Child rated ADHD, Group 2 vs Group 1 at 14y (d=-0.03, P=0.321), Group 2 vs Group 1 at 23y (d=0.02, P=0.758), Group 3 vs Group 1 at 14 y (d=0.34, P=*0.042), Group 3 vs Group 1 at 23y (d=-0.10, P=0.579). ADHD, attention-deficit/hyperactivity disorder.

The delayed brain and neurocognitive development in Group 3 led us to ask whether these adolescents had increased risk for neuropsychiatric disorders. Consistent with the improvements of neurocognition, we observed decreased attention-deficit/hyperactivity disorder (ADHD) symptoms, but increased depression symptoms in Group 3 (Fig. 2b and Supplementary Table 4). This indicated that although neurocognitive abilities in Group 3 exhibited pronounced improvement during adolescence, this was not necessarily the case for mental disorder symptoms. Furthermore, consistent with their worsened neurocognitive performances, we observed increased depression symptoms in Group 2 at the last follow-up compared to baseline, (Fig. 2b and Supplementary Table 4). The continuously worsened neurocognition and mental health problems in Group 2 indicated biological, social and mental disadvantages among these adolescents.

Given the slightly different patterns of GMV development for males and females²⁸, we conducted the analyses stratified by sex following the same workflow. Results of group clustering largely overlapped with the original analyses (Supplementary Table 7). In general, the sex-stratified analyses revealed similar patterns of neurocognition and mental health symptoms among three groups of adolescents. However, differences of neurocognition among these groups were manifested more for risk-taking and impulsive behaviors in males, while for spatial working memory in females (Supplementary Table 8). Besides, increase of the depressive symptoms in Group 2 was only observed in males, and increase of the depressive symptoms in Group 3 was only observed in females.

In addition, we compared the genetic liability to major neurodevelopmental disorders and related traits, including ADHD, autism spectrum disorder (ASD), educational attainment (EA) and intelligence (IQ), by calculating the corresponding polygenic risk scores (PRS) for each adolescent. Group 3 had higher PRS for ADHD than both Group 1 (P_adj = 0.007) and Group 2 (P_adj = 0.017), while Group 2 was not statistically different from Group 1 (P_adj = 0.42). We did not observe significant differences among the three groups in terms of the PRS of ASD, EA and IQ (Supplementary Table 9). The higher genetic liability of ADHD in Group 3 led us to ask whether genetic variants could explain the delayed neurodevelopment and neurocognitive performances in this group.

Genetic variation and epigenetic changes contribute to structural neurodevelopment

To better understand the genetic basis of structural neurodevelopment, we conducted genome-wide association studies (GWAS) for Group 3 versus Groups 1/2 using group-reweighted GMV as the proxy-phenotype among 7,662 adolescents in ABCD, since GWAS was under-powered for the IMAGEN study due to limited sample size³⁹. Group-reweighted GMV was derived and used as the proxy phenotype because GMV developmental patterns could not be estimated in ABCD due to limited age range. This continuous phenotype represented one’s tendency of being in Group 3 relative to Groups 1/2, or in other words, one’s propensity of having delayed brain development. Specifically, we began by calculating the ROI-specific weight in discriminating Group 3 (relative to Groups 1/2) in IMAGEN adjusting for potential confounders, and applying these weights to corresponding ROIs in ABCD to obtain the Group3-reweighted GMV, which was then used as the proxy phenotype in the Group 3 GWAS (Methods). The Group3-reweighted GMV showed negative correlation with neurocognition in ABCD, indicating the validity of using Group3-reweighted GMV as appropriate proxy for delayed neurodevelopment (Fig. 3a and Supplementary Table 10). Similarly, Group2-reweighted GMV was calculated and used as the proxy phenotype in the Group 2 GWAS.

• Download figure
• Open in new tab

Fig. 3 Genome-wide association study (GWAS) identified two significant loci associated with delayed neurodevelopment in Group 3.

(a) Correlation between Group3-reweighted GMV and neurocognition (Supplementary Methods) in ABCD (n=11,101) indicated the validity of using the proxy phenotype for delayed neurodevelopment in the following GWAS. GDT, Game of Dice Task; DDT, Dealy Discounting Task; PVT, Picture Vocabulary Test; Flanker, Flanker Inhibitory Control and Attention Test; List, List Sorting Working Memory Test; DCCS, Dimensional Change Card Sort Test; Pattern, Pattern Comparison Processing Speed Test; PSMT, Picture Sequence Memory Test; Reading, Oral Reading Recognition Test; FluidCog, fluid cognition; CrystalCog, crystallized cognition; TotalCog, total cognition. *** P<0.001; ** P<0.01; * P<0.05. (b) GWAS Manhattan plot for Group3-reweighted GMV in the ABCD population (n=7,662) . Group3-reweighted GMV was calculated for each adolescent (Methods) and used as the proxy phenotype for delayed neurodevelopment. Multiple SNPs on chromosome 6 achieved genome-wide significant effects (P<5×10^-8). SNPs on chromo-some 6 were mapped to the intronic region of CENPW. Results from gene based association analysis (Supplementary Fig. 11) confirmed the significant effect of CENPW on delayed neurodeveloopment. Box plot in (c) showed that CENPW score of delayed neurodevelopment was higher in Group 3 (n=60) compared to Group 1 and 2 (n=1,338) (two-sided t-test: P=0.028). (d) indicated that CENPW score of delayed neurodevelopment was negatively correlated with baseline (BL) neurocognitive performance, and became non-significant at the last follow-up (FU). Here, Worse indicated higher CGT Delay aversion score, lower CGT risk adjustment score, longer CGT Delibration time and SST GoRT. CGT Delay aversion, BL (r=0.09, *P_adj=0.027), FU3 (r=0.07, P_adj=0.239); CGT Delibration time, BL (r=0.08, *P_adj=0.027), FU3 (r=-0.02, P_adj=0.983); CGT risk adjustment, BL (r=-0.08, *P_adj=0.027), FU3 (r=0.022 P_adj=0.983); SST GoRT, BL (r=-0.06, *P_adj=0.038), FU3 (r=-0.03, P_adj=0.472). CGT, Cambridge Gambling Task; SST GoRT, reaction time for ‘Go’ trials in Stop Signal Task. (c-d) confirmed the relationship between CENPW and delayed neurodevelopment identified in (b).

One locus showed genome-wide significant effects in the Group 3 GWAS (Fig. 3b and Supplementary Table 11). The lead single-nucleotide polymorphisms (SNP), rs9375442 (β = 0.52, P = 1.02×10^-8) on chromosome 6, is an intronic variant located on CENPW (Supplementary Fig. 9). CENPW is a protein-coding gene involved in the packaging of telomere ends and cell cycle mitotic^{40, 41}, and increased CENPW expression in progenitors could lead to decreased cortical volume and cognitive function by altering neurogenesis or increasing apoptosis⁴². Other variants on these genes were reported to be associated with cortical surface area and brain volume^43–48 (Supplementary Fig. 10), general cognitive ability^49–51 and physical growth^52–56. Gene-based association analysis confirmed the identification of CENPW (Supplementary Fig. 11). Next, we conducted validation of the Group 3 GWAS back in IMAGEN. We began by calculating the polygenic risk score for SNPs (N_SNP = 4) residing in CENPW (CENPW score) and across the whole genome (PRS) that are associated with Group3-reweighted GMV for each adolescent in IMAGEN, tested for the differences of PRS among these groups, and correlated the PRS with neurocognition and behavioral risk factors. Consistent with the Group 3 GWAS, we observed higher CENPW score in Group 3 relative to Groups 1/2 (Fig. 3c) and positive correlations between CENPW score and improvement of neurocognition and conduct problems (Fig. 3d). Similar results were obtained for PRS (Supplementary Fig. 12).

No genome-wide significant SNPs were identified in the Group 2 GWAS (Supplementary Fig. 13). However, the large overlap between the neurodevelopmental patterns and homogeneous genetic liability in Group 2 with Group 1 led us to reason that the differences of neurocognitive performances between Group 1 and 2 were quantitative (rather than qualitative) and might subject to the effects of environmental exposure. To test this, we performed epigenome-wide association studies (EWAS) in IMAGEN (Methods) using group label as the phenotype. A significant hypermethylation site cg06064461 (β = 25.40, P = 4.24×10^-8) (Fig. 4a) was identified and mapped to ATF2 and MIR933 on chromosome 2. ATF2 encodes a transcription factor of the activator protein-1 family, is ubiquitously expressed in the brain and was found to be associated with both neurodegeneration and neurogenesis^57–59. MIR933 shares a common promoter with ATF2 and offers neuroprotection against neurodegenerative diseases by regulating brain-derived neurotrophic factor (BDNF)⁶⁰. To validate the EWAS results, we correlated the methylation of cg06064461 with estimated GMV trajectory and peak GMV in Groups 1/2, and calculated the mediation effects of cg06064461 methylation in the adverse environment - neurodevelopment pathway. Consistent with the EWAS results, positive correlation between cg06064461 methylation and total GMV trajectory (r = 0.14, P < 0.001) (Fig. 4b) and negative correlation between cg06064461 methylation and peak GMV (r = -0.07, P = 0.020) (Fig. 4c) were observed. Although none mediation effects of cg06064461 methylation on the environment – neurodevelopment pathway showed statistical significance after correcting for multiple testing (Fig. 4d and Supplementary Table 12), nominal significance was observed for family affirmation, where higher level of family affirmation was associated with higher peak GMV through reduced cg06064461 methylation (Fig. 4e). Family affirmation was defined as behaviors that the parent puts in place to support or help children in various situations or to show them approval and affection and refers to parent-child relationship⁶¹. These results indicated that environmental exposure could contribute to disadvantaged neurodevelopment and neurocognition by inducing epigenetic changes of neurogenesis-related genes.

• Download figure
• Open in new tab

Fig. 4 Epigenome-wide association study (EWAS) identified significant signals associated with lowered neurodevelopment in Group 2.

(a) EWAS Manhattan plot in the IMAGEN population. Group 2 (n=463) (relative to Group 1, n=446) status was used as the phenotype, adjusting for potential confounders (Methods). One hypermethylated site cg06064461 achieved genome-wide significant effect (P<5×10^-8,P_adj<0.05) and was mapped to ATF2 and MIR933 on chromosome 2. (b-c) Validation of EWAS results in IMAGEN (n=909). cg06064461 methylation was positively correlated with total GMV trajectory (b; r=0.14, ***P<0.001) and negatively correlated with peak gray matter volume (GMV) (c; r=-0.07, *P=0.020), adjusting for potential confounders (Methods). (d) Proportion of the mediation effects through cg06064461 methylation in the environmental exposure - peak GMV pathway, adjusting for potential confounders. Environmental factors were sorted by P values of the mediation effects. Although none mediation effects of cg06064461 methylation showed statistical significance after correcting for multiple testing, nominal significance was observed between family affirmation and peak GMV. Childexp, child’s experience of family life; FamStress, family stressors; CTQ, Childhood Trauma Questionnaire. (e) Mediation model was conducted to analyze the direct and indirect effect of family affirmation on peak GMV, with cg06064461 methylation as the mediator. Results showed that cg06064461 methylation significantly mediate the relationship between family affirmation and peak GMV (β=0.005, *P=0.048).

Genetic variation had limited effects on the cognitive, mental health and socio-economic outcomes in mid-to-late adulthood

Both genetic vulnerability and structural neurodevelopment are well-established to have profound impact on one’s physical, social and mental well-being in mid-to-late adulthood^62–65. However, the neurobiological mechanisms through which the long-term effects of genetic variation could be manifested remain largely unknown. We conclude our study by testing in the UK Biobank whether, and to what degree, polygenic risk for delayed neurodevelopment could have impact on the long-term brain structure, cognition, social-economic outcomes and mental well-being.

Motivated by the Group 3 GWAS results, we first calculated the PRS and CENPW score of delayed neurodevelopment for each participant in UK Biobank and then correlated them with outcomes of interest. Both PRS and CENPW score were approximately normally distributed and negatively associated with total GMV among this population (Fig. 5a). Next, we inspected the association of PRS and CENPW score with GMV in multiple brain regions in UK Biobank and found that inferior temporal, fusiform, middle temporal, medial orbitofrontal and lateral orbitofrontal areas were among the most correlated ROIs with PRS of delayed neurodevelopment, and lateral orbitofrontal, caudal middle frontal, rostral middle frontal, insula and superior frontal areas were among the most correlated ROIs with CENPW score (Fig. 5b and Supplementary Tables 13-14). These were consistent with the worse spatial working memory among participants with higher PRS of delayed neurodevelopment and CENPW score (Supplementary Table 15). Finally, we conducted non-superiority tests of the correlation coefficients and found that correlations between PRS of delayed neurodevelopment and CENPW score and all outcomes of interest were smaller than 0.05 (Fig. 5c and Supplementary Figs. 14-15). This indicated that although polygenic risks contributed to delayed neurodevelopment during adolescence, their long-term influences on the cognitive, mental health and socio-economic outcomes were limited once neurocognitive abilities were fully developed.

• Download figure
• Open in new tab

Fig. 5 Genetically-predicted neurodevelopment had limited impact on socio-economic, cognitive and mental health outcomes in mid-to-late adulthood.

(a) Correlation between polygenic risk score (PRS) of delayed neurodevelopment and total gray matter volume (GMV) for participants in UK Biobank (n=337,199). Marginal distributions of PRS and total GMV were both normal. PRS and CENPW score both showed negative correlation with total GMV(r=-0.07, ***P<0.001 for PRS and r=-0.07, ***P<0.001 for CENPW score). PRS were averaged over different P value thresholds. (b) Correlation between averaged PRS of delayed neurodevelopment and CENPW score and regional GMV for participants in UK Biobank. Inferior temporal (r=-0.07, ***P_adj<0.001), fusiform (r=-0.06, ***P_adj<0.001), middle temporal (r=-0.06, ***P_adj<0.001), medialorbitofrontal (r=-0.06, ***P_adj<0.001) and lateralorbitofrontal (r=-0.06, ***P_adj<0.001) were among the ROIs with the strongest correlation with PRS, while lateralorbitofrontal (r=-0.06, ***P_adj<0.001), caudalmiddlefrontal (r=-0.05, ***P_ad-j<0.001), rostralmiddlefrontal (r=-0.05, ***P_adj<0.001), insula (r=-0.05, ***P_adj<0.001) and superiorfrontal (r=-0.05, ***P_adj<0.001) were ROIs having the strongest correlation with CENPW score. These were consistent with the results that participants with higher PRS of delayed neurodevelopment also had worse performance in spatial working memory in UK Biobank. (c) Inferiority test of the correlation between averaged PRS and socio-economic, cognitive and mental health outcomes (Supplementary Methods) indicated that polygenic risk of delayed neurodevelopment had limited effect on the long-term socio-economic, cognitive and mental health outcomes. Full results were displayed in Supplementary Fig. 14. Similar results were observed between CENPW score and these long-term outcomes (Supplementary Fig. 15). IMD, Indices of Multiple Deprivation; Scot, Scotland; Edu, the highest educational level; IQ, intelligence.

Participants

Genomic, neuroimaging, environmental exposure, behavioral and mental health related data used to identify adolescent neurodevelopmental patterns were obtained from the IMAGEN study. Individuals with GMV beyond 4 interquartile ranges (IQRs) in any ROI were considered as outliers and were excluded from the analyses. After applying the exclusion criteria, 1543 adolescents with at least two structural MRI scans from 14 to 23 years old were included in the analyses (Supplementary Table 16). The average number of structural MRI scans per participant was 2.63, with 974 adolescents having a total of 3 scans (at 14, 19 and 23 years, respectively) and 569 adolescents having a total of 2 scans (384 at 14y and 19y, 147 at 14y and 23y and 38 at 19y and 23y). In addition, genotyping data used in GWAS, validation of GWAS and investigation of the long-term impact were obtained from ABCD, IMAGEN and UK Biobank, respectively. A total of 11,760 participants aged between 9 and 11 years old from ABCD were included, with the average number of structural MRI scans per adolescent 1.68. Further, 502,409 participants aged between 37 and 73 years old from UK Biobank were included in the long-term analyses of structural brain development. Demographics and baseline characteristics of participants from the three large population cohorts were summarized in Supplementary Table 17. A full description of all population cohorts used in the analyses can be found in Supplementary Methods.

Data preprocessing

In brief, quality-controlled processed T1-weighted neuroimaging data were obtained from ABCD, IMAGEN, Human Connectome Project Development (HCP-D), HCP Young Adult (HCP-YA) and the Philadelphia Neurodevelopmental Cohort (PNC). Assessment of regional morphometric structure were extracted by FreeSurfer v6.0 using Desikan-Killiany (h.aparc) atlas for cortical regions, and ASEG atlas for subcortical regions. Quality check was performed according to FreeSurfer reconstruction QC measures. Detailed description of data collection and preprocessing is provided in Supplementary Methods.

Estimation of GMV developmental trajectory

GMV trajectory in each of the 44 ROIs was estimated for each adolescent using linear mixed effect regression model (lme4 1.1-31 package) (since at most three structural MRI scans were available for each adolescent, only random slope model could be robustly estimated). Empirical Bayes estimate of the random slope was extracted for each adolescent. Intracranial volume (ICV), sex, handedness and imaging site were used as covariates to adjust for potential confounding.

Principal component analysis (PCA) and group clustering

Dimension reduction via PCA was performed to individual GMV trajectories estimated in the 44 ROIs. The first 15 principal components (Supplementary Table 1), which explained 80% of the total variation, were used in the multivariate k-means clustering. The optimal number of clusters was selected based on Elbow method with the constraint that each cluster contain at least 4% of the overall population.

Permutation test

Permutation was conducted by shuffling the estimated GMV trajectory in each ROI simultaneously and re-performing the dimension reduction and multivariate clustering repeatedly over 1000 times. P value was calculated as the proportion of Between-cluster Sum of Squares / Total Sum of Squares ratio greater than the estimated ratio in the original sample across all 1000 permutations.

Comparison of GMV trajectory among groups

Pairwise comparisons of GMV developmental trajectories in each ROI among the three groups were conducted via t test. The top 5 ROIs with the largest absolute t values were selected as the top distinguishing ROIs between the corresponding groups (Supplementary Table 3). Cohen’s ds (calculated using effectsize 0.8.3 package) for these regions were provided in Fig. 1d and Supplementary Fig. 16.

Estimation of age and region-specific GMV development among groups

To illustrate the region-specific GMV development in an extended time frame ranging from late childhood to early adulthood, external neuroimaging data from several population cohorts were incorporated. This includes a total of 21,826 participants comprising of 11,811 participants aged 9-14y with 19,587 scans in ABCD, 652 participants aged 5-22y in HCP-D, and 1,587 participants aged 8-23y in PNC study. Since cubic model could not capture GMV trajectory beyond 23y and quadratic model could not utilize data before 9y, we used a reference curve estimated from cross sectional studies (HCP-D + PNC) in estimating the region-specific GMV developmental curve over 5-25y. Distance between GMV in IMAGEN and that in the reference population in the corresponding ROI was used as the dependent variable in the quadratic linear mixed effect model with random intercept and slope, adjusting for ICV and site. Empirical Bayes estimates of the random effects for each group were added to the population averaged estimates to yield the group-specific developmental curve.

Estimation of group-specific developmental curve of total GMV in IMAGEN

A two-stage estimating procedure was adopted. Optimal model was selected among a series of polynomial mixed effect models using likelihood ratio test. First, population ICV developmental curve over 5-23y was estimated using the above-mentioned population neuroimaging data. Quadratic linear mixed effect models with random intercept at the individual and study level were fitted. To estimate the developmental curve of total GMV in the reference population (ABCD + HCP + PNC), cubic model adjusting for ICV was selected with random intercepts at the individual and study level. To estimate the developmental curve of total GMV in the ABCD and IMAGEN population, cubic model adjusting for ICV was selected with random intercept and slope at the individual level. Empirical Bayes estimates of the random effects were extracted and averaged in each group. Population ICV estimated in stage 1 was used to fit group-specific curves. The 5^th and 95^th percentile of the group-specific total GMV were calculated as the 95% confidence interval at each age.

Estimation of peak total GMV in IMAGEN

To estimate the peak total GMV in the IMAGEN population, 1,113 participants aged 22-38y in HCP-YA study were added to the reference population. A similar two-stage estimating procedure was used and the optimal model was selected based on Bayesian information criterion (BIC) and likelihood ratio test. First, population developmental curve of ICV was estimated using mixed effect regression model (nlme 3.1-160 package) with random intercept and slope. Basis function involving centered age was determined as the natural spline with 5 degrees of freedom. Interaction effects between age and sex, sex and study were also included in the regression model. Estimated ICV at each age was retained for the following analysis. Next, linear mixed effect regression model with random intercept and slope was fitted for total GMV. Basis function involving centered age was determined to be B spline with 12 degrees of freedom. Interaction effects between age and sex, sex and study were included in the regression model. Peak total GMV was defined as the highest total GMV one can achieve during brain maturation.

Comparisons of environmental burden, neurocognition, behavior and mental disorder

To assess whether environmental burden, neurocognition, behavioral risk factors and mental symptoms differ by groups, we analyzed their longitudinal measurements at 14y, 16y, 19y and 23y in IMAGEN. Personal traits, including personality, temperament and characters, were obtained from the NEO Five-Factor Inventory (NEO-FFI) and temperament and character inventory (TCI-R). Environmental burden, including prenatal exposures (parental smoking, maternal drinking and maternal medical problems during pregnancy), birthweight, stressful life events, child trauma experiences, child’s experience of family life and family stressors, were obtained from Pregnancy and Birth Questionnaire (PBQ), life-events questionnaire (LEQ), Childhood Trauma Questionnaire (CTQ), and Family Stress Scale and Family Life Questionnaire from development well-being assessment interview (DAWBA). Neurocognitive performances were obtained from Cambridge Neuropsychological Test Automated Battery (CANTAB) tests, Monetary-Choice Questionnaire (KIRBY) and Stop Signal Task (SST) results. Behavioral assessments, including conduct problems and substance use, were obtained from strengths and difficulties questionnaire (SDQ), European school survey project on alcohol and drugs (ESPAD), Fagerstrom test for nicotine dependence (FTND). Mental health conditions, including ADHD and depression, were obtained from self-rated development well-being assessment interview (DAWBA), where Attention-deficit/hyperactivity disorder (ADHD) score was additionally calculated using parent-rated interview. A detailed description of these assessment instruments is provided in Supplementary Methods. Generalized linear models adjusting for sex, handedness and ICV were used for comparing these tests at baseline and at the last follow-up visit (19y for Pattern recognition memory (PRM), Affective Go-No Go (AGN) and Rapid visual information processing (RVP); 23y for all other tests). For child-rated ADHD and depression score, baseline scores were also included as covariates. Intra-Extra Dimensional Set Shift (IED) test was only available at age 23y, and parent-rated ADHD score was only available at 14y and 16y. Cohen’s d was calculated for each measurement after regressing out the covariates. False discovery rate (FDR) was used to correct for multiple testing within scales.

Quality control of genomic data

In this study, we performed stringent QC standards using PLINK 1.90. Individuals with > 10% missing rate and single-nucleotide polymorphisms (SNPs) with call rates < 95%, minor allele frequency < 0.1%, deviation from the Hardy-Weinberg equilibrium with P < 1E-10 were excluded from the analysis. For ABCD, we only selected subjects with self-reporting White ancestral origins using the public release 3.0 imputed genotype data, which was imputed with the HRC reference panel⁸⁷. Considering that ABCD is oversampled for siblings and twins, we randomly selected one participant within each family. For IMAGEN, details about preprocessing of genomic data can be found in previous reports ⁸⁸ and data was imputed with the HapMap3 reference panel⁸⁹. For UKB, we selected subjects that were estimated to have recent British ancestry and have no more than ten putative third-degree relatives in the kinship table using the sample quality control information provided by UKB. For more details, please refer to the official document for genetic data of the UKB (http://www.ukbiobank.ac.uk/scientists-3/genetic-data/). After quality control, we obtained a total of 4,244,228 SNPs and 7,662 participants in ABCD, 5,966,316 SNPs and 1,982 participants in IMAGEN, and 616,339 SNPs and 337,199 participants in UKB.

Calculation of genetic liability

For each individual, polygenic risk scores (PRS) of ADHD, ASD, EA and IQ were calculated based on the public GWAS summary statistics^{51, 90–92} using PRSice v2.3.3. For ADHD, ASD and IQ, optimal p-value thresholds were determined based on the best-fit R2 using parent-rated psychiatric scores for ADHD and ASD, and the total WISCIV score. For EA, variants were selected using a P value threshold from 5e-08 to 1 with a step of 5e-05 and an average score under each P value threshold was calculated.

GWAS and validation

Since it was difficult to estimate individual GMV developmental trajectory in ABCD with limited number of structural MRI scans per participant and limited age range, we calculated the group-reweighted GMV as a proxy phenotype. There are several underlying assumptions in this calculation. Firstly, it assumes that all brain regions exhibit a comparable linear change from childhood to adolescence. Secondly, it assumes that the participants from ABCD and IMAGEN are drawn from a homogeneous population. Once again, we only included individuals in ABCD with self-report White ancestral origins. ROI-specific loading contributing to group classification (Group 2 vs Group 1, Group 3 vs Group 1, Group 3 vs Group 2, and Group 3 vs Groups 1/2) were obtained by regressing baseline GMV in 44 ROIs adjusting for age, sex, handedness and site in IMAGEN. Logistic regression model was used as the classification model and top 10 ROIs with the largest loadings were used to calculate the group-reweighted GMV in ABCD. Since results remained similar when comparing Group 3 vs Group 1 and when comparing Group 3 vs Group 2 (Supplementary Figs. 17-18), we combined Group 1 and 2 for increased statistical power, and performed the GWAS to investigate the genetic variations associated with Group 3 vs Groups 1/2 (delayed brain development). GWAS was conducted in the White population adjusting for sex, scanner effect and top 20 PCA components using Plink 2 (Supplementary Fig. 19). To ensure the validity of group reweighted phenotype, we correlated the Group-3 reweighted GMV and neurocognitive assessments (Game of Dice Task, Delay Discounting Task and NIH Toolbox) in ABCD (Supplementary Methods). Gene-based association analysis was conducted via MAGMA (version 1.10) using raw genomics data with the same covariate adjustment. To validate the GWAS results, polygenic risk scores for SNPs residing in CENPW (referred to as CENPW score) and across the whole genome (referred to as PRS) were calculated. Four SNPs were obtained by clumping within 250 kb upstream and downstream of CENPW (chr6:126339789-126483320) using Plink 2. PRS was calculated using PRSice using the most predictive P threshold for group-reweighted baseline GMV (r=0.06, P=0.033). Distribution of PRS between Group 3 vs Groups 1/2 and correlation coefficients between PRS and neurocognition, behavior and mental disorder at age 14y and 23y were obtained. FDR was used for multiple tests correction within scales.

EWAS, gene-specific methylation analysis and results validation

EWAS was performed among 909 adolescents in Group 1 and Group 2 in IMAGEN. Methylation data were collected using the Illumina Infinium HumanMethylation450 BeadChip. Locus-specific genome-wise methylation analysis was conducted and beta values at each Autosomal CpG site between groups was compared using logistic regression adjusting for sex, experimental batches (recruitment center and acquisition wave), the first two principal components of methylation composition and the first four principal components of estimated differential cell counts. Statistical significance was set as false discovery rate (FDR) adjusted p-value < 0.05. Next, we aimed to investigate the association between CpG site and gene methylation with environmental factors of interest. We conducted mediation analyses (sem function in the lavaan 0.6-12 package) and estimated the total effect of childhood experiences on estimated peak GMV and the indirect effect mediated by cg06064461 hypermethylation. Sex, batches effects, methylation composition components and differential cell count components were included as covariates. Total, direct and indirect effect and their standard deviations were estimated using 1000-iterated nonparametric bootstrap approach. False discovery rate (FDR) was used to correct for multiple testing within scales.

Long-term impacts of neurodevelopment in UK Biobank

Socio-economic, cognitive and mental health outcomes were obtained at baseline visit among participants in UK Biobank. Socioeconomic conditions were assessed by household income, jobs involved in physical activity and Indices of Multiple Deprivation (IMD) (education, employment and income scores). Cognition was assessed by fluid intelligence and the highest educational attainment. The educational level was divided into four ordinal categories: (1) College or University degree; (2) A levels/AS levels, NVQ or HND or HNC, other professional qualifications or equivalent; (3) O levels/GCSEs, CSEs or equivalent; (4) None of the above. Mental health was assessed by diagnosis of mental disorders (anxiety, depression and mania), summary score of neuroticism and self-reported mental symptom appearances. A detailed description of assessment instruments used in the analysis can be found in Supplementary Methods. To estimate the long-term effect of delayed neurodevelopment, we calculated CENPW score and RPS according to the results of Group 3 GWAS and correlated these scores with outcomes of interest after regressing out the age effect at recruitment, site and gender. It should be noted that these scores only reflect a genetic predicted risk for delayed brain development. Given the large age gap between participants in UK Biobank and IMAGEN, it is challenging to disentangle the long-term impacts of neurodevelopment from those due to potential environmental confounding in mid-to-late adulthood. Therefore, this analysis only serves to explore the potential long-term influence of genetically predicted delayed neurodevelopment and does not account for potential confounding due to environmental factors. Similarly, we assume the homogeneity of study populations between IMAGEN and UK Biobank. For PRS calculation, we used P value thresholds from 5E-08 to 1 with a step of 5E-05 and calculated an average PRS score for each individual. Due to the large sample size and easily-obtainable statistical significance, inferiority tests against 0.05 were conducted against the null hypothesis that the absolute correlation coefficient was less than 0.05.