## Brain anatomy revealed in startling detail

The mammalian cerebral cortex is an enormously complex network of neuronal processes that are long and thin, branching, and extremely densely packed. This high packing density has made the reconstruction of cortical neuronal networks challenging. Motta *et al.* used advanced automated imaging and analysis tools to reconstruct with high spatial resolution the morphological features of 89 neurons and their connections in the mouse barrel cortex. The reconstruction covered an area more than two orders of magnitude larger than earlier neuroanatomical mapping attempts. This approach revealed information about the connectivity of inhibitory and excitatory synapses of corticocortical as well as excitatory thalamocortical connections.

*Science*, this issue p. eaay3134

## Structured Abstract

### INTRODUCTION

The brain of mammals consists of an enormously dense network of neuronal wires: the axons and dendrites of nerve cells. Their packing density is so high that light-based imaging methods have so far only been able to resolve a very small fraction of nerve cells and their interaction sites, the synapses, in mammalian cortex. Recent advances in three-dimensional (3D) electron microscopy allow researchers to image every nerve cell and all chemical synapses in a given piece of brain tissue, opening up the possibility of mapping neuronal networks densely, not just sparsely. Although there have been substantial advances in imaging speed, the analysis of such 3D image data is still the limiting step. Therefore, dense reconstructions of cortical tissue have thus far been limited to femtoliter-scale volumes, keeping the systematic analysis of axons, neuronal cell bodies and their dendrites of different types, and the dense connectome between them out of reach.

### RATIONALE

Image analysis has made decisive progress using artificial intelligence–based methods, but the resulting reconstructions of dense nerve tissue are still too error-prone to be scientifically meaningful as is. To address this, human data analysis has been integrated into the generation of connectomes and it is the efficiency of this human–machine data analysis that now determines progress in connectomics. We therefore focused on efficiency gains by: (i) improving the automated segmentation quality, (ii) analyzing the automated segmentation for locations of likely errors and directing the human work to these locations only, and (iii) optimizing human data interaction by helping annotators to immediately understand the problem to be solved, allowing fast, in-browser parallel data flight, and by minimizing latency between annotator queries. With this, close to 100 student annotators solved hundreds of thousands of reconstruction problems within just 29 s each, including all preparation and transition time.

### RESULTS

We reconstructed 2.7 m of neuronal wires densely in layer 4 of mouse somatosensory cortex within only ~4000 invested human work hours, yielding a reconstruction ~300 times larger than previous dense cortical reconstructions at ~20-fold increased efficiency, a leap for the dense reconstruction of connectomes. The resulting connectome between 6979 presynaptic and 3719 postsynaptic neurites with at least 10 synapses each, comprising 153,171 synapses total, was then analyzed for the dense circuit structure in the cerebral cortex. We found that connectomic data alone allowed the definition of inhibitory axon types that showed established principles of synaptic specificity for subcellular postsynaptic compartments, but that at scales beyond ~5 μm, geometric predictability of the circuit structure was low and coarser models of random wiring needed to be rejected for dense cortical neuropil. A gradient of thalamocortical synapse density along the cortical axis yielded an enhanced variability of synaptic input composition at the level of single L4 cell dendrites. Finally, we quantified connectomic imprints consistent with Hebbian synaptic weight adaptation, obtaining upper bounds for the fraction of the circuit that could have undergone long-term potentiation.

### CONCLUSION

By leveraging human–machine interaction for connectomic analysis of neuronal tissue, we acquired the largest connectome from the cerebral cortex to date. Using these data for connectomic cell-type definition and the mapping of upper bounds for the learned circuit fraction, we establish an approach for connectomic phenotyping of local dense neuronal circuitry in the mammalian cortex, opening the possibility for the connectomic screening of nervous tissue from various cortices, layers, species, developmental stages, sensory experience, and disease conditions.

## Abstract

The dense circuit structure of mammalian cerebral cortex is still unknown. With developments in three-dimensional electron microscopy, the imaging of sizable volumes of neuropil has become possible, but dense reconstruction of connectomes is the limiting step. We reconstructed a volume of ~500,000 cubic micrometers from layer 4 of mouse barrel cortex, ~300 times larger than previous dense reconstructions from the mammalian cerebral cortex. The connectomic data allowed the extraction of inhibitory and excitatory neuron subtypes that were not predictable from geometric information. We quantified connectomic imprints consistent with Hebbian synaptic weight adaptation, which yielded upper bounds for the fraction of the circuit consistent with saturated long-term potentiation. These data establish an approach for the locally dense connectomic phenotyping of neuronal circuitry in the mammalian cortex.

The cerebral cortex of mammals houses an enormously complex intercellular interaction network implemented with neuronal processes that are long and thin, branching, and extremely densely packed. Early estimates indicated that 4 km of axons and 400 m of dendrites are compressed into a cubic millimeter of cortical tissue (*1*). This high packing density of cellular processes has made the locally dense mapping of neuronal networks in the cerebral cortex challenging.

So far, reconstructions of cortical tissue have been either sparse (*2*–*7*) or restricted to small volumes of up to 1500 μm^{3} (*8*–*10*). Consequently, the detailed network architecture of the cerebral cortex is unknown. Particular open questions are to what degree local neuronal circuits are explainable by geometric rules alone (*1*, *2*, *11*–*13*) and on which spatial scales cortical connectivity is only explainable by innervation preferences beyond such geometric models (*5*, *8*, *9*, *14*, *15*). Similarly, although numerous cortical neuronal cell types have been described based on protein expression, morphology, and electrophysiological characteristics (*16*), and these have been shown to have particular synaptic target patterns (*17*), the inverse question—whether, at the level of the dense cortical circuit, axons represent a continuum of synaptic preference or a set of distinct innervation paradigms that would allow for a purely connectomic cell type definition [as has been successful in the retina (*18*, *19*)]—is still open. Next, at the level of synaptic input to the primary dendrites of cortical excitatory cells, it is not known whether the typically three to 10 primary dendrites of a cortical neuron that leave the cell body homogeneously sample the available excitatory and inhibitory synaptic inputs or if there is an enhanced heterogeneity of synaptic input composition, making it possible to exploit the numerous mechanisms that have been discussed for the nonlinear integration of local synaptic inputs (*20*–*23*). Finally, whereas the change of synaptic weights in response to electrical and sensory stimulation has been widely studied (*24*–*28*) and connectomic data consistent with LTP have been described (*29*, *30*), the fraction of a given cortical circuit that is plausibly shaped by processes related to Hebbian learning under undisturbed conditions is still unknown.

