Inhibitory specificity from a connectomic census of mouse visual cortex

In mammalian cortex, information processing involves a diverse population of neurons distributed across six layers in an arrangement described as a cortical column³. Cell types are a central concept for understanding how the columnar network is organized¹. Originally classified on the basis of morphology⁴, cortical cell types have been increasingly characterized by transcriptomic, molecular, electrophysiological and functional properties as well^5,6,7,8,9. Excitatory neurons make up almost 90% of neocortical neurons¹⁰ and vary not only across cortical layers but also by long-range projection targets¹¹. Inhibitory neurons, although much fewer in total number, have at least as much diversity as excitatory neurons in a single region^6,7,8, offering the potential for highly selective control of cortical activity.

Determining how fine-scale cell type definitions are reflected in synaptic connectivity remains difficult. Most of our understanding of inhibitory connectivity is based not on individual cell types but on cardinal subclasses based on marker genes parvalbumin (PV), somatostatin (SST), vasoactive intestinal polypeptide (VIP) and Id2, each with shared developmental, functional and synaptic properties^12,13,14,15. Within these cardinal subclasses, individual cell types can be highly diverse^6,7,16 and functionally distinct¹⁷, but little is known about connectivity for most cell types. Although some studies have observed largely unspecific connectivity onto nearby cells¹⁸, others have found examples of selective targeting of subpopulations of excitatory cells based on the layer¹⁹ or long-range axonal projection of target cells^20,21. It is not known whether such selectivity is common or rare relative to unspecific connectivity. Likewise, basic organizational properties remain unclear: for example, which excitatory neurons receive inhibition from the same interneurons.

To date, physiological^22,23 or viral²⁴ approaches to measuring connectivity are still challenging to scale to the full diversity of potential cell type interactions. In smaller model organisms like Caenorhabditis elegans²⁵ and Drosophila melanogaster^26,27, dense reconstruction using large-scale electron microscopy (EM) has been instrumental for discovering cell types and their connectivity. In mammalian cortex, technical limitations on EM volume sizes have meant that similar studies could not examine complete neuronal arbours, making the link between cellular morphology and connectivity difficult to address^28,29,30. However, recent advances in data generation and machine learning have helped to produce EM datasets at the scale of a cubic millimetre, making circuit-scale cortical EM volumes now possible².

In this study, we used a millimetre-scale EM volume of mouse primary visual cortex (VISp)² to reconstruct the anatomy and synaptic connectivity for a continuous population of 1,352 neurons in a column spanning from layer 1 to white matter. The scale of this data, combined with the resolution provided by EM, led us to ask how morphological cell types relate to the synaptic connectivity of inhibitory neurons. Inspired by classical neuroanatomical methods, we classified inhibitory neurons into connectivity-based subclasses largely aligned with molecular subclasses and developed a new classification of excitatory neurons using morphological and synaptic properties, capturing features that were not clear from morphology alone. By analysing the synaptic output of inhibitory neurons at both the single-cell and subclass levels, we found that inhibitory neurons exhibited widespread target specificity and identified groups of interneurons with similar subclass-specific targeting but with different compartmental targeting. Our data not only identified a new class of disinhibitory specialist but also indicate an organizing principle for inhibitory connectivity that is complementary to, but distinct from, cell types: diverse groups of inhibitory neurons that are positioned to collectively control activity of the same target populations with remarkable precision.

To measure synaptic connectivity and neuronal anatomy for a large neuronal population, we used a serial section transmission EM volume of mouse visual cortex acquired as part of the broader MICrONS project². Specifically, we analysed a volume of mouse visual cortex spanning 523 × 1,100 × 820 µm (anteroposterior × mediolateral × depth), covering pia to white matter and including parts of VISp and higher-order visual areas (Fig. 1a–c). Importantly, these dimensions were sufficient to capture the entire dendritic arbour of typical cortical neurons (Fig. 1a) at a resolution capable of resolving ultrastructural features such as synaptic vesicles (Fig. 1b). Convolutional networks generated an initial autosegmentation of all cells, segmented nuclei, detected synapses and assigned synaptic partners². Owing to reduced alignment quality near the edge of tissue, segmentation began about 10 µm from the pial surface and continued into white matter.

a, The millimetre-scale EM volume is large enough to capture complete dendrites of cells across all layers. Neurons shown are a random subset of the volume, with a single example at right for clarity. b, The autosegmented EM data show ultrastructural features such as membranes, synapses and mitochondria. Scale bar, 500 nm. c, Top view of EM data with approximate regional boundaries indicated. The yellow box indicates the 100 µm × 100 µm column of interest. Scale bar, 200 µm. d, All soma locations in the column coloured by cell class. Scale bar, 100 µm. e, Example neurons from along the column. Note that anatomical continuity required adding a bend in deeper layers. f, Proofreading workflow by cell class. g, Cell density for column cells along cortical depth by cell class. Scale bar, 200 µm. h, Input synapse count per micrometre of depth across all excitatory (purple) and inhibitory (green) column cells along cortical depth by target neuronal cell class. Scale bar, 200 µm. i, All excitatory dendrites, with arbours of cells with deeper somata coloured darker. Same orientation as in d. Scale bar, 200 µm. j, Number of input synapses for each excitatory neuron as a function of soma depth. k, All inhibitory dendrites, as in j. l, Number of input synapses for inhibitory neurons, as in k. m, Axons of inhibitory neurons, as in j. n, Number of output synapses for inhibitory neurons, as in k. VISrl, rostrolateral visual area; VISal, anterolateral visual area; Exc, excitatory; inh, inhibitory; non, non-neuronal; syn, synapses; WM, white matter.

Source Data

Full size image

To generate an unbiased sample of cells across all layers, we selected all cells whose soma fell within a 100 × 100-µm-wide column from pia to white matter and centred on the VISp portion of the volume (Fig. 1c–e). This ___location was chosen to be far from dataset edges to avoid truncated arbours as much as possible. To follow a continuous population of neurons, the column bends in lower layer 5, defined such that the apical dendrites of deep layer cells would be intermingled with the cell bodies of superficial cells (Fig. 1d,e and Methods). This trajectory was also followed by primary axons of superficial cells and the translaminar axons of inhibitory neurons, indicating that this bend is shared across cell types.

We classified all 1,886 cells in the column as excitatory neurons, inhibitory neurons or non-neuronal cells on the basis of morphology (Fig. 1d). For neurons, we performed extensive manual proofreading—more than 46,000 edits in all (Fig. 1f), guided by computational tools to focus attention on potential error locations (Methods). We selected a proofreading strategy to efficiently measure the connectivity of inhibitory neurons across all possible target cell types. Proofreading of excitatory neurons aimed to reconstruct complete dendritic arbours, combining both manual edits and computational filtering of false axonal merges onto dendrites (Methods), and for inhibitory neurons we reconstructed both complete dendritic arbours and extensive (but incomplete) axonal arbours.

Consistent with previous reports¹⁰, excitatory cell densities varied between layers, while inhibitory neurons and non-neuronal cells were more uniform (Fig. 1g). Dendritic reconstructions included the locations of a total of 4,490,649 synaptic inputs across all cells. Synaptic inputs onto excitatory cell dendrites were more numerous in layers 1–4 compared to layers 5–6 (Spearman correlation of synapse count with depth: r = −0.92, P = 1.3 × 10⁻¹¹), whereas inputs onto inhibitory cells were relatively uniform across depths (Fig. 1h; Spearman correlation of synapse count with depth: r =−0.06, P = 0.76).

Reconstructions captured rich anatomical information for individual cells across all layers. Excitatory cell dendrites (Fig. 1i) typically had thousands of synaptic inputs, with laminar differences in total synaptic input per cell (analysis of variance for layer effect: F = 82.9, P = 1.6 × 10⁻⁴⁸, Fig. 1j). Typical inhibitory neurons had 10³–10⁴ synaptic inputs (Fig. 1k,l) and 10²–10⁴ outputs (Fig. 1m,n) but did not show strong laminar patterns (analysis of variance for layer effect: F = 0.72, P = 0.53). Collectively, inhibitory axon reconstructions had 427,294 synaptic outputs. Attempts were made to follow every main inhibitory axon branch, but for large inhibitory arbours not every tip was reconstructed to completion; axonal properties should be treated as a lower bound. Comparing to a subset of neurons where reconstruction aimed for completeness³¹, we estimate that typical axonal reconstructions captured 50–75% of their total synaptic output compared to exhaustive proofreading.

Molecular expression is a powerful organizing principle for inhibitory neurons, with four cardinal subclasses having distinct connectivity rules, synaptic dynamics and developmental origins¹⁴. However, EM data have no direct molecular information, nor do simple rules map morphology to molecular identity. Classical neuroanatomical studies often used the postsynaptic compartments targeted by an inhibitory neuron as a key feature of its subclass^12,32: for example, distinguishing soma-targeting basket cells from dendrite-targeting Martinotti cells.

