INTRODUCTION

Schizophrenia (SCZ) is a devastating brain disorder and severe form of psychotic illness with a complex phenotype presentation including hallucinations, delusions, thought disorder and blunted or reduced affect [1]. According to the stress-vulnerability model, onset and course of SCZ result from multiple factors [2]: On the one hand, there is an inherited risk of developing the disease as demonstrated by twin [3] and family studies [4,5]. Even though the contribution of individual allelic variants to the development of SCZ is not fully understood, genome-wide association studies (GWAS) have detected numerous common alleles of small effect [6] that overall explain up to half of the genetic variance [7]. On the other hand, a range of environmental factors can influence the expression of psychotic symptoms. In addition to cannabis abuse, growing up in an urban environment or belonging to a minority, exposure to early life adversity can also critically impair the development of the brain [8] and promote the development of later psychopathology [9]. For example, sexual abuse as a child, disruption of attachment relations or chronic victimization has been associated with hallucinations and paranoid symptoms in adulthood [10–12].

Structural MRI has augmented our understanding of SCZ as complex developmental brain disorder by longitudinal cortex mapping in adolescent patients with SCZ [13], or meta-analyses of grey matter (GM) morphology, demonstrating widespread cortical thinning and reduced volumes of subcortical structures, including the hippocampus, amygdala, thalamus and nucleus accumbens (NAcc) as well as reduced intracranial volume [14]. White matter changes are similarly multiregional, and like GM reductions and ventricular enlargement, intensify with longer disease duration [15,16]. Typical aberrant task responses localize to frontotemporal areas, including the dorsolateral prefrontal cortex (DLPFC), orbitofrontal cortex and left superior temporal gyrus [16]. Pivotally, these abnormalities are detectable in very similar areas during resting state fMRI (rsfMRI) that is acquired in the absence of any specific cognitive constraints. RsfMRI functional connectivity (FC) measures can be calculated by a broad range of methods directed towards both regional and multiregional FC properties [17]. A meta-analysis of heterogeneous rsfMRI studies suggests that medial PFC and anterior cingulate cortex are particularly sensitive to changes in SCZ [18]. Beyond these cortical areas, FC is also altered in subcortical areas, and particularly increases have been reported for the NAcc [19–22].

The NAcc is part of the mesolimbic pathway and the main projection site of dopaminergic neurons in the ventral tegmental area (VTA)[23] with additional inputs from other brain areas including the thalamus, hippocampus, amygdala and the PFC [24]. Its involvement in SCZ is closely linked to the dopamine hypothesis in which altered dopaminergic neurotransmission is conceptualized as causal to psychotic symptoms [25]: Mesolimbic dopamine projections are supposedly hyperactive, through presynaptic D2-(auto)receptors and postsynaptic D1-/D2-receptors, resulting in positive symptoms such as delusions and hallucinations; on the contrary, mesocortical dopamine projections, predominated by D1-receptors, are hypoactive, resulting in negative symptoms such as blunted affect and stereotyped thinking [25,26]. Of relevance to our approach, dopamine release in the NAcc, through several mechanisms, increases the local blood oxygen level dependent (BOLD) signal which implies that BOLD fMRI can be used to investigate dopaminergic dysregulation [27]. Increased FC of the NAcc to association cortices has been related to hallucinations in SCZ, especially their degree of complexity [28]. The NAcc is also a main target of antipsychotic drugs that block D2-receptors in this area [29] which makes this structure a primary region of interest.

Regarding genetic influences, ventral striatal seed FC has proven heritable in two childhood samples [30], with additive genetic influences ranging between 19% and 67% for six connected cortical regions. This provides a framework in which NAcc connectivity could be considered an intermediate phenotype (or potentially endophenotype [31]) that captures one (patho-)physiological process on which genetic and environmental factors could converge. Aggregating multiple risk polymorphisms into a polygenic risk score (PGRS) is one effective approach to parameterize the concurrence of many small genetic effects. Several studies have combined a PGRS of SCZ with task fMRI and rsfMRI (for review, see Dezhina et al. [32]). Here, more consistent findings were correlations between the PGRS and ventral striatal activation during choice behaviour [33] and reward processing [34,35]. Wang et al. [36] investigated PGRS effects on resting state FC of the insula, finding both positively and negatively associated regions. One study has linked specific dopaminergic risk allele load to putamen volumes and found an interaction with CA [37].

More recently, neuroimaging correlates of CA in patients with SCZ have been meta-analyzed based on 15 studies [38]: Patients with a history of childhood trauma showed lower total cerebral grey matter volume, particularly in the prefrontal cortex, reduced white matter integrity in specific tracts, and altered FC in a network comprising the amygdala, anterior cingulate cortex (ACC), precuneus/posterior cingulate cortex and the temporoparietal junction. FC analysis revealed an association between early life stress and altered amygdala/DLPFC and amygdala/rostral ACC connectivity that mediated stress response patterns in adulthood [39]. A relevant connection between early life stress and the mesocorticolimbic dopaminergic system has been detected using FC analysis of the VTA, identifying age-related attenuation of VTA/hippocampal connectivity [40]. Altered amygdala-striatal connectivity has been identified as a specific trait in male cocaine users with a history of CA [41], and impaired GABA signalling in the NAcc has been observed in animal models of early-life adversity [39].