We used dense connectomic reconstruction to quantitatively address these questions about the formational principles of a dense cortical circuit.

## Results

We acquired a three-dimensional (3D) EM dataset from upper layer 4 of primary somatosensory cortex of a 28-day-old mouse (Fig. 1, A to D, likely located within a barrel, see supplementary materials) using serial block-face electron microscopy [SBEM (*31*); dataset size: 61.8 × 94.8 × 92.6 μm^{3}; voxel size: 11.24 × 11.24 × 28 nm^{3}]. For dense reconstruction (Fig. 1, E to H), we 3D aligned the images and applied a sequence of automated analyses [SegEM (*32*), SynEM (*33*), ConnectEM, and TypeEM; Fig. 2, supplementary materials and methods, and table S2], followed by focused manual annotation (FocusEM). We reconstructed 89 neurons that had their cell body in the dataset (Fig. 1, E and F). These neurons constituted only 2.6% of the total path length (69 mm; Fig. 1G). To reconstruct axons, which constitute most of the wiring in the dense circuit (1.79 m, 66.6%, Fig. 1H), we applied a scalable distributed annotation strategy that identified locations of uncertainty in the automated reconstruction, which were then resolved by targeted manual annotation. To reduce the required manual annotation time, it was critical to obtain an automated reconstruction with low error rates, to use efficient algorithms for identifying locations for focused manual inspection (queries), and to minimize the time spent per user query. For this (Fig. 2A), we developed artificial intelligence–based algorithms that evaluated the EM image data and convolutional neural network (CNN)–filtered versions of the image data in the surrounding of interjunctions between segmented pieces of neurites (Fig. 2B). Together with classifiers that computed the probability of volume segments belonging to an axon, a dendrite, a spine head, or a glial process (using, among others, shape features; Fig. 2C), this allowed us to automatically connect parts of dendrites, attach spine heads to dendritic shafts (by a greedy stepwise agglomeration initiated at the spine head, Fig. 2D; 58.9% of spine heads unaffected by the dataset boundary were automatically attached), and reconstruct parts of axons. Similarly, synapses were automatically detected by evaluating pre- and postsynaptic volumes at neurite interfaces [Fig. 2E and figs. S2 to S5 (*33*); for shaft synapses, additional CNN-based classifiers for vesicle clouds and mitochondria were used]. To manually correct remaining errors in axons (Fig. 2, F to H), we detected ending locations of automatically reconstructed axon pieces (Fig. 2F) and directed user queries to these locations. For this, we used an egocentric 3D image display mode [“flight mode,” Fig. 2G (*34*)] and oriented the user annotation along the axis of the neurite for which a local annotation (“query”) was requested (movie S2). Together with data preloading, this yielded a low-latency, targeted neurite annotation in which individual user queries took 29.4 s to resolve (traveled path length per query: 5.49 μm). These queries could be easily distributed among 87 annotators. Similarly, we detected locations of likely mergers between axons (Fig. 2H, “chiasmata”) and directed user queries to reconnect the chiasma exits along actual axons. Using this scalable annotation architecture, we obtained a dense reconstruction of 2.69 m of neuronal processes (Fig. 1, G and H) with a total investment of 3981 human work hours, ~10 times faster than a recent dense reconstruction in the fly larval brain (*35*) (Fig. 2, I and J), ~20 times faster than the previous dense reconstruction in the mammalian retina (*18*), and ~25 times faster than the previous dense reconstruction in mammalian cortex (*9*) (Fig. 2, I and J). To quantify remaining reconstruction error rates in this dense neuropil reconstruction, we measured the remaining errors in a set of 10 randomly chosen axons and found 12.8 errors per millimeter of path length (of these, there were 8.7 continuity errors per millimeter; see materials and methods). This is indistinguishable from the error rates previously found in fast human annotations (*18*, *34*, *36*).

We obtained a connectome (Fig. 3) between 34,221 presynaptic axonal processes and 11,400 postsynaptic processes [6979 × 3719 connectivity matrix (Fig. 3E) when restricted to those pre- and postsynaptic neurites that established at least 10 synapses]. Among the postsynaptic processes, we classified *n* = 169 apical dendrites (ADs) that traversed the dataset along the cortical axis without connection to one of the neuronal cell bodies in the dataset (Fig. 2A), 246 smooth dendrites (SDs, Fig. 2B), 80 somata, 116 axon initial segments (AISs; Fig. 2C), and 89 proximal dendrite (PD) trees connected to a soma in the dataset (movie S1; note that some of these neurons also had ADs that were classified as PDs and not included in the AD definition above; see materials and methods and tables S1 and S2).

### Connectomic definition of axon types

We investigated whether, based solely on connectomic information (Fig. 3), we could extract the rules of subcellular innervation preference described for inhibitory axons in the mammalian cortex (*17*) and if such synaptic target preference could also be found for excitatory axons. We first measured the preference of each axon for innervating dendritic spine heads versus dendritic shafts and other targets (Fig. 4, A and B) because, in the mammalian cortex, most axons of inhibitory interneurons (INs) preferentially innervate the dendrites’ shafts or neuronal somata (*17*) and most excitatory glutamatergic axons preferentially innervate the spine heads of dendrites (*1*). The fraction of primary spine synapses per axon (out of all synapses of that axon) accordingly allowed the identification of spine-preferring, likely excitatory axons with at least 50% primary spine innervations (*n* = 5894 axons) and shaft-preferring, likely inhibitory axons with <20% primary spine innervations (*n* = 893 axons, or 13.2% of all axons; for exceptions to this rule and control measurements, see the supplementary materials and tables S1 and S2).

We then determined for each of the subcellular synaptic target classes defined above (Figs. 3 and 4C) the per-synapse innervation probability that would best explain whether an inhibitory axon establishes at least one synapse onto each of these targets. These inhibitory “single-hit” binomial innervation probabilities were 4.2% (somata), 17.8% (PD), 4.9% (SD), 3.3% (AD), and 0.5% (AIS) (Fig. 4D). We then computed the expected distribution of synapses per axon made onto each target class assuming the double-hit, triple-hit, etc., innervation probabilities are the same as the probability to establish at least one synapse onto that target. When comparing these target distributions with the measured distributions of synapses per axon onto each target class (Fig. 4E), we found that inhibitory axons established enhanced preference for somata (*p* = 2.4 × 10^{−34}, *n* = 893, one-sided Kolmogorov–Smirnov test), PDs (*p* = 6.0 × 10^{−77}), ADs (*p* = 2.5 × 10^{−4}), and to a lesser degree for SDs (*p* = 1.7 × 10^{−3}, table S1), but no enhanced preference for AISs in L4 (*p* = 0.648). AISs were synaptically innervated by 0.172 input synapses per micrometer of AIS length, but these innervations were not made by axons with an enhanced preference for AISs, unlike in supragranular and infragranular layers (*37*).

