Neocortical layer 4 in adult mouse differs in major cell types and circuit organization between primary sensory areas

Layer 4 (L4) of mammalian neocortex plays a crucial role in cortical information processing, yet a complete census of its cell types and connectivity remains elusive. Using whole-cell recordings with morphological recovery, we identified one major excitatory and seven inhibitory types of neurons in L4 of adult mouse visual cortex (V1). Nearly all excitatory neurons were pyramidal and almost all Somatostatin-positive (SOM+) neurons were Martinotti cells. In contrast, in somatosensory cortex (S1), excitatory cells were mostly stellate and SOM+ cells were non-Martinotti. These morphologically distinct SOM+ interneurons correspond to different transcriptomic cell types and are differentially integrated into the local circuit with only S1 cells receiving local excitatory input. Our results challenge the classical view of a canonical microcircuit repeated through the neocortex. Instead we propose that cell-type specific circuit motifs, such as the Martinotti/pyramidal pair, are optionally used across the cortex as building blocks to assemble cortical circuits.

To support our expert classification, we fully reconstructed a subset of neurons from each 142 inhibitory type (n=92 in total) and trained a regularized logistic regression classifier to 143 discriminate between each pair of inhibitory cell types (see Methods). We used 2D density maps 144 and a set of morphometric statistics (Fig. S3) as predictors 43 . Across all 21 pairs, the average 145 cross-validated classification accuracy was 0.92, with most pairs discriminated almost perfectly 146 ( Fig. 2a, left). However, SC/HEC and SC/NGC pairs showed only ~0.65 classification accuracy. 147 Visualisation of this dataset with t-SNE (Fig. 2a, right) indicated that SC/HEC and SC/NGC 148 types were partially overlapping, as well as BC/DBC. Overall, this analysis suggests that while 149 most morphological classes can be very well discriminated, some may be partially overlapping. To further explore variability in electrophysiological properties between cell types, we 166 characterized the firing pattern of a subset of neurons (n=235) using 13 automatically extracted 167 electrophysiological features (Fig. S4). Most features exhibited strong differences between cell 168 types (Fig. S5). Two-dimensional visualisation of this dataset using t-SNE (Fig. 2b) showed that 169 all four PV + cell types overlapped in one group of electrophysiologically similar FS neurons, 170 while the other four types (PYR, NGC, BPC, and MC) each had distinct firing patterns. We 171 confirmed this using pairwise classification with regularized logistic regression (Fig. 2b): the 172 average cross-validated classification accuracy between the FS types was only 0.67, while the 173 average accuracy across all other pairs was 0.98. 174 175 V1 differs from S1 in major L4 cell types 176 In contrast to V1 L4, stellate cells are known to be abundant in S1 L4 of rats and mice [1][2][3][4][5][6][7] . To 177 confirm this, we recovered L4 excitatory cells in S1 (n=24 mice, including n=5 Scnna1-Cre/Ai9) 178 with the same method as in V1. We found that indeed 76.6% (85/111) of the recovered spiny 179 neurons did not have a clear apical dendrite and were thus classified as stellate cells (Fig. 3B) The NMCs also differed in their firing pattern from MCs recorded in V1: they had a higher 195 maximal firing rate, a lower AP width, and a lower membrane time constant (Fig. 4a and Fig.  196 S5). This resembles the FS firing pattern, and one previous study called NMCs "quasi-FS" 13 . 197 Comparison of electrophysiological features between MCs, NMCs, and FS cells revealed that 198 NMCs were "in between" MCs and FS cells in terms of their firing patterns and intrinsic 199 membrane properties (Fig. S5, S6). 200

