Abstract

Introduction:

The heritability of white matter (WM) has been a central focus of diffusion neuroimaging and genetic research. Twin studies have reported high heritability estimates for WM integrity, primarily based on fractional anisotropy (FA). However, other diffusion tensor imaging (DTI) metrics—mean diffusivity (MD), axial diffusivity (AD), and radial diffusivity (RD)—remain less explored.

Methods:

In the present twin study, we assessed WM heritability using FA, MD, AD, and RD metrics. Tract-Based Spatial Statistics was applied to DTI data from 81 healthy twin pairs (33 monozygotic, 48 dizygotic) recruited through the population-based Italian Twin Registry. Twin correlations and genetic and environmental variance components were estimated for each DTI index across brain regions.

Results:

Monozygotic twins exhibited higher correlations than dizygotic twins across most brain regions; however, a few regions showed an unexpected inverse pattern. Notably, while several regions demonstrated strong genetic influences, others showed no detectable genetic contribution, suggesting a substantial role for environmental factors in shaping WM characteristics in those areas.

Conclusion:

These findings revealed regional variability in WM heritability and challenged the assumption of uniform genetic influence, highlighting the importance of considering early environmental factors and supporting the development of more nuanced models of WM development.

1 Introduction

White matter (WM) is the main connective tissue which facilitates efficient communication between disparate regions of the human brain, underpinning a myriad of cognitive and behavioral functions (Bassett et al., 2011). Advancements in neuroimaging, particularly diffusion tensor imaging (DTI), have enabled researchers to non-invasively probe the microstructural integrity of WM tracts, offering insights into the neural substrates of various neuropsychiatric conditions and cognitive processes (Tournier et al., 2011). Among the DTI-derived metrics, fractional anisotropy (FA)—the degree of anisotropy within axonal fibers—has been extensively utilized to assess WM integrity, with studies consistently reporting moderate to high heritability estimates for FA across different age groups. For instance, research leveraging data from the Human Connectome Project and the ENIGMA-DTI consortium has demonstrated heritability estimates for FA ranging from 53 to 90%, underscoring a substantial genetic influence on WM microstructure (Kochunov et al., 2015).

In the last two decades, DTI has been employed to investigate the genetic and environmental influences on WM tissues in twin populations (Blokland et al., 2012). Twin studies indeed provide an excellent tool for analyzing the intricate relationship between genetic and environmental influences on traits of interest. By comparing monozygotic (MZ) and dizygotic (DZ) twins, researchers can estimate the heritability of various brain traits, including susceptibility to psychiatric disorders (Brambilla et al., 2014; Colli et al., 2024), the influence of genetic factors on behavior (Maggioni et al., 2024), and patterns of functional connectivity activation (Tassi et al., 2023). From this perspective, investigating WM heritability in twins provides valuable insights into the interplay between genetic and environmental influences on brain connectivity. This interaction underpins individual differences in cognition and behavior, as well as susceptibility to psychopathological conditions (Fagnani and Stazi, 2020; Maggioni et al., 2020). Twin studies have consistently demonstrated that DTI metrics are significantly influenced by genetic factors, underscoring their fundamental role in WM development. In neonatal populations, research has reported substantial heritability estimates for FA, radial diffusivity (RD) and axial diffusivity (AD)—which reflects diffusion perpendicular and parallel to axonal fibers, respectively—across several WM tracts, including fornix, cingulum and corpus callosum (Lee et al., 2015). These findings suggest that genetic factors play a critical role in the early development of WM microstructure, as supported by previous evidence (Gao et al., 2019), which has highlighted how the genetic architecture of WM is complex and highly polygenic (Arnatkeviciute et al., 2021; Zhao et al., 2021a). It is not surprising that large-scale genome-wide association studies (GWAS) have identified numerous genetic loci associated with variations in WM microstructure (Sha et al., 2023). For example, a study analyzing data from over 17,000 UK Biobank participants found that DTI parameters are substantially heritable across all WM tracts, with a mean heritability of 48.7%. Moreover, the study identified 213 independent significant single-nucleotide polymorphisms (SNPs) associated with 90 DTI parameters, many of which have been implicated in cognitive and mental health traits (Zhao et al., 2021b).

