The influence of polygenic risk for bipolar disorder on neural activation assessed using fMRI.

Genome-wide association studies (GWAS) have demonstrated a significant polygenic contribution to bipolar disorder (BD) where disease risk is determined by the summation of many alleles of small individual magnitude. Modelling polygenic risk scores may be a powerful way of identifying disrupted brain regions whose genetic architecture is related to that of BD. We determined the extent to which common genetic variation underlying risk to BD affected neural activation during an executive processing/language task in individuals at familial risk of BD and healthy controls. Polygenic risk scores were calculated for each individual based on GWAS data from the Psychiatric GWAS Consortium Bipolar Disorder Working Group (PGC-BD) of over 16 000 subjects. The familial group had a significantly higher polygene score than the control group (P=0.04). There were no significant group by polygene interaction effects in terms of association with brain activation. However, we did find that an increasing polygenic risk allele load for BD was associated with increased activation in limbic regions previously implicated in BD, including the anterior cingulate cortex and amygdala, across both groups. The findings suggest that this novel polygenic approach to examine brain-imaging data may be a useful means of identifying genetically mediated traits mechanistically linked to the aetiology of BD.

Introduction

Evidence from family, twin and adoption studies has consistently shown a significant genetic contribution to the aetiology of bipolar disorder (BD) with heritability estimates of between 0.6 and 0.8.[1, 2] The identification of susceptibility genes for BD, from both candidate gene and genome-wide association studies (GWAS), means there are now a number of single-nucleotide polymorphisms (SNPs) that are robustly associated with risk of BD.[2, 3] Individually these variants account for a small proportion of the heritability of BD, suggesting that each risk allele may have a small additive and cumulative effect across the genome. Evidence that liability to BD is due to the accumulation of many common genetic variants is consistent with the threshold-liability model of mental illness,[4] and with the observation that traits linked to BD are expressed to a greater degree in unaffected family members than unrelated individuals.[5]

In order to test the polygenic model of BD, two studies have examined whether a polygenic risk score derived from a discovery data set can be used to predict BD in an independent target sample.[6, 7] Both studies concluded that the polygenic contribution could explain a small but significant portion of the total variation in liability, and both demonstrated that a polygenic risk score could be used to classify cases significantly better than chance in the target data set (maximum effect size: R2=0.028, P=1.7 × 10−9). This proportion of explained variance is expected to increase with larger samples and with more dense SNP genotyping arrays that capture a higher proportion of total genetic variation.[8]

In order to better understand the genetic contributions to BD, previous imaging studies have investigated the structural and functional effects of individual SNPs in affected and unaffected individuals, for example, CACNAIC, DISC1, BDNF, DGKH and NRG1,[9, 10, 11, 12, 13] or the interaction of a small number of SNPs.[14, 15, 16, 17] Although these studies have been informative regarding the effects of specific SNPs, a polygenic model is likely to capture a greater proportion of the shared genetic architecture of both the BD and imaging phenotypes. It is therefore likely to identify brain regions linked to the genetic aetiology of BD that have broader applicability than the effects of individual alleles.