When performing the same analysis for excitatory axons (Fig. 4F), we found clear target preference for ADs (*p* = 2.5 × 10^{−34}, Fig. 4F), SDs (*p* = 7.6 × 10^{−25}), and PDs (*p* = 1.3 × 10^{−169}). By contrast, thalamocortical (TC) axons [detected using the criteria reported in (*38*); see fig. S6 and materials and methods] indicated a target preference for PDs (*p* = 2.5 × 10^{−31}), but not for ADs (*p* = 0.019) or SDs (*p* = 0.723). To determine the fraction of inhibitory and excitatory axons that had an unexpectedly high synaptic preference for one (or multiple) of the subcellular target classes, we applied the false detection rate criterion used for the determination of significantly expressed genes [*q* value (*39*); see materials and methods] and obtained lower bounds on the fractions of axons in the tissue that preferentially innervate the various subcellular target classes (Fig. 4G; at least 58.0% of inhibitory and 24.4% of excitatory axons). Inhibitory axons (Fig. 4H), but not excitatory axons (Fig. 4I), showed higher-order innervation preferences, indicating that at the level of the dense cortical circuit, synaptic target preferences established by axons were not a continuum but allowed cell-type classification without the need for measurements of neuronal morphology, electrical activity, protein expression, or transcription levels.

### Geometric sources of synaptic innervations

Could these local connectivity rules have been derived solely from the geometry of axons and dendrites? We first quantified the overall relation between the spatial distribution of axons and dendrites and the establishment of synapses between them (Fig. 5). One paradigm, originally proposed by Peters (*11*), states that TC axons entering a certain cortical tissue volume would sample the available cortical dendrites for synaptic innervation according to their relative prevalence in the tissue (*1*). This model (Fig. 5A) predicted the TC innervation of most cortical dendrites rather well, with the exception of smooth dendrites [an exception reported by White et al. (*14*)] and the enhanced TC innervation of PDs of layer 4 cells (*21*). When applied to corticortical excitatory and inhibitory axons (Fig. 5A), we found that this model predicted excitatory innervation of most spiny dendrites rather well, but again failed to predict innervation of SDs and the proximal bias of inhibitory synapses. Because this model [which has been most widely used for circuit inference (*2*, *12*)] implicitly accounts for the density of synapses along the presynaptic axons, it was capable of capturing the increased synapse density of TC axons (Fig. 5A). A simpler variant of the Peters model (*1*, *15*) (Fig. 5B), which uses the density of pre- and postsynaptic path length as basis for the synaptic innervation prediction, failed at predicting the TC innervation but captured the corticocortical innervation of spiny dendrites (Fig. 5B). We then analyzed whether a Peters model normalized for postsynaptic synapse density (Fig. 5C) would better capture synaptic innervation and found that, in fact, the dendritic model was a far better predictor of synaptic innervation (compare Fig. 5, C and B). This indicated that SDs and ADs sampled synaptic input according to the relative path length of the presynaptic axons (Fig. 5C). We then investigated whether a Peters model accounting for pre- *and* postsynaptic synapse densities would improve the innervation prediction (Fig. 5D). In this model, both the output and the input of cortical excitatory neurites were properly predicted, but the suppressed innervation of SDs and ADs by TC axons and the proximal bias of inhibitory axons was not. Notably, none of the Peters models could account for this proximal bias of inhibitory synapses [Fig. 5D; for other failures of Peters predictions, see, e.g., (*3*, *6*, *8*, *9*)].

More recently, the Peters model has been investigated for the close proximity between axons and dendrites on the scale of few micrometers (*8*, *9*) and concluded poor (*8*) or absent (*9*) geometric predictability of synaptic innervation. We used our larger dense reconstruction to investigate the geometric prediction over a substantially broader spatial scale from 1 to ~30 μm and accounted for inhibitory axons, excitatory axons, and postsynaptic target types (Fig. 5, E to H). We measured whether the postsynaptic membrane surface available within a certain radius *r*_{pred} around a given axon (Fig. 5E) would be a predictor of synaptic innervation for that given axon. We measured the available membrane surface belonging to the five subcellular target classes around all 6979 axons (Fig. 5F) and used a linear multinomial regression model to predict synaptic innervation from these data (Fig. 5G). Then, we computed the coefficient of determination (*R*^{2}) reporting the fraction of axonal synaptic innervation variance that could be explained purely based on the geometrical information (Fig. 5H; for details, see the materials and methods). In fact, for small spatial scales of 1 to 5 μm, the membrane surface available around an axon was a rather good predictor of synaptic innervation from excitatory axons (range, 16 to 90%, Fig. 5H; less so for inhibitory axons: range, 23 to 79%).

Would this imply that axonal and dendritic proximity at the single-axon level can be used to infer synaptic connectivity in the cortex (*13*)? We found that for the spatial alignment scales that can be achieved in light-microscopy–based neuron reconstructions from multiple animals (10 to 20 μm), predictability dropped substantially (Fig. 5H), making circuit inference by an emulation of growth processes based on light-microscopically aligned data (*13*, *40*) implausible.

### Subcellular synapse placement

We used our dense reconstruction to study the spatial distribution of synapses along somata and dendrites in the cortical neuropil. The density of TC synapses had a substantial dependence on cortex depth (Fig. 6, A to D): the absolute density of TC synapses in the volume increased by ~93% over 50 μm cortex depth (Fig. 6, A and B); the TC excitatory synapse fraction TC/(TC+CC) (where CC is cortico-cortical) increased by 82.6%, corresponding to an absolute increase in the TC synapse fraction of 5.8% per 50 μm cortex depth (Fig. 6D). This gradient was consistent with light-microscopic analyses of TC synapses indicating a decrease of TC synapse density from lower to upper L4 (*41*). Neither the inhibitory nor the corticocortical excitatory synapse densities showed a comparable spatial profile (Fig. 6C).