Here we evaluated whether at the level of rsfMRI, CA and polygenic risk for SCZ have effects in the direction of a pattern seen in established SCZ, and if these influences interact. For this purpose, we applied two complementary FC analysis techniques: First, functional connectivity density (FCD) was calculated on a voxel-by-voxel: FCD, by an averaging or counting step, aggregates information from the entire FC network of a voxel back into a single scalar [42–44]. FCD maps of the whole brain were calculated and regional values from anatomically segmented NAcc area extracted. FCD is sensitive to anatomically diffuse or multiregional FC effects, yet due its averaging principle cannot localize the anatomical sources of its effects. Thus, NAcc seed analysis was added to backtrack anatomical sources of FC changes with the NAcc. Conceptually, both techniques are complementary and were employed in 253 healthy subjects to investigate the effects of polygenic risk for schizophrenia and its interaction with childhood adversity on resting state FC measures, particularly of the NAcc. Comparisons between healthy subjects and patients were performed in an endophenotype framework to replicate mainly FCD effects in schizophrenia and allow for effect size comparisons.

MATERIALS AND METHODS

The here presented fMRI data originated from four imaging genetics studies were acquired at the MPIP between 2014 and 2018. Study protocols were in line with the Declaration of Helsinki and approved by the local ethical review committee (IST: 121-14, 19 May 2014; TMEM: 174-13, 20 June 2013; BeCOME: 350-14, 9 June 2015 PSC: 084-14, 11 August 2014). The Imaging Stress Test (“IST”) study recruited 59 healthy subjects of which 53 (six genotyping QC failures) were included in this study. Second, an imaging genetics study directed towards effects of anxiety-related genes (ref. to as “TMEM”) recruited a total of 160 healthy subjects of which 136 were included here after ten general exclusions (e.g., incidental findings, abortion of the study protocol), 13 genetics-related (e.g., relatedness with other subjects, ethnicity outliers) and one rsfMRI specific dropout. The Biological Classification of Mental Disorders (“BeCOME”) study represents an ongoing deep phenotyping study of the MPI of Psychiatry that recruits healthy subjects and patients with currently symptomatic psychiatric disease mainly of the depression/anxiety spectrum. We included 38 healthy control subjects for which genotyping had been completed by the end 09/2018. The PsyCourse study (“PSC”) represents an ongoing multicenter German Research Foundation funded longitudinal biological psychiatric study for which healthy controls, patients with unipolar and bipolar depression as well as patients with SCZ are included [45]. From this study we included 26 healthy controls and 26 patients with a clinically confirmed diagnosis of SCZ (N = 15, hereof 3 secondary exclusions due to motion artefact) or schizoaffective disorder (N = 11). Inclusion and exclusion criteria of the four studies are listed in Supplemental Table S1.

All measurements were performed on a 3 Tesla scanner (Discovery MR750, General Electric, Milwaukee, USA) using a 32-channel head coil (MR Instruments Inc., Minneapolis, MN, USA). For each subject, functional resting-state whole-brain echo-planar imaging (EPI) data were acquired along with anatomical images and an AC-PC orientation of the slices. Three of four studies also included a single high-contrast EPI image. Subjects were instructed to lie still and fixate their gaze on a cross back projected on a monitor during the acquisition. Functional MR scanning parameters and length of the measurement varied slightly (Supplemental Table S2). In all studies, the rsfMRI acquisition was performed after anatomical sequences to allow for a general habituation phase regarding the MRI scanning environment. The anatomical high resolution image was uniformly acquired across all four studies (T1-weighted, sagittal FSPGR BRAVO, TR 6.2 ms, TE 2.3 ms, slice thickness 1.0 mm, 196 slices, field of view 25.6 × 25.6 cm², effective resolution 1 × 1 × 1 mm³).

All participants completed self-rating questionnaires including demographic information and childhood trauma experience. Childhood trauma was quantified with the German version 34-item Childhood Trauma Questionnaire (CTQ) [46,47] that classifies five types of adverse childhood and adolescence experience: physical abuse, emotional abuse, sexual abuse, emotional neglect and physical neglect. Each subscale consists of five items each rated on a five-level Likert scale. Severity classes per category can be built according to validated ranges [46,48,49]. In this study, we focused on the sum score of all five subscales. Besides these five subscales, CTQ includes a three-item minimization/denial validity scale to detect underreporting of maltreatment. In order to specify psychotic symptoms and quantify their severity, SCZ patients were clinically rated using the Positive and Negative Syndrome Scale (PANSS) [50].

