A Comparison of Heritability Maps of Cortical Surface Area and Thickness and the Influence of Adjustment for Whole Brain Measures: A Magnetic Resonance Imaging Twin Study

Abstract

Understanding the genetic and environmental contributions to measures of brain structure such as surface area and cortical thickness is important for a better understanding of the nature of brain-behavior relationships and changes due to development or disease. Continuous spatial maps of genetic influences on these structural features can contribute to our understanding of regional patterns of heritability, since it remains to be seen whether genetic contributions to brain structure respect the boundaries of any traditional parcellation approaches. Using data from magnetic resonance imaging scans collected on a large sample of monozygotic and dizygotic twins in the Vietnam Era Twin Study of Aging, we created maps of the heritability of areal expansion (a vertex-based area measure) and cortical thickness and examined the degree to which these maps were affected by adjustment for total surface area and mean cortical thickness. We also compared the approach of estimating regional heritability based on the average heritability of vertices within the region to the more traditional region-of-interest (ROI)-based approach. The results suggested high heritability across the cortex for areal expansion and, to a slightly lesser degree, for cortical thickness. There was a great deal of genetic overlap between global and regional measures for surface area, so maps of region-specific genetic influences on surface area revealed more modest heritabilities. There was greater inter-regional variability in heritabilities when calculated using the traditional ROI-based approach compared to summarizing vertex-by-vertex heritabilities within regions. Discrepancies between the approaches were greatest in small regions and tended to be larger for surface area than for cortical thickness measures. Implications regarding brain phenotypes for future genetic association studies are discussed.

Individual differences in cortical surface area are likely to reflect important developmental processes (Dickerson et al., 2009; Ostby et al., 2009) and appear to be strongly related to behavioral measures (Koscik et al., 2009; Schwarzkopf et al., 2011). Several recent studies have begun to address the genetic underpinnings of human variation in regional surface area measurements. Our group found high heritabilities for total surface area (Kremen et al., 2010) and moderate regional surface area heritabilities (Eyler et al., 2011). The degree of genetic contribution to regional surface area was reduced after controlling for total surface area, suggesting that some of the genetic variance in these parcellations is related to genes that influence the global size of the cortical surface (Eyler et al., 2011). Genetic influences on surface area were pronounced in frontal and parietal cortex, whereas surface area in medial temporal lobe regions appeared to have stronger environmental than genetic influences, particularly after adjustment for total area (Eyler et al., 2011).

Cortical thickness is another structural measure that has demonstrated sensitivity to disease, changes over the course of development, and relationships to cognitive performance (Dickerson et al., 2009; Fjell et al., 2009; Karama et al., 2011). Others and we have found that individual differences in cortical thickness were under substantial genetic control (Joshi et al., 2011; Kremen et al., 2010; Rimol et al., 2010; Schmitt et al., 2008). Our studies have also shown that the genetic influences on surface area were distinct from the genetic influences for cortical thickness: the genetic correlation between surface area and thickness measures was near zero (Panizzon et al., 2009). In addition, unlike for surface area, controlling for total thickness had only a modest effect on regional heritabilities for thickness measures, suggesting unique genetic influences on regional cortical thickness that are not shared with the genetic factors that determine overall thickness of the cortex (Eyler et al., 2011). Winkler et al (2010) found similarly low genetic correlations between surface area and thickness in a family pedigree study, and also observed regional heritabilities of moderate size, even after adjustment for global measures.

To expand on our previous work and better understand regional variability in the genetic determinants of cortical surface area, we used a regional measure of surface area (i.e., areal expansion; see Methods) that was spatially unconstrained by traditional boundaries based on anatomical features, and estimated genetic and environmental variance components using a large sample of twins. The advantage of this method over those used previously to examine regional heritability of surface area is that it allows us to find patterns of variability that cross the boundaries contained in pre-determined atlas parcellation systems. Specifically, heritability could be inaccurately estimated if multiple genetic or environmental sources contribute to the surface area of a region determined through anatomical landmarks or cellular features, because those regions may or may not form genetically meaningful units (Chen et al., 2011; Rimol et al., 2010). Furthermore, greater measurement error in the determination of surface area for some small cortical regions could lead to underestimates of true regional heritability.

