Single-Cell RNA Sequencing.

Mice were collected at time points E13.5, E14.5, E17.5, E18.5, P1, and P5. Adrenal glands were isolated and dissociated into single cells. Adrenal glands were collected in PBS containing 1 mg/mL collagenase/dispase (Roche) and sheared using a 1-mL pipette on ice until single-cell suspension was observed. The enzymatic digestion was stopped by the addition of basic medium (Advanced DMEM/F12 [Gibco] containing N2 [1×; Gibco], B27 [1×; Gibco], Hepes [10 mM; Gibco], penicillin–streptomycin [100 U/mL; Gibco], Glutamax [1×; Gibco], N-acetyl cysteine [1.25 mM; Sigma]), and the suspension was centrifuged at 250 × g before being resuspended in PBS. All centrifugation steps were done at 4 °C, with the samples kept on ice for the duration of the experiment.

Single-cell suspensions were stained using 1 µM DAPI (Sigma). Viable single cells were sorted based on forward/side-scatter properties and DAPI staining (gating strategy shown in SI Appendix, Fig. S7) using FACS (FACSJazz; BD Biosciences) into 384-well plates (Bio-Rad) containing 10 µL of mineral oil (Sigma) and 50 nL of RT primers. Libraries were prepared using the Sort-seq method (18) and subjected to paired-end sequencing with 75-bp read length using the Illumina NextSeq500 sequencer.

Sequencing data were processed using the Sharq pipeline as described (49). Mapping was performed using STAR, version 2.6.1, on the Genome Reference Consortium Mouse Build 38. Read assignment was with featureCounts, version 1.5.2, using a gene annotation based on GENCODE, version M14.

Transcripts mapping to the mitochondrial genome and external RNA controls (ERCCs) were removed from all cells. A total of 410 cells forming two distinct clusters of erythroid lineages was identified based on high levels of hemoglobin complex genes and removed. To avoid contamination, all cells expressing more than 25% erythroid marker genes (SI Appendix, Fig. S8 and Table S2) were also removed. Finally, cells with mitochondrial-encoded transcripts over 50% of nuclear ones were removed, as well as cells with less than 1,000 nuclear-encoded transcripts or 334 genes. Subsequently, unique transcript counts were normalized using sctranform (50) and analyzed using the Seurat R package, version 3.1.5 (51). A Pearson residuals variance of 2 was used as a cutoff for identification of the variable genes. The following genes were removed from the list of variable genes in order to avoid biases in the cell clustering: genes associated with sex (Xist, Tsix, and Y chromosome-specific genes), cell cycle phase, dissociation stress (heat shock and chaperone proteins according to GO:0006986), and activity (ribosomal protein genes according to GO:0022626), as described before (SI Appendix, Table S3) (52). The first 25 principal components were used to calculate dimensionality reduction using the UMAP technique. Clustering was performed using the same principal components, a resolution of 2 and the Louvain algorithm.

The resulting 23 clusters were annotated on a single-cell basis using SingleR 1.2.4 (20) with the MouseRNAseqData reference dataset, consisting of 358 bulk RNA-seq samples of sorted cell populations (SI Appendix, Fig. S1A). Marker genes were used for further refining the cell annotations (SI Appendix, Fig. S1B). The 23 clusters were partitioned in five groups based on correlation. Differential expression analysis was performed with the FindMarkers Seurat (51) function for groups, clusters, and their combination using the Wilcoxon test with 1.5-fold change expression cutoff and 5% Bonferroni multiple testing corrected statistical significance cutoff. All genes passing the above criteria are shown in Dataset S3. Each marker gene was assigned to a unique group/cluster, based on higher overall expression.

Heatmaps were generated by heatmaply, version 1.1 (53). Columns represent cells ordered by cluster and developmental stage. Rows represent, in order, the top 20 differentially expressed genes of each group, followed by the top 20 differentially expressed genes in each of the clusters. The mean-centered expression of each gene–cell combination is represented in color according to the legend. Dot plots show gene expression changes across groups or clusters shown in the y axis. The size of the dot encodes the percentage of cells expressing each gene and the color encodes the log differential expression levels, as calculated by the DotPlot function in Seurat.

For the cortex and medulla group, trajectory graphs and pseudotime projections were generated with monocle3, version 0.2.2 (23), using code described in the SeuratWrappers package of the Seurat developers (https://github.com/satijalab/seurat-wrappers) and specifying E13.5 and E14.5 SCPs as the root node.