How is the synaptic TC gradient mapped onto the dendrites of L4 neurons along the cortex axis (Fig. 6, E to G)? One possibility is that the TC synapse gradient is used to enhance the variability of synaptic input composition between different primary dendrites of the L4 neurons such that a neuron’s dendrites pointing upward toward the pia would sample relatively less TC input than dendrites pointing toward the white matter. Alternatively, mechanisms to establish synaptic target preference (such as those reported in Fig. 4) could be used to counterbalance this synaptic gradient and equilibrate the synaptic input fractions on the differently oriented dendrites. Our analysis showed that, in fact, even at the level of single primary dendrites, TC input fractions were 1.28-fold higher for dendrites pointing upward toward the cortical surface versus downward toward the white matter (Fig. 6, F and G; TC input fractions of each dendrite were corrected for the entire neuron’s TC input fraction; for this analysis, see the materials and methods). We then investigated whether this differential composition of the excitatory inputs is accompanied by different compositions of the inhibitory input synapses (Fig. 6, H to L). We found that the fraction of TC input to a neuron’s dendrites was anticorrelated to the fraction of inhibitory synapses that originated from AD-preferring inhibitory axons (Figs. 6I and 4), both at the level of the input to L4 neurons and at the level of single primary dendrites of L4 neurons (Fig. 6, I and J). The effect was absent for all other synapse classes, most notably the soma-preferring inhibitory axons (Fig. 6K; see discussion).

### Connectomic mapping of the plasticity-consistent circuit fraction

The concept of Hebbian plasticity, thought to be at the core of experience-dependent changes of synaptic weights in the brain, makes predictions about the temporal evolution of synaptic weights in multiple synaptic contacts between the same pre- and postsynaptic neurons (AADD joint synapses; Fig. 7, A and B): Because Hebbian synaptic plasticity is dependent on the electrical activity of the pre- and postsynaptic neurons, which in a first approximation can be assumed to be similar at joint synapses, long-term potentiation (LTP) predicts joint synapses to become stronger and relatively more similar in weight (especially if synaptic weight saturates) and long-term depression (LTD) predicts joint synapses to become weaker and relatively more dissimilar in weight (but more similar if synaptic weights saturate; Fig. 7A). Models of LTP and LTD thus make particular predictions about the temporal evolution of joint synaptic weights, and the mapping of synaptic weights and synaptic weight similarity in the connectome allows the quantification of upper bounds on the fraction of the circuit that can have undergone such particular patterns of weight change before the connectomic experiment (we denote those synapse pairs for which such patterns of weight change occurred to a sufficient degree as “having undergone LTP/LTD”; see discussion).

We set out to leverage our large connectomic dataset (*n* = 5290 excitatory joint synaptic pairs onto spines; Fig. 7C) to map the relation between synaptic size and synaptic size similarity in joint synapse pairs [Fig. 7E; for visualization, the figure reports relative synaptic size *dis*similarity on the x-axis; for the utilization of the axon–spine interface area (Fig. 7D) as an indicator of synaptic weight (*42*, *43*)]. These data would allow us to determine upper bounds on the plasticity-consistent fraction of the circuit beyond the previous finding that in joint synapse pairs, synaptic size is more similar than for randomly shuffled synapse pairs (*9*, *29*, *30*, *44*).

Synaptic size similarity in joint synapse pairs showed a broad distribution (Fig. 7E). When comparing this distribution with the synaptic size and synaptic size similarity distribution obtained from a random assignment of the same synapses into “random pairs” (Fig. 7F and fig. S7, C and D), we observed that the population of oversimilar synapse pairs (Fig. 7F) was split into a region of oversimilar and large synapses (mean synaptic size 0.23 to 1.19 μm^{2}; 16 to 20% of all joint synapse pairs are found in this region; the above-random synapse pairs constitute 3.6 to 3.9% of all joint synapse pairs; see fig. S7, C and D, and the materials and methods for details of the region definition and statistics), and oversimilar and small synapses (mean synaptic size 0.06 to 0.2 μm^{2}; 15 to 19% of all joint synapse pairs were found in this region; 3.0 to 3.4% of all joint synapse pairs were above random in this region). With this information, we obtained upper bounds on the fraction of the circuit that can have undergone LTP and LTD with weight saturation (compare Fig. 7, F and A).

To what degree was the observed synaptic weight similarity a result of subtypes of neurons establishing differently sized synapses? Although the quantification of the upper bounds of the plasticity-consistent circuit fraction would remain unaffected, we could use this more detailed analysis to understand whether the plasticity-consistent circuit fractions were specific to types of neuronal connections.

First, we considered the possibility that certain presynaptic cell types made consistently larger or consistently smaller synapses (Fig. 7G). In this case, the distribution of synaptic weight similarity for same-axon different-dendrite (AADd) synapse pairs would also show a bias toward more similarly sized synapses. However, we found no such evidence (Fig. 7H), excluding cell-type-specific synapse size of either presynaptic (axonal) or postsynaptic (dendritic) origin as the cause of the observed oversimilar synapse pairs.

Next, we separated those connections established by TC axons from those made by the remaining excitatory (i.e., CC) axons (Fig. 7I and fig. S7, A and B). We found an excess of oversimilar synapse pairs in the TC connections as well, with 7 to 16% of pairs found in a region of overly similar and large synapses (i.e., an upper bound of 16% on LTP). The region of overly similar and small synapse pairs, however, only comprised 2 to 7% of joint synapse pairs. This remaining number of overly similar small synapse pairs could in fact be induced by the overly similar large synapse pairs (see the supplementary materials). At 28 days of age, ~3 weeks after the proposed critical period during which LTP can be induced in TC connections (*45*, *46*), a fraction of up to 16% of joint synapse pairs was still consistent with previous episodes of LTP that led to stabilized potentiated synapse pairs at dendritic spines (*47*, *48*), but 84% were not.

Repeating these analyses for other combinations of pre- and postsynaptic neurite types (Fig. 7J), we found upper bounds for LTP and LTD of ~10 to 20%. For each of these subtype-specific connections, we could then again analyze whether any purely presynaptic or purely postsynaptic subtype within the already type-selected connections (corresponding to squares in the table of Fig. 7J) could be the cause of the observed synapse similarity. For example, the connections from corticortical axons onto spiny L4 neurons (*49*) showed no evidence for presynaptic axonal subtypes yielding oversimilar synapses (Fig. 7K; for additional controls of these findings, see the supplementary materials; fig. S7, A and B; and table S1).

Together, these results provided a first quantitative upper bound on the fraction of the circuit consistent with previous episodes of saturated Hebbian synaptic plasticity leading to strengthening or weakening of synapses (a “connectomic fingerprint” of the maximum possible plasticity fraction of the circuit) and excluded obvious cell-type-based connection strength differences as the origin of these observations. Because these results were obtained from brains of untrained animals and were not the result of electrical or other stimulation (“plasticity induction”), these data may represent an unbiased screening of upper bounds of plasticity traces in local cortical circuits, for which the dense connectomic mapping was essential.