In the current study, our first aim was to map heritability of areal expansion and cortical thickness within the same sample, using identical methods, in order to reveal spatial patterns of the relative influence of genetic versus environmental factors in these measures without the constraints of macroanatomical structural boundaries. Our second aim was to explore the influence of statistical control for global surface area and thickness on heritability maps. This analysis would help to reveal the degree to which the genetic factors that influence overall surface area and thickness are shared by regional measures. Our third aim was to compare two methods of determining regional heritability of surface area and cortical thickness within a sulcal-based parcellation system (Desikan et al., 2006): one approach first estimated heritability at each vertex and then averaged within a region; the second, more traditional approach first determined the surface area and thickness of each region of interest (ROI) and then estimated heritability of those values. These approaches are likely to differ in the degree of variability in heritability estimates between regions because of the influence of spatial averaging on measurement error and, thus, on heritability estimates. It is a well-known psychometric property that measurement error for an overall test score is less than that of any individual item; similarly, in measures of brain features, measurement error decreases with greater summation or averaging of surface features. In the case of vertex-based estimates, the amount of spatial averaging is fixed because it is imposed by the processing method across all vertices on the surface (a process referred to as spatial smoothing). For the ROI-based measures, the degree of averaging or summation varies considerably from region to region, increasing with increasing region size. Because measurement error is included in the unique environmental variance component of the twin model, heritability estimates should increase with increased region size for ROI-based calculations, but not for vertex-based calculations.

Methods

Participants

The Vietnam Era Twin Study of Aging (VETSA) project has been described previously (Kremen et al., 2006). The VETSA sample was drawn from the Vietnam Era Twin (VET) Registry (Goldberg et al., 2002), a sample of male-male twin pairs born between 1939 and 1957 who had both served in the United States military between 1965 and 1975. The study sample is not a VA or patient group; the majority of individuals were not exposed to combat. For this analysis, 474 individual VETSA participants were included. Of those, 406 were paired (i.e., 203 twin pairs): 110 monozygotic (MZ) and 93 dizygotic (DZ) pairs. Zygosity for 92% of the sample was determined by analysis of 25 satellite markers that were obtained from blood samples. For the remainder of the sample, zygosity was determined through a combination of questionnaire and blood group methods (Eisen et al., 1989).

Mean age of the MRI participants was 55.8 years (SD = 2.6, range = 51–59), mean years of education was 13.9 (SD = 2.1), and 85.2% were right-handed. Most participants (74.9%) were employed full-time, 4.2% were employed part-time, and 11.2% were retired. There were 88.3% non-Hispanic white, 5.3% African-American, 3.4% Hispanic, and 3.0% ‘other’ participants. Self-reported overall health status was as follows: excellent (14.8%); very good (36.5%); good (37.4%); fair (10.4%); and poor (0.9%). Demographic characteristics of the VETSA MRI sample did not differ from the entire sample, and are comparable to U.S. census data for similarly aged men (Centers for Disease Control and Prevention, 2003; National Center for Disease Statistics, 2003). There were no significant demographic differences between MZ and DZ twins.

All participants gave informed consent to participate in the research and the study was approved by the institutional review boards of the University of California, San Diego (UCSD), Boston University, and the Massachusetts General Hospital (MGH).

Image Acquisition

