https://orcid.org/0000-0002-6342-2951
Abstract
The human visual cortex is a heterogeneous entity that has multiple subregions showing substantial variability in their functions and connections. We aimed to identify genes associated with resting-state functional connectivity (rsFC) of visual subregions using transcriptome-neuroimaging spatial correlations in discovery and validation datasets. Results showed that rsFC of eight visual subregions were associated with expression measures of eight gene sets, which were specifically expressed in brain tissue and showed the strongest correlations with visual behavioral processes. Moreover, there was a significant divergence in these gene sets and their functional features between medial and lateral visual subregions. Relative to those associated with lateral subregions, more genes associated with medial subregions were found to be enriched for neuropsychiatric diseases and more diverse biological functions and pathways, and to be specifically expressed in multiple types of neurons and immune cells and during the middle and late stages of cortical development. In addition to shared behavioral processes, lateral subregion associated genes were uniquely correlated with high-order cognition. These findings of commonalities and differences in the identified rsFC-related genes and their functional features across visual subregions may improve our understanding of the functional heterogeneity of the visual cortex from the perspective of underlying genetic architecture.
Introduction
The human visual cortex, primarily located in the occipital lobe, is a heterogeneous entity that has multiple subregions showing substantial variability in their functions. According to classical Brodmann’s map of the human cortex, the occipital lobe could be parcellated into three subregions including the primary visual cortex (BA17) and visual association cortex (BA18 and BA19) (Zilles and Amunts 2010). It is well established that the primary visual cortex predominantly involves receiving visual information from the lateral geniculate nucleus and then delivering information to other visual regions for further complex processing after initial integration (Prasad and Galetta 2011). There is also growing evidence that the primary visual cortex is engaged by higher level cognitive processes (Roelfsema and de Lange 2016) via interactions with widespread cortical regions through numerous feedback projections (Clavagnier et al. 2004; Markov et al. 2011; Muckli and Petro 2013; Qiao et al. 2020). The visual association cortex is usually segregated into two anatomically and functionally distinct pathways: a ventral occipitotemporal pathway subserving object perception (“what”), and a dorsal occipitoparietal pathway subserving object localization (“where”) and visuomotor control (“how”) (Creem and Proffitt 2001; Freud et al. 2016).
It is generally accepted that information processing in the brain involves complex interactions among broadly distributed regions. Resting-state functional magnetic resonance imaging (rs-fMRI) (Logothetis et al. 2001; Rosen and Savoy 2012) allows researchers to characterize the interactions by providing a useful resting-state functional connectivity (rsFC) measure reflecting the temporal correlation of blood-oxygen-level-dependent (BOLD) signals between spatially distinct brain regions (Biswal et al. 1995; Yeo et al. 2011). By leveraging this technique, previous studies have demonstrated that different visual subregions have distinct rsFC patterns (Shen et al. 2010; Tosoni et al. 2015; Park et al. 2018), supporting the concept that different functions of visual subregions can be attributable to variations in their connection profiles. Moreover, the rsFC patterns of visual subregions have been proved to be highly stereotyped across individuals (Mueller et al. 2013), implying a conserved molecular program responsible for their development. Nevertheless, the molecular mechanisms underlying the rsFC variations across visual subregions have yet to be determined. More broadly, as there is a diverse and convergent body of data showing selective functional dysconnectivity of visual subregions in multiple neuropsychiatric conditions (e.g., schizophrenia (Reavis et al. 2020), autism (Villalobos et al. 2005; Borowiak et al. 2020), depression (Le et al. 2017; Su et al. 2018), anxiety (Liao et al. 2010; Arnold Anteraper et al. 2014), and migraine (Burke et al. 2020)), elucidating such mechanisms may not only have theoretical implications for further understanding disease neuropathology, but also have clinical relevance for developing new prevention and treatment programs targeted to underlying molecular substrates.
With the recent availability of comprehensive, brain-wide gene expression atlases such as the Allen Human Brain Atlas (AHBA) (Hawrylycz et al. 2012) comes the unprecedented capacity to link microscopic molecular function to macroscopic brain organization. The integration of brain-wide transcriptomic and imaging data has given rise to the nascent field of neuroimaging transcriptomics, whose goal is to offer new insights into the genetic basis of various neuroimaging phenotypes (Fornito et al. 2019). Two most commonly used approaches in this field are transcription-neuroimaging spatial association analysis and gene co-expression analysis. The former refers to relating spatial variations in gene expression to spatial variations in a given neuroimaging phenotype. The latter is to relate region-to-region correlations of gene expression to structural or functional connectivity measures between pairs of brain regions. In this context, there have been recent attempts to unveil the transcriptional architecture accounting for rsFC in healthy populations using these approaches (Hawrylycz et al. 2015; Richiardi et al. 2015; Krienen et al. 2016; Vertes et al. 2016; Anderson et al. 2018; Zhu et al. 2021). These efforts, in combination with the state-of-the-art rsFC analysis methods and data-mining techniques, have identified genes with spatial profiles of expression that track anatomical variations in rsFC measures as well as have established the correspondence between inter-regional gene co-expression and inter-regional rsFC. Despite this previous work, there is a paucity of research examining the transcriptional correlates underlying the rsFC patterns of different visual subregions. A full explanation of these underpinnings might benefit from an elaborate transcriptome-neuroimaging correlation analysis in independent discovery and validation datasets using tailored methodologies including standardized gene expression data processing and rsFC measurement based on more fine-grained segments of the visual cortex.
Toward this goal, rs-fMRI data were obtained from a discovery dataset and two independent cross-scanner, cross-race validation datasets. Brain gene expression data from the AHBA dataset were processed using a newly proposed standardized pipeline (Arnatkeviciute et al. 2019). Seed-based rsFC was calculated based on 22 fine-grained occipital visual subregions derived from the Human Brainnetome Atlas (Fan et al. 2016), a new brain atlas constructed using a connectivity-based parcellation framework. Then, we identified genes associated with rsFC of each visual subregion using transcriptome-neuroimaging spatial correlation analyses, followed by performance of gene functional annotation. A schematic overview of the study design and analysis pipeline is provided in Figure 1. We hypothesized there would be commonalities and differences in the identified genes and their functional features across visual subregions.
A schematic overview of the study design and analysis pipeline. First, we defined 22 subregions for the bilateral occipital visual cortex based on the Human Brainnetome Atlas. For a given subregion, we calculated its rsFC to each brain tissue sample by computing Pearson’s correlation coefficients between BOLD time courses of this subregion and each sample, resulting in a subregion-samples rsFC vector for each subject. Then, we obtained a sample × subject rsFC matrix for each visual subregion. Second, we used the newly proposed pipeline to process gene expression data, resulting in normalized expression data of 5013 genes for 1280 tissue samples. Third, we performed transcriptome-neuroimaging spatial correlation analyses to identify genes associated with rsFC of each visual subregion. In the discovery experiment, cross-sample Pearson’s correlations between gene expression and rsFC were performed in a gene-wise manner, yielding 5013 correlation coefficients for each subject. Next, genes with significant correlations in more than 90% of subjects were used for further validation analyses. In the validation experiment, the same analysis procedure was conducted for the CNP and SALD datasets, respectively. Only genes consistently identified in both discovery and validation experiments were considered as rsFC-related genes. Finally, a range of annotation analyses were carried out for the identified rsFC-related genes including functional enrichment, tissue-specific expression, cell type-specific expression, temporal-specific expression, and PPI. Abbreviations: rsFC: resting-state functional connectivity; BOLD, blood-oxygen-level-dependent; AHBA, Allen Human Brain Atlas; CNP, Consortium for Neuropsychiatric Phenomics; SALD, Southwest University Adult Lifespan Dataset; PPI, protein–protein interaction.
Materials and Methods
Participants
Our study included a discovery dataset along with two independent cross-scanner and cross-race validation datasets. The discovery participants were healthy adults of Chinese Han and right handedness, recruited from the local universities and community through poster advertisements. Exclusion criteria included neuropsychiatric or severe somatic disorder, a history of head injury with loss of consciousness, pregnancy, MRI contraindications, and a family history of psychiatric illness among first-degree relatives. This study was approved by the ethics committee of The First Affiliated Hospital of Anhui Medical University. Written-informed consent was obtained from all participants after they had been given a complete description of the study. The validation samples were from two publically available datasets: the Consortium for Neuropsychiatric Phenomics (CNP, https://openneuro.org/datasets/ds000030/versions/1.0.0) (Poldrack et al. 2016), and the Southwest University Adult Lifespan Dataset (SALD, http://dx.doi.org/10.15387/fcp_indi.sald) (Wei et al. 2018). Of note, we solely selected the healthy adults from the cross-disorder CNP dataset. Full details regarding the two validation samples (e.g., ethics, informed consent, inclusion and exclusion criteria, among others) have been described in the data descriptor publications (Poldrack et al. 2016; Wei et al. 2018). To rule out the potential influence of neurodevelopment and neurodegeneration, all the participants were restricted to an age range of 18–60 years. Additionally, participants with poor image quality or excessive head motion during scanning were excluded. This brought the final samples used in this study to 361 in the discovery dataset, 103 in the CNP dataset, and 329 in the SALD dataset. Details of the demographic data of the three datasets are presented in Supplementary Table 1 in the Supplementary Materials.
Image Acquisition
MRI data of the discovery sample were obtained using the 3.0-Tesla General Electric Discovery MR750w scanner, and those of the validation samples were acquired using the 3.0-Tesla Siemens Trio scanners. Details of the resting-state fMRI protocols for the three datasets can be found in Supplementary Table 2 in the Supplementary Materials. Of note, during the fMRI scanning, the participants from the discovery and SALD datasets were asked to close their eyes, while the participants from the CNP dataset were asked to open their eyes. There is empirical evidence that rsFC of visual regions are different between eyes open versus closed conditions (Xu, Huang, et al. 2014a; Agcaoglu et al. 2019). Therefore, we used fMRI datasets collected in two eyes conditions to exclude the influence of eyes conditions.
Image Processing
Resting-state BOLD data were preprocessed using Statistical Parametric Mapping (SPM12, http://www.fil.ion.ucl.ac.uk/spm) and Data Processing & Analysis for Brain Imaging (DPABI, http://rfmri.org/dpabi) (Yan et al. 2016). The first several volumes (discovery: 10, CNP: 5, SALD: 10) for each participant were discarded to allow the signal to reach equilibrium and the participants to adapt to the scanning noise. The remaining volumes were corrected for the acquisition time delay between slices. Then, realignment was performed to correct the motion between time points. Head motion parameters were computed by estimating the translation in each direction and the angular rotation on each axis for each volume. All participants’ BOLD data were within the defined motion thresholds (i.e., translational or rotational motion parameters less than 2.0 mm or 2.0°). We also calculated frame-wise displacement (FD), which indexes the volume-to-volume changes in head position. Several nuisance covariates (the linear drift, the estimated motion parameters based on the Friston-24 model, the spike volumes with FD > 0.5, the white matter signal, and the cerebrospinal fluid signal) were regressed out from the data. The datasets were then band-pass filtered using a frequency range of 0.01–0.1 Hz. In the normalization step, individual structural images were firstly co-registered with the mean functional image; then the transformed structural images were segmented and normalized to the Montreal Neurological Institute (MNI) space using the diffeomorphic anatomical registration through the exponentiated Lie algebra (DARTEL) technique (Ashburner 2007). Finally, each filtered functional volume was spatially normalized to the MNI space using the deformation parameters estimated during the above step and resampled into a 3-mm cubic voxel.
Definition of Occipital Visual Subregions
Subregions of the occipital visual cortex were defined according to a recent connectivity-based parcellation study using multimodal neuroimaging techniques (Fan et al. 2016). These subregions are parcellated based on their specific anatomical and functional connectivity patterns. In each hemisphere, the medial occipital cortex was divided into caudal lingual gyrus (cLinG), rostral cuneus gyrus (rCunG), caudal cuneus gyrus (cCunG), rostral lingual gyrus (rLinG) and ventromedial parietooccipital sulcus (vmPOS), and the lateral occipital cortex into middle occipital gyrus (mOccG), middle temporal visual area (V5/MT+), occipital polar cortex (OPC), inferior occipital gyrus (iOccG), medial superior occipital gyrus (msOccG) and lateral superior occipital gyrus (lsOccG). Thus, we defined a total of 22 subregions for the bilateral occipital visual cortex (Fig. 2).
Illustration of occipital visual subregions. Abbreviations: L, left; R, right; cLinG, caudal lingual gyrus; rCunG, rostral cuneus gyrus; cCunG, caudal cuneus gyrus; rLinG, rostral lingual gyrus; vmPOS, ventromedial parietooccipital sulcus; mOccG, middle occipital gyrus; V5/MT+, middle temporal visual area; OPC, occipital polar cortex; iOccG, inferior occipital gyrus; msOccG, medial superior occipital gyrus; lsOccG, lateral superior occipital gyrus.
Brain Gene Expression Data Processing
Brain gene expression information was obtained from the downloadable AHBA dataset (http://www.brain-map.org) (Hawrylycz et al. 2012, 2015). The dataset was derived from six human post-mortem donors (Supplementary Table 3 in the Supplementary Materials). Custom 64 K Agilent microarrays were used to measure the expression of more than 20 000 genes at 3702 spatially distinct brain tissue samples. We used the newly proposed pipeline to process gene expression data for transcriptome-neuroimaging associations in this study (Arnatkeviciute et al. 2019). Specifically, we first updated the probe-to-gene annotations based on the latest available information from the National Center for Biotechnology Information (NCBI) using the Re-Annotator package (Arloth et al. 2015). With intensity-based filtering, we excluded probes that did not exceed the background noise in at least 50% of samples across all donors. As multiple probes were used to measure expression level of a single gene, we further used the RNA-seq data as a reference to select probes. After excluding genes that do not overlap between RNA-seq and microarray datasets, we calculated the correlations between microarray and RNA-seq expression measures for the remaining genes. After excluding probes with low correlations (r < 0.2), a representative probe for a gene was selected based on the highest correlation to the RNA-seq data. In this study, we only included the tissue samples in the left cerebral cortex. For one, all six donors had expression data in the left hemisphere, but only two donors had samples in the right hemisphere. For another, the inclusion of subcortical samples might introduce potential biases because of the great differences in gene expression between cortical and subcortical regions (Hawrylycz et al. 2012). To account for potential between-sample differences and donor-specific effects in gene expression, we performed both within-sample cross-gene and within-gene cross-sample normalization by using the scaled robust sigmoid normalization method. Differential stability (DS) is a measure of consistent regional variation across donor brains. Previous research has reported that genes with high DS scores demonstrate more consistent spatial expression patterns between donors and are enriched for brain-related biological function (Hawrylycz et al. 2015). As gene expression conservation across subjects is a prerequisite for the transcriptome-neuroimaging spatial correlations, we only selected genes with relatively more conserved expression patterns for analysis. To realize this goal, we ranked the genes by their DS values and chose the 50% of the highest DS genes for the main analysis. After these processing procedures, we obtained normalized expression data of 5013 genes for 1280 tissue samples, resulting in a sample × gene matrix of 1280 × 5013.
Calculating rsFC of Occipital Visual Subregions to Tissue Samples
To calculate rsFC of a given subregion to brain tissue samples from the AHBA, we drew a spherical region (radius = 3 mm) centered at the MNI coordinate of each tissue sample and extracted the mean BOLD time course of voxels within each sphere, followed by calculation of Pearson’s correlation coefficient between the mean time courses of this subregion and each sample in the brain. Then, the Pearson’s correlation coefficients were converted to Fisher’s Z-scores to improve normality, resulting in a subregion-samples rsFC vector for each subject. Finally, we obtained a sample × subject rsFC matrix for each occipital visual subregion.
Identifying Genes Related to rsFC of Occipital Visual Subregions
Genes associated with rsFC of each occipital visual subregion were determined separately using the following spatial correlation analyses. In the discovery experiment, cross-sample (1280 samples) Pearson’s correlations between gene expression and rsFC were performed in a gene-wise manner (5013 genes), yielding 5013 correlation coefficients for each subject. Next, genes with significant correlations (Bonferroni correction for the number of genes, P < 0.05/5013 = 9.974 × 10−6) in more than 90% of subjects were used for further validation analyses. In the validation experiment, the same analysis procedure was conducted for the CNP and SALD datasets, respectively. Only genes consistently identified in both discovery and validation experiments were considered as rsFC-related genes. For ease of interpretability, we pooled genes associated with the homotopic subregions (e.g., the left and right cLinG) for subsequent analysis.
To further test whether the number of the identified rsFC-related genes was significantly greater than random level, a permutation test was pursued to establish the significance of our results. As both transcriptional and rsFC data are spatially autocorrelated, that is, nearby anatomical regions tend to have more similar patterns of gene expression and rsFC than spatially distant regions, the standard non-parametric null (i.e., randomly shuffling the sample labels) is strongly violated by the spatial autocorrelation of brain maps, yielding increased family-wise error rates (Markello and Misic 2021). To address this issue, we used a spatially-constrained null model proposed by Burt et al. (2020) to conduct the permutation test. This method is implemented in an open-access, Python-based software package, BrainSMASH: Brain Surrogate Maps with Autocorrelated Spatial Heterogeneity (https://github.com/murraylab/brainsmash). It can simulate volumetric surrogate brain maps that preserve the spatial autocorrelation using three-dimensional Euclidean distance between regions. To correct the spatial autocorrelation in both transcriptional and brain imaging data, we used this method to generate spatial autocorrelation-preserving surrogate maps for each gene and rsFC of each subregion. These surrogate maps were used to re-identify the rsFC-related genes using exactly the same method as described above. Then, we repeated this procedure 1000 times and recorded the number of genes identified in each test to build a null distribution. Finally, we compared the number of genes identified using the real data with this null distribution to determine whether our results were different from random.
Gene Enrichment Analysis
For each subregion, a range of enrichment analyses were carried out for the identified rsFC-related genes. First, these genes were functionally annotated using the ToppGene portal (https://toppgene.cchmc.org/) (Chen et al. 2009). Gene ontology (GO) was used to determine their biological functions including molecular functions (MFs), biological processes (BPs), and cellular components (CCs). The pathway and disease databases were used to determine the biological pathways and related diseases. Second, we used online tissue-specific expression analysis (TSEA) (http://genetics.wustl.edu/jdlab/tsea/) and cell type-specific expression analysis (CSEA) tools (http://genetics.wustl.edu/jdlab/csea-tool-2/) (Dougherty et al. 2010) to conduct tissue, cell type, and temporal specific expression analyses. These specific expression analyses could help to determine the specific tissues, cortical cell types, and developmental stages in which the rsFC-related genes were overrepresented. A specificity index probability (pSI = 0.05) was used to determine how likely a gene was to be specifically expressed (Xu, Wells, et al. 2014b). Finally, we examined the overlap between rsFC-related genes found in the current study and those reported in previous rsFC studies (Richiardi et al. 2015; Anderson et al. 2018), using 20 737 genes with unique Entrez IDs in the AHBA as the background list. Fisher’s exact tests were used to assess the significance of the above-mentioned enrichment analyses. Multiple testing was corrected using the Benjamini and Hochberg method for false discovery rate (FDR-BH correction) with a corrected P value of 0.05.
Capturing the Behavioral Relevance of Genes Related to rsFC of Occipital Visual Subregions
To further explore the behavioral relevance of the identified rsFC-related genes, we examined the associations of gene expression with behavioral domains via the Neurosynth (https://neurosynth.org/), a well-validated and publicly available platform for meta-analysis of functional neuroimaging literature (Yarkoni et al. 2011). The Neurosynth database provides a range of activation maps of behavioral terms capturing conceptually distinct aspects of human behavior. Cross-sample spatial correlation analyses were performed between expression measures of the identified rsFC-related genes and activation values in each behavioral term map. For each behavioral term, we obtained a set of correlation coefficients corresponding to the number of the rsFC-related genes. A positive correlation coefficient indicates that a brain region with higher gene expression tends to show greater neural activation, while a negative correlation coefficient means that a brain region with lower gene expression tends to show greater neural activation. Thus, both positive and negative correlation coefficients indicate that a gene contributes to a behavioral term. To avoid biases due to offset, we averaged the absolute values of these correlation coefficients (|r|mean) to index the extent to which this set of genes was linked to each behavioral term. Consequently, the behavioral terms were ordered based on their |r|mean and those with higher |r|mean were identified to characterize the behavioral relevance of genes associated with rsFC of each visual subregion. Here, a threshold of |r|mean > 0.2 was used for visualization and interpretation.
Protein–Protein Interaction Analysis
Protein–protein interaction (PPI) analysis was performed by using STRING v11.0 (https://string-db.org/) to construct PPI networks and identify hub genes from the rsFC-related genes. Hub genes were defined by the top 10% of the node degree in the PPI networks with a highest confidence interaction score of 0.9. Besides, the Human Brain Transcriptome database (http://hbatlas.org/) was used to characterize the spatial–temporal expression trajectory of hub genes with the highest node degree.
Sensitivity Analysis
Considering the potential influences of different processing methods for gene expression and brain imaging data, we conducted several sensitivity analyses to test for robustness of our results. First, we chose the 50% of the highest DS genes to focus our efforts on genes with relatively more conserved expression patterns across six donors in the main analysis. To assess the effect of different DS threshold selections, we repeated our analysis with two other DS cutoff thresholds (top 40% and 60%). Second, since global signal regression (GSR) has long been a controversial topic in rs-fMRI analyses (Murphy and Fox 2017), we re-calculated rsFC of each visual subregion using BOLD data with GSR and then repeated the entire analyses.
Results
The rsFC Patterns of Occipital Visual Subregions
The rsFC maps of occipital visual subregions in one discovery dataset and two validation datasets are depicted in Supplementary Figure 1 in the Supplementary Materials. Each subregion showed similar rsFC patterns across the discovery and validation datasets. While different visual subregions exhibited distinct rsFC patterns, they consistently had strong connections with the sensory systems including the visual, somatosensory, and auditory areas.
Genes Related to rsFC of Occipital Visual Subregions
For each visual subregion, we performed gene expression-rsFC spatial correlations analyses at the individual level in the discovery and validation datasets. Only genes with significant correlations (P < 0.05, Bonferroni corrected) in more than 90% of subjects in each dataset were considered as rsFC-related genes. As a result, we found eight sets of genes associated with rsFC of eight visual subregions (see Supplementary File 1 for the pooled genes and Supplementary File 2 for genes related to rsFC of left and right visual subregions separately). Specifically, expression measures of 52 genes were correlated with rsFC of the cLinG, 971 with the rCunG, 30 with the cCunG, 406 with the rLinG, 501 with the vmPOS, 95 with the mOccG, 2 with the V5/MT+, and 172 with the msOccG. Moreover, the spatially-constrained permutation tests showed that none out of 1000 permutations resulted in more rsFC-related genes than those identified using the real data for these subregions (Pperm < 0.001), indicating that our results were different from random. Scatter plots of the representative correlations between gene expression and rsFC of each visual subregion in the discovery dataset are shown in Supplementary Figure 2 in the Supplementary Materials.
Gene Functional Enrichment
To characterize the biological functions, pathways, and diseases of the rsFC-related gene sets, we performed functional enrichment analyses using the ToppGene portal. The results of enrichment analyses for seven visual subregions with exception of the non-significant mOccG are listed in Supplementary File 3 and are depicted in Figure 3. Overall, the genes related to rsFC of different visual subregions showed distinct enrichment results. Specifically, the genes related to rsFC of the cLinG were enriched for MFs including oxysterol binding; for pathways including antiarrhythmic drug pathways such as flecainide pathway and antihypertensive drug pathways such as felodipine pathway; and for diseases including ophthalmic diseases such as optic nerve hypoplasia and neuropsychiatric diseases such as epilepsy and psychosis (Fig. 3A). The genes related to rsFC of the rCunG were enriched for MFs including nuclear receptor activity and rho GDP-dissociation inhibitor activity, BPs including nervous system development, regulation of neurotransmitter levels, synaptic signaling and neurogenesis, and CCs including synapse, neuron projection, neuronal cell body, dendrite and axon; for pathways including neuronal system, glutamatergic synapse, and signaling pathways such as calcium signaling pathway and cAMP signaling pathway; and for diseases including neuropsychiatric diseases such as schizophrenia, impaired cognition, alcoholic intoxication, and bipolar disorder (Fig. 3B). The genes related to rsFC of the cCunG were enriched for MFs including oxysterol binding, sulfur carrier activity, nitric-oxide synthase inhibitor activity, sodium channel regulator activity, calcium activated cation channel activity, and protein kinase C activity (Fig. 3C). The genes related to rsFC of the rLinG were enriched for MFs including rho GDP-dissociation inhibitor activity, BPs including regulation of neurotransmitter levels and nitric oxide metabolic process, and CCs including synapse (Fig. 3D). The genes related to rsFC of the vmPOS were enriched for MFs including rho GDP-dissociation inhibitor activity, BPs including protein autophosphorylation, and CCs including synapse, neuron projection, neuronal cell body, axon, intrinsic component of plasma membrane, GABA-ergic synapse and growth cone membrane; and for pathways including calcium signaling pathway, cAMP signaling pathway, melanogenesis, neuronal system, neuroactive ligand-receptor interaction, and ion channel transport (Fig. 3E). The genes related to rsFC of the V5/MT+ were primarily enriched for BPs including detection of visible light; and for pathways including inactivation, recovery, regulation of the phototransduction cascade (Fig. 3F). The genes related to rsFC of the msOccG were mainly enriched for MFs including ion channel activity (Fig. 3G).
Functional enrichment of the genes related to rsFC of occipital visual subregions. For each bubble chart corresponding to a visual subregion, the x-axis denotes the rich factor and the y-axis denotes items from the GO, pathway and disease databases. The rich factor refers to the ratio of the number of rsFC-related genes annotated to the item to the number of all genes annotated to the item. The bubble size represents the number of genes overlapping with those belonging to each item, and the bubble color represents the −log10(P) with the P value corrected by the FDR-BH method. Abbreviations: rsFC, resting-state functional connectivity; GO, gene ontology; cLinG, caudal lingual gyrus; rCunG, rostral cuneus gyrus; cCunG, caudal cuneus gyrus; rLinG, rostral lingual gyrus; vmPOS, ventromedial parietooccipital sulcus; V5/MT+, middle temporal visual area; msOccG, medial superior occipital gyrus.
Tissue, Cell Type, and Temporal Specific Expression
The genes related to rsFC of the rCunG, rLinG, vmPOS, mOccG, and msOccG showed specific expression in brain tissue (Fig. 4A and Supplementary File 4, sheet “Tissues”). For gene sets that were specifically expressed in brain tissue, we further investigated the cortical cell types and developmental stages in which these gene sets were specifically expressed. As a consequence, we found that the genes related to rsFC of the rCunG, rLinG, vmPOS, and msOccG were specifically expressed in multiple types of neurons and immune cells, while those of the mOccG were only specifically expressed in Pnoc+ neurons (Fig. 4B and Supplementary File 4, sheet “Cell types”). The temporal-specific expression analyses showed that the genes related to rsFC of the rCunG, rLinG, and vmPOS were preferentially expressed during middle late childhood, adolescence, and young adulthood; those of the mOccG were preferentially expressed during middle late childhood; and those of the msOccG were preferentially expressed during neonatal early infancy and middle late childhood (Fig. 4C and Supplementary File 4, sheet “Developmental stages”).
Specific expression of the genes related to rsFC of occipital visual subregions. (A) Tissue-specific expression. The inner location indicates a greater −log10(P) (FDR-BH corrected) which represents a higher statistical significance. Colored fonts represent tissues in which genes were significantly specifically expressed. (B) Cell type-specific expression. (C) Temporal-specific expression. Y axis is −log10(P); if −log10(P) > 5, then −log10(P) = 5. The red dashed line represents the statistical significance threshold of P < 0.05 (FDR-BH corrected). *P < 0.05, FDR-BH corrected. Abbreviations: rsFC, resting-state functional connectivity; cLinG, caudal lingual gyrus; rCunG, rostral cuneus gyrus; cCunG, caudal cuneus gyrus; rLinG, rostral lingual gyrus; vmPOS, ventromedial parietooccipital sulcus; mOccG, middle occipital gyrus; V5/MT+, middle temporal visual area; msOccG, medial superior occipital gyrus; Myeli, myelinating oligodendrocytes; Proge, oligodendrocyte progenitor cells; Immu, immune cells, Astro, astrocytes.
Overlap with rsFC-Related Genes in Previous Studies
Fisher’s exact tests revealed that rsFC-related genes found in the current study significantly overlapped with those reported in previous rsFC studies (P < 0.05, FDR-BH corrected; Supplementary Tables 4 and 5 in the Supplementary Materials), indicating some degree of rsFC specificity of our findings.
Behavioral Relevance of Genes Related to rsFC of Occipital Visual Subregions
By linking gene expression with behavioral domains via the Neurosynth, we found that genes related to rsFC of visual subregions within both the medial and lateral occipital cortex (the rCunG, rLinG, vmPOS, and msOccG) were correlated with shared behavioral terms including sensation, language, motion, perception, emotion, and attention; in addition, those within the lateral occipital cortex (the msOccG) were correlated with unique behavioral terms such as personality, semantic memory, imagery, and social cognition, among others (Fig. 5). Remarkably, behavioral terms pertaining to visual processes showed the strongest correlations, corroborating the prominent visual involvement of the identified genes.
Correlations between rsFC-related genes and behavioral terms from the Neurosynth. Four visual subregions were selected because their rsFC-related genes could construct PPI networks with statistical significance. The rCunG, rLinG, and vmPOS are within the medial occipital cortex and the msOccG is within the lateral occipital cortex. The coordinate values represent the mean absolute correlations between expression measures of the identified rsFC-related genes and activation values in each behavioral term map. Abbreviations: rsFC, resting-state functional connectivity; PPI, protein–protein interaction; rCunG, rostral cuneus gyrus; rLinG, rostral lingual gyrus; vmPOS, ventromedial parietooccipital sulcus; msOccG, medial superior occipital gyrus.
PPI Networks and Hub Genes
PPI analyses demonstrated that genes related to rsFC of four visual subregions (rCunG, rLinG, vmPOS, and msOccG) could construct PPI networks with statistical significance (Fig. 6). Genes with top 10% of the node degree in each PPI network were identified as the hub genes. As such, there were 34, 11, 14, and 2 hub genes involved in the PPI networks constructed by the gene sets associated with the rCunG, rLinG, vmPOS, and msOccG, respectively (Fig. 7). We also characterized the spatial–temporal expression trajectory of four hub gene with the highest node degree for each visual subregion (GNB4 for the rCunG, HIST1H2BK for the rLinG, ADCY7 for the vmPOS, and ADAMTS3 for the msOccG) (Fig. 6).
PPI networks. PPI networks with statistical significance were constructed by the genes related to rsFC of four visual subregions including rCunG (A), rLinG (B), vmPOS (C), and msOccG (D). The subregions, whole PPI networks, PPI subnetworks centered at the hub genes with the highest node degree, and spatial–temporal specific expression curves of the hub genes with the highest node degree are shown from the left to right panels. The P value denotes the statistical significance of how likely the proteins encoded by the input genes are connected to construct a network. Abbreviations: PPI, protein–protein interaction; rsFC, resting-state functional connectivity; rCunG, rostral cuneus gyrus; rLinG, rostral lingual gyrus; vmPOS, ventromedial parietooccipital sulcus; msOccG, medial superior occipital gyrus.
Circos map showing hub genes in the PPI networks. The orange, purple, pink and blue tracks represent 34, 11, 14, and 2 hub genes involved in the PPI networks constructed by the gene sets associated with the rCunG, rLinG, vmPOS and msOccG, respectively. For each track, the outermost ring shows the symbols of the hub genes. The second circle represents −log10(P) (Bonferroni corrected) of the transcriptome-neuroimaging spatial correlations for each subject (a dot or triangle). The significant P value is shown as pink dot (positive correlation) or blue triangle (negative correlation), and the non-significant P value is shown as gray. The third circle represents the percentage of subjects with consistently significant correlations in each dataset sorted clockwise as the discovery, CNP, and SALD datasets. The gray lines from the outside to inside represent the percentages of 90%, 92%, 94%, 96%, 98%, and 100%, respectively. The innermost ring represents the average correlation coefficients between gene expression and rsFC of all subjects from the three datasets, with a darker red (blue) color indicating a greater positive (negative) correlation coefficient. Abbreviations: PPI, protein–protein interaction; rsFC, resting-state functional connectivity; rCunG, rostral cuneus gyrus; rLinG, rostral lingual gyrus; vmPOS, ventromedial parietooccipital sulcus; msOccG, medial superior occipital gyrus; CNP, the Consortium for Neuropsychiatric Phenomics; SALD, the Southwest University Adult Lifespan Dataset.
Sensitivity Analysis
To determine the effect of different DS threshold selections, we used two other DS cutoff thresholds (top 40% and 60%) during the brain gene expression data processing to obtain normalized expression measures of 4010 and 6016 genes, respectively. By repeating the gene expression-rsFC association analysis, we found substantial overlaps between the identified genes in the main analyses and those in the sensitivity analyses based on the thresholds of 40% (overlap ratio: 81.08–100%) and 60% (overlap ratio: 99.14–100%) (Supplementary Tables 6 and 7 in the Supplementary Materials and Supplementary File 5). In addition, using BOLD data with GSR led to a considerable reduction in the identified rsFC-related genes, while most of the identified significant genes (84.97–100%) were included in those found in our main analyses (Supplementary Table 8 in the Supplementary Materials and Supplementary File 6).
Discussion
The present study is the first to investigate the molecular mechanisms underlying different rsFC patterns of occipital visual subregions by applying transcriptome-neuroimaging correlation analyses to gene expression data from the AHBA dataset and rs-fMRI data from three cross-scanner, cross-race datasets. We found rsFC of eight visual subregions were associated with expression measures of eight gene sets, which were specifically expressed in brain tissue and showed the strongest correlations with visual behavioral processes. More importantly, there was a significant divergence in the identified rsFC-related genes and their functional features between medial and lateral visual subregions. First, there were more genes found for medial visual subregions than lateral visual subregions. Second, contrasting with those of lateral visual subregions enriched for relatively specialized biological functions and pathways, genes related to rsFC of medial visual subregions were enriched for more diverse biological functions and pathways. Third, medial visual subregions associated genes were enriched for multiple neuropsychiatric diseases, while this enrichment was absent for lateral visual subregions associated genes. Fourth, medial visual regions associated genes were specifically expressed in multiple types of neurons and immune cells, yet lateral visual region associated genes were specifically expressed in a single type of neurons. Fifth, medial visual subregions associated genes were preferentially expressed during the middle and late stages of cortical development, while lateral visual subregions associated genes were preferentially expressed during the early and middle stages. Finally, although genes associated with both medial and lateral visual subregions were correlated with shared behavioral terms, lateral visual subregion associated genes were correlated with unique high-order cognitive terms. Collectively, these findings confirm our hypothesis of commonalities and differences in the identified genes and their functional features across visual subregions.
Previous studies have documented that the synchronized patterns of gene expression could contribute to the formation and maintenance of brain functional connectivity (Richiardi et al. 2015; Zhu et al. 2021). In support of this notion, our results demonstrated that the transcriptional profiles of genes were associated with the rsFC patterns of visual subregions. Furthermore, the rsFC-related genes found in the current study showed significant overlaps with those reported in prior studies, indicating some degree of rsFC specificity of these findings. However, for some visual subregions within the lateral occipital cortex (the OPC, iOccG, and lsOccG), we found no significant spatial correlations between their rsFC and gene expression. The reasons for this are unknown and are likely to be multifactorial. One plausible explanation lies in the fact that these subregions belong to the visual association cortex. Empirical evidence suggests higher individual variability in rsFC in heteromodal association cortex and lower variability in unimodal cortex (Mueller et al. 2013). Thus, one may speculate that the muted gene expression-rsFC association may arise from highly variable rsFC of these visual association regions.
In the current study, one important finding is the significant divergence in the identified rsFC-related genes between medial and lateral visual subregions. First, there were more genes found for medial visual subregions than lateral visual subregions. The medial primary visual cortex is thought to entail basic neural processes including receiving visual information and providing first-order processing, while the lateral visual association cortex entails complex neural processes. Our results are consistent with the previous hypothesis that brain structures subserving basic processes may be under strong genetic control, while the association cortex may be less genetically influenced during development, enabling more variable impact of environmental factors that ultimately leads to the diversity of neural connections beyond their genetic determination (Brun et al. 2009). Second, contrasting with those of lateral visual subregions enriched for relatively specialized biological functions and pathways, genes related to rsFC of medial visual subregions were enriched for more diverse biological functions and pathways. For example, in terms of cell communication, the msOccG associated genes were enriched for ion channel activity, such as voltage-gated potassium channel (coded by KCNA3, KCNAB3, and KCNMB4), sodium channel (coded by SCN3B) and calcium channel (coded by CACNA1H). These ion channels are crucial mediators of neuronal excitability and are actively implicated in cellular and molecular signaling pathway (Eijkelkamp et al. 2012; Shah and Aizenman 2014; Kumar et al. 2016). By contrast, medial visual subregions (the rCunG, cCunG, rLinG, and vmPOS) associated genes were enriched not only for ion channel activity, but also for other cell communicational types including synaptic signaling, an important information exchange manner between neurons (Biederer et al. 2017), and nitric oxide, a highly reactive gas with considerable diffusion power that can serve as a multifunctional signaling molecule (Picon-Pages et al. 2019). With respect to biological processes, the V5/MT+ associated genes were only enriched for detection of visible light, while the rCunG associated genes were enriched for nervous system development and neurogenesis. At the molecular level, medial visual subregions associated genes were involved in molecular functions that affect neuronal development such as nuclear receptor activity and cAMP signaling pathway. Nuclear receptors are a family of transcription factors (Evans and Mangelsdorf 2014; Lazar 2017) that can regulate the downstream gene transcription (Sun and Xu 2020). For instance, the nuclear receptor subfamily 3 group C member 1 (NR3C1) encodes the glucocorticoid receptor (Weikum et al. 2017), which can be combined with glucocorticoids to alter neuronal excitability (Groeneweg et al. 2012), regulate dendritic spine development and plasticity (Liston and Gan 2011), and remodel the architecture and function of dendrite and synapse (Chen et al. 2017). cAMP signaling affects a large number of the developmental processes in neural system including cell differentiation, axon outgrowth, response to guidance molecules, or modulation of neuronal connectivity (Nicol and Gaspar 2014). Furthermore, previous literature has highlighted the pivotal role of cAMP signaling in shaping the neuronal connectivity of visual system (Seol et al. 2007; Nicol and Gaspar 2014). Finally, medial visual subregions (the cLinG and rCunG) associated genes were enriched for neuropsychiatric diseases including epilepsy, schizophrenia, impaired cognition, alcoholic intoxication and bipolar disorder, while this enrichment was absent for lateral visual subregions associated genes. Of note, among subregions within the medial occipital cortex, we found that the rCunG had the maximum rsFC-related genes enriched for the most diverse biological functions and pathways as well as multiple neuropsychiatric diseases. From the view of brain imaging, the cuneus is a prominent functional hub in the human brain (Tomasi and Volkow, 2011a, 2011b), and its activity and connectivity have been shown to be implicated in many neuropsychiatric diseases such as depression (Marwood et al. 2018; Ma et al. 2019; Zhu et al. 2020), autism (Rolls et al. 2020), cognitive impairment (Xu et al. 2020), and alcohol dependence (Wang et al. 2018). Our present observation emphasizes the critical role of the cuneus in the human visual cortex from the perspective of molecular mechanisms underlying its rsFC.
Tissue-specific expression analysis showed that the genes related to rsFC of most visual subregions showed specific expression in brain tissue, which confirms the credibility of our results. Despite specific expression in brain tissue, gene sets related to rsFC of medial and lateral visual regions diverged on cell types and developmental stages for which they were enriched. Medial visual regions (the rCunG, rLinG, and vmPOS) associated genes were specifically expressed in multiple types of neurons and immune cells, yet lateral visual subregion (the mOccG) associated genes were specifically expressed in a single type of neurons (Pnoc+ neurons). The overarching specific expression in neurons is congruent with prior findings that underscore the importance of neurons in mediating the genetic effects on rsFC (Richiardi et al. 2015; Anderson et al. 2018). In addition, it is largely known that immune cells such as microglial cells have the potential to regulate the development, structure, and function of neuronal network in the central nervous system (Garaschuk and Verkhratsky 2019). Accordingly, it is reasonable to assume that immune cells may exert their indirect effects on rsFC of some medial visual subregions via regulating neurons. In temporal-specific expression analysis, we found that medial visual subregions (the rCunG, rLinG, and vmPOS) associated genes were preferentially expressed during the middle and late stages of cortical development, a period when the brain undergoes dramatic remodeling and also a time when neuropsychiatric conditions often emerge (Paus et al. 2008). This result echoes the above-described finding that medial visual subregions associated genes were enriched for neuropsychiatric diseases. By contrast, lateral visual subregions (the mOccG and msOccG) were preferentially expressed during the early and middle stages of cortical development.
By linking gene expression with behavioral domains via the Neurosynth, interestingly, we found that behavioral terms pertaining to visual processes showed the strongest correlations with genes related to rsFC of occipital visual subregions, confirming the pronounced visual involvement of these identified genes. More importantly, we observed that genes associated with medial and lateral visual subregions differed in their related behavioral terms. Genes related to rsFC of visual subregions within both the medial and lateral occipital cortex (the rCunG, rLinG, vmPOS, and msOccG) were correlated with shared behavioral terms including sensation, language, motion, emotion and attention. In line with this finding, pilot work has established links between these shared behavioral terms and neural activity in the visual cortex (Kastner and Ungerleider 2000; Harris and Mrsic-Flogel 2013; Marques et al. 2018; Kanjlia et al. 2019; Kragel et al. 2019). Critically, lateral visual subregion (the msOccG) associated genes were uniquely correlated with high-order cognitive processes such as personality, semantic memory, imagery and social cognition, in favor of the hypothesis that functional variability in association cortex has important implications for the evolution of higher-order cognitive abilities (Mueller et al. 2013).
The visual subregions have been defined according to different parcellation strategies, such as the Human Brainnetome Atlas based on connectional architecture (Fan et al. 2016), the Brodmann atlas based on microscale cytoarchitecture (Zilles and Amunts 2010), and the PALS-12 atlas based on structural MRI volume (Van Essen 2005). The aim of the current study was to explore the molecular mechanisms underlying rsFC of different visual subregions. Thus, voxels within a subregion should have identical connectivity profiles, while different subregions should have distinct connectivity patterns. Accordingly, we preferred the connectivity-based Human Brainnetome Atlas in our analysis.
The GSR has been considered a controversial topic in rs-fMRI analyses (Murphy and Fox 2017). In the current study, we investigated genes associated with rsFC of visual subregions using BOLD data without and with GSR, respectively. Results showed a considerable reduction in the identified rsFC-related genes using data with GSR relative to those without GSR. As the global signal contains both non-neuronal confounds (e.g., head motion and physiological artifacts) and neuronal component, the application of GSR may remove not only global confounds, but also physiologically relevant neuronal signals. This has led to some speculation that the loss of neuronal information may hamper the identification of genes related to these neuronal activities, resulting in a reduction of the rsFC-related genes.
Several limitations should be noted in our research. First, the gene expression and brain imaging data were collected from different individuals. Adopting DS thresholding to focus our analyses on genes with more conservative expression profiles, coupled with strong evidence of conserved gene expression across human populations (Zeng et al. 2012; Hawrylycz et al. 2015), makes it feasible to conduct transcriptome-neuroimaging spatial correlation analysis. However, genes with great individual transcriptional differences would be inevitably missed. Second, we only included the tissue samples in the left cerebral cortex due to limited gene expression data in the right hemisphere and great gene expression differences between cortical and subcortical regions. However, the reduced tissue samples along with hemisphere and region selections may have biased our results. Third, the correlational nature of the analyses does not resolve causality. Future animal experiments are needed to validate and further interpret our preliminary results. Finally, although the AHBA dataset offers an excellent resource to investigate the transcriptome-neuroimaging relationship, we must consider that this gene expression information comes from only six individuals in light of the scarcity of donor brains.
In conclusion, this transcriptome-neuroimaging association study demonstrated that rsFC of different occipital visual subregions were modulated by distinct gene sets. Moreover, the gene sets associated with medial and lateral visual subregions diverged on their numbers and functional features including biological functions and pathways, diseases, cell type and temporal specific expression, and behavioral relevance. These findings may improve our understanding of the functional heterogeneity of the human visual cortex from the perspective of underlying genetic architecture.
Funding
National Natural Science Foundation of China (81801679, 81771817, and 82071905).
Notes
We thank the Allen Institute for Brain Science founders and staff who supplied the brain expression data. We also thank the subjects who contributed to this study. Conflict of Interest: The authors declare no conflict of interests.
Author Contributions
C. Z.: Conceptualization, Methodology, Investigation, Formal analysis, Visualization, Data curation, Writing-original draft; H. C.: Conceptualization, Methodology, Investigation,Visualization, Data curation; X. X.: Visualization, Data curation; Q. L.: Visualization, Data curation; X. L.: Visualization, Data curation; W. Z.: Visualization, Data curation; Y. Q.: Supervision, Data curation; J. Z.: Conceptualization, Methodology, Investigation, Supervision, Resources, Writing-review and editing; Y. Y.: Conceptualization, Supervision, Project administration; Resources, Funding acquisition.
