Introduction

Stress-related disorders, such as posttraumatic stress disorder (PTSD), are a function of altered hypothalamic–pituitary–adrenal (HPA) axis activity. Traumatizing events can negatively impact HPA responses,1 including enhancing glucocorticoid receptor (GR) sensitivity.2 Activating the HPA axis via stressors increases cortisol response in individuals with PTSD, while inhibition of arousal via dexamethasone administration results in enhanced cortisol suppression.3 It has been proposed that PTSD involves enhanced negative-feedback inhibition, which contributes to ongoing hyperarousal in PTSD patients.4

Relevant to the role of GRs in the stress response is the polymorphism FKPB5. The GR co-chaperone protein FK506 binding protein 5 (FKBP5) mediates GR activity, including the regulation of ligand binding, receptor activation and transcriptional processes.5 A number of FKBP5 polymorphisms have been found to have a role in stress responsivity. For instance, a history of childhood adversity interacts with rs1360780 to predict cortisol response.6 In response to a psychosocial stressor, healthy individuals with homozygous FKBP5 minor alleles have insufficient cortisol recovery and higher self-reported anxiety (for example, rs4713916, rs3800373, rs1360780, rs4713916, rs4713902 and rs3800374 (refs 7, 8)). The severity of childhood abuse also appears to predict adulthood PTSD via interactions with FKBP5 polymorphisms (for example, rs9470080, rs926158, rs3800373, rs1360780, rs9470080 (refs 9, 10)). From an epigenetic perspective, people with a risk allele of FKBP5 experience distinct chromatin conformations following childhood abuse, and this altered structure leads to greater GR receptivity.11 Taken together, these studies implicate an association between stress-related psychiatric conditions and FKBP5; however, the mechanisms underpinning this putative link have not been fully elucidated.

Few studies have investigated neural activation associated with FKBP5 allelic variations. One study found that one allelic variation of FKBP5 (rs1360780) was associated with altered hippocampal shape, greater attentional bias to threat and increased hippocampal activation in response to threat.12 It has also been reported that FKBP5 risk alleles interact with emotional neglect to predict dorsal amygdala reactivity to threat,13 FKBP5 gene expression is reduced in the amygdala in suicide victims,14 and PTSD patients have lower FKBP5 gene expression and decreased hippocampal and medial orbitofrontal cortex.15 Little is understood about how these risk alleles may affect neural functioning within otherwise healthy individuals.

Recent attention has turned to major inter-relationships between neural networks and how distinctive patterns of connectivity between networks underpin psychiatric conditions. Several key networks have been identified, and notably some of these are central to psychopathological stress responses, including PTSD. The salience network (SN), which monitors for and detects salient stimuli and involves the insula and amygdala, is overactivated in PTSD.16 The default mode network (DMN) includes the ventromedial prefrontal cortex (vmPFC), hippocampus and posterior cingulate cortex, and is linked with internal monitoring and autobiographical memory updating;17 further, PTSD is associated with reduced hippocampal and vmPFC activation suggesting disturbances in this basic regulatory state.18 Another major network is the central executive or frontoparietal (FP) network, and involves the dorsolateral prefrontal cortex and lateral parietal networks, and is associated with executive control and working memory.19

Resting-state connectivity provides a useful means to investigate the connections between neural networks because it is free of the potentially confounding influences of task-related activity and provides insight into spontaneous inter-relationships between brain regions.20 This approach focuses on low-frequency (0.01–0.10 Hz) fluctuations of the BOLD response during resting states—reflecting baseline activity of the brain.21 Resting-state studies of fear circuitry disorders, and in particular PTSD, have found altered connectivity in the DMN21 and the SN.22 Notably, there is decreased connectivity between these networks in PTSD.23