## Discussion

Using FocusEM, we obtained the first dense circuit reconstruction from the mammalian cerebral cortex at a scale that allowed the analysis of axonal patterns of subcellular innervation, ~300 times larger than previous dense reconstructions from cortex (*9*). Inhibitory axon types preferentially innervating certain postsynaptic subcellular compartments could be defined solely on the basis of connectomic information (Figs. 3 and 4). In addition to inhibitory axons, a fraction of excitatory axons also exhibited such subcellular innervation preferences (Fig. 4). The geometrical arrangement of axons and dendrites explained only a moderate fraction of synaptic innervation, revoking coarse random models of cortical wiring (Fig. 5). A substantial TC synapse gradient in L4 gave rise to an enhanced heterogeneity of synaptic input composition at the level of single cortical dendrites (Fig. 6), which was accompanied by a reduced innervation from AD-preferring inhibitory inputs. The consistency of synapse size between pairs of axons and dendrites signified fractions of the circuit consistent with saturated synaptic plasticity, placing an upper bound on the “learned” fraction of the circuit (Fig. 7). FocusEM allowed the dense mapping of circuits in the cerebral cortex at a throughput that enables connectomic screening.

### Synaptic input composition along L4 dendrites

Our finding of a covariation of enhanced TC inputs to L4 excitatory cells with reduced direct inhibitory input from AD-preferring INs (Fig. 6, H to K) could be interpreted in the context of a disinhibitory circuit described previously (*50*, *51*). Taking into account the preferential targeting of ADs and of soma-preferring parvalbumin (PV)–positive INs by somatostatin (SST)–positive INs, this could imply that SST-IN–based disinhibition can enhance TC input by silencing perisomatic PV inputs recruited by feedforward inhibition (*52*) and concomitantly reducing the direct inhibitory component from SST INs. In any case, this finding of per-dendrite input variation points to a circuit configuration in which TC input variability is enhanced between neurons of the same excitatory type in cortical layer 4, and furthermore provides evidence for a per-dendrite synaptic input composition of enhanced heterogeneity.

### Connectomic traces of plasticity

We interpreted the joint synapse data (Fig. 7) in terms of upper bounds of synapse pairs that could have undergone certain models of plasticity. Although this analysis detects those synapse pairs that were exposed to saturating plasticity (i.e., the possible plasticity event led to a final weight state of both synapses), an alternative interpretation is a dynamic circuit in which at any given point in time, only a fraction of synapses has expressed saturated plasticity, whereas other (or all) synapses are in the process of undergoing plastic changes. We expect that more elaborate plasticity models of entire circuits will also make testable predictions that are accessible by connectomic snapshot experiments as shown here.

### Outlook

The presented methods and results open the path to the connectomic screening of nervous tissue from various cortices, layers, species, developmental stages, sensory experiences, and disease conditions. The fact that even a small piece of mammalian cortical neuropil contains a high density of relevant information so rich as to allow the extraction of possible connectomic signatures of the “learnedness” of the circuit makes this approach a promising endeavor for the study of the structural setup of mammalian nervous systems.

## Materials and Methods

### Animal experiments

A wild-type (C57BL/6) male mouse was transcardially perfused at postnatal day 28 under isoflurane anesthesia using a solution of 2.5% paraformaldehyde and 1.25% glutaraldehyde (pH 7.4) following the protocol in (*53*). All procedures followed the animal experiment regulations of the Max Planck Society and were approved by the local animal welfare authorities (Regierungspräsidien Oberbayern and Darmstadt).

### Tissue sampling and staining

The brain was removed from the skull after 48 hours of fixation and sliced coronally using a vibratome. Two samples were extracted using a 1-mm biopsy punch (Integra Miltex, Plainsboro, NJ) from a 1-mm-thick slice 5 mm from the front of the brain targeted to layer 4 in the somatosensory cortex of the right hemisphere. The corresponding tissue from the left hemisphere was further sliced into 70-μm-thick slices followed by cytochrome oxidase staining, indicating the location of the coronal slice to be in barrel cortex.

Next, the extracted tissue was stained as described previously (*53*). Briefly, the tissue was immersed in a reduced osmium tetroxide solution (2% OsO_{4}, 0.15 M CB, 2.5 M KFeCN), followed by a 1% thiocarbohydrazide step and a 2% OsO_{4} step for amplification. After an overnight wash, the sample was further incubated with 1.5% uranyl acetate solution and a 0.02 M lead(II) nitrate solution. The sample was dehydrated with propylene oxide and EtOH, embedded in Epon Hard (Serva Electrophoresis GmbH, Germany), and hardened for 48 hours at 60°C.

### 3D electron microscopy experiment

The embedded sample was placed on an aluminum stub and trimmed such that the tissue was directly exposed on all four sides of the sample. The sides of the sample were covered with gold in a sputter coater (Leica Microsystems, Wetzlar, Germany). Then, the sample was placed into an SBEM setup [(*31*), Magellan scanning electron microscope, FEI Company, Hillsboro, OR, equipped with a custom-built microtome courtesy of W. Denk]. The sample was oriented so that the radial cortex axis was in the cutting plane. The transition between L4 and L5A was identified in overview electron microscopy (EM) images by the sudden drop in soma density between the two layers (Fig. 1C). A region of size 96 × 64 μm^{2} within L4 was selected for imaging using a 3 × 3 image mosaic, a pixel size of 11.24 × 11.24 nm^{2}, an image acquisition rate of 10 MHz, a nominal beam current of 3.2 nA (thus a nominal electron dose of 15.8 e^{–}/nm^{2}), an acceleration voltage of 2.5 kV, and a nominal cutting thickness of 28 nm. The effective data rate, including overhead time spent during motor movements for cutting and tiling, was 0.9 MB/s. A total of 3420 image planes were acquired, yielding 194 GB of data.

### Image alignment

After 3D EM dataset acquisition, all images were inspected manually and marked for imaging artifacts caused by debris present on the sample surface during imaging. Images with debris artifacts were replaced by the images at the same mosaic position from the previous or subsequent plane. First, rigid translation-only alignment was performed based on the procedures in (*53*). The following modifications were applied. When shift vectors were obtained that yielded offsets of >100 pixels, these errors were iteratively corrected by manually reducing the weight of the corresponding entry in the least-squares relaxation by a factor of 1000 until the highest remaining residual error was <10 pixels. Shift calculation of subsequent images in cutting direction was found to be the most reliable measurement and was therefore weighted 3-fold in the weighted least-squares relaxation. The resulting shift vectors were applied (shift by integer voxel numbers) and the 3D image data were written in KNOSSOS format (*34*, *36*). For further improvement, subimage alignment was applied (see the supplementary materials and methods).