In the current study we generated a polygenic risk profile score for BD based on all available Psychiatric GWAS Consortium summary data for BD (https://pgc.unc.edu/). We then sought to determine the impact of this score on brain activation during performance of an executive function/language task. This task has been shown to differentiate patients with BD, and those at increased familial risk, from healthy controls.[18, 19] Moreover, it probes neuropsychological deficits in executive function, verbal initiation and verbal fluency observed in patients with BD.[20]

Neuroimaging studies of BD patients have generally implicated regions of the limbic network, including the amygdala and anterior cingulate cortex, and regions of the prefrontal cortex.[21] Functional imaging studies of BD have reported hyperactivity of limbic regions, in line with theories indicating a heightened emotional reactivity in the disorder.[21, 22, 23, 24] There is also evidence that unaffected relatives, who share a proportion of genetic material with the affected patient, display similar patterns of structural and functional abnormality to that seen in patients.[25, 26, 27, 28, 29, 30, 31, 32, 33, 34] Indeed, our previous study of familial high-risk individuals versus controls demonstrated a significant difference in amygdala activation between the groups using the same task.[19] Studies in relatives also allow the investigation of genetic effects against a genetic and (shared) environmental background similar to patients, but in the absence of psychotropic medication and other illness-associated confounding factors. In the current study we hypothesised, based on these previous findings, that the polygenic score derived for BD would be associated with increased activation of fronto-limbic circuitry underlying emotion regulation.

Methods and materials

Study population

The two groups examined in this study comprised a group of individuals at familial risk of BD and comparison group of healthy controls.[19] We note that individuals at familial risk of BD will also be at enhanced risk for other affective disorders and the findings should be viewed in this context. In brief, individuals at high genetic risk of BDI (BDI) because of a close family history of the disorder and control subjects were recruited as part of the Scottish Bipolar Family Study as has been described in detail elsewhere.[19, 35] To identify high-risk participants, caseloads of psychiatrists across Scotland were searched for individuals diagnosed with BDI. Diagnoses were confirmed with the Structural Clinical Interview for DSM-IV-TR Axis I Disorders (SCID-I)[36] or the symptom checklist of the Operational Criteria (OPCRIT).[37] Subsequently, subjects affected by BD were asked to identify a first or second-degree relative aged 16–25 years not suffering from the disorder. These unaffected individuals were invited to participate in this study providing that they had at least one first degree, or two-second degree relatives suffering from BDI. Control subjects with no personal history of BD or family history of a mood disorder in first-degree relatives were identified from the personal contacts of the bipolar high-risk subjects. Only unrelated individuals were included in the current analysis. Exclusion criteria for all groups included a personal history of major depression, mania or hypomania, psychosis, substance dependence, an IQ <70 or clinical diagnosis of learning disability, or any major neurological disorder or history of head injury that included loss of consciousness, and any contraindications to MRI. A total of 87 bipolar high risk and 71 controls provided genetic information, and a subgroup comprising 73 bipolar high risk and 52 controls provided suitable fMRI data. All participants provided written informed consent and the study was approved by the multicentre research ethics committee for Scotland.

Genotyping and derivation of polygenic scores

Genomic DNA was extracted from venous blood. Genotyping was conducted at the Wellcome Trust Clinical Research Facility, Edinburgh, United Kingdom (www.wtcrf.ed.ac.uk) and used the Illumina OmniExpress 730 K SNP array (San Diego, CA, USA). SNPs were excluded where the minor allele frequency was <1%, if the call rate was <95% or if the χ2-test for Hardy–Weinberg Equilibrium was <1 × 10−3. Strand ambiguous SNPs were also removed. The resulting SNP set was then used to calculate multidimensional scaling (MDS) components to assess for population stratification and adjust for this in later analyses. The data were then imputed to HapMap version 3 in mach software[38] and then converted back to plink (map/ped) format for later analysis.[39]

Summary results from the most recent international GWAS of 7481 individuals with BD and 9250 controls was obtained from PGC-BD. Details of the methods used by the consortium are given elsewhere.[7] Polygenic scores were calculated according to the methods by Purcell et al.[6] Four lists of significant SNPs were generated from the PGC-BD association data at significance thresholds of P=0.5, P=0.1, P=0.05 and P=0.01. These were then used to select common SNPs from our imaging-GWAS data sets resulting in four separate files that contained the genotype of each individual. Our primary analyses concerned those SNPs from PGC-BD that met a significance level of P=0.5 or less as this was the level that most efficiently discriminated individuals with and without BD in an independent test set.[6] Full details of findings for the other thresholds are presented in Supplementary Information (Supplementary Table 2). In order to identify polygenic effects due to independent SNPs in linkage equilibrium with one another, each SNP set was then pruned using a published method based on the variance inflation factor.[39] Linkage equilibrium-based SNP pruning was conducted in a sliding window of 50 SNPs with each calculation performed iteratively by moving the window by five SNPs in the 3′ direction. SNPs were conservatively selected on the basis of a variance inflation factor of 2 or less. Finally, these four SNP sets were then scored using the sum of the number of reference alleles multiplied by the natural logarithm of the odds ratio across the whole genome.

All analyses were performed in Plink[39] with the exception of imputation to HM3 and data manipulation, which were performed in Mach (www.sph.umich.edu/csg/abecasis/MACH) and R, respectively.[38, 40]

Clinical assessments

All participants were interviewed by one of two experienced psychiatrists (AMM, JES) using the SCID[36] to confirm the absence of any lifetime axis I disorders. Current manic and depressive symptoms were rated using the Young Mania Rating Scale[41] and Hamilton Rating Scale for Depression (HAM-D).[42] Estimates of temperamental variations in minor affective symptoms were assessed using the Temperament Evaluation of Memphis, Pisa, Paris and San Diego Autoquestionnaire (TEMPS-A),[43] providing measures of cyclothymic, depressive, hyperthymic, irritable and anxious temperament. Statistical analysis of demographic characteristics was conducted using independent t-tests or χ2-tests. For the clinical assessments and measures of temperament, comparison of groups was conducted using Mann–Whitney U-tests.

Experimental paradigm

Subjects performed the verbal initiation section of the Hayling Sentence Completion Test[44] in the scanner.[45] This is considered an extension of the verbal fluency task, where constraint is based on sentence context rather than letter or semantic category. Briefly, subjects were shown sentences with the last word missing and asked to silently think of an appropriate word to complete the sentence and press a button when they had done so. The task has four levels of difficulty, according to the range of suitable completion words suggested by the sentence context. This design allows a standard subtraction analysis (sentence completion versus baseline) and a parametric analysis (examining areas of increasing activation with increasing task difficulty). Sentences were presented in blocks of fixed difficulty. Each block lasted 40 s and included eight sentences. Sentences were presented for a period of 3 s followed by a fixation cross for 2 s. The baseline condition consisted of viewing a screen of white circles on a black background for 40 s. The order of the blocks was pseudo-random, and each block was repeated four times using different sentences.

Immediately after scanning, subjects were given the same sequence of sentences on paper and requested to complete each sentence with the word they first thought of in the scanner. ‘Word appropriateness' scores were determined from the word frequency list of sentence completion norms.[46] A score of one was given to the most frequently produced word in the word frequency list, a score of two for the next most frequently produced word and so on.

Scanning procedure

Imaging was carried out at the Brain Imaging Research Centre (BIRC) for Scotland on a GE 1.5 T Signa scanner (GE Medical, Milwaukee, WI, USA). The functional imaging protocol consisted of axial gradient-echo planar images (EPI) (TR/TE=2000/40 ms; matrix =64 × 64; field of view (fov)=24 cm) acquired continually during the experimental paradigm. Twenty-seven contiguous 5 mm slices were acquired within each TR. Each EPI acquisition was run for 404 volumes the first four of which were discarded. The T1 sequence yielded 180 contiguous 1.2 mm coronal slices (matrix=192 × 192; fov=24 cm; flip angle=8°). Visual stimuli were presented using a screen (IFIS, MRI Devices, Waukesha, WI, USA) placed in the bore of the magnet.

Image processing and analysis

The EPI and T1 images were reconstructed into nifti format (Mayo Foundation, Rochester, MN, USA) using DICOM convert functions available in SPM5 (Statistical Parametric Mapping: The Wellcome Department of Cognitive Neurology and collaborators, Institute of Neurology, London) running in Matlab (The MathWorks, Natick, MA, USA). Images were pre-processed using standard protocols available in SPM5. All EPI images were realigned to the mean volume in the series. Functional images were then normalised according to standard co-registration procedures using the individual's structural scan. Finally, all realigned and normalised images were smoothed with an 8 × 8 × 8 mm3 full width half maximum (FWHM) Gaussian filter.

First-level statistical analysis was performed using the general linear model approach. At the individual subject level the data were modelled with four conditions corresponding to the four difficulty levels each modelled by a boxcar convolved with a synthetic haemodynamic response function. Estimates of the subject's movement during the scan were entered as ‘covariates of no interest'. Before fitting the model, the participant's data was filtered in the time domain using high pass filter (128 s cutoff) and serial correlations were accounted for by using the first order autoregressive model. Contrasts were constructed to examine all four levels of sentence completion difficulty versus baseline, and areas of increasing activation with increasing task difficulty (the parametric contrast).

Second-level analysis

For each contrast of interest (sentence completion versus baseline and parametric effects), one contrast image per subject was entered into a second-level random effects analysis. The polygene score for each individual was entered into a full factorial model as a single regressor per group, hence adjusting for any main group effects. As described above the three MDS factors were entered as additional ‘nuisance' covariates in the analysis to control for population stratification. Contrasts were then constructed to examine polygene × group-interaction effects. Positive and negative associations across the groups between the polygene score and activation were then explored.

Statistical maps were thresholded at a level of P<0.001 (uncorrected) and regions were considered significant at a cluster level of P<0.05, corrected for multiple comparisons. All coordinates are quoted in Montreal Neurological Institute (MNI) convention (http://www.mni.mcgill.ca) and images are overlaid onto standard brain in MNI space using Mango software package (http://ric.uthscsa.edu/mango). On the basis of our prior hypothesis, small volume corrections were applied for the amygdala created using the WFU PickAtlas[47, 48] as in our previous study.[19]

Results

Demographic, clinical, temperament and behavioural measures

Subject details are presented in Table 1. There were no significant differences between the groups in terms of age, gender, handedness or NART IQ. The groups differed in terms of the polygenic score for BD, with the familial group scoring significantly higher than controls (P=0.04). For the clinical and temperament measures, the groups differed on measures of depression (as determined by the HAM-D, P=0.03), and cyclothymic temperament (P=0.03), with the high-risk group scoring highest for both measures as reported previously in this sample.[19] The polygene score did not correlate with NART IQ, or any of the clinical measures as described above, for the groups combined, or analysed separately. There were no significant differences in the within-scanner measures of reaction time or post-scan word appropriateness scores between the groups for those who underwent fMRI (demographics for these individuals are supplied as Supplementary Material). These measures also did not correlate with the polygene score, for the groups combined, or analysed separately. Both groups demonstrated the typical gradation in these behavioural measures according to task difficulty demonstrating compliance with the task.[19]

Table 1

Participants details

	Controls (n=71)		Bipolar high risk (n=87)		t/Z	P
Demographics
Mean age (years) (s.d.)	20.66	(2.40)	20.89	(2.88)	0.52	0.60
Gender (M:F)	36:35	—	43:44	—	0.03a	0.87
Handedness (R:L+mixed)	66:5	—	79:8	—	0.24a	0.62
Mean NART IQ (s.d.)	109.19	(7.42)	108.41	(9.00)	0.58	0.56
Polygene score
Polygene score	0.2290	(0.05)	0.2466	(0.06)	2.05	0.04
Clinical measuresb (median (interquartile range))
Young Mania Rating Scale score	0	(0)	0	(0)	1.05	0.30
Hamilton Depression Scale score	0	(1)	0	(2)	2.19	0.03
Temperament measuresb
TEMPS-A (median (interquartile range))
Cyclothymia	1	(3.5)	2	(3.5)	2.17	0.03
Depressive	0	(2)	0	(1)	0.06	0.95
Irritability	1	(2)	1	(2.5)	1.57	1.16
Hyperthymia	3	(3)	2	(3)	0.84	0.40
Anxious	1	(1.5)	0	(2)	0.55	0.58
Total score	7	(9)	8	(9)	0.77	0.44
Behavioural measures
Word appropriateness scores	3.04	(0.54)	2.91	(0.54)	1.38	0.17
Reaction time	2473	(619)	2558	(631)	0.75	0.46

χ2-test.

Mann–Whitney tests, median and interquartile range presented for skewed variables. Task-related behavioural measures only for those who underwent imaging protocol (n=52, n=73 for controls and high risk, respectively).

Task-related brain activation patterns

The controls demonstrated the expected patterns of brain activation and behavioural responses indicating subjects were performing the task appropriately in the scanner.[19, 45, 49] Regions activated for the sentence completion versus baseline contrast included the left medial and lateral prefrontal regions, left lateral temporal cortex, sub-cortical structures, left lateral parietal cortex, occipital lobes bilaterally and right cerebellum (see Supplementary Figure 1). Regions where there was reduced activation during the task (baseline greater than sentence completion) included the anterior cingulate cortex and amygdala. For the parametric contrast, all groups demonstrated activation in similar areas including left lateral and medial prefrontal cortex, left lateral temporal cortex and right cerebellum. The main patterns of activation and main group differences have been presented previously.[19]

Polygenic score and brain activation

Sentence completion versus baseline

For the sentence completion versus baseline contrast, there were no statistically significant polygene × group-interaction effects. There was, however, a significant positive association between the polygene score and activation in the anterior cingulate cortex (P<0.001, KE =477, Z=4.30, x=2, y=40, z=−6, r=0.32 corrected for multiple comparisons at the whole brain level,) and in the right amygdala (P=0.02, KE=40, Z=3.39, x=18, y=−6, z=−14, r=0.30 applying small volume correction for the amygdala) across the groups (see Figures 1 and 2). These polygenic effects were present over and above any main group effects, as group was taken into account in the statistical model. There were no regions of significant negative associations with the genetic score.

Figure 1

Associations between polygene score and activation. Depicts regions of increasing activation with increasing polygenic score for BD (a) in subgenual anterior cingulate and (b) amygdala. Images are overlaid onto standard brain in MNI space using Mango software package (http://ric.uthscsa.edu/mango). Map represents F-statistic images thresholded equivalent to P uncorrected=0.001, see methods for further details.

Figure 2

Scatter graphs of associations between the polygene score and activation in significant clusters. Depicts graphs of extracted data for clusters of significant association in (a) anterior cingulate and (b) amygdala. Triangles represent control subjects, cross hairs represent bipolar high-risk individuals. There were no significant groups × polygene score interaction effects. Overall Pearson's r-values were r=0.32 for the sgACC and r=0.30 for amygdala. For the groups separately, controls r=0.50, r=0.45, and for bipolar high-risk r=0.20, r=0.18 for the sgACC and amygdala, respectively.

Parametric contrast

For the parametric contrast, there were no significant positive or negative associations between activation and polygene score across the groups and no evidence for interaction effects between the groups.

Group effects controlling for polygene score

Based on our previous results demonstrating increased amygdala activation in the bipolar high-risk individuals versus controls, we further explored the data to determine whether there was any attenuation of the previously significant group differences. We found that although there were differences in amygdala activation occurring in the same direction as we reported in the previous paper,[19] these effects were attenuated and no longer statistically significant for either the sentence completion versus baseline, or for the parametric contrast. These results demonstrate that differences in polygenic load may have mediated our previous finding of increased amygdala activation in the high-risk group.

Population stratification

We also examined the relationship between the MDS components included in the model addressing population stratification and the extracted values from the clusters of activation associated with the polygene score, as above (the sgACC and amygdala). There were no significant correlations between these measures indicating that the above findings were not confounded by population stratification.

Discussion

In the current study we derived an overall polygenic risk profile score for BD in those at familial risk and healthy controls, and examined the collective impact of these risk variants on brain activation. We report a significant difference in polygene scores between the groups, where the familial group demonstrated higher profile scores for BD than controls. We did not find any evidence for interaction effects between polygene score and group, but we did demonstrate a significant positive association between an increasing risk score for BD and activation in the subgenual portion of the anterior cingulate cortex (sgACC, BA 32), and in the amygdala. Importantly, we were able to show these effects of overall risk on brain activation in the absence of confounding disease and medication effects. There was also no evidence for significant relationships between activation in these clusters and the MDS components indicating that these results were not significantly influenced by population stratification effects.

Neuroimaging and neuropathological studies have consistently implicated the sgACC in mood disorders, including BD.[21, 24, 50, 51, 52] Neuroimaging studies indicate differences in the activation patterns of the ACC in both early and later stages of the illness, and in those at familial risk.[53, 54] Functional imaging studies in healthy controls have suggested that this region is involved in the regulation of emotional behaviour,[50] and activation in this region has also been demonstrated to be associated with severity of depressive symptoms[55] and with sadness.[56] The sgACC has strong anatomical connections to the amygdala, which is similarly implicated in the pathophysiology of BD.[22, 23, 57, 58] The amygdala has an important role in emotional processing and in the attribution of emotion and motivational salience to stimuli. Amygdala dysfunction has been found in BD irrespective of mood state, in adulthood and adolescence, and in those at familial risk.[23, 59] Both of these regions have also been previously implicated in imaging studies of individual genetic loci associated with risk to the disorder, such as CACNA1C.[60, 61]

The sgACC and amygdala form part of a network traditionally implicated in neurobiological models of BD. These models suggest that in BD there is increased emotional reactivity along with an imbalance in the neural systems responsible for regulating emotional states.[22, 23, 62, 63] Regions of the ventral-limbic pathway, which includes the ACC and amygdala, are specifically proposed to be involved in the ability to identify the emotional significance of stimuli, the production of affective states and regulation of emotional responses. The complimentary dorsal system involving the prefrontal cortex is proposed to sub-serve more cognitive aspects of the regulation of affective states.[22, 23] It is suggested therefore that abnormalities in these functional systems controlling emotional reactivity might underlie deficits in emotional regulation and predispose individuals to development of the disease. The current findings are therefore highly consistent with both the literature and neurobiological models of the disorder and indicate an effect of cumulative genetic risk on these regions.

In our previous imaging study we reported that the bipolar high-risk cohort failed to show normal inhibition of the amygdala in response to increasing cognitive load.[19] The sgACC was also found to demonstrate task-related deactivation (see Supplementary Figure 1). The current findings therefore suggest that there is a failure to demonstrate normal inhibition of these regions in the context of a cognitive task, which relates to the degree of genetic risk for BD. A failure to inhibit activation of these limbic regions may therefore represent a heritable feature of the disorder and may underlie the abnormally sensitized response to non-emotional stimuli seen in the disorder. As we had previously reported a significant difference in amygdala activation between the groups,[19] we further explored the data to determine whether this difference would still be present if the polygene score was controlled for in the statistical model. We report that the group difference in the amygdala was no longer significant, strengthening the suggestion that there is an important contribution of polygenetic effects to amygdala activation in this context.

The polygenetic score for BD was significantly higher in our familial BD group than in the control sample, supporting the liability threshold model for BD. It should be considered, however, that because of the age of the individuals in the current study, the familial group will likely contain individuals who will remain well, along with individuals who will go on to develop psychiatric illness. We could speculate that high-risk individuals who subsequently become unwell may carry an increased number of risk alleles, whereas those who remain well carry protective alleles. Further longitudinal clinical investigation of these individuals will allow determination of these issues. Another important finding is that the polygene score was associated with activation in these ventral-limbic regions across both study groups with no evidence of significant interaction effects. The current findings suggest therefore that polygenic risk load exerts an effect on neural activation independent of group status. This indeed may have been anticipated as the polygene score is considered to capture a sizable proportion of variation in genetic liability.[6] This is also what might be expected under an additive model, as the control participants in our sample also carried a variable number of risk-alleles. On studying effects of individual SNPs, we previously reported significant gene by group-interaction effects in individuals at familial risk and controls.[13] This we interpreted as relating to the presence of an existing genetic background in the familial sample, consistent with the current findings.

In terms of future directions, polygenic risk profiling will enable further studies of psychiatric disorder and their shared genetic architecture with measures of brain structure and function, leading ultimately to the identification of shared variants and biological pathways. Using polygenetic risk profiling, it may also be possible to subsequently model complex psychiatric disorders in vitro, selecting individuals from the extremes of the distribution in order to identify the cellular process that underlie the imaging findings and disease biology.

The main limitations of the current study relate to sample size, which although substantial in terms of functional imaging, is relatively small in terms of genetic studies. Also, owing to the diversity of functional imaging paradigms applied across the field and due to practicalities in recruiting such a large familial sample, we were not able to perform a replication analysis of these findings. The findings should therefore be viewed considering these limitations.

In conclusion, we demonstrate that activation in regions known to be involved in the regulation of mood, and implicated in theoretical models of BD, are associated with overall genetic risk for the illness. These results suggest that hyper-activation in these regions may indeed be a trait-related feature of the disorder with an overlapping genetic architecture with BD. The advantage of this approach is that the polygenic score takes account of cumulative risk from multiple weakly associated genomic loci, and examines their overall association with neuroimaging measures, which to our knowledge have not previously been explored in BD. Previously, genetic models based on family studies would have determined risk to the disorder based on proximal distance to the primary affected family member hence all first-degree relatives would be assumed to have the same genetic risk. In the current study, however, polygenetic risk scores vary depending on each individual's profile of contributing risk alleles and this may therefore capture more of the phenotypic variation due to genetic factors than other approaches. The current findings therefore contribute to the understanding of deficits seen in patient groups by suggesting a mechanistic link between trait features of the disorder and genetic aetiology.