Inspired by this approach, we used the targeting properties of inhibitory neurons to assign cells to anatomical subclasses (Fig. 2a). For all excitatory neurons, we divided the dendritic arbour into four compartments: soma, proximal dendrite (less than 50 µm from the soma), apical dendrite and distal basal dendrite (Fig. 2b; see Extended Data Fig. 1 and Methods for apical classification). Inhibitory cells were treated as a fifth target compartment. For each inhibitory neuron, we measured the distribution of synaptic outputs across compartments (Fig. 2c). We also included two measures of how a cell distributes its synapses onto individual targets: (1) the fraction of all synapses that were part of a multisynaptic connection and (2) the fraction of synapses in a multisynaptic connection that were close together along the axon (‘clumped’; Fig. 2d). We used a distance threshold of 15 µm, about a quarter of the circumference of a typical cell body, and measurements were robust to the exact value (Extended Data Fig. 2). We use the term ‘connection’ to indicate a pre- and postsynaptic pair of cells connected by one or more distinct synapses and ‘multisynaptic connection’ for a connection with at least two synapses. We trained a linear classifier on the basis of expert annotations of the four cardinal subclasses for a subset of inhibitory neurons and applied it to all cells (Fig. 2d and Extended Data Fig. 3).

a, Example of an inhibitory axon making synaptic outputs (green dots) onto specific locations on a target pyramidal cell (purple). b, Dendritic compartment definitions for excitatory neurons. c, Cartoon definition for a multisynaptic connection (left) and the synapses in the multisynaptic connection considered ‘clumped’ along the presynaptic axon (right). d, Targeting features for all inhibitory neurons, measured as fraction of synapses onto column cells (for fraction clumped only: synapses in multisynaptic connections). e, Relationship between anatomical connectivity categories (top), typical associated classical cell categories (middle) and anatomical examples (bottom) of the inhibitory subclasses. Dendrite is darker, axon lighter. f, Adjacency matrix for inhibitory neurons. Each dot represents a connection from a presynaptic to a postsynaptic cell, with dot size proportional to synapse count. Dots are coloured by presynaptic subclass and ordered by subclass, connectivity group (Fig. 5) and soma depth. g, Standard model of inhibition of inhibition between molecular subclasses. h, Mean number of synaptic inputs a postsynaptic cell received from all cells of a given presynaptic subclass. i, Potential InhTC targets. j, Synaptic output fraction each InhTC (columns) places onto target subclasses (rows). InhTCs are clustered into two subtypes: one targeting DistTCs (InhTC^dist) and another targeting PeriTCs (InhTC^peri). k, Connectivity diagram for InhTC^peri suggested by data. l, Morphology of example InhTC^dist. m, Morphology of all InhTC^peri. n, Median synapse size (arbitrary units measuring voxels in segmented cleft) from InhTC^dist (left) and InhTC^peri (right) onto inhibitory subclasses. Error bars indicate 95% confidence interval. T-test P-values indicated: *, P < 0.05; ***, P < 0.005 after Holm–Sidak correction. o, Distribution of synapses per connection for InhTC^peri and InhTC^dist onto their preferred and non-preferred targets. Scale bars, 500 µm. CCK, cholecystokinin; frac, fraction; multisyn, multisynaptic; no, number of.

Source Data

Full size image

We named each subclass on the basis of its dominant anatomical property: perisomatic targeting cells (PeriTC) that primarily target soma or proximal dendrites, distal dendrite targeting cells (DistTC) that primarily target distal basal or apical dendrites, sparsely targeting cells (SparTC) that make few multisynaptic connections and inhibitory targeting cells (InhTC) that primarily target other inhibitory neurons (Fig. 2e). Typical examples of each subclass correspond roughly to classical or molecular subclasses (Fig. 2e), but there is not a one-to-one match^12,14. PeriTCs would include soma-targeting cells from multiple molecular subclasses (for example, both PV and CCK⁺ basket cells)³³. DistTCs would include SST⁺ Martinotti and non-Martinotti cells but also any neuron that strongly targets apical dendrites. InhTCs align well with disinhibitory specialist VIP neurons. The SparTC subclass included both neurogliaform cells and all layer 1 interneurons, indicating that it largely contained cells from the Id2 class¹⁵. Note that some cell types, such as chandelier cells, had no examples in the column, and some column cells did not fall into clear classical categories.

Numerous studies have identified a standard architecture for the inhibition of inhibition at the subclass level³⁴: PV neurons inhibit other PV neurons, SST neurons inhibit all other subclasses (but not themselves), and VIP neurons inhibit SST neurons (Fig. 2h). Variations on this broad pattern have been found; for example, VIP⁺ neurons have been shown to target both SST and PV cells³⁵, but little is known about the relationship between these connections and individual cells. The EM data contained 9,235 synapses between pairs of inhibitory neurons across 3,569 distinct connections (Fig. 2g and Extended Data Fig. 4), allowing us to examine whether single-cell resolution offered new insights into circuit organization.

To validate reconstructions and labels, we first measured inhibitory connectivity at the level of cardinal subclasses. As a proxy for presynaptic influence of a subclass, we computed the average number of synapses between all neurons from each presynaptic subclass onto each inhibitory neuron and averaged them within postsynaptic subclass. The five expected subclass-level connections aligned with the five strongest connections measured from EM (Fig. 2i), on the basis of the approximate correspondence (Fig. 2e). Both cardinal subclass identification and neuronal reconstructions were thus consistent with established connectivity.

At the level of individual cells, however, the data showed new connectivity patterns. We focused on InhTCs, ‘disinhibitory specialists’ that almost exclusively target other inhibitory neurons rather than excitatory cells (mean: 82% of synaptic outputs). For each InhTC, we computed its distribution of synaptic outputs across inhibitory subclasses (Fig. 2j,k). VIP-positive disinhibitory specialists in visual cortex have been shown to preferentially target SST cells^23,34, and thus we expected InhTCs would largely target DistTCs.

As expected, synaptic output was principally onto DistTCs for 21 of 29 InhTCs (mean of 74% of those synapses onto inhibitory neurons), a group we denoted InhTC^Dist (Fig. 2k). Single-neuron consideration of InhTC^Dist connectivity showed striking laminar organization, with InhTC^Dist in layers 2–4 targeting those DistTCs in layers 4 and 5 but not those in layer 2/3 (Extended Data Fig. 5). Those DistTCs in layer 2/3 made few synapses onto InhTC^Dist in return. Interestingly, layer 2/3 DistTCs typically targeted excitatory neurons in upper (‘layer 2’) but not lower (‘layer 3’) layer 2/3, indicating that InhTC-mediated disinhibition differs across layer 2/3 pyramidal cells.

Unexpectedly, we also found a second population of disinhibitory specialists. This smaller group of InhTCs (8 of 29) specifically targeted PeriTCs (mean of 82% of those synapses onto inhibitory neurons), and hence we called them InhTC^Peri (Fig. 2k). Although InhTC^Dist had bipolar or multipolar dendrites and were concentrated in layers 2–4 (Extended Data Fig. 5), consistent with typical VIP neurons (Fig. 2m), InhTCs^Peri all had multipolar dendrites and were distributed across layers (Fig. 2n). The eight InhTC^Peri in the column targeted 56 of 58 PeriTCs with a mean of 10.5 net synapses per target cell, indicating that this connectivity probably includes basket cells from PV and other molecular subclasses (Extended Data Fig. 4). We next asked if InhTC^Peri receive reciprocal inhibition from PeriTCs in analogy to the reciprocal inhibition between VIP and SST cells^23,34. However, we found few reciprocal synapses from PeriTCs back onto InhTC^Peris but numerous inhibitory inputs from DistTCs (Extended Data Fig. 5), suggesting a new pathway on the standard inhibitory diagram (Fig. 2k).

The targeting preference of InhTCs was seen across several aspects of their connectivity. We first looked at a measure of synapse size on the basis of the automatic synapse detection (Methods). InhTC^Dist → DistTC synapses had a median size 44% larger than those onto other inhibitory subclasses (Fig. 2o). Similarly, the median InhTC^Peri → PeriTC synapse was 69% larger than synapses onto other inhibitory subclasses (Fig. 2o). In addition, the mean number of synapses per unique connection was significantly higher between InhTC^Dist → DistTC compared to other targets (Fig. 2p) (3.6 versus 1.6 synapses per connection; P = 1.5 × 10⁻¹⁰, Student’s t-test) and between InhTC^Peri → PeriTC compared to other targets (3.1 versus 1.5 synapses per connection; P = 1.1 × 10⁻⁵, Student’s t-test). The ___location of synapses onto preferred targets was similar for the two InhTC subgroups, with a median distance from soma of 83.5 µm (InhTC^Dist) and 86.2 µm (InhTC^Peri) and no significant difference in distribution (Kolmogorov–Smirnov test, P = 0.25) (Extended Data Fig. 5). Taken together, both InhTC^Dist and InhTC^Peri express their distinct targeting through increased synapse count, larger synapses and more synapses per connection.

Although inhibitory neurons have frequently been described as having dense, non-specific connectivity onto nearby neurons³⁶, many studies have shown examples not only of layer-specific connectivity³⁷ but also of selectivity in spatially intermingled excitatory subpopulations^20,21,38. It is unclear the degree to which inhibition is specific, and in general, the principles underlying which excitatory neurons are inhibited by which inhibitory neurons are not well understood.