However, most studies investigating WM heritability in twin populations present several limitations that may affect the generalizability and interpretation of their findings (Videtta et al., 2024). First, research samples are often heterogeneous, frequently including unpaired twins or siblings, which can confound heritability estimates. Second, the reliance on specific datasets may reduce the representativeness of study populations. For example, commonly used datasets such as Vietnam Era Twin Study of Aging (VETSA) (Kremen et al., 2013) and Older Australian Twins Study (OATS) (Sachdev et al., 2009) are drawn from highly specific cohorts, including war veterans and older twin adults, respectively. Such sampling constraints may also introduce gender imbalances; notably, the VETSA dataset comprises only male twin pairs (Kremen et al., 2013). Additionally, while FA is widely used as a measure of WM integrity, other DTI metrics—such as mean diffusivity (MD), which reflects average water diffusion within axonal fibers, as well as AD and RD—remain underexplored in the context of heritability, despite offering complementary insights into WM microstructure.

In light of these limitations, the present study aims to estimate WM heritability across all major DTI indices (FA, MD, AD, RD) in a representative sample of healthy Italian twins. To strengthen current evidence, we combined Tract-Based Spatial Statistics (TBSS) with classical twin modeling, focusing on data from paired twin samples to ensure methodological consistency and demographic representativeness.

2 Method2.1 Participants

A total sample of 121 pairs (242 twins) was selected from the twin cohorts participating in the “WhyMe?” and in the “SPES” Projects (Italian Ministry of Health, grants n. RF-2011-02352308, and n. GR-2010-2316745, respectively). Twins were recruited through the population-based Italian Twin Registry (Medda et al., 2019) in two different centers, Scientific Institute IRCCS Eugenio Medea of Bosisio Parini (Lecco, Italy) and Scientific Institute IRCCS Eugenio Medea of Pasian di Prato (Udine, Italy). Inclusion criteria were: (i) healthy twins; (ii) availability of Magnetic Resonance Imaging (MRI) scan; (iii) quality of DTI images acceptable for analyses. Unpaired twins and those pairs in which only one member met inclusion criteria were excluded. From 121 pairs (242 twins), 6 pairs (12 twins) were excluded due to psychiatric disorders in one or both twins and 19 pairs (38 twins) due to missing MRI data. An additional 15 pairs (30 twins) were ruled out due to low-quality DTI data, leading to a final sample of 81 pairs (162 twins): 47 pairs from the Scientific Institute IRCCS Eugenio Medea of Bosisio Parini (Lecco, Italy), and 34 pairs from the Scientific Institute IRCCS Eugenio Medea of Pasian di Prato (Udine, Italy).

The study was approved by the local ethical committee according to the Declaration of Helsinki and its subsequent versions (World Medical Association, 2013). Participants provided a written informed consent. For minors (<18 years), an assent was obtained, and a legal tutor provided informed consent.

2.2 MRI acquisition

Each participant underwent an MRI acquisition either at Scientific Institute IRCCS Eugenio Medea of Bosisio Parini (Lecco, Italy) or Scientific Institute IRCCS Eugenio Medea of Pasian di Prato (Udine, Italy). Two different MRI scanners were employed: Philips Achieva D-Stream 3T and Philips Achieva 3T. Acquisition protocols were identical across scanners. Diffusion images were acquired axially with anterior–posterior direction and one no-weighted image was acquired at b = 0 s/mm2: spin-echo sequencing, 64 directions, b = 1,000 s/mm2, echo time (TE) = 76 ms, repetition time (TR) = 8,895 ms, voxel size = 1.5 × 1.5 × 1.5, number of slices = 65, matrix = 138 × 138, field of view (FOV) = 207 × 207, flip angle = 90 degrees.

2.3 Tract-based spatial statistics