In this study, we investigated the extent to which FKBP5 moderates the connectivity within and between networks implicated in PTSD by assessing resting-state activity. To examine the impact of FKBP5 on resting state, we focus here on a fully healthy sample who vary in their distribution of FKBP5 risk alleles; this design permits delineation of the contribution of FKBP5 and the confounding influences of psychopathology on resting-state activity. Although different studies have reported various risk alleles of FKBP5, we focus on four single-nucleotide polymorphisms (SNPs) in the FKBP5 gene region that have been repeatedly associated with increased risk for enhanced stress responses. The selection of these specific alleles (and their justification) were: rs3800373 (increased FKBP5 protein levels, impaired HPA negative feedback following stressors, psychopathology risk following trauma;9, 24 rs9296158 (impaired negative feedback of the HPA axis, association with PTSD and depression;9, 25 rs1360780 (reduced basal cortisol levels, risk for PTSD and depression, attentional bias to threat;10, 12, 24, 26 and rs9470080 (lower basal cortisol level, increased risk for PTSD, and increased amygdala reactivity.9, 13, 27

Taken together, this accumulating evidence of FKBP5 risk alleles suggests that carriers of these alleles are more at risk of stress-related disorders, as well as the core cognitive and neural features associated with these disorders, because these alleles modulate GR sensitivity.28 On the basis that glucocorticoid system activity underpins stress and drives heightened emotional responsivity and considering that stress-related conditions are underpinned by distinctive resting-state networks, including the SN,16 we hypothesized that high-risk FKBP5 carriers would have a discernibly different resting-state profile in networks involved in emotional processing, such as the SN.

Materials and methods

Participants

Following written informed consent approved by the University of New South Wales Human Ethics Committee, 204 undergraduate psychology students of white European heritage provided saliva samples in return for course credit. Samples were collected using the Oragene DNA collection kit (DNA Genotek, Ottawa, ON, Canada) and genomic DNA was extracted. On the basis of associations between SNPs and indices of stress response,9, 10 (rs3800373, rs9296158, rs1360780 and rs9470080) were selected for genotyping. Genotypes were determined using iPLEX Gold primer extension followed by mass spectrometry analysis on the Sequenom MassARRAY system (Sequenom, San Diego, CA, USA) by the Australian Genome Research Facility (http://www.agrf.org.au/). Genotype frequencies in the total sample were as follows: rs9296158 AA=18, AG=88, GG=98; rs1360780 TT=19, CT=86, CC=99; rs3800373 GG =15, GT=85, TT=104; rs9470080 TT=24, CT =84, CC=96. The distributions of rs4713916, rs1360780, rs3800373 and rs9470080 did not differ from Hardy–Weinberg equilibrium (P=0.72, P=0.69, P=0.72 and P=0.66, respectively). Two risk groups were defined based on participants’ extreme FKBP5 risk haplotype values in order to test whether those with high risk versus those with no risk allele values demonstrated differential resting-state neural act. The low-risk group consisted of carriers that were heterozygous for zero high-risk alleles. The high-risk group represented carriers with four or more high-risk alleles, including a combination of hetero- and homozygous high-risk alleles. The final sample comprised 18 high risk (6 males and 12 females) and 19 low risk (7 males and 12 females). The sample size was determined on the basis of prior studies that have demonstrated differential resting-state patterns in groups of varying stress levels.23 On the day of scanning, these participants were administered the Beck Depression Inventory29 to assess depression levels, the Trait version of the State-Trait Anxiety Inventory30 to index trait anxiety and the Difficulties in Emotion Regulation Scale31 to index emotional distress and emotion dysregulation.

MRI data acquisition

One 5-min resting-state run for each subject was acquired using a Siemens Magnetom Trio 3T scanner located at the Advanced Research and Clinical High-field Imaging (ARCHI) facility of the University of Sydney. Functional scans were acquired transaxially using gradient-echo echo-planner imaging with the following parameters: repeat time 2.0 s; echo time 35ms; field of view 24 cm; acquisition matrix 64 × 64; flip angle 70°; voxel size 3.75 × 3.75 × 5 mm3; 29 axial slices; slice thickness 4 mm; gap 1 mm. Subjects were instructed to keep their eyes open during the functional scan and stare passively at a foveally presented fixation cross, as this is suggested to facilitate network delineation compared to eyes-closed conditions.32 Head motion during scanning was restrained using foam pads inserted on each side. In total, 155 whole-brain volumes of functional data including five ‘dummy’ scans were acquired. All participants were judged as awake and alert at the start and conclusion of scanning.

fMRI data analysis

Data pre-processing

Pre-processing was done using SPM5 software (http://www.fil.ion.ucl.ac.uk/spm/ software/spm5). Each subject’s functional images were first inspected visually for scanner artifacts and gross anatomical abnormalities, and then re-oriented so that the origin of the image lay within 3 cm of the anterior commissure. The initial 5 ‘dummy’ images were discarded to remove longitudinal equilibration effects. Images were then realigned using INRIAlign—a motion correction algorithm unbiased by local signal changes and corrected for slice time variation using the middle slice as the reference frame. Then these images were spatially normalized to a common stereotactic space using the Montreal Neurological Institute (MNI) EPI template and spatially smoothed with a Gaussian kernel of 10 mm3 full width at half maximum. Following spatial normalization, the data (originally acquired at 3.75 × 3.75 × 5 mm3) were slightly subsampled to 3 × 3 × 3 mm3, resulting in 53 × 63 × 46 voxels. Exclusion criteria regarding head motion was that subjects were excluded from analysis if their functional scans showed extreme head motion (translation>2 mm and rotation>1°). No participant was excluded from the analysis due to head movement.

ICA and identification of resting-state networks

Group spatial independent component analysis (ICA; Calhoun et al.33) was used to decompose data from all participants into components using the GIFT software (version 2.0e; http://mialab.mrn.org/software/gift/). The independent number of components in the data was determined using the ‘minimum description length’ (MDL) criteria. Twenty-six components was predicted by MDL, after which the data from all participants were concatenated and reduced to 26 temporal dimensions using standard principal component analysis, followed by an independent component estimation using the Infomax algorithm which attempts to minimize mutual information.33 The Infomax ICA algorithm wasrepeated 10 times in ICASSO (http://research.ics.tkk.fi/ica/icasso/) and resulting components were clustered to estimate the reliability of the decomposition.34 Following the group decomposition, single subject time courses and spatial maps (SMs) were then back-reconstructed using GICA and calibrated using z-scores. Components were visually inspected for artifacts.

Selection of networks of interest and visualization

Networks for analysis were chosen on the basis of three conditions: (1) they had to exhibit activations in gray matter combined with low spatial overlap with known vascular, ventricular motion and susceptibility artifacts; (2) be dominated by low-frequency fluctuations, that is, the ratio of integral of spectral power in low-frequency range (between 0.01 and 0.10 Hz) to high-frequency range (between 0.15 and 0.25 Hz) has to be greater than 5,35 and 3) represent the three core neurocognitive networks: the DMN, the FP network, and the SN36 by encompassing the primary regions of these networks. To visualize the spatial maps of a component, all subjects’ maps for that particular component were entered into a random-effect analysis model (1 sample t-test in SPM5). Brain regions were considered to be within each network if they met a height threshold of P<0.00001 corrected for multiple comparisons using the family-wise error (FWE) correction and an extent threshold of 50 voxels.

Second-level analyses of networks of interest

Translations and rotations of the head estimated during realignment stage of pre-processing were regressed out from the time course of each of the chosen networks in each subject and then these time courses were used to investigate differences between groups in three outcome variables: (1) network’s power spectra; (2) network’s SM; and (3) functional connectivity between networks (FNC). In order to focus on the subset of voxels which most likely represent the network and show strong and consistent activation across subjects component spatial maps were thresholded. Thresholding was based on the distribution of voxelwise t-statistics. Only voxels exceeding the threshold on the right tail: t>mean+4s.d. were selected to represent a network. These thresholded spatial maps were subsequently compared between groups. To determine differences between groups in network’s spectral power, a network’s time course was transformed to 129 spectral bins using fast Fourier transform and power in each bin was compared. For comparing differences between groups in between networks connectivity, first FNC was calculated following the procedure described by Jafri et al.37 For this, the time course data associated with the selected components were detrended, despiked and filtered usinga fifth-order Butterworth low-pass filter with a high-frequency cutoff of 0.15Hz.35 Then for each subject, correlations between pairwise combinations were calculated using Pearson’s correlations and resultant r-values were transformed to z-scores using Fisher’s transformation as a standardized measure of FNC. In order to determine differences between groups in SMs, spectral power, and FNC, MANCOVAN utility (within GIFT) was used.35 Primary aim of MANCOVAN is to uses a multivariate model selection strategy which reduces the total number of statistical tests performed and facilitates testing predictors on the response matrices (in this case it contains SM, power spectra and FNC) as a whole. It performs backward selection by testing whether each predictor in the model explains variability in the multivariate response using a multivariate analysis of covariance. The multivariate results of dimension-reduced responses, however, are difficult to interpret. Therefore, univariate tests were performed on each of the covariates of interest retained in the final model to understand the nature and extent of the relationship between covariates and response variables. An α level of 0.05 was used for all analyses and in addition to that univariate results were corrected for multiple comparisons using false discovery rate.38 The design matrix included high- vs low-risk allele group membership as a covariate of interest. A post hoc analysis was performed to investigate distribution of spectral power in the frontotemporoparietal (FTP) network over different frequency bands in the two allele groups. For this, entire frequency range (0.01–0.25 Hz) was divided into six equally spaced frequency bins (Bin1:0.01–0.04 Hz; Bin2:0.05–0.08 Hz; Bin3:0.09–0.12 Hz; Bin4:0.13–0.16 Hz; Bin5:0.17–0.20 Hz; Bin6: 0.21–0.25 Hz) and spectral power in each bin was calculated for each group. Significant differences between groups were determined by using two sample t-tests and correcting the results for multiple comparisons using false discovery rate.38

Results

Participant characteristics

Table 1 presents participant characteristics. There were no significant differences between low-risk and high-risk allele groups in age, depression, anxiety or emotional distress levels. Groups also did not differ in head motions.

Table 1 Participant characteristics
Full size table

Resting-state patterns

Out of the 26 components initially identified, six components were chosen for further analysis. In the present study, these networks were identified on the basis of their visual similarity to the networks identified from other resting-state studies (Jafri et al.;37 Allen et al.;35 Das et al., 201439, 40). These networks were highly stable (reliability index>0.89) as determined by ICASSO. Primary regions with each network are provided in Table 2. The first two networks represent the left and right FP networks (LFP and RFP), and the second set of three networks different manifestations of the core default mode network (DM1, DM2 and DM3).

Table 2 Primary brain regions within each network
Full size table

The sixth network contained the primary regions of the SN; specifically, the frontoinsular cortex and the dorsomedial prefrontal cortex/anterior cingulate cortex (DMPFC/ACC). This network also contained a temporoparietal region, therefore, in this study we referred to this network as a FTP network. Spatial maps of these chosen networks are shown in Figure 1.

Figure 1
figure 1

Six components/networks of interest. DM, default mode; FTP, frontotemporoparietal; LFP, left frotoperietal; RFP, right frotoperietal.

Full size image

The MANCOVAN yielded a significant effect (P<0.05) of group only in spectral power of the FTP network. Univariate tests also confirmed this finding after correcting for multiple comparisons. The FKBP5 high-risk allele group displayed less spectral power in the low frequency range than the FKBP5 low-risk allele group (Figure 2). Differences between groups in distribution of spectral power in the FTP network are shown in Figure 3. This analysis revealed a different pattern of fluctuation between the allele groups. Specifically, the low-risk allele group had significantly more power in the low frequency range (<0.1 Hz) compared to the high-risk allele group, who had significantly more power in higher frequency bins (>0.15 Hz). This difference was apparent only in the FTP network.

Figure 2
figure 2

Left panel depicts the significance and direction of group effect as a function of frequency for each component, displayed as the −(t)log10(p). Significant differences (P<0.01, corrected) between groups were observed only in the FTP network. It was found that in the low frequency range FKBP5 high-risk allele group had less spectral power than low-risk allele group. Right panel shows bar plots of the average β-values for group. β-Values were averaged over frequency bands with effects of the same directionality where test statistics exceeded the false discovery rate (FDR) threshold (P<0.05). The color of the bar is proportional to the fraction of contributing frequency bins; the absence of a bar indicates that multivariate tests did not identify this as an important predictor and therefore univariate tests were not performed or the univariate results did not survive significant threshold. FTP, frontotemporoparietal.

Full size image
Figure 3
figure 3

Distribution of spectral power of the frontotemporoparietal network over six equally spaced frequency bins (Bin1:0.01–0.04 Hz; Bin2:0.05–0.08 Hz; Bin3:0.09–0.12 Hz; Bin4:0.13–0.16 Hz; Bin5:0.17–0.20 Hz; Bin6: 0.21–0.25 Hz) in both high-risk and low-risk allele groups. High-risk allele group displayed significant less spectral power in low frequency range and greater power (but this did not survive multiple comparison correction) in high frequency range.

Full size image

Discussion

The major finding of this study was that carriers of the high-risk alleles of the FKBP5 polymorphism showed less spectral power than low-risk carriers, and specifically had more power in higher frequency bins in a resting-state network involving frontoinsular cortex, dorsal medial prefrontal cortex/anterior cingulate cortex and the temporoparietal cortex. Most noteworthy in this pattern is that these differences exist between participants who are otherwise healthy, do not differ in terms of any psychological indices related to emotional disturbance, and only differ according to their FKBP5 allelic status.

The current finding of a distinct resting state in the FTP network in those with different risk alleles of the FKBP5 polymorphism highlights a potential mechanism for why this genetic profile renders these individuals more at risk of PTSD following exposure to a traumatic event.9, 10 The mechanism by which FKBP5 is thought to render individuals with the risk alleles vulnerable to developing stress disorders when exposed to adversity is via impaired negative feedback of the glucocorticoid response, and thereby impoverished return to baseline of stress-induced cortisol levels.28 Altered desynchronization in resting-state networks that are implicated in salience monitoring, emotion processing and regulation may represent a vulnerability factor among those with FKBP5 risk alleles, such that when exposed to the stressful experiences of trauma, these networks function distinctively in association with glucocorticoid activation, thereby contributing to ongoing stress responses.

The finding that high-risk allele participants had less spectral power in the low frequency range than low-risk allele participants may suggest that the high-risk participants have less functional connectivity between regions in the FTP. Low-frequency spectral power is associated with functional connectivity in healthy participants,41 and so it is possible that carriers of the high-risk alleles possess less connectivity in the network involved in emotional processing. In this context, it is interesting to speculate how this finding may be interpreted in terms of current understanding of the implicated network. The DMPFC/ACC are both essential to emotion processing,42 and are active during appraisal of threat,43 attentional biases to threat44 and induced anxiety.45 The DMPFC is linked to other emotional and cognitive processing regions (including the ACC), suggesting that the DMPFC has a role in integrating emotion and cognition.46 One meta-analysis of brain networks implicated the role of DMPFC in emotional evaluation.47 The ACC has been associated across numerous studies with emotion regulation,48 and dorsal sub-sections are more active during appraisal and evaluation of emotional content.49 Impaired ACC functioning has been observed in emotion provocation tasks in PTSD populations.50 Both the frontoinsular cortex and the ACC are involved in the SN, which is implicated in detecting stimuli51 and activity is heightened in PTSD.16 The frontoinsular cortex is also associated with sustaining emotional awareness.52 Further, there is evidence of increased functional connectivity between the DMPFC and temporoparietal junction during social emotion processing.53 The temporoparietal cortex has been shown to be increased in response to threatening stimuli,54 and is involved with mentalizing,55 which can result in hypervigilence to threat evident in anxiety. Taken together, the distinctive FTP network resting-state pattern observed in the high-risk allele group accords with convergent evidence that these neural networks underpin threat detection, emotion processing and appraisal regulation.

In terms of potential limitations of the current study, we note the small sample size. We did not conduct a power analysis prior to the study simply because of the paucity of previous studies that have evaluated resting-state activity according to the risk alleles; instead we used as a guide the findings of distinct networks in stress-related conditions in comparative studies that had sample sizes smaller than ours.22 We also note that we focused on healthy participants, and so the findings cannot be generalized to clinical populations. However, it was our intention to initially isolate the effect of the FKBP5 risk alleles and in order to limit confounds; the lack of a clinical sample is a potential strength of this study in that the findings are not complicated by psychiatric morbidity and medication. We also note further clarification is required from studies that link resting-state analyses with fMRI studies that utilize provocation tasks. It would be interesting to integrate resting state and fear provocation paradigms, along with measures of cortisol, to determine the extent to which the FTP network identified in this study is associated with elevated neural fear responses in groups with FKBP5 risk alleles. Further, longitudinal studies are needed that assess the FTP network in participants with different FKBP5 risk alleles to index how they respond clinically and neutrally following exposure to traumatic events. We also note that ‘it is possible that because no strategies were employed to diminish potential physiological confounds such as the use of concurrent recordings of physiological changes and the use of removal strategies to exclude non-neuronal variables, these may ultimately have affected the findings. However, the use of an ICA framework, which manages to partially separate physiological noise is likely to have limited the impact of this potential confounds.35 These limitations notwithstanding, this study does provide novel insight into a potential neural basis for the increased risk for posttraumatic stress reactions in carriers of the FKBP5 risk alleles.

Recognizing that the actual mechanisms reflected by resting-state MRI are not well articulated,56 the current findings nonetheless point to distinct neural patterns in healthy people with risk alleles known to infer risk for fear-related disorders after exposure to adversity.