To address these questions, we first anatomically characterized excitatory neurons subclasses in the EM data. Previous approaches to data-driven clustering of excitatory neuron morphology used dendritic shape^8,9, but the EM data also has the ___location and size of all synaptic inputs (Fig. 3a). We reasoned that such synaptic features would help characterize the landscape of excitatory neurons, because synapses directly reflect how neurons interact with one another. We assembled a suite of 29 features to describe each cell, including synapse properties such median synapse size, skeleton qualities such as total branch length and spatial properties characterizing the distribution of synapses with depth (Fig. 3b and Methods). The synapse-detection algorithm did not distinguish between excitatory and inhibitory synapses, and thus all synapse-based measures include both types of synapses. We performed an unsupervised consensus clustering of these features (Fig. 3b–d), identifying 18 ‘morphological types’ (M-types; Methods).

a, Morphology (black) and synapse (cyan dots) properties were used to extract features for each excitatory neuron, such as this layer 2/3 pyramidal cell. b, Heatmap of Z-scored feature values for all excitatory neurons, ordered by anatomical cluster (Fig. 5) and soma depth. Anatomical properties were tip length, tortuosity, dendritic and somatic synapse counts, total path length, radial extent, median synapse distance from soma, somatic and dendritic synapse sizes, dynamic range of synapse size, shallowest and deepest ranges of synapse depth, range of synapse depths, linear synapse density and dendritic radius. All synapse measures use synaptic inputs only. See Methods for detailed feature descriptions. c, Uniform manifold approximation and projection (UMAP) of neuron features coloured by anatomical cluster. Inset shows number of cells per cluster. d, Example morphologies for each cluster. Scale bar, 500 µm. e, Soma depth of cells in each anatomical cluster. f, Median linear density of input synapses across dendrites by M-type. g, Median synapse size (Methods). In f and g, coloured dots indicate single cells; black dots and error bars indicate a bootstrapped (n = 1,000) estimate of the median and 95% confidence interval. a.u., arbitrary units.

Source Data

Full size image

To relate this landscape to known cell types, we compared M-type classifications to expert labels of layer and long-range projection type (intratelencephalic/intracortical (IT); extratelencephalic/subcortical-projecting (ET), near-projecting (NP), corticothalamic (CT))³⁹. Each layer contained several M-types, some spatially intermingled and others separating into subdomains in the layer (Fig. 3e). M-types were named by the dominant expert label (Extended Data Fig. 6), with M-types in the same layer being ordered by projection subclass and average soma depth. For clarity, we use the letter ‘L’ in the name of M-types (which may include some cells outside the given layer) and the word ‘layer’ to refer to a spatial region. Upper and lower layer 2/3 emerged as having distinct clusters, which we denoted ‘L2’ and ‘L3’, respectively. Layer 6 had the most distinct M-types, broadly split into two categories: those with short or inverted apical dendrites (L6short), consistent with IT subclasses; and those with tall apical and narrow basal dendrites (L6tall), consistent with CT subclasses⁸. It was not possible to unambiguously label some layer 6 neurons as either IT or CT on the basis of anatomy alone, but 99% (n = 142 of 143) of manually assigned CT cells fell into one of the L6tall M-types.

Most M-types had visually distinguishable characteristics (Fig. 3d and Extended Data Fig. 2), but in some cases subtle differences in skeleton features were differentiated by stark differences in synaptic properties. For example, the two layer 2 M-types are visually similar, although L2a had a 29% higher overall dendritic length (L2a, 4,532 µm; L2b, 3,510 µm). However, L2a cells had 80% more synaptic inputs than L2b cells (L2a, 4,758; L2b, 2,649), a 40% higher median synapse density (L2a, 1.04 synapses per micrometre; L2b, 0.72 synapses per micrometre) (Fig. 3f,g) and a wider distribution of synapse sizes (Extended Data Fig. 1). Median synapse size turned out to differ across M-types, often matching layer transitions (Fig. 3g and Extended Data Fig. 1). Strikingly, L5 NP cells were outliers across synaptic properties, with the fewest total dendritic inputs, lowest synaptic input density and among the smallest synapses (Fig. 3g,h). Excitatory M-types thus differed not only in morphology but also in cell-level synaptic properties like total synaptic input and local properties like synapse size.

Subtype definitions based on structural properties may or may not be meaningful to cortical circuitry. If different M-types received input from different inhibitory populations, it would indicate potential for other circuit differences as well. Having classified inhibitory subclasses and excitatory M-types, we thus analysed how inhibition is distributed across the landscape of excitatory neurons.

The column reconstructions included 70,884 synapses from inhibitory neurons onto excitatory neurons (Fig. 4a). PeriTCs and DistTCs were by far the dominant source of inhibition, with individual cells having as many as 2,118 synapses onto excitatory cells in the column (mean PeriTC, 581 synapses per presynaptic cell; mean DistTC, 596 synapses per presynaptic cell), whereas SparTCs and InhTCs made far fewer synapses per presynaptic cell (mean SparTC, 74 synapses; mean InhTC, 16 synapses; Fig. 4b). Inhibition was distributed unequally across M-types (Fig. 4c). Much of this difference was related to differences in overall synaptic input. Across M-types, synaptic input at the soma, which is almost completely inhibitory, was strongly correlated (r = 0.96, P = 5 × 10⁻¹⁰) with net synaptic input onto dendrites, which is primarily excitatory (Extended Data Fig. 7). Notably, this structural balance of dendritic and somatic input also remained significant across individual cells in 16 of 18 M-types.

a, Connectivity from all inhibitory neurons (columns) onto all excitatory neurons (rows), sorted by M-type and soma depth. Dot size indicates net number of synapses observed. b, Net synapses onto column cells for each inhibitory subclass. Black dots indicate median; bars show 5% confidence interval. c, Mean net synapses per target cell from each inhibitory subclass onto each excitatory M-type. d, Spearman correlation of PeriTC and DistTC net input onto individual cells, measured in each M-type. Bars indicate 95% confidence interval based on bootstrapping (n = 2,000). Stars indicate M-types significantly different from zero with a P value < 0.05 after Holm–Sidak multiple test correction. e,f, Pearson correlation of connectivity density between excitatory M-types on the basis of PeriTCs (e) and DistTCs (f). Dotted lines indicate groups of cells roughly in a layer.

Source Data

Full size image

Similarly, synaptic input from PeriTC and DistTC was also typically balanced onto individual cells for each M-type. We examined the number of PeriTC and DistTC inputs onto individual excitatory neurons for each M-type and found significant positive correlation for 12 of 18 M-types (Fig. 4d), indicating coordinated amounts of inhibitory synaptic inputs across the entire arbour of target cells. M-types in upper layers had particularly heterogeneous amounts of inhibitory input, with L2b cells receiving 60% fewer synapses from intracolumnar interneurons as spatially intermingled L2a cells (L2b, 37.7 ± 0.27 synapses; L2a, 94.8 ± 0.58 synapses), whereas L3b cells had nearly as many intracolumnar inhibitory inputs as much larger L5 ET cells. All layer 6 M-types had relatively few intracolumnar inhibitory inputs compared to upper layers (Fig. 4c). However, note that the columnar sampling only reflects local sources of inhibition and does not capture the net effect of potentially wider or narrower spatial domains of inhibitory integration between layers.

Individual inhibitory neurons often targeted several M-types, indicating that certain combinations could be inhibited together. For each inhibitory neuron, we computed the connection density onto each M-type: that is, the fraction of cells in the column that received synaptic input from it (Fig. 4e). To measure the structure of co-inhibition, we computed the correlation of inhibitory connection density between M-types across PeriTCs and DistTCs separately (Fig. 4f). A high correlation would indicate that the same inhibitory neurons that connected more (or less) to one M-type also connect more (or less) to another, whereas zero correlation would indicate independent sources of inhibition between M-types.

These correlations showed several notable features of the structure of inhibition across layers. In superficial cortex, the layer 2 and layer 3 M-types were strongly correlated in the layer but had modest correlations between layers, indicating largely different sources of inhibition. Layer 4 M-types, in contrast, were all highly correlated with one another. Layer 5 M-types were more complex and suggested largely non-overlapping sources of inhibition, particularly among neurons with different long-range projection targets. Layer 6 inhibition was virtually independent from other layers, with DistTC connectivity also distinct between IT-like L6short cells and CT-like L6tall cells. Collectively, both layer and projection subclass were key factors in shaping co-inhibition. Importantly, most cotargeting relationships were consistent for both PeriTC and DistTC output, indicating that cardinal inhibitory subclasses distribute their output across excitatory neurons with similar patterns of connectivity.

How do individual neurons distribute their output to produce the patterns of inhibition described above? To compare patterns of output, for every inhibitory neuron, we measured the fraction of synaptic outputs made onto each M-type (Fig. 5a). This normalized synaptic output budget reflected factors such as the number of synapses per connection and the number of potential targets but was not strongly affected by partial arbours. We performed a consensus clustering (Methods), identifying 18 ‘motif groups’, sets of cells with similar patterns of output connectivity (Fig. 5a and Extended Data Fig. 8). Although this measurement only included synapses with cells in the column, interneurons made more than four times more synapses onto cells outside the column than within (Extended Data Fig. 8s,t). To check whether these results would hold with data outside the column, we used a prediction of neuronal M-types on the basis of perisomatic features and trained on column M-type labels⁴⁰. We found that within-column and predicted dataset-wide synaptic output budgets were highly correlated (Pearson r = 0.90), confirming that the columnar sampling provided a good estimate of overall neuronal connectivity (Extended Data Fig. 8u,v).