Genotyping was performed using Illumina OmniExpress ([30] all TMEM participants and [51] all IST participants; separate runs) and Illumina GSA arrays [52] (all BeCOME, all PSC participants; one run). After genotyping, we applied a stringent quality control (QC) and removed all SNPs with a call rate below 98%, a minor allele frequency below 1% or a p-value for deviation from Hardy-Weinberg-Equilibrium below 1 × 10⁻⁵. We also excluded individuals with a call rate below 98%. Within each cohort, we performed multidimensional scaling analysis (MDS) on the LD-pruned genotypes in PLINK [53]. After this QC, we imputed missing genotypes using shapeit2 [54] and impute2 [55]. We used the 1000 Genomes Phase II dataset as reference sample. After imputation, we filtered for SNPs with an info-score above 0.8. In the final dataset, 8847825, 8256847 and 7734296 SNPs were available for batches (1)–(3). The SCZ PGRS was calculated based on a genome-wide association analysis for schizophrenia [6]. At first, best-guessed genotypes (based on a probability threshold of 0.9) were LD-clumped based on the genome-wide summary statistics, afterwards dosage data from the clumped SNPs were extracted and used for the calculation of the risk scores. The PGRS was calculated by summing up the risk alleles of those SNPs that presented genome-wide significance in the mentioned GWAS weighted with the respective log(OR). All risk score calculations were performed with PLINK. Resulting scores were scaled to the number of SNPs available for PGRS calculation (193 for TMEM and IST; 195 for BeCOME and PSC); scaled scores were then Z-transformed (PGRS_Z), separately for the three genotyping batches (1)–(3) to control for batch-specific ranking shifts. Regarding more lenient p^T such as 0.05, we considered the resulting PGRS not sufficiently reliable due to batch-related differences of the PGRS distributions also after the scaling step.

