Spatial transcriptomics reveals human cortical layer and area specification

Human sample preparation and acquisition

Research performed on samples of human origin was conducted with approval of the institutional review board of Boston Children’s Hospital. De-identified human fetal tissue samples were collected after written informed consent was acquired from two sources. (1) Samples were received from the Beth Israel Deaconess Medical Center Department of Pathology, with a maximum post-mortem interval of 12 h. Tissue was transported in Hibernate-E medium on ice to the laboratory for research processing. Cortical tissue was then dissected into coronal pieces of 1–3 cm² cross-section area and 0.5–1 cm thickness. The dissection was performed under the supervision of a neuropathologist to provide annotation of cortical areas based on anatomical location and features. Tissue pieces were then directly frozen in liquid nitrogen and stored at −80 °C. (2) Banked fresh-frozen tissues were obtained from the University of Maryland Brain and Tissue Bank through the NIH NeuroBioBank. Only samples with postnatal intervals of <12 h were used. Cortical area annotation was provided by NeuroBioBank. Tissue was shipped overnight in dry ice and stored at −80 °C. For both sources, only samples with no neurological anomalies were analysed. Samples were screened for RNA quality by collecting 50-µm-thick cryosections, isolating total RNA and measuring the RNA integrity number using an Agilent 4200 TapeStation system, and RNA samples with integrity numbers of ≥6.5 were used in the study. Details of the samples used are summarized in Supplementary Table 2.

Annotation of cortical areas

