Introduction

The hippocampus carries notable significance concerning its phylogenetic origins^1-4 and its engagement in a broad range of cognitive functions and behaviors beyond its initial implication in memory and spatial navigation^5-12. It was traditionally assumed that the hippocampus undergoes its most significant development in infancy and early childhood and achieves full maturity by the end of this period¹³. Yet, the advances in high-quality magnetic resonance imaging (MRI) data in recent years have allowed researchers to discern that the hippocampus undergoes protracted growth beyond early childhood, which may play an important role in the maturation of sophisticated cognitive functions during youth¹⁴. Notably, this development was primarily characterized by changes affecting subregions along the cytoarchitectural and anterior-posterior (i.e. long-axis) gradients of hippocampus^15-20, highlighting the importance of understanding the organizational principles of the hippocampus to map its structural and functional development, and its role in cognitive maturation.

Despite such progress, theories centered solely on the hippocampus fail to fully explain its extensive functions ^13,21, highlighting the need for a systems-level perspective that embeds the hippocampus within brain-wide networks^22-25. Previous findings have illuminated the complex, multidimensional organizational principles of hippocampal-cortical connectivity gradients, which underpin key hippocampal functions in local information processing, large-scale neural network integration and the mapping of behavioral phenotypes^26,27. Accumulating research has further revealed the potential of the hippocampal-cortical connectivity and its gradients in accounting for the protracted developmental profile of cognition during youth^28-32. Nevertheless, a comprehensive elucidation of this developmental process still requires integrating the intricate, multidimensional organizational architecture of the hippocampus itself.

Moreover, these organizational principles are likely sculpted by multiscale factors—from macroscopic geometric constraints to molecular genetic mechanisms—that collectively modulate hippocampal structure and function, and ultimately drive the reorganization of its connectivity gradients. At the macroscopic level, geometric constraints play a crucial role in organizing hippocampal activity³³, fundamentally influencing its functional topography. Moving to the mesoscopic scale, myelin content has been found to exhibit a significant medial-lateral gradient along the hippocampus²⁶. Myelination during nurturing process is identified as a key factor in neural plasticity and circuit refinement, affecting neural activity synchronization^34-38. Finally, at the molecular level, emerging evidence reveals systematic variations in gene expression along the anatomical long-axis of the hippocampus^39,40, underscoring the importance of genetic factors in determining both its structural organization and functional properties. Importantly, the prominence of cortical hierarchical gradient was found to be established during adolescence^41,42 and may support the maturation of varies higher-order cognitions, including executive functions^42,43. However, the mechanisms underlying the association between hippocampal gradient organization and cortical functional maturation, as well as how multiscale factors relate to their developmental trajectories, remain largely unknown.

We hypothesized that the hippocampus’s multidimensional gradient organization undergoes significant reorganization during youth, shaped by multiscale factors including geometry, myelination, and gene expression, and that this reorganization is associated with the concurrent maturation of cortical functional organization and the development of high-order cognition. To test this hypothesis, as illustrated in Figure 1, we first derived the triple gradient organization of the hippocampus through analysis of hippocampal-cortical connectomes and geometric modes. We then projected the hippocampal gradients onto the cortex to elucidate their roles in cortical functional processing. Subsequently, we examined the developmental trajectories of this gradient organization and their impact on cortical hierarchy maturation, episodic memory, and executive function. Furthermore, we identified the key drivers of hippocampal gradient reorganization, including myelination and transcriptomic factors. Our findings were replicated across two independent datasets, ensuring the robustness of our results. This comprehensive study integrates molecular, cellular, geometric, and functional analyses to shed light on the hippocampal triple gradient organization, its development, and its role in cortical hierarchy maturation and high-order cognitive development in youth.

Figure 1.

Results

Our primary findings were obtained from a large-scale MRI dataset sourced from the Human Connectome Project Development (HCP-D)⁴⁴, encompassing 652 typically developing participants aged 5-21 years (351 females). The validation results were obtained from two longitudinal multimodal MRI datasets (Enhanced Nathan Kline Institute - Rockland Sample (NKI-RS) and Children School Functions and Brain Development Project (CBD, Beijing Cohort)) (Supplementary Fig. 1). Structural and resting-state functional MRI (fMRI) underwent HCP minimal preprocessing⁴⁵ with an additional denoising process for fMRI (see Methods). To reduce the partial volume effect and adapt to the diverse variants of hippocampal folding observed among individuals^26,46, we employed an innovative deep learning technique integrated with a topological constraint approach (HippUnfold⁴⁷) to segment the hippocampus and generate its mid-thickness surfaces (483 vertices) featuring coordinates of two intrinsic orthogonal geodesic axes (posterior-anterior (P-A) and proximal-distal (P-D) axes)(Figure 2a; Supplementary Fig. 2; see Methods).

Figure 2.

Triple gradient organization of the hippocampus

We first derived the group-level hippocampal-cortical functional topography in youth. We constructed hippocampal-cortical FCs between each hippocampal vertex and each cortical region (defined by the Glasser atlas⁴⁸), and then applied a nonlinear dimension reduction algorithm to obtain the hippocampal functional gradient (See Methods and Figure 2a). The group-level hippocampal gradient was initially derived from the averaged FCs across all participants in the HCP-D dataset (Figure 2b).

