Cells with over 200,000 mapped mRNA reads were clustered iteratively based on WGCNA gene modules as described previously (Tasic et al., 2016; Thomsen et al. 2016). Due to strong temporal gene expression variation across time points, we clustered independently for cells at each age group. Within each age group, we assessed whether gene modules were represented by one batch. Modules that showed significant batch effect, or corresponded to cell cycle states, ribosome rRNAs, ribosomal proteins, immediate early response genes, mitochondrial genes, RNA binding, and mRNA quality (as assessed by mRNA mapping percentage of the single cells) were removed from downstream analysis, as our aim was to identify cell types independent of these factors. The cells were then clustered based on eigengenes from remaining gene modules using hierarchical clustering with “Ward’s” method. The number of clusters were selected dynamically so that every pair of clusters can be distinguished by at least 5 differentially expressed genes (2-fold change, adjusted Pvalue < 0.01) whose sum of –log10(adjusted Pvalue) is greater than 40 (each gene’s contribution capped at 20). We produced cluster heatmaps at
