Participants
The dataset comprises individuals from the Behavioral Brain Research Project of Chinese Personality (BBP), which launched in September 2019 (still in progress) to recruit participants from that year’s freshmen at Southwest University, Chongqing, China. All procedures involving human subjects were approved by the Ethics Committee of the Southwest University. Written informed consents were obtained from all subjects. In total, 1059 participants completed the cross-sectional neuroimaging protocol, the schizotypal personality questionnaire-Brief (SPQ-B), as well as basic demographical variables, i.e., headedness, age, gender. Psychiatric history of the participants was collected via self-report. Forty-six participants were excluded due to excessive head motion during scanning (with a mean framewise displacement (FD) larger than 0.25 mm), resulting in a final sample of 1013 participants (mean age = 19.11 years old, SD = 0.95; 336 males and 677 females; 90.33% right handedness; mean FD = 0.09, SD = 0.04), see Table S1 for details.
Questionnaire Assessment
The schizotypal personality questionnaire-Brief (SPQ-B) is a 22-item self-report scale that measures schizotypal traits (36). These items can be further grouped into three factors, i.e., cognition-perceptual, interpersonal and disorganization factors. SPQ scale modeled on the DSM-III-R criteria for schizotypal personality disorder. SPQ-B has been widely used to screen for psychometrically defined schizotypal personality disorder and identify people who may be at an elevated risk for prodromal mental illness in the general population (37). The Chinese version of the SPQ-B used in the present study has also been shown to possess good reliability (test-retest reliability = 0.82; (38)). Consistent with previous study in the literature (34, 39), we used the SPQ total score to characterize the schizotypal trait. In the present study, the internal consistency for the whole scale is evidenced by a Cronbach α coefficient of 0.85.
The measurement of family socioeconomic status (SES) consists of two questions answered by the participants: annual family income (1–7 level, Table S1) and highest educational qualification of each parent. We converted the highest educational qualification into years taken to complete the degree. The averaged educational years of both parents was used. Consistent with previous study (40), we used the averaged z-scores for family annual income and averaged parental education length to index SES.
MRI Acquisition And Image Data Preprocessing
All structural and resting-state functional MRI data were acquired on a 3T Siemens Primsa-fit scanner with a standard 32 channel head coil located at the MRI Center of Southwest University. All participants were instructed to stay awake, keep their eyes open, fixate on the displayed crosshair, and move their head as little as possible. See supplementary materials for detailed scanning parameters.
All preprocessing steps were consistent with our previous study(41), See supplementary materials for details. To account for artifact of head motion, we excluded participants with > 0.25 mean framewise displacement (FD) and regarded mean FD as a covariate in statistical analysis.
Connectivity Gradient Mapping
We used the diffusion map embedding approach (42), implemented in BrainSpace Toolbox (43) to construct the connectome gradient for each participant. First, we constructed the parcellation-based functional connectivity matrix using Fisher’s Z-transformed Pearson’s correlations. We used Schaefer's 1000 cortical parcellation (44) and Brainnetome 36 subcortical parcellation (45), resulting in 1036 subregions. We excluded the cerebellum because not all participants’ cerebellum was fully covered during fMRI scanning. For each subject, we obtained a FC matrix (1036*1036), then retained the top 10% of connections per row retained, and further converted the thresholded connectivity matrix into cosine distance-based similarity matrix (24, 25). Finally, diffusion map embedding approach, a nonlinear and unsupervised dimensionality reduction algorithm, was adopted to capture spatial axes (gradient components) in connectivity variation across different areas.
Schizotypy-related Connectome Hierarchy
In order to implement statistical analysis across participants, we performed iterative Procrustes rotation to align the gradients across participants (46). For each subject, four global topographic metrics were used to quantify gradient properties for the two connectome gradients 1 and 2: (1) explanation ratio, which characterizes the percentage of connectivity variance explained by a given gradient (20); (2) gradient range, which characterizes the difference between the greatest positive and negative value in a given gradient (24); (3) gradient variation, which represents the standard deviation of the gradient scores (47); (4) gradient dispersion, which characterizes the summed square of the Euclidean distance from each brain nodal to the centroid of the gradient space determined by the first two gradients (47). Considering the SPQ total score did not follow a normal distirbution (Fig.S1, Kolmogorov-Smirnov test, p < 0.01), the association between connectome gradient metrics (global topography and regional gradient value) and schizotyoy were assessed by using a permutation analysis of general linear model (48) with age, gender, handedness and head motion as covariates. Consistent with previous studies (24, 25), Z-normalization was conduct to gradient map before the regional level association analysis. For global gradient metrics, the threshold of statistical significance was set to p < 0.05. For regional gradient score maps, the threshold of statistical significance was set to FDR q < 0.05.
Meta-analytic Decoding Of The Function Of Significant Connectome Hierarchy Using NeuroSynth
We used NeuroSynth (49), a large-scale database-informed meta-analytic approach to decode the functional properties associated with hierarchy alterations in schizotypy. The top 20 terms (excluding anatomical terms for brain regions) exhibiting the highest correlations with schizotypy-related positive or negative hierarchy associative regions were extracted. These cognitive terms were visualized on a word cloud. The font size of these terms in each word cloud is proportional to their correlation strength with corresponding meta-analytic neuro-maps included in Neurosynth. Considering its nature of exploratory and inverse inference, it should be noted that, the corresponding functional properties should be interpreted cautiously.
Association Between Schizotypy-related And Schizophrenia-altered Connectome Hierarchy
To further assess the similarity between schizotypy-related and schizophrenia-altered connectome hierarchy, we calculated the spatial correlation between schizotypy-related connectome hierarchy and our previous case-control study of schizophrenia-altered cortical connectome hierarchy (Schizophrenia patients, n = 96; Health control, n = 122; (25)). To keep the same dimension, we extracted the regional average T-statistical values for Schaefer's 1000 cortical parcellations from our previous study and computed the Pearson correlation. The significance levels were determined by using permutation tests based on spherical rotations (50), to account for spatial autocorrelation of schizotypy-related connectome hierarchy map (5,000 times).
Estimation Of Regional Gene Expression
First, we downloaded the Allen Human Brain Atlas (AHBA, http://human.brain-map.org), a whole-genome expression atlas of the human brain created by the Allen Institute for Brain Sciences using 6 donors 24 to 57 years of age (51). Then, we followed Allen Human Brain Atlas (AHBA) processing pipeline with the default recommended settings to obtain brain regional gene expression. Briefly, preprocessing for the gene expression microarray data included the following steps: verification of the probe-to-gene annotations, probe selection, tissue sample assignment to each donor’s cortical parcellation with 1000 regions, normalization across donors, gene-set filtering. After that, for each region, gene expression was obtained by calculating mean of all cohorts in a region, resulting in 1,000 × 10,027 regional transcription matrices. See Arnatkevic̆iūtė, Fulcher (52) for more details. Considering AHBA dataset only includes right hemisphere for two donors, in the subsequent analysis, we only used gene expression in left hemisphere, i.e., the region × gene matrix size was 500 × 10,027.
Association between schizotypy-related connectome hierarchy and gene expression profiles
We used Partial Least Squares (PLS) regression to associate the regional schizotypy-related connectome hierarchy to the gene expression measurements for all 10, 027 genes (53, 54). In the PLS regression model, gene expression date was regarded as the predictor variables to predict the regional schizotypy-related connectome hierarchy (response variables). The first PLS component (PLS1) is the linear combination of the predictor variables (i.e., gene expression) that can explain most of the variance in the response variables (i.e., regional schizotypy-related connectome hierarchy). In order to evaluate the significance of the PLS1, we adopted permutation testing based on spherical rotations, to account for spatial autocorrelation of schizotypy-related connectome hierarchy map (5,000 times) to determine whether the variance explained by PLS1 was significantly greater than that expected by chance (50). The variability (standard error) of each gene’s weight on PLS1 was assessed by bootstrapping, and then we divided the weight of each gene by its bootstrap-estimated standard error to obtain the corrected weight for each gene.
We ranked the genes according to their corrected weights, which represent their contribution to the PLS1. Both the ascending and descending gene lists were enrolled for the following gene enrichment analysis. We used GOrial (http://cbl-gorilla.cs.technion.ac.il/, Eden, Navon (55)) to identify the enriched Gene Ontology terms for each ranked list of genes. Consistent with previous studies, we regarded terms with Benjamini-Hochberg false discovery rate (FDR)-corrected q-value less than 0.05 as significant enrichment.
In addition, we directly tested whether there was a significant correlation between PLS1 weights of all genes and the differential expression values in postmortem brain tissue measurements of messenger RNA from case-control studies of schizophrenia. To complement this, we used 2 lists of genes that were transcriptionally dysregulated: (1) 1) the list of genes reported by Gandal, Haney (56) in the prefrontal and parietal brain regions in schizophrenia (n = 159); and 2) the list of genes reported by Fromer, Roussos (57) in the dorsolateral prefrontal cortex in schizophrenia (n = 258). Correlations were calculated using genes that were common across postmortem and PLS1 gene lists. We used repeated random relabeling of genes to test the significance of correlation. To test whether the PLS1 weights gene captured the specificity of schizophrenia, we also used the same procedures to data of differential gene expression from case-control studies of bipolar disorder (BD), autism spectrum disorder (ASD) and major depressive disorder (MDD), which was also reported by Gandal, Haney (56).
Control Analysis
We validated schizotypy-related connectome hierarchy by considering several control analyses. First, considering global signal regression (GSR) is still controversial in the rsfMRI image preprocessing, we repeated core analyses with GSR in the control analyses. Second, we used a different parcellation scheme to validate our results. Specifically, 1078 regions were used to construct connectivity matrix, including 1024 cerebral regions with uniform sizes (58) and 54 subcortical regions (59).
Data And Code Availability
The statistical maps of the associations between schizotypy and connectome hierarchy, and gene list revealed by PLS will be provided once published at GitHub (https://github.com/DeboDong/Hierarchy_schizotypy). The code for spatial permutation testing can be found at https://github.com/frantisekvasa/rotate_parcellation. Connectome gradients were constructed by using BrainSpace toolbox (https://github.com/MICA-MNI/BrainSpace). we followed Allen Human Brain Atlas (AHBA) processing pipeline (https://github.com/BMHLab/AHBAprocessing) with the default recommended settings to obtain brain regional gene expression. The PLS code to associate schizotypy-related connectome hierarchy with gene expression can be found at https://github.com/SarahMorgan/Morphometric_Similarity_SZ. We used GOrial (http://cbl-gorilla.cs.technion.ac.il/) to implement the Gene Ontology enrichment analysis and visualization. Brain maps were visualized using BrainNet Viewer (60), https://www.nitrc.org/projects/bnv/.