### Methods description for software code

All routines described in the following are available as software at https://gitlab.mpcdf.mpg.de/connectomics/L4dense, which is the relevant reference for the exact sequence of processing steps applied. The following descriptions and the more detailed ones in the supplementary materials and methods are aimed at pointing to the key algorithmic steps rather than enumerating all detailed computations.

### Workflow for dense circuit reconstruction

The workflow for volume reconstruction of the acquired 3D EM volume (Fig. 2 and fig. S1) was as follows. We first detected blood vessels and cell bodies using automated heuristics, and then processed the remaining image volume using machine-learning-based image segmentation [CNN and watershed as described in SegEM (*32*)]. The result of this processing was 15 million volume segments corresponding to pieces of axons, dendrites, and somata (volume: 0.0295 ± 0.3846 μm^{3}, mean ± SD). We then constructed the neighborhood graph between all these volume segments and computed the properties of interfaces between directly adjacent volume segments. On the basis of these features, we trained a connectivity classifier (ConnectEM; Fig. 2, A and B) to determine whether two segments should be connected (along an axon or a dendrite or a glial cell) or if they should be disconnected. Using the SynEM classifier (*33*), we determined whether an interface between two disconnected processes corresponded to a chemical synapse and, if so, which was the presynaptic and which was the postsynaptic neurite segment (see below for more details). We furthermore trained a set of classifiers (TypeEM; Fig. 2C) to compute for each volume segment the probability of being part of an axon, a dendrite, a spine head, or a glia cell (precision and recall were 91.8 and 92.9% for axons, 95.3 and 90.7% for dendrites, 97.2 and 85.9% for astrocytes, and 92.6 and 94.4% for spine heads, respectively; see table S2).

### Cell body–based neuron reconstruction

We next reconstructed those neurons that had their cell bodies in the tissue volume (Fig. 1, E and F, cell gallery in movie S1; *n* = 125 cell bodies; of these, 97 were neuronal and of these 97, 89 were reconstructed with dendrites in the dataset). For this, we used a set of simple growth rules for automatically connecting neurite pieces on the basis of the segment-to-segment neighborhood graph and the connectivity and neurite type classifiers (fig. S1, “automated agglomeration”; see the supplementary materials and methods). As a result, we obtained fully automated reconstructions of the neuron’s soma and dendritic processes. With a minimal additional manual correction investment of 9.7 hours for 89 cells (54.5 mm dendritic and 2.1 mm axonal path length), the dendritic shafts of these neurons could be reconstructed without merge errors, but there were 37 remaining split errors, at 87.3% dendritic length recall (table S2). This reconstruction efficiency compares favorably to recent reports of automated segmentation of neurons in 3D EM data from the bird brain obtained at ~2-fold higher imaging resolution (*54*), which reports soma-based neuron reconstruction at an error rate of beyond 100 errors per 66 mm dendritic shafts at lower (68%) dendritic length recall with a similar resource investment (see the supplementary materials and methods).

In addition to the dendritic shafts, the dendritic spines constitute a major fraction of the dendritic path length in cortical neuropil (Fig. 1G). Using our spine head classifier (part of the TypeEM classifiers; Fig. 2C), we found 415,797 spine heads in the tissue volume, which is a density of 0.784 per μm^{3} (0.98 per μm^{3} of neuropil, when excluding somata and blood vessels). To connect these to the corresponding dendritic shafts, we trained a spine neck continuity algorithm that was able to automatically attach 58.9% of these spines (evaluated in the center of the dataset at least 10 μm from the dataset border), yielding a dendritic spine density of 0.672 per μm dendritic shaft length [comparable to spine densities in the bird brain (*55*)]. However, in mammalian cerebral cortex, the density of spines along dendrites is even higher (at least 1 per μm dendritic shaft length). The remaining spine heads were then attached to their dendritic shafts by seeding manual reconstructions at the spine heads and asking annotators to continue along the spine necks to the dendritic shafts. This annotation was performed in the “orthogonal mode” configuration of webKnossos (*34*), in which the annotator viewed three orthogonal image planes to decide where to continue the respective spine neck [as in KNOSSOS (*36*)]. The annotation of all remaining spine necks consumed an additional 900 hours of human work for the attachment of 98,221 spines, resulting in a final overall spine density of 0.959 per μm dendritic shaft length.

### Dense tissue reconstruction

The reconstruction of neurons starting from their cell bodies, however, was not the main challenge. Rather, the remaining processes, axons and dendrites not connected to a cell body within the dataset and densely packed in the tissue constitute ~97% of the total neuronal path length in this volume of cortex (Fig. 1G). To reconstruct this vast majority of neurites (Fig. 1H), we first used our connectivity and neurite type classifiers (ConnectEM and TypeEM, respectively; Fig. 2) to combine neurite pieces into larger dendritic and axonal agglomerates (“automated agglomeration,” fig. S1 and supplementary materials and methods). Then, we took those agglomerates that had a length of at least 5 μm (*n* = 74,074 axon agglomerates), detected their endings that were not at the dataset border, and directed focused human annotation to these endings (“queries,” Fig. 2, F and G).

For human annotation, we used an egocentric directed 3D image data view (“flight mode” in webKnossos), which we had previously found to provide maximized human reconstruction speed along axons and dendrites in cortex (*34*). Here, however, instead of asking human annotators to reconstruct entire dendrites or axons, we only queried their judgment at the endings of automatically reconstructed neurite parts. To make these queries efficient, we made three additions to webKnossos: (i) we oriented the user along the estimated direction of the neurite at its ending, reducing the time the user needs to orient within the 3D brain tissue; (ii) we dynamically stopped the user’s flight along the axon or dendrite whenever another of the already reconstructed neurite agglomerates had been reached; and (iii) we preloaded the next query while the user was annotating (Fig. 2, F and G). Movie S2 illustrates this annotation process for cases of splits and mergers, respectively. Note that the user was able to switch quickly to the next query and, based on its 3D orientation, spent little time orienting in the tissue at the new location. With this, the average user interaction time was 21.3 ± 36.1 s per query, corresponding to an average of 5.5 ± 8.8 μm traveled per query. In total, 242,271 axon-ending queries consumed 1978 paid-out work hours (i.e., including all overheads, 29.4 s per query).

