To analyze transcriptional profiles associated with major laminar and areal axes of cortical organization, laser microdissection (LMD) was used to selectively isolate individual cortical layers in ten discrete areas of the neocortex from two male and two female adult rhesus monkeys. As shown schematically in , these areas spanned primary sensorimotor cortices (S1, M1, A1, and V1), higher order visual areas (V2, MT and TE), and frontal cortical areas (DLPFC, OFC and ACG). In each cortical region, samples were isolated from layers definable on the basis of lightly stained Nissl sections used for the sample preparation, taking care to avoid layer boundaries. In most areas, 5 layers were isolated (L2, L3, L4, L5 and L6), although in M1, OFC and ACG no discernible L4 could be isolated. 8 layers were sampled in V1 () to include the functionally specialized and cytoarchitecturally distinct sublayers of L4 (4A, 4B, 4Ca and 4Cb). For a non-neocortical comparator data set, samples were also isolated from subfields of the hippocampus (CA1, CA2, CA3 and dentate gyrus) and from the magno-, parvo- and koniocellular layers of the dorsal lateral geniculate nucleus (LGN). Collectively, the selected regions allowed for interrogation of differences in gene expression between cortical areas and layers located distal or proximal to each other, and from regions that comprise specific functional types or streams. Representative pre- and post-cut images from each structure are shown in Figure S1
, and stereotaxic locations of sampled cortical regions in Table S1
. RNA was isolated from LMD samples, and 5ng total RNA per sample was amplified to generate sufficient labeled probe for use on Affymetrix rhesus macaque microarrays.
Experimental Paradigm for Profiling Macaque Cortical Regions and Layers
Multiple analytical methods were used independently to identify the most robust patterns of gene expression. Principle component analysis (PCA) can often illustrate the major organizational features of microarray data sets (Colantuoni et al., 2011
), and we initially applied it to the whole sample set comprising 225 cortical, hippocampal and thalamic samples across all 52,865 probes. A significant proportion of the variance was accounted for by the first three components (12.5%, 8.7% and 6.8% respectively; Figure S2
). As shown in , samples from major structures (cortex, hippocampus, and thalamus) cluster together, have highly distinct molecular signatures and appear well-segregated. Considering the cortical samples alone, the first three components accounted for a similar proportion of variance (13.6%, 8.5% and 6.6% respectively), and plotting samples by areal or laminar class revealed striking organization along two orthogonal axes reflecting the areal () and laminar () dimensions of the neocortex. Remarkably, the spatial relationships between neocortical samples are recapitulated by the transcriptional relationships between samples. Samples align in a rostral to caudal orientation by cortical area along the first principle component horizontally in , and appear tightly clustered in their native laminar order along the second principle component vertically in .
Major Organizational Features of the Rhesus Neocortical Transcriptome
To identify differentially expressed genes, three-way ANOVA of the cortical data set identified large numbers of probes that vary between cortical regions (6,170 at p<10-12
), layers (4,923) and individual animals (2,347; and Table S2
). Importantly, there was a high degree of overlap between the sets of genes varying by cortical region and layer, suggesting that a substantial proportion of the genes differentiating cortical areas vary within specific cortical layers. Gene set analysis of both areal and laminar genes showed enrichment for genes associated with axonal guidance signaling and ephrin receptor signaling, synaptic long term potentiation (LTP) and neuronal activities (Table S2
). Gene expression patterns associated with gender and individual animals were also identified by ANOVA (Figure S2
), and individual-associated differences were enriched with genes related to metabolism, mitochondria and antigen presentation (Table S2
). Gender-specific gene expression was observed both on sex and autosomal chromosomes (Fig S2
), and there was significant overlap (p<10-9
) between the individual-related genes identified here and gender-related genes identified in human brain (Kang et al., 2011
We next applied WGCNA to identify sets, or modules, of highly co-expressed genes by searching for genes with similar patterns of variation across samples as defined by high topological overlap (Zhang and Horvath, 2005
). Applied to the entire set of neocortical samples, WGCNA revealed a series of gene modules (named here as colors) related to different features of the data set (, also and ). Gene assignment to modules and gene ontology analysis for the whole cortex network are shown in Table S3
. The majority of these modules correlated with laminar and regional patterns as described below. Several modules were related to gender and individual differences, as previously observed in humans (Oldham et al., 2008
). In the lightyellow
module was strongly enriched in male versus female samples (upper panel), while the grey60
module was selectively lowest in samples originating from one particular animal. The top (hub) genes in the lightyellow
module were on the Y chromosome, including the putative RNA helicaseDDX3Y
and the 40S ribosomal protein RPS4Y1
Robust Transcriptional Signatures of Cortical Laminar Structure
Molecular Signatures of Cortical Regions
Cortical laminar variation
The most striking features were the robust molecular signatures associated with different cortical layers. As shown in , a wide variety of transcriptional patterns were associated with individual cortical layers or subsets of layers. For example, ANOVA of laminar expression in all cortical regions and clustering of these genes identified large gene sets enriched in specific subsets of (generally proximal) layers (). Notably, the majority of these laminar patterns are consistent across different cortical areas, reflecting conserved laminar and cellular architecture across the cortex. Gene set analysis suggests these layer-associated clusters are associated with neuronal function, including neuronal activity, LTP/LTD, calcium, glutamate and GABA signaling ( and Table S4
). Consistent with functional studies of superficial layer synaptic plasticity, genes and pathways involved in LTP and calcium signaling were most represented in L2 and L3. Pathways related to cholesterol metabolism were enriched in deeper layers, likely reflecting the greater proportion of oligodendrocytes closer to the underlying white matter. Similarly, many of the gene modules identified through WGCNA of all cortical samples were correlated with specific cortical layers (). By ANOVA-based clustering and WGCNA, proximal layers showed the strongest correlations, with superficial L2 and L3 highly correlated with one another, and the deeper L4-6 highly correlated as well (dendrograms in ).
Individual layers showed highly specific gene expression signatures. Layer-enriched expression patterns were identified by searching for genes with high correlation to layer-specific artificial template patterns (Lein et al., 2004
). shows cohorts of genes with remarkably layer-specific expression that was relatively constant across all cortical areas. These observations demonstrate the specificity of the laminar dissections with minimal interlaminar contamination, and also the constancy of laminar gene expression across the neocortex.
WGCNA gene modules derived from the whole cortex network also showed highly layer-enriched expression, demonstrating the robustness of our findings. For example, the black
module contains genes enriched in superficial L2 (hub genes plotted in , top row). While some layer-specific genes could be identified by targeted analyses, the dominant patterns were more complex, with most network modules being associated with combinations of layers, typically proximal to one another. For example, individual modules were enriched in L2-4 (salmon
), L3-5 (greenyellow
), L4-5 (royalblue
) and in a gradient increasing from L2 to L6 (red
). This tendency for coexpression between adjacent layers is also apparent in the heatmap representation of gene clusters in . Gene ontology (GO) analysis of these modules provides some insight into their functional relevance (Table S3
). The greenyellow
module was enriched for genes associated with axons and neuron projections, potentially related to long-range pyramidal projection neurons in L3 and L5. The red
module showed enrichment in genes associated with myelination, consistent with the presence of oligodendrocytes in deep layers, and this module was highly correlated with oligodendrocyte-associated gene networks in other studies (data not shown, (Oldham et al., 2008
Interestingly, the expanded L4 of V1 displayed a distinct signature from the rest of L4 (see top of middle box in ). To explore this further, we performed ANOVA and WGCNA selectively on samples from V1 ( and Tables S6 and S7
). A comparison between V1 ANOVA-derived laminar differential expression and membership in whole cortex WGCNA modules is in Table S8
. Similar to the whole cortex analysis, robust clusters and network modules were associated with individual cortical layers. As shown in the unsupervised hierarchical 2D clustering of ANOVA results in , individual samples from each layer cluster together, and neighboring cortical layers are most similar to one another. Interestingly, L4A clusters with more superficial layers, while L4B, L4Ca and L4Cb display a distinct transcriptional pattern, most easily seen by the dendrograms based on ANOVA and network analysis in .
To investigate whether layer specificity of gene expression may relate to selective patterns of connectivity, we examined the relationship between thalamocortical inputs and their targets in V1. L4Ca and L4Cb receive input selectively from magnocellular (M) and parvocellular (P) divisions of the LGN, respectively. Hypothesizing that there may be substantial shared gene expression patterns selective for specific pairs of connected neurons, we searched for genes that were differentially expressed between the thalamic inputs and between the cortical targets. 1,002 probes were differentially expressed between L4Ca and L4Cb (t-test, p<.01), and 825 probes between M and P. Surprisingly, these gene sets did not significantly overlap (13/1827, p=.08). Although the possibility certainly exists that specific ligand-receptor pairs are associated with this selective connectivity, it would appear that the specificity of these connections is not associated with specific large-scale correlated gene expression patterns.
To validate the specificity of the microarray findings and test hypotheses about laminar enrichment based on ANOVA and WGCNA, we examined a set of genes displaying layer-enriched patterns using in situ hybridization (ISH) in areas V1 and V2 (). Overall the laminar specificity of gene expression and variations between cortical areas predicted by microarrays were confirmed by cellular-level analysis, and illustrate the high information content of layer-specific expression profiling and gene specificity of the microarray probesets. For example, GPR83 is selectively expressed in L2 of all cortical areas, both by microarray and ISH analysis (). Laminar specificity was confirmed for RORB (L3-5; ), PDYN (L4-5; ), CUX2 (L2-4; ), and SV2C (L3-4 enriched; ). Specificity for deep cortical layers was prominent, as shown for PDE1A (L5-6; ), NR4A2 (L5-6; ), COL24A1 (L6; ) and RXFP1 (L5-6; ). Differences in laminar specificity were sometimes apparent between V1 and V2 (generally V1 versus all other areas); CUX2 was expressed in L2 through L4Cb in V1 but more limited to L2 and L3 in V2 (), and SV2C was highest in L4B in V1, but highest in L3 in area V2 ().
Layer-specific Gene Expression in Visual Cortical Areas
Cortical areal variation
Both ANOVA () and WGCNA analysis () identified gene clusters enriched in specific subsets of cortical regions. As illustrated in the dendrograms from both methods, the strongest relationships between cortical areas were based on areal proximity rather than functionally connectivity. For example, the caudal visual areas V1, V2 and MT showed highly correlated patterns of gene expression, while the functionally related but distal visual region TE had greater transcriptional similarity to its proximal neighbor A1 in temporal cortex. Strong relationships were observed for the adjacent primary motor and sensory cortices M1 and S1, and for the frontal DLPFC and OFC regions. Differentially expressed genes showed enrichment in specific subsets of (generally proximal) cortical areas (), generally related to neuronal development and function (axon guidance, neuronal activities, LTP/LTD; Table S9
). Areal expression also had a strong laminar signature, easily visualized by grouping these ANOVA-derived genes by cortical layer (Figure S3
Parallel relationships between cortical areas were observed by WGCNA demonstrating the robustness of these observations (), with individual gene modules showing enrichment in specific cortical regions (). Module eigengenes revealed additional patterning, including rostrocaudal gradients and laminar components to areal patterning. For example the tan module (, upper left) reflected a caudal low to rostral high patterning enriched in deep L5 and L6. Another gene module (purple, upper right) had an opposite gradient from high caudal to low rostral, in this case enriched in L3 and L4. Other modules were more area-specific: in V2, MT, DLPFC and OFC (blue) or lowest in V1, V2 and MT, with enrichment in L2 and L3 (pink).
Individual genes showed a wide range of areal patterns reflecting the modules described above, as well as patterns related to individual cortical areas or combinations of areas. Example gene patterns derived from the clustering analyses above, as well as analysis of the genes showing maximal cross-area fold changes, are shown in and . A large cohort of genes displayed rostrocaudal gradients. For example, MET, PVALB
were expressed most strongly in caudal V1 and decreased moving rostrally. Typically this gradient expression also had a laminar component. For example, MET
, which has been associated with autism (Campbell et al., 2006
) was expressed preferentially in L4 of V1, expanded to include L2-4 in temporal areas, and had no detectable expression in frontal cortices (), consistent with observations in fetal human brain (Mukamel et al., 2011
). Conversely, WFDC1
() was expressed most strongly in frontal cortical areas (e.g. DLPFC) and lowest in caudal V1, with expression primarily restricted to L2.
Diversity of Regional Gene Expression Patterns
V1-selective Macaque Patterns of Gene Expression and Cross-species Comparisons
Many genes showed area-selective expression (). NEFH, an intermediate filament heavy chain subunit generally expressed in long range projection neurons, was selectively enriched in L5 of M1 where the longest, spinally-projecting neurons, Betz cells, are located. Other genes were selectively enriched in specific regions of frontal cortex, including ACG (CALML4, IGFBP5, and LXN) or DLPFC and OFC (CD53). ISH showed selective enrichment of IGFBP5 in ACG compared to DLPFC and OFC (). Conversely, ISH analysis of CD53 () showed enrichment in dorsolateral, ventrolateral and ventromedial cortex relative to dorsomedial cortex and ACG.
Comparative analysis of laminar and areal expression patterns
The laminar and areal patterning observed in macaque was then compared to homologous structures in human and mouse. This phylogenetic comparison is shown in , made possible by the availability of mouse and human ISH data in the Allen Mouse Brain Atlas (http://mouse.brain-map.org
) and Allen Human Brain Atlas (http://human.brain-map.org
) generated using the same ISH technology platform as the macaque data. Some macaque data was also derived from the NIH Blueprint Non-Human Primate Atlas (www.blueprintnhpatlas.org
). Most genes in above showed laminar expression patterns that were highly conserved between macaque and human, with lesser conservation in mouse. Laminar patterns were generally conserved between primates and mice in visual cortex for CUX2, RORB, RXFP1
(left panels in ), although RXFP1
in mouse showed some regional differences between visual and somatosensory cortices not apparent in the macaque (). Other genes with highly laminar patterns showed major differences between species, indicated by blue arrowheads in . PDYN
() was enriched in excitatory neurons in L4 and L5 in rhesus and human, but in neurons with a broader laminar distribution that co-label with GAD1
indicating expression in GABAergic interneurons in mice (data not shown). The synaptic vesicle protein SV2C
was predominantly enriched in superficial L3 pyramidal neurons in most cortical areas in primates (e.g. V2), while it is fairly selective for deep L5 pyramidal neurons in mice. NR4A2
was expressed selectively in deep L5 and L6 in macaque and human, and in both L6 and L2/3 in mice ().
Cross-species Comparison of Laminar Gene Expression Patterns Between Adult Macaque, Human and Mouse Neocortex
While laminar distributions between V1 and V2 were highly conserved between macaque and human, differences were also noted. PDYN was expressed in L4B in macaque V1 in addition to the dominant L4Cb/5 expression (), while in human only L 4Cb/5 expression was observed. Similarly, faint L2 expression in addition to L5/6 enrichment was observed for PDE1A in macaque but not human (). In both cases, only the fainter expression in macaque was not observed in human, leaving the possibility that these differences relate less to true biological differences than to detection sensitivity on postmortem human tissues as compared to rapidly frozen rhesus macaque specimens.
The most robust patterns of areal specificity, both in terms of numbers of genes and their relative fold differences, were related to the highly specialized area V1 (). Both selective enrichment in V1 and selective lack of expression in V1 were observed, with a sharp boundary corresponding to the cytoarchitectural boundary observed by Nissl staining. This areal patterning was typically restricted to particular cortical layers as well. Some of this selective expression related to the expanded input L4 in V1. For example, ASAM, VAV3 and ESRRG were enriched primarily in L4 of V1 (). However, selective enrichment or decreased expression was seen in all cortical layers, including L2 and L3 (MEPE and RBP4; ), L5 (HTR2C; ), and L6 (CTGF, SYT6 and NPY2R; ).
The V1-selective patterning appeared to be highly conserved between macaque and human, while significant differences were observed between primates and mice (). For example, the enrichment of SYT6 () and NPY2R () in L6 of V1 relative to V2 was conserved between macaque and human, as was absence in L5 of V1 for HTR2C (). NPY2R expression showed a completely different pattern in mice, restricted to sparse (presumably GABAergic) neurons scattered across the cortex. Conversely, for both SYT6 and HTR2C, laminar restriction to L6 and L5 respectively was conserved in mice, but with no selective enrichment or lack of expression in V1. Thus these V1-specific gene expression differences correlate with primate-specific cytoarchitectural and functional specialization, rather than with the functional sensory modality subserved by visual cortex.