The elbow point in the explained variance plot (Extended data Fig. 1a) indicates that the first five gradients collectively account for a substantial portion of the variance in hippocampal FCs (72.9% and 71.5% for the left and right hemispheres, respectively. The first principal gradient (G1) accounts for 39.5% (left hippocampus) and 38.6% (right) of the variance and is strongly linearly correlated with P-A coordinates (left: r² =0.823, right: r²=0.784, p<0.001, corrected for spatial autocorrelation, as were all subsequent analyses; see Methods). These findings are illustrated in Figure 2b. We thus termed this gradient as the linear long-axis gradient (ILG). The second gradient (G2) explains 16.6% (left) and 15.5% (right) of the variance, showing a significant quadratic relationship with P-A coordinates (left: r² =0.524, right: r² =0.474, p<0.001, Figure 2c). This gradient hence was named the quadratic long-axis gradient (qLG). The third gradient (G3) explains 7.7% and 8.3% of the variance for the left and right hemispheres, respectively, displaying a significant cubic association with P-D coordinates (left: r² =0.210, right: r² =0.320, p<0.001, Figure 2d). This gradient exhibits similarities to the previously established cytoarchitectonic variations within the hippocampus^4,27. We hence termed it the cubic iso-allocortical gradient (cIG). Noted that no significant evidence shows the relationship between the triple hippocampal gradients and the alternative anatomical axis (Extended Data Fig. 2). We refrained from analyzing G4 and G5 due to the lack of clarity regarding their biological significance (Extended data Fig. 1b).

To investigate the age-related changes of this triple gradient organization, we derived the averaged functional gradients within four age-specific groups (5-9, 10-13, 14-17, 18-21 years, Extended Data Fig. 3a and b). The range of ILG expanded in early childhood but contracted in adolescence. The range of qLG consistently expanded through these stages, while the cIG remained relatively stable.

Additionally, we examined the developmental trajectories of explained variance and distribution characteristics (including standard deviation, range, skewness, and kurtosis) for the triple principal gradients using generalized additive models (GAMs) (Figure 1e, f and Extended Data Fig. 4). We fit the model with a smooth term for age as a fixed effect; sex, scanning site, and in-scanner head motion as linear covariates; and hemisphere as a random effect (see Methods). Our results revealed a progressive dominance of the , as indicated by an increase in explained variance (Δ Adjusted R² = 0.0192, p = 1.15 × 10⁻⁵). Furthermore, the standard deviation of ILG changed significantly, with an increase during childhood (5-11 years), a decrease during adolescence (12-17 years), and subsequent stability in adulthood (18-21 years) (Δ Adjusted R² = 0.0127, p = 2.04 × 10⁻⁴). The qLG demonstrated increasing differentiation as evidenced by the standard deviation (Δ Adjusted R² = 0.0062, p=0.003), without significant variation in its explained variance (Δ Adjusted R² = 0.0005, p=0.583). The cIG displayed a decrease of dominance as indicated by the explained variance (Δ Adjusted R² = 0.0047, p=0.011) and a relatively stable differentiation (Δ Adjusted R² = 0.0011, p=0.144). Analysis of the other characteristics of the triple functional gradients are presented in Extended Data Fig. 4.

Reorganization of hippocampal triple gradients in youth is linked to cortical hierarchy maturation and cognitive functions

Due to the profound developmental changes in the cortex during youth, we further investigated the link between the triple functional organization of the hippocampus and cortical functional organization. We initially constructed whole-brain maps of cortical contributions to each hippocampal gradient by projecting them to the cortex (Figure 2a & b; see Methods), forming group-averaged projection maps (Figure 3c, and Extended Data Fig. 5a). The Extended Data Fig. 5b & d illustrate how projection values differentiate distinct cortical-hippocampal functional connectivity (FC) profiles. Specifically, cortical regions with the highest (top 5%) and lowest (bottom 5%) projection exhibited markedly segregated FC patterns along the corresponding hippocampal gradients, underscoring that projection strength serves as a key dimension segregating hippocampal-cortical interaction modes.

Figure 3.

Intriguingly, when we tested the associations between the hippocampal projection patterns and the first ten corticocortical functional gradients⁴⁹, we found that these projection maps showed a striking correspondence with distinct cortical intrinsic functional organization (Extended Data Fig. 5c & e). We found a significant correlation between the projection patterns of ILG and cortical modulation-representation (MR) gradient corresponding to the 3^rd corticocortical functional gradient (r =−0.760, p<0.001). RMG differentiates attention and frontoparietal systems (modulation) from the sensorimotor and default mode systems (representation); The qLG significantly correlated with cortical functional hierarchy (FH) that corresponds to the 1^st corticocortical functional gradient (r =0.775, p<0.001). FH denotes the cortex organization where information flows from simpler low-level to complex high-level processing; While the reflects the internal-external gradient (IEG) that corresponds to the 5^th corticocortical functional gradient (r=0.469, p=0.007). IEG differentiates limbic, default mode, and ventral attention systems (internally oriented processing) from visual, dorsal attention, and frontoparietal systems (externally oriented, goal-directed processing). The associations between hippocampal projection maps and other cortical functional gradients can be seen in Supplementary Fig. 3.

Given the observation that hippocampal gradients closely correlated with the cortical hierarchy in youth, we further investigated whether the development of hippocampal gradient is linked to the maturation of the cortical hierarchy. We calculated the correlations between individual gradient projections and the adult cortical FH map and examined the developmental pattern of this association using GAMs. Significantly, we observed a strengthening link between triple gradient projections and the cortical FH map as individuals aged (ILG: Δ Adjusted R² = 0.008, p=0.003; qLG: Δ Adjusted R² = 0.020, p= 3.75 × 10⁻⁷; cIG: Δ Adjusted R² = 0.026, p= 4.51 × 10⁻⁷. See Figure 3d).

We also characterized the maturational changes in the hippocampal gradient projections in individual cortical regions. The developmental trajectories of the lLG projection and corresponding significant age effects (FDR-corrected p<0.05) are depicted in Extended Data Fig. 5f. Consistent with the development of the lLG shown in Figure 2e, f, the development of lLG projection also presents two inflection points, which coincided with the transition into adolescence (around age 12) and into adulthood (around age 18). Moreover, changes in lLG projection were particularly prominent in regions associated with the frontoparietal and ventral attention systems (Figure 3e, showing only regions that remained significant after FDR correction), indicating significant changes in their FC differentiation pattern along the hippocampal lLG. Furthermore, the range of qLG projection exhibited a consistent increase throughout the entire period (see the middle panel in Extended Data Fig. 5f), in line with the findings from the qLG range shown in Extended Data Fig. 4. We observed significant developmental changes in qLG projection predominantly within the ventral attention system (Figure 3f). In addition, the range of cIG projections exhibited a consistent decrease (see the right panel in Extended Data Fig. 5f), consistent with the findings from the explained variance changes shown in Figure 2e. The significant developmental changes of the cIG projections were in sensory-motor and default mode systems (Figure 3g).

This reorganization of triple hippocampal gradients, which correlates with increasingly complex hierarchical cortical processing, may link to the maturation of episodic memory and executive function (EF)—a set of top-down mental processes potentially associated with cortical hierarchical organization⁴³. We conducted 10-fold cross-validated LASSO-PCR model (the penalty parameter was selected via nested 5-fold cross-validation to identify the optimal value) using hippocampal gradient projections to predict the cognitions, while correcting for potential confounding variables including age, sex, scanning site, and in-scanner head motion, and applied FDR correction to all resulting p-values (12 gradient-cognition pairs) (Figure 3h). Interestingly, episodic memory is best predicted by lLG (r = 0.2259, FDR-corrected p < 0.001), followed by qLG (r = 0.119, FDR-corrected p = 0.019) and cIG (r = 0.139, FDR-corrected p = 0.008). In contrast, cognitive flexibility is most strongly linked to qLG (r = 0.447, p < 0.001), with lLG (r = 0.430, FDR-corrected p < 0.001) and cIG (r = 0.251, FDR-corrected p < 0.001) also contributing. While working memory is primarily associated with cIG (r = 0.197, FDR-corrected p = 0.004) over lLG (r = 0.081, FDR-corrected p = 0.127) and qLG (r = 0.049, FDR-corrected p = 0.354). Inhibitory control is similarly best predicted by cIG (r = 0.218, FDR-corrected p < 0.001), followed by lLG (r = 0.137, FDR-corrected p = 0.013) and qLG (r = 0.074, FDR-corrected p = 0.169).

Collectively, these results underscore the important and distinct role that the reorganization of hippocampal triple gradients during youth plays in the maturation of cortical hierarchy—particularly within frontoparietal and ventral attention systems—and in the development of cognitive functions, including episodic memory and executive function.

Progressive relaxation of geometric constraints on hippocampal gradients in youth

We aimed to investigate whether geometric constraints characterized by Laplace eigenmodes shape the triple dominant functional organization of the hippocampus. We derived the geometric eigenmodes of the hippocampus for each participant using the Laplace–Beltrami operator. Figure 4 a, b, and c showcase the group-averaged first three eigenmodes are significantly correlated with the triple principal functional gradients (lLG : left, r =0.90, right, r =0.89, p<0.001; qLG: left, r =0.72, right, r =0.67, p<0.001; cIG: left, r =0.48, p=0.002, right, r =0.38, p=0.014).

Figure 4.

Building on the above findings of significant geometry-function associations, we further employed GAMs to fit the age-related trajectory of coupling between hippocampal geometric modes and their closely correlated hippocampal gradients (Figure 4d, e, and f). The magnitude and direction of age effects were also derived. Interestingly, the geometry-function coupling of the λ and cIG display a continuous decrease during this period (lLG: Δ Adjusted R² =0.055, p<1 ×10⁻¹⁶; cIG: Δ Adjusted R² =0.003, p=0.020), whereas the coupling of the qLG display a decline with a slight fluctuation (Δ Adjusted R² =0.007, p=0.026).

Development of hippocampal triple gradients parallel myelin maturation

We first aimed to investigate whether regional myelin content, as measured by T1w/T2w intensity⁵¹, is associated with the organization of the hippocampal triple gradients. Our analysis revealed that the group average T1w/T2w intensities exhibited a significantly stronger correlation with cIG than with lLG and qLG (cIG: left, r² =0.155, p=0.039, right, r² =0.161, p=0.042; lLG: left, r² =0.008, p=0.017, right, r² =0.031, p=0.007; qLG: r² =0.013, p=0.063. Figure 5a, b).

Figure 5.

Considering the pivotal role of myelination in regulating and restraining neural plasticity^36,52, our investigation sought to elucidate the potential relationship between the maturation of myelin content and the development of hippocampal triple gradient profiles. We initially partitioned the group-averaged T1w/T2w intensity (a structural MRI measure sensitive to myelin content; Figure 5a) into nine uniform bins to ensure consistency. Then, within each bin, we examined the developmental trajectories of mean myelin content and mean gradient values for each of the triple gradients. Consistent with our expectations, myelin content in all bins increases continuously throughout childhood and adolescence, reaching a relatively stable level when individuals enter adulthood. Furthermore, the age effect demonstrates a gradient along the myelin axis, i.e., hippocampal regions with higher myelin content exhibited more pronounced age-related changes in myelin content. To distinguish between positive and negative age effects, the Δ Adjusted R² of GAM fits that decreased with age were assigned a negative sign for visualization (Figure 5c; Extended Data Table 1 provides the Δ Adjusted R² and FDR-corrected p value for each bin). (Figure 5c, Extended Data Table 1 provides the Δ Adjusted R² and FDR-corrected p value for each bin). Interestingly, the developmental variability of the hippocampal triple gradients also displayed a prominent pattern aligned with the myelin axis (see Figure 5d as well as Extended Data Table 1). We identified a significant spatial correspondence in developmental variability between the refinement of myelin content and the hippocampal triple functional gradients. Notably, these findings remained consistent regardless of the number of bins used (see Extended Data Fig. 6 and Supplementary Tables 1–4).

Transcriptomic substrates of the hippocampal triple gradients across development

We initially conducted a transcriptomic association analysis to explore the presence of gene expression variations along the triple hippocampal functional organization. We obtained normalized gene expression data from 58,692 probes, extracted from 173 hippocampal samples of six deceased human donors provided by the Allen Human Brain Atlas. We then employed a LASSO-PCR algorithm to predict gradient value using gene expression profile (see Methods), through 10 repeated ten-fold cross-validation (CV). We found a significant association between gene expression and the lLG (lLG: left, r² =0.580, p=0.001, right, r² =0.239, p=0.002; qLG: left, r² =0.068, p=0.179, right, r² =0.033, p=0.252; cIG: left, r² =0.011, p=0.334, right, r² =0.012, p=0.346. Figure 6a).

Figure 6.

To delve into the biological processes and molecular functions associated with the important gene set responsible for predicting the hippocampal triple gradients (details see Methods), we performed enrichment analyses for biological pathways and gene ontology utilizing Metascape⁵³. The left panels within Figure 6b, c and d show the most significant enrichment term clusters derived from the identified important genes, with FDR-corrected p values below 0.01. We adopted enrichment networks by connecting enriched terms with Kappa similarities exceeding 0.3 to depict similarities between term clusters and redundancies within clusters. The most significant gene set linked to the lLG prediction exhibits enrichment in terms associated with neuroactive signaling (e.g. neuroactive ligand-receptor interaction, chemical synaptic transmission, calcium ion binding, and behaviors elicited by internal or external stimuli), regulation of anatomical structure morphogenesis and neural development, and regulation of stress hormones (i.e., positive regulation of cortisol secretion). Furthermore, the most significant gene set linked to qLG displays enrichment in terms correlated with neurodevelopment and synaptic plasticity (e.g. synapse organization, brain development, neuron projection development, synaptic signaling, and neurotrophins signaling pathway), homeostatic mechanisms (regulation of cellular response to stress), and protein localization. We did not obtain significant enrichment terms correlated with cIG. Detailed information on gene annotation and enrichment for the triple gradients can be found in Supplementary Data 1–3, respectively.

To ascertain whether the noteworthy gene set related to the triple gradients exhibits enrichment specifically during the developmental period under investigation within the hippocampus, the identified important gene sets were compared against developmental expression profiles from the BrainSpan dataset (http://www.brainspan.org/) by employing the developmental-specific expression analysis (SEA) tool⁵⁴. The results underscored that the important gene expression for lLG exhibited enrichment within the hippocampus during mid/late childhood, adolescence, and young adulthood. This enrichment was notably accentuated in the mid/late childhood period (depicted in the right panel of Figure 6b), which aligned with the fastest developmental period in the hippocampal lLG (see left panel in Figure 4e and f). Comparable results were also obtained from the developmental enrichment analysis of the important gene list for qLG (right panel of Figure 6c). We did not obtain significant enrichment developmental periods correlated with cIG (right panel of Figure 6d).

These results suggest that the molecular processes —encompassing neurodevelopment, neuroactive signaling, and stress hormone regulation—likely contribute to the development of hippocampal dual long-axis gradients in youth.

Replications of hippocampal triple gradients and development

Thus far we have demonstrated the fundamental role and developmental significance of the hippocampal triple gradient organization, we aimed to further investigate whether these organizational principles generalize to two independent datasets: the CBD and NKI-RS. The CBD dataset includes 300 healthy children (140 females, aged 6–13 years, 478 total scans), while the NKI-RS dataset comprises 240 healthy youth (109 females, aged 6–21 years, 363 scans). Our findings confirm that the functional gradients and geometric modes observed in the HCP-D dataset are well replicated in both datasets (Extended Data Fig. 7a & b).

We further assessed the generalizability of the developmental trajectories of the hippocampal triple gradients and their cortical projections across independent datasets. Our primary focus was the NKI-RS dataset, given its broader age range compared to the CBD dataset, which also encompasses the age range of the HCP-D dataset. The results showed that the developmental trajectories of the hippocampal triple gradients were largely replicated in the NKI-RS dataset (Extended Data Fig. 7c, d). Specifically, while not reaching statistical significance, the lLG exhibited a similar trend to that observed in the HCP-D dataset (Δ Adjusted R² = 0.007, p = 0.07 for standard deviation). Consistent with our findings, the qLG demonstrated relatively stable explained variance (Δ Adjusted R² = 0.001, p = 0.47), while the cIG showed a significant decline in explained variance (Δ Adjusted R² = 0.022, p = 0.0007) and a relatively stable standard deviation (Δ Adjusted R² = 0.001, p = 0.61). For the age range in the CBD dataset, we observed a significantly increasing trajectory (Δ Adjusted R² = 0.001, p = 0.02) in the standard deviation of the qLG, consistent with the findings from the HCP-D dataset within the corresponding age range (Supplementary Fig. 4).

Notably, we observed an increasing association between hippocampal triple gradient projections and the cortical functional hierarchy with age in NKI-RS dataset (lLG : Δ Adjusted R² = 0.004, p = 0.06; qLG: Δ Adjusted R² = 0.014, p = 0.008; cIG: Δ Adjusted R² = 0.016, p = 0.006), which aligns with the findings from the HCP-D dataset, confirming the consistency of this developmental relationship across independent samples (Extended Data Fig. 7e).

Together, these results show that the hippocampal triple gradients, their developmental trajectories, and increasing alignment with the cortical functional hierarchy are robustly replicated in the CBD and NKI-RS datasets, confirming the generalizability of these organizational principles in youth.

Discussion

The human hippocampus and cortex depend on a set of fundamental organizing principles to process information, regulate behavior, and facilitate various intricate functions. Our findings enrich the comprehension of the hippocampal organizational principles in terms of hippocampal-cortical connectivity, hippocampal geometric eigenmodes, myelin pattern, gene expression, as well as their distinct relationships with cortical functional processing in youth and their relevance to the maturation of cortical functional hierarchy and cognition. Specifically, the hippocampal geometric modes significantly constraints the triple hippocampal functional gradients, underscoring the fundamental role of geometry in shaping hippocampal spontaneous functional activity, and potentially by influencing wave dynamics^33,55,56. Additionally, we found that the gene expression pattern can predict linear long-axis gradient best compared to the other two gradients while the myelin content predicts cubic iso-allocortical gradient best (Figure 5b and Figure 6a). Since hippocampal linear long-axis gradient is more evolutionarily conserved and thus it may depend more on shared gene expression pattern^4,39. While cIG reflects the part of hippocampal-cortical connectome influenced by cytoarchitectonic differentiation^4,27. In contrast, hippocampal quadratic long-axis gradient showed a prominent association with cortical functional hierarchy (Figure 3c and Extended Data Fig. 5a, c and e). Given the greater dominance of hierarchical processing in humans⁴¹ and macaques⁵⁷ in contrast to other species such as marmosets⁵⁸ or rodents⁵⁹, quadratic long-axis gradient is probably not as evolutionarily conserved as the linear long-axis gradient. This intriguing finding also calls into question the understanding of the hippocampus as an evolutionarily conserved structure with respect to its structural and functional organization.

Hippocampal gradients sculpt cortical function and cognition

Additionally, we demonstrated that the maturation of hippocampal triple gradients is closely linked to the refinement of the cortical hierarchy, primarily through the reorganization of its connectivity with the frontoparietal and ventral attention systems. In line with our observations, recent evidence indicates that the cortical functional organizations that our hippocampal gradients interface with—specifically the sensory-association (SA, closely related with cortical functional hierarchy) and modulation-representation (MR) gradients—undergo profound reorganization during youth, shifting from a sensory-dominated architecture in infancy to a highly differentiated associative and controlled system in early adulthood⁶⁰. Additionally, these developmental shifts are primarily driven by changes in the ventral attention system⁶¹. Moreover, attention and frontoparietal control systems play a crucial role in coordinating bottom-up and top-down information flow^56,62-64. Furthermore, the triple hippocampal-cortical gradients take distinct contributions to episodic memory and executive functions, thereby highlighting their differential influence on cognitive processes. Specifically, the hippocampal long-axis gradient (aligned with cortical MR) predominantly correlates with episodic memory, consistent with a model in which the integration of sensory and default-mode representations alongside frontoparietal and attentional modulation relates to memory encoding and retrieval processes^18,65-70. Additionally, cubic iso-allocortical gradient (mapped to cortical internal–external processing) relates more to working memory and inhibitory control, reflecting coordination between internal limbic/ventral attention and external frontoparietal/dorsal attention systems^71-75. Finally, the quadratic long-axis gradient (associated with cortical functional hierarchy) best predicts the cognitive flexibility, suggesting frontoparietal-attention specialization and cross-hierarchy communication support task shifting.^43,76. This role of hippocampal-cortical gradients in scaffolding executive control is further underscored by lifespan evidence showing that the maturation of these cortical gradients—particularly the differentiation of the MR and SA axes—peaks during adolescence and is strongly associated with the emergence of high-level cognitive performance⁶⁰.

Developmental decoupling of hippocampal geometry and function

The geometric eigenmodes derived from both neocortical and non-neocortical structures have recently been demonstrated to be a more concise and precise representation of their macroscale functional activity compared to the connectome-based method³³. Furthermore, multiple recent studies have demonstrated that wave dynamics constrained by geometry and distance-dependent horizontal connectivity (i.e. dense short-range connections that decreases approximately exponentially with distance⁷⁷) may have a dominant influence on spatiotemporal brain activity^33,55,56. Therefore, the reduction trajectories in geometry-function coupling between hippocampal geometry and dual long-axis functional gradient suggest the development of long-range and non-horizontal connectivity and that the dominance of wave dynamics within the hippocampus may decrease during its maturation process. We inferred that this decoupling process may enable efficient communication both within and outside the hippocampus.

Myelination regulates developmental plasticity in hippocampal gradients

We found that developmental changes in the hippocampal triple functional gradients closely align with spatial patterns of myelination during youth, and the molecular process related to neurodevelopment, stress hormone regulation, and neuroactive signaling. Myelination—a key regulator of neural plasticity^36,37,52,78—exhibits spatial variations that parallel the neurodevelopmental axis of hippocampal gradients during this sensitive period. Notably, recent studies have suggested the potential existence of critical periods with heightened plasticity for hippocampus-dependent learning and memory, which have been proposed to be fundamental to cognitive development^79-82. In the primary sensory cortex, myelin formation and associated signaling via neurite outgrowth inhibitor (Nogo) receptors have been identified as pivotal mechanisms that may contribute to constraining or terminating these critical periods^36,52. In line with these findings, our results suggest that myelin not only regulates the developmental plasticity of hippocampal functional gradients but may also represent a candidate factor in orchestrating the timing and regulation of critical period plasticity during hippocampal maturation in youth. Future studies could consider combining pharmacological or chemogenetic approaches with neuroimaging techniques to establish connections between cellular- or molecular-level plasticity-regulating mechanisms and noninvasive human neuroimaging^83-85. This integrated approach can provide further insights into specific neurobiological mechanisms that may underlie putative critical periods in the development of the human hippocampus.

Molecular processes underpin development of hippocampal gradients

We demonstrated the potential role of distinct neurodevelopment process in forming the hippocampal dual long-axis gradient (i.e., linear and quadratic long-axis gradient) in youth. Specifically, gene expression related to the macroscopic patterning and structural formation (i.e., tissue morphogenesis) of the hippocampus shows a graded pattern along the linear long – axis gradient. Furthermore, processes governing the precise circuit wiring within this structure—namely, neuron projection development and synapse organization (including formation and pruning)—are most prominently enriched in the genes that define the quadratic long-axis gradient. This gradient may therefore emerge as a consequence of a transcriptional gradient optimized for experience-dependent circuit refinement. This fine-grained plasticity underlies complex cognition and may offer evolutionary advantages, consistent with the integration role of quadratic long-axis gradient within the cortical functional hierarchy and its protracted developmental trajectory. Notably, we found that the regulation of stress hormones also displayed variations along the linear long-axis gradient. This finding suggests a potential association between cortisol regulation and this hippocampal gradient, shedding light on the hippocampus’s role in stress-related cortisol activity^86,87. Since stress is linked to impaired neuroplasticity^88-90 and neurodevelopment^91-93—as well as its relevance to disorders like depression^94,95—future research could focus on how stress-related conditions or early-life adversity disrupt hippocampal gradient organization. Furthermore, mental training programs, known to reduce stress^96,97, could be investigated for their effects on the hippocampal gradient. In contrast, our linear modeling approach did not yield a robust genetic signature underlying the hippocampal cubic iso-allocortical gradient, suggesting that its organizational principle may be governed by more complex, nonlinear patterns of gene regulation.

Limitations and future directions

We replicated the main findings in hippocampal triple gradient organization, its cortical projections, development and the association with the maturation of cortical functional hierarchy on two independent longitudinal datasets. This replication underscores the robustness and generalizability of these findings across diverse populations and methodologies. We noted that the absence of pubertal hormone measures in the HCP-D dataset limits our ability to assess their impact on hippocampal triple gradients, highlighting the need for future studies incorporating hormonal assessments. By showing that the hippocampus is already organized into a distinctive triple-gradient architecture in youth and continues to undergo developmental refinement, this study provides a critical framework for tracing the emergence of hippocampal organization in early life and elucidating its role in shaping cognition and behavior.

Materials and methods

Participants and brain imaging dataset

HCP-D dataset (discovery)

For the primary scientific objectives of this study, we acquired multimodal MR images (T1w, T2w, and rs-fMRI) from 652 typically developing participants in the latest release of the Human Connectome Project (HCP) Development dataset. Following rigorous quality control procedures⁹⁸, 601 participants (330 female; age range = 5.58–21.92 years) were retained for final analysis. Participants were scanned at four sites using 3 Tesla Siemens Prisma platforms. Structural scans encompassed high-resolution MPRAGE T1-weighted images (TR/TI = 2,500/1,000 ms, TE = 1.8/3.6/5.4/7.2 ms, flip angle = 8°) at a resolution of 0.8 mm isotropic. A variable-flip-angle turbo-spin-echo T2-weighted sequence (TR/TI = 3,200/564 ms, turbo factor = 314) was also conducted with the same resolution. The incorporation of multiband acquisitions facilitated a subsecond temporal resolution for all functional images, characterized by a voxel size of 2.0 mm isotropic, TR/TE of 800/37 ms, and a flip angle of 52°. The participants were instructed to view a small white fixation crosshair on a black background during fMRI scanning.

Cognitive function test

We employed episodic memory and multiple cognitive measures related to EF, including cognitive flexibility, inhibitory control, and working memory, in the latest release of the HCP-D dataset. All these performance tests were conducted using NIH Toolbox (https://www.nihtoolbox.org/domain/cognition/). Specifically, episodic memory was estimated by Picture Sequence Memory Test in 472 participants by having participants recall sequences of illustrated objects and activities presented. Scoring was based on the number of correctly ordered adjacent picture pairs. Cognitive flexibility was measured by the Dimensional Change Card Sort Test, involving 463 participants. In this test, participants engaged in a bivalent matching task, wherein they matched test pictures (e.g., yellow balls and blue trucks) to target pictures. Initially, matching was based on one dimension (e.g., color), followed by a shift to the other dimension (e.g., shape) after several trials. The final score combined accuracy and reaction time. Accuracy (≤80%) alone determined the score; higher accuracy (>80%) added a speed component based on median reaction time, with faster responses earning more points. Additionally, inhibitory control performance was assessed with the Flanker Inhibitory Control and Attention Test in 462 participants. Here, participants focused on a target stimulus while ignoring flanking stimuli. It was scored identically to the cognitive flexibility test. Finally, working memory performance was measured using the List Sorting Working Memory Test, which evaluates both information storage and processing (manipulation) and included 473 participants. This test necessitated the rapid recall and arrangement of various visually and orally presented stimuli. The total score was the sum of items correctly recalled and sequenced. The detailed task description and scoring interpretation can be found in the documentation for the HCP-Aging Lifespan 2.0 Release and NIH Toolbox website.

CBD and NKI-RS dataset (replication)

The CBD dataset includes 300 children aged 6-14 years (140 females, 478 total scans, 165 children with one scan, 92 children with two scans, and 43 children with three scans; the scan interval was approximately one year) This data were collected by the Children School Functions and Brain Development Project (CBD, Beijing Cohort). Study procedures were approved by the Ethics Committee of Beijing Normal University, and written informed consent was obtained from all participants or their parents/guardians. The NKI-RS dataset comprises 240 typically developing children and adolescents aged 6–21 years (109 females; 363 total scans). All participants in these two datasets were cognitively normal, with no history of neurological disorders, mental health conditions, head injuries, or chronic physical illnesses. Imaging modalities included T1w images and rs-fMRI. For the CBD dataset, structural scans were acquired using a T1-weighted sequence on a Siemens Prisma scanner (TR/TE/TI = 2530/2.98/1100 ms, flip angle = 7°) with 1 mm isotropic resolution. Resting-state fMRI data were obtained with an echo-planar imaging sequence (TR/TE = 2000/30 ms, flip angle = 90°, slice thickness/gap = 3.5/0.7 mm) over an 8-minute scan duration, yielding 240 volumes in total. For the NKI-RS dataset, structural images were acquired using a T1-weighted MPRAGE sequence on a Siemens Tim Trio 3T scanner (TR/TE/TI = 2500/3.5/1200 ms, flip angle = 8°) with 1 mm isotropic voxel size. Resting-state fMRI data were collected with a multiband echo-planar imaging sequence (TR/TE = 645/30 ms, flip angle = 90°) at 3 mm isotropic resolution, comprising 900 time points. All scans underwent the same stringent quality-assurance procedures⁹⁹. Due to the limited spatial resolution and the relatively short duration of fMRI scans, as well as the short age range, we employed these data to validate our results at the group level.

Image preprocessing

Images from all three datasets underwent HCP minimum preprocessing⁴⁵ with several modifications to adapt to the specific developmental dataset.

1) Structural MRI

Images underwent gradient distortion correction and anterior commissure-posterior commissure (AC-PC) alignment, followed by brain extraction. Subsequently, a rigid body transformation with a boundary-based registration cost function¹⁰⁰ was employed to coregister the T1w and T2w images. To correct for bias fields, the square root of the product of T1w and T2w images was used, with nonbrain tissues thresholded out. The preprocessed images were subjected to nonlinear registration onto a standard template (MNI). Cortical surfaces were extracted in native space using FreeSurfer 6.0-HCP¹⁰¹, with a modification using T2w images to refine the pial surface obtained from T1w images. The individual native cortical surfaces were registered first to standard template space, then to FreeSurfer’s standard surface atlas (fsaverage), and finally to the Conte69 template¹⁰². Additionally, a T1w/T2w cortical myelin map was computed.

2) rs-fMRI

Images were registered to standard template space in a one-step spline resampling by concatenating a set of relevant transformations, including gradient distortion correction, distortion correction in the phase encoding direction, rigid body motion correction, two-step registration to the native T1w space with a combination of rigid body and boundary-based registrations, and nonlinear T1w-to-MNI registration. Following this, the processing steps involved the removal of the bias field calculated from the structural image, brain extraction, and normalization of whole-brain intensity. Thereafter, the volume time series in standard template space were mapped onto native cortical surfaces using a partial volume-weighted ribbon-constrained mapping algorithm⁴⁵, and then the signals on the surface were resampled and registered to the fsaverage and Conte69 template.

The subsequent preprocessing steps aimed at reducing structured noise, head motion artifacts, and physiological confounds. This included linear detrending, nuisance regression (incorporating motion parameters, white matter, cerebrospinal fluid, and global signals), and band-pass temporal filtering (0.01–0.08 Hz). Additional dataset-specific denoising steps (ICA-AROMA¹⁰³) were applied to further address residual motion and spatial noise.

Hippocampus segmentation and mid-thickness surface generation

We employed HippUnfold (v1.2.1)⁴⁷ (https://hippunfold.readthedocs.io), which integrated a deep learning method with a topological constraint strategy. This tool enabled the creation of mid-thickness surfaces featuring coordinates of two geodesic axes. Specifically, for each of the three geodesic axes—namely, inner-outer (I-O), posterior-anterior (P-A), and proximal-distal (P-D) axes—a Laplace field was calculated by solving Laplace’s equation across the hippocampal gray matter, spanning from source boundary tissue (representing one axis end) to the opposite sink boundary (the other end of this axis). Take the P-A axis as an example, considering the hippocampus as having cool and warm spots at the posterior and anterior ends, respectively. Solving the Laplace equation is like predicting how heat would evenly spread across the hippocampus, creating a smooth gradient from posterior to anterior, which topologically and intuitively maps the P-A coordinates.

In this framework, the outer boundary mirrors the cortical pial surface, and the inner boundary mirrors the white matter surface, thereby enabling the derivation of the mid-thickness surface of the hippocampus through the application of I-O coordinates (see Supplementary Fig. 2 for sagittal and coronal views of the three generated surfaces). This extracted surface also encompassed the P-A and P-D coordinates, allowing for the reasonable positioning of vertices along these two geodesic axes. Importantly, this mesh surface maintained a one-to-one correspondence of vertices across the hippocampi of distinct individuals. The vertex spacing was set at 2 mm, a choice that aligns well with the requirements for future fMRI signal sampling. Furthermore, due to its distinctive topology, the DG was delineated as a separate surface using a method analogous to the one employed for the remainder of the hippocampus, as elucidated above. The amalgamation of these procedures resulted in the creation of two separate mid-thickness surfaces corresponding to each hippocampus in native space. The two surfaces were composed of 483 vertices, with a specific allotment of 64 vertices dedicated to the DG. Notably, segmentations and surfaces were visually quality controlled for all participants.

To guarantee the segmentation accuracy and reliability in our study, we undertook rigorous quality control procedures with two steps derived from the recommendation of the HippUnfold and a further step of hippocampus volume analysis:

(1). We automatically evaluated hippocampus segmentation accuracy by computing the dice coefficient between HippUnfold-generated hippocampal masks and those obtained through diffeomorphic registration to a standard template. All Dice coefficient values exceeded 0.7, confirming the high sensitivity of the HippUnfold method to detect detailed hippocampal structures.

(2). We conducted visual quality control checks for the volumetric outputs of hippocampal segmentations and the topology or integrity of hippocampal surfaces across all participants. Snapshot images from individual participants’ segmentations were meticulously reviewed to ensure that no segmentation anomalies remained and that all segmentations met our standards of quality

Hippocampal feature mapping and gradient computation

Blood oxygen level-dependent (BOLD) time series and T1w/T2w intensities were mapped to each of the 483 vertices on the mid-thickness surfaces. Before the mapping of BOLD signals, fMRI images underwent a transformation into their native space. For optimal precision, we employed a same time series mapping approach akin to the method used for the cortical surface, known as the partial volume-weighted ribbon-constrained mapping algorithm. This approach leveraged both the outer and inner boundary surfaces to define the fMRI voxels situated within the hippocampal gray matter ribbon. The signal value attributed to each surface vertex represented the weighted mean of voxels entirely or partially contained within the hippocampal gray matter ribbon. In cases of partial voxels, their weighting was determined based on their fractional volume within the ribbon. These feature mapping procedures were performed utilizing the Workbench Command tool (v1.5.0). It is important to highlight that utilizing mid-thickness surfaces for feature sampling offers several advantages: 1) it can reduce the influence of the partial volume effects on the hippocampal boundaries. 2) it ensures heightened sensitivity and spatial precision in localizing BOLD signal sources within the hippocampus compared to volume-based methods. This advantage stems from the ability of mid-thickness surface to accommodate diverse variants in hippocampal folding⁴⁶, leading to improved localization accuracy. 3) This technique can increase statistical power due to the demonstrated superior alignment across distinct individuals resulting from surface registration^104,105.