a, Distribution of synaptic output for all interneurons, clustered into motif groups with common target distributions. Each row is an excitatory target M-type, each column is an interneuron, and colour indicates fraction of observed synapses from the interneuron onto the target M-type. Only synapses onto excitatory neurons are used to compute the fraction. Neurons are ordered by motif group and soma depth. Bar plots along top indicate number of synapses onto column cells, with colour showing subclass (as in d). Bar plots along right indicate number of cells in target M-type. b,c, Morphology of all cells in group 4 (b) and group 13 (c), with colours as in a. Scale bar, 500 µm. d, Soma depth and subclass for cells in each motif group. Scale bar, 200 µm. e, Net synaptic output distribution across M-types for each motif group. f, Synaptic input for each M-type from each motif group as a fraction of all within-column inhibition. g, Schematic of motif group connectivity in upper layers. h, Schematic of motif group connectivity in Layer 5.

Source Data

Full size image

Each motif group represented a collection of cells that targeted the same pattern of excitatory cell types. Although some motif groups focused their output onto single excitatory M-types (such as group 9) or layers (such as group 7), others spanned broadly (such as group 6). However, motif groups were not simply individual cell types. Motif groups (Fig. 5b,c) showed diversity in both individual cell morphology and connectivity subclasses (Extended Data Figs. 3 and 8). Indeed, 15 of 18 groups (comprising 156 of 163 cells) included neurons from at least two subclasses, often aligned in cortical depth (Fig. 5d).

To summarize the relationship between motif groups and M-types, we computed both the average output fraction from each motif group onto each M- type (Fig. 5e) and the average input fraction of within-column inhibitory synaptic inputs onto a given M-type from each motif group (Fig. 5f). Input fraction often followed output fraction for particularly strong connections, but not always. For example, although group 3 more strongly targeted M-types in layer 3 than layer 2, it still contributed a substantial fraction of all inhibitory layer 2 input. In addition, we found that dominant connections for motif groups had both high connectivity density and several synapses per connection (Extended Data Fig. 9), properties that indicate a strong functional role in the circuit.

Inhibitory circuits were organized differently in upper layers compared to layers 5 and 6. In layers 2–4, each excitatory M-type received strong inhibition from 2–3 motif groups with overlapping combinations of targets, some specific in layers and others that cross layer boundaries (Fig. 5g). In contrast, most motif groups targeted only single M-types in layer 5, although in some cases they also targeted cells in other layers (Fig. 5h). Connectivity patterns in layer 6 included clear examples of IT-specific and putatively CT-specific cells, similar to layer 5 projection subclasses, but also had cells, particularly PeriTCs, that targeted widely in layer 6.

Cell type specificity, how concentrated the output onto targets, is typical among interneurons described here. However, specificity can arise in many ways. Different neurons have varying dendritic and axonal morphologies, synaptic densities and compartment preferences that constrain potential interactions^41,42. In addition, they can exhibit cell type selectivity, which we define as forming synapses with particular cell types more or less than might be expected on the basis of other factors such as axon/dendrite overlap.

To differentiate the effects of different contributing factors, for each interneuron we assembled information on morphology (Fig. 6a), synaptic connectivity and how its output is distributed across compartments (Fig. 6b) or excitatory M-types (Fig. 6c) in the column. We developed a selectivity index by comparing observed synaptic connectivity to a null model ignoring M-type but capturing compartment preference and postsynaptic factors such as number of synapses a cell typically receives and the spatial heterogeneity of potential targets⁴² (Fig. 6d). Because many of these factors are correlated with target cell properties, such a null model aimed to address confounding between the M-type label itself and those structural properties that affect connectivity more generally: for example, whether cells with higher input synapse density receive more inhibitory inputs irrespective of their M-type.

a, Example inhibitory neuron (cell ID 303085). Axon in blue, dendrite in red. Scale bar, 200 µm. b, Distribution of synaptic outputs across target compartments for the cell in a. c, Distribution of synaptic outputs across M-types (bar length) and compartments (bar colours) for the cell in a. d, Selectivity index values for the cell in a, measured as the ratio of observed synapse count to median shuffled synapse count for a null model as described below. Error bars indicate 95th percentile interval. Coloured dots (blue, low; orange, high) indicate significant differences (two-sided shuffle, P < 0.05) relative to the shuffle distribution after Holm–Sidak multiple test correction. e, As a baseline synapse distribution for null models, all synaptic inputs onto all cells in the column were binned by compartment, depth and M-type. (See Extended Data Fig. 10 for more details.) f, Shuffled connectivity for the cell in a was computed by sampling from the baseline synapse distribution with the observed depth and compartment bins and counting synapses onto each M-type across all bins (n = 1,000 shuffles). Example shuffle values for L3a (top) and L4a (bottom) M-types versus observed synapses are shown. g, Selectivity index for all cells in motif group 5. Non-significant values are assigned a value of 1. The cell in a is highlighted by a black box. h, Direction of the median cell’s selectivity index from each motif group onto each M-type. Orange indicates more connected, blue less connected. Connections where the median selectivity index was non-significant are indicated with a dot. i–l, Compact cell connectivity cards encapsulating anatomy (left), M-type target distribution (middle, bar length), compartment targeting (middle, bar colours as in d) and selectivity index (right, as in g) for four example neurons: an L5ET-specific basket cell (i), a deep-layer-specific upper layer neuron (j), a translaminar basket cell (k) and a translaminar layer 6 interneuron (l). Full connectivity cards for all cells can be found in Supplementary File 1. Scale bars, 200 µm. SI, selectivity index.

Source Data

Full size image

We computed a global baseline distribution of all synaptic inputs onto all column neuron dendrites, binned by cortical depth (20-µm-deep bins), M-type and target compartment (Fig. 6e and Extended Data Fig. 10). For each interneuron, we computed a shuffled output distribution across M-types by repeatedly sampling connectivity from the baseline distribution, matching the synapse depth and compartment targeting distributions for that cell’s outputs. For each connection from an interneuron onto a potential target M-type, we defined the selectivity index as the ratio of observed connectivity to the median of the distribution of shuffled connectivity (n = 10,000 repeats; Fig. 6f), reflecting the amount of cell-type-dependent selectivity beyond the factors included in the null model. Although this sampling included both excitatory and inhibitory synapses, previous studies^43,44 and our data indicate that excitatory and inhibitory inputs are proportional to one another, even at the level of individual cells (Extended Data Fig. 7).

Because motif groups had common specificity by definition, we asked if cells in motif groups had common patterns of selectivity as well (Fig. 6g). We computed the median selectivity index for each M-type/motif group pair, setting non-significant selectivity index values to 1 (Fig. 6h). We found that although 17 of 18 motif groups showed consistent positive or negative selectivity for some targets, in many cases highly specific connectivity was not associated with increased selectivity. For example, group 1 is highly specific to layer 2 targets but surprisingly did not show consistent positive selectivity for them (Extended Data Fig. 10). Examination of each constraint of the null model—synapse abundance of different targets, presynaptic compartment specificity and presynaptic depth distribution—suggested that this was because of group 1 axons having a narrow spatial distribution of axons that strongly overlapped layer 2 targets, which for many cells was sufficient to explain their connectivity (Extended Data Fig. 10). This effect was more pronounced for PeriTCs, which tended to target a more compact spatial ___domain with less overlap between M-types. In contrast, for DistTCs, the increased spatial overlap of distal and apical dendrites of different M-types required more selectivity to explain their connectivity (Extended Data Fig. 10). Collectively, this suggests that to achieve specific targeting patterns, interneurons both project their axons to precise spatial domains and selectively favour or disfavour making synapses with specific targets, with the relative contribution of these factors differing across cell types.

Here we generated a detailed map of neuronal structure and inhibitory connectivity in a column of visual cortex using EM. Using synaptic properties in addition to traditional morphological features, we found a collection of excitatory M-types with distinct patterns of inhibitory input, demonstrating that anatomical distinctions are reflected in the local inhibitory circuit. We use the term ‘motif groups’ to describe this organization of inhibitory neurons—a diverse collection of cells, extending beyond the concept of cell types, that target specific combinations of M-types’ perisomatic and dendritic compartments. The distribution of inhibitory motif groups also offers insights into the functional relationships of excitatory cell types. In layers 5 and 6, each projection subclass (IT, ET, NP and CT) had a collection of inhibitory cells for which they were the predominant target. This affords the network the potential to individually control each projection subclass via selective inhibition both at the soma and across dendrites, potentially with different inhibitory types active under different network conditions and behavioural states.

Cell connectivity cards