Diffusion-weighted images were preprocessed using tools from the FMRIB Software Library (FSL; Jenkinson et al., 2012). For each participant, diffusion data were corrected for eddy current-induced distortions and simple head motion using FSL’s eddy_correct, applying an affine registration of all diffusion-weighted volumes to the first b0 image. Non-brain tissue was removed using the Brain Extraction Tool (BET) with a fractional intensity threshold of 0.1, and the resulting brain mask was used for subsequent tensor fitting. Because the diffusion preprocessing was performed using eddy_correct, which does not include slice-wise motion correction or outlier replacement, residual motion effects, particularly relevant in pediatric and adolescent samples, could not be fully excluded. A diffusion tensor model was fitted at each voxel to compute the maps of diffusion scalar metrics: FA, MD, AD. Maps of RD were computed using FSL’s fslmaths by dividing the sum of the second (L2) and third (L3) eigenvalue maps by 2. Participants’ DTI maps were processed using Tract-Based Spatial Statistics (Smith et al., 2006). In detail, nonlinear registration was applied to align the images to the FMRIB58 FA template. Then, the registered images were averaged and thresholded at 0.2 to create a mean FA skeleton representing the central trajectories of WM tracts. Each participant’s DTI data were then projected onto this skeleton.

Regional diffusion metrics were extracted from the TBSS skeleton using a custom Python script based on Nibabel and NumPy (see Supplementary materials). Skeletonised diffusion maps (FA, MD, AD, RD) generated by TBSS were stored as 4D images in standard space, with the fourth dimension indexing participants. The Johns Hopkins University (JHU) ICBM-DTI-81 white matter atlas was used to define regions of interest (Mori et al., 2008). For each atlas-defined region, voxel indices corresponding to the region label were identified. For each participant, diffusion values were extracted from the skeletonised maps at those voxel locations. To ensure that only voxels belonging to the TBSS skeleton were included, zero-valued voxels (i.e., voxels outside the skeleton projection) were excluded from the computation. Mean diffusion values were then calculated across the remaining voxels within each region, yielding one regional value per participant per metric. Regions defined bilaterally in the atlas were analysed separately for left and right hemispheres. Six regions (middle cerebellar peduncle, pontine crossing tract, bilateral medial lemniscus, and bilateral inferior cerebellar peduncle) were excluded due to systematic quality issues and poor skeleton representation. Diffusion data quality was assessed following preprocessing and TBSS projection through visual inspection of corrected diffusion volumes, tensor-derived maps, and skeletonised images. Datasets showing severe motion artefacts, signal dropout, or excessive noise that compromised reliable tensor estimation or skeleton projection were excluded. Based on these criteria, 15 twin pairs were removed from the analysis. Exclusion rates did not differ between monozygotic and dizygotic twins, reducing the likelihood of zygosity-related bias.

Finally, ComBat harmonization (Fortin et al., 2017) was performed at the regional level on diffusion values extracted from the TBSS skeleton, rather than at the voxelwise level, to mitigate scanner-related site effects. Scanner model was specified as the batch variable, while age and sex were included as biological covariates to be preserved, ensuring that developmental and sex-related variability in WM microstructure was retained.

Figure 1 reports the schematic workflow of the diffusion MRI preprocessing and analysis pipeline.

Overview of the diffusion MRI preprocessing and analysis pipeline. Diffusion-weighted images were preprocessed using FSL, followed by TBSS-based skeletonization. Regional diffusion metrics were extracted from the skeleton using the JHU ICBM-DTI-81 atlas, harmonized across sites using ComBat (preserving age and sex), and analyzed using twin modeling.

2.4 Twin modeling