The cortical surface was parcellated into 360 regions according to the Glasser atlas⁴⁸, and we computed the averaged time series within each region. We identified the hippocampal-cortical functional connectome gradients by performing the following steps. First, we correlated hippocampal and cortical time series and transformed them with Fisher’s z-transformation to obtain approximately normally distributed correlation coefficients. Then, we obtained the group-averaged z-transformed connectivity and inversely transformed it to correlation values. Diffusion embedding was performed to obtain the intrinsic connectome gradients based on the connectivity matrix using the BrainSpace tool (v0.1.10)¹⁰⁶ with default parameters in MATLAB (r2018a). Left/right hippocampal gradients underwent Procrustes alignment. We also derived individual gradients for each hippocampus based on the individual connectivity matrix.

Cortical projections of hippocampal gradients

The cortical projections for specific hippocampal gradients were computed by matrix production between cortical-hippocampal FCs and gradient values,

where Pⁱ ∈ R^360×1 represents the cortical projection corresponding to i-th order hippocampal gradient, Gⁱ ∈ R^483×1, and FC ∈ R^360×483 denotes the cortical-hippocampal FC. To clarify, the projection value for a particular cortical region r can be represented as the weighted sum of the gradient across all hippocampal vertices (equation (4)), where these weights were determined by the cortical-hippocampal FCs. Consequently, the cortical projection indicated the extent of similarity or alignment between the cortical-hippocampal FCs and the corresponding hippocampal gradient.

