Genetic analysis of cortical thickness and fractional anisotropy of water diffusion in the brain.

OBJECTIVES: The thickness of the brain's cortical gray matter (GM) and the fractional anisotropy (FA) of the cerebral white matter (WM) each follow an inverted U-shape trajectory with age. The two measures are positively correlated and may be modulated by common biological mechanisms. We employed four types of genetic analyses to localize individual genes acting pleiotropically upon these phenotypes.
METHODS: Whole-brain and regional GM thickness and FA values were measured from high-resolution anatomical and diffusion tensor MR images collected from 712, Mexican American participants (438 females, age = 47.9 ± 13.2 years) recruited from 73 (9.7 ± 9.3 individuals/family) large families. The significance of the correlation between two traits was estimated using a bivariate genetic correlation analysis. Localization of chromosomal regions that jointly influenced both traits was performed using whole-genome quantitative trait loci (QTL) analysis. Gene localization was performed using SNP genotyping on Illumina 1M chip and correlation with leukocyte-based gene-expression analyses. The gene-expressions were measured using the Illumina BeadChip. These data were available for 371 subjects.
RESULTS: Significant genetic correlation was observed among GM thickness and FA values. Significant logarithm of odds (LOD ≥ 3.0) QTLs were localized within chromosome 15q22-23. More detailed localization reported no significant association (p < 5·10(-5)) for 1565 SNPs located within the QTLs. Post hoc analysis indicated that 40% of the potentially significant (p ≤ 10(-3)) SNPs were localized to the related orphan receptor alpha (RORA) and NARG2 genes. A potentially significant association was observed for the rs2456930 polymorphism reported as a significant GWAS finding in Alzheimer's disease neuroimaging initiative subjects. The expression levels for RORA and ADAM10 genes were significantly (p < 0.05) correlated with both FA and GM thickness. NARG2 expressions were significantly correlated with GM thickness (p < 0.05) but failed to show a significant correlation (p = 0.09) with FA. DISCUSSION: This study identified a novel, significant QTL at 15q22-23. SNP correlation with gene-expression analyses indicated that RORA, NARG2, and ADAM10 jointly influence GM thickness and WM-FA values.

Introduction

The human cerebrum is a complex, multi-compartmental structure whose anatomy and function are influenced by individual genetic variations (Thompson et al., 2001; Glahn et al., 2007). Understanding the genetic mechanisms that control inter-subject cerebral variability is critical for deciphering the brain’s normal and pathological function. Recent discoveries of genetic factors that can increase the likelihood of developing neurodegenerative disorders such as dementias has emphasized the need to identify genes that influence integrity of cerebral tissue (Glahn et al., 2007; Meyer-Lindenberg, 2010; Stein et al., 2010a). Magnetic resonance imaging (MRI) facilitates the in vivo assessment of cerebral atrophy and structural integrity. Statistical genetic methods have been developed to measure the genetic modulation of inter-subject variability. By combining neuroimaging with genetics clinical investigators have promoted a better understanding of cortical variability (Walsh, 2000; Farnham et al., 2004; Gaitanis and Walsh, 2004; Jones et al., 2004; Edenberg et al., 2005; Klein et al., 2005; Kamarajan et al., 2006; Brouwer et al., 2010; Kochunov et al., 2010a; Joshi et al., 2011) and neurodegenerative processes (Biffi et al., 2010; Joshi et al., 2011).

We hypothesized that some of the genetic factors responsible for the integrity of cortical gray matter (GM) are also responsible for the integrity of cerebral white matter (WM). To test this hypothesis we identified genetic sources of shared covariance in indices of GM and WM integrity measured in a well-characterized population of Mexican Americans (Mitchell et al., 1996; Olvera et al., 2010). Integrity of the GM and WM compartments was assessed using the thickness of cortical GM and the fractional anisotropy (FA) of cerebral WM. GM thickness is calculated as the distance from the outer cortical surface to the inner cortical WM–GM boundary (Fischl and Dale, 2000), or a related symmetric measure (Lerch and Evans, 2005; Aganj et al., 2009). GM thickness is an indirect measure of a complex cortical architecture that is sensitive to cortical myelination, synaptic pruning, and cell density (Huttenlocher and Dabholkar, 1997). Integrity of cerebral WM was assessed using FA of water diffusion, measured from diffusion tensor imaging data. FA describes the directional selectivity of the random diffusion of water molecules (Beaulieu, 2002). WM–FA values are sensitive to various tissue properties including myelin content, intra-voxel axonal crossing, and axonal fiber density and diameter (Beaulieu, 2002). Though gray and WM measures are highly heritable (h2 = 0.4–0.8; p < 10−6) the extent to which these different measures are influenced by overlapping genetic factors is unknown (Brouwer et al., 2010; Kochunov et al., 2010a; Winkler et al., 2010; Chiang et al., 2011).

Gray matter thickness and FA values follow an inverted U-shape trajectory with age (Sowell et al., 2003; Gogtay et al., 2004; Salthouse, 2009; Kochunov et al., 2010b). The two measurements have a positive linear relationship across the lifespan (Kochunov et al., 2011). Age-related changes in cerebral myelination are thought to be the putative biological mechanisms responsible for the shape of their age-trajectories and the relationship between them (Sowell et al., 2003; Gogtay et al., 2004). Histological studies that compared GM thickness measured from MRI with those from the corresponding myelin-stained sections observed that the MR signal was sensitive to intra- and inter-layer myelination and the intra- and inter-layer density of glial cells. Likewise, inter-subject differences in the regional WM FA values have been reported to be predominantly due to the regional inter-subject differences in myelination levels and density of the glial cells (Song et al., 2005; Madler et al., 2008). Age-related changes in regional axonal myelin, reflected by changing FA values in WM, are likely to be similar to age-related changes in cortical myelin (Abe et al., 2002; Gao et al., 2009).

In our current study we attempted to localize chromosomal sites and identify individual genes that jointly influence these quantitative traits using bivariate genetic analyses. For correlated traits bivariate genetic analyses can reveal genetic factors that contribute jointly to both traits and can greatly increase the power of genetic discovery compared to univariate analyses (Amos and Laing, 1993; Schmitz et al., 1998; Amos et al., 2001). We performed analyses in this well-characterized population of Mexican Americans where we have demonstrated significant heritability (h2 = 0.4–0.6; p < 10−6) of cortical GM thickness (Winkler et al., 2010) and FA values (Kochunov et al., 2010a).