Sagittal T1-weighted Magnetization Prepared RApid Gradient Echo (MPRAGE) images (two per case) were acquired on Siemens 1.5 Tesla scanners (241 at UCSD; 233 at MGH). Scan parameters were: TI = 1000 ms, TE = 3.31 ms, TR = 2730 ms, flip angle = 7 degrees, slice thickness = 1.33 mm, voxel size 1.3 × 1.0 × 1.3 mm. Data were reviewed for quality, registered, and averaged to improve signal-to-noise. Of the 493 scans available at the time of these analyses, quality control measures excluded 0.6% (3 cases) due to scanner artifact, and 3% (16 cases) due to inadequate image processing results (e.g., poor contrast caused removal of non-brain to fail). The resultant 474 available cases included 203 twin pairs (406 individuals) that were used in the present study.

Image Processing

As in our previous work (Eyler et al., 2011; Kremen et al., 2010), the cortical surface was reconstructed using methods based on the publicly available FreeSurfer software package (Dale et al., 1999; Desikan et al., 2006; Fischl & Dale, 2000; Fischl et al., 1999; Fischl et al., 2004). Processing began with correction for variation in image intensity due to magnetic field inhomogeneities, creation of a normalized intensity image, and removal of the skull (non-brain). Preliminary segmentation using a connected components algorithm was performed, and interior holes in the components representing white matter were filled, resulting in a single filled volume for each cortical hemisphere. The resulting surface was covered with a polygonal tessellation and smoothed to reduce metric distortions. To obtain a representation of the gray/white boundary, a refinement procedure was applied, and the resulting surface was deformed outwards to obtain an explicit representation of the pial surface. Once generated, the cortical surface model was manually reviewed and edited for technical accuracy. Minimal manual editing was performed by applying standard, objective editing rules. Maps were spatially smoothed using iterative nearest neighbor smoothing with 2,819 iterations.

Each individual's map was placed into a common coordinate system using a non-rigid, high-dimensional, spherical averaging method to align cortical folding patterns (Fischl et al., 1999). This procedure provides accurate matching of morphologically homologous cortical locations across subjects on the basis of each individual's anatomy while minimizing metric distortion. The maps thus produced are not restricted to the voxel resolution of the original images and allow for submillimeter spatial resolution (Kremen et al., 2008). Estimates of cortical area were obtained by computing the area of each triangle in the standardized, spherical atlas space surface tessellation, when mapped into the individual subject space. This provides point-by-point estimates of the relative areal expansion or compression from the individual subject space to the atlas space for each location in atlas space. A standard bivariate twin model was then fitted separately to each of the uniformly-distributed surface locations (vertices). Cortical thickness at each of these vertices was also calculated as distance from the pial surface to the white matter surface along a line that is oriented perpendicular to the local white matter surface.

Statistical Analysis

Based on our previous findings of minimal common environmental influences on surface area (Eyler et al., 2011) and cortical thickness (Kremen et al., 2010; Rimol et al., 2010), we used a twin model that estimated contributions of additive genetic effects (A) and individual-specific environmental effects (E) to the variance in areal expansion or cortical thickness at each vertex. The variance-covariance patterns were examined by fitting models with Mx, a maximum-likelihood-based structural equation modeling program (Neale et al., 2003). We sought to map ‘unadjusted’ genetic and environmental effects at each cortical location, which includes those genetic and environmental effects shared with total surface area, as well as map estimates of ‘adjusted’ genetic and environmental contributions to areal expansion or cortical thickness at a particular location that are unique to this location. To accomplish this, bivariate twin models (using both vertex-based areal expansion and total surface area measures, and both vertex-based and total cortical thickness measures) were used to estimate the genetic and environmental contributions to the total phenotypic variance at each vertex. The unique genetic contributions to areal expansion or cortical thickness at each cortical location (adjusted heritability) were estimated by using the bivariate model to account for genetic covariances with the global measures.

The traditional approach to determining the heritability of a cortical region does not measure vertices individually; rather, the heritability of an ROI's total surface area or average thickness is estimated, and this value is entered into the model for estimating variance components. Here, we compared this traditional approach with regional heritability estimates determined by calculating the heritability of each vertex and then averaging the heritability estimates across all vertices that fall within a region in standard space. We first compared the inter-regional variability between these two approaches using Levene's test for equivalence of variances.