Derivation of hippocampal geometric eigenmodes

We derived the geometric eigenmodes for each hippocampus through the following process, which is adapted from Pang, et al. ³⁴: constructing a Laplace–Beltrami operator (LBO) based on the individual hippocampal mesh and subsequently solving the associated eigenvalue problem. To provide a more unified representation, a tetrahedral mesh was utilized to account for the complete three spatial dimensions of the hippocampus. Specifically, we first employed FreeSurfer’s mri_mc function, which implemented a marching-cubes algorithm, to generate a 2D surface of hippocampal boundary. This was achieved by tessellating the volumetric hippocampal mask at a resolution of 2 mm. Subsequently, we utilized Gmsh software (https://gmsh.info/) to transform the 2D boundary surface into a 3D tetrahedral mesh.

Subsequently, we formulated the Laplace-Beltrami operator (LBO) using this 3D tetrahedral mesh as a foundation. The LBO adeptly captures the spatial relationships and curvature between adjacent vertices in a localized manner. With the LBO established, we proceeded to solve the eigenvalue problem associated with it,

where Δ is the LBO and ϕ = {ϕ₁, ϕ₂, …} is the set of orthogonal geometric eigenmodes with a corresponding set of eigenvalues λ = [λ₁, λ₂, …]. These eigenvalues were sequentially organized based on the spatial frequency inherent to the spatial patterns of each respective mode, such that ϕ₁ represents the mode characterized by the lowest frequency. Note that the first eigenmode is a constant function as λ₁ is approximately equal to zero, and we discarded this mode for subsequent analysis. The general definition of the LBO is encompassed in equation (6), where denotes the component of the Riemannian metric tensor and x_i denotes the local coordinates on the surface. This tensor underpins the local geometric properties and distance measurements intrinsic to the manifold. We utilized the LaPy Python library (v0.6.0)^107,108 to extract the hippocampus’s geometric eigenmodes. Finally, we interpolated the eigenmodes from the tetrahedral mesh onto the mid-thickness surface, facilitating the comparison with hippocampal functional gradients.

Transcriptomic association analysis of the hippocampal gradients

1) Preparation of the human gene expression data