The cortical areas of tissue sections used in MERFISH analyses were annotated on the basis of the reference atlas for human fetal brain developed by the Allen Institute for Brain Science (https://atlas.brain-map.org)³⁰. For GW20 and GW22 samples, dissection of the tissue was performed under the supervision of a neuropathologist to distinguish the frontal lobe, the parietal lobe, the occipital lobe and the temporal lobe. The relative AP, dorsal–ventral and medial–lateral location of each piece of dissected tissue was recorded by notes and photographs. For GW15, GW18 and GW21 samples obtained from the NIH NeuroBioBank, annotation of cortical lobes was provided by NeuroBioBank staff. For more specific annotation of cortical areas, the GW15 samples were compared with the GW15 reference atlas (https://atlas.brain-map.org/atlas?atlas=138322603), and the GW18, GW20, GW21 and GW22 samples were compared with the GW21 reference atlas (https://atlas.brain-map.org/atlas?atlas=3). The relative AP locations of tissue sections were matched with the reference atlas, and cytoarchitectural features such as early sulci were used as reference landmarks. The transition between the PMC and the M1 lacks clear definition in mid-gestation; therefore, some tissues analysed were annotated as PMC/M1. By contrast, the calcarine sulcus was morphologically well defined in GW15–GW22 samples, which enabled the identification of the V1 and the V2. For GW34 and adult samples from NeuroBioBank, the specific Brodmann’s area corresponding to the postnatal human cortex was determined by NeuroBioBank staff and was provided when tissues were requested. Definition of the specific Brodmann’s area is possible only for GW34 and adult samples because by this age, the major sulci and gyri structures have formed, which is not the case for GW15–GW22 periods. Annotation of the cortical area for all samples is summarized in Supplementary Table 2.

Immunohistochemistry and microscopy

Cryosections (10 µm thick) of human fetal tissue were placed on Superfrost slides (Fisher) for immunohistochemistry. Preparation of fetal macaque tissue samples are described separately, but immunostaining for NPY followed the same procedure. The samples were first fixed in 4% paraformaldehyde (PFA) for 15 min at room temperature, followed by permeabilization with 0.5% Triton-X in PBS for 1 h and blocked with blocking solution of 10% donkey serum in PBS and 0.05% Triton-X (PBST) for 30 min. Primary antibodies diluted 1:500 in blocking solution were applied to the sections overnight at 4 °C. Primary antibodies used were rat anti-CTIP2 (Abcam, ab18465), mouse anti-SATB2 (Abcam, ab9244), rabbit anti-TBR1 (Abcam, ab31940), rabbit anti-NPY (Abcam, ab30914) and goat anti-synaptophysin (R&D Systems, AF5555-SP). After washing with PBST for a minimum of 5 times, secondary antibodies and DAPI (1:2,000) diluted in blocking solution were applied to the sections for 1–4 h at room temperature or overnight at 4 °C. Secondary antibodies used were AlexaFluor 488-conjugated, 555-conjugated or 647-conjugated donkey antibodies (Invitrogen) used at 1:500 dilution. Finally, sections were washed with PBST for a minimum of 5 times before mounting with Vectashield Vibrance Antifade mounting medium (Vector Labs, H-1700-10). Images were captured using a Zeiss LSM 980 confocal microscope. The z stack function was used to image at 10 µm thickness, and the tile-stitching function was used. Sample images were prepared in ImageJ software (Fiji 2.15.1) and Adobe Photoshop (m.2506).

MERFISH gene-panel selection and probe construction

We initially designed a panel of 300 genes (Supplementary Table 3) that was used for most of the experiments. The panel contained 39 canonical markers for major cell types in the human fetal cortex, including markers for ENs and INs, IPCs and neural progenitor cells; 32 were validated cortical-layer markers in the adult human cerebral cortex from a previous study³⁸. Thirty genes associated with autism spectrum disorders were selected from the SFARI Gene Portal (https://gene.sfari.org/database/human-gene/)⁶⁴, and one putative regulator long noncoding RNA for each SFARI gene was included (30 in total). The remaining genes on the panel were obtained from the top enriched cluster markers in published single-cell RNA SMART-seq data of mid-gestation human fetal cerebral cortex⁷. Single-cell clustering and marker identification based on differential expression were performed in the original study⁷. The top 20 marker genes based on the fold change for EN clusters and the top 5–10 marker genes for other clusters were manually curated. The genes that did not overlap with previous selected categories were included in the gene panel, which resulted in a total of 300 genes. All genes in the panel were included for subsequent analyses.

We also designed a new 960 gene panel that contained all 300 genes from the first gene panel and expanded. This panel was only used for the UMB1031-O1, UMB1759-O1 and UMB1932-O1 experiments. The following criteria were used for including additional genes: (1) genes clinically associated with malformation of the cerebral cortex; (2) genes exhibiting high expression variance from our snRNA-seq data of GW20 human fetal cortex; (3) DEGs between the V1 and V2 at GW20 identified by our Visium data; (4) genes involved in postnatal V1 specification in mice, taken from a published scRNA-seq analysis⁵¹; (5) new markers for interneuron subtypes identified in developing macaque brain, taken from a scRNA-seq analysis⁶⁵; (6) additional classical interneuron-subtype markers; (7) marker genes for cell types found in the human fetal GE, taken from a scRNA-seq analysis⁶⁶; and (8) marker genes for cell types found in the human fetal thalamus, taken from a scRNA-seq analysis⁶⁷.

Merscope-encoding probes for the 300 and 960 genes were constructed by Vizgen using a commercial pipeline. Each of the 300 and 960 genes was assigned a unique binary barcode drawn from a 22-bit, Hamming-Distance-4, Hamming-Weight-4 encoding scheme⁹ (Supplementary Table 4). Fifteen extra barcodes (for the 300 gene panel) and 40 extra barcodes (for the 960 gene panel) were included as blank barcodes, which were not assigned to any genes to provide a measure of the false-positive rate in MERFISH, as previously described⁹.

MERFISH imaging

MERFISH experiments was performed using a Vizgen Merscope system. Sample preparation was performed according to manufacturer’s instructions (‘Merscope Fresh and Fixed Frozen Tissue Sample Preparation User Guide’, number 91600002). In brief, fresh-frozen tissues were embedded in OCT and sectioned into 10-µm-thick sections using a cryostat (Leica) and adhered to Merscope slides (Vizgen, 1050001) placed in a 6-cm Petri dish. For GW15–GW22 samples, slides were kept inside the cryostat maintained at −15 °C for 30 min to allow the section to dry and firmly adhere to the glass slides before fixation. For GW34 samples, slides were kept inside the cryostat for 5 min and then transferred to room temperature for 15 min before fixation. Slides were then fixed in 4% PFA and permeabilized by 70% ethanol overnight, with Parafilm sealing the Petri dish to prevent evaporation. Slides were then treated in a Merscope photobleacher instrument for autofluorescence quenching for 3 h. Slides were stored overnight or up to 1 week before proceeding to the next step. The encoding probe mix was added directly on top of the tissue section for hybridization at 37 °C for 36 h in a humidified incubator. After probe hybridization, sections were fixed again using formamide and embedded in gel. After gel embedding, tissue samples were cleared using a clearing mix solution supplemented with proteinase K for 24–48 h at 37 °C until no visible tissue was evident in the gel. After clearing, sections were stained for DAPI and PolyT and fixed with formamide before imaging. No additional cell boundary staining was used. Reagents used for these steps were included in Merscope sample preparation kit (Vizgen 10400012).

The MERFISH imaging process was performed according to the Merscope Instrument Preparation guide (number 91500001). In brief, an imaging kit was thawed in a 37 °C water bath for 45 min, activated and loaded into the Merscope instrument. The flow chamber was then assembled, fluidics were primed and the flow chamber filled with liquid. A low-resolution image for the DAPI and PolyT stains was taken at ×10 magnification, and a region of interest (ROI) was manually drawn in the Merscope software, followed by automated image acquisition and fluidic control in the Merscope instrument. For each section, a ROI of up to 1 cm² area was imaged as 2,000–2,500 tiles at ×40 magnification, and images were collected in the 750 nm, 650 nm and 560 nm channels for the readout probes and in the 488 nm and 405 nm channels for PolyT and DAPI staining, respectively. Seven z stacks were captured over a thickness of 10 µm. After imaging, image processing and transcript decoding were performed using Merscope proprietary software. The transcript matrix with spatial coordinates and the stitched tiled DAPI images acquired were transferred for subsequent processing and single-nucleus segmentation. Vizgen Merscope Visualizer (v.2.1.2595.1) software was used to visualize each experiment to determine whether the MERFISH run was successful before systematic analysis was performed.

Nucleus segmentation and MERFISH data processing

Automated segmentation was performed on the DAPI channel using a custom CellPose model^23,24. The model was initialized with the CellPose (2.0) ‘cyto’ weights, then trained for 300 epochs with a learning rate of 0.1 and a weight decay of 1 × 10^–4 using 145 manually segmented images for training and 3 for testing. All images, before training and running segmentation, were filtered using a difference of Gaussians filter, which consisted of a positive Gaussian having zero standard deviation (all weight at zero) and a negative Gaussian with a standard deviation of 20 pixels. When running the CellPose segmentation, for testing and in the Vizgen pipeline, the cell diameter was set to 55 pixels (average from training data), the flow threshold was set to 0.5, the cell probability threshold was set to –3 and the minimum mask size was 500 pixels.

The Vizgen Merscope output consisted of 7 planes evenly spaced across 10 µm. Given the wide point-spread function and the small spacing, we used the maximum projection image over the first six planes as input to the CellPose segmentation algorithm, using the CellPose parameters and filtering as specified above. Because nucleus staining (DAPI) was used, we dilated the masks by 10 pixels to approximate the cytoplasmic area of the cells. These processing steps were added to the Vizgen Post-processing Tool (v.1.2.0); the modified code is available from GitHub (https://github.com/carsen-stringer/vizgen-postprocessing). The Vizgen processing pipeline uses the CellPose masks to define ROIs and then assigns the RNA transcripts to each ROI to return a cell-by-gene matrix. In Extended Data Fig. 1b, we quantified the segmentation accuracy using the average precision metric, which is the number of true positives divided by the total number of true positives, false negatives and false positives²⁴. The true positives were defined as segmented ROIs that matched ground-truth ROIs at or above a defined intersection-over-union threshold. The false negatives were the ground-truth ROIs that were missed, and the false positives were the predicted ROIs that did not match any ground-truth ROIs.

To compare the transcriptomic profiles between MERFISH cells segmented with and without dilation, we performed a test (Extended Data Fig. 1d,e). The same experiment (FB080-O1c) was segmented twice using our CellPose and Vizgen Post-processing Tool pipeline, with the dilation function disabled or enabled, respectively. The cells segmented with and without dilation were combined, preprocessed and clustered into 30 unsupervised clusters using sklearn.cluster.KMeans(). The resulting clusters were visualized using UMAP in Extended Data Fig. 1d. To further compare the source distribution of cells in each cluster, we calculated and visualized the proportions of dilated and non-dilated cells in each cluster as a bar plot in Extended Data Fig. 1d. We also compared the distribution of transcript count in each cluster (Extended Data Fig. 1e).

MERFISH data quality control and integrated hierarchical clustering

After performing nucleus segmentation, for each experiment, we filtered out all cells with a total transcript count below the tenth percentile specific to that experiment. We then normalized each cell by its total transcript count using scanpy.pp.normlize_total(). Transcript counts were then log-transformed and z score-normalized using scanpy.pp.log1p() and scanpy.pp.scale(), respectively. Once preprocessed, we integrated the gene-expression data from all samples. Subsequently, hierarchical clustering was performed on the integrated dataset to identify groups of cells with similar gene-expression profiles. Initially, all cells were grouped into eight H1 clusters at the first hierarchy (cell class). After this step, each H1 cluster was subclustered into 5 preliminary H2 clusters, which resulted in 40 preliminary H2 clusters (cell type). Each H2 cluster was further subdivided into 5 preliminary H3 clusters to form a total of 200 preliminary H3 clusters (cell subtypes). The clustering results at each hierarchy were obtained using sklearn.cluster.KMeans(). Marker genes were then identified using scanpy.tl.rank_gene_groups() and scanpy.get.rank_genes_group_df(). For a given H1 cluster, genes with a log₂[fold change] of at least 0.25 compared with the cells in all other H1 clusters were denoted as marker genes for that cluster. Marker genes for each H2 or H3 cluster were identified using the log₂[fold change] between that cluster and the other four clusters belonging to the same H1 or H2 group, respectively.

Cluster annotations were conducted sequentially from H1 to H3 on the basis of the marker gene list and the spatial distribution of each cluster. The eight H1 cell classes were annotated on the basis of expression of canonical markers. Some preliminary H2 and H3 clusters were manually merged. Because ECs were not of particular interest in our study and our gene panel has limited relevant genes for ECs, all EC clusters were merged into one cluster. Similarly, astrocytes and oligodendrocyte (OPC) clusters were merged into three H2 clusters, and subsequently five H3 clusters. H3 clusters for IPCs derived from each H2 IPC cluster were highly similar in spatial distribution and gene expression; therefore, these were also merged. Finally, we noted that a few clusters exhibited aberrant spatial distribution that reflected technical artefacts during the imaging process. These clusters, easily recognizable because of their exclusive localization at the edge of certain tissue sections or surrounding bubbles in the tissue section, were the accidental result of tissue hydrogel detaching from the surface during imaging. We removed cells from these artefact clusters from all subsequent analyses.

After transcript count filtering, artefact cluster removal and supervised merging, we analysed 15,927,370 cells from the initial set of experiments, which were clustered into 8 H1 clusters, 33 H2 clusters and 114 H3 clusters (Supplementary Table 5). H2 cell types were annotated to reflect their spatial distribution. For EN-IT, EN-ET and EN-Mig clusters, their approximated layer enrichment at GW34 were denoted. For RG and IPCs, their enriched laminar structures (such as the VZ, the SVZ and the IZ) were denoted in cluster annotation. In addition, if a cluster was predominantly present only in GW15 samples, it was denoted with ‘early’, whereas a cluster predominantly present only in GW34 was denoted with ‘late’. H3 EN clusters were similarly annotated to reflect their layer enrichment and their areal distribution enrichment (A for anterior, P for posterior, T for temporal lobe, and so on). In subsequent set of experiments that were separately analysed, 2,635,057 cells were analysed using the same quality-control procedures.

Correlation between technical and biological replicates

To verify the reproducibility of our MERFISH experiments, we computed the mean expression for each gene and visualized the correlation between the mean expression levels from two replicates as scatterplots (Extended Data Fig. 1i,j). Pearson correlation coefficients were calculated to assess the agreement. Horizontal and vertical dashed lines marked the mean expressions of blank genes for each replicate. Technical replicates analysed included UMB1117-F1a and UMB1117-F1b, UMB1117-F2a and UMB1117-F2b, FB080-F2a and FB080-F2b, FB080-P1a and FB080-P1b, FB080-O1a and FB080-O1b, FB080-O1b and FB080-O1c, and FB080-O1a and FB080-O1c. Biological replicates included FB080-F1 and FB121-F1, FB080-F2a and FB121-F2, FB080-P1a and FB121-P1, FB080-P2 and FB121-P2, FB080-T1 and FB121-T1, and FB080-O1a and FB121-O1. The mean and standard deviation of correlation coefficients were computed and visualized in Extended Data Fig. 1k.

Cell-label-transfer analysis

Reference-based cell-type annotation was conducted using SingleR (v.1.8.1)⁶⁸ to transfer cell-type labels from the reference data to the testing data with the following parameters: “genes”: “de”; “de.method”: “classic”; “fine.tune”: TRUE; “prune”: TRUE. The five separate analyses are summarized below. (1) For Fig. 1f and Extended Data Fig. 3d, previous fetal human cortex scRNA-seq data⁷ were used as a reference to transfer labels to our MERFISH data. In the reference, we excluded ‘MGE-’ (cells from dorsal cortical specimens) and ‘unknown’ clusters and aggregated all ‘nIN’ subclusters into one and denoted it as ‘IN’. (2) For Extended Data Fig. 3b, previous scRNA-seq data⁸ were used as a reference to transfer cell-type labels to the MERFISH data. (3) For Extended Data Fig. 3e, previous scRNA-seq data⁶ were used to transfer cell-type labels, whereby ‘CombinedCluster-Final’ in the metadata was used to define reference cell types. Because the reference dataset from ref. ⁶ covered more gestational ages and cortical areas and non-cortical structures, it contains cell types for which counterparts are not expected to be found in our MERFISH dataset. Thus, we added a filter to hide rows with maximal correspondence fractions of <0.05 to save space on the graph. (4) For Extended Data Fig. 6d, cells from excitatory neuron (Exc) clusters in a published adult human cortex scRNA-seq dataset³⁵ were extracted and served as a reference to transfer cell-type labels to EN-ET and EN-IT cells in our MERFISH data. For simplicity, the 56 original EN clusters in the adult human cortex dataset were merged into 14 groups on the basis of their shared layer and molecular marker identities. (5) For Extended Data Fig. 18d, cells from the EN clusters in our snRNA-seq data were used as a reference to transfer labels to EN cells or spots in our MERFISH and Visium data separately. The pruned labels returned from SingleR were used to annotate the target cells. Correspondences between raw and transferred cell types in the target data were visualized using heatmaps.

Quantitative spatial analysis of MERFISH data

For the majority of tissues analysed by MERFISH, one to three fan-shaped regions were manually drawn for quantitative spatial analysis, referred to as ‘annotated area A, B or C’ in Supplementary Table 2. The fan-shaped regions were selected on the basis of relative geometrical uniformity, avoiding anatomical structures such as sulci, which may complicate location quantification. Areas with tissue section tears and bubbles were also avoided. For each tissue section analysed by MERFISH, spatial graphs of H1 cell-class distribution were generated. On the basis of these graphs, each fan-shaped region was made of hand-traced vector lines marking the apical surface and basal surface and were connected by two straight lines using the pen tool in Adobe Photoshop software. The basal borders were defined by the pial surface of the cortex. For GW15 and GW20 samples, the apical borders were drawn at the ventricular surface. For GW22 samples, the apical borders were drawn in the oSVZ as the VZ could not fit in the imaged area owing to the larger tissue size. For GW34 samples, the apical surfaces were drawn at the approximated transition between layer 6b and the white matter.

In each fan-shaped region, laminar structures of the developing cerebral cortex were further manually defined with respect to the distribution and morphology of cell classes. The CP was defined as the condensed layer of excitatory neurons (EN-ITs and EN-ETs) near the pial surface. The thin layer above the CP to the pial surface was defined as the MZ. The SP was defined as the layer containing only EN-ETs, below the CP, and with lower cell density than the CP. The IZ was defined as the relatively cell-sparse layer predominantly composed of EN-Mig cells below the SP. The oSVZ was defined as the region composed of a mixture of EN-Mig cells, IPCs and RG below the IZ. The VZ was defined as the condensed layer of RG showing clear ventricular morphology at the apical surface. The thin layer above the VZ composed of a high density of IPCs was defined as the iSVZ. The boundaries between each neighbouring laminar structure were manually traced using the pen tool in Photoshop, which established an annotated vector mask for each fan-shaped region.

The vector mask was then imported into the computational analysis object by overlapping the spatial coordinates of cells with the vectors. Therefore, each cell in the fan-shaped region was also assigned with its laminar structure (VZ, iSVZ, oSVZ, IZ, SP, CP or MZ). For each annotated region, three images were exported from Photoshop: (1) a mask containing the boundary of the fan-shaped region along with the lines delineating the boundaries of the annotated structures; (2) a mask containing only the boundary of the fan-shaped region; and (3) a mask containing the boundary of the fan-shaped region, where the lines marking the apical surface and basal surface were replaced by straight lines. The initial objective was to determine the coordinates of the apical and basal surfaces. Using image (3), we applied cv2.goodFeaturesToTrack() to identify the four corners of the fan shape. Using the orientation of the annotated region (up, down, left and right), corner coordinates were assigned to their respective ends of the apical and basal surfaces. We calculated the equations of the lines connecting the corners of the apical and basal surfaces, and the point at which the lines intersected (O) was determined. To identify all points comprising the apical and basal surfaces, we first set all pixels equal to 0 in the mask of image (1) that were in a 3-pixel radius of each of the corners. We then used cv2.connectedComponents() to categorize each of the four lines in the mask into different components. The component identifiers for the apical and basal surfaces were then determined on the basis of the previously calculated corner locations. To extract the cells in a specific layer, we first subtracted image (2) from image (1) to eliminate the boundaries of the fan-shaped region, including the apical and basal surfaces. With only the layer-defining boundaries remaining, we applied cv2.connectedComponents() to assign each boundary to a different component. The component identifiers corresponding to each of the layer boundaries were determined using the orientation of the fan-shaped region. To calculate the mask boundary for a given layer, we used cv2.line() to draw lines between the corresponding end points on the top and bottom boundaries for that layer. We then created the mask by applying cv2.binary_fill_holes() to the boundary. Using the coordinates of the cells located in the annotated region, we then extracted all cells with coordinates contained in the layer-specific mask.

For each cell in the annotated fan-shaped region, its RH was calculated to represent its laminar location from the apical to basal surfaces. The RH calculation normalizes for the tilting of the fan-shape created by the asymmetrical morphology of the human cerebral cortex. Similarly, for a cell in the CP of an annotated fan-shaped region, its CD was calculated to quantify its relative laminar location in the CP. To calculate the RH for a given cell, we used cv2.line() to draw a line from point O to the cell’s location at point C. Using the mask containing this line, we identified intersection points with the apical and basal surfaces (designated as C₃ and C₁, respectively). The Euclidean distances between C₃ and C₁, and between C and C₁ were then computed. The RH of the cell was then calculated as the ratio of these distances. To calculate the CD, we first extracted all cells located in the CP of the fan-shaped region. Then, for each cell located in the CP, we used cv2.line() to draw a line connecting O and the location of the cell, C. Using the mask containing this line, we identified the points where the line intersected with the top and bottom boundaries of the cortical plate (denoted by C₃ and C₂, respectively). We then calculated the Euclidean distances between C₃ and C, and between C and C₂. The ratio of these distances defined the CD of the cell. For a given annotated region, the RH violin plot was constructed using the H2 cell-type annotations for all cells in that region. The CD violin plot was created using the H3 cell-type annotations for all EN-IT, EN-ET and EN-Mig cells in the CP of that region. For any clusters with fewer than 50 cells in that region, we represented the RH or CD distribution with individual cell dots instead of a violin plot.

Cortical layers in the CP of an annotated region were computationally defined by the distribution of CD values for specific H3 EN subtypes. Although the laminar location of most H3 EN subtypes were consistent across gestational ages, some clusters showed varying abundance with gestational age, which prevented us from using the exact same set of H3 EN subtypes for the definition of cortical layers in samples from all gestational weeks analysed. For GW15, EN-IT-L4-1, EN-ET-L5-1 and EN-IT-L6-1 were used to defined layers 3/4, layer 5 and layer 6, respectively. For GW20, EN-IT-L2/3-A1, EN-IT-L4-1, EN-ET-L5-1 and EN-IT-L6-1 were used to defined layers 2/3, layer 4, layer 5 and layer 6 respectively. For GW22, EN-L2-1, EN-IT-L3-A, EN-IT-L4-1, EN-ET-L5-1 and EN-IT-L6-1 were used to define layer 2, layer 3, layer 4, layer 5 and layer 6, respectively. For GW34, EN-L2-4, EN-IT-L3-late, EN-IT-L4-late, EN-ET-L5-1 and EN-IT-L6-late were used to define layer 2, layer 3, layer 4, layer 5 and layer 6, respectively. The borders between adjacent layers, for example, between layer 5 and layer 6, were calculated as the mean CD value between the lower 25% of layer 5 cells and the upper 75% of layer 6 cells. Therefore, the CD value for each cell in the CP was used to assign a cortical layer identity for that cell.

Spatiotemporal expression patterns of genes

For each gestational week and cortical area, we isolated all cells in each annotated layer using their coordinates and the corresponding layer-defining masks. Cells in the CP were assigned to specific cortical layers on the basis of their CD values. Subsequently, for each gene, we calculated the mean expression across all cells in each annotated layer for a given gestational age. We then standardize the expression for that gene to the [0, 1] range across all layers and gestational weeks. The results are visualized in the summary expression heatmaps (Supplementary Data). To construct z score expression spatial graphs for a given cortical area, the normalized expression values for all cells in the area were extracted, and gene-expression values were z score-normalized using scanpy.pp.scale(). Gene-expression heatmaps were then plotted using scanpy.pl.embedding().

Cell number proportion and enrichment quantifications

The two panels at the bottom of Extended Data Fig. 9a illustrate (1) the proportion of cells in a H3 cell type for each cortical area (the PFC, the PMC/M1, the Par, the Temp and the Occi) and (2) the cell number of a H3 cell type for each gestational week (GW15, GW20, GW22 and GW34). For the cortical area proportion plot in (1), the proportion (r) of cells for each cortical area (CA) was calculated as the total number of cells for each CA for the H3 cell type divided by the total number of cells for that H3 cell type across all the CAs: \(r_\rmCA,\rmH3=\frac\,N_\rmCA,\rmH3N_\rmH3\) where N_CA,H3 is the number of cells for a CA at H3 and N_H3 is the total number of cells at H3. For the number of cells plot in (2), the number of cells for each gestational week at each H3 subcluster is shown.

The results of spatiotemporal enrichment analysis of EN subtypes (EN-ETs, EN-ITs and EN-Mig cells) in Extended Data Fig. 9a and RG, IPC and IN subtypes in Extended Data Fig. 9b were done using the normalization process described below. The total number of cells for each sample–region (SR) pair at H1 and H3 were counted separately for each gestational week. The number of cells of SR pairs at each H3 cell type was then divided by the number of cells of SR pairs at the corresponding H1 cell type. \(r_\rmSR,\rmH3=\frac\,N_\rmSR,\rmH3N_\rmSR,\rmH1\) where N_SR,H1 is the number of cells of a SR pair at H1 and N_SR,H3 is the number of cells of a SR pair at H3. The ratio above was then divided by the maximum number of cells of SR pairs at H3. \(\rmNormalized(r_\rmSR,\rmH3)=\fracr_\rmSR,\rmH3R_\rmH3R_\rmH3\) is the maximum of r_SR,H3 among all the SR pairs at that H3 cell type.

For Extended Data Fig. 2c, the proportions of all H2 cell types in each SR pair was calculated. The number of cells of SR pairs at each H2 cell type was divided by the total number of cells in that SR. \(r_\rmH2,\rmSR=\frac\,N_\rmH2,\rmSRN_\rmSR\), where N_H2,SR is the number of cells of a SR pair at H2 and N_SR is the total number of cells of an SR pair. Figures 4f and 5b were calculated in a similar manner.

For UMAP of EN subtype compositions (Extended Data Fig. 9c), all EN-ET and EN-IT cells from each annotated cortical region were extracted. The relative proportion of the EN-ET and EN-IT H3 subtypes was calculated for each cortical area. We computed the UMAP of the matrix of relative EN H3 cell-type proportions, providing a two-dimensional representation for each cortical area based on the EN subtype composition.

Identification of DEGs in anterior and posterior clusters

In Extended Data Fig. 10a, five pairs of anteriorly enriched and posteriorly enriched neuronal subtype H3 clusters were compared on the basis of their DEGs: EN-IT-L2/3-A2 versus EN-IT-L3-P, EN-IT-L4-A versus EN-IT-L4-late, EN-IT-L4/5-1 versus EN-IT-L5/6-P, EN-ET-L6-A versus EN-ET-L6-P, and EN-ET-SP-A versus EN-ET-SP-P1. The DEGs were detected using the scanpy.tl.rank_genes_groups() function with the default t-test from the Scanpy_1.8.2 in Python (v.3.10). The top ten and bottom ten genes were selected separately on the basis of their log[fold change], and the log₁₀[P] was used to break any ties. The top ten genes represent the most upregulated genes in the DEG analysis, whereas the bottom ten genes represent the most downregulated genes. In the bubble plots, the mean expressions for each gene across all the cells were calculated and represented by the spot colour. The percentage of expressed cells was calculated by dividing the number of cells expressed for each gene by the total number of cells expressed across all the genes, represented by the spot size.

Identification of area-enriched genes

To identify area-enriched genes (Fig. 3h), we analysed EN-ET and EN-IT cells from GW20 and GW22 samples in the annotated fan-shaped regions. Cells in the PFC and the PMC/M1 were grouped as anterior regions and cells in the Par and the V2 were grouped as posterior regions. We identified anteriorly and posteriorly enriched genes by t-tests using the scanpy.tl.rank_gene_groups function. The top ten genes with the highest anterior enrichment and posterior enrichment were selected separately on the basis of their P values. Moreover, the top five genes with the most enrichment exclusively in the Temp were detected using the same methods. EN-Mig cells, RG and IPCs from GW20 and GW22 samples were analysed using a similar method in Extended Data Fig. 10d,e.

Clustering by individual gestational age and single-cell significant hierarchical clustering

To compare neuronal subtypes at different gestational ages, we used an alternative clustering strategy in Figs. 4, 5 and Extended Data Fig. 11, for which only cells from the same gestational age were clustered together. Using the preprocessed data, all cells belonging to samples from the same gestational week were combined. Subsequently, hierarchical clustering was conducted for each gestational week independently. Analogous to the method applied for the integrated analysis, for each time point, we first clustered the cells into 8 H1 and 40 preliminary H2 clusters using sklearn.cluster.KMeans(). H1 and H2 clusters were annotated using the same strategy as for integrated clustering described above.

The single-cell significant hierarchical clustering (scSHC) pipeline (v.0.1.0) was used to determine in an unbiased manner the number of significant H3 subclusters⁴⁰. scSHC adopts a hypothesis-testing approach, recursively dividing the cells into two groups and testing the significance of each split with an adjusted threshold to control the family-wise error rate. It defines Ward linkage as test statistics and uses parametric bootstrapping to estimate its null distribution. If the P value at a given node is larger than the adjusted threshold, the two subsequent branches are merged. For each gestational age, we randomly sampled 500,000 cells and ran scSHC in each H2 cluster. To ensure model robustness, we increased the number of cells used for null distribution estimation from 1,000 in the original method to one-third of the cluster size. Furthermore, to avoid overly small subclusters, if one subcluster comprised fewer than 10% of the cluster size, we stopped further splitting. For each run of scSHC, the following parameters were used: alpha, the family-wised error rate, was fixed at 5 × 10^–4; num_features, the number of genes used as features, was set to 300 (all MERFISH panel genes); num_PCs, the number of top principal components retained for gene expression dimension reduction, was set to 30. The H3 subclusters generated by scSHC are summarized in Supplementary Table 7. Experiments that were subsequently added, including those using the expanded 960 gene panel, were individually clustered using the same method, as illustrated in Supplementary Table 9.

scSHC subcluster correspondence analysis across gestational ages

To evaluate the transcriptomic correspondence of EN-ET and EN-IT scSHC subclusters across different ages, we applied XGBoost (v.2.0.3)⁴¹, a distributed gradient-boosted decision-tree-based classification method. For each adjacent pair of gestational stages, we trained a XGBoost classifier to learn cluster labels from gene-expression data in the earlier gestational stage. Consequently, we used the classifier to classify cells from the later gestational stage. The correspondence between the original clusters and the classifier-assigned labels in the later gestational age dataset were used to map clusters between ages. The general classification workflow is described below and was applied to each adjacent gestational age pair.

Let X_E represent the earlier gestational stage dataset grouped into N_E clusters, and X_L denote the later stage dataset grouped into N_L clusters. Two gene-expression matrices are normalized and log-transformed. The main steps are as follows:

1. 1.

   We trained a multiclass XGBoost classifier on X_E using all 300 genes as features. For clusters with fewer than 15,000 cells, we upsampled by bootstrapping to 15,000 cells to make the dataset more balanced. XGBoost classifier parameters were set to the following values: “objective”: “multi:softmax”; “eval_metric”: “mlogloss”; “num_class”: \(N_\rmE\); “eta”: 0.2; “max_depth”: 20; “subsample”: 0.6; “num_boost_round”: 1000.

2. 2.

   For validation, we randomly sampled 80% of cells in each cluster of X_E to train the classifier and predicted the cluster labels for the remaining 20% of cells. For each testing cell k, the classifier returns a cluster a.

3. 3.

   Assignment probability vector \(p\in \mathbbR^N_\rmE\) and the final cluster label is assigned as cluster (k) = argmax_i p_i. For EN-ET cells, the classifiers achieved over 90% accuracy (96.17% between GW15 and GW20, 90.40% between GW20 and GW22, and 92.98% between GW22 and GW34), whereas for EN-IT cells, the accuracy was more than 83% (93.65% between GW15 and GW20, 88.96% between GW20 and GW22, and 83.09% between GW22 and GW34).

4. 4.

   For prediction, we re-trained the classifier on 100% of cells in X_E and applied it to X_L to obtain the predicted label for each cell. The overall correspondence flows between clusters across gestational ages were visualized using Sankey diagrams (Extended Data Fig. 11).

To identify genes exhibiting upregulation or downregulation across gestational ages in one layer-based EN group, we calculated the expressed cell fraction and mean expression for each gene at each time point. We then fit gene-specific linear-regression models for each of these two metrics, using time points as input, for which the time points were encoded as 1, 2, 3 and 4 for GW15, GW20, GW22 and GW34, respectively. Genes were ranked according to the regression coefficient, with the top ten positive ones indicating upregulated genes and the top ten negative ones indicating downregulated genes. The analyses were repeated for each EN group, and the expression patterns of identified genes were visualized in dot plots (Extended Data Fig. 12a). Cluster-specific expression levels of identified genes are shown in Extended Data Fig. 12b.

To pinpoint genes that drive subcluster specification, we first created a pseudobulk for each cluster in the same EN group and gestational age by averaging cell expression. Subsequently, if there was more than one cluster in the same EN group and gestational age pair, we computed the pseudobulk expression variance for each gene. Next, we selected genes with expression variances greater than 1 in at least one EN group and gestational age pair and showed their expression variances using a heatmap (Extended Data Fig. 12d). We also picked the top ten genes with the highest expression variance in each EN group and gestational age pair and visualized their expression patterns as dot plots (Extended Data Fig. 12e,f).

Visium spatial transcriptomic analysis

Visium spatial gene expression analysis (10x Genomics) was performed according to the manufacturer’s user guide and demonstrated protocol (‘Visium CytAssist Spatial Gene Expression for Fresh Frozen–Methanol Fixation, H&E Staining, Imaging & Destaining, CG000614, Rev A; Visium CytAssist Spatial Gene Expression Reagent Kits User Guide, CG000495, Rev D’). In brief, 10-µm-thick cryosections adjacent to the sections used for MERFISH analyses were adhered to Superfrost slides. Two samples were used for Visium analysis: FB080-O1 (a GW20 occipital cortex sample containing the V1 and V2) and FB121-F1 (a GW20 prefrontal cortex sample). FB080-O1 was coded as A1 and FB121-F1 was coded as D1 in the data-processing and analyses steps.

The tissues were fixed in chilled methanol at −20 °C for 30 min and stained with H&E. The slides were mounted with 85% glycerol in water. H&E staining was imaged at ×20 magnification using a Zeiss Axioscan 7 microscope. After imaging, the slides were placed in Visium tissue slide cassettes with a 6.5 mm gasket and destained with 0.1 N HCL for 15 min at 42 °C. Probe hybridization was performed overnight at 50 °C on a thermal cycler using the provided slide adaptor, followed by post-hybridization buffing. A probe ligation mix was then added into the slide gasket and incubated for 1 h at 37 °C, followed by a post-ligation wash. RNA digestion and tissue removal were performed using a CytAssist instrument (10x Genomics) before probe extension and probe elution. Pre-amplification, SPRIselect and library construction were performed according to the user’s guide. GEX post-library construction quality control was performed using an Agilent 4200 TapeStation System.

Visium libraries were sequenced on a NovaSeq6000 according to recommended parameters by 10x Genomics to the depth of 214–235 million PE150 (paired-end reads with a read length of 150 base pairs) reads. Demultiplexed reads were processed using Space Ranger (v.2.1.0) and Spaceranger-count with default parameters and manual alignment. Visium data analysis was conducted separately on two slices: A1 (FB080-O1) and D1 (FB121-F1). For visualization and clustering, Seurat_5.0.1 and R (v.4.3.1, 2023-06-16) were used on the BCH compute nodes. SCTransform normalization was performed with the ‘spatial’ assay parameter. Optimal clusters were generated, and clusters with negligible cells were removed. Marker identification was performed using FindAllMarkers analysis from the Seurat package using Wilcoxon rank-sum test with the following parameters: (1) only return positive markers; and (2) only test genes that are detected in a minimum fraction of 1% of cells in either of the two populations. This was the default. After these results, we filtered for markers grouped by clusters with an average log₂[fold change] greater than 1. We generated different visualizations using top genes from this analysis. R (v.4.1.2) along with Seurat_4.1.1 were used for both exploratory analyses and final visualization on a local machine. Spatial feature plots were generated to visualize the expression of genes of interest in the tissue sections.

Ferret tissue preparation and immunohistochemistry

All procedures on ferrets were performed under protocols approved by the Institutional Animal Care and Use Committee at Boston Children’s Hospital. Ferrets (Mustela putorius furo) were obtained from Marshall BioResources and were housed in a vivarium under a 12-h light–dark cycle. Food and water were available ad libitum.

Ferrets were deeply anaesthetized with ketamine (50 mg ml^–1) and xylazine (20 mg ml^–1) by intraperitoneal injection and transcardially perfused with 0.9% NaCl solution (w/v) followed by 4% PFA in PBS. Brains were dissected out, post-fixed overnight in 4% PFA in PBS at 4 °C, cryoprotected by serial overnight incubations with 15% and 30% sucrose in PBS and sectioned frozen at 50 µm on a sliding microtome (Leica SM2010 R). Free-floating brain slices were permeabilized with 0.25% Triton X-100 in PBS for 4 times for 15 min each and then blocked in a solution containing 5% BSA (Sigma-Aldrich, 05470), 0.3% Triton X-100, and 10% normal donkey serum (Abcam, ab7475) for 2 h. Brain slices were then incubated overnight at 4 °C with primary antibodies on an orbital shaker.

The next day, slices were washed for 4 times for 15 min each in PBS and then incubated with secondary antibodies for 2 h at room temperature on an orbital shaker. All primary and secondary antibodies were diluted in a solution containing 1% BSA, 0.3% Triton X-100 and 5% normal donkey serum. The following antibodies were used: mouse anti-RORB (1:500, R&D Systems, PP-N7927-00), rabbit anti-NPY (1:500, Abcam, ab30914), rat anti-CTIP2 (1:500, Abcam, ab18465), Alexa Fluor 488 donkey anti-mouse (1:500, Invitrogen), Alexa Fluor 555 donkey anti-rabbit and Alexa Fluor 647 donkey anti-rat. Slices were counterstained with DAPI solution. Slices were mounted on microscope glass slides and coverslipped with Vectashield Vibrance antifade mounting medium. Imaging of ferret brain slices was performed using a LSM980 confocal microscope at 8-bit depth with a ×10 objective and 0.6 digital zoom. Tile scan images were acquired and stitched.

Fetal macaque tissue preparation and single-molecule fluorescent in situ hybridization

Preserved cryosections of fetal rhesus macaque monkeys at E93 and E76 were shared from a previous study⁴⁶. Macaque fetal brains were dissected and immersed in 4% PFA overnight at 4 °C. Fixed brain blocks were immersed in step-gradients of sucrose–PBS up to 30% for 2–3 days at 4 °C, then embedded in OCT and frozen at −80 °C. Sagittal sections (for E76) and coronal sections (for E93) were cut at 25 μm using a Leica CM3050S cryostat and stored at −80 °C until use. Immunostaining of fetal macaque samples was performed as described above, but with an extra step of antigen retrieval before permeabilization. Antigen retrieval was performed by heating slides in a steam cooker for 30 min in Retrievagen A (pH 6.0) solution (BD Biosciences), followed by two washes of PBS wash. Single-molecule fluorescent in situ hybridization was performed on macaque sections using a RNAscope Multiplex Fluorescent V2 assay (ACD). A custom RNAscope probe for macaque NPY was designed and made by ACD (RNAscope Probe, Mmu-NPY). The manufacturer’s standard protocol for multiplex fluorescent in situ hybridization using a HybEZ Hybridization system was followed. Whole tissue mRNA in situ hybridization imaging was performed on a Zeiss LSM 980 confocal microscope with image tiling and stitching, z stack of 10 µm and 2× averaging at ×20 magnification.

Re-analysis of a macaque scRNA-seq dataset

scRNA-seq data (Macaque.dev.rds Seurat object) for macaque fetal cortex⁴⁶ were downloaded from http://resources.sestanlab.org/devmacaquebrain/. The standard workflow for visualization and clustering was executed using Seurat_5.0.1 and R (v.4.3.1). To refine the analysis of ENs, we used the cell_subtype annotations from the pre-processed dataset. Specifically, we subclustered cells labelled as ‘enIPC’ and ‘Excitatory neurons’. This iterative subclustering process continued until we identified clusters exclusive to the V1. In these V1-specific clusters, we assessed the expression patterns of NPY, PDZRN4 and ABI3BP. Dimensionality reduction plots were generated to visualize the subclusters, and DGE was performed using the FindAllMarkers function to identify key markers defining each subcluster. To evaluate the representation of the V1 among the identified subclusters, contribution plots were generated, quantifying the proportion of cells from the V1 in each cluster.

V1-enriched and V2-enriched gene set analysis

To identify sets of genes that are enriched in the V1 or V2 at different ages (GW15, GW20, GW34 and adult), we conducted DEA between V1-specific and V2-specific clusters in EN-ITs and EN-ETs separately at GW15, GW20, GW34 and adult stages (Fig. 5d–g and Supplementary Table 11). Lists of V1-specific and V2-specific clusters are provided below:

• V1-specific clusters in EN-ITs. GW20: ’EN-IT-L3-c1 and ’EN-IT-L3/4-c4. GW34: EN-IT-L3/4-c2, EN-IT-L3/4-c3, EN-IT-L3/4-c4, EN-IT-L3/4-c5, EN-IT-L3/4-c6, EN-IT-L3/4-c7, EN-IT-L4-2-c2 and EN-IT-L4-2-c3. Adult: EN-L4c-c1, EN-L4-c0, EN-L4-c1, EN-L4-c2, EN-L4-c3 and EN-L3-2-c2.

• V2-specific clusters in EN-ITs. GW20: ’EN-IT-L3-c4 and ’EN-IT-L3/4-c1. GW34: EN-IT-L3/4-c1, EN-IT-L3-c0 and EN-IT-L4-2-c1. Adult: EN-L4c-c0 and EN-L3-2-c1.

• V1-specific clusters in EN-ETs. GW20: ’EN-ET-L5/6-c5 and ’EN-ET-SP-2-c2. GW34: EN-ET-L6-c1 and EN-ET-L6-V1-c3. Adult: EN-L2/3-2-c0.

• V2-specific clusters in EN-ETs. GW20: ’EN-ET-L5/6-c4 and ’EN-ET-SP-2-c4. GW34: EN-ET-L5-c1. Adult: EN-L2 | 3-1-c0.

For GW15, no V1-specific or V2-specific clusters were found, and DEA was conducted across the annotated V1–V2 regions. The DEGs were detected using the scanpy.tl.rank_genes_groups() function with t-tests. Genes with adjusted P < 0.05 and log[fold change] > 0.5 were selected as significant DEG sets. For each annotated V1 or V2 region at every time point, the log[fold change] for the combined expression of the DEG set was calculated, and the mean log[fold changes] were visualized as heatmaps.

Nucleus isolation and snRNA-seq

Crysoections (100 µm thick) consecutive to the ones collected for MERFISH and Visium analyses were collected in a 1.5 ml tube. Nuclei were isolated from the cryosection as previously described but with minor modifications⁶⁹. In brief, one cryosection per sample was resuspended in 1 ml homogenization buffer with additives (10 mM Tris buffer pH 8.0, 250 mM sucrose, 25 mM KCl, 5 mM MgCl₂, 0.1% Triton X-100, 0.1 mM DTT, 1× complete, mini, EDTA-free protease inhibitor cocktail (Roche 11836170001) and 25 µl Protector RNAse inhibitor (Roche 3335399001, 0.2 U µl^–1)) and transferred to a 7 ml douncer and dounced 10 times with a ‘tight’ pestle. Homogenized nuclei were spun for 10 min at 900g at 4 °C, then washed once with blocking buffer (1× PBS pH 7.4 and 1% BSA), centrifuging for 5 min at 400g at 4 °C. All centrifugations were done in a bucket centrifuge. Nuclei were resuspended in 300 µl blocking buffer with Protector RNAse inhibitor (Roche 3335399001, 1 U µl^–1) and DAPI (final concentration 1 µg ml^–1) and passed through a 40 µm filter. DAPI-positive nuclei (13,500 nuclei per 10x reaction) were fluorescence-activated nuclei sorted directly into Chromium Next GEM Single Cell 3′ Reagent Kit v3.1 GEM master mix (Step 1.1) minus RT Enzyme C. RT Enzyme C was then added and reactions loaded into a Chromium Next GEM Chip and snRNA-seq libraries were generated according to the Chromium Next GEM Single Cell 3′ Reagent Kit v.3.1 manual with 11 PCR cycles for cDNA amplification and 11–12 PCR cycles used for final amplification of gene expression libraries. Five reactions were performed for FB080-O1 and three reactions were performed for FB121-F1 samples. Three reactions were performed on FB080-O2, a consecutive tissue block posterior to FB080-O1, near the occipital pole of the brain. snRNA-seq libraries were sequenced on a NovaSeq6000 according to the recommended parameters of 10x Genomics to the depth of 909 million PE150reads for FB080-O1, 541 million PE150 reads for FB080-O2 and 484 million PE150 reads for FB121-F1.

snRNA-seq processing, clustering and analysis

Demultiplexed reads were processed in Cell Ranger (v.7.2.0) with cellranger-count and default parameters, which resulted in 41,740 estimated cells for FB080-O1, 23,001 estimated cells for FB080-O2 and 26,161 estimated cells for FB121-F1. For visualization and clustering, Seurat_5.0.1 and R (v.4.3.1) were used. Clustering optimization was achieved through iterative application of the FindClusters algorithm to attain optimal cluster resolution. Feature plots were constructed for genes of interest, and integration into the cell-by-gene tool was verified. Contribution plots were generated to observe the contribution of cell types to the overall dataset. To find DEGs between V1-enriched (EN-IT-L4-V1) and V2-enriched (EN-IT-UL-2) clusters, we performed Wilcoxon rank-sum tests at log₂[fold change] > 0.5 and Bonferroni-corrected P < 0.05. Only test genes that were detected in a minimum fraction of 25% of cells in either of the two populations were considered, filtering the genes that were infrequently expressed and only returning positive markers. To verify that we captured DEGs between the V1 and V2, we also identified DEGs between V1 and V2 regions in Merscope, which revealed overlap only between V1-upregulated genes (50%) or V2-upregulated genes (73%) but not vice versa (0%). For GO analysis, we used Bioconductor (v.3.19) to perform GO analyses in R (https://bioconductor.org/packages/release/bioc/html/motifmatchr). Specifically, DEGs upregulated in the V1 or V2 were both tested for GO enrichments using a background of all genes tested for differential expression (odds ratio > 2 and FDR < 0.01).

ENVI imputation from scRNA-seq to Merscope

We used ENVI (v.0.1.0)⁴⁷ on our MERFISH and snRNA-seq data to expand the MERFISH gene panel. ENVI leverages a conditional variational autoencoder to integrate spatial and snRNA-seq data, which can simultaneously impute missing gene expression for genes that are not included in the spatial transcriptomic data and transfer spatial information onto the snRNA-seq data. We selected the top 1,000 highly variable genes from snRNA-seq, along with the following genes of interest for imputation: ABI3BP, PDZRN4, FLRT2, TAFA2, NR1D1, IL1RAP, CCBE1, THSD7B, TRPC6, CHRM2, LUZP2, LRP1B, LRRTM4, HDAC9, FBXL7, DTNA, SYNDIG1, SDK1, LMO3, TRIQK, UNC13C, CNTNAP2, KCNIP4, PDZRN3, DLX6, DLX6-AS1, ADARB2, ERBB4, NRXN3, DLX2, ZNF536, PRKCA, THRB, TSHZ1, PBX3, MEIS2, CALB2, CDCA7L, SYNPR, SP8, CASZ1 and FOXP4. We ran ENVI with the following key parameters: “k_nearest”: 100; “num_cov_genes”: 50; “num_layers”: 3; “num_neurons”: 1024; “latent_dim”: 512; “spatial_dist”: “pois”; “sc_dist”: “nb”; “cov_dist”: “OT”; “prior_dist”: “norm”; “spatial_coeff”: 1; “sc_coeff”: 1; “cov_coeff”: 1; “kl_coeff”: 0.3.

Confidence interval analysis for imputed gene expression

To quantify the imputation uncertainty, we used transcript imputation with spatial single-cell uncertainty estimation (TISSUE) (v.1.0.1)⁴⁸ to generate the 95% confidence intervals for imputation values. TISSUE assumes that neighbouring cells of the same cell type share similar gene-expression profiles. It defines the cell-centric variability measure and adopts a conformal inference framework to calibrate uncertainties for new expression predictions. We ran TISSUE with the following key parameters: “grouping_method”: “kmeans_gene_cell”, “k”: 4, “k2”: 2, “alpha_level”: 0.05. After obtaining the confidence intervals, we truncated the lower bounds to 0 and computed interval widths. These widths were then visualized for each gene using box plots (Extended Data Fig. 18a).

Constellation plot analysis

The constellation plot was generated using the ENVI imputation results. We selected all RG, IPC, EN-Mig, EN-ET and EN-IT cells from the V1 and V2 of FB080-O1c and the PFC of FB121-F1. On the basis of the H2 annotations, we combined all EN-ET cells from layers 5 and 6 into a single node (EN-ET L5/6) and did the same for layers 5 and 6 EN-IT cells (EN-IT L5/6). Other EN-ET and EN-IT cells were categorized into nodes according to their H2 annotations, whereas RG and IPC cells were classified by their H1 annotations. Cells were further subclassified on the basis of the sample and annotated region in which they were located. The positions of the cell types on the constellation plot were determined using the median UMAP coordinates of all cells of each type. Node sizes were scaled according to the number of cells of each type. To establish edges between pairs of cell types, we first performed principal component analysis on the ENVI gene expression matrix. Using the top 50 principal components, we then identified the 15 nearest neighbours for each cell by applying the Ball Tree algorithm. For any given pair of cell types, A and B, we added a connection from A to B if a cell in A had a nearest neighbour in B. We then calculated the total number of connections from A to B (n_A→B) and the total number of connections from A to all types (n_A). An edge was drawn between A and B on the constellation plot if both n_A→B/n_A and n_B→A/n_u exceeded 0.02. The edge from A to B was coloured according to the proportion of connections from A to B relative to the total connections from A (n_A→B/n_A).

Cell–cell communication analysis

We used the CellChat package (v.1.6.1) to decipher cell–cell communication networks using its extensive database of known human ligand–receptor interactions. To assess the communication probability, we followed the methodologies outlined in the original publication⁴⁹. This method applies the law of mass action to the average expression levels of ligands and receptors across annotated cell groups and evaluates their significance using the random permutation statistical test. Interactions with P < 0.05 were considered significant and those identified in fewer than ten cells in each cell group were discarded to enhance the robustness of our findings. For this study, we analysed all cell types available in our dataset but specifically focused on EN-IT-L4-V1, EN-IT-UL-1, EN-IT-UL-2, EN-IT-UL-3, EN-IT-UL-4 and EN-IT-UL-5 as receivers in the NRXN pathway. Visualization of the communication networks and pathways was performed using the native plotting functions in CellChat, which facilitated a comprehensive representation of the cellular interactions in our dataset.

Interactive web-based browser for MERFISH data

To create a web browser of our MERFISH data, we used the WebAtlas pipeline, which leverages Vitessce (v.1.0.15) to visualize spatial transcriptomic data, in the cloud-optimized Zarr format and the data are uploaded to Amazon Web Service (AWS) S3 buckets for easy access and retrieval⁷⁰. WebAtlas retrieves the data directly from S3 using its web portal hosted by the Wellcome Sanger Institute, which enables seamless integration and efficient visualization by leveraging cloud storage. WebAtlas provides interactive visualizations including UMAP plots, cell-type annotation and gene-expression distributions, all powered by AWS S3 cloud storage. Examples of visualization can be found at the links below:

FB121-F1 (GW20 PFC experiment): https://webatlas.cog.sanger.ac.uk/dev/index.html?theme=dark&config=https://f1-fb121.s3.us-east-1.amazonaws.com/Merfish_F1_FB121/0.5.2/F1_FB121-scRNAseq-config.json.

FB123-F1 (GW22 PFC experiments): https://webatlas.cog.sanger.ac.uk/dev/index.html?theme=dark&config=https://f1-fb123.s3.us-east-1.amazonaws.com/Merfish_F1_FB123/0.5.2/F1_FB123-scRNAseq-config.json.

Statistics and reproducibility

Replicates and statistical tests are described in the figure legends and Methods. Sample sizes were not predetermined using statistical methods. Tissue samples were not randomized, nor were the investigators blinded during collection as no subjective measurements were taken. Data for snRNA-seq, Visium and MERFISH were collected from all available samples and no randomization was necessary. To identify DEGs between clusters, Wilcoxon rank-sum tests were performed.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.