Materials and Methods

Subjects

The Genetics of Brain Structure and Function study recruited, among others, active participants from the San Antonio Family Heart Study (SAFHS; Mitchell et al., 1996) for brain imaging and neurocognitive assessment (Olvera et al., 2010). Analyzes were performed on data from 712 subjects in total (274 men/438 women). The Mexican American individuals are from large extended pedigrees selected randomly from the community. Subjects ranged in age from 19 to 85 years of age (47.9 ± 13.2 years) and were members of 73 families (9.7 ± 9.3 individuals/family; range 2–49). Subjects were excluded if they had MRI contraindications, history of neurological illnesses, stroke, or any other major neurological event. All experiments were performed with IRB approval from the University of Texas Health Science Center at San Antonio (UTHSCSA). All participants provided written informed consent on forms approved by the Institutional Review Boards at the University of Texas Health Science Center San Antonio and Yale University.

MR imaging

Imaging data were collected using a Siemens 3T Trio scanner located at the Research Imaging Institute, UTHSCSA.

T1-weighted imaging

This study used an MRI protocol specifically optimized for GM thickness measurements (Kochunov and Davis, 2009). The average GM thickness of the human cerebral cortex is approximately 2.5 mm, but varies regionally from 1.8–3.2 mm. High spatial resolution is necessary for accurate mapping of inter-subject GM thickness differences (Fischl and Dale, 2000). Our protocol was designed to collect data to resolve the cortical ribbon across the cortex using an isotropic spatial resolution of 0.8 mm (voxel volume = 0.5 mm3). T1-weighted contrast was achieved using a magnetization prepared sequence with an adiabatic inversion contrast-forming pulse (scan parameters: TE/TR/TI = 3.04/2100/785 ms, flip angle = 11°). A retrospective motion correction technique (Kochunov et al., 2006) was used to reduce artifacts related to subject motion.

Diffusion tensor imaging

Diffusion tensor images were collected using a single-shot, echo-planar, single refocusing spin-echo, T2-weighted sequence with a spatial resolution of 1.7 mm × 1.7 mm × 3.0 mm. The sequence parameters were: TE/TR = 87/8000 ms, FOV = 200 mm, axial slice orientation with 50 slices and no gaps, 55 isotropically distributed diffusion-weighted directions, two diffusion weighting values (b = 0 and 700 s/mm2) and three b = 0 images. These parameters were calculated using an optimization technique that maximizes the contrast to noise ratio for FA measurements (Jones et al., 1999).

Data processing

Processing of T1-weighted images

The processing of the T1-weighted images is summarized in Figure 1; see (Kochunov et al., 2007) for details. Briefly, the processing consisted of removing non-brain tissues, global spatial normalization, and radio frequency (RF) inhomogeneity correction. Non-brain tissues, such as skin, muscle, and fat, were removed using an automated skull stripping procedure and images were corrected for RF inhomogeneity (Smith et al., 2004). Next, images were imported into a freely distributed, structural analysis package, BrainVisa and processed using its cortical extraction and parcellation pipelines, as described by Kochunov et al. (2005). This pipeline extracts the pial and GM/WM interface surfaces, performs extraction, labeling, and verification of sulcal surfaces as described by Mangin et al. (2004). It segments the cortical landscape into 14 cortical regions using the primary sulcal structures (Cachia et al., 2003). For the frontal lobe, these regions are included: superior (S. Frontal), intermediate (M. Frontal), inferior frontal (I. Frontal), and pre-central gyri (Pre. Central). Parietal lobe was segmented into the post-central gyrus (Post. Central) and superior and inferior parietal lobules (S. and I. Parietal; Table 2). The temporal lobe was segmented into superior (S. Temp), intermediate and inferior temporal gyri (M. I. Temp), and the fusiform gyrus (Fusiform). Occipital lobe was segmented into lateral (L. Occipt) and medial occipital areas (M. Occipt). The limbic lobe was segmented into the cingulate and hippocampal gyri.

Figure 1

T1-weighted image processing pipelines. A T1-weighted image is skull-stripped, globally spatially normalized, and RF-inhomogeneity corrected (A). Next, cerebral hemispheres and cerebellum and identified and tissue classified (B). Cortical surfaces for GM and WM are calculated (C,D) and a homotopic erosion operation and crevasse detector are used to reconstruct sulcal surface as the medial surface of the two opposing gyral banks (E). The sulcal identification pipeline uses a congregation of 500 artificial neural network-based pattern classifiers to identify (F) sulcal landmarks, and to perform gyral segmentation of the cortex (G). GM thickness for 14 cortical regions (I), was calculated as the distance between the pial (C) and GM/WM interface surfaces (H).

Table 2

Phenotypic, corrected for age and sex, correlation coefficients between regional GM thickness values, and by tract FA measurements.