We then explored the relationship between the size of each region and the magnitude of the discrepancy between the heritability estimates from the two approaches. Importantly, the degree of spatial averaging for the vertex-based estimates is constant; the effective resolution that resulted from the applied smoothing kernel was approximately 5,580 vertices. For the ROI-based estimates, the amount of spatial averaging or summation varied considerably between regions, because they ranged in size from 336 vertices (left frontal pole) to 13,062 vertices (left superior frontal gyrus). We expected that regional differences in size would influence measurement error, and thus heritability, for the ROI-based but not the vertex-based approach. The relationship of the ratio of heritability estimates from the two approaches to region size would therefore follow a predictable pattern, such that the ratio would increase as ROI size increases. To examine this pattern, we plotted the relationship between the ratio of ROI-based and vertex-based heritabilities and region size (in number of vertices) and fit the following function: Ratio = 1/(1 + b/ROI\Size) (see Appendix for derivation). The parameter ‘b’ estimates the error variance as a fraction of total variance; higher b means that the curve more quickly reaches a point where the two methods yield similar heritability estimates. This function will fit the data well if the discrepancy in heritability estimates between the two approaches is driven primarily by a direct effect of ROI size on measurement error in the ROI-based approach.

We would also expect a somewhat stronger association between heritability and ROI size for surface area than for cortical thickness, due to the combination of two factors that can vary across regions and affect measurement error: degree of spatial averaging and inaccuracy (or variability) of boundary placement. Although degree of spatial averaging (decreasing with decreasing ROI size) impacts both thickness and surface area measures, inaccuracy of boundary placement (which has a proportionally greater impact on smaller ROIs) poses a particular problem for estimates of surface area because the magnitude of the measured area is directly affected by the placement of the boundary. For cortical thickness (which is averaged across the ROI on an individual basis), as well as for both measures when calculated by averaging across a region in standard space, boundary placement only affects which points are included in the average.

Total surface area was computed by calculating the sum of areal expansion measures across all vertices; mean cortical thickness was calculated by averaging the cortical thickness measure across all vertices. All bivariate models included the effects of site of data collection (MGH or UCSD) and age as fixed effects on the means.

Results

Continuous Maps of Areal Expansion and Cortical Thickness Heritability and Effect of Adjustment for Global Measures

Before adjustment for total surface area, the high-resolution map of heritability of areal expansion at each point on the cortex showed very high estimates with some regional variability (Figure 1). The degree of genetic correlation with total surface area was pronounced in all regions, suggesting a large overlap in the genes that influence total and regional surface area (Figure 1, middle row). Because of this overlap, heritabilities were lower after adjustment for total surface area (Figure 1, bottom row), but considerable region-specific genetic contributions still were evident in some regions.

FIGURE 1 Top Row: Map of heritability for areal expansion at each point on the cortex that is unadjusted for total surface area. Middle Row: Map of the genetic correlation between total surface area and areal expansion at each point on the cortex. Bottom Row: Map of heritability for areal expansion at each point on the cortex adjusted for the genetic contributions shared with total surface area. Color scales are adjusted for each row to encompass the range of mapped values.

Unadjusted maps of the heritability of cortical thickness showed moderate to high genetic contributions to this anatomical measure (Figure 2, top row). The genetic correlations between vertex-based cortical thickness and average thickness across the entire cortex were moderate and somewhat lower than those for areal expansion (Figure 2, middle row). After adjustment for global thickness, regional heritabilities for cortical thickness were reduced somewhat (Figure 2, bottom row).

FIGURE 2 Top Row: Map of heritability for cortical thickness at each point on the cortex that is unadjusted for total cortical thickness. Middle Row: Map of the genetic correlation between total cortical thickness and thickness at each point on the cortex. Bottom Row: Map of heritability for cortical thickness at each point on the cortex adjusted for the genetic contributions shared with total cortical thickness. Color scales are adjusted for each row to encompass the range of mapped values.