Twin modeling analyses were performed with the Mx Software (Neale et al., 2006). Univariate biometric twin models incorporating additive genetic (A), either shared environmental (C) or non-additive genetic (D), and unshared (individual-specific) environmental (E) components of variance (Neale and Cardon, 1992) were applied to DTI metrics for the 42 WM regions by using sex- and age-adjusted correlations as input data for MZ and DZ groups. The A component represents the additive effects of all gene variants (i.e., alleles) that influence the trait, without interactive effects; the D component represents interactions between alleles at the same fixed chromosomal site (i.e., locus) (“dominance”) or at different loci (“epistasis”); the C component represents the effects of environmental factors that are shared by the twins within the family (e.g., rearing environment, family socio-economic status, parental behaviors) or in the womb (e.g., hormonal exposures); the E component represents the effects of environmental factors that are unique to an individual (e.g., lifestyles, relations with peers, infections), including measurement error. Only three of these four variance components can be simultaneously estimated because C and D components are confounded when analyzing data from MZ and DZ twins reared together. Given the limited power of the study, it was deemed prudent to report the estimated genetic and environmental components of variance for each of the regions and DTI indices under the fully parameterized (i.e., either ACE or ADE) models only, without looking for more parsimonious sub-models, which might penalize magnitude-relevant, yet not significant components. For the same reason, only point estimates of correlations in MZ and DZ pairs and of variance components were interpreted, because large standard errors often resulted in non-significant p-values, even in the case of sizeable correlation differences and variance components. The genetic component (G), referred to within the text and depicted in Figure 2, indicates the additive (“narrow-sense,” A) heritability under the ACE model, or the additive plus non-additive (“broad-sense,” A + D) heritability under the ADE model.

Genetic and environmental variance components for each diffusion tensor imaging metrics. (A) FA genetic and environmental variance components. (B) MD genetic and environmental variance components. (C) AD genetic and environmental variance components. (D) RD genetic and environmental variance components. C, Shared environment; G, Heritability [equal to additive (“narrow-sense”, A) heritability under the ACE model, or to additive plus non-additive (“broad-sense”, A + D) heritability under the ADE model]; E, unshared environment; gCC, Genu of Corpus Callosum; bCC, Body of Corpus Callosum; sCC, Splenium of Corpus Callosum; FX_R, Fornix Right; FX_L, Fornix Left; CST_R, Corticospinal Tract Right; CST_L, Corticospinal Tract Left; SCP_R, Superior Cerebral Peduncle Right; SCP_L, Superior Cerebral Peduncle Left; CP_R, Cerebral Peduncle Right; CP_L, Cerebral Peduncle Left; ALIC_R, Anterior Limb of Internal Capsule Right; ALIC_L, Anterior Limb of Internal Capsule Left; PLIC_R, Posterior Limb of Internal Capsule Right; PLIC_L, Posterior Limb of Internal Capsule Left; RIC_R, Retrolenticular part of Internal Capsule; RIC_L, Retrolenticular part of Internal Capsule Left; ACR_R, Anterior Corona Radiata Right; ACR_L, Anterior Corona Radiata Left; SCR_R, Superior Corona Radiata Right; SCR_L, Superior Corona Radiata Left; PCR_R, Posterior Corona Radiata Right; PCR_L, Posterior Corona Radiata Left; PTR_R, Posterior Thalamic Radiation Right; PTR_L, Posterior Thalamic Radiation Left; SS_R, Sagittal Stratum Right; SS_L, Sagittal Stratum Left; EC_R, External Capsule Right; EC_L, External Capsule Left; CgC_R, Cingulum (Cingulate gyrus) Right; CgC_L, Cingulum (Cingulate gyrus) Left; CgH_R, Cingulum (Hippocampus) Right; CgH_L, Cingulum (Hippocampus) Left; FX/ST_R, Fornix (cres)/Stria terminalis Right; FX/ST_L, Fornix (cres)/Stria terminalis Left; SLF_R, Superior Longitudinal Fasciculus Right; SLF_L, Superior Longitudinal Fasciculus Left; SFOF_R, Superior Fronto-Occipital Fasciculus Right; SFOF_L, Superior Fronto-Occipital Fasciculus Left; UF_R, Uncinate Fasciculus Right; UF_L, Uncinate Fasciculus Left; TAP_R, Tapetum Right; TAP_L, Tapetum Left.

2.5 Statistical analysis

Statistical significance for gender and age differences was assessed using a Python script. A chi-squared (χ2) test was applied for gender, and an unpaired t-test for age. A significance threshold of p < 0.05 was used for both tests.

3 Results3.1 Sociodemographic data