Although motif groups described the broad organization of groups of cells, individual interneurons showed fascinating but idiosyncratic structural properties. To concisely convey individual cell properties, we summarized morphology and connectivity into ‘connectivity cards’ (Fig. 6i–l). Individual cards can show unique features that were not clear from groups alone, such as extreme specificity (Fig. 6i) or different patterns of translaminar connectivity (Fig. 6j–l). An atlas for all interneurons can be found in Supplementary File 1.

Sublaminar inhibitory specificity

This question of inhibitory specificity has been perhaps best studied in layer 5, with its highly distinct ET and IT excitatory projection suclasses⁴⁵. For dendrite-targeting cells, precise genetic targeting of layer 5 SST subtypes identified distinct cell types that targeted ET versus IT cells²¹. In addition, developmental perturbation altering ET or IT neurons has suggested that they have different perisomatic input as well, with PV cells preferentially targeting ET and cholecystokinin neurons targeting IT^46,47,48. Here, ET cells received input from a larger number and diversity of inhibitory cells than layer 5 IT, despite being less numerous, indicating that as the primary subcortical output cell, ET neurons have a larger and more diverse inhibitory network than IT cells. ET cells were also frequently involved in translaminar circuits, with several examples of both ascending and descending translaminar PeriTCs and ascending DistTCs that targeted both L2/3 neurons and ETs but not ITs, indicating bidirectional pathways for coordinated inhibition. In addition, layer 5 IT and NP cells had distinct collections of inhibitory neurons. Projection-specific inhibition was also found in layer 6 between putative IT and CT neurons. In contrast to layer 5, however, there was a combination of both projection-specific and broad layer 6 inhibition. These connectivity patterns afford the cortical network the potential to selectively inhibit distinct projection subclasses, potentially with different cell types active under different network states or using different plasticity rules.

Even in layers 2–4, with only IT cells, there was significant sublaminar specificity. The differential inhibition of layer 2 versus layer 3 cells suggests that are functionally distinct subnetworks with independent modulation. This could mirror depth-dependent differences in intracortical projection patterns⁴⁹, similar to prefrontal cortex, where amygdala-projecting layer 2 cells receive inhibition that selectively avoids neighbouring cortical-projecting cells^20,38. Another possibility is that they are well posed to differentially modulate top-down versus sensory-driven activity^50,51, as layer 3 receives more sensory thalamic input than layer 2 (ref. ¹¹). More generally, the distinct inhibitory environments of upper and lower layer 2/3 have been observed across cortex, from primary sensory areas⁵² to higher-order association areas⁵³, indicating that they may reflect a more general functional specialization.

Basket cell disinhibition

Many VIP interneurons preferentially inhibit other inhibitory neurons. Such disinhibitory VIP neurons have been shown to strongly target SST neurons across cortical areas^34,35 and, to a lesser extent, fast spiking or PV neurons^35,54. Here we found two classes of disinhibitory specialists with distinct and specific targets: one preferring putative SST cells and one preferring basket cells.

Fast spiking or PV basket cells are inhibited by many sources, including other PV cells^23,34, SST cells⁵⁵ and even neurogliaform cells⁵⁶. However, the basket-targeting disinhibitory specialists differ from these other pathways in their specificity—not only do they distribute most synaptic output onto basket cells instead of any other inhibitory or excitatory targets but they do so with larger synapses and more synapses per connection. This highly specific targeting offers an intriguing pathway to control basket-cell-mediated excitatory gain or synchrony without significantly affecting other neuronal populations. Determining what conditions cause these cells to be active will be important for understanding their functional effect. Future experiments will also be required to determine their molecular subclass.

Limitations

The principal concern is the generalizability of data, because it comes from a single animal, in one ___location near the edge of VISp, and has at most a few examples per cell type. However, companion work from the same dataset focusing on several examples of morphologically defined cell types shows consistent target preferences^31,57, and our data also agree with recent functional measurements of type-specific connectivity of SST cells²¹. We thus expect that the broad connectivity results will apply generally, although it will be important to measure the variability across individual cells, distinct animals and locations in cortex.

This study also only considered cells and connectivity in a narrow range of distances and limited volume. If cells change their connectivity with distance, as has been seen in excitatory neurons⁵⁸, this would bias the observed connectivity distributions. Extending a similar analysis across a much wider extent will be important for building a complete map of inhibitory cell types and more firmly establishing the nature of inhibitory motif groups.

Multimodal cell typing with EM

The M-types found here from morphology and synaptic properties generally agree with approaches from morphology alone⁵⁹ or in combination with other modalities^8,9, in particular distinguishing cells in upper layer and lower layer 2/3 and differentiating between projection subclasses. Sublaminar variation is also found in transcriptomic studies, with several excitatory clusters in upper layer 2/3 in VISp^60,61 and variation in other layers, although it is not clear if they correspond precisely to the M-types observed here.

To facilitate subsequent analysis of anatomy, connectivity and ultrastructure, all EM data, segmentations, skeletons and tables of synapses and cell types are available are available via the MICrONS-Explorer website². However, making the best experimental use of EM data will require linking EM to genetic tools. Patch-seq, which generates combined electrophysiology, transcriptomic and morphological data, was used in a companion study to quantitatively link particular Martinotti cells from EM to specific transcriptomic subtypes³¹. At present, however, transcriptomic clusters often have diverse morphologies and probably diverse connectivity^9,16. This suggests that the process of linking structural and molecular datasets should aim to become bidirectional, not only decorating EM reconstructions with transcriptomic information but also using EM to identify cell types with distinct connectivity and analysing Patch-seq data to identify distinguishing transcriptomic markers or collecting more examples.

This dataset was acquired, aligned and segmented as part of the larger MICrONS project. Methods underlying dataset acquisition are described in full detail elsewhere^2,62,63,64, and the primary data resource is described in a separate publication². We repeat some of the methodological details for the dataset here for convenience.

Animal preparation for EM

All animal procedures were approved by the Institutional Animal Care and Use Committee at the Allen Institute for Brain Science or Baylor College of Medicine. Neurophysiology data acquisition was conducted at Baylor College of Medicine before EM imaging. Afterwards the mice were transferred to the Allen Institute in Seattle and kept in a quarantine facility for 1–3 days, after which they were euthanized and perfused. All results described here are from a single male mouse, age 64 days at onset of experiments, expressing GCaMP6s in excitatory neurons via SLC17a7-Cre and Ai162 heterozygous transgenic lines (recommended and generously shared by Hongkui Zeng at the Allen Institute for Brain Science; JAX stock 023527 and 031562, respectively). Two-photon functional imaging took place between P75 and P80 followed by two-photon structural imaging of cell bodies and blood vessels at P80. The mouse was perfused at P87.

Tissue preparation

After optical imaging at Baylor College of Medicine, candidate mice were shipped via overnight air freight to the Allen Institute. Mice were transcardially perfused with a fixative mixture of 2.5% paraformaldehyde, 1.25% glutaraldehyde and 2 mM calcium chloride, in 0.08 M sodium cacodylate buffer, pH 7.4. A thick (1,200 µm) slice was cut with a vibratome and post-fixed in perfusate solution for 12–48 h. Slices were extensively washed and prepared for reduced osmium treatment based on the protocol of ref. ⁶⁵. All steps were performed at room temperature, unless indicated otherwise. The first osmication step involved 2% osmium tetroxide (78 mM) with 8% v/v formamide (1.77 M) in 0.1 M sodium cacodylate buffer, pH 7.4, for 180 min. Potassium ferricyanide 2.5% (76 mM) in 0.1 M sodium cacodylate, 90 min, was then used to reduce the osmium. The second osmium step was at a concentration of 2% in 0.1 M sodium cacodylate for 150 min. Samples were washed with water and then immersed in thiocarbohydrazide for further intensification of the staining (1% thiocarbohydrazide (94 mM) in water, 40 °C, for 50 min). After washing with water, samples were immersed in a third osmium immersion of 2% in water for 90 min. After extensive washing in water, Walton’s lead aspartate (20 mM lead nitrate in 30 mM aspartate buffer, pH 5.5, 50 °C, 120 min) was used to enhance contrast. After two rounds of water wash steps, samples proceeded through a graded ethanol dehydration series (50%, 70%, 90% w/v in water, 30 min each at 4 °C, then 3 × 100%, 30 min each at room temperature). Two rounds of 100% acetonitrile (30 min each) served as a transitional solvent step before proceeding to epoxy resin (EMS Hard Plus). A progressive resin infiltration series (1:2 resin:acetonitrile (for eample, 33% v/v), 1:1 resin:acetonitrile (50% v/v), 2:1 resin acetonitrile (66% v/v) and then 2 × 100% resin, each step for 24 h or more, on a gyrotary shaker), was done before final embedding in 100% resin in small coffin moulds. Epoxy was cured at 60 °C for 96 h before unmoulding and mounting on microtome sample stubs. The sections were then collected at a nominal thickness of 40 nm using a modified ATUMtome (RMC/Boeckeler⁶²) onto six reels of grid tape^62,66.

Transmission EM imaging