Vertex-Based Versus ROI-Based Heritability Estimates

To compare vertex-based heritability results to those calculated using regions of interest that are determined by gyral features, we applied a standard parcellation system (Desikan et al., 2006) in atlas space to the maps and averaged the vertex-wise heritability estimates (unadjusted and adjusted) within each cortical region. We found that there was greater inter-regional variation in heritability estimates when calculated using the ROI-based method compared to the vertex-based method (Table 1) for unadjusted, F(130) = 69.7, p < .001, and adjusted surface area F(130) = 10.2, p = .002, and for unadjusted F(130) = 35.6, p < .001, and adjusted cortical thickness F(130) = 9.7, p = .002. Further, discrepancies between vertex-based and ROI-based heritabilities (as measured by the ratio of the two) were greatest for small regions (Figure 3). The curvilinear associations fit the predicted function based on the influence of spatial averaging on measurement error (surface area R ² = .66; cortical thickness R ² = 0.29). The curves for the adjusted heritabilities reached an asymptote sooner than the unadjusted curves (not shown). Presumably, the adjustments remove the global positive correlations; thus, the spatial averaging associated with increased ROI size can no longer increase heritability beyond some relatively low level. Nevertheless, the magnitude of the relationship to ROI size was greater for surface area than for cortical thickness in both cases. For example, it can be seen from the red lines in Figure 3 that, for a small region (2,000 vertices), the ROI-based approach will yield heritability estimates for surface area that are 43% the size of those found with the vertex-based approach. For cortical thickness, the same small region will yield ROI-based heritabilities that are much closer to the size of the vertex-based estimates (80%).

TABLE 1 Unadjusted and Adjusted Heritability Estimates for Areal Expansion/Surface Area and Cortical Thickness in Regions of the Desikan-Killiany Cortical Parcellation System

Values are calculated by averaging vertex-based estimates and by basing the estimates on measures in each participant's regions of interest (ROI-based).

FIGURE 3 Plot of the ratio of unadjusted ROI (region of interest)-based to vertex-based heritabilities within 132 cortical parcellations against the size of each ROI in number of vertices. Surface area is indicated with solid circles and cortical thickness with open squares. Fit line for the equation Ratio = 1/(1+b/ROI_Size) is shown as a solid line for surface area and a dashed line for cortical thickness. Red lines illustrate the predicted ratio of ROI-based to vertex-based heritabilities for a small region (2000 vertices) for surface area (solid line) and cortical thickness (dashed line).

Discussion

Consistent with our previous findings, the current results showed that genetic influences play a large role in individual differences in both surface area and cortical thickness. Using vertex-based measures of areal expansion, we created heritability maps that revealed only modest regional variation. Unadjusted cortical thickness heritabilities were generally lower than for areal expansion. This is consistent with one previous study that reported heritabilities for both surface area and cortical thickness as measured with sulcal-based regions of interest (Winkler et al., 2010) and with our own study that reported lobar heritabilities for both measures in the same sample (Panizzon et al., 2009).

When examining the relationship of regional to global measures of surface area, we found large genetic correlations between total and regional surface area across the entire cortex. Examination of maps of the region-specific genetic contributions to local areal expansion after accounting for shared genetic variance with total surface area reveal that heritability was moderate overall, and slightly more variable across the cortex than seen in unadjusted maps. These findings of a strong genetic overlap between global and regional surface area measures are similar to what has been observed in previous analyses using sulcal-based regions of interest (Eyler et al., 2011; Winkler et al., 2010). For cortical thickness, the genetic correlations with total mean thickness were weaker, so adjustment for total thickness made less of a difference to the regional heritabilities. This is consistent with our reported lobar findings (Eyler et al., 2011) and similar to findings from a family study (Winkler et al., 2010), although adjustment was for both total surface area and total cortical thickness in that analysis.

