Participants and clinical assessment

This study protocol was approved by the ethics committees of our institution (No.2021082) and registered at https://clinicaltrials.gov/ (Identifier: ChiCTR2100048109). Written informed consent was obtained from all participants’ guardians. From November 2020 to June 2023, 68 pediatric GHD were prospectively recruited. Meanwhile, 45 typically developing controls (TDs) matched for age and gender was recruited. Detailed clinical and demographic data of the participants are shown in Table 1. Finally, 61 pediatric GHD were enrolled, and 7 patients were excluded due to image artifacts (n = 3) and significant registration errors (n = 4), three TDs were excluded due to image artifacts (n = 1) and significant registration errors (n = 2) (Fig. 1A).

Table 1 Demographic and clinical characteristics of children with GHD and HCsFig. 1

Flowchart of selection of GHD group (A) and the primary analytical process of gray matter MBNs (B) in the current study. (i) 3D T1-weighted imaging and (ii) preprocessing (segment, normalize, modulate, and smooth); (iii) nodes are defined according to the automated anatomical labelling (AAL-116) atlas; (iv) edges are defined according to the combined Euclidean distance method; (v) an individual similarity matrix is obtained; (vi-vii) network properties are calculated and analyzed

For children with GHD, age, gender, height, weight, body mass index (BMI), serum IGF-1, adrenocorticotropic hormone (ACTH), cortisol, and thyroid stimulating hormone (TSH) were obtained from the medical records. Furthermore, GHD children underwent two provocation tests. Blood samples were collected at time 0 and after 30, 60, 90, 120 minutes after intravenous bolus injection of pyridostigmine combined with levodopa. The GH peak was recorded after the provocation tests. The Achenbach's child behavior cheek list (CBCL) was also assigned. Regarding TDs, age, gender, height, weight, and BMI were also recorded.

The inclusion criteria of pediatric GHD: 1) short stature, less than the third percentile or below 2 standard deviations of mean age-matched population height; 2) less than 10 μg/L peak serum GH level with at least two provocative stimulations; 3) no adrenocorticotropic hormone deficiency, hypoglycemia, thyroid-related diseases, and familial genetic and metabolic diseases; 4) right-handedness. The exclusion criteria of each participant were: 1) combined with other mental disorders, personality disorders, or psychotropic drug dependence; 2) challenging to cooperate during MRI examination, and the image quality is too poor for image analysis; 3) a history of other brains organic and metabolic diseases; 4) other MRI contraindications.

Image acquisition and preprocessing

All participants underwent sagittal three-dimensional T1 imaging with a 3.0-T MR imaging system (SIGNA Pioneer GE Healthcare, WI, USA) using a 32-channel phased-array head coil. The head was stabilized with cushions and earplugs. Images were acquired using the fast spoiled gradient recalled echo (FSPGR) sequence, with the following parameters: repetition time (TR) = 8.6 ms, echo time (TE) = 3.3 ms, flip angle = 12°, 188 sagittal slices with slice thickness = 1 mm with no slice gap, a field of view= 256 × 256 mm2, and data matrix = 256 × 256.

Automated segmentation of the whole brain based on 3D T1-weighted images was processed with the CAT12 toolbox (http://www.neuro.uni-jena.de/cat/) within the SPM12 environment (http://www.fil.ion.ucl.ac.uk/spm/software/spm12/) running under MATLAB R2019b (MathWorks). The preprocessing steps involved spatial normalization to the Montreal Neurologic Institute (MNI) space and segmentation. Modulated GM images were resliced to a 2 mm isotropic voxel size and spatially smoothed using a 3D Gaussian kernel with an FWHM of 6 mm, which was chosen in line with previous studies [24, 21]. T1-weighted were downsampled from the raw 1 mm to 2 mm isotropic voxel size, consistent with previous individual-level morphological brain network studies [21, 24]. This resolution balances anatomical detail preservation with computational efficiency for whole-brain network analysis.

Construction of individual morphological similarity networks

The Anatomical Automatic Labeling (AAL-116) atlas [25] was applied to define network nodes or brain regions, and each hemisphere was divided into 45 anatomical regions of interest (ROIs) and cerebellum was divided into 26 anatomical ROIs. Next, the approach named Multivariate Euclidean Distances (MEDs) [23] was performed to estimate the inter-regional morphological similarities between each of the 6670 pairs of the 116 cortical, subcortical and cerebellum regions derived from each individual GMV Map. We calculated Euclidean Distances based on gray matter volumetric values rather than spatial coordinates. This approach assumes that regions with similar gray matter volumes across subjects are more likely to be structurally connected. The details of MEDs were previously described [23], and we outline the key aspects of this algorithm here for clarity and completeness:

For the \(k\) th subject, each pair of anatomical regions \((X,Y)\) from the AAL template was computed using the combined Euclidean distance \(_(X,Y)\), defined as follows:

$$e_k\left(X,Y\right)=\frac\left(\frac2\sum\nolimits_^\sum\nolimits_^_2-\frac1\sum\nolimits_^\sum\nolimits_^_2-\frac1\sum\nolimits_^\sum\nolimits_^_2\right)$$

(1)

Here \(X=\_, \dots , _\}\) and \(Y=\_, \dots , _\}\), where \(x\) and \(y\) denote vertices in regions \(X\) and \(Y\), respectively. \(_\) and \(_\) are the numbers of vertices in \(X\) and \(Y\). The Euclidean distance is computed by the 2-norm (\(_\)).

