The Impact of Genome-Wide Supported Schizophrenia Risk Variants in the Neurogranin Gene on Brain Structure and Function

The neural mechanisms underlying genetic risk for schizophrenia, a highly heritable psychiatric condition, are still under investigation. New schizophrenia risk genes discovered through genome-wide association studies (GWAS), such as neurogranin (NRGN), can be used to identify these mechanisms. In this study we examined the association of two common NRGN risk single nucleotide polymorphisms (SNPs) with functional and structural brain-based intermediate phenotypes for schizophrenia. We obtained structural, functional MRI and genotype data of 92 schizophrenia patients and 114 healthy volunteers from the multisite Mind Clinical Imaging Consortium study. Two schizophrenia-associated NRGN SNPs (rs12807809 and rs12541) were tested for association with working memory-elicited dorsolateral prefrontal cortex (DLPFC) activity and surface-wide cortical thickness. NRGN rs12541 risk allele homozygotes (TT) displayed increased working memory-related activity in several brain regions, including the left DLPFC, left insula, left somatosensory cortex and the cingulate cortex, when compared to non-risk allele carriers. NRGN rs12807809 non-risk allele (C) carriers showed reduced cortical gray matter thickness compared to risk allele homozygotes (TT) in an area comprising the right pericalcarine gyrus, the right cuneus, and the right lingual gyrus. Our study highlights the effects of schizophrenia risk variants in the NRGN gene on functional and structural brain-based intermediate phenotypes for schizophrenia. These results support recent GWAS findings and further implicate NRGN in the pathophysiology of schizophrenia by suggesting that genetic NRGN risk variants contribute to subtle changes in neural functioning and anatomy that can be quantified with neuroimaging methods.


SI 1.1 Antipsychotic history
Antipsychotic history was collected as part of the psychiatric assessment using the PSYCH instrument [1], and cumulative and current antipsychotic exposure was calculated using the chlorpromazine (CPZ) conversion factors of Woods et al. [2]. To calculate cumulative dose years, the following formulas were applied:

SI 1.2 Sternberg item recognition paradigm
Sternberg Item Recognition Paradigm (SIRP). For each task block, participants were presented with a prompt, ''Learn'', and then following a 0.5 s delay, they were shown a memory set comprised of one, three, or five digits for 6 s. This was followed by a ''probe epoch'', which lasted 38 s, and consisted of a series of 14 probe digits presented for 1.1 s with a jittered intertrial interval of ≤1.6 s. Participants used a button box to indicate whether each probe digit was a member of the memory set (''target'') or not (''foil''). They were instructed to respond as quickly and accurately as possible and were given a bonus of 5 cents for each correct response.
This bonus was provided after completion of the scan. Within each block, half of the items were targets and the other half were foils. For each participant, target-and foilbutton responses were randomly assigned to the right or left thumbs. The stimuli were projected onto a screen positioned on the head coil. Each of 3 runs contained 2 blocks of each of the 3 load conditions, presented in pseudorandom order, with the blocks of each load condition alternating with fixation (baseline) resting periods. Each run lasted 6 min.
The stimuli and responses were presented and collected using E-prime software (EPrime v1.1, Psychology Software Tools, Inc., Pittsburg, PA). Four participants (three patients and one control) were excluded from further analysis, because they completed a block with less than a 75% accuracy rate and/or with more than 6 probes not answered within a block. Please see [3] for more information and a schematic depiction of the SIRP protocol.