However, we had to account for a second kind of reconstruction error, so-called mergers, which can originate from the original segmentation, the agglomeration procedure, or erroneous flight paths from human queries (Fig. 2H). To detect such mergers, we started with the notion that most of these merger locations will yield a peculiar geometrical arrangement of a 4-fold neurite intersection once all neurite breaks have been corrected (Fig. 2H, “chiasma”). Because such chiasmatic configurations occur rarely in branching neurites, we directed human focused annotation to these locations. First, we automatically detected these chiasmatic locations using a simple heuristic to detect locations at which axon-centered spheres intersected more than three times with the axon [Fig. 2H, *n* = 55,161 chiasmata; for approaches to detect such locations by machine learning, see (*56*, *57*)]. Then, we positioned the user queries at a certain distance from the chiasma location pointing inward (Fig. 2H) and used a set of case distinctions to query a given chiasma until its configuration had been resolved (see the supplementary materials and methods for details). Chiasma annotation consumed an additional 1132 work hours [note that the detection of endings and chiasmata was iterated eight times for axons (see the supplementary materials and methods) and that, in a final step, we also detected and queried 3-fold neurite configurations to remove remaining mergers].

### Synapse detection, types of postsynaptic targets, and connectome reconstruction

Given the reconstructed pre- and postsynaptic neurites in the tissue volume, we then went on to extract their connectome. For this, we used SynEM (*33*) to detect synapses between the axonal presynaptic processes and the postsynaptic neurites (Fig. 2E).

We trained a dedicated interface classifier for nonspine synapses using training data containing only shaft and soma synapses (Figs. S2 to S5; see the supplementary material and methods). This classifier also used four additional texture filters compared with SynEM in (*33*), which originated from the voxelwise predictions of a multiclass CNN trained on synaptic junctions, vesicle clouds, mitochondria, and a background class (Fig. 2E).

Because we were interested in analyzing the subcellular specificity of neuronal innervation, we had to also classify which of the postsynaptic membranes belonged to cell bodies; to classify spiny dendrites as belonging to excitatory cells and smooth dendrites as belonging to INs; and to detect AISs and those dendrites that were likely ADs of neurons located in deeper cortical layers. We developed semiautomated heuristics to detect these subcellular compartments (Fig. 3, A to D; see the supplementary materials and methods for details).

### Definition of excitatory and inhibitory axons

We used the fraction of primary spine synapses per axon (out of all synapses of that axon, only axons with at least 5 μm path length and at least 10 synapses were analyzed), which had a peak at ~80% (Fig. 4, A and B), to identify spine-preferring, likely excitatory axons with at least 50% primary spine innervations. Similarly, we identified shaft-preferring, likely inhibitory axons with <20% primary spine innervations. Together, this yielded 6449 axons with clear shaft or spine preferences. For the remaining *n* = 528 axons with primary spine innervations >20% and <50%, we first wanted to exclude remaining mergers between excitatory and inhibitory axons (that would yield intermediate spine innervation rates) and split these axons at possible merger locations (at least 3-fold intersections). Of these, 338 now had at least 10 synapses and spine innervation rates <20% or >50%. The remaining *n* = 192 axons (2.75% of all axons with at least 10 synapses) were not included in the following analyses. This together yielded *n* = 5894 excitatory and *n* = 893 inhibitory axons in our data. For additional controls, see the supplementary materials.

TC axons were defined following parameters described previously (*38*) (see the supplementary materials and methods).

### Analysis of subcellular synaptic target preference

First, we assumed that all synapses of a given axon class have the same probability to innervate a particular postsynaptic target class (as above). We then inferred this single-hit innervation rate for all combinations of presynaptic axon classes and postsynaptic target classes by determining the probability that best explains whether an axon innervated the target class under a binomial model. The optimized binomial model was then used together with the measured number of synapses of each axon to calculate the expected distribution of target innervation rates. A one-sided Kolmogorov–Smirnov test was used to search for the existence of a subpopulation with an increased target innervation rate. To identify those axons that innervated a given target class beyond chance (Fig. 4G), we computed the probability *p*^{(}^{t}^{)}_{meas,}* _{i,k}* of finding at least the measured fraction of synapses onto target

*t*for each axon

*i*from axon class

*k*. The

*p*values were also calculated for the expected distribution of target innervation rates and combined with

*p*

^{(}

^{t}^{)}

_{meas}

*to estimate the*

_{,i,k}*p*-value threshold

*p̂*

^{(}

^{t}^{)}

*at which the false discovery rate*

_{k}*q*(

*39*) crosses 20%. Eighty percent of the axons with

*p*

^{(}

^{t}^{)}

_{meas,}

_{i,k}< p̂^{(}

^{t}^{)}

*innervate target*

_{k}*t*with a rate above the single-hit innervation probability and are thus considered to be

*t*preferring.

For the analysis of second-order innervation preference (Fig. 4, H and I), we reported the fraction of synapses onto target τ by *t*-preferring axons of class *k* after removal of synapses onto *t*. This innervation rate was compared against the fraction of synapses onto target τ by all axons of class *k*.

### Geometrical predictability analysis

Peters’ rule (*1*) stipulates that synapses between classes of axons and dendrites are established in proportion to the prevalence of these classes. One variant of Peters’ rule considered (Fig. 5B) makes the prediction that the fraction of synapses from axon class A onto target class T is the product of *p*_{A} and *q*_{T}, where *p*_{A} is the proportion of axonal path length made up by class A, and *q*_{T} is the proportion of dendritic path length (excluding spines) made up by class T. The measured synapse fractions were compared against the predictions by calculating the ratio of observed to predicted synapse fractions.

Other formulations evaluate these predictions independently for each axon class (Fig. 5A) or each dendrite class (Fig. 5C).

Finally, to assess the effect of incorporating explicit knowledge about the synapse densities of different axon and dendrite classes, a fourth variant of Peters’ rule (Fig. 5D) was considered in which the predicted synapse fraction from axon class A onto target class T is the product of *p*′_{A} and *q*′_{T}, the overall fractions of synapses originating from A and innervating T, respectively.

How much additional information about the neuropil composition around an axon helps to predict its postsynaptic targets was assessed as follows. For each axon, we determined the total surface area of the target classes that were contained within the cylinder of radius *r*_{pred} around the axon (Fig. 5E) and compared it with the actually innervated target fraction of each axon (Fig. 5, E and F). We then analyzed the correlation between the availability of the target surfaces and the established synapses on these target classes (Fig. 5G).