We collected post-mortem gene expression data from the Allen Human Brain Atlas. In summary, a total of 3702 tissue samples were gathered, including samples extracted from both hemispheres of two human brain donors and the left hemisphere of four additional donors. Microarray analysis and preprocessing were performed on each sample, leading to the quantification of gene expression across 58,692 probes. In line with preceding studies employing this dataset, the emphasis lies in discerning common patterns of human gene expression. Through the combination of samples from six donors, donor-specific characteristics including age and sex were regressed using linear models. Subsequently, standardized residuals were derived to remove individual variance from each probe. For the identification of hippocampal samples, we specifically chose those labeled as CA1 field, CA2 field, CA3 field, CA4 field, Subiculum, and Dentate Gyrus, encompassing both hemispheres. This selection constituted a total of 188 samples. Following this, the native hippocampal mid-thickness surfaces for all individuals underwent transformation into MNI space and were subsequently averaged to generate the population’s template mid-thickness surfaces. Upon contrasting the positions of the selected samples with the template surfaces, 15 samples were found to have MNI coordinates with a minimum distance from the hippocampal mid-thickness surface exceeding 6 mm and were excluded, leaving a total of 173 hippocampal samples. We obtained the nearest vertex on the template mid-thickness surfaces for each of the 173 hippocampal samples. The group-averaged gradient value was then associated with each sample.

