The regional variation of laminar thickness in the human isocortex is related to cortical hierarchy and interregional connectivity

The human isocortex consists of tangentially organized layers with unique cytoarchitectural properties. These layers show spatial variations in thickness and cytoarchitecture across the neocortex, which is thought to support function through enabling targeted corticocortical connections. Here, leveraging maps of the 6 cortical layers based on 3D human brain histology, we aimed to quantitatively characterize the systematic covariation of laminar structure in the cortex and its functional consequences. After correcting for the effect of cortical curvature, we identified a spatial pattern of changes in laminar thickness covariance from lateral frontal to posterior occipital regions, which differentiated the dominance of infra- versus supragranular layer thickness. Corresponding to the laminar regularities of cortical connections along cortical hierarchy, the infragranular-dominant pattern of laminar thickness was associated with higher hierarchical positions of regions, mapped based on resting-state effective connectivity in humans and tract-tracing of structural connections in macaques. Moreover, we show that regions with similar laminar thickness patterns have a higher likelihood of structural connections and strength of functional connections. In sum, here we characterize the organization of laminar thickness in the human isocortex and its association with cortico-cortical connectivity, illustrating how laminar organization may provide a foundational principle of cortical function.

based on their expression intensity relative to background noise [4], such that probes with intensity less than the background in ≥ 50% or 25% of samples across donors were discarded.The threshold of ≥ 50% was used for correlated gene expression and developmental gene enrichment analyses, and ≥ 25% was used for creating transcriptomics maps of neuronal subtypes, as it led to removal of fewer genes associated with these subtypes.When multiple probes indexed the expression of the same gene, we selected and used the probe with the most consistent pattern of regional variation across donors (i.e., differential stability) [5], calculated with: where  is Spearman's rank correlation of the expression of a single probe, p, across regions in two donors  and , and N is the total number of donors.Here, regions correspond to the structural designations provided in the ontology from the AHBA.The MNI coordinates of tissue samples were updated to those generated via non-linear registration using the Advanced Normalization Tools (ANTs; https://github.com/chrisfilo/alleninf).Samples were assigned to brain regions by minimizing the Euclidean distance between the MNI coordinates of each sample and the nearest surface vertex.Samples where the Euclidean distance to the nearest vertex was more than 2 standard deviations above the mean distance for all samples belonging to that donor were excluded.To reduce the potential for misassignment, sample-to-region matching was constrained by hemisphere and gross structural divisions (i.e., cerebral cortex, subcortex/brainstem, and cerebellum, such that e.g., a sample in the left cerebral cortex could only be assigned to an atlas parcel in the left cerebral cortex) [3].If a brain region was not assigned a sample from any donor based on the above procedure, the tissue sample closest to the centroid of that parcel was identified independently for each donor.The average of these samples was taken across donors, weighted by the distance between the parcel centroid and the sample, to obtain an estimate of the parcellated expression values for the missing region.All tissue samples not assigned to a brain region in the provided atlas were discarded.
Inter-subject variation was addressed by normalizing tissue sample expression values across genes using a robust sigmoid function [6]: where <> is the median and IQR is the normalized interquartile range of the expression of a single tissue sample across genes.Normalized expression values were then rescaled to the unit interval: /01234 =  *+,-− min( *+,-) max( *+,-) − min( *+,-) Gene expression values were then normalized across tissue samples using an identical procedure.Samples assigned to the same brain region were averaged separately for each donor and then across donors, yielding a regional expression matrix.We then excluded the right hemisphere regions due to the sparsity of samples and large number of regions with no expression data.

Transcriptomics maps of excitatory and inhibitory neuronal subtypes
The RNA sequencing of isolated cortical neuronal nuclei and their data-driven clustering has previously identified 16 neuronal subtypes, consisting of eight excitatory and eight inhibitory subtypes [7].We used the list of subtype-specific genes and estimated the spatial distribution of each subtype by taking the weighted average mRNA expression of genes associated with the subtype, weighting each gene by the fraction of its positive values using exon-only derived transcripts per million [8].Of note, few of the genes associated with each neuronal subtype had no eligible expression data in the processed AHBA dataset and were omitted from the calculation of subtypespecific gene expression (Table I).
We observed significant correlations of the principal axis of laminar thickness covariance (LTC G1), asymmetry-based hierarchy and laminar-based hierarchy maps with the spatial variation of excitatory and inhibitory neuronal subtypes based on gene expression (Fig I) [1,2,7].This finding was in line with previous work on the association of cortical hierarchy and microstructure with cortical microcircuits [8][9][10].