The parallel imaging pipeline used in this study⁶² used a fleet of transmission electron microscopes that had been converted to continuous automated operation. It was built on a standard JEOL 1200EXII 120 kV transmission electron microscope that had been modified with customized hardware and software, including an extended column and a custom electron-sensitive scintillator. A single large-format CMOS (complementary metal–oxide–semiconductor) camera outfitted with a low-distortion lens was used to grab image frames at an average speed of 100 ms. The autoTEM was also equipped with a nano-positioning sample stage that offered fast, high-fidelity montaging of large tissue sections and a reel-to-reel tape translation system that locates each section using index barcodes. During imaging, the reel-to-reel GridStage moved the tape and located the targeting aperture through its barcode and acquired a 2D montage. We performed quality control on all image data and reimaged sections that failed the screening.

Image processing

Volume assembly

The volume assembly pipeline is described in detail elsewhere^63,64. Briefly, the images collected by the autoTEMs are first corrected for lens distortion effects using nonlinear transformations computed from a set of 10 × 10 highly overlapping images collected at regular intervals. Overlapping image pairs are identified in each section, and point correspondences are extracted using features extracted using the scale-invariant feature transform. Montage transformation parameters are estimated per image to minimize the sum of squared distances between the point correspondences between these tile images, with regularization. A downsampled version of these stitched sections is produced for estimating a per-section transformation that roughly aligns these sections in three dimensions. The rough aligned volume is rendered to disk for further fine alignment. The software tool used to stitch and align the dataset is available on GitHub (https://github.com/AllenInstitute/render-modules). To fine align the volume, we needed to make the image processing pipeline robust to image and sample artefacts. Cracks larger than 30 um (in 34 sections) were corrected by manually defining transforms. The smaller and more numerous cracks and folds in the dataset were automatically identified using convolutional networks trained on manually labelled samples using 64 × 64 × 40 nm³ resolution images. The same was done to identify voxels containing tissue. The rough alignment was iteratively refined in a coarse-to-fine hierarchy⁶⁷ using an approach based on a convolutional network to estimate displacements between a pair of images⁶⁸. Displacement fields were estimated between pairs of neighbouring sections and then combined to produce a final displacement field for each image to further transform the image stack. Alignment was refined first using 1,024 × 1,024 × 40 nm³ images and then 64 × 64 × 40 nm³ images. The composite image of the partial sections was created using the tissue mask previously computed.

Segmentation

The image segmentation pipeline is fully described in ref. ⁶³. Remaining misalignments were detected by cross-correlating patches of image in the same ___location between two sections after transforming into the frequency ___domain and applying a high-pass filter. Combining with the tissue map previously computed, a ‘segmentation output mask’ was generated that sets the output of later processing steps to zero in locations with poor alignment. Using previously described methods⁶⁹, a convolutional network was trained to estimate intervoxel affinities that represent the potential for neuronal boundaries between adjacent image voxels. A convolutional network was also trained to perform a semantic segmentation of the image for neurite classifications, including (1) soma + nucleus, (2) axon, (3) dendrite, (4) glia and (5) blood vessel. Following the methods described in ref. ⁷⁰, both networks were applied to the entire dataset at 8 × 8 × 40 nm³ in overlapping chunks to produce a consistent prediction of the affinity and neurite classification maps, and the segmentation output mask was applied to predictions. The affinity map was processed with a distributed watershed and clustering algorithm to produce an oversegmented image, where the watershed domains are agglomerated using single-linkage clustering with size thresholds^71,72. The over-segmentation was then processed by a distributed mean affinity clustering algorithm^71,72 to create the final segmentation.

For synapse detection and assignment, a convolutional network was trained to predict whether a given voxel participated in a synaptic cleft. Inference on the entire dataset was processed using the methods described in ref. ⁷⁰ using 8 × 8 × 40 nm³ images. These synaptic cleft predictions were segmented using connected components, and components smaller than 40 voxels were removed. A separate network was trained to perform synaptic partner assignment by predicting the voxels of the synaptic partners given the synaptic cleft as an attentional signal⁷³. This assignment network was run for each detected cleft, and coordinates of both the presynaptic and postsynaptic partner predictions were logged along with each cleft prediction.

For nucleus detection², a convolutional network was trained to predict whether a voxel participated in a cell nucleus. Following the methods described in ref. ⁷⁰, a nucleus prediction map was produced on the entire dataset at 64 × 64 × 40 nm³.

Column description and cell classes

The column borders were found by manually identifying a region in primary visual cortex that was far from both dataset boundaries and the boundaries with higher-order visual areas. A 100 × 100 µm box was placed on the basis of layer 2/3 and was extended along the y axis of the dataset.

While analysing data, we observed that deep layer neurons had apical dendrites that were not oriented along the most direct pia-to-white-matter direction, and we adapted the definition of the column to accommodate these curved neuronal streamlines. Using a collection of layer 5 ET cells, we placed points along the apical dendrite to the cell body and then along the primary descending axon towards white matter. We computed the slant angle as two piecewise linear segments, one along the negative y axis to lower layer 5 where little slant was observed, and one along the direction defined by the vector-averaged direction of the labelled axons. We believe the slant to be a biological feature of the tissue and not a technical artefact for several reasons:

1. 1.

   The curvature is not aligned to a sectioning plane or associated with shearing or other distortion in the imagery, making it unlikely to be a result of the alignment process.

2. 2.

   Blood vessel segmentation does not show a large, correlated distortion in deep layers, making it unlikely to be a result of mechanical stress on the tissue (https://ngl.microns-explorer.org/#!gs://microns-static-links/mm3/blood_vessels.json). Moreover, it is unclear why such stress would affect only layer 5b and below.

3. 3.

   Individual examples of neurons with slanted morphologies can be found among single-cell reconstructions in the literature: for example, several descending bipolar VIP interneurons and layer 6 pyramidal cells in ref. ¹⁶. It is not possible to determine whether these individual cases correspond to a larger population of correlated arbours, but it suggests these morphologies are not atypical.

4. 4.

   Similar curvature has been observed in other large EM datasets from visual cortex (data not shown) and light level morphological reconstructions, particularly among layer 6 pyramidal cells.

Using these boundaries and previously computed nucleus centroids², we identified all cells in the columnar volume. Coarse cell classes (excitatory, inhibitory and non-neuronal) were assigned on the basis of brief manual examination and rechecked by subsequent proofreading and automated cell typing⁴⁰. To facilitate concurrent analysis and proofreading, we split all false merges connecting any column neurons to other cells (as defined by detected nuclei) before continuing with other work.

Proofreading

Proofreading was performed primarily by five expert neuroanatomists using the CAVE infrastructure⁷⁴ and a modified version of Neuroglancer⁷⁵. Proofreading was aided by on-demand highlighting of branch points and tips on user-defined regions of a neuron on the basis of rapid skeletonization (https://github.com/AllenInstitute/Guidebook). This approach quickly directed proofreader attention to potential false merges and locations for extension, as well as allowed a clear record of regions of an arbour that had been evaluated.

For dendrites, we checked all branch points for correctness and all tips to see if they could be extended. False merges of simple axon fragments onto dendrites were often not corrected in the raw data because they could be computationally filtered for analysis after skeletonization (see below). Detached spine heads were not comprehensively proofread, and previous estimates place the rate of detachment at roughly 10–15%. Using this method, dendrites could be proofread in about 10 min per cell.

For inhibitory axons, we began by ‘cleaning’ axons of false merges by looking at all branch points. We then performed extension of axonal tips until either their biological completion or data ambiguities, particularly emphasizing all thick branches or tips that were well-suited to project to new laminar regions. For axons with many thousands of synaptic outputs, we followed many but not all tips to completion once primary branches were cleaned and established. For smaller neurons, particularly those with bipolar or multipolar morphology, most tips were extended to the point of completion or ambiguity. Axon proofreading time differed significantly by cell type not only because of differential total axon length but also because of axon thickness differences that resulted in differential quality of autosegmentations, with thicker axons being of higher initial quality. Typically, inhibitory axon cleaning and extension took 3–10 h per neuron.

The lack of segmentation in the top 10 µm of layer 1 truncates some apical tufts and limited reconstruction quality of layer 1 interneurons. For those excitatory neurons with extensive apical tufts, particularly layer 2 and L5ET cells, the reconstructions here might miss both distinguishing characteristics and sources of inhibitory input in that region. Similarly, axons in deep layer 6 were generally less complete because of alignment quality in white matter.

Manual cell subclass and layer labels

Expert neuroanatomists further labelled excitatory and inhibitory neurons into subclasses. Layer definitions were based on considerations of both cell body density (in analogy with nuclear staining) supplemented by identifying kinks in the depth distribution of nucleus size near expected layer boundaries⁴⁰.

For excitatory neurons, the categories used were Layer 2/3-IT, Layer 4-IT, Layer 5-IT, Layer 5-ET, Layer 5-NP, Layer 6-IT, Layer 6-CT and Layer 6b (‘L6-WM’) cells. Excitatory expert labels did not affect analysis but were used as the basis for naming morphological clusters. Layer 2/3 and upper Layer 4 cells were defined on the basis of dendritic morphology and cell body depth. Layer 5 cells were similarly defined by cell body depth, with projection subclasses distinguished by dendritic morphology following ref. ⁸ and classical descriptions of thick (ET) and thin-tufted (IT) cells. Layer 5 ET cells had thick apical dendrites, large cell bodies, numerous spines and a pronounced apical tuft, and deeper ET cells had many oblique dendrites. Layer 5 IT cells had more slender apical dendrites and smaller tufts, fewer spines and fewer dendritic branches overall. Layer 5 NP cells corresponded to the ‘Spiny 10’ subclass described in ref. ⁸; these cells had few basal dendritic branches, each very long and with few spines or intermediate branch points. Layer 6 neurons were defined by cell body depth, and some cells were able to be further labelled as IT or CT by human experts. Layer 6 pyramidal cells with stellate dendritic morphology, inverted apical dendrites or wide dendritic arbours were classified as IT cells. Layer 6 pyramidal cells with small and narrow basal dendrites, an apical dendrite ascending to Layer 4 or Layer 1 and a myelinated primary axon projecting into white matter were labelled as CT cells.

For inhibitory neurons, manual cell typing considered axonal and dendritic morphology as well as connectivity. Cells that primarily contacted soma or perisomatic regions were labelled as basket cells. Cells that made arbours that extended up to layer 1 or formed a dense plexus and primarily targeted distal dendrites were labelled as putative SST cells. Cells that remained mostly in layer 1 or had extensive arbourization and many non-synaptic boutons were labelled as putative Id2 or neurogliaform cells. Finally, cells with a bipolar dendritic morphology or a multipolar dendritic morphology and output onto other inhibitory neurons were labelled as putative VIP cells. Several cells, particularly in layer 6, had an ambiguous subclass assignment, typically when their connectivity was not basket-like but their morphology was also not similar to upper layer Martinotti or non-Martinotti cells.

Skeletonization

To rapidly skeletonize dynamic data, we took advantage of the PyChunkedGraph data structure that collects all supervoxels belonging to the same neuronal segmentation into 2 × 2 × 20 µm ‘chunks’ with a unique ID and precisely defined topological adjacency with neighbouring chunks of the same object. Each chunk is called a ‘level 2 chunk’, and the complete set of chunks for a neuron and their adjacency we call the ‘level 2 graph’ on the basis of its ___location in the hierarchy of the PyChunkedGraph data structure²⁷. We precompute and cache a representative central point in space and the volume and the surface area for each level 2 chunk and update this data when new chunks are created because of proofreading edits. Using the level 2 graph and assigning edge lengths corresponding to the distance between the representative points for each vertex (that is, each level 2 chunk), we run the TEASAR⁷⁶ algorithm (10 µm invalidation radius) to extract a loop-free skeleton. Each of the level 2 vertices removed by the TEASAR algorithm is associated with its closest remaining skeleton, making it possible to map surface area and volume data to the skeleton. Typical edges between skeleton vertices are about 1.7 µm, and new skeletons can be computed de novo in about 10 s, making them useful for analysis over length scales of tens of micrometre or larger.

To represent the cell body, a further vertex was placed at the ___location of the nucleus centroid, and all vertices within an initial radius and topologically connected to centroid were collapsed into this vertex with associated data mapping. The radius was determined for each neuron separately by consideration of the volume of each cell body. A companion work⁴⁰ computed the volume of each cell body, and we generated an effective radius on the basis of the sphere with the same volume. To ensure that our values captured potentially lopsided cell bodies, we padded this effective radius by a further factor of 1.25. Skeletons were rooted at the cell body, with ‘downstream’ meaning away from soma and ‘upstream’ meaning towards soma. Each synapse was assigned to skeleton vertices on the basis of the level 2 chunk of its associated supervoxel. For each unbranched segment of the skeleton (that is, between two branch points or between a branch point and end point), we computed an approximate radius r on the basis of a cylinder with the same path length L and total volume V associated with that segment: \((r=\sqrt{V/{\rm{\pi }}L})\).

Axon/dendrite classification

To detect axons, we took advantage of the skeleton morphology, the ___location of presynaptic and postsynaptic synapses and the clear segregation between inputs and outputs of cortical neurons. For inhibitory cells, we used synapse flow centrality⁷⁷ to identify the start of the axon as the ___location of maximum paths along the skeleton between sites of synaptic input and output. Two inhibitory neurons had two distinct, biologically correct axons after proofreading (cell IDs 258362 and 307059). For these cells, we ran this method twice, masking off the axon found after the first run, to identify both. For excitatory neurons that did not have extended axons, there were often insufficient synaptic outputs on their axon for this approach to be reliable. Excitatory neurons with a segregation index⁷⁷ of 0.7 (on a scale with 0 indicating random distribution of input and output synapses and 1 indicating perfect input/output segregation) or above were considered well-separated, and the synapse flow centrality solution was used. For cells with a segregation index less than 0.7, we instead looked for branches near the soma with few synaptic inputs. Specifically, we took identified all skeleton vertices within 30 µm of the cell body and looked at the distinct branches downstream from this region. For each branch, we computed the total path length and the total number of synaptic inputs to get a linear input density. Branches with both a path length more than 20 µm and an input density less than 0.1 synaptic inputs per micrometre were labelled as being axonal and filtered out of subsequent analysis.

We further filtered out any remaining axon fragments merged onto pyramidal cell dendrites using a similar approach. We identified all unbranched segments (regions between two branch points or between a branch point and end point) on the non-axonal region of the skeleton and computed their input synapse density. Starting from terminal segments (that is, those with no downstream segments), we labelled a segment as a ‘false merge’ if it had an input density less than 0.1 synaptic inputs per micrometre. This process iterated across terminal segments until all remaining had an input density of at least 0.1 inputs per micrometre. Falsely merged segments were masked out of the skeleton for all analysis.

Excitatory dendrite compartments

We assigned all synaptic inputs onto excitatory neurons to one of four compartments: soma, proximal dendrite, distal basal dendrite and distal apical dendrite. The most complex part was distinguishing the basal dendrite from the apical dendrite. Although easy in most cases for neurons in layer 3–5 because of the consistent nature of apical dendrites being single branches reaching towards layer 1, this is not true everywhere. In upper layer 2/3, cells often have several branches in layer 1 equally consistent with apical dendrites, and in layer 6 there are often cells with apical dendrites that stop in layer 4 and that point towards white matter or even that lack a clear apical branch entirely. To objectively and scalably define apical dendrites, we built a classifier that could detect between zero and three distinct apical branches per cell. Following the intuition from neuroanatomical experts, we used features on the basis of the branch orientation, ___location in space, relative ___location compared to the cell body and branch-level complexity. Specifically, we trained a random forest classifier to predict whether a skeleton vertex belonged to an apical dendrite on the basis of several features: depth of vertex, depth of soma, difference in depth between soma and vertex, vertex distance to soma along the skeleton, vertex distance to farthest tip, normalized vertex distance to tip (between 0 and 1), tortuosity of path to root, number of branch points along the path to root, radial distance from soma, absolute distance from soma and angle relative to vertical between the vector from soma to vertex. We aggregated predictions in each branch by summing the log-odds ratio from the model prediction, with the net log-odds ratio saturating at ±200. Finally, for each branch i with aggregated odds ratio R_i, we compare branches to one another via a soft-max operation: \({S}_{i}=\exp \left({R}_{i}/50\right)/\,{\sum }_{j}\exp ({R}_{j}/50)\). Branches with a maximum tip length of less than 50 µm were considered too short to be a potential apical dendrite and excluded from consideration and not included in the denominator. Branches with both R_i > 0 (evidence is positive towards being apical) and S_i > 0.25 were defined to be apical. Note that the soft-max was chosen to allow multiple apical branches if they had similar aggregated odds ratios, which was found to be necessary for upper layer pyramidal neurons. Training data were selected from an initial 50 random cells, followed by a further 33 cells chosen representing cases where the classifier did not perform correctly. Performance on both random and difficult cells had an F1-score of 0.9297 (86 true positives, 599 true negatives, 2 false positives and 11 false negatives) on the basis of leave-one-out cross validation, with at least one apical dendrite correctly classified for all cells.

Compartment labels were propagated to synapses on the basis of the associated skeleton vertices. Soma synapses were all those associated with level 2 chunks in the soma collapse region (see the section ‘Skeletonization’). Proximal dendrites were those outside of the soma but within 50 µm after the start of the branch. Distal basal synapses were all those associated with vertices more distant than the proximal threshold but not on an apical branch. Apical synapses were all those associated with vertices more distant than the proximal threshold and on an apical branch.

Inhibitory feature extraction and clustering

Many classical methods of distinguishing interneuron classes are based on how cells distribute their synapses across target compartments. Following proofreading, expert neuroanatomists attempted to classify all inhibitory neurons broadly as ‘basket cells’, ‘SST-like cells’, ‘VIP-like cells’ and ‘neurogliaform/layer 1’ cells on the basis of connectivity properties and morphology. Although 150 cells were labelled on this basis, a further 13 neurons were considered uncertain (primarily in layer 6), and in some cases manual labels were low confidence. To classify inhibitory neurons in a data-driven manner, we thus measured four properties of how cells distribute their synaptic outputs:

1. 1.

   The fraction of synapses onto inhibitory neurons.

2. 2.

   The fraction of synapses onto excitatory neurons that are onto soma.

3. 3.

   The fraction of synapses onto excitatory neurons that are onto proximal dendrites.

4. 4.

   The fraction of synapses onto excitatory neurons that are onto distal apical dendrites.

   Because the fraction of synapses targeting all compartments sums to one, the last remaining property, synapses onto distal basal dendrites, was not independent and thus was measured but not included as a feature. Inspection of the data suggested two more properties that characterized synaptic output across inhibitory neurons:

5. 5.

   The fraction of synapses that are part of multisynaptic connections, those with at least two synapses between the same presynaptic neuron and target neuron.

6. 6.

   The fraction of multisynaptic connection synapses that were also within 15 µm of another synapse with the same target, as measured between skeleton nodes. Note that we evaluated the robustness of this parameter and found that intersynapse distances from 5 to more than 100 µm have qualitatively similar results (Extended Data Fig. 2).

Using these six features, we trained a linear discriminant classifier on cells with manual annotations and applied it to all inhibitory cells. Differences from manual annotations were treated not as inaccurate classifications but rather as a different view of the data.

Excitatory feature extraction and clustering

To characterize excitatory neuron morphology, we computed features based only on excitatory neuron dendrites and soma. The features were as follows:

1. 1.

   Median distance from branch tips to soma per cell.

2. 2.

   Median tortuosity of the path from branch tips to soma per cell. Tortuosity is measured as the ratio of path length to the Euclidean distance from tip to soma centroid.

3. 3.

   Number of synaptic inputs on the dendrite.

4. 4.

   Number of synaptic inputs on the soma.

5. 5.

   Net path length across all dendritic branches.

6. 6.

   Radial extent of dendritic arbour. We define ‘radial distance’ to be the distance in the same plane as the pial surface. For every neuron, we computed a pia-to-white-matter line, including slanted region in deep layers, passing through its cell body. For each skeleton vertex, we computed the radial distance to the pia-to-white-matter line at the same depth. To avoid any outliers, the radial extent of the neuron was defined to be the 97th percentile distance across all vertices.

7. 7.

   Median distance to soma across all synaptic inputs.

8. 8.

   Median synapse size of synaptic inputs onto the soma.

9. 9.

   Median synapse size of synaptic inputs onto the dendrites.

10. 10.

    Dynamic range of synapse size of dendrite synaptic inputs. This was measured as the difference between 95th and fifth percentile synapse sizes.

11. 11.

    Shallowest extent of synapses, on the basis of the fifth percentile of synapse depths.

12. 12.

    Deepest extent of synapses, on the basis of the 95th percentile of synapse depths.

13. 13.

    Vertical extent of synapses, on the basis of the difference between 95th and fifth percentile of synapse depths.

14. 14.

    Median linear density of synapses. This was measured by computing the net path length and number of synapses along 50 depth bins from layer 1 to white matter and computing the median. A linear density was found by dividing synapse count by path length per bin, and the median was found across all bins with non-zero path length.

15. 15.

    Median radius across dendritic skeleton vertices. To avoid the region immediately around the soma from having a potential outlier effect, we only considered skeleton vertices at least 30 µm from the soma.

Three more sets of features used component decompositions. To more fully characterize the absolute depth distribution of synaptic inputs, for each excitatory neuron, we computed the number of synapses in each of 50 depth bins from the top of layer 1 to surface of white matter (bin width approximately 20 µm). We Z-scored synapse counts for each cell and computed the top six components using SparsePCA. The loadings for each of these components on the basis of the net synapse distribution were used as features.

To characterize the distribution of synaptic inputs relative to the cell body instead of cortical space, we computed the number of synapses in 13 soma-adjusted depth bins starting 100 µm above and below the soma. As before, synapse counts were Z-scored, and we computed the top five components using SparsePCA. The loadings for each of these components were used as further features.

To characterize the relationship with branching to distance, we measured the number of distinct branches as a function of distance from the soma at ten distances, every 30 µm starting at 30 µm from the soma and continuing to 300 µm. For robustness relative to precise branch point locations, the number of branches were computed by finding the number of distinct connected components of the skeleton found in the subgraph formed by the collection of vertices between each distance value and 10 µm towards the soma. We computed the top three singular value components of the matrix on the basis of branch count versus distance for all excitatory neurons, and the loadings were used as features.

All features were computed after a rigid rotation of 5 degrees to flatten the pial surface and translation to set the pial surface to 0 on the y axis. Features on the basis of apical classification were not explicitly used to avoid ambiguities on the basis of both biology and classification.

Using this collection of features, we clustered excitatory neurons by running phenograph⁷⁸ 500 times with 95% of cells included each time. Phenograph finds a nearest-neighbourhood graph on the basis of proximity in the feature space and clusters by running the Leiden algorithm for community detection on the graph. Here we used a graph on the basis of ten nearest neighbours and clustered with a resolution parameter of 1.3. These values were chosen to consistently separate layer 5 ET, IT and NP cells from one another, a well-established biological distinction. A coclustering matrix was assembled with each element corresponding to the number of times two cells were placed in the same cluster. To compute the final consensus clusters, we performed agglomerative clustering with complete linkage on the basis of the coclustering matrix, with the target number of clusters set by a minimum Davies–Bouldin score and a maximum Silhouette score. Clusters were then named on the basis of the most frequent manually defined cell type in the cluster and reordered on the basis of median soma depth. The labelling of cells as layer 2 and layer 3 was formed on the basis of soma depth and a morphology with a relatively flat morphology, often with no distinct apical trunk, although often apical-tuft-like branches emitted directly from the cell body. The L2c subclass was ambiguously defined between the two categories, with cells that had a distinct apical trunk but with connectivity and other properties that seemed more similar to layer 2 subclasses.

To compute the importance of each feature for each M-type, for each M-type we trained a random forest classifier to predict whether a cell belonged to it using scikit-learn⁷⁹. Because the classes were strongly imbalanced, we used SMOTE resampling to oversample datapoints from the smaller class. We used the Mean Decrease in Impurity metric, which quantifies how often a given feature was used in the decision tree ensemble.

Inhibitory connectivity and selectivity

To measure intracolumnar inhibitory connectivity, we first restricted synaptic outputs to the axon of each inhibitory neuron, as we have not observed any correctly classified synaptic outputs on dendritic arbours in this dataset. One cell with fewer than 30 synaptic outputs was omitted because of insufficient size. All remaining synaptic outputs across all interneurons were then filtered to include only those that target cells in the column, unless otherwise specified. Each output synapse was also labelled with the target skeleton vertex, dendritic compartment and M-type of the target neuron on the basis of the compartment definitions above.

For the inhibitory motif group clustering, for each interneuron we first computed the number of synapses across each excitatory M-types in the column. This synaptic output budget was then normalized per cell to generate a vector for each neuron with elements ranging from zero to one. Normalized synaptic output budgets were oversegmented using k-means (k = 20) with Euclidean distances 500 times, and a matrix of coclustering frequency—that is, the number of times two cells were put in the same k-means cluster—between individual cells was computed. Final M-types were found through agglomerative clustering with complete linkage of the coclustering matrix, scanning from two to 25 output clusters and selecting a final value of 18 on the basis of silhouette score and Davies–Bouldin score.

For measuring the synaptic output budget across cell types across the dataset (that is, inside and outside the column), we used a hierarchical classifier on the basis of a collection of perisomatic features that was trained on the data-driven clustering from the column sample⁴⁰. Only synapses onto object segmentation associated with a single nucleus and a cell type classification were used. Although most of these other targets were not proofread, estimates on the basis of proofread neurons indicate that 99% of non-proofread input synapses are accurate⁴⁰.

To measure inhibitory selectivity in the column, we compared the M-type distribution of its synaptic outputs to the M-type distribution of synaptic inputs according to a null model accounting for cell abundance, synapse abundance and depth. We first generated a baseline distribution of all 4,504,935 somatic or dendritic synaptic inputs to all column cells, where each synapse was associated with a precise depth, target compartment and M-type. We discretized synaptic inputs into 50 depth bins spanning pia to white matter, each covering approximately 20 µm and each of the five compartments: soma, proximal dendrite, basal dendrite, apical dendrite or inhibitory neuron. For each interneuron, we similarly discretized its synaptic output into the same bins, compartments and M-types. To generate a randomized output distribution preserving both observed depth and compartment distributions, we randomly picked synapses from the baseline distribution with the observed depth bins and compartment targets but without regard to M-type. We computed 10,000 randomized distributions per interneuron. To get a selectivity index, we compared the observed number of synapses onto a given M-type to the median of the number of synapses from the shuffle distribution. To get a significance for the selectivity index for a given M-type, we directly computed the two-sided P value of the observed number of synapses relative to the shuffle distribution for that M-type. P values were corrected for several comparisons using the Holm–Sidak method in each interneuron for those M-types with non-zero potential connectivity. Selectivity was only measured in the column because we did not generate compartment labels for unproofread dendrites outside of the column.

On connectivity cards, we also show a similar selectivity index on the basis of compartment rather than M-type. In that case, the shuffled distribution preserves observed depth and M-type output distributions but not compartments.

Reporting summary

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