2) Model fitting to predict hippocampal gradients

We employed a LASSO-PCR approach here using the scikit-learn package (version 0.21.3) in Python (version 3.8.16). This method effectively decreased the data dimensions and maintained gene co-expression networks while permitting sparse feature selection. We performed principal component analysis on the normalized gene expression data X, which was a matrix of dimensions 173(sample) × 58692(probe). The outcome of this process was the transformed data T, which then took the form of a 173(sample)×173(component score) matrix.

Then, we fit the Lasso linear regression model to predict the hippocampal gradient values based on the transformed data T. Then, we used the equation

where G_i is the i-th order hippocampal gradient, α = [α₁, α₂, …]^T, α_i is the regression coefficient of component score i, β = [β₁, β₂, …]^T, and β_i is the impact of probe i. The estimated impact of individual probes can thus be calculated based on the estimated using equation (10). Noted that the penalty parameter was selected via nested 5-fold cross-validation to identify the optimal value. We employed ten repeated 10-fold CV to estimate the model generalizability.

3) Model feature deconstruction to identify important genes

Regarding the unreliability of interpreting global feature weights derived from a Lasso model, due to the potential for significant reshuffling of feature importance resulting from feature inclusion or exclusion¹⁰⁹, we systematically eliminated the top 50 probes (positive-associated) and the bottom 50 probes (negative-associated) with the highest weights. Following this, we retrained the model and evaluated the CV accuracy of the updated model. This process was iterated until all 58,692 probes were removed. As a control, we repeated this same process iteratively removing 100 random probes instead of the 100 most important probes. To ensure methodological parity, we conducted this control procedure in parallel with the primary process. We visually examined the change in CV accuracy throughout successive rounds of probe removal. Inflection points were pinpointed at rounds where the CV accuracy experienced a decline and did not subsequently rebound. Notably for λ prediction, removing the initial set of 200 probes resulted in a distinct and unrecoverable decline in CV accuracy, underscoring the significance of this gene set for the model’s performance. Conversely, the stepwise removal of sets of 100 random probes resulted in a gradual and intermittent decline in accuracy, culminating at its lowest point solely upon the exclusion of a substantial majority of probes. Furthermore, because the feature deconstruction process for predicting qLG and cIG failed to robustly identify significant genes, we employed univariate correlation analysis between each probe and these two gradient components. This approach identified 600 probes associated with qLG and 70 probes associated with cIG.