rp	GCC	BCC	SCC	Cingulum	CR	EC	IC	FO	SLF
S. frontal	0.21: p = 10 − 7	0.21: p = 10 − 7	0.11: p = 10−3	0.19: p = 10 − 6	0.19: p = 10 − 6	0.18: p = 10 − 6	0.07: p = 0.1	0.2: p = 10 − 7	0.17: p = 10 − 5
M. frontal	0.19: p = 10 − 6	0.21: p = 10 − 7	0.11: p = 10−3	0.19: p = 10 − 6	0.2: p = 10 − 7	0.2: p = 10 − 6	0.06: p = 0.1	0.21: p = 10 − 7	0.18: p = 10 − 5
I. frontal	0.16: p = 10 − 5	0.16: p = 10 − 4	0.12: p = 10 − 3	0.17: p = 10 − 6	0.18: p = 10 − 5	0.22: p = 10 − 7	0.08: p = 0.04	0.19: p = 10 − 6	0.17: p = 10 − 5
Post central	0.21: p = 10 − 7	0.15: p = 10 − 4	0.07: p = 0.1	0.27: p = 10 − 11	0.28: p = 10 − 12	0.31: p = 10 − 14	0.25: p = 10 − 10	0.23: p = 10 − 7	0.30: p = 10 − 14
Pre central	0.18: p = 10 − 6	0.16: p = 10 − 4	0.02: p = 0.7	0.21: p = 10 − 7	0.23: p = 10 − 7	0.28: p = 10 − 12	0.17: p = 10 − 5	0.17: p = 10 − 5	0.22: p = 10 − 7
S. parietal	0.31: p = 10 − 14	0.22: p = 10 − 7	0.13: p = 10−3	0.35: p = 10 − 18	0.35: p = 10 − 18	0.38: p = 10 − 21	0.32: p = 10 − 15	0.26: p = 10 − 10	0.36: p = 10 − 19
I. parietal	0.28: p = 10 − 12	0.22: p = 10 − 7	0.09: p = 0.05	0.34: p = 10 − 17	0.32: p = 10 − 15	0.35: p = 10 − 18	0.29: p = 10 − 13	0.25: p = 10 − 10	0.34: p = 10 − 17
S. temp	0.22: p = 10 − 7	0.17: p = 10 − 7	0.01: p = 0.7	0.25: p = 10 − 9	0.25: p = 10 − 9	0.29: p = 10 − 12	0.2: p = 10 − 7	0.23: p = 10 − 7	0.27: p = 10 − 11
M. I. temp	0.31: p = 10 − 14	0.23: p = 10 − 7	0.07: p = 0.1	0.33: p = 10 − 16	0.32: p = 10 − 15	0.36: p = 10 − 19	0.28: p = 10 − 12	0.27: p = 10 − 11	0.34: p = 10 − 17
Cingulate	0.27: p = 10 − 11	0.27: p = 10 − 11	0.17: p = 10 − 5	0.22: p = 10 − 7	0.18: p = 10−5	0.19: p = 10 − 6	0.04: p = 0.3	0.2: p = 10 − 7	0.16: p = 10 − 5
M. occipt	0.29: p = 10 − 13	0.25: p = 10 − 9	0.1: p = 0.01	0.34: p = 10 − 17	0.32: p = 10 − 16	0.36: p = 10 − 19	0.27: p = 10 − 11	0.25: p = 10 − 9	0.33: p = 10 − 16
L. occipt	0.29: p = 10 − 12	0.23: p = 10 − 7	0.1: p = 0.01	0.31: p = 10 − 14	0.31: p = 10 − 14	0.33: p = 10 − 16	0.26: p = 10 − 10	0.25: p = 10 − 9	0.31: p = 10 − 14
Fusiform	0.32: p = 10 − 15	0.23: p = 10 − 7	0.06: p = 0.1	0.3: p = 10 − 13	0.29: p = 10 − 13	0.34: p = 10 − 17	0.25: p = 10 − 10	0.26: p = 10 − 10	0.31: p = 10 − 14
Hippocampal	0.15: p = 10 − 4	0.15: p = 10 − 4	0.16: p = 10 − 5	0.17: p = 10 − 5	0.05: p = 0.06	0.11: p = 6·10−3	0.01: p = 1	0.11: p = 0.01	0.04: p = 0.3