We then computed *R*^{2} using the following model. For all axons of a given type, we used the fraction of target innervations and fractional surface availabilities in a surround of radius *r*_{pred} to find the optimal multivariate linear regression parameters. To estimate best-case geometric predictability, we then calculated *R*^{2} as 1 minus the ratio between the squared residuals of the regression and the synaptic variance on the same axons used for parameter optimization. Here, we corrected for the variance introduced by the finite number of synapses per axon: we used the axons’ fractional surface availabilities within *r*_{pred} and absolute synapse numbers to calculate the expected binomial variance and subtracted it from the squared residuals.

This analysis made several assumptions that were in favor of a geometrical explanation of synaptic innervation [therefore, the conclusions about a minimal predictability (Fig. 5H) are still upper-bound estimates]. It was assumed that the number of synapses for a given axon was already known; in most settings, only average synapse rates are known for a given circuit. It also assumed that a precise knowledge of the axonal trajectory and the surrounding target surface fractions were available; again, this is usually only available as an average on the scale of *r*_{pred} of several tens of micrometers.

To relax the assumption of complete knowledge about target availabilities, we repeated the above *R*^{2} analysis for a model in which the predicted fractional innervation of a target is the fractional surface availability of that target.

The computational routine used can be found at https://gitlab.mpcdf.mpg.de/connectomics/L4dense in +connectEM/+Connectome/plotGeometricPredictability.m.

### Synapse-size consistency analysis

To determine the consistency of primary spine synapses between a given axon–dendrite pair, we calculated the axon–spine interface area (ASI) (*42*, *43*) of a synapse as the total contact area between the corresponding axon and spine head agglomerates. For axon–dendrite pairs connected by exactly two primary spine synapses, we then calculated the coefficient of variation (CV) of the ASI areas by CV = 2^{1/2}(ASI_{1} – ASI_{2})/(ASI_{1} + ASI_{2}), where ASI_{1} and ASI_{2} are the larger and smaller of the two ASI areas, respectively. To avoid false same-axon, same-dendrite (AADD) pairs caused by remaining merge errors in the axon reconstruction, this analysis was performed only after splitting axons at their branch points. The measured distribution of CV values was compared against the CV values obtained by randomly drawing pairs from all AADD synapses and against the CV values of observed AADd synapse pairs and pairs from different axons onto the same dendrite (AaDD) and from different axons onto different dendrites (AaDd; Fig. 7H). To test whether AADD primary spine synapse pairs are more similar in size than pairs in the control conditions, a one-sided Kolmogorov–Smirnov test was used. We calculated the decimal logarithm of the average ASI area (in square micrometers) and the CV of the ASI areas of each synapse pair to map the size-similarity plane (Fig. 7, F and I). The kernel density estimate of the observed distribution was compared against the distribution expected from random pairs (5000 Monte Carlo samples; fig. S7C) to identify statistically significantly overrepresented regions. Contour lines show the intersection of the significance regions for *p*-value thresholds of 0.5% and 5% (Fig. 7, E, F, I, and fig. S7, C and D), with the convex hull around the set of all data points. The fraction of data points contained within a contour was used as the upper bound on the fraction of connections consistent with saturated Hebbian plasticity (Fig. 7J).

### Statistical methods

The following statistical tests were performed (in order of presentation in the figures):

The existence of axon subpopulation with unexpectedly high synapse rate onto a given target class was tested using the one-sided Kolmogorov–Smirnov test (Fig. 4, E and F). Axons belonging to a given target-preference class were identified on the basis of the false detection rate criterion [*q* = 20% (*39*)] (Fig. 4G).

The degree to which synaptic variance is explainable by geometry-based models was evaluated using *R*^{2} (Fig. 5H). Binomial variance was corrected for by subtracting the surface fraction-based expected binomial variance from the squared residuals.

*F* tests were used to evaluate synaptic gradients as function of cortical depth (Fig. 6, B and D) or dendritic orientation (Fig. 6, F and G). For correlation of the TC input fraction with other synaptic input fractions along dendrites, the inhibitory input fraction and seven target-preferential inhibitory and excitatory synapse types were tested. AD-preferring inhibitory synapses were the only ones with significant and substantial correlation (Pearson’s correlation after Bonferroni’s correction for *n* = 8 multiple tests). The correlation was also significant at the soma level (Pearson’s correlation). Both correlations were also significant using Spearman’s rank correlation.

The four variants of Peters’ rule (Fig. 5, A to D) were compared using a likelihood-ratio test based on the following multinomial model. It was assumed that the pre- and postsynaptic classes of each synapse in the connectome were sampled either after the path-length fractions of these classes (*p*_{A} and *q*_{T}) or after the product of the path length and a class-specific likelihood-maximizing relative synapse density. Wilks’ theorem was used to compute the corresponding *p* values.

To test whether the axon–spine interface areas of a given spine–synapse pair configuration were more similar than randomly sampled pairs, a one-sided Kolmogorov–Smirnov test was used (Fig. 7, H and K).

## Supplementary Materials

science.sciencemag.org/content/366/6469/eaay3134/suppl/DC1

Materials and Methods

Supplementary Text

Figs. S1 to S7

Tables S1 and S2

Captions for Movies S1 and S2

This is an article distributed under the terms of the Science Journals Default License.

## References and Notes

**Acknowledgments:**We thank Y. Hua and K. M. Harris for discussions; A. Gour, Y. Hua, P. Laserstein, and H. Schmidt for discussing unpublished observations; B. Cowgill for experimental contributions; E. Eulig, R. Hesse, C. Schramm, and M. Zecevic for code contributions; C. Guggenberger for compute infrastructure; Scalable Minds for implementation of the fast-swapping user queries in webKnossos; and annotators (see the supplementary materials) for reconstruction work.

**Funding:**Funding was provided by the Max Planck Society.

**Author contributions:**M.H. conceived, initiated, and supervised the project; K.M.B. performed experiments; P.H. provided analysis methods; A.M., M.B., K.M.B., and B.S. developed and performed analyses with contributions by all authors; and A.M. and M.H. wrote the manuscript with contributions by all authors.

**Competing interests:**The authors declare no competing financial interests.

**Data and materials availability:**The raw EM data are available for browsing at https://wklink.org/9276 (see data links in the figure legends). All data and code are publicly available at https://L4dense2019.brain.mpg.de and in updated form at https://gitlab.mpcdf.mpg.de/connectomics/L4dense. Less than 2% of the 3D EM image dataset using the initial image alignment was previously published (

*32*–

*34*,

*58*).