Preprocessing of fMRI data was based on inhouse software written in Matlab R2017a (MathWorks Inc., Natick, Massachusetts, USA) scripts combined with SPM (version SPM12) (https://www.fil.ion.ucl.ac.uk/spm/software/spm12) and FSL functions (https://fsl.fmrib.ox.ac.uk/fsl/fslwiki) that were identically applied to all four subsamples using a Linux cluster. The steps performed per subject were as follows: (1) The first four EPI volumes were excluded from the time series analysis to prevent T1-equilibration taking influence on the slice timing correction process. (2) Realignment was performed using rigid body motion correction with the first image as reference. Images were interpolated (5th spline interpolation). In addition, FSL based rigid body motion was performed on a copy of the timeseries to calculate root-mean-squared (RMS) intensity difference of volume N to volume N + 1 (ref. to as DVARS)[56] and defining outliers by default rules (values larger than 75th percentile plus 1.5 times the interquartile range). The resulting matrix of dummy regressors marking critical images of the time series was saved for later denoising. (3) Slice timing correction considering study specific number of slices, TR and slice order was performed using SPM12. (4) Spatial normalisation was supported by segmenting a specifically acquired strongly T2-weighted single EPI (for IST, TMEM and BeCOME) or an average of the first four T1-unequilibrated images (for PSC) using the Unified Segmentation algorithm of SPM12 [51]. For both types of segmentation driving input images, a good separation between GM, white matter (WM) and cerebrospinal fluid (CSF) could be achieved. Resulting GM and WM probability maps were entered as dual input into the iterative DARTEL normalisation algorithm (six generation of IXI templates; DARTEL default settings)[51] to produce a flowfield with warping information between native and MNI space. The time series was then spatially transformed using this transformation matrix and interpolated to a resolution of 3 × 3 × 3 mm³ (for FCD) or to 2 × 2 × 2 mm³ (for seed analysis). Spatially normalized versions of GM, WM and CSF probability maps, also at both resolutions, were generated. (5) A brain extraction mask was created (FSL brain extraction tool) and applied to the time series. (6) Denoising followed the CompCor strategy [57] in a multiple linear regression framework, including each five PCA components of WM (explained variance ~45%) and CSF (explained variance ~50%) in addition to six motion coefficients and their temporal derivatives, and the DVARS-based binary matrix. Resulting residuals were spatially smoothed (6 × 6 × 6 mm³) and added to the intercept for further analysis. (7) Temporal bandpass filtering (infinite impulse response filter)(range 0.008–0.100 Hz) was applied. (8) Eventually, the time series of all four substudies were temporally interpolated to a common TR of 2.25 ms using standard Matlab function (interp1, piecewise cubic spline). From this time series, the first five images and a varying (substudy specific) number or right-end images were skipped to gain a uniform 136 images (representing 5 min and 10 s) used as uniform time window. (9) Quality control involved (a) visualization of an exemplary preprocessed image of the time series overlaid with the normalized GM and WM segment, and (b) exclusion of subjects with more than 10 abrupt motions between two directly consecutive volumes exceeding 2 mm (translation) or 0.02 rad (rotation). The latter criterion marked a boundary line between the bulk of subjects and three SCZ patients with extreme outlier characteristics that were excluded.

High-resolution T1-weighted images (T1WI) were processed in two ways to be used in MNI space along with EPI time series: First, cortical reconstruction and volumetric segmentation was performed using the Freesurfer image analysis suite (http://surfer.nmr.mgh.harvard.edu, version 6.0). In parallel, DARTEL was applied to gain a transformation matrix to the MNI space. The discretized segmentation result image of FreeSurfer was then spatially transformed to MNI space and interpolated to 3 × 3 × 3 mm³ (for use on FCD maps) and 2 × 2 × 2 mm³ for seed analysis (see below). For both versions, an intensity cutoff of >0.6 was applied to re-binarize the mask after interpolation. Supplemental Figure S1 shows a thresholded (>0.6) group average of the MNI normalized NAcc of all study subjects at 2 × 2 × 2 mm³. The concept of fusing fMRI and T1WI pipelines in MNI space through DARTEL to combine mapping analyses in MNI space with regional analyses has proven useful before [58].

FCD accounts for the average FC between the index voxel and every other voxel in a predefined analysis volume. In-house Matlab software executed on a Linux cluster was used to calculate FCD based on the 3 × 3 × 3 mm³ preprocessed time series. The analysis space for each individuum was defined by a thresholded canonical GM probability mask (containing 47351 voxels) multiplied with the individual analysis mask of the residualisation model, resulting in an average of 43199 (SD 1590) voxels. Using intermediate storage of 500 × 500 cross correlation submatrices, the (absolute) FC values (Pearson correlation coefficient), after Fisher’s Z transformation, were averaged, resulting in one FCD map per subject. These maps were smoothed (FWHM 3 × 3 × 3 mm³) to increase the sensitivity and improve distribution properties of the residuals. Proof-of-concept group averages of controls and patients showed a posterior midline hub as expected (Supplemental Figure S2)[42]. A global FCD value (FCD_global; average of all voxels with a value >0.001 (identical with the 2nd level masking threshold)) was calculated for the later use as covariate. Further, average FCD (FCD_NAcc) of the left and right NAcc were recalculated from all individuals’ binarized 3 × 3 × 3 mm³ NAcc mask.

Seed-to-brain analyses were performed on the isometric 2 mm version to allow for higher accuracy of the seed placements and recruited networks. Before the actual NAcc analysis, a seed analysis using the anterior and posterior midline hubs of the default mode network (as gained from an independent sample of 25 subjects with EEG-validated wakefulness [59]) was performed in order to validate expected FC patterns. First level models contained an intercept and the regressor of interest (time series of the seed region), and resulting (unstandardized) beta coefficients were forwarded to a second level model covarying for age, gender and three source sample dummy variables. Inspection of the default mode map confirmed positively correlated areas, and anticorrelated areas in bilateral insular and inferior parietal areas, suggesting efficient global signal correction and mild, yet no inflated anticorrelations. Second level NAcc seed patterns revealed strongest FC in (orbito-)medial prefrontal cortex areas, as reported [59] (Supplemental Figure S3).

Group comparisons between patients and controls were performed for FCD maps and for seed based NAcc connectivity (left, right, bilateral): For FCD effects, age, sex and global FCD were used as nuisance covariates, in addition to three sample source dummy covariates. Age differed between patients and controls (Table 1, T = −5.654, p < 0.001) while no differences of the male/female distribution was found (Pearson Chi-Square 0.190, p = 0.663); no group difference was found for FCD_global after correction for age, sex and sample source (ANCOVA, F = 0.746, p = 0.388). Further, the relationship between focal (bilateral NAcc) and FCD_global differed between groups (partial correlation 0.532 for controls and 0.735 for patients, respectively; group-by-FCD_global effect [multiple linear regression] t = 6.966, p < 0.001). Therefore, at the voxel level, covariates age, sex and global FCD were modelled and centered separately per group. We analyzed the group main effect and, through conjunction analysis [60], combined group and group-by-covariate effects. For NAcc seed analyses, the procedure was identical, yet without a global correction variable.

TABLE 1

Table 1. Demographic and clinical characteristics of participants.

Generally, testing of polygenic risk or childhood adversity (CA) effects or interactions of these, were performed in the sample of healthy subjects (N = 253). The CA score (CTQ sum score) was log- and Z-transformed (CTQ_LOG-Z), and an interaction term (referred to as INTER_Z) calculated by multiplication with the PGRS_Z. Multiple linear regression (predictors age, sex, sample source, CTQ_LOG-Z, PGRS_Z, INTER_Z) was the general approach to investigate influences of CA, the PGRS and their interaction. For better readability, effects are referred to as CA, PGRS or their interaction in the following sections. Using SPSS (Release 18.0.0, SPSS Inc., Chicago, Illinois), this model was first applied to extracted FCD_global and FCD_NAcc (with FCD_global as covariate). FCD effects of the three regressors of interest were also analyzed at the voxel level (3 × 3 × 3 mm³) within a multiple linear regression model. Multiple test correction of clusters collected at p_voxel < 0.001 was performed for (i) a whole brain mask (51395 voxels), (ii) a study specific mask (FWE_voxel < 0.05, k = 10) gained from the SCZ/controls conjunction analysis (182 voxels), and (iii) for a bilateral NAcc mask (70 voxels). To understand, which functional NAcc connections are influenced by CA, genetic effects or their interaction, NAcc seed analyses were studied using the same second level setup.

To explore if FCD of the NAcc—when considered as endophenotypic marker of SCZ on which genetic and potentially environmental risk project—could be useful as risk indicator, we stratified controls into three risk levels, using a 75th percentile split point: low risk (PGRS and CA low), medium risk (either PGRS or CA high) and high risk (both PGRS and CA high). These groups were compared with the patient group (ANCOVA, 4-level factor risk, adjusted for age, sex, sample source, risk-by-age and risk-by-FCD_global).

For analyses in the anatomical space within the general linear model (GLM) framework of SPM, significance criteria were clusters collected at p < 0.001 followed by whole brain correction of cluster p-values for family-wise error (FWE). Significant result clusters (p_cluster.FWE < 0.05) were anatomically labelled using the Automated Anatomical Labelling (AAL) toolbox. For the analysis of FCD_global and regional FCD values of the NAcc, respective models were tested in SPSS (Release 18.01, SPSS Inc., Chicago, Illinois, USA) as described and p < 0.05 was accepted as significance criterium.

RESULTS

Table 1 summarizes demographic details of the healthy subjects (N = 253) and the patient sample (N = 23) included in the fMRI and genetic analyses. Participants of the original studies that lacked either CTQ information (incomplete or missing data), genetic data, or failed to pass the fMRI quality control are not listed. Diagnoses were defined by the attending psychiatrist during the hospital stay and represented the clinical diagnosis at discharge. The mean age differed between patients and controls by about 8.9 years (p < 0.001); no differences of the male/female distributions were detected. CTQ of the healthy subjects (N = 253) were mostly in the none to minimal severity class (63.6–94.1%), followed by low to moderate (2.4–26.5%), moderate to severe (2.0–7.1%) and severe to extreme severity (0.4–3.6%)(Supplemental Table S4). The distribution was very similar to one reported from a large German cohort [49] with about half of subjects scoring in the low to moderate severity class or higher in at least one of five categories. CTQ sum values were higher in the patient group by about 9 points (p < 0.001). The distribution of subjects over the number of CTQ categories with at least low to moderate severity (range 0 to 5 categories) differed between SCZ and controls (Pearson Chi Square 19.77, p < 0.001), confirming a CA/SCZ association [12]. Of the healthy subjects, 120 had an increased CTQ score in at least one category. Disease duration was about 12.4 years on average, though, with a large range between 0 (acutely diagnosed) and 43 years. 96% of patients were treated with antipsychotics (both atypical and typical) at the time of fMRI. Positive and negative symptoms in the SCZ patients were moderately high as expected for a sample with longer going disease under treatment.

Proof-of-concept Visualization of group average FCD maps were as expected from previous reports [42](Supplemental Figure S2). Higher FCD was found in five whole brain significant clusters collected at p_voxel < 0.001 (Supplemental Table S3) localizing to the right caudate, right posterior thalamus, and mid-/posterior cingulate and precuneus area (Supplemental Figure S4). Strong group-by-FCD_global effects (SCZ > controls) mapped to the bilateral subcortical nuclei including the NAcc and bilateral thalamus and extending to the anterior cingulate cortex (Supplemental Figure S4B). The global conjunction of both effects confirmed the bilateral subcortical nuclei including the NAcc, posterior cingulate/precuneus, right temporal cortical and mesencephalic clusters.

FCD_global analyzed as dependent variable showed effects of the sample source, yet no effect of CA, PGRS or PGRS-by-CA. Multiple regression analysis focusing on FCD_NAcc and corrected for FCD_global revealed significant effects for one sample source (IST study, T = 2.025, standardized beta (beta_s) = 0.114, p = 0.046), FCD_global (T = 9.760, beta_s = 0.5614, p < 0.001), the PGRS (T = 2.351, beta_s = 0.119, p = 0.020) and the PGRS-by-CA interaction (T = 2.453, beta_s = 0.126, p = 0.015). No effect of CA was detected (p = 0.215). Covariages age, sex, other sample source and MDS components were all not significant. The same four significant predictors emerged when a stepwise multiple linear regression was performed. Here, the PGRS increased the explained variance by 1.5% and the interaction term by 1.1%. An alternative statistical approach exemplifiesd the interaction effect: restricting the analysis to subjects with at least one rated childhood trauma domain (N = 120) resulted in an increase of the explained variance of the PGRS for the FCD_NAcc from 1.1% to 4.6%.

Explorative whole brain FCD analysis of PGRS effects revealed a single positively correlated cluster with k > 10 localizing to the right NAcc and caudate (k = 24, whole brain: p_cluster.FWE = 0.349; striatal mask: p_cluster.FWE = 0.015). The CA and PGRS-by-CA interaction regressors revealed no robust clusters. The global conjunction (positive contrasts on PGRS and interaction) revealed one cluster in the left NAcc (k = 31, whole brain: p_cluster.FWE = 0.076; striatal mask: p_cluster.FWE = 0.003).

NAcc seed analysis was used to dissect if the PGRS and CA-by-PGRS effects on FCD stem from specific brain areas. Here, the bilateral NAcc seed showed positive correlation with the PGRS in the right middle/superior occipital gyrus and middle temporal and parahippocampal gyrus (Table 2). The PGRS-by-CA interaction mapped to an occipital cluster (p_cluster.FWE = 0.024) exactly co-localizing with the PGRS effect. The conjunction of both contrasts mapped to six clusters with a bilateral occipito-parietal distribution that were non-overlapping with the normative positive and negative NAcc connectivity network (Figure 1).

TABLE 2

Table 2. NAcc seed analysis: effect of PGRS and PGRS-by-CA interaction.

FIGURE 1

Figure 1. PGRS effects at the FCD level and the NAcc seed analysis level. (A) FCD analysis: Note right striatal area with positive correlation with the PGRS (yellow) and left striatal area with conjunction effect (positive correlations with PGRS and PGRS-by-CT)(red), both significant within the volume-of-interest mask (see results section for details). (B) NAcc seed analysis: Combined PGRS and PGRS-by-CA effects (red) depicted against the background of the normative NAcc connectivity (positive NAcc seed correlation in white, negative NAcc seed correlation in blue, both p_voxel.FWE < 0.05). Both analyses revealed no main effect of CT.

As for the FCD_NAcc, also at the FCD and NAcc seed level, PGRS effects were investigated separately in subjects with at least one rated childhood trauma domain (N = 120) and the remaining group (N = 133, no trauma). A PGRS effect on FCD was only detected in the trauma group (k = 20, whole brain: p_cluster.FWE = 0.460; striatal mask: p_cluster.FWE = 0.030). Similarly, PGRS effects on NAcc seed connectivity were only found for the trauma group (Supplemental Figure S5).

FCD_NAcc values were higher in patients (F = 6.769, p = 0.010, partial eta² = 0.025) compared with healthy controls (ANCOVA with covariates age, sex, FCD_global and sample sources). Application of the same model using a 4-level factor risk (low (N = 145), medium (N = 91) and high risk (N = 17) in healthy subjects, and patients (N = 25)) showed a significant main effect of risk (F = 2.746, p = 0.043). Post-hoc pairwise comparisons confirmed differences between SCZ and the low (p = 0.009) and medium (p = 0.010) risk group, yet not between SCZ and high risk (p = 0.169) (Figure 2A). This pattern was confirmed when using the ratio between FCD_NAcc/FCD_global instead of a covariate approach (Figure 2B). In either model, re-separating the medium risk group into high CTQ or high PGRS did not change the result pattern.

FIGURE 2

Figure 2. FCD of the NAcc for three risk PGRS/CA risk levels compared with SCZ. (A) Bars represent mean (1 SEM) values of FCD of the bilateral NAcc adjusted for age, sex, source sample and FCD_global. (B) Same comparison, yet using the ratio between FCD_NAcc and FCD_global as dependent variable. Note intermediate position of the high risk group that overlaps with the SCZ group. Numbers printed on the bars represent the number of subjects within the respective group. Low risk: both CTQ and PGRS ≤ 75th percentile; medium risk: either CTQ or PGRS > 75th percentile; high risk: both CTQ and PGRS > 75th percentile.

DISCUSSION

We employed rsfMRI analysis to compare the effects of established SCZ with effects of polygenic risk for SCZ and its interaction with CA in healthy individuals. One benefit of combining case/control comparisons and imaging genetics analyses in one framework is the minimization of methodological sources of inconsistency that represent an acknowledged problem in rsfMRI. Our first set of analyses replicated previously reported increases in FCD in SCZ mainly in subcortical nuclei including the NAcc. In our second set of analyses, we found that higher polygenic risk for SCZ correlated positively with FCD in the bilateral NAcc and with specific NAcc connections in the secondary visual association cortices, and that polygenic effects were stronger in subjects with higher levels of CA.

Datasets with rsfMRI can be easily acquired and provide a rich source for analyzing spontaneously occurring functional patterns. Still, as the initial analysis space is high-dimensional (e.g., ~25.000 native space voxels with grey matter BOLD signal over ~150 time points, resulting in about 4.5 × 10⁶ raw data points) and as the space of possible (mostly complementary) analyses techniques is also large [61], a decision was made to use FCD as an entry technique. FCD aggregates connectivity information of one voxel to all other voxels in a single value and can be understood as sensitive, yet anatomically nonspecific screening technique [42,44,62]. In our analysis of 23 patients with SCZ or schizoaffective disorder, mostly treated with neuroleptics (96%) and with a moderate level of current symptoms, we found abnormally increased FCD in subcortical nuclei and posterior midline areas, and in addition, increased voxel/whole brain FCD correlations in the ventral striatum and mesencephalon. Conjunction analysis affirmed that particularly the NAcc presented such hyperconnectivity. Both the directionality and anatomical pattern matches two prior FCD reports on SCZ [36,22], even under consideration of technical differences: Global FCD (as used here) considers the connectivity between an index voxel and all other voxels, whereas local FCD restricts the voxels of interest to those adjacently connected to the index voxel [42]. Zhuo et al. [22,63] employed unweighted global FCD to compare 95 patients against 93 controls, counting voxels with a correlation coefficient larger than 0.6 and scaling these to the average count in GM. We emulated this counting approach post-hoc and found strong count increases in the subcortical nuclei (data not shown), supporting that our results are not specific to the decision weighted vs. unweighted FCD [62,36] focused on local FCD changes and mapped those to bilateral subcortical nuclei, (para-)hippocampal areas, and less prominently, to frontotemporal regions. Of note, long-range FCD (derived from subtracting local from global FCD) was increased in the bilateral thalami and pallidum, confirming that both short- and long-range connectivity of the basal ganglia are increased in SCZ. In terms of clinical specificity, both our sample and the one analyzed by Zhuo et al. [22,63] and Liu et al. [21] comprised patients treated with antipsychotics, not allowing to fully tease apart disease from medication effects. One study on drug-naive, first episode SCZ patients reported reduced FC between the NAcc and the anterior cingulate cortex [64], yet results on the full band width FCD are not reported. Another rsfMRI study on young, unmedicated patients with SCZ reported FCD effects from specific frequency bands, yet not analyzing the full band width [65].

Building on this replication of subcortical FCD increases in SCZ, we turned to the endophenotypic properties of FCD, particularly NAcc connectivity, in healthy subjects. By excluding patients, we prevented that clinical variability, such as differences in disease duration, acuity levels or medication, influences the imaging genetics results. In our healthy sample (aged 18–66 years), polygenic risk for SCZ, as calculated from a meta-GWAS based set of SNPs, correlated positively with FDC_NAcc, explaining about 1.1% of the variance. This is a first indication that rsfMRI NAcc connectivity is under specific genetic control and could reflect an intermediate pathophysiological level sensitive to inherited risk. Seed connectivity of the ventral striatum has proven heritable in two childhood samples [30] while heritability for FCD of the NAcc has not been determined to our knowledge. The association with the PGRS was weak, yet, it should be considered that the PGRS (when built on collection threshold of p^T < 0.05) explained about 7% of schizophrenia liability whereby half of the variance was attributable to the small group of whole genome wide significant SNPs. This strong weight of the most strongly associated SNPs had led us to construct the PGRS in this way although other thresholds and stratification techniques have been employed for imaging phenotypes. Regardless of the threshold, both the rsfMRI marker and the genetic marker are complex compound markers, integrating more than 100 genetic loci and FC of the NAcc to the whole brain. It is thus likely that correlations between specific NAcc connections and specific SCZ-related genes might be stronger. Filtering the genetic space for more NAcc-associated functions, e.g., genes regulating dopamine release and reception, could be a next step. Equally, focussing on the bivariate connectivity between the NAcc and the ventral tegmental area (VTA) could increase circuitry specificity.

The influence of marker selection on the fMRI side is also exemplified by a recent negative report on the association between a whole genome based PGRS for SCZ and rsfMRI [52,66]. Here, the marker was the Eucledian distance of an individual from the group’s centroid in terms of functional network connectivity in 50-dimensional space, i.e., the (temporal) relationships between time courses of resting state networks. This marker reflects brain dysfunction at a very high level of abstraction, yet proved sensitive to SCZ. Possibly, its disease specificity is lower than assumed, or differences in generating the PGRS (P^T < 0.0001) influenced the result.

We also searched the entire FCD voxel space, identifying an isolated correlated cluster of 24 voxels in the right NAcc/caudate that was not significant at the whole brain level, yet within the striatal mask, particularly within the group of subjects with at least one rated childhood trauma domain. Anatomically, our result resembles reports on PGRS effects on striatal responses during reward anticipation [34,35]. In conjunction with reported genome wide associations with SCZ at dopamine receptor D2 [6] and pathway analysis indirectly pointing to a role of internalization of the dopamine transporter during development [6,67], we speculate that fMRI is sensitive to ventral striatal functions indirectly reflected in the PGRS. NAcc seed analysis, in turn, revealed that the PGRS was positively correlated with NAcc connectivity to visual association cortex, in a pattern resembling the dorsal and ventral visual stream [68]. This resonates directly with reported increased connectivity between the NAcc and visual association cortex in SCZ patients with visuoauditory hallucinations [28]. It appears that PGRS correlates at the rsfMRI level do not only converge with increased FCD in SCZ, but also with NAcc connections specifically related to positive symptoms.

Beyond the direct genetic effect, we found that CA as measured by self-reports interacted with polygenic risk into the same effect direction. For FCD_NAcc the interaction term significantly increased the explained variance by 1.1%. At the NAcc seed level, the interaction term projected to a cluster overlapping with the genetic main effect which is noteworthy as both regressors were orthogonal. Reflecting the same pattern, the amount of variance explained by the PGRS differed strongly between subjects with no rated childhood trauma domain (N = 133) compared with subjects with at least one rated childhood trauma domain (N = 120). This pattern was also found at the FCD mapping and the NAcc seed level (Supplemental Figure S5). No main effect of CA was observed at the FCD or NAcc seed level which calls into attention that for mediation of some environmental factors, high genetic risk could be a prerequisite. We see a potential link between our result and epidemiological data that report hallucinations (of different types) to be disproportionally prevalent in patients with SCZ that have a history of childhood abuse [12]. This relationship exists across diagnostic borders (e.g., bipolar disorder) and also among non-patients: for example, bullying was found associated with visual hallucinations in the non-clinical general population [69]. A conceptually similar report revealed an interaction of CA with dopaminergic risk alleles on the volume of the left putamen in a combined SCZ/control sample [37]. So far, no other reports on interactions of CA with SZ-relevant genetic risk conditions exist for human functional striatal measures. Direct CA effects on the mesolimbic dopaminergic system as measured by rsfMRI have been found in 43 previously trauma exposed children and adolescents [40], and amygdala-striatal FC changes were reported in male drug cocaine users with a history of CA [41]. In experimental animals, early life-adversity was identified as factor influencing GABAergic functions of the NAcc [70]. The stronger polygenic effects we found in subjects with higher CA levels might point to epigenetic mechanisms that, as a consequence of CA, may lead to long-term changes of gene expression and subsequent functional adaptations.

Interesting analogies exists between subcortical hyperconnectivity (as mapped by FCD) in SCZ and cannabis abuse [71]. The epidemiological link between cannabis abuse, particularly when occurring at a younger age and over a longer period, and SCZ is a matter of a long debate [72–74]. Of value for this discussion, it was reported that patients with SCZ and cannabis use disorder showed reduced NAcc/prefrontal connectivity that was alleviated by cannabis [75]. These analogies suggest that NAcc or VTA based rsMRI connectivity markers, particularly when applied to sufficiently phenotyped population cohorts, may be useful to characterize individuals with regard to abnormalities in dopaminergic and thus SCZ relevant circuits. To exemplify the potential of FCD_NAcc as a “functional risk marker”, we quantified that a “double risk” subgroup of 17 subjects with both high (>75th percentile) polygenic risk and CA showed FCD_NAcc values in the range of the SCZ patients, whereas “no risk” or “single risk” was accompanied by FCD levels significantly different from patients. Serial studies, ideally in adolescents or young adults, are needed that investigate the actual predictive value of this marker.

A first and major limitation of our study is in heterogeneous characterization of the healthy subjects. For example, symptom scales that reflect the continuum between schizotypy, schizotypic personality and SCZ [76], cognitive tests, personality scores or other psychometric variables were not consistently available in all subjects. This hampered to probe with sufficient power if the identified fMRI risk pattern or combined CA/genetic risk translate to such liability traits as would be expected from a endophenotypic model. A second limitation is that due to the small size the heterogeneous patient sample, collapsing SCZ and schizoaffective disorder, more specific clinical correlates of NAcc connectivity could not be identified. In this line, a higher age in our patients complicated statistical modelling of group vs. group-by-age effects, particularly as some ventral striatal responses such as loss aversion are age-dependent [77,78] This theoretical limitation is attenuated to some degree as, statistically, no age correlation with global or voxelwise FCD in 253 healthy subjects was found, and as group comparisons mainly aimed at verifying a previously reported SCZ-related FCD pattern. Third, MR acquisition differences between our four substudies and different genotype batches may have reduced our statistical power in the imaging genetics part. Still, meta-analysis of FCDNAcc resulted in practically exchangabe effects compared with mega-analysis. Fourth, though post-acquisition motion correction was performed at usual standards, the degree of motion that needed correction was higher in the SCZ group, leaving the possibility of an influence on FC(D) group comparisons. Last, the assessment of CA using self-rating instruments may have introduced a memory bias in patients—yet, in this work, the reported CTQ difference between patients and controls has not influenced the main results from the healthy group.

CONCLUSIONS

We report that polygenic risk for SCZ explained a significant part of the variance of FC density of the NAcc in healthy subjects and that high self-reported CA intensified this relation. More specifically, increased connectivity between the NAcc and visual association cortices emerged as one underlying pattern. Subjects with high genetic and environmental risk approached FCD levels overlapping with SCZ, strengthening the role of NAcc FC as a potential endophenotype of SCZ. More generally, self-reported CA levels emerged as modifier of imaging genetic correlations which might help to understand result heterogeneity in imaging genetics.

SUPPLEMENTARY MATERIALS

Supplementary File 1 (PDF, 6.55 MB)

AUTHOR CONTRIBUTIONS

Conceptualization: PGS, CE; Methodology: PGS, DC, MC; Formal Analysis: CE, DC, PGS, BMM; Data Curation: CE, YP, PGS, FJ, DC, JM, VS; Writing-Original Draft Preparation: CE, PGS; Writing-Review & Editing: all coauthors.

ACKNOWLEDGMENTS

We are grateful to all participants of our studies, particularly patients recruited in the PsyCourse study curing their hospital stay at the hospital of the MPI of Psychiatry. We thank our MRI technicians Rosa Schirmer, Ines Eidner and Anna Hetzel for their longstanding support of MRI acquisition and data management throughout all four studies. We thank Alina Tontsch and the team of the Max Planck Institute of Psychiatry Biobank for sample processing. We thank all members of the BeCOME study team for developing the BeCOME study logistics and supporting the BeCOME recruitment and measurements. We thank Konstantin Strauch for valuable discussions of statistical aspects. We are grateful to Benno Pütz and Herbert Steininger for their enduring IT support.

FUNDING

CE, PGS and EB receive funding by the German Research Foundation (DFG, SA 1358/2-1 and Bl 806/6-1) for the PsyCourse study. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

CONFLICTS OF INTEREST

The authors declare that there is no conflict of interest regarding the publication of this paper. The authors declare that they have no competing interests.

How to Cite This Article

Eberle C, Peterse Y, Jukic F, Müller-Myhsok B, Czamara D, Martins J, et al. Endophenotype Potential of Nucleus Accumbens Functional Connectivity: Effects of Polygenic Risk for Schizophrenia Interacting with Childhood Adversity. J Psychiatry Brain Sci. 2019;4:e190011. https://doi.org/10.20900/jpbs.20190011