The large heritability estimates we observed for both regional surface area and cortical thickness using a vertex-based approach contrasts with the more moderate values reported in previous studies by our group (Eyler et al., 2011; Kremen et al., 2010) and by others (Schmitt et al., 2009; Winkler et al., 2010). Greater heritability for vertex-based results may be due to the degree of spatial averaging applied in the vertex-based approach. In the present study, we directly compared the average vertex-based heritability estimate within each region to the heritability estimate calculated in the traditional ROI approach. Inter-regional differences in heritability were greater in the ROI-based approach, and the degree of discrepancy between the approaches was dependent on the ROI size. Specifically, the curvilinear form of this relationship was consistent with the function expected due to an effect of degree of spatial averaging (which varies in the ROI-based approach and is constant and large in the vertex-based approach) on measurement error.

These results suggest that relative levels of spatial averaging, both by initial image processing and by methods that average data points on the surface maps, should be considered when interpreting inter-regional and inter-study differences in magnitudes of heritability. For very small regions, heritability is likely to be underestimated by an ROI-based approach. When the size of regions varies substantially, regional differences in heritability are likely to be overestimated. As predicted, we did also find a somewhat stronger association of heritability with ROI size for surface area than for cortical thickness. As suggested, this difference is consistent with the fact that inaccuracy of boundary placement will have a far greater effect on heritability estimates for surface area than on heritability estimates for cortical thickness.

Implications

Our results have implications for investigators who seek brain phenotypes that are likely to be associated with genetic polymorphisms. First, it is clear that investigators should determine whether they are searching for genes that are important only for determining the size of particular regions or for genes whose influences may be both regional and global. Given the high genetic correlation between areal expansion measures and total surface area, this issue is of particular relevance for surface area investigations. While in certain regions this distinction may not matter much, a number of regions appear to have relatively fewer unique genetic contributions and might not be good candidate phenotypes for association studies seeking to identify regionally-specific genes that determine surface area.

Also, if there is an interest in examining genetic associations or influences on particular regions from a given parcellation system, it appears that calculating heritabilities or associations at each point and then summarizing into mean values within a region may be preferable, at least for small regions. There are, however, some limitations to averaging heritabilities. There would be subsequent difficulties calculating confidence intervals, but there still may be instances where the vertex-based approach is preferred for understanding regional genetic contributions. At the least, interpretations of how heritability of surface area and thickness varies from region to region should be made while keeping in mind the possible effects of spatial averaging and boundary inaccuracy. For measures of area, the difference between the sensitivity of ROI-based and vertex-based approaches to variation in region size may be particularly striking because the placement of the boundaries between regions in each individual will have a direct and strong impact on the measured surface area of the region in the ROI-based approach.

Limitations

There were some minor limitations to our study that should be noted. First, our sample consisted of only middle-aged men, so the results may not generalize to women or other age groups. Second, although our vertex-based approach gave good qualitative information about patterns of heritability, we were not able to determine which observed regional variations were statistically meaningful. As a guide for future studies and generation of hypotheses, however, maps have the advantage of not being restricted by a priori boundaries that may not be genetically meaningful. Future studies will use patterns of genetic covariation in measures of surface area and cortical thickness to determine parcellation systems that are most relevant for genetic investigations.

Summary

In conclusion, use of spatially-unconstrained maps of areal expansion and cortical thickness is a powerful method to examine spatial patterns of the contribution of genetic and environmental influences. Even if particular regions based on functional or anatomical characteristics are of interest, a vertex-based approach may be the preferred first step for genetic studies. Further investigation of genetic patterning for areal expansion and cortical thickness using measures of genetic correlation among vertices is warranted, and will inform the search for particular genes or particular environmental conditions that influence cortical structure and the effect of development and disease on these measures.