Multiscale neurobiological correlates of human neuroticism.

Neuroticism is a heritable personality trait associated with negative emotionality; however, we know little regarding the association between the microscale and macroscale neurobiological substrates of human neuroticism. Cross-scale correlation analysis may provide such information. In this study, voxel-wise neuroimaging-neuroticism correlation analyses consistently showed a positive correlation between neuroticism and functional connectivity density (FCD) in the ventral striatum in 274 young Chinese adults. Partial least squares regression analysis showed that the FCD-neuroticism correlation map was significantly spatially correlated with gene expression profiles in each of six donated human brains. Neuroticism-related genes derived from the six donors consistently showed significant enrichment in the chemical synaptic transmission, circadian entrainment, long-term potentiation, inflammatory mediator regulation of transient receptor potential channels, and amphetamine addiction pathways. The protein-protein interaction analysis revealed four hub genes involved in the above pathways, including G protein subunit gamma 10, 5-hydroxytryptamine receptor 2C, prodynorphin, and calcium/calmodulin-dependent protein kinase II alpha. By combining multiscale correlation analyses and functional annotations, this study advances our understanding of the genetic and neural substrates of human neuroticism and emphasizes the importance of striatal functional properties in human neuroticism.

INTRODUCTION

Neuroticism is a heritable personality trait reflecting an individual's negative emotional reactivity and has been frequently related to mood disorders (Mincic, 2015). An investigation of the neurobiological substrates of human neuroticism may shed light on the neurobiological predispositions for neuroticism‐related mood disorders. The neurobiological substrates of human neuroticism have been explored from two aspects. First, candidate gene studies and hypothesis‐free genome‐wide association studies (GWAS) (microscale‐personality) have identified hundreds of neuroticism‐related genes and genetic variants. Early candidate genetic association studies have identified several genes and genetic variants associated with neuroticism, such as serotonin transporter polymorphism, brain‐derived neurotrophic factor, and catechol‐o‐methyltransferase (Sanchez‐Roige, Gray, MacKillop, Chen, & Palmer, 2018; Sen et al., 2003; Stein, Fallin, Schork, & Gelernter, 2005). More recent GWAS and meta‐analyses of neuroticism have found more consistent significant loci on chromosomes 1, 8, 9, 11, 15, 17, 18, and 22 and in neurodevelopmental process and neurotransmitter pathways (Nagel, Jansen, et al., 2018). Second, neuroimaging–neuroticism association studies (macroscale‐personality) have revealed significant correlations between neuroticism and structural and functional properties of the human brain mainly in brain regions involving emotional processing, including limbic regions (amygdala and hippocampus) and neocortical regions (temporal, parietal, and occipital brain areas) (Canli et al., 2001; Servaas et al., 2013; Wright et al., 2006). However, we know little regarding the association between the microscale and macroscale neurobiological substrates of human neuroticism, which is critically important to create neurobiological pathways of microscale–macroscale‐neuroticism.

At least two strategies can be used in transcription‐neuroimaging association studies. One is to investigate the consistency of interindividual variations between transcription and neuroimaging profiles and the other is to study the consistency of interbrain region variations between them (spatial correlation across brain regions). The former requires a very large sample of subjects with both transcriptional and neuroimaging data, but no such data are available. The latter can be performed in a single individual (such as any one of the Allen Human Brain Atlas [AHBA] donors) who has both neuroimaging and gene expression data with enough spatial coverage. Based on the fact that some genes show high consistency in cortical regional expression patterns across donors in the Allen microarray data (Fu et al., 2020; Hawrylycz et al., 2015), the spatial expression patterns of these genes in any human brain can be approximately represented by those in AHBA donated brains (Liu et al., 2019). At least for these conserved genes, we can investigate the spatial correlation between gene expression from the AHBA donors and neuroimaging metrics from other individuals. In this study, we used a measure of spatial expression stability (SES) to identify genes with conserved transcriptional profiles and performed leave‐one‐donor‐out cross‐validation analysis to validate the transcription‐neuroimaging association results. Although our method may miss genes with great interindividual variation but that are associated with neuroticism‐related neuroimaging phenotypes, our analyses can still provide insights into how spatial patterns of gene expression relate to regional variations of neuroimaging phenotypes to bridge the gap between macroscale neuroimaging measures and microscale transcriptional profiles in terms of neuroticism. By investigating sample‐wise spatial correlations between gene expression in the postmortem human brains from the AHBA and neuroticism‐related neuroimaging statistical map from a different imaging dataset, one can identify genes whose transcriptional profiles are associated with neuroticism‐related macroscale neuroimaging properties.

In this study, we performed comprehensive analyses to investigate the neurobiological substrates of human neuroticism. We first performed voxel‐wise neuroimaging–neuroticism correlation analyses to investigate the associations of FCD and gray matter volume (GMV) with neuroticism scores in 274 healthy Chinese Han young adults. FCD mapping was selected because it is ultrafast, comprehensive, and hypothesis‐free method for quantitatively assessing functional connectivity at a voxel‐wise level. Then, we performed genome‐wide SES analysis to identify genes with highly consistent transcriptional patterns across brain regions in six AHBA brains. Third, a multivariate partial least squares (PLS) regression analysis was used to investigate the consistent transcription‐neuroimaging associations to identify genes whose transcriptional profiles were significantly associated with neuroticism‐related neuroimaging properties. Finally, we performed a series of gene enrichment analyses to further characterize the biological functions of the significant PLS genes. We expected to identify new biological pathways and hub genes associated with human neuroticism and to provide new insights on the neurobiological substrates of human neuroticism.

METHODS

Subjects