The Min-Max normalization was performed to minimize possible bias in different ranges of different subjects. The Min-Max normalization between regions X and Y of the \(k\) th subject is computed as follows:

$$_\left(X,Y\right)=\frac_\left(X,Y\right)-_}_-_}$$

(2)

Where \(_\) and \(_\) are the minimum and maximal value in the combined Euclidean distance of the \(k\) th subject, respectively.

In the last step, to obtain the morphological similarity, the value of \(_\left(X,Y\right)\) should be converted to a similarity measurement using the following equation;

$$_\left(X,Y\right)=\text(-_\left(X,Y\right))$$

(3)

Finally, a 116×116 MBNs of each subject was obtained. The values of the edges range from 0 to 1, and 1 represents identical morphological feature distributions in the two AAL regions. A flowchart of the construction of individual-level grey matter MBNs is presented in Fig. 1B.

Network analysis

Network properties were calculated using GRETNA toolbox (http://www.nitrc.org/projects/gretna/) [26] in MATLAB. To ensure the thresholder networks were estimable with sparse properties and small-world index was > 1.0 [27], the minimum and maximum sparsity values were determined. Then, the threshold range was set as 0.05 < S < 0.40 with an interval of 0.05. At each sparsity level, the topologic profiles of brain networks at both global and nodal levels were calculated. Global network profiles included the clustering coefficient (Cp), characteristic path length (Lp), normalized clustering coefficient (γ), normalized characteristic path length (λ), small-world parameters (σ), global efficiency (Eglob), local efficiency (Eloc), and nodal network topological profiles including nodal efficiency (\(_\)), nodal degree (\(_\)), and nodal betweenness (\(_\)). Considering the network sparsity dependent network characteristics, the area under the curve (AUC) of each global profile (Eglob, Eloc, Cp, Lp, \(\gamma\), \(\lambda\), \(\sigma\)) and nodal profile (\(_\), \(_\), \(_\)) across a range of interested densities were calculated as the summarized scalar for each measure, denoted as \(_^\), \(_^\), \(_^\), \(_^\), \(^\), \(^\), \(^\), \(_^\), \(_^\), \(_^\), respectively.

Network-based statistics analysis

To identify the differences in brain network connectivity between GHD and TD groups, the network-based statistics (NBS) method [28] was also used. The NBS has become increasingly popular in recent years for network-level statistical analyses in neuroimaging studies [21, 24], as it offers greater sensitivity in detecting connected components of altered connectivity compared to traditional mass-univariate approaches. To perform NBS analysis, in current study, we first examined whole-brain networks to identify nodes showing significant between-group differences (P < 0.05, uncorrected) in at least one centrality measure (degree, efficiency, or betweenness), which resulted in a set of suprathreshold connections. Within the set of suprathreshold connections, we identified topologically connected components using a breadth-first search algorithm [28]. The size of each identified component was compared against a null distribution obtained through permutation testing (5000 permutations). This step controlled for multiple comparisons at the component level. Significant components (PFDR < 0.05) were visualized and interpreted as subnetworks showing between-group differences in connectivity [28].

Statistical analysis

SPSS v21.0 (IBM Corp., Armonk, New York) was used to perform statistical analysis. The Shapiro-Wilk test was used to evaluate the normality of the data for continuous variables. Qualitative variables were compared by Chi-squared tests, and quantitative variables were compared using two-tailed independent-sample t-tests. A P < 0.05 was set as statistically significant.

Analysis of covariance (ANCOVA) was used to compare between-group differences of the AUC of each network metric (including global and nodal metrics) with diagnosis as fixed factors, age, and gender added to the model as covariates, respectively. The Benjamini-Hochberg false discovery rate (BHFDR) correction was applied to multiple comparisons. Finally, partial correlation analysis was used to examine relationships between significant network metrics and clinical variables, controlling for age and gender as confounding variables (P < 0.05).

Reproducibility analyses

Similar network analysis was repeated for reproducibility analysis with an additional Harvard Oxford atlas with 112 brain regions (HOA-112 atlas, Table S1) [29] to evaluate the potential effects of different parcellation schemes.