201
To further investigate the differences between MCs in V1 and NMCs in S1, we used the Patch-202 seq 29-31 technique which combines patch-clamp recordings with single cell transcriptomics. 203 Using n=6 SOM-Cre/Ai9 mice, we sequenced RNA of SOM-Cre + neurons in L4 of V1 and S1 204 (n=42 in V1 and n=35 in S1 after quality control), and also in L5 of each area (n=17 and n=16 205 respectively). We obtained on average 1.1 million reads per cell (median; mean±SD on a log10 206 scale: 6.0±0.3) and detected 6.4±1.6 thousand (mean±SD) genes per cell (Fig. S7). We 207 mapped these cells to a large transcriptomic cell type dataset 24 that contained 21 somatostatin 208 types with 2880 neurons from V1 and ALM. The quality of the mapping was equally good for V1 209 and S1 cells (Fig. S7), suggesting that the V1+ALM dataset is a reasonable reference for S1 210 interneurons. This is in agreement with the idea that inhibitory transcriptomic cell types are 211 shared across cortical regions 24,25 . Three cells (excluded from the counts given above and from 212 further analysis) had fast-spiking firing pattern, did not express SOM, and mapped to Pvalb Reln 213 Itm2a transcriptomic type, likely corresponding to the basket cells that we found labeled in the 214 SOM-Cre line (Fig. 3). All other cells mapped to Sst transcriptomic types. To answer this question, we looked at electrophysiological features that were most different 234 between SOM + interneurons in V1 and S1 (Cohen's d>1: input resistance, membrane time 235 constant, maximum firing rate, and AP width) and found that for two of them (input resistance 236 and membrane time constant) V1 cells belonging to the Hpse type had values more similar to 237 the S1 cells than to the V1 cells from the Calb2 type (Fig. 4d). This suggests that 238 electrophysiologically, V1 Hpse MC cells are in between V1 Calb2 MC cells and S1 NMC cells. 239

240
The relationship between gene expression and electrophysiological features can be visualized 241 using the sparse reduced-rank regression analysis that we have recently introduced 46 . This 242 technique aims to reconstruct all the electrophysiological features using a two-dimensional 243 13 projection of the expression levels of a small set of genes (Fig. 4d). The optimal number of 244 genes was found using cross-validation (see Methods). This analysis supports our conclusion 245 that V1 Hpse MCs are "in between" Calb2 MCs and NMCs in terms of electrophysiology. 246 Interestingly, this analysis also showed that some of the cells assigned to the Tac1 and Mme 247 types had a distinct fast-spiking-like firing pattern which was different from firing patterns of MCs 248 and NMCs (but was not as sustained as the proper FS pattern). These three SOM + 249 transcriptomic cell types have recently been identified in Tasic  The L5 SOM + cells that we sequenced in both areas mostly mapped to a different set of 279 transcriptomic types than the L4 SOM + cells, but there were no apparent differences between 280 S1 and V1 in terms of transcriptomic cell types (Fig. 4b). 281 282 Connectivity among excitatory and SOM + neurons in L4 of V1 vs. S1 283 So far, we have described major differences in the morphology, electrophysiology, and 284 transcriptomic signatures of excitatory neurons and SOM + interneurons in L4 between V1 and 285 S1. We next investigated whether there are differences in their connectivity profiles as well, 286 using simultaneous multi-cell patch-clamp recordings. We found that certain connectivity 287 patterns between them are very similar across both areas (Fig. 5). First, the connection 288 probabilities among excitatory cells were low in both areas (1.0%, 7/701 in V1; 2.5%, 3/122 in 289 S1). Second, the connection probabilities between SOM + cells were also low in both areas (0%, 290 0/68 in V1; 3.8%, 2/52 in S1). Third, the connection probabilities from SOM + cells to excitatory 291 cells were high in both areas (21.1%, 30/142 in V1, 26.6%, 17/64 in S1). In addition, despite 292 their low connectivity via chemical synapses, both MCs in V1 and NMCs in S1 were similarly On the other hand, we found a striking area-specific difference in connection probabilities from 296 excitatory to SOM + neurons. In S1, NMCs received facilitating synaptic connections from local 297 excitatory neurons (12.5%, 8/64), in line with previous studies in young rodents 14,16 . In contrast, 298 we did not find any connections (0%, 0/142) from local excitatory neurons to MCs in V1 299 (p=0.0002, Fisher's exact text). This was also in stark contrast to MCs in L2/3 and L5 of adult 300 mouse V1, which receive strong facilitating synaptic inputs from local PYRs in the same layers 301 28 (see Discussion for further considerations). 302 303 In addition, we tested the connectivity in V1 L4 including BCs (Fig. S10). We found that BCs 304 followed the same connectivity rules as previously found in other layers 28,48 . PYRs connected to 305 BCs with probability 12.5% (38/303), MCs inhibited BCs with probability 32.6% (15/46), and BCs 306 inhibited each other (36.7%, 75/204), MCs (13.0%, 6/46) and PYRs (25.7%, 78/303). All of 307 these connection patterns have also been reported in S1 L4 of young mice 23 . We also found 308 that BCs were electrically coupled to each other with probability 27.5% (28/102) but were never 309 electrically coupled to MCs (0/23), in agreement with previous findings that gap junctions exist 310 between inhibitory cells of the same type 49 . 311 312 Notably, the connection probability between PYRs in V1 L4 was very low, consistent with our 313 previous work in other layers in adult animals 28 , but in contrast to the findings in young and 314 juvenile rodents 50,51 . To confirm that this low connectivity reflects an age effect, we measured 315 the connectivity between PYRs in V1 L4 at different ages (P15-20 and P30-40, n=5 each) using 316 Scnn1a-Cre/Ai9 mice. We found that the connection probability monotonically decreased with 317 age (Fig. S11): from 13.2% in P15-20 (15/114) to 5.1% in P30-40 (8/156) to 1.0% (7/701) 318 reported above for the P55+ mice with median age P71. This is in agreement with the recent 319 study that found 6.3% (20/315) connection probability in V1 L4 of P47±6 mice 52 . 320 321 When measuring connectivity in S1 L4, no special care was taken to ensure that the tested cells 322 were within the same barrel. At the same time, it is known that cells in S1 L4 preferentially make 323 intra-barrel connections 2,3 . To address this concern, we performed a separate series of 324 experiments using n=8 Scnn1a-Cre/Ai9 mice to test intra-barrel connectivity of excitatory 325 neurons. We used the tdTomato fluorescence signal to detect the barrels during patch clamp 326 recordings 53 and performed cytochrome oxidase staining in a subset of slices to confirm that 327 the fluorescence signal reliably corresponded to barrel boundaries 3 (see Methods and Fig.  328 S12). The measured connection probability was 5.2% (5/104) which was larger than the value 329 reported above (2.5%, 3/122) but not significantly different from it (p=0.48, Fisher's exact test). 330 Both estimates are substantially lower than the existing estimates of intra-barrel connectivity 331 obtained in young animals (30--35%) 2,3,54 which is in line with the decrease in local excitatory 332 connectivity with age discussed above for V1 (see also question, why pyramidal cells in rodent V1 L4 remain pyramidal, whereas L4 excitatory cells in 361 rodent S1 and in V1 of other non-rodent species develop into stellate cells. We suggest one 362 hypothesis below when discussing their circuit organization. 363

364
We found that all non-fast-spiking SOM + neurons in V1 L4 are Martinotti cells (MCs), which is 365 also in contrast to S1 L4 where almost all SOM + neurons are non-Martinotti 13,14 . Using Patch-366 seq, we showed that SOM + MCs in V1 L4 and SOM + NMCs in S1 L4 correspond to two different There is a broad agreement between their types from L4 and our types. There are also some 376 differences: they split abundant types (e.g. PYRs and BPCs) into multiple narrow sub-clusters, 377 while at the same time missing some rare types such as HECs. Cbln4 respectively), although distinct, were close and partially overlapping in the t-SNE 384 visualisation (Fig. 4b). Consistent with this, Tasic et al. 24

also found intermediate cells between 385
the "core" members of these two clusters. We showed that, electrophysiologically, MCs and 386 NMCs also form a continuum (Fig. 4d, Fig. S6) with all electrophysiological features having 387 unimodal distributions (Fig. 4c). This is in agreement with the findings of Naka et al. 45 who 388 demonstrated an electrophysiological continuum between NMCs in S1 L4 and MCs in S1 L5. correspond to MCs based on our data. We hypothesize that these cells may be "latent NMCs", 397 present in V1, but failing to develop a NMC morphology due to the nearly complete absence of 398 stellate cells in V1. Tasic et al. 24 showed that the majority of transcriptomic inhibitory types are 399 20 shared between two very different cortical areas (V1 and ALM). Our findings demonstrate that 400 this does not necessarily imply that morphological types are also all shared. 401 402 Using Patch-seq, we also performed single-cell RNA-sequencing of a small number of L5 SOM + 403 cells in both S1 and V1. Morphologically, almost all SOM + cells in V1 L5 (except some fast-404 spiking cells) 28,45 and the majority of SOM + cells in S1 L5 21,28,45,61 are known to be MCs. We 405 found that L5 SOM + cells had electrophysiological features similar to L4 MCs (Fig. 4d), but 406 mostly mapped to a different set of transcriptomic clusters than the L4 SOM + cells (Fig. 4b). types (m-types) as electrophysiological ones (e-types). We only obtained the transcriptomic 419 information for SOM + neurons, but found out that MCs in V1 L4 could belong to two different 420 transcriptomics types (t-types), one of which corresponded to NMCs in S1; inside V1, the cells 421 from these two t-types had slightly different electrophysiology (Fig. 4). An integrative definition of 422 cell type in V1 should take this into account. 423 424 21 Importantly, the SOM-Cre line does not label neurons in exact correspondence with 425 transcriptomic classes. We found that ~8% of V1 L4 neurons labeled in the SOM-Cre line were 426 fast-spiking cells with the morphology of basket cells (Fig. 3), in agreement with previous reports 427 28,41,42 . In our patch-seq experiments we found three sequenced SOM-Cre + neurons that were 428 fast-spiking and mapped to Pvalb transcriptomic types. We did not detect SOM (zero read 429 count) in either of these three cells, suggesting that they likely had transiently expressed it 430 during development, as hypothesized by Hu et al. 41 Interestingly, all three cells mapped to the 431 same Pvalb type: Pvalb Reln Itm2a. 432 433 Circuit organization in L4: V1 vs. S1 434 In terms of connectivity, both MCs in V1 and NMCs in S1 avoid connecting to each other (apart 435 from forming gap junctions; Fig. S8), and project to excitatory population in L4. Moreover, the 436 axonal morphologies of these two cell types seemed to match the respective dendritic 437 morphologies of their excitatory neuronal targets. In V1, axons of L4 MCs primarily projected to 438 L1 where they are potentially able to synapse onto the tuft of L4 PYRs, similar to the pattern 439 described in other cortical layers 63,64 . In S1, by contrast, axons of L4 NMCs were more 440 localized, matching the more compact dendritic structure of stellate cells. This observation is in 441 line with previous findings that the excitatory identity controls the survival and wiring of local 442 interneurons 65,66 . We suggest that the difference in the morphology of SOM + neurons between 443 these two cortical areas might be a result of the difference in dendritic organization of the 444 targeted excitatory neurons. Consistent with this, in S1 L5 where the principal excitatory cells 445 are pyramidal, their inhibitory input comes from L5 MCs 45 . We hypothesize that the reshaping 446 of excitatory neurons' apical dendrites in S1 L4 during development, which depends on the 447 sensory input 7 , could be followed by the corresponding reshaping of SOM + neurons. It will be 448 interesting to test whether this MC/pyramidal and NMC/stellate paring exists in other cortical 449 areas and other species. On the other hand, while we found that SOM + cells receive inputs from local excitatory neurons 452 in S1 L4, in agreement with previous studies 14,23 , we did not detect connections from L4 PYRs 453 to L4 MCs in V1. SOM + MCs in other layers are known to receive facilitating excitatory inputs 454 from local principal neurons in both S1 67 In addition to examining the connectivity among PYRs and MCs in V1 L4, we tested the 463 connectivity of BCs, another major cell type in L4 (Fig. S10). We found that BCs in L4 followed 464 the same connectivity rules as described for basket cells in other layers 28 and in younger 465 animals 48 : BCs inhibit other BCs, MCs, and pyramidal cells, and are inhibited by MCs and 466 excited by PYRs. All of these connection patterns have also been reported in S1 L4 in young 467 mice 23 , suggesting that the circuitry wiring involving PV + cells is roughly conserved between 468 these two areas and across age. 469 470 Finally, we found very low connection probability between PYRs in V1 L4, which was consistent 471 with the findings in V1 L2/3 and V1 L5 of adult mice 28 , but much lower than what was reported 472 in young animals 70,71 . We directly showed that this difference in connection probability among 473 excitatory neurons is due to the age of the animal (Fig. S11). 474 475 Summary 476 In conclusion, we confirmed the difference in morphology of L4 principal cells and revealed a 477 difference in morphology of L4 SOM + interneurons in V1 and S1 of adult mice. In each area, the 478 morphology of SOM + interneurons matched that of the excitatory neurons, suggesting that one 479 of them might adapt to another. Furthermore, we found differences in the connections from 480 excitatory neurons to SOM + interneurons, suggesting a different functional role of SOM + 481 interneurons in different cortical areas. In addition, we found that there is no one-to-one match 482 between the morphological and the transcriptomic types of SOM + interneurons, highlighting the 483 need of multi-modal profiling of cell types in the neocortex. Our results support the view that 484  Fig. S1) were performed using wild-type (n=24), Viaat-Cre/Ai9 (n=47), 502 Scnna1-Cre/Ai9 (n=5 for V1 and n=5 for S1), SOM-Cre/Ai9 (n=14 for V1 and n=19 for S1), VIP- We followed the pipeline that we recently benchmarked in Laturnus et al. 43 . As predictors for 653 pairwise classification we used morphometric statistics and density maps. Due to the very high 654 dimensionality of the density maps, we reduced them to 10 principal components (for cross-655 validation, PCA was computed on each outer-loop training set separately, and the same 656 transformation was applied to the corresponding outer-loop test set). This makes the final 657 feature dimensionality equal to 36. 658 659 For classification, we used logistic regression regularized with elastic net. Regularization 660 parameter alpha was fixed to 0.5, which is giving equal weights to the lasso and ridge penalties. 661 We used nested cross-validation to choose the optimal value of the regularization parameter 662 lambda and to obtain an unbiased estimate of the performance. The inner loop was performed 663 using the civisanalytics Python wrapper around the glmnet library 77 that does K-fold 664 31 cross-validation internally. We used 5 folds for the inner loop. We kept the default setting which 665 uses the maximal value of lambda with cross-validated loss within one standard error of the 666 lowest loss (lambda_best) to make the test-set predictions: 667 LogitNet(alpha=0.5, n_splits=5, random_state=42) 668 Note that the default behavior of glmnet is to standardize all predictors. The outer loop was 10 669 times repeated stratified 5-fold cross-validation, as implemented in scikit-learn by 670 RepeatedStratifiedKFold(n_splits=5, n_repeats=10, random_state=43) 671 Outer-loop performance was assessed via test-set accuracy. The resting membrane potential and the input resistance were computed differently for the 691 standard patch-clamp/morphology recordings and for the Patch-seq recordings, because of the 692 differences in the stimulation protocol between these two sets of experiments. In the Patch-seq 693 experiments, the current clamp value before each current stimulation was fixed at 0 pA for all 694 cells. Consequently, we computed the resting membrane potential as the median membrane 695 voltage before stimulation onset. Input resistance for each hyperpolarizing stimulation was 696 calculated as the ratio of the maximum voltage deflection to the corresponding current value. 697 We took the median of all hyperpolarizing currents as the final input resistance value. In 698 contrast, in the standard patch-clamp experiments, the current clamp before current stimulation 699 was not always fixed at 0 pA. For that reason we used linear regression (for robustness, random 700 sample consensus regression, as implemented in scikit-learn) of the steady state 701 membrane voltage onto the current stimulation value to compute the input resistance 702 (regression slope) and the resting membrane potential (regression intercept) (Fig. S1D). For this 703 we used five highest hyperpolarizing currents (if there were fewer than five, we used those 704 available). 705 706 To estimate the rheobase (minimum current needed to elicit any spikes), we used robust 707 regression of the spiking frequency onto the stimulation current using the five lowest 708 depolarizing stimulation currents with non-zero spike count (if there fewer than five, we used 709 those available) (Fig. S1D). The point where the regression line crosses the x-axis gives the 710 rheobase estimate. We restricted the rheobase estimate to be between the highest current 711 clamp value eliciting no spikes and the lowest current clamp value eliciting at least one spike. In 712 the rare cases when the regression line crossed the x-axis outside of this interval, the nearest 713 edge of the interval was taken instead as the rheobase estimate. 714 715 33 The action potential (AP) threshold, AP amplitude, AP width, afterhyperpolarization (AHP), 716 afterdepolarization (ADP), and the first spike latency were computed as illustrated in Fig. S1C, 717 using the very first spike fired by the neuron. AP width was computed at the half AP height. 718

719
The adaptation index (AI) is defined as the ratio of the second interspike interval to the first one 720 (Fig. S1B). We took the median over the five lowest depolarizing current stimulation that elicited 721 at least three spikes (if fewer than five were available, we used all of them). For the t-SNE visualization (Fig. 2B), we log-transformed the AI values because this feature had 732 a strongly right-skewed distribution (Fig. S1). We also excluded ADP and latency; ADP 733 because it was equal to zero for most neurons and rare cells with non-zero values appeared as 734 isolated subpopulations in the t-SNE representation, and latency because it had high outliers 735 among the FS types, also yielding isolated subpopulations. The remaining 11 features were z-736 scored and exact (non-approximate) t-SNE was run with perplexity 15 and random initialisation 737 with seed 42 using scikit-learn implementation: 738 TSNE(perplexity=15, method='exact', random_state=42) 739 740