Bolded values are significant after corrected for multiple (.

Skeletonized, average FA values are shown on the population average FA image. WM-tract labels for nine major tracts are taken from the Johns Hopkins DTI WM atlas. The average FA values were calculated for the following tracts: CR, corona radiata; SLF, superior longitudinal fasciculus, cingulate; EC, external capsule; IC, internal capsule, FO, fronto-occipital and the genu, body, and splenium of corpus callosum. The population-wise average FA values are color coded on the skeleton of the tract. The location of the axial slices (the Z-coordinate) is provided in the Talairach reference frame. The three segment of corpus callosum are shown on the mid-sagittal slice (Talairach X = 0).

Measurements of cortical gray matter thickness

Gray matter thickness is commonly defined as the distance from the outer cortical surface to the inner cortical WM–GM boundary (Figures 1D,H). Computationally, cortical thickness is determined by measuring the distance between two polygonal meshes, one representing the pial cortical surface and the other the white–gray interface. Multiple GM thickness measurement algorithms have been proposed to deal with the complex topography of the cerebral cortex. We used a GM thickness measurement tool distributed as a plug-in for BrainVisa. This tool uses a “normal-average” algorithm that has been reported to be a good compromise between accuracy and performance (Lerch and Evans, 2005). The GM thickness is measured as the Euclidean distance from an inner mesh vertex to the outer mesh along the direction normal to the inner mesh polygon. The measurement is repeated, inward, along the direction normal the outer mesh and the two distances are averaged (Kochunov et al., 2008). The consistency of the distance measurements is ensured by verifying that the line connecting the two surfaces does not intersect other polygons, inner or outer, along the way. GM thickness measurements were averaged for individual cortical areas for both hemispheres; the whole-brain GM thickness measurement was obtained by averaging all gyral GM thickness measurements.

Processing of diffusion-weighted images

Details for the processing of DTI scans are discussed elsewhere (Kochunov et al., 2007). In short, the tract-based spatial statistics (TBSS) software (Smith et al., 2006) was used for multi-subject analysis of FA images. FA images were created by fitting the diffusion tensor to the raw diffusion data (Smith, 2002). All FA images were non-linearly aligned to a group-wise, minimal-deformation target (MDT) brain (Kochunov et al., 2001). Next, individual FA images were averaged to produce a group-average anisotropy image. This image was used to create a group-wise skeleton of WM tracts, which encodes the medial trajectory of the WM fiber-tracts. Finally, FA values from each image were projected onto the group-wise skeleton of WM structures. This step helps account for residual misalignment among individual WM tracts. FA values are assigned to each point along a skeleton using the peak value found within a 20-mm distance perpendicular to the skeleton. The FA values vary rapidly perpendicular to the tract direction but very slowly along the tract direction. By assigning the peak value to the skeleton, this procedure effectively maps the center of individual WM tracts onto the skeleton.

Tract-based analysis

Whole-brain average FA values were calculated by averaging values for the entire WM skeleton. Average FA values were calculated for nine major WM tracts (Table 1) as described elsewhere (Kochunov et al., 2010b). Briefly, the population-based, 3D, DTI cerebral WM tract atlas developed at the John Hopkins University (JHU) and distributed as part of the FSL software package (Wakana et al., 2004) was used to calculate population average diffusion parameter values along the nine major WM tracts (Table 1; Figure 3). The JHU atlas was non-linearly aligned to the MDT brain and labels for individual tracts were transferred to the MDT brain using nearest-neighbor interpolation. Per-tract average values were calculated by averaging the values along the tracts for both hemispheres. The overall average FA values were calculated by averaging values for the entire WM skeleton.

Table 1

Cerebral white matter tracts.

Tract	Fiber type	Connections
Genu of corpus callosum (GCC)	C	Cerebral hemispheres
Body of corpus callosum (BCC)	C	Cerebral hemispheres
Splenium of corpus callosum (SCC)	C	Cerebral hemispheres
Cingulum	A	Cingulate gyrus/hippocampus
Corona radiata (CR)	P	Cortical/subcortical
External capsule (EC)	A	Frontal/temporal/occipital
Internal capsule (including thalamic radiation; IC)	P	Subcortical/brainstem/cortex
Superior/inferior fronto-occipital fasciculi (FO)	A	Frontal/parietal/occipital
Superior longitudinal fasciculus (SLF)	A	Frontal/temporal/occipital

Figure 3

Bivariate linkage analysis results for the GM thickness of the superior parietal lobule and FA values of the body of corpus callosum.

Genotyping

Linkage markers

Details of the genotyping procedure are described elsewhere (Kammerer et al., 2003) After DNA was extracted from lymphocytes, fluorescently labeled primers from the MapPairs Human Screening set (versions 6 and 8; Research Genetics, Huntsville, AL, USA) and PCR were used to amplify 417 microsatellite markers spaced at approximately 10-cM intervals across 22 autosomes. An automated DNA sequencer (ABI Model 377 with Genescan and Genotyper software; Applied Biosystems, Foster City, CA, USA) was used. The average heterozygosity index for these markers was approximately 0.76. The sex-averaged marker map was confirmed by deCODE genetics and markers not on this map were placed by interpolation based on physical location (Goring et al., 2007).

Association markers

The DNA for all subjects was processed using one million single nucleotide polymorphism (SNP) genotypes using the Illumina Human1M-Duo BeadChip microarray. All SNP genotypes were checked for Mendelian consistency using the program SimWalk2 (Sobel and Lange, 1996). Maximum likelihood techniques that account for pedigree structure were used to estimate allelic frequencies, taking pedigree structure into account (Boehnke, 1991). To avoid excluding individuals with missing SNP data from the analysis, we employed likelihood-based imputation with the MERLIN software package (Abecasis et al., 2002; Burdick et al., 2006).

Gene-expression measurements

Gene-expression measurements, collected 17 years ago (from 1991 to 1995), were available for 371 (150/221 males/females; average age = 47.8 ± 12.9 years) of our subjects. The details of gene-expression measurements are discussed elsewhere (Goring et al., 2007). Total RNA was isolated from lymphocyte samples using a modified procedure of the QIAGEN RNeasy 96 protocol for isolating total RNA from animal cells using spin technology (QIAGEN, Inc.). Its integrity was verified and anti-sense RNA was synthesized, amplified and purified. Hybridization of anti-sense RNA was carried out using Illumina’s BeadChip 6 × 2 protocol. Samples were scanned on the Illumina BeadArray 500GX Reader. Expression levels were z-normalized to make the expression phenotypes comparable among individuals and across transcripts (Goring et al., 2007). This minimized the influence of overall signal levels, which may reflect RNA quantity and quality rather than true individual differences. The z-normalization was performed on the abundance values of all transcripts within individuals using decile percentage bins of transcripts, grouped by average log-transformed raw signals across individuals. This was followed by a linear regression against the individual-specific average log-transformed raw signal and its squared value. Finally, each transcript’s residual expression scores were normalized by using an inverse Gaussian transformation across individuals to ensure the data normality (Goring et al., 2007).

Genetic analyses

Four sets of analyzes were performed to ascertain the magnitude of shared genetic variance between the GM thickness and FA and localize genetic factors that jointly influence whole-brain average and regional measurements. First, we formally estimated the significance of the phenotypic, genetic, and environmental correlation between both traits with a bivariate correlation analysis (Almasy et al., 1997; Almasy and Blangero, 1998). Further, we used bivariate whole-genome linkage analyses to formally localize chromosomal regions that jointly influenced the traits that showed a significant phenotypic correlation. Finally, gene-level localization was performed using SNP association and gene-expression analyses to identify putative genes located within chromosomal regions located under quantitative trait loci (QTL) peaks. All analyzes were performed using the SOLAR software package with age, sex, age × sex, age2, age2 × sex included as covariates.

Quantitative trait linkage analysis

Quantitative genetic analyses were performed using a variance-components, maximum likelihood method implemented in SOLAR (Almasy and Blangero, 1998). Bivariate quantitative trait linkage analyses of GM–FA trait pairs were performed to localize potentially pleiotropic QTLs (Almasy et al., 1997). Bivariate QTL analysis was chosen over the simpler, univariate linkage analysis on an arithmetic combination of two traits approach, because bivariate genetic analyses can greatly increase the power of genetic discovery compared to the univariate approach (Amos et al., 2001). Linkage significance was assessed by comparing the likelihood of a model with QTL-specific partial heritabilities estimated simultaneously for both traits to a null model in which both partial heritabilities were constrained to zero. Significance p-values from the 2 degree of freedom (df) bivariate linkage tests were converted to 1 df-equivalent LOD scores to aid comparison with univariate linkage results. We chose LOD scores of 2.0 and 3.0 as the thresholds for “suggestive” (likely to occur one or fewer times by chance in a genome scan) or significant scores (genome-wide p-value = 0.05), respectively.

Single nucleotide polymorphism association analyses

Single nucleotide polymorphism association analyses were performed for chromosomal regions located under the QTL peaks. SNP genotype values were represented as 0, 1, or 2 copies of the minor allele (or, for missing genotypes, the weighted covariate based on imputation). The significance of association was calculated by comparing the model that included the SNP value with the null models that incorporated standard polygenic effects and fixed covariates. The SNP genotype values were included as covariates in the variance-components mixed models for measured genotype analyses (Boerwinkle and Sing, 1986). The critical p-value, corrected for multiple comparisons, was computed based on the effective number as suggested in (Moskvina and Schmidt, 2008).

Results

Phenotypic correlation analyses between whole-brain average GM thickness and FA values reported that the two traits were significantly and positively correlated (rP = 0.27; p < 10−7). The positive sign of the phenotypic correlation coefficient suggested that the same genetic factors associated with higher FA values were linked to progressively higher GM thickness values. There are 126 possible GM–FA pairs between the 14 regional GM thickness values and nine regional FA values (Table 2). 114 out of 126 GM–FA correlation analyses showed a significant (p < 0.05) phenotypic correlation (rp), and 101 of the 126 were significant when corrected for multiple comparisons (p < 0.0004, Bonferroni correction for 126 tests).

Bivariate whole-genome QTL analysis of whole-brain average GM thickness and FA value produced a single peak of suggestive (LOD = 2.56) significance on the q-arm of chromosome 15 (ch15:54) located at 15q22–23. Regional, bivariate QTL analysis was limited to the chromosome 15 and performed for 101 trait pairs that showed significant phenotypic correlation (Table 2). This analysis produced 13 significant (LOD > 3.0) and 51 suggestive (LOD > 2.0) linkage peaks on the chromosome 15q22–23 (Table 3). The highest LOD (LOD = 4.51) score was observed for the GM thickness values of the superior parietal lobule and FA values of the body of CC (Figures 3 and 4, Table 3, top). Assuming that a LOD ≥ 3.0 criterion controls the risk of one or more false positives over all chromosomes at level 0.05, then 101 analyses should generate no more than 5.05 false positives, on average; thus 13 positive results is more than twice the number of positive results expected by chance. The locations of significant LOD peaks were flanked by markers rs12901270 and rs1574119 (15q22–23; Table 4; Figure 4).

Table 3

Logarithm of Odds (LOD) values and location of the peak [chromosome: location (cM)] for traits.

LOD	GCC	BCC	SCC	Cingulum	CR	EC	IC	FO	SLF
S. frontal	0.94 (15:54)	2.07 (15:55)	N/A	1.54 (15:51)	1.53 (15:53)	1.54 (15:54)	N/A	1.55 (15:54)	1.08 (15:51)
M. frontal	1.42 (15:54)	2.34 (15:55)	N/A	2.04 (15:54)	1.76 (15:52)	1.8 (15:53)	N/A	1.54 (15:54)	1.38 (15:52)
I. frontal	1.82 (15:54)	LOD < 2.00	N/A	2.18 (15:53)	1.74 (15:52)	1.69 (15:53)	N/A	1.54 (15:54)	1.57 (15:52)
Post central	1.67 (15:52)	2.52 (15:54)	N/A	1.34 (15:54)	0.88 (15:44)	1.12 (15:53)	1.07 (15:107)	1.40 (15:48)	1.11 (15:53)
Pre central	1.77 (15:54)	2.95 (15:56)	N/A	1.65 (15:54)	1.08 (15:44)	0.96 (15:54)	0.65 (15:107)	1.17 (15:48)	0.96 (15:53)
S. parietal	3.34 (15:54)	4.50 (15:54)	N/A	3.03 (15:54)	2.94 (15:52)	2.97 (15:54)	2.95 (15:54)	2.87 (15:54)	3.04 (15:54)
I. parietal	2.72 (15:54)	3.63 (15:54)	N/A	2.49 (15:54)	2.13 (15:52)	2.21 (15:62)	2.42 (15:54)	1.19 (15:48)	2.21 (15:53)
S. temp	2.06 (15:54)	2.97 (15:55)	N/A	1.79 (15:54)	1.73 (15:52)	1.74 (15:53)	1.53 (15:53)	2.21 (15:49)	1.95 (15:54)
M. I. temp	2.40 (15:54)	4.03 (15:54)	N/A	2.70 (15:54)	2.51 (15:54)	2.62 (15:54)	2.41 (15:54)	2.34 (15:54)	2.60 (15:54)
Cingulate	2.81 (15:54)	3.43 (15:57)	2.23 (15:65)	2.68 (15:65)	2.12 (15:63)	2.23 (15:54)	N/A	1.8 (15:63)	1.89 (15:54)
M. occipt	3.41 (15:54)	4.04 (15:55)	N/A	2.85 (15:54)	2.47 (15:54)	2.64 (15:54)	2.60 (15:54)	2.27 (15:54)	2.48 (15:54)
L. occipt	3.19 (15:54)	4.17 (15:54)	N/A	2.65 (15:54)	2.62 (15:54)	2.63 (15:54)	2.45 (15:54)	2.29 (15:54)	2.69 (15:53)
Fusiform	2.17 (15:52)	3.92 (15:55)	N/A	2.62 (15:55)	1.60 (15:54)	2.24 (15:54)	2.26 (15:54)	2.38 (15:54)	2.30 (15:54)
Hippocampal	2.44 (15:57)	3.58 (15:57)	1.59 (15:55)	2.92 (15:57)	N/A	N/A	N/A	N/A	N/A

Bolded values indicated a significant (LOD > 3.0) QTL. N/A are reported for trait pairs that did not have a significant phenotypic correlation coefficient (.

Figure 4

The LOD peak was flanked by markers rs12901270 and rs1574119 and covered 5.6 Mb of genome (top). The chromosomal region under the linkage peak harbored 20 genes with 22 gene expression measurements (middle). Results (-log10 of the p-value) of the association analysis for 1565 SNPs located under the linkage peak are plotted along the x-axis (bottom). The red line corresponds to the significant evidence value of p = 5·10−5. The green line corresponds to the suggestive evidence value of p = 1·10−3.

Table 4

SNP name (p) and (genetic association) for the highest and all suggestively (.

SNP	GCC	BCC	Cingulum	SLF
S. parietal	RS460458(0.0001) (intergenic)	RS460458 (0.0002) (intergenic)	RS7183502(0.0001) (intergenic)	RS7183502(0.0003) (intergenic)
	RS2456930(0.0003) (intergenic)	RS1545540 (0.0003) (intergenic)		RS17302400(0.0004) (intergenic)
		RS7174111 (0.0004) (intergenic)		rs12438866 (0.0008) (RORA)
		RS1869486 (0.0005) (RORA)
		RS10459628 (0.0005) (NARG2)
		RS12910118 (0.0006) (intergenic)
		RS1869485(0.0007) (intergenic)
		RS12594972 (0.001) (RORA)
I. parietal	N/A	RS1869485(0.001) (NARG2)	N/A	N/A
		RS7174111(0.0003) (intergenic)
		RS17270188 (0.0003) (intergenic)
		RS1545540 (0.0003) (intergenic)
		RS10459628 (0.0003) (NARG2)
		RS11071529 (0.0005) (NARG2)
		RS1869486 (0.001) (RORA)
M. I. temp	N/A	RS1545540 (6·10−5) (intergenic)	N/A	N/A
		RS1869486 (7·10−5)(RORA)
		RS7174111(9·10−5) (intergenic)
		RS1869485 (0.0001) (intergenic)
		RS10459628 (0.0002) (NARG2)
		RS11071529 (0.0002) (NARG2)
		RS12910118 (0.0002) (intergenic)
		RS460458(0.0004) (intergenic)
		RS4561404 (0.0006) (intergenic)
		RS12909890 (0.0007) (RORA)
		RS2456930 (0.0009) (Intergenic)
		RS7169359 (0.0009) (intergenic)
		RS12594972 (0.0009) (RORA)
Cingulate	N/A	RS4620904 (0.001) (intergenic)	N/A	N/A
M. occipt	RS17270188(0.001) (RORA)	Rs10459628(0.0003 (NARG2)	N/A	N/A
		RS1545540 (0.0003) (intergenic)
		RS7174111(0.0004) (intergenic)
		RS1869485 (0.0008) (NARG2)
		RS12910118 (0.0008) (intergenic)
L. occipt	RS7183502(0.002) (intergenic)	RS7174111(0.0003 (NARG2)	N/A	N/A
		RS1545540 (0.0005) (intergenic)
		RS12910118 (0.0006) (intergenic)
		RS7174111(0.0007) (intergenic)
		RS10459628 (0.0008) (NARG2)
Fusiform	N/A	RS1545540(8·10−5 (NARG2)	N/A	N/A
		RS1869486 (9·10−5)(RORA)
		RS7174111(0.0001) (intergenic)
		RS10459628 (0.0002) (NARG2)
		RS1869485 (0.0003) (intergenic)
		RS12910118 (0.0003) (intergenic)
		RS11071529 (0.0005) (NARG2)
		RS460458 (0.0005) (intergenic)
		RS4561404 (0.0006) (Intergenic)
		RS16940027 (0.0009) (Intergenic)
Hippocampal	N/A	RS16940027 (0.0001) (Intergenic)	N/A	N/A
		RS8041394 (0.0003) (GTF2A2)
		RS8038077 (0.0005) (RORA)
		RS1693534 (0.0007)(RNF111)
		RS7882 (0.0008) (GTF2A2)

N/A is reported for trait pairs that did not have a significant (LOD > 3.0) QTL.

Bivariate SNP association analyses were performed for the thirteen trait pairs that showed significant QTL to achieve finer localization of the linkage region. A total of 1565 SNPs were extracted from the region of 5.6 Mb (Figure 4 middle). The effective number of SNPs, corrected for linkage disequilibrium, was 985 with a Bonferroni corrected significance criterion of 5·10−5. No SNP met this criterion of significance. The highest significance (6·10−5) was observed for the intergenic rs154554. This SNP is located in the immediate vicinity (103 kb pairs) of the p-end of the NMDA receptor regulated 2 (NARG2) gene (Table 4). The highest significant SNPs for the remaining 12 trait pairs, were localized to NARG2 gene (N = 5), retinotopic acid receptor (RAR)–related orphan receptor alpha (RORA) gene (N = 1), and the rest (N = 6) were intergenic.

Next, we performed two post hoc analyses to identify potential clustering patterns of the SNP associations. We identified the 22 suggestively (p ≤ 10−3) significant SNPs across all thirteen trait pairs (Figure 3; Table 5). None was in linkage disequilibrium with others based on the r2 > 0.80 thresholds. Nine of the 22 polymorphisms (41%) were localized to two neighboring genes: NARG2 (N = 3) and RORA (N = 6; Table 5; Figure 3, bottom). Two (rs8041394, p = 0.0003 and rs7882, p = 0.0001) were localized to the General Transcription Factor IIA, 2 gene (Table 5; Figure 3, bottom). One (rs1693534, p = 0.0007) was localized to the Ring Finger Protein 111 (RNF11) gene (Figure 3, bottom). Among the remaining (N = 7) intergenic polymorphisms, we identified rs2456930 (p = 0.0003; Figure 3, bottom). This polymorphism was reported as a highly associated (p = 3.1·10−7) finding resulting from genome-wide search by Stein and colleagues for bilateral temporal lobe volume in which each copy of the G allele correlated with reductions of temporal lobe volume (Stein et al., 2010a). In the current study, this SNP was observed as the second most significant association in the analysis of the GM thickness of the superior parietal lobule and the FA of the genu of corpus callosum. Residual (corrected for age and sex) GM thickness and FA values for subjects with homozygous minor (AA) allele were 3.8 and 3.4% higher than for subjects with GG allele (Figure 5).

Table 5

Suggestively (.

Name	Location	−log10(p)	Gene
rs16940027	58593730	3.95	Intergenic
rs1693534	59281962	3.14	RNF111
rs17302400	59824987	3.39	Intergenic
rs8041394	59949178	3.56	GTF2A2
rs7882	59951573	3.01	GTF2A2
rs7183502	60334878	3.49	Intergenic
rs12910118	60558741	3.63	Intergenic
rs7169359	60633676	3.01	Intergenic
rs4561404	60697836	3.24	Intergenic
rs1545540	60699608	4.21	Intergenic
rs7174111	60706446	4.02	Intergenic
rs11071529	60709904	3.68	NARG2
rs10459628	60710269	3.78	NARG2
rs1869485	60746605	3.89	NARG2
rs12594972	60794649	3.00	RORA
rs1869486	60794807	4.14	RORA
rs17270188	60799909	3.52	RORA
rs8038077	61004647	3.27	RORA
rs12909890	61017977	3.11	RORA
rs12438866	61068520	3.02	RORA
rs2456930	62687339	3.56	Intergenic
rs460458	62692778	3.74	Intergenic

*Location is reported based on the 37.1p map.

Figure 5

Residual (corrected for age and sex) GM thickness and FA values for subjects with grouped by genotype for the rs2456930 polymorphism. The allele frequencies for this population were AA = 14.6%, AG = 48.3%, and GG = 36.9%.

Next, we identified the frequency at which polymorphisms appeared in the list of top-10 most significant polymorphisms for the thirteen trait pairs (Table 6). This list was composed of 46 unique polymorphisms (Table 6). None of the SNPs were in linkage disequilibrium with others based on the r2 > 0.80 threshold. Three polymorphisms (rs1545540, rs1869485, and rs717411) had the highest co-occurrence (occurred for 11 out of 13 trait pairs). Two were intergenic (rs1545540 and rs7174111) and one (rs1869485) resided with the NARG2 gene. Both intergenic (rs1545540 and rs7174111) polymorphisms were proximal (12 Kbt and 5 Kbt downstream, respectively) to the NARG2 gene. We further analyzed the genetic association of the 130 most significant polymorphisms. Polymorphisms localized to two genes: RORA and NARG2 accounted for about 22% of the polymorphisms (Table 7).

Table 6

Frequency repeats by SNPs within .

SNP	Number of repeats	Average p	Location*	Minor allele	Frequency	Gene
RS1545540	11	0.001	60699608	G	0.25	Intergenic
RS1869485	11	0.002	60746605	G	0.31	NARG2
RS7174111	11	0.001	60706446	T	0.24	Intergenic
RS10459628	10	0.001	60710269	T	0.25	Intergenic
RS12910118	8	0.001	60558741	G	0.18	intergenic
RS1869486	6	0.001	60794807	T	0.18	RORA
RS460458	6	0.001	62692778	G	0.10	intergenic
RS11071529	5	0.001	60709904	G	0.31	intergenic
RS7174111	5	0.002	60706446	T	0.24	Intergenic
RS8031866	5	0.002	57445040	A	0.28	Intergenic

*Location is reported based on the 37.1p map.

Table 7

Analysis of repeats by gene within top-10 significant association for each of the thirteen traits.

Gene	Count	Average p	SNP
RORA	14	0.001	RS1869486, RS8038077, RS10431796, RS10431796, RS12909890, RS341401
NARG2	14	0.002	RS1869485, RS1107153, RS17270146, RS11071537
GTF2A2	4	0.002	RS1871500, RS7882, RS8041394, RS8041394
ANAX2	1	0.005	RS7182242
RS2456930	1	0.002	RS2456930

Finally, we tested 22 gene-expression measurements from 20 genes residing under the linkage peak. These data were collected 17 years ago and were available for 371 subjects. Bivariate genetic correlation analysis was used to calculate phenotypic correlation between gene-expression values and whole-brain GM thickness and FA values separately. Phenotypic (rP) correlation values for two transcripts of the RORA gene and one transcript of disintegrin and metalloproteinase domain-containing protein 10 (ADAM10) were significant (p < 0.05) for both GM thickness and FA values (Table 8).

Table 8

Correlation coefficient between the expression level and whole-brain average GM thickness and FA values.

Transcript name	rP WB–GM	rP WB–FA	Gene	Location
GI_19743904-A	0.008	0.02	RORA	15q22.2
GI_4557250-S	0.009	0.02	ADAM10	15q22
GI_19743904-I	0.01	0.02	RORA	15q22.2
ρP Significant for GM only
GI_37202122-S	0.04	0.09	NARG2	15q22.2
ρP Significant for FA only
GI_4757855-S	0.14	0.04	BNIP2	15q22.2

A transcript for NARG2 gene was correlated with the GM thickness and its correlation coefficient with FA approached significance (p < 0.10). However, these correlations did not reach significance (p = 0.002) once corrected for multiple (N = 22) comparisons. Additionally, the ADAM10 expressions were significantly and positively correlated with the NARG2 expressions (r = 0.44; p < 10−5). No further significant (r = < 0.10) correlations were observed among other transcripts.

Discussion

The life-long relationship between cortical GM thickness and WM FA values (Kochunov et al., 2011) suggests that shared genetic factors might be influencing their age-trajectories. We examined these trends in a large, well-characterized sample of Mexican–American participants from the Genetics of Brain Structure and Function study and demonstrated a significant shared variance between them. We used bivariate genetic analyses to localize chromosomal regions and to identify individual genes involved in this relationship. The bivariate linkage analysis of the whole-brain average GM thickness and FA values produced a suggestive locus on the chromosome 15q22–23 near the previously reported univariate QTL for the FA values (Kochunov et al., 2010a). Bivariate linkage analyses of regional GM and FA trait pairs produced thirteen significant loci in the same region. This chromosomal region of 5.6 Mb harbors over 30 genes, many of which are critical to normal brain development. Reports on individuals with chromosomal deletions in this region have described severe congenital defects such as holoprosencephaly (Machado et al., 2010) and cerebral dysmorphia (Koivisto et al., 1999; Lalani et al., 2006). A microdeletion of 1 Mb of genome in the immediate vicinity to the rs1574119 resulted in cerebral dysmorphia, autism, mental retardation, and developmental delay (Smith et al., 2000). Additionally, a study in 175 families of Latino descent reported this region as a locus of suggestive (LOD = 2.11) significance for harboring genes that act as risk factors for schizophrenia (Escamilla et al., 2009).

We attempted to identify individual genes using SNP associations and correlations with the gene-expression analyses. The bivariate association analysis for the 1565 SNPs uniformly covering this region did not return a definitive finding. It produced several useful results. 41% of the potentially (p < 10−3) significant SNPs were localized to two proximal genes NARG and RORA. These two genes span 14% of the genome between markers rs12901270 and rs1574119 and therefore received a higher (two to three times) number of associations than would be expected by chance. Examining the top-10 most significant SNPs for each of the 13 trait pairs found that 22% (28 out of 130) of the polymorphism were also localized to the NARG2 or RORA genes. Correlation with gene-expression values identified three transcripts whose expression levels were significantly (r < 0.05) correlated with whole-brain average measurements for both traits. Two of these transcripts were from the RORA gene. The third transcript was from the ADAM metallopeptidase domain 10 (ADAM10) gene. The transcript for NARG2 gene showed a significant correlation with GM thickness and a suggestive (p = 0.09) correlation with the FA values.

We consider RORA as the likely candidate for pleiotropic association between FA and GM thickness. RORA was discovered as the site of the “staggerer” mutation in a mouse model that is characterized by cerebellar ataxia and profound neurodegeneration (Hamilton et al., 1996; Dussault et al., 1998). RORA encodes a protein that acts as a constitutive activator of transcription and a receptor for glucocorticoids (Jarvis et al., 2002). It is expressed in many tissues including cortical, subcortical and cerebellar neurons and neuroglial cells (Jarvis et al., 2002). In neuronal and glial cells the RORA transcription factor exerted both neuroprotective and anti-inflammatory functions. It provides neuroprotection by increasing the expression of antioxidant proteins, which makes neurons less sensitive to apoptotic stimuli (Boukhtouche et al., 2006). In glial cells its expression is up-regulated by pro-inflammatory cytokines and loss of its function leads to pro-inflammatory activation in astrocytes. Epidemiological genetic analyses implicated this gene in several disorders including: depression (Garriock et al., 2009; Lavebratt et al., 2009; Le-Niculescu et al., 2009; McGrath et al., 2009; Terracciano et al., 2010), Parkinson’s disease (Edwards et al., 2010), age-related macular degeneration (Schaumberg et al., 2010), and attention deficit hyperactivity disorder (Neale et al., 2008). Its lymphocytic expression levels in autistic individuals were lower than in healthy controls and correlated with the severity of the disorder (Hu et al., 2006, 2009; Nguyen et al., 2010; Sarachana et al., 2010).

We consider NARG2 as a potential candidate gene. Relatively little is know about this gene. It codes the N-methyl-d-aspartate (NMDA) receptor, a glutamate-gated ion channel (Sugiura et al., 2001). It is highly expressed in the neonatal brain and plays an important role in brain development by regulating neuronal survival, migration, and proliferation (Sugiura et al., 2004). Though NARG2 is located next to RORA on the q-terminal, the lack of correlation between NARG2 and RORA indicates that NARG2 modulation of GM thickness is unrelated to its proximity to RORA.

Notably, the gene-expressions for the α-secretase, ADAM10 gene were correlated with both GM thickness and FA values. ADAM10 is an α-secretase capable of anti-amyloidogenic proteolysis of the amyloid precursor protein. This gene is considered a novel treatment target for Alzheimer’s disease (AD; Kim et al., 2009). In vitro studies demonstrated that α-secretase is responsible for the non-amyloidogenic cleavage of Alzheimer’s plaque protein. Up-regulation of its expression can prevent intracellular and extracellular accumulation of the toxic β-amyloid species (Aβ; Jarvis et al., 2010). In vivo animal studies in Alzheimer mouse models showed that blocking ADAM10 expression, or its trans-synaptic transport, led to a rapid increase in the concentration of Aβ plaque and increased brain atrophy (Epis et al., 2010; Wang et al., 2010; Zhang et al., 2010). Reduced ADAM10 trans-synaptic trafficking rates were also demonstrated in the hippocampal tissue samples from AD patients (Marcello et al., 2010). Two ADAM10 polymorphisms were associated with increased risk of developing a late-onset AD (Kim et al., 2009). These polymorphisms significantly attenuated α-secretase activity of (>70% decrease), and elevated Aβ levels (1.5–3.5-fold).

Our study replicated an association between rs2456930 and brain integrity (Stein et al., 2010b). Stein and colleagues identified this SNP as a significant GWAS finding in 742 participants in the Alzheimer disease neuroimaging initiative (ADNI). The genotypic frequencies for this SNP in our subjects (AA = 14.6%, AG = 48.3%, GG = 36.9%) were nearly identical to those reported by Stein and colleagues (AA = 13.0%, AG = 49.5%, GG = 37.5%; Stein et al., 2010a). The magnitudes and the direction of the genetic effects for were also similar. In both populations the homozygous minor allele (AA) exerted neuroprotective properties. The AA carriers had higher GM thickness and FA values (3.8 and 3.4%, respectively) than carriers of the homozygous major allele (GG). Stein et al. (2010a) reported a 3.13% higher volume of the temporal lobe for AA versus GG carriers. This polymorphism is found within a “gene desert” over 200 kb away from known genes. It is located within a 10 kb of a region with both promoter and enhancer histone markings which characterizes it as an area with potential regulatory function (Schones and Zhao, 2008). It potentially has high functional relevance, even though the gene target is not yet known. Our finding suggests the neuroprotective effects of this polymorphism can be observed in normal aging and therefore, are probably not directly associated with AD pathology.

Our analysis of a population of healthy aging Mexican Americans has identified a novel, significant QTL at the chromosome 15q22–23. SNP association and correlation with gene-expression analyses identified three genes located under the linkage peak, RORA, NARG2, and ADAM10, showed potential for jointly influencing cortical thickness and cerebral WM integrity. Genetic factors that are jointly influencing two different imaging traits may be acting through neuroprotective effects. The scientific literature indicates that RORA and ADAM10 genes are considered potential treatment targets for neurodegenerative disorders, and two mutations of ADAM10 gene are being considered as candidates for AD susceptibility. Localization of these genes reached only trend significance and further deep sequencing and functional variants analyses will be required for the true identification.

Limitations

This study was performed in Mexican Americans, a population with significant Native American admixture. If relatively rare variants are involved in the determination of quantitative variability, we may expect considerable differences in the localization of the most important genetic loci across populations (Blangero et al., 2003). Linkage studies of such complex phenotypes cannot be used to exclude genetic regions for important QTLs. Therefore, the lack of concordance cannot be interpreted as evidence against the hypothesis that a QTL exists in a particular genomic region.

The gene-expression correlation analyses were limited by two factors. First, these were available for about half (N = 371, 150/221 males/females; average age = 47.8 ± 12.9 years) of our subjects. Second, the collection of gene-expression data preceded the acquisition of brain images by 17 years. This prospective analysis approach may limit future replications studies.

Another limitation of this analysis is the use of Bonferroni correction for multiple comparisons. Permutation analysis can be used to infer empirical p-values and provide less conservative correction for multiple comparisons than that obtained with Bonferroni. However, the variance-components estimation in large extended pedigrees is not trivial and sometimes may fail to converge as permutation does not take family structures into the account. Hence permutation would require significant additional effort, and since the Bonferroni thresholds are conservative, we chose to use that methodology.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.