We recruited 323 healthy Chinese Han young adults with strong right‐handedness (157 males and 166 females; age: 22.7 ± 2.5 years, range: 18–31 years) under the approval of the Medical Research Ethics Committee of Tianjin Medical University General Hospital. All recruited subjects signed an informed consent form and were strongly right‐handed assessed as the revised Chinese edition of the Edinburgh Handedness Inventory (Oldfield, 1971). The final imaging analysis included 274 subjects (126 males and 148 females; age: 22.8 ± 2.4 years; years of education: 15.6 ± 2.1 years) after excluding 49 subjects who did not fill out a personality assessment form (n = 23) or had bad image quality, including metal artifacts (n = 8), brain abnormalities (n = 2), or excessive head motion (n = 16). Other exclusion criteria included histories of psychiatric or neurological illness, drug or alcohol abuse, and contraindications for MRI examinations.

Personality assessments

Negative emotionality traits (Mincic, 2015), including neuroticism, trait anxiety, harm avoidance and depression were assessed by the Eysenck Personality Questionnaire (EPQ) (Eysenck, 1991), the State–Trait Anxiety Inventory (Spielberger & Gorsuch, 1983), the Tridimensional Personality Questionnaire (Cloninger, Przybeck, & Svrakic, 1991), and Beck Depression Inventory‐II (Beck, Steer, & Brown, 1996), respectively. We included the assessments of trait anxiety, harm avoidance, and depression to reduce the influence of these traits on neuroticism‐related analysis because of the presence of strong correlations between neuroticism and these traits (Mincic, 2015). Another personality trait‐extraversion, which is thought to be an independent personality dimension according to Eysenck's personality theory, was also assessed by the EPQ.

We transformed the raw neuroticism score into a T score (T = 50 + 10 × (raw score‐mean)/SD) according to the Chinese population norm of the EPQ (Qian, Wu, Zhu, & Zhang, 2000). The T scores of neuroticism were used in the following neuroimaging‐personality association analyses.

MRI data acquisition

All subjects were scanned using a 3.0‐Tesla MR system (Signa HDx, General Electric, Milwaukee, WI). Resting‐state functional images were acquired with a single‐shot gradient echo planar imaging sequence. The parameters were as follows: repetition time (TR) = 2,000 ms; echo time (TE) = 30 ms; field of view (FOV) = 240 mm × 240 mm; matrix = 64 × 64; flip angle = 90°; slice thickness = 4 mm, no gap; 40 interleaved transverse slices; and 180 volumes. Sagittal 3D T1‐weighted images were collected using a brain volume sequence with the following scan parameters: TR/TE = 8.1/3.1 ms; inversion time = 450 ms; flip angle = 13°; FOV = 256 mm × 256 mm; matrix = 256 × 256; slice thickness = 1 mm, no gap; and 176 sagittal slices. All participants were instructed to move as little as possible, keep their eyes closed, think of nothing, and refrain from sleeping during the fMRI scans. After scanning, a questionnaire was used to assess their waking state during the fMRI scans. If the participant was not awake, the fMRI data were discarded, and the participant was scanned again.

fMRI data preprocessing

The resting‐state fMRI data were preprocessed using SPM12 (www.fil.ion.ucl.ac.uk/spm) with the following steps. The first 10 volumes of each functional time series were removed to allow the signal to reach equilibrium and the participants to adapt to the scanning noise. Acquisition time delays were corrected between slices for the remaining 170 volumes, and then these volumes were realigned to the first volume to correct head motions. We further excluded subjects with translational or rotational motion more than 2 mm or 2°, respectively. Moreover, we calculated the frame‐wise displacement (FD), an index indicating volume‐to‐volume changes in head position (Power, Barnes, Snyder, Schlaggar, & Petersen, 2012). We obtained the FD from the derivatives of the rigid body realignment estimates from the realignment of head motion (Power et al., 2012; Power, Barnes, Snyder, Schlaggar, & Petersen, 2013). These functional images were spatially normalized to the Montreal Neurological Institute (MNI) standard space using two‐step normalization. We first coregistered individual structural images (high‐resolution T1‐weighted images) to the mean motion‐corrected functional image by linear transformation. The coregistered structural images were then segmented into gray matter (GM), white matter (WM), and cerebrospinal fluid (CSF), and the GM was nonlinearly coregistered to the MNI space. Finally, the motion‐corrected functional volumes were normalized to the MNI standard space by using the parameters estimated during the nonlinear coregistration of the structural images. Subsequently, we resampled the normalized functional volumes into a voxel size of 3 × 3 × 3 mm3. At last, functional images were band‐pass filtered with a frequency range of 0.01–0.1 Hz and several nuisance covariates (including six parameters of head motion and the average BOLD signals of the CSF and WM) were regressed out.

FCD calculation and analysis

The FCD value of each voxel was calculated using an in‐house script according to a method described previously (Tomasi & Volkow, 2011). Functional connections were calculated using Pearson's linear correlation, a functional connection was identified if the correlation coefficient between any two voxels was greater than 0.6 (Tomasi & Volkow, 2012). The FCD of a specific voxel x 0 was defined as the total number of functional connections between x 0 and all other voxels. To minimize unwanted effects from susceptibility‐related signal‐loss artifacts, we applied a GM mask to restrict the calculation of the FCD to voxels only in the GM regions with a signal‐to‐noise ratio >50 (Tomasi & Volkow, 2012). To improve the normality of the distribution, the grand mean scaling of FCD was obtained by dividing by the mean FCD value of the qualified voxels in the whole brain. Finally, the normalized FCD maps were spatially smoothed with a 6 × 6 × 6 mm3 full‐width‐at‐half‐maximum (FWHM) Gaussian kernel. In our study, to keep the results reliable, we also used threshold r‐values = .5 and .7 to calculate the FCD.

GMV calculation and analysis