4) Enrichment analysis

Gene annotation, biological pathway, and gene ontology enrichment analysis were conducted using the Metascape webtool⁵³ (www.metascape.org). The input species were specified as Homo sapiens, and the standard parameters were employed. This analytical process centered around the gene sets significantly linked to the triple gradients. Gene sets that exhibited a statistically significant overrepresentation among the members of the input gene list were identified. These highlighted gene sets were considered potential sources of neurobiological insights into hippocampal functional gradients. Additionally, to determine whether the identified gene set demonstrated enrichment within the hippocampus during the specified developmental phase, we conducted a developmental expression analysis. The genes pinpointed through the preceding analysis were cross-referenced with developmental expression profiles sourced from the BrainSpan Atlas of the Developing Human Brain (http://www.brainspan.org/), which contains transcriptomic data from 29 developing brain samples spanning from 8 post-conceptual weeks (PCW) to 40 years of age. The SEA (Specific Expression Analysis) tool, developed by the Dougherty Lab⁵⁴, was employed to perform comparative quantitative analysis, identifying candidate genes that were enriched in specific developmental stages or brain regions across multiple profiles from the given gene list. The significance of the overlap between the identified gene list and lists of transcripts enriched in a particular developmental stage or brain region was then quantified using Fisher’s exact test with Benjamini-Hochberg correction.

Statistical testing of spatial alignment

When testing the statistical significance of the spatial alignment between different hippocampal or cortical maps, it is inappropriate to employ a standard correlation significance test due to the potential bias introduced by spatial auto-correlation within MRI data^110,111. In this context, we employed the variogram matching approach introduced in Burt et al.¹¹². This methodology generates new surrogate hippocampal or cortical maps with spatial autocorrelation resembling that of the input data. The algorithm conceptually involves two main steps: first, the values within a target hippocampal or cortical map are randomly permuted, and then the permuted map undergoes smoothing and rescaling to recover the original spatial autocorrelation structure. Specifically, the smoothing procedures were executed using a distance-dependent kernel, enhancing the similarity of brain features in closely situated regions. Following the application of smoothing to the permuted map, the brain map underwent rescaling to align with the variogram of the target brain map. The variogram, serving as a concise assessment of autocorrelation within spatial data, was used to quantify pairwise feature variation relative to distance. To implement this method, we computed the geodesic distances between all vertex pairs on the hippocampal mid-thickness surface, or between pairs of regions on the cortical mid-thickness surface. Then we utilized the BrainSpace tool¹⁰⁶ to generate surrogate brain maps. In all instances of statistical testing for spatial alignment within our study, we produced 1000 surrogate brain maps and then refitted the model to generate the null distribution.

Analysis of developmental effects

All developmental effects on HCP-D dataset were studied using GAMs with penalized regression splines to capture nonlinear trends in the data while mitigating the risk of complex over-fitting estimates. This model was implemented using the mgcv package (v1.8-42) in R (version 4.2.3). The GAM is a generalized linear model (GLM) in which the linear predictor is a combination of smooth functions of covariates and conventional linear terms. Our specific model configuration incorporated a smooth term for age as a fixed effect. Linear covariates included sex, scanning site, and in-scanner head motion (mean framewise displacement, mFD). Furthermore, the model incorporated the hemisphere as a random effect. To maintain a balanced trade-off between flexibility and computational efficiency, four basis-functions were selected as the upper permissible limit for flexibility within the age smooth term across all models. The optimization process for smoothing parameters was executed through the utilization of restricted maximum likelihood. To quantify the age effect, we computed the change in adjusted R²(Δ Adjusted R², or ΔAdj. R²) between the full GAM model and the reduced model excluding the age term. Finally, the confidence band displayed in all figures of developmental trajectories is associated with a 95% confidence level. Noted that developmental effects in the CBD and NKI-RS datasets were analyzed using generalized additive mixed models (GAMMs), which incorporated both individual participants and hemisphere as random effects.

Data availability

Segmentations and mid-thickness surfaces of the hippocampus, the sampled BOLD time series on the generated mid-thickness surfaces, and cortical parcels for the HCP-D, CBD, and NKI-RS datasets are available upon request. The HCP-D data, after meeting eligibility requirements, are accessible at https://humanconnectome.org/study/hcp-lifespan-development. The NKI-RS data are accessible at https://fcon_1000.projects.nitrc.org/indi/pro/nki.html. All codes are available at https://github.com/debinz/hippocampus_cortex_gradient_youth.

Supplementary material

Supplementary Tables

Supplementary Table 1.

Supplementary Table 2.

Supplementary Table 3.

Supplementary Table 4.

Extended Data

Extended Data Fig. 1.

Extended Data Fig. 2.

Extended Data Fig. 3.

Extended Data Fig. 4.

Extended Data Fig. 5.

Extended Data Fig. 6.

Extended Data Fig. 7.

Extended Data Table 1.

Supplementary Figures

Supplementary Fig. 1.

Supplementary Fig. 2.

Supplementary Fig. 3.

Supplementary Fig. 4.

Acknowledgements

The work is supported by the STI 2030 - the major projects of the Brain Science and Brain-Inspired Intelligence Technology (2021ZD0200500). Shuyu Li is supported by NSFC (32271146) and the Startup Funds for Top-notch Talents at Beijing Normal University. Yong He is supported by NSFC (82021004). Qiongling Li is supported by NSFC (82202245). Qi Dong is supported by NSFC (31521063). Sha Tao, Yanpei Wang, Daoyang Wang, Mingming Hu, and Zhiying Pan are supported by NSFC (31521063) & Beijing Municipal Science & Technology Commission (Z15110000391512). Xi-Nian Zuo receives the Start-up Funds for Leading Talents at Beijing Normal University. Debin Zeng is supported by the China Postdoctoral Science Foundation (2025M772874).

Additional information

Author contributions

Shuyu L., X.N. Z., Yong H., and D. L. conceptualized the work. Shuyu L. and X.N. Z. supervised the study. D. Z. contributed to hippocampal segmentation and mid-thickness surface generation. Yirong H. and X. D. conducted the quality control of segmentations and surfaces of the hippocampus for HCP-D data. Shaoxian L. and S. B. conducted this quality control for CBD data. Z. Y. and Y. Z. conducted this quality control for NKI-RS data. D. Z. and Q. L. performed MR image preprocessing for CBD data. D. Z. and Q. L. conducted structural and functional feature mapping and data analysis. D. Z., X.N. Z., Shuyu L., and Q. L. conducted data interpretation. D. Z. and Q. L. performed data visualization. L. S., and D. Z. performed quality control of MR images for HCP-D data. T. X., L. S., and D. Z. conducted this quality control for NKI-RS data. W. M., Y. W., Shuping. T., J. G., S. Q., Sha T., Q. D. and Yong H. contributed to MR images acquisition of CBD data. Sha T., D. W., Y. W., M. H., and Z. P. ran the cohort and conducted all arrangements of CBD data. Y. X. and L. S. organized these data. X. L., T. Z., X. C. and T. L. performed the quality control of MR images for CBD data. D. Z., X.N. Z., Yong H., Shuyu L., Q. L. and Yirong H. wrote the paper.

Funding

Ministry of Science and Technology of the People's Republic of China (MOST) (2021ZD0200500)

• Qi Dong

MOST | National Natural Science Foundation of China (NSFC) (32271146)

• Shuyu Li

Beijing Normal University (BNU) (the Startup Funds for Top-notch Talents)

• Shuyu Li

MOST | National Natural Science Foundation of China (NSFC) (82021004)

• Yong He

MOST | National Natural Science Foundation of China (NSFC) (31521063)

• Qiongling Li

• Yanpei Wang

• Qi Dong

• Mingming Hu

• Daoyang Wang

• Sha Tao

• Zhiying Pan

Beijing Municipal Science and Technology Commission, Adminitrative Commission of Zhongguancun Science Park (北京市科学技术委员会) (Z15110000391512)

• Sha Tao

• Mingming Hu

• Zhiying Pan

• Yanpei Wang

• Daoyang Wang

Beijing Normal University (BNU) (Start-up Funds for Leading Talents)

• Xi-Nian Zuo

China Postdoctoral Science Foundation (中国博士后科学基金会职) (2025M772874)

• Debin Zeng

Additional files

Supplementary Data 1.

Supplementary Data 2.

Supplementary Data 3.