The sample consisted of 81 pairs (162 twins), including 78 males and 84 females, with a mean age of 16.7 ± 6.1 years (range: 9–32 years) (Supplementary Table S1). Of these, 33 were MZ pairs (18 males, 15 females; mean age: 16.8 ± 6.4 years) and 48 were DZ pairs (21 males, 27 females; mean age: 13.3 ± 2.5 years), of which 38 were same-sex pairs and 10 were opposite-sex pairs.

Comparing MZ and DZ pairs, no differences were found in gender (χ2 = 3.9, p = 0.27) and age (t = 0.14, p = 0.88) (Supplementary Table S1).

3.2 Twin correlations

Twin correlations were performed on MZ pairs and DZ pairs. Substantial correlations were observed across all DTI metrics, with MZ twins generally showing higher correlations than DZ twins in most brain regions (Table 1). Specifically, the mean correlations for MZ and DZ pairs, respectively, were: (i) FA: 0.71 (range: −0.07 to 0.91) and 0.49 (range: 0.07–0.65); (ii) MD: 0.58 (range: 0.15–0.89) and 0.37 (range: 0.01–0.57); (iii) AD: 0.56 (range: 0.08–0.85) and 0.39 (range: 0.06–0.57); and (iv) RD: 0.65 (range: 0.30–0.91) and 0.42 (range: 0.05–0.57).

RegionTwins correlationsGenetic and environmental variance componentsMZDZp-value*ModelACDA + DGETotalFAgCC0.810.370.0006ADE0.65–0.170.810.810.191.00bCC0.850.610.0126ACE0.470.38––0.470.151.00sCC0.890.540.0002ACE0.700.19––0.700.111.00FX0.440.240.1578ACE0.420.03––0.420.561.00CST_R−0.020.070.3386ACE0.000.03––0.000.971.00CST_L0.220.120.3262ACE0.210.01––0.210.781.00SCP_R0.220.320.3201ACE0.000.28––0.000.721.00SCP_L0.410.160.1193ADE0.22–0.190.410.410.591.00CP_R0.730.600.1494ACE0.270.47––0.270.271.00CP_L0.780.500.0168ACE0.560.22––0.560.221.00ALIC_R0.810.660.0696ACE0.310.50––0.310.191.00ALIC_L0.780.660.1429ACE0.240.54––0.240.221.00PLIC_R0.880.610.0018ACE0.550.34––0.550.121.00PLIC_L0.860.570.0026ACE0.580.28––0.580.141.00RIC_R0.740.490.0417ACE0.490.24––0.490.261.00RIC_L0.750.450.0191ACE0.600.15––0.600.251.00ACR_R0.830.400.0006ADE0.75–0.080.830.830.171.00ACR_L0.750.300.0022ADE0.45–0.300.750.750.251.00SCR_R0.780.520.0202ACE0.540.25––0.540.221.00SCR_L0.850.460.0008ACE0.770.08––0.770.151.00PCR_R0.870.520.0007ACE0.690.18––0.690.131.00PCR_L0.790.560.0268ACE0.480.32––0.480.211.00PTR_R0.790.440.0059ACE0.700.09––0.700.211.00PTR_L0.820.580.0194ACE0.470.35––0.470.181.00SS_R0.760.530.0439ACE0.460.30––0.460.241.00SS_L0.720.510.0705ACE0.420.30––0.420.281.00EC_R0.680.660.4346ACE0.040.64––0.040.321.00EC_L0.880.570.0011ACE0.620.26––0.620.121.00CgC_R0.600.640.3861ACE0.000.63––0.000.371.00CgC_L0.710.600.2028ACE0.220.49––0.220.291.00CgH_R0.550.450.2866ACE0.200.35––0.200.451.00CgH_L0.650.370.0495ACE0.560.09––0.560.351.00FX/ST_R0.710.560.1438ACE0.300.41––0.300.291.00FX/ST_L0.650.620.4339ACE0.050.60––0.050.351.00SLF_R0.770.560.0513ACE0.420.35––0.420.231.00SLF_L0.790.640.0941ACE0.300.49––0.300.211.00SFOF_R0.800.650.0882ACE0.290.51––0.290.201.00SFOF_L0.810.510.0099