SI 1.3 Image acquisition and processing
The T1-weighted structural brain scans at each of the four sites were acquired with a coronal gradient echo sequence: TR=2530 ms for 3T, TR=12ms for 1.5T; TE=3.79 for 3T, TE=4.76ms for 1.5T; TI=1100 for 3T; Bandwidth=181 for 3T, Cross-site calibration and reliability of these acquisition sequences for each scanner and the experimental set up for functional imaging as well as potential site and scanner differences were investigated prior to the study [4,5]. Additionally, test-retest reliability of functional and structural imaging data from ten MCIC subjects, who were all scanned at all four sites, was analyzed. Results showed that, even with different scanner manufacturers and field strengths, activation variability due to site differences is small compared to variability due to subject differences [4][5][6].
Structural data were analyzed using FreeSurfer. Within the validated and complex surface reconstruction algorithm for structural brain scans, white matter segmentations were produced, topological defects in the surface were automatically corrected and the white and gray matter boundary was tessellated. After intensity normalization, the gray/white matter and the gray/cerebrospinal fluid borders were detected at the location where the greatest shift in intensity defines the transition to the other tissue class [7,8]. Final surfaces were used to calculate cortical thickness at each vertex on the tessellated surface as the closest distance from the gray/white boundary to the gray/cerebrospinal fluid boundary [9]. Segmentation and surface reconstruction quality were assured by manual inspection of all raw MRI volumes, segmented volumes in three planes and pial as well as inflated volumes.
Functional images were registered to a high-resolution T1 image of the same subject (using a new algorithm called Boundary-Based Registration [10]) and to the standard space defined by the MNI-152 atlas. We did this by first registering the T1 images to the standard brain using FLIRT [11,12] and then composing the functional-to-T1 and T1-to-standard registrations. We then fit a general linear model to the fMRI time course at each voxel in a whole brain model to estimate the average activation during the three loads of the probe condition in all trials. The magnitude of each Contrast Of Parameter Estimate (COPE), along with an estimate of its variability derived from model residuals, was passed to a second-level fixed effects analysis to combine COPE's from separate runs, yielding a composite T-statistic map for each contrast of interest for each subject.
Quality assurance steps for functional images included checks for whole-brain coverage of brain masks, motion and global mean intensity outlier timepoints, alignment of structural and functional scans, and registration problems (Epi to T1 and T1 to template). Outlier time frames in each fMRI data time series (detected using the artefact detection tools (ART) [13]) were defined by: (i) Global mean image intensity that differed by more than 3 standard deviations from the mean of the entire series of time frames in a scan, (ii) Displacement due to motion by more than 1 mm in the x, y or z direction relative to the previous time frame or (iii) Rotation due to motion by more than 0.1 rad around any of the three axes relative to the previous time frame.
We removed the outlier time frames through the use of nuisance regressors in the linear model. In the case of runs where more than 15% of the time frames were flagged as outliers, the entire run was dropped from the analysis or the subject had to be excluded.

SI 1.6 Monte-Carlo simulation
The Monte-Carlo simulation included the following steps: (1) An initial vertexwise threshold (VWT) was set to p=0.05 to form spatially contiguous areas of association (referred to as a cluster).
(2) The likelihood that a finding (cluster) of this size and magnitude (difference in thickness as specified by the VWT) would appear by chance, i.e. when using repeated random sampling, was tested using Monte-Carlo simulation with 10,000 repeats. This resulted in cluster-wise probabilities (CWP), which are reported using p-values throughout the results section.

SI 1.7 Estimates of cortical thickness and activation indices
We estimated the actual cortical thickness in millimeter (mm) in the identified cluster by labeling the cluster region as a region-of-interest (ROI), mapping it back on to each individual subject's unfolded surface by applying the same algorithm that morphed each subject's unfolded surface to the average spherical surface representation and extracting the average thickness for each individual.
Dorsolateral prefrontal cortex (DLPFC) ROI was derived from FreeSurfer cortical parcellations as described previously [14]. We obtained indices of activation for the DLPFC using the Contrast Of Parameter Estimates (COPEs) obtained from the second-level fixed-effects analysis for each participant. We applied an additional functional mask, based on the COPE of all loads (load 1, load 3 and load 5) versus fixation exceeding a threshold of Z = 2.3 and extracted the mean percent signal change (mean %∆), defined as the mean COPE of all working memory loads.

SI 1.8 Power Calculations
SI Figure S1 displays the effect size-by-power relationship for the functional analyses. Rs12541 accounted for 8.9% of total variance of mean cluster activation (R 2 change after accounting for the effect of diagnosis and tesla field strength), corresponding to f 2 =0.09 and a power of 99%.