Correlated gene expression matrix
The correlated gene expression matrix was computed by correlating the regional gene expression patterns between the regions and across all genes, followed by its Z-transformation.Further supporting shared genetic effects between regions with similar laminar thickness, we observed a correlation of r = 0.24 (p spin = 0.010) between LTC and the correlated gene expression matrix.This matrix shows the correlation of gene expression between regions and across all genes, based on transcriptomics data obtained from the Allen Human Brain Atlas (AHBA; Fig II) [1,2].

Fig II. Correlated gene expression in association with laminar thickness covariance.
The gene expression data from Allen Human Brain Atlas (left) was used to create correlated gene expression (CGE) matrix, showing inter-regional similarity of gene expression profiles (middle), which was compared to the laminar thickness covariance matrix (LTC; right).The data and code needed to generate this figure can be found in https://zenodo.org/record/8410965.

Partial least squares regression of LTC G1 with gene expression maps and developmental enrichment analysis
We used partial least squares regression (PLS) to identify genes that show high alignment in their regional expression variability with the LTC G1 within the left hemisphere.The scikit-learn package (https://scikit-learn.org/stable/)[11] was used to perform single-component PLS between the parcellated map of LTC G1 on one side, and the regional expression of all genes on the other side.We then ranked the genes based on the absolute value of their weights and selected the top 500 genes, which were subsequently divided into genes with positive and negative weight.Next, we fed the list of top genes with positive and negative alignment to LTC G1 to the online cell-type specific expression analysis (CSEA) developmental expression tool (http://genetics.wustl.edu/jdlab/cseatool-2/)[12], which involves comparison of the genes against developmental expression profiles from the BrainSpan Atlas of the Developing Human Brain (http://www.brainspan.org), and for each developmental stage and brain structure reported the inverse log of Fisher's Exact p-values corrected using Benjamini-Hochberg approach and based on the specificity index (pSI) threshold of 0.05.
We identified distinct sets of genes expressed preferentially at the opposite ends of the LTC G1 in adulthood based on the AHBA data, and investigated their developmental enrichment using the BrainSpan dataset [12].We observed distinct spatiotemporal patterns for the expression of each set of genes associated with the rostral versus caudal ends of the LTC G1.Specifically, within the cerebral cortex, genes preferentially expressed at the caudal end of LTC G1 are mainly expressed in mid fetal and post-pubertal stages, while genes with higher expression at the rostral end of LTC G1 were expressed in late fetal and pre-pubertal stages of development (Fig III ).

Fig III.
Developmental enrichment of genes expressed in alignment with the principal axis of laminar thickness covariance.The average expression of top genes with regional expressions aligned to the principal axis of laminar thickness covariance (LTC G1), over-expressed rostrally (blue) or caudally (yellow) in the left hemisphere, and their spatiotemporal developmental enrichment.The data and code needed to generate this figure can be found in https://zenodo.org/record/8410965.The clip-arts are based on Ref. [13] published under CC BY 4.0 (https://creativecommons.org/licenses/by/4.0/).

Fig I .
Fig I. Microcircuitry in association with cortical hierarchy and the principal axis of laminar thickness covariance.Spatial distribution of excitatory (a) and inhibitory (c) neuronal subtypes based on gene expression (left hemisphere), and their correlation with the principal axis of laminar thickness covariance (LTC G1) as well as asymmetry-based and laminar-based hierarchy (b, d).Bar length shows the correlation coefficient and its color represents the level of statistical significance from white (p variogram, FDR > 0.05) to black (p variogram, FDR < 0.001).The data and code needed to generate this figure can be found in https://zenodo.org/record/8410965.