All structural images were carefully checked slice by slice to ensure image quality and exclude gross anatomical abnormalities. The structural MR images were segmented into GM, WM, and CSF using SPM12, which is thought to reveal the density or concentration of GM, WM, and CSF (Ashburner & Friston, 2000). Following segmentation, a GM template was generated from the entire image dataset using diffeomorphic anatomical registration through the exponentiated Lie algebra technique (Ashburner, 2007). After initial affine registration of the GM template to the tissue probability map in the MNI space, the GM images were nonlinearly warped to the GM template in the MNI space with a resolution of 1.5 mm3. The nonlinear components derived during spatial normalization were used to modulate the GM value of each voxel. According to the manual of the VBM8 toolbox (http://www.neuro.uni-jena.de/vbm/segmentation/modulation/), GMV is defined as the modulated output of the GM value of each voxel corrected for intracranial volume. Finally, to compensate for residual between‐subject anatomical differences, the GMV images were smoothed using a Gaussian kernel of 4 × 4 × 4 mm3 FWHM. Then, the preprocessed GMV maps were used for further analysis.

Personality‐imaging (FCD and GMV) association analysis

With SPM12 software, multiple regression analysis was used to investigate voxel‐wise correlations between the FCD values in the whole brain GM and personality scores (neuroticism and extroversion, respectively) at a connection threshold of r = .6 while controlling for the effects of mean FD, age, sex, educational years, trait anxiety, harm avoidance, and depression. To make the results reliable, we used the same model of multiple regression analysis between personality scores and FCD for two additional connection thresholds: r = .5 and .7. Because trait anxiety and harm avoidance are most likely highly correlated with neuroticism, to test whether removing these parts of the variance would bias the results, we performed spatial correlation between the uncorrected voxel‐wise FCD‐neuroticism correlation map after controlling for trait anxiety and harm avoidance and the uncorrected voxel‐wise FCD‐neuroticism correlation map without controlling for trait anxiety and harm avoidance. In addition, we validated the FCD‐neuroticism regression analysis by using Human Connectome Project (HCP) data from the WU‐Minn HCP Consortium. Detailed information on data acquisition and preprocessing can be found in the supplementary materials.

For personality score‐FCD correlation analysis, we obtained the significant clusters related to personality scores. Then, we used connection probability maps to show the functional connectivity patterns of these significant clusters. The connection probability maps were generated individually for each voxel in the significant clusters and then calculated to generate probability maps at the group level. Specifically, all the preprocessing steps were the same as those for the FCD calculation. Then, at the individual level, for a significant cluster in the FCD‐personality analysis, we generated a functional connection map for each voxel in this cluster for each subject and binarized each connection of each voxel using the same connection threshold as that used for the FCD calculation. Then, at the group level, we calculated the number of connections of each voxel at the same connection threshold and divided by the number of all subjects to generate a connection probability map of each voxel in this cluster for this connection threshold. The connection probability maps could provide information on which functional connections may contribute to the FCD‐personality correlation.

In addition, we also performed voxel‐wise regression analysis to determine the association between GMV and personality scores while controlling for the effects of age, gender, years of education, trait anxiety, depression, and harm avoidance. Nonstationary cluster extent correction was utilized to correct inhomogeneities of local smoothness which might result in misestimating cluster size, producing false‐positive results (Kurth, Gaser, & Luders, 2015). Multiple comparisons were corrected using the same method and threshold as in the FCD analysis (i.e., voxel‐level uncorrected p < .001 and cluster‐level corrected p < .05).

Brain gene expression dataset

The Allen Institute for Brain Science provides the normalized microarray expression data for six donated human brains (five males and one female) with an average age of 42.5 years (https://human.brain-map.org/). A total of 3,702 tissue samples were collected from the cortex, subcortex, cerebellum, and brainstem in the six brains and only cerebral samples (cortex and subcortex) were included in our analysis. Each tissue sample was measured using a custom Agilent 8 × 60 K array chip, containing 58,692 probes for more than 20,000 genes. More importantly, each postmortem brain underwent structural MRI scanning prior to dissection to determine anatomical structures, which was registered to the MNI standard coordinate space, and thus we could obtain MNI coordinates for each sample, providing an opportunity to bridge the gene expression dataset and imaging dataset. Although the AHBA did not include resting‐state fMRI data and personality assessments, we were able to investigate sample‐wise spatial correlation between gene expression data from the AHBA and personality‐related imaging metrics from our imaging dataset. Details of the data normalization, probe filtration, annotations from probes to genes are provided in the supplementary methods.

SES calculation

Differential stability has been proposed as a correlation‐based metric for measuring the reliability of the expression patterns of the differentially expressed genes across brain structures in the six donated brains (Hawrylycz et al., 2015). Using the same line of thinking, we proposed SES to measure the reliability of the spatial expression patterns of the genes of interest across brain structures in the six donated brains. For a certain gene, we calculated the Pearson correlation coefficient of gene expression across common brain regions between each pair of donated brains; and then defined the SES of this gene as the average Pearson correlation coefficient over all 15 donated brain pairs. The bootstrapped 95% confidence interval (CI) of SES was used to estimate the probability of the true mean. The common brain structure was defined by using an automatic anatomical labeling (AAL) template (Tzourio‐Mazoyer et al., 2002); we matched tissue samples to AAL regions by their MNI coordinates, and averaged the gene expression values of the samples in the same AAL region. Then we ranked the genes by their SES values from high to low and used the top half of the high‐SES genes to perform the following PLS regression analysis.

In the initial analysis, we used SES = 0.6 to obtain the genes with stable spatial expression across the six donors, and then we used two different SES cutoff thresholds (0.5 and 0.7) to ensure that our initial results were reliable. Last, we performed leave‐one‐donor‐out cross‐validation analysis to validate the SES calculation and the following analysis. Details of the validation of the SES calculation are provided in the supplementary methods.

Transcription‐imaging association using PLS regression analysis

PLS regression, a multivariate regression analysis that can predict a set of dependent variables from a set of independent variables, was used to investigate the association between the gene expression of high‐SES genes and personality‐related imaging metrics. PLS regression analysis is best applied when the number of predictors obviously exceeds the number of observations. However other multivariate methods, such as canonical correlation analysis (CCA) and multiset CCA (MCCA), require that the number of predictors cannot exceed the number of observations, because CCA can obtain weight vectors that are not uniquely defined, leading to overfitting issues (Grellmann et al., 2015). Moreover, PLS regression is also applied for highly interdependent or multicollinear independent variables (Krishnan, Williams, McIntosh, & Abdi, 2011). Here, we used the SIMPLS algorithm (De Jong, 1993) to investigate the predictor matrix and response vector, which was also utilized by a previous study (Romero‐Garcia, Warrier, Bullmore, Baron‐Cohen, & Bethlehem, 2019). PLS regression was used to identify significant genes contributing to the personality‐related imaging metrics. Generally, the main steps of PLS regression analysis include constructing predictor matrix and response variable vectors, regression analysis, and statistical significance assessment of the PLS components. Below we take the neuroticism‐FCD correlation T map as an example to describe the steps of the analysis.

First, the matrix of the sample‐wise high‐SES gene expression in the AHBA dataset (number of samples in each donor with gene expression data × number of high‐SES genes) were defined as the independent variables, and the vector of sample‐wise T values of the FCD‐neuroticism correlation map in our imaging dataset (number of samples in each donor with gene expression data × one neuroticism‐FCD T map) was designated the dependent variable. There is only one voxel‐wise FCD‐neuroticism correlation map generated by regression analysis between FCD and neuroticism in the 274 subjects. In this map, each voxel has a statistical t‐value reflecting the correlation strength between the FCD values of this voxel and the neuroticism scores in the 274 subjects. For each donated brain, each tissue sample has an MNI coordinate in standard space. We then drew a spherical region of interest (radius = 6 mm) centered at the MNI coordinate of this tissue sample on the voxel‐wise FCD‐neuroticism correlation map. The t‐values of voxels within the sphere were averaged to represent the t‐value of this tissue sample, by which a sample‐wise FCD‐neuroticism correlation matrix was obtained for this donor. By repeating these steps, we can obtain sample‐wise FCD‐neuroticism correlation matrices for other five donors. In this case, the sample‐wise gene expression matrix (n × p1) comprised the predictor variables, and the sample‐wise FCD‐neuroticism correlation T map vector comprised the response variables (n × p2), where p1 indicates the number of high SES genes (p1 = 1,751 when the SES cutoff = 0.6) used in main analysis, p2 indicates the number of imaging measures (here, the neuroticism‐related FCD T map) used in the main analysis, and n indicates the number of brain samples with gene expression data used in our main analysis. Specifically, in the PLS regression of Donor 1, the size of the predictor matrix was (246 × 1,751) and the size of the response vector was (246 × 1). In the PLS regression of Donor 2, the size of the predictor matrix was (328 × 1,751) and the size of the response vector was (328 × 1). In the PLS regression of Donor 3, the size of the predictor matrix was (296 × 1,751) and the size of the response vector was (296 × 1). In the PLS regression of Donor 4, the size of the predictor matrix was (325 × 1,751) and the size of the response vector was (325 × 1). In the PLS regression of Donor 5, the size of the predictor matrix was (641 × 1,751) and size of the response vector was (641 × 1). In the PLS regression of Donor 6, the size of the predictor matrix was (549 × 1,751) and the size of the response vector was (549 × 1).

Second, the PLS identified components from the gene expression matrix that had maximum covariance with the response variables in the FCD‐neuroticism correlation T map vector. The PLS components were ranked by the covariance between the gene expression matrix and FCD‐neuroticism correlation T map vector, so the first several PLS components will transform the covariance between the high‐dimensional data matrices into the best low‐dimensional representation. In our case, we first conducted individual PLS regressions on the high‐SES gene expression matrix in each of the six donated brain expression datasets to find linear combinations of genes with significant predictive ability for the response variables (six sample‐wise FCD‐neuroticism correlation T‐maps).

Third, we subsequently used the bootstrapping resampling method (resampling with replacement of the sample label in each donor) to test the goodness of fit of the estimated PLS component by repeating the analysis 1,000 times after shuffling the sample labels assigned to the FCD‐neuroticism correlation T‐map. The error in the PLS weights for each gene was also estimated using bootstrapping and then we ranked the genes according to their contributions to each PLS component calculated as the ratio of the weight of each gene to its bootstrapped SE. Furthermore, we used significant genes with both negative and positive weights in the following functional annotation analysis as suggested by (Romero‐Garcia, Warrier, Bullmore, Baron‐Cohen, & Bethlehem, 2018). Overlapping genes among the six PLS regressions that significantly contributed to the corresponding PLS component for all six donors were used for subsequent functional annotation analysis.

Functional annotations

With the Enrichr online toolset (http://amp.pharm.mssm.edu/Enrichr/), the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway was used to identify biological pathways of surviving genes in the significant PLS components; and gene ontology (GO) gene sets were used to find biological processes in the GO of these genes. Cell‐type specific expression analysis (CSEA) (http://genetics.wustl.edu/jdlab/csea-tool-2/) was used to disclose the cell‐type enrichment of these genes under a particular specificity index threshold (Xu, Wells, O'Brien, Nehorai, & Dougherty, 2014). We then performed enrichment analysis with Fisher's exact test for surviving genes in the significant PLS components using the 275 neuroticism‐related GWAS significant genes (Nagel, Watanabe, Stringer, Posthuma, & van der Sluis, 2018) and the 8 eQTL significant genes in brain tissue from GTEx database (Nagel, Jansen, et al., 2018). Protein–protein interaction (PPI) hub protein analysis was performed by STRING v11.0 (https://string-db.org/) to identify the hub genes in the significant PLS components. Hub genes were defined by the top half of the node degree in the PPI network at a medium confidence interaction score and consistently showed reliability for different SES thresholds. Finally, spatiotemporal gene expression analysis was used to characterize the gene expression trajectory of hub genes by using the Human Brain Transcriptome database (http://hbatlas.org/).

A schematic diagram of the analysis step is provided in Figure 1.

FIGURE 1

Schematic diagram of the methodology. (a) The SES of one certain gene is defined as the average Pearson correlation coefficient over 15 pairs of 6 donated brains (b) Sample‐wise FCD‐neuroticism correlation matrices are obtained via extracting the average T value of one voxel‐wise FCD‐neuroticism correlated T map within each spherical region of interest (ROI, radius = 6 mm) in each of six donated brains respectively. Different colored spheres indicate samples in six donated brains. (c) Six individual PLS regressions are conducted on high SES gene expression matrix in each of six donated brain expression datasets to find the linear combinations of genes with the significant predictive ability for the response variables (six sample‐wise FCD‐neuroticism correlation T‐maps). (d) Overlapped significant genes among six PLS regressions that contributed to one certain PLS component would be used for following (e) functional annotation analysis (GO, KEGG, cell‐type specific enrichment analysis, PPI and spatiotemporal expression analysis). (f) The reliability and specificity analysis includes reliability of FCD‐neuroticism correlations, validation of SES calculation, validation of PLS regression and gene enrichment analysis, specificity analysis of neuroticism‐imaging correlation analysis, specificity analysis of neuroticism‐related transcription‐neuroimaging associations. FCD, functional connectivity density; GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; PLS, partial least squares; PPI, Protein–protein interaction; SES, spatial expression stability

Statistical analyses

Demographic and behavioral data were analyzed using the Statistical Package for the Social Sciences (SPSS 20.0, Chicago, IL). Descriptive data are shown in the mean ± SD and data range.

For all voxel‐wise statistical analyses, multiple comparisons were corrected using family‐wise error (FWE) cluster‐level correction with a corrected threshold at p < .05 and an initial voxel‐level p < .001.

The significance of the explained variance of each PLS component was assessed by the permutation test (p < .05). To correct winners curse bias, a false discovery rate (FDR) inverse quantile transformation correction was used to obtain genes that survived after FDR correction of p < .05 (Bigdeli et al., 2016). All these enrichment analyses were considered significant after Benjamini–Hochberg FDR correction (p < .05).

RESULTS

Initial findings

Personality assessment

Personality assessment is shown as follows as the mean ± SD and range: neuroticism raw score (10.81 ± 5.38, range: 1–24); neuroticism T score (50.04 ± 10.03, range: 31.77–74.52); trait anxiety (36.96 ± 8.45, range: 20–75); harm avoidance (14.79 ± 6.42, range: 1–32); and depression (7.79 ± 6.87, range: 0–52).

Neuroimaging–neuroticism correlations

In 274 healthy young adults, we used the GMV as a structural imaging measure and the FCD as a functional imaging measure to investigate the voxel‐wise correlations between neuroimaging measures and neuroticism scores (p < .05, FWE cluster level corrected, voxel p < .001). Although there was no significant correlation between neuroticism and the GMV, we found a significant positive correlation between neuroticism and the FCD in the ventral striatum (peak MNI coordinate: x = −6, y = 6, z = −6, peak T (1,265) = 4.53, cluster size = 130 voxels, Figure 2a). The connection probability map showed that this significant cluster was functionally connected with brain regions in the anterior cingulate cortex (ACC)–striatum–thalamus circuit (Figure 2b).

FIGURE 2

Correlations between neuroticism and FCD and connectivity pattern of significant clusters. (a) Brain region with a correlation between neuroticism and FCD at connection threshold of r = .6. (b) Functional connection probability map of significant cluster (ventral striatum) at connection threshold of r = .6 in 274 participants. FCD, functional connectivity density

Genome‐wide SES

We used SES to measure the reliability of the expression patterns of the 20,737 genes across brain structures in the six donated brains. The high SES genes represent those with highly consistent expression patterns across the brains. Based on these genes, we identified canonical genes that were associated with the neuroticism‐related FCD statistical map (T‐map). Using a dichotomous cutoff value of SES = 0.6, 1,751 genes were obtained from 20,737 genes with higher SES values (Figure 2a), which were then used in the PLS regression analysis.

Gene expression associated with neuroticism‐related FCD map

With gene expression data from AHBA brains and the neuroticism‐FCD correlation T‐maps from the 274 healthy young adults, we explored the sample‐wise spatial correlations between gene expression and the neuroticism‐related FCD T‐maps. Here, PLS regression analysis was performed between the gene expression of each donated brain and the neuroticism‐related FCD T‐map separately to discover consistent findings across the six AHBA brains. The top two PLS components significantly explained the variances (PLS1: 27.84–34.29%, PLS2: 5.75–11.10%) of the neuroticism‐related FCD T‐map after a permutation test (p < .05). Although the PLS1 scores calculated from each donated brain were positively correlated with the neuroticism‐related FCD T‐map (Figure S1), the observed associations could be mainly explained by the gene expression differences between cortical and subcortical samples (Hawrylycz et al., 2012). Therefore, the PLS2 scores calculated from each donated brain may truly represent the positive association with the neuroticism‐related FCD T‐map (r = .28–.42, Figure 3b). Specifically, brain regions with higher PLS2 scores also have dense connections between neuroticism and FCD (Figure 3b). Using FDR correction (p < .05), 161 overlapping genes that significantly contributed to the PLS2 component for all six donors were used for subsequent functional annotation analysis. Each PLS2 gene weight and its bootstrapped standard error for each donor at a threshold SES = 0.6 are shown in Supplementary file S1.

FIGURE 3

Relationship between significant genes in PLS2 and neuroticism‐related FCD map. (a) Gene ranks according to SES and number of genes at different SES threshold. Then, 1,751 genes at threshold SES = 0.6, 3,236 genes at threshold SES = 0.5, 740 genes at threshold = 0.7. (b) Positive correlations between PLS2 scores and T values of neuroticism‐related FCD map in each donor at threshold SES = 0.6. (c) Positive correlations between PLS2 scores and T values of neuroticism‐related FCD map in each donor at threshold SES = 0.5. (d) Positive correlations between PLS2 scores and T values of neuroticism‐related FCD map in each donor at threshold SES = 0.7. Pearson correlation coefficients, 95% confidence interval (in parentheses) and p values are provided in the bottom of the respective panels. Amy, amygdala; BF, basal forebrain; CgG, cingulate gyrus; Cl, claustrum; dStriatum, dorsal striatum; dThalamus, dorsal thalamus; Epi, epithalamus; FC, frontal cortex; FCD, functional connectivity density; GP, globus pallidus; Hippo, hippocampus; Hypo, hypothalamus; Ins, insula; OC, occipital cortex; PC, parietal cortex; PG, parahippocampal gyrus; PiC, piriform cortex; PLS, partial least squares; SES, spatial expression stability; sThalamus, subthalamus; TC, temporal cortex; vStriatum, ventral striatum; vThalamus, ventral thalamus

Functional annotations

The PLS2 genes were enriched in 17 pathways in the KEGG database in the initial findings. However, only seven pathways were consistently significant in validation analysis (p < .05, FDR correction, Figure 4). Enriched pathways included serotonergic synapse, cholinergic synapse, dopaminergic synapse, circadian entrainment, long‐term potentiation, inflammatory mediator regulation of transient receptor potential (TRP) channels, and amphetamine addiction. In the GO database of biological processes, the PLS2 genes showed significant enrichment in the biological process of chemical synaptic transmission. The cell‐type specific analysis showed that the PLS2 genes exhibited a significant enrichment in neurons (p < .05, FDR correction, Figure 4), but not in any other cell types including astrocytes, oligodendrocytes and immune cells (p > .01, FDR correction). We identified the PLS2 genes that were enriched for neuroticism‐related GWAS significant genes (p = .025, FDR correction) and neuroticism related eQTL significant genes in the brain tissue from the GTEx database (p = .046, FDR correction). To further identify the hub genes in the enriched pathways and biological processes, we performed PPI analysis and found that PLS2 gene‐coded proteins formed a PPI network (Figure 5) with four hub genes, including G protein subunit gamma 10 (GNG10), 5‐hydroxytryptamine receptor 2C (HTR2C), prodynorphin (PDYN), and calcium/calmodulin‐dependent protein kinase II alpha (CAMK2A). These four genes were involved in the enriched pathways and biological processes and consistently ranked in the top half in terms of node degree. In spatiotemporal gene expression analysis, HTR2C and PDYN were highly expressed in the striatum, hippocampus, and amygdala, and the temporal expression curves of CAMK2A and HTR2C exhibited a gradual increase from the embryo stage and reached a plateau between the child and young adult stage. However, GNG10 is not available in the Human Brain Transcriptome database.

FIGURE 4

Enrichment analysis for PLS2 genes. KEGG pathway, GO, GWAS, eQTL, cell types enrichment analyses for significant genes in PLS2 associated with neuroticism‐related FCD map. The same color of a bar graph in each panel indicates the same enriched class. BP, biological process; eQTL, expression quantitative trait loci; FCD, functional connectivity density; GO, gene ontology; GWAS, genome‐wide association analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes; PLS, partial least squares; TRP, transient receptor potential

FIGURE 5

Hub genes in the PLS2 component. Protein–protein interaction network visualization (a) and spatiotemporal expression curves of the HTR2C, CAMK2A, and PDYN (b). The description of Period 1–15 in spatiotemporal expression curves was provided in http://hbatlas.org/. AMY, amygdala; CBC, cerebellar cortex; HET, heterozygotes; HIP, hippocampus; MD, mediodorsal nucleus of the thalamus; NCX, neocortex; STR, striatum

Reliability and specificity analyses

Reliability of FCD‐neuroticism correlations

To ensure that the results of the FCD‐neuroticism correlations derived from the connection threshold of 0.6 were reliable, we repeated the voxel‐wise correlation analysis using connection thresholds of 0.5 and 0.7. Similar to the main findings, we found that neuroticism was only positively correlated with FCD in the ventral striatum under connectivity thresholds of 0.5 (peak MNI coordinate: x = −6, y = 6, z = −6, peak T (1,265) = 4.54, cluster size = 107 voxels, Figure S2a) and 0.7 (peak MNI coordinate: x = −9, y = 6, z = −6, peak T (1,265) = 3.94, cluster size = 107 voxels, Figure S2b). We found a highly significant spatial correlation (r = .91, p < .001) between the uncorrected voxel‐wise FCD‐neuroticism correlation map with and without controlling for trait anxiety and harm avoidance (Figure S3), indicating that removing the variance from trait anxiety and harm avoidance may not obviously bias our results. With HCP Caucasian data, we found a significant positive correlation between neuroticism and FCD in the ventral striatum after FWE correction (p < .001, peak MNI coordinate: x = −12, y = 10, z = −22, peak T (1, 55) = 5.18, cluster size = 258 voxels). The statistical T map of the Chinese data was overlaid with the T map of HCP data with initial voxel‐level p < .001 (Figure S4).

Validation of SES calculation

The SES calculated after leaving each donor out showed high correlation with the SES calculated for all six donors (r leave‐donor1‐out = .9805, 95% CI (0.9800–0.9811), p < .001; r leave‐donor2‐out = .9760, 95% CI (0.9754–0.9767), p < .001; r leave‐donor3‐out = .9837, 95% CI (0.9833–0.9842), p < .001; r leave‐donor4‐out = .9861, 95% CI (0.9857–0.9865), p < .001; r leave‐donor5‐out = .9814, 95% CI (0.9809–0.9819), p < .001; r leave‐donor6‐out = .9866, 95% CI (0.9862–0.9869), p < .001). We also found PLS2 could significantly explain the variance (5.29–11.07%) of the neuroticism‐related FCD T‐map after the permutation test (p < .05) (Figure S5). Functional annotation analysis revealed similarly consistent enriched pathways, in the PPI network (Figure S6).

The performance of average Kendall tau was closely related to average the Pearson correlation (r = .9312, 95% CI [0.9293–0.9329], p < .001). We also provided the CI of the SES value of each gene with a high SES at the 0.6 threshold in Supplementary file S2.

Validation of PLS regression and gene enrichment analysis

We used two different SES cutoff values (0.5 and 0.7) to perform genome‐wide SES analysis and obtained consistent expressed genes across the six donors (3,236 genes for SES = 0.5 and 740 genes for SES = 0.7). Then, these genes were also used in the subsequent PLS regression separately. The procedures of the PLS regression using different SES thresholds were the same as our initial analysis in the calculation of variance explained by the PLS components and correlation analysis between PLS scores and the neuroticism‐related neuroimaging map. We also found that only PLS2 truly displayed a positive association with the neuroticism‐related FCD T‐map (Figure 3c,d). In the gene enrichment analysis, PLS2 genes also had a significant enrichment for similar GO terms, pathways, cell‐types, GWAS significant genes, eQTL significant genes in the brain tissue of the GTEx database (except for PLS2 genes at an SES threshold 0.5, p = .071, FDR correction), and PPI analysis with initial enrichment analysis (Figures 4 and 5).

Specificity analysis of neuroticism‐imaging correlation analysis

We performed similar voxel‐wise correlation analysis between FCD and extroversion, which is one of the two independent personality traits initially proposed by Eysenck (Bech, Lunde, & Moller, 2012). Using the same statistical threshold of p < .05, FWE cluster level corrected, voxel p < .001 as in FCD‐neuroticism correlations, we did not find any significant correlation between FCD and extroversion. To compare the distribution differences in brain regions whose FCD values correlated with neuroticism and extroversion, we used a lenient threshold (voxel‐wise p < .01, uncorrected, cluster size = 10) to depict the correlation maps between FCD and neuroticism and extroversion. We found completely different distributions between brain regions whose FCD values were correlated with extroversion and those correlated with neuroticism (Figure S7).

Specificity analysis of neuroticism‐related transcription‐neuroimaging associations

The PLS2 identified in the transcription‐neuroimaging correlations for neuroticism was used to investigate whether there was also a significant correlation between PLS2 scores and the extroversion‐FCD correlation T‐map. However, we only found much weaker correlations in three donors using the one‐tailed Fisher‐z test (r = −.12 to .11; p donor1 < 1 × 10−4, Z donor1 = 6.26; p donor2 = 6 × 10−4, Z donor2 = 3.23; p donor3 = 3 × 10−4, Z donor3 = 3.41; p donor4 = 5 × 10−4, Z donor4 = 3.27; p donor5 < 1 × 10−4, Z donor5 = 3.95; p donor6 = 2 × 10−4, Z donor6 = 3.61) (Figure S8).

DISCUSSION

To our knowledge, this is the first study that combines cross‐scale correlation analyses and up‐to‐date functional annotations to investigate the neurobiological correlates of human neuroticism. We found that the functional connectivity of the ventral striatum was most strongly associated with human neuroticism and the neuroticism‐related functional connectivity was related to the expression of the GNG10, HTR2C, PDYN, and CAMK2A hub genes of the chemical synaptic transmission, circadian entrainment, long‐term potentiation, inflammatory mediator regulation of TRP channels, and amphetamine addiction pathways in cortical and striatal neurons.

There was no significant neuroticism‐GMV correlation, although one previous study showed a negative association between the GMV of the left superior frontal gyrus and neuroticism (Lu et al., 2014). One of the probable reasons is that the striatum is closely related to human neuroticism and striatal functional connectivity could better account for the individual differences in neuroticism than striatal volume. Another important reason is that the multiple comparison correction method used in a previous study is insufficient for reducing false‐positive rates (Silver, Montana, Nichols,, & Alzheimer's Disease Neuroimaging, 2011). However, we found that individuals with high neuroticism scores had high FCD values in the ventral striatum, which was reliable across different connection thresholds. In the original Eysenck's personality theory (Bech et al., 2012), extraversion and neuroticism were identified as two possible personality traits with distinct neural and genetic substrates. To prove neuroimaging findings specific to neuroticism, we performed a voxel‐wise extraversion‐FCD correlation analysis and found that brain regions with extraversion‐FCD correlations were completely different from those with neuroticism‐FCD correlations, which is in agreement with previous findings of the differences between the neuroticism‐ and extraversion‐related functional connectivity of the amygdala and brain network topology (Aghajani et al., 2014). Our findings suggest that the functional connectivity of the ventral striatum is specifically associated with human neuroticism, which is also consistent with prior findings reporting that individuals with high neuroticism scores demonstrate enhanced activity in the ventral striatum in response to reward stimuli (Schaefer, Knuth, & Rumpel, 2011), and that neuroticism is positively correlated with dopamine D2 receptor density in the striatum (Lee et al., 2005). Since neuroticism has been considered as a risk factor for depressive and anxiety disorders (Griffith et al., 2010), its neural substrate, ventral striatum functional connectivity may be associated with these disorders, which is supported by the predictive value of this connectivity for future risk in adolescents for depressive disorder (Pan et al., 2017). Functional connectivity probability mapping revealed that the ventral striatum was functionally connected with brain regions of the ACC–striatum–thalamus circuit, which is responsible for mediating motivated behavior and lesions of this circuit manifest as apathy (Mayer & Richard, 2001). More importantly, abnormal connectivity in this circuit can be used to distinguish the neurophysiological subtypes of depression and deep stimulation of this circuit can relieve the core symptoms of depression and anxiety disorders (Drysdale et al., 2017; Figee et al., 2013). These findings indicate that the functional connectivity of the ventral striatum may be a neural substrate bridging neuroticism and mood disorders, which provides a neurobiological explanation for the increased risk of neuroticism in depressive and anxiety disorders (Griffith et al., 2010).

Using PLS regression, we identified a genetic component (PLS2) for each donor that could significantly explain the variance of the neuroticism‐related FCD map, which showed reliability in gene sets extracted with different SES cutoff values. To ensure that the obtained genetic component was specific to the neuroticism‐related FCD map, we also performed sample‐wise spatial correlation analysis between the transcription profiles of PLS2 and the extraversion‐related FCD map and found weaker correlations (r = −.12 to .11) than those obtained with the neuroticism‐related FCD map (r = .28–.42). These findings suggest that the PLS2 genes are specific to the neuroticism‐related FCD map.

The PLS2 genes were consistently enriched in the serotonergic synapse, cholinergic synapse, dopaminergic synapse, circadian entrainment, long‐term potentiation, inflammatory mediator regulation of TRP channels, and amphetamine addiction of KEGG pathway and chemical synaptic transmission of biological process ontology and in cortical and striatal neurons in cell‐type specificity analysis, indicating that these genes may affect neuroticism‐related functional connectivity by regulating trans‐synaptic neurotransmitter transmission, including serotonin, choline and dopamine. In addition, it may also work on circadian entrainment, long‐term potentiation, TRP channels, and the amphetamine addiction pathway to influence neuroticism‐related FCD. Interestingly, the serotonergic, dopaminergic synapse, and circadian entrainment pathways have been previously implicated in neuroticism‐related GWAS and molecular imaging studies (Farde, Plaven‐Sigray, Borg, & Cervenka, 2018; Ferguson et al., 2018; Nagel, Jansen, et al., 2018). These results were further supported by enrichment results obtained using neuroticism‐related GWAS and eQTL genes. Although the pathways of cholinergic synapse, long‐term potentiation, TRP channels, and amphetamine addiction have not previously been directly associated with neuroticism, many studies have shown that these pathways are involved in several neuropsychiatric disorders, such as anxiety, depression and addiction (Conrad & Winder, 2011; Naziroglu & Demirdas, 2015; Salo et al., 2011; Scarr, Gibbons, Neo, Udawela, & Dean, 2013), suggesting that our study may advance the neurobiological knowledge of neuroticism. In the PPI analysis of the PLS2 genes associated with neuroticism‐related FCD, we consistently observed significant enrichment of four PPI hub proteins (GNG10, HTR2C, PDYN, and CAMK2A).

The 5‐HT2C receptor encoded by HTR2C is a subtype of 5‐HT receptor that binds the endogenous neurotransmitter serotonin and regulates dopamine release in the striatum, prefrontal cortex, hippocampus, hypothalamus, and amygdala (Alex, Yavanian, McFarlane, Pluto, & Pehek, 2005). In addition, 5‐HT2C receptors have been associated with several behaviors and psychological states, including mood, anxiety, and depressive states (Heisler, Zhou, Bajwa, Hsu, & Tecott, 2007; Millan, 2005). PDYN is an opioid polypeptide hormone encoded by PDYN, which is involved in chemical signal transduction and is highly enriched in the striatum. Moreover, individuals with mood disorders exhibit a reduction in the mRNA expression of PDYN in the amygdala (Hurd, 2002). The protein encoded by CAMK2A is a member of the Ca2+/calmodulin‐dependent protein kinase II subfamily and plays a critical role in hippocampal long‐term potentiation (Lisman, Yasuda, & Raghavachari, 2012) and amphetamine action (Steinkellner et al., 2015). In addition, the interaction of CAMK2A function with serotonin receptor 1A during the postnatal period results in the development of anxiety‐related behavior (Lo Iacono & Gross, 2008). Another hub gene GNG10 encodes a member of the gamma‐subunit of a modulator of the G‐protein‐coupled receptor signaling pathways, which has been little investigated except for melanoma (Cardenas‐Navia et al., 2010).

The AHBA provided the only high‐resolution gene expression dataset derived from brain samples; thus, we were able to investigate the association between macroscale MRI data and microscale gene expression. However, there are two aspects that should be further investigated. First, whole genome‐wide genotyped data should be collected in future studies, and we could use polygenic risk scores for significant genes derived from AHBA gene expression analysis to study the relationship with behavior and imaging metrics. Fortunately, the establishment of world's largest sample of the Chinese Han imaging genetics cohort, the Chinese Imaging Genetics will provide opportunity to identify these associations (Xu et al., 2020). Second, bioinformatic analysis should be expanded to biological validation, such as behavioral analysis of knockout mice, which could further identify causal relationships rather than only putative associations.

However, several limitations should be noted. First, we employed an eye‐closed rfMRI paradigm, which may have significantly affected functional connectivity. Previous studies indicated that the eye‐open paradigm demonstrates significantly greater resting‐state activity in visual and attention networks and significantly lower activity in the sensorimotor network (Wei et al., 2018). Second, the gene expression dataset and imaging dataset used in our study came from different individuals. To reduce the impact of this limitation, we first found robust neuroticism‐related imaging associations in our fMRI dataset and investigated transcription‐imaging associations in each individual gene expression dataset using consistently expressed genes.

CONCLUSION

In this study, we found that human neuroticism was correlated with functional connectivity in the ventral striatum and the neuroticism‐FCD correlations could be specifically explained by a major genetic component implicated in the chemical synaptic transmission, circadian entrainment, long‐term potentiation, inflammatory mediator regulation of TRP channels, and amphetamine addiction pathways in neurons. In the PPI analysis, we identified several genes (GNG10, HTR2C, PDYN, and CAMK2A) as hub genes associated with the neuroticism‐related FCD map. These findings may improve our understanding of the genetic and neural substrates of human neuroticism.

CONFLICT OF INTEREST

The authors declare no conflicts of interest.

AUTHOR CONTRIBUTIONS

C.Y. and T.J. designed the experiment. Q.X. and F.L. performed the experiments. Q.X., W.Q., and F.L. analyzed the data. Q.X. and C.Y. drafted and revised the manuscript. All authors discussed the results.

Appendix S1: Supporting Information

Click here for additional data file.

Supplementary file S1 Each PLS2 gene weight and its bootstrapped standard error in each donor at threshold SES = 0.6

Click here for additional data file.

Supplementary file S2 The SES and its bootstrapped CI of high SES genes at threshold 0.6.

Click here for additional data file.