Skip to main content

Thank you for visiting nature.com. You are using a browser version with limited support for CSS. To obtain the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles and JavaScript.

  • Article
  • Published:

The single-cell transcriptional landscape of mammalian organogenesis

Nature volume 566, pages 496–502 (2019)Cite this article

Subjects

Abstract

Mammalian organogenesis is a remarkable process. Within a short timeframe, the cells of the three germ layers transform into an embryo that includes most of the major internal and external organs. Here we investigate the transcriptional dynamics of mouse organogenesis at single-cell resolution. Using single-cell combinatorial indexing, we profiled the transcriptomes of around 2 million cells derived from 61 embryos staged between 9.5 and 13.5 days of gestation, in a single experiment. The resulting ‘mouse organogenesis cell atlas’ (MOCA) provides a global view of developmental processes during this critical window. We use Monocle 3 to identify hundreds of cell types and 56 trajectories, many of which are detected only because of the depth of cellular coverage, and collectively define thousands of corresponding marker genes. We explore the dynamics of gene expression within cell types and trajectories over time, including focused analyses of the apical ectodermal ridge, limb mesenchyme and skeletal muscle.

This is a preview of subscription content, access via your institution

Access options

Buy this article

  • Purchase on SpringerLink
  • Instant access to full article PDF
Buy now

Prices may be subject to local taxes which are calculated during checkout

Fig. 1: sci-RNA-seq3 enables profiling of 2,072,011 cells from 61 mouse embryos across 5 developmental stages in a single experiment.
Fig. 2: Identifying the major cell types of mouse organogenesis.
Fig. 3: Identification and characterization of epithelial cell subtypes and the limb AER.
Fig. 4: Characterization of ten major developmental trajectories present during mouse organogenesis.
Fig. 5: UMAP visualization of individual major trajectories.
Fig. 6: Resolving cellular trajectories in myogenesis.

Similar content being viewed by others

Data availability

The sci-RNA-seq3 protocol and all data have been made freely available, including through a cell-type wiki to facilitate their ongoing annotation by the research community (http://atlas.gs.washington.edu/mouse-rna). The data generated in this study can be downloaded in raw and processed forms from the NCBI Gene Expression Omnibus under accession number GSE119945.

References

  1. Kojima, Y., Tam, O. H. & Tam, P. P. L. Timing of developmental events in the early mouse embryo. Semin. Cell Dev. Biol. 34, 65–75 (2014).

    Article  CAS  Google Scholar 

  2. Tam, P. P. L. & Loebel, D. A. F. Gene function in mouse embryogenesis: get set for gastrulation. Nat. Rev. Genet. 8, 368–381 (2007).

    Article  CAS  Google Scholar 

  3. Dickinson, M. E. et al. High-throughput discovery of novel developmental phenotypes. Nature 537, 508–514 (2016).

    Article  CAS  Google Scholar 

  4. Meehan, T. F. et al. Disease model discovery from 3,328 gene knockouts by The International Mouse Phenotyping Consortium. Nat. Genet. 49, 1231–1238 (2017).

    Article  CAS  Google Scholar 

  5. Wagner, D. E. et al. Single-cell mapping of gene expression landscapes and lineage in the zebrafish embryo. Science 360, 981–987 (2018).

    Article  ADS  CAS  Google Scholar 

  6. Briggs, J. A. et al. The dynamics of gene expression in vertebrate embryogenesis at single-cell resolution. Science 360, eaar5780 (2018).

    Article  Google Scholar 

  7. Farrell, J. A. et al. Single-cell reconstruction of developmental trajectories during zebrafish embryogenesis. Science 360, eaar3131 (2018).

    Article  Google Scholar 

  8. Mayer, C. et al. Developmental diversification of cortical inhibitory interneurons. Nature 555, 457–462 (2018).

    Article  ADS  CAS  Google Scholar 

  9. Lescroart, F. et al. Defining the earliest step of cardiovascular lineage segregation by single-cell RNA-seq. Science 359, 1177–1181 (2018).

    Article  ADS  CAS  Google Scholar 

  10. Han, X. et al. Mapping the Mouse Cell Atlas by Microwell-seq. Cell 172, 1091–1107 (2018).

    Article  CAS  Google Scholar 

  11. The Tabula Muris Consortium, Quake, S. R., Wyss-Coray, T. & Darmanis, S. Transcriptomic characterization of 20 organs and tissues from mouse at single cell resolution creates a Tabula Muris. Preprint at https://www.biorxiv.org/content/10.1101/237446v2 (2018).

  12. Amini, S. et al. Haplotype-resolved whole-genome sequencing by contiguity-preserving transposition and combinatorial indexing. Nat. Genet. 46, 1343–1349 (2014).

    Article  CAS  Google Scholar 

  13. Adey, A. et al. In vitro, long-range sequence information for de novo genome assembly via transposase contiguity. Genome Res. 24, 2041–2049 (2014).

    Article  CAS  Google Scholar 

  14. Cusanovich, D. A. et al. Multiplex single cell profiling of chromatin accessibility by combinatorial cellular indexing. Science 348, 910–914 (2015).

    Article  ADS  CAS  Google Scholar 

  15. Vitak, S. A. et al. Sequencing thousands of single-cell genomes with combinatorial indexing. Nat. Methods 14, 302–308 (2017).

    Article  CAS  Google Scholar 

  16. Ramani, V. et al. Massively multiplex single-cell Hi-C. Nat. Methods 14, 263–266 (2017).

    Article  CAS  Google Scholar 

  17. Cao, J. et al. Comprehensive single-cell transcriptional profiling of a multicellular organism. Science 357, 661–667 (2017).

    Article  ADS  CAS  Google Scholar 

  18. Mulqueen, R. M. et al. Highly scalable generation of DNA methylation profiles in single cells. Nat. Biotechnol. 36, 428–431 (2018).

    Article  CAS  Google Scholar 

  19. Cao, J. et al. Joint profiling of chromatin accessibility and gene expression in thousands of single cells. Science 361, 1380–1385 (2018).

    Article  ADS  CAS  Google Scholar 

  20. Rosenberg, A. B. et al. Single-cell profiling of the developing mouse brain and spinal cord with split-pool barcoding. Science 360, 176–182 (2018).

    Article  ADS  CAS  Google Scholar 

  21. La Manno, G. et al. RNA velocity of single cells. Nature 560, 494–498 (2018).

    Article  ADS  Google Scholar 

  22. Wolock, S. L., Lopez, R. & Klein, A. M. Scrublet: computational identification of cell doublets in single-cell transcriptomic data. Preprint at https://www.biorxiv.org/content/10.1101/357368v1 (2018).

    Article  Google Scholar 

  23. Qiu, X. et al. Reversed graph embedding resolves complex single-cell developmental trajectories. Nat. Methods 14, 979–982 (2017).

    Article  CAS  Google Scholar 

  24. Yang, A. et al. p63 is essential for regenerative proliferation in limb, craniofacial and epithelial development. Nature 398, 714–718 (1999).

    Article  ADS  CAS  Google Scholar 

  25. McQualter, J. L., Yuen, K., Williams, B. & Bertoncello, I. Evidence of an epithelial stem/progenitor cell hierarchy in the adult mouse lung. Proc. Natl Acad. Sci. USA 107, 1414–1419 (2010).

    Article  ADS  CAS  Google Scholar 

  26. Cichorek, M., Wachulska, M., Stasiewicz, A. & Tymińska, A. Skin melanocytes: biology and development. Postepy Dermatol. Allergol. 30, 30–41 (2013).

    Article  Google Scholar 

  27. Tomihari, M., Hwang, S.-H., Chung, J.-S., Cruz, P. D. Jr & Ariizumi, K. Gpnmb is a melanosome-associated glycoprotein that contributes to melanocyte/keratinocyte adhesion in a RGD-dependent fashion. Exp. Dermatol. 18, 586–595 (2009).

    Article  CAS  Google Scholar 

  28. Varjosalo, M. & Taipale, J. Hedgehog: functions and mechanisms. Genes Dev. 22, 2454–2472 (2008).

    Article  CAS  Google Scholar 

  29. Strähle, U., Lam, C. S., Ertzer, R. & Rastegar, S. Vertebrate floor-plate specification: variations on common themes. Trends Genet. 20, 155–162 (2004).

    Article  Google Scholar 

  30. Holmes, G. P. et al. Distinct but overlapping expression patterns of two vertebrate slit homologs implies functional roles in CNS development and organogenesis. Mech. Dev. 79, 57–72 (1998).

    Article  CAS  Google Scholar 

  31. Akle, V. et al. F-spondin/spon1b expression patterns in developing and adult zebrafish. PLoS ONE 7, e37593 (2012).

    Article  ADS  CAS  Google Scholar 

  32. Zeisel, A. et al. Molecular architecture of the mouse nervous system. Cell 174, 999–1014 (2018).

    Article  CAS  Google Scholar 

  33. Hartman, B. H., Durruthy-Durruthy, R., Laske, R. D., Losorelli, S. & Heller, S. Identification and characterization of mouse otic sensory lineage genes. Front. Cell. Neurosci. 9, 79 (2015).

    Article  Google Scholar 

  34. Szenker-Ravi, E. et al. RSPO2 inhibition of RNF43 and ZNRF3 governs limb development independently of LGR4/5/6. Nature 557, 564–569 (2018).

    Article  ADS  CAS  Google Scholar 

  35. Cai, X. et al. Tbx20 acts upstream of Wnt signaling to regulate endocardial cushion formation and valve remodeling during mouse cardiogenesis. Development 140, 3176–3187 (2013).

    Article  CAS  Google Scholar 

  36. Miller, R. A., Christoforou, N., Pevsner, J., McCallion, A. S. & Gearhart, J. D. Efficient array-based identification of novel cardiac genes through differentiation of mouse ESCs. PLoS ONE 3, e2176 (2008).

    Article  ADS  Google Scholar 

  37. Petit, F., Sears, K. E. & Ahituv, N. Limb development: a paradigm of gene regulation. Nat. Rev. Genet. 18, 245–258 (2017).

    Article  CAS  Google Scholar 

  38. Guo, Q., Loomis, C. & Joyner, A. L. Fate map of mouse ventral limb ectoderm and the apical ectodermal ridge. Dev. Biol. 264, 166–178 (2003).

    Article  CAS  Google Scholar 

  39. Lewandoski, M., Sun, X. & Martin, G. R. Fgf8 signalling from the AER is essential for normal limb development. Nat. Genet. 26, 460–463 (2000).

    Article  CAS  Google Scholar 

  40. Gerdes, J., Schwab, U., Lemke, H. & Stein, H. Production of a mouse monoclonal antibody reactive with a human nuclear antigen associated with cell proliferation. Int. J. Cancer 31, 13–20 (1983).

    Article  CAS  Google Scholar 

  41. Bergman, D., Halje, M., Nordin, M. & Engström, W. Insulin-like growth factor 2 in development and disease: a mini-review. Gerontology 59, 240–249 (2013).

    Article  CAS  Google Scholar 

  42. Trapnell, C. et al. The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nat. Biotechnol. 32, 381–386 (2014).

    Article  CAS  Google Scholar 

  43. McInnes, L., Healy, J. & Melville, J. UMAP: Uniform Manifold Approximation and Projection for dimension reduction. Preprint at https://arxiv.org/abs/1802.03426 (2018).

  44. Alexander Wolf, F. et al. Graph abstraction reconciles clustering with trajectory inference through a topology preserving map of single cells. Preprint at https://www.biorxiv.org/content/10.1101/208819v2 (2018).

  45. Braun, T. & Gautel, M. Transcriptional mechanisms regulating skeletal muscle differentiation, growth and homeostasis. Nat. Rev. Mol. Cell Biol. 12, 349–361 (2011).

    Article  CAS  Google Scholar 

  46. Comai, G., Sambasivan, R., Gopalakrishnan, S. & Tajbakhsh, S. Variations in the efficiency of lineage marking and ablation confound distinctions between myogenic cell populations. Dev. Cell 31, 654–667 (2014).

    Article  CAS  Google Scholar 

  47. Halperin-Barlev, O. & Kalcheim, C. Sclerotome-derived Slit1 drives directional migration and differentiation of Robo2-expressing pioneer myoblasts. Development 138, 2935–2945 (2011).

    Article  CAS  Google Scholar 

  48. Heimberg, G., Bhatnagar, R., El-Samad, H. & Thomson, M. Low dimensionality in gene expression data enables the accurate extraction of transcriptional programs from shallow sequencing. Cell Syst. 2, 239–250 (2016).

    Article  CAS  Google Scholar 

  49. Osterwalder, M. et al. Enhancer redundancy provides phenotypic robustness in mammalian development. Nature 554, 239–243 (2018).

    Article  ADS  CAS  Google Scholar 

  50. Dickel, D. E. et al. Ultraconserved enhancers are required for normal development. Cell 172, 491–499.e15 (2018).

    Article  CAS  Google Scholar 

  51. Kraft, K. et al. Deletions, inversions, duplications: engineering of structural variants using CRISPR/Cas in mice. Cell Rep. 4, S2211–S1247 (2015).

    Google Scholar 

  52. Buenrostro, J. D., Giresi, P. G., Zaba, L. C., Chang, H. Y. & Greenleaf, W. J. Transposition of native chromatin for fast and sensitive epigenomic profiling of open chromatin, DNA-binding proteins and nucleosome position. Nat. Methods 10, 1213–1218 (2013).

    Article  CAS  Google Scholar 

  53. Renaud, G., Stenzel, U., Maricic, T., Wiebe, V. & Kelso, J. deML: robust demultiplexing of Illumina sequences using a likelihood-based approach. Bioinformatics 31, 770–772 (2015).

    Article  CAS  Google Scholar 

  54. Dobin, A. et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, 15–21 (2013).

    Article  CAS  Google Scholar 

  55. Anders, S., Pyl, P. T. & Huber, W. HTSeq–a Python framework to work with high-throughput sequencing data. Bioinformatics 31, 166–169 (2014).

    Article  Google Scholar 

  56. Wolf, F. A., Angerer, P. & Theis, F. J. SCANPY: large-scale single-cell gene expression data analysis. Genome Biol. 19, 15 (2018).

    Article  Google Scholar 

  57. Zheng, G. X. Y. et al. Massively parallel digital transcriptional profiling of single cells. Nat. Commun. 8, 14049 (2017).

    Article  ADS  CAS  Google Scholar 

  58. Wolock, S. L., Lopez, R. & Klein, A. M. Scrublet: computational identification of cell doublets in single-cell transcriptomic data. Preprint at https://www.biorxiv.org/content/10.1101/357368v1 (2018).

  59. Pliner, H. et al. Chromatin accessibility dynamics of myogenesis at single cell resolution. Preprint at https://www.biorxiv.org/content/10.1101/155473v1 (2017).

  60. Kuleshov, M. V. et al. Enrichr: a comprehensive gene set enrichment analysis web server 2016 update. Nucleic Acids Res. 44, W90–W97 (2016).

    Article  CAS  Google Scholar 

  61. Levine, J. H. et al. Data-driven phenotypic dissection of AML reveals progenitor-like cells that correlate with prognosis. Cell 162, 184–197 (2015).

    Article  CAS  Google Scholar 

  62. Mao, Q., Wang, L., Tsang, I. & Sun, Y. Principal graph and structure learning based on reversed graph embedding. IEEE Trans. Pattern Anal. Mach. Intell. 39, 2227–2241 (2016).

    Article  Google Scholar 

  63. Mao, Q., Yang, L., Wang, L., Goodison, S. & Sun, Y. SimplePPT: a simple principal tree algorithm. In Proc. 2015 SIAM International Conference on Data Mining (eds Venkatasubramanian, S. & Ye, J.) 792–800 (SIAM, 2015).

  64. Moran, P. A. P. Notes on continuous stochastic phenomena. Biometrika 37, 17–23 (1950).

    Article  MathSciNet  CAS  Google Scholar 

  65. Li, D. et al. Formation of proximal and anterior limb skeleton requires early function of Irx3 and Irx5 and is negatively regulated by Shh signaling. Dev. Cell 29, 233–240 (2014).

    Article  CAS  Google Scholar 

Download references

Acknowledgements

We thank members of the Shendure and Trapnell labs, especially D. Cusanovich, R. Daza, G. Findlay, A. McKenna, H. Pliner and V. Ramani, as well as L. McInnes, D. Beier, N. Ahituv and S. Tapscott for helpful discussions and feedback; M. Zager for major contributions to the website; R. Hunter, and R. Rualo at the Transgenic Resources Program of University of Washington and N. Brieske and A. Stiege at the Max Planck Institute for Molecular Genetics for their assistance; S. Geuer for the Fndc3a probe. M.S. was supported by a grant from the Deutsche Forschungsgemeinschaft (SP1532/3-1). This work was funded by the Paul G. Allen Frontiers Group (Allen Discovery Center grant to J.S. and C.T.), grants from the NIH (DP1HG007811 and R01HG006283 to J.S.; DP2 HD088158 to C.T.), the W. M. Keck Foundation (to C.T. and J.S.). J.S. is an Investigator of the Howard Hughes Medical Institute.

Reviewer information

Nature thanks Alistair Forrest, Peter Sims, Patrick Tam and the other anonymous reviewer(s) for their contribution to the peer review of this work.

Author information

Author notes

  1. These authors contributed equally: Junyue Cao, Malte Spielmann

Authors and Affiliations

  1. Department of Genome Sciences, University of Washington, Seattle, WA, USA

    Junyue Cao, Malte Spielmann, Xiaojie Qiu, Xingfan Huang, Andrew J. Hill, Cole Trapnell & Jay Shendure

  2. Molecular and Cellular Biology Program, University of Washington, Seattle, WA, USA

    Junyue Cao & Xiaojie Qiu

  3. Department of Computer Science, University of Washington, Seattle, WA, USA

    Xingfan Huang

  4. Max Planck Institute for Molecular Genetics, RG Development & Disease, Berlin, Germany

    Daniel M. Ibrahim & Stefan Mundlos

  5. Institute for Medical and Human Genetics, Charité Universitätsmedizin Berlin, Berlin, Germany

    Daniel M. Ibrahim & Stefan Mundlos

  6. Illumina, San Diego, CA, USA

    Fan Zhang, Lena Christiansen & Frank J. Steemers

  7. Brotman Baty Institute for Precision Medicine, Seattle, WA, USA

    Cole Trapnell & Jay Shendure

  8. Allen Discovery Center for Cell Lineage Tracing, Seattle, WA, USA

    Cole Trapnell & Jay Shendure

  9. Howard Hughes Medical Institute, Seattle, WA, USA

    Jay Shendure

Authors
  1. Junyue Cao
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Malte Spielmann
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Xiaojie Qiu
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Xingfan Huang
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Daniel M. Ibrahim
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Andrew J. Hill
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Fan Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Stefan Mundlos
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Lena Christiansen
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. Frank J. Steemers
    View author publications

    You can also search for this author in PubMed Google Scholar

  11. Cole Trapnell
    View author publications

    You can also search for this author in PubMed Google Scholar

  12. Jay Shendure
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

J.C. developed techniques and performed sci-RNA-seq3 experiments with assistance from M.S., F.Z., L.C. and F.J.S. M.S. performed embryo collection and in situ hybridization validations with assistance from D.M.I. and S.M. J.C. and C.T. performed computational analyses with assistance from M.S., X.Q. and A.J.H. X.Q. and C.T. developed Monocle 3. X.H. developed the website with assistance from J.C. J.S. and C.T. supervised the project. J.S., C.T., J.C. and M.S. conceived the project and wrote the manuscript.

Corresponding authors

Correspondence to Cole Trapnell or Jay Shendure.

Ethics declarations

Competing interests

L.C., F.Z. and F.J.S. declare competing financial interests in the form of stock ownership and paid employment by Illumina. One or more embodiments of one or more patents and patent applications filed by Illumina may encompass the methods, reagents and data disclosed in this manuscript. Some work in this study may be related to technology described in the following exemplary published patent applications: WO2010/0120098 and WO2011/0287435.

Additional information

Publisher’s note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Extended data figures and tables

Extended Data Fig. 1 Performance and quality-control-related analyses for sci-RNA-seq3.

a, Comparison of fixation conditions in human HEK-293T cells. PFA-fixed nuclei yielded the highest numbers of UMIs. Cell number: n = 21 for fresh nuclei, 17 for frozen nuclei, 32 for PFA-fixed cells and 31 for PFA-fixed nuclei. b, Tn5 transposomes loaded only with N7 adaptor (cell number, n = 13 cells) increased UMI counts by over 50%, relative to the standard Nextera Tn5 (cell number, n = 11), in human HEK-293T cells. c, Bar plot showing the number of reverse transcription wells used for each of 61 mouse embryos. d, Histogram showing the distribution of raw sequencing reads from each PCR well in sci-RNA-seq3. e, Scatter plot of mouse (NIH/3T3) versus human (HEK-293T) UMI counts per cell. f, g, Box plot showing the number of UMIs and purity (proportion of reads mapping to the expected species) per cell from HEK-293T (cell number n = 7,943) and NIH/3T3 cells (cell number, n = 10,914). At a sequencing depth of 23,207 reads per cell, we observed a median of 5,461 UMIs per HEK-293T cell and 5,087 UMIs per NIH/3T3 cell, with 3.9% and 2.9% of reads per cell mapping to incorrect species, respectively. h, Box plot comparing the number of UMIs per cell (downsampled to 20,000 raw reads per cell) for sci-RNA-seq3 (cell number, n = 689 for HEK-293T and 997 for NIH/3T3) versus sci-RNA-seq (cell number, n = 47 for HEK-293T and 120 for NIH/3T3). i, Correlation (Pearson’s correlation) between gene expression measurements in aggregated profiles of HEK-293T from sci-RNA-seq3 nuclei versus sci-RNA-seq cells. j, Scatter plot showing correlation between number of reverse transcription wells used and number of cells recovered per embryo. k, Box plot showing the number of genes and UMIs detected per cell. l, Box plot showing the number of UMIs detected per cell from embryos across five developmental stages. Cell number: n = 152,120 for E9.5; 378,427 for E10.5; 615,908 for E11.5; 475,047 for E12.5; 437,150 for E13.5. m, Histogram showing the distribution of the cell doublet score for the actual mouse embryo data versus doublets stimulated by Scrublet. n, Scatter plot of the number of cells profiled per reverse transcription well and the detected doublet-cell ratio. Blue line shows the linear regression. The detected doublet-cell rate was modestly correlated with number of cells profiled per well during reverse transcription (Spearman’s ρ = 0.35). o, Scatter plot of unique reads aligning to Xist (female-specific) versus chrY transcripts (male-specific) per mouse embryo. Sex assignments of individual embryos inferred from these data. p, Bar plot showing the number of male and female embryos profiled at each developmental stage. q, t-SNE of the aggregated transcriptomes of single cells derived from each of 61 mouse embryos results in 5 tightly clustered groups perfectly matching their developmental stages (embryo number, n = 61). r, Pseudotime trajectory of pseudobulk RNA-seq profiles of mouse embryos (embryo number, n = 61); identical to Fig. 1c, but coloured by pseudotime. s, The E10.5 embryos were ordered by pseudotime. The 3 earliest versus 3 latest (in pseudotime) E10.5 embryos are shown in photographs, and appear to potentially be morphologically distinct. Notably, the distinct colouring of E10.5 embryos positioned earlier versus later in developmental pseudotime is potentially due to different levels of haemoglobin. For all box plots: thick horizontal lines, medians; upper and lower box edges, first and third quartiles, respectively; whiskers, 1.5 times the interquartile range; circles, outliers.

Extended Data Fig. 2 Identifying the major cell types and cell-composition dynamics during mouse organogenesis.

a–e, t-SNE visualization of mouse embryo cells from different developmental stages, with sampling 10,000 cells per stage and colouring by embryo ID: E9.5 (a), E10.5 (b), E11.5 (c), E12.5 (d), E13.5 (e). We consistently observe that cells derived from independent embryos at the same time point are similarly distributed. f, The same t-SNE as Fig. 2a is shown, with subsets of cells highlighted. The first panel only shows cells from E9.5 embryos, and cells from subsequent developmental stages are progressively added. g, Box plot showing the number of UMIs detected per cell for major cell types (cell number n for each cell type is listed in Supplementary Table 3). Thick horizontal lines, medians; upper and lower box edges, first and third quartiles, respectively; whiskers, 1.5 times the interquartile range; circles, outliers. h, t-SNE visualization of a randomly sampled 100,000 cells coloured by expression level of Hbb-bh1 (top) or Fndc3c1 (bottom). ‘High’ indicates cells with UMI count for Hbb-bh1 >3 or Fndc3c1 >1. i, Bar plot showing the number of marker genes in each major cell type, defined as differentially expressed genes (5% FDR) with a >twofold (green) or >fivefold (red) expression difference between first- and second-ranked cell types. j, Left, t-SNE visualization of a randomly sampled 100,000 cells coloured by expression level of Shh (top) or Tox2 (bottom). Right, WISH images of Shh (top) or Tox2 (bottom) in embryos. n = 5. ‘High’ indicates cells with UMI count for Shh> 0 or Tox2 >1. Arrow, site of gene expression. k, Bar plot showing the number of cells profiled for each cell type, split by development stage. l, Heat map showing the estimated relative number of each cell type (rows) in 61 mouse embryos (columns). An estimate of the absolute cell number per cell type per embryo was calculated by multiplying the proportion that cell type contributed to a given embryo by the estimated total number of cells at that development stage. For presentation, these estimates are normalized in each row by the maximum estimated cell count for that cell type across all 61 embryos. Embryos are sorted left-to-right by developmental pseudotime. m, Line plot showing the estimated relative cell numbers for primitive erythroid and definitive erythroid lineages, calculated as in l. Dashed lines show relative expression of marker genes for primitive erythroid (Hbb-bh1) and definitive erythroid (Hbb-bs) major cell types. Data points for individual embryos were ordered by development pseudotime and smoothed by the LOESS method.

Extended Data Fig. 3 Louvain clustering and t-SNE visualization of subclusters of the each of 38 major cell types.

As cell-type heterogeneity was readily apparent within many of the 38 clusters shown in Fig. 2a, we adopted an iterative strategy, repeating Louvain clustering on each main cell type to identify subclusters. After subclusters dominated by 1 or 2 embryos were removed and highly similar subclusters were merged, a total of 655 subclusters were identified. (also termed ‘subtypes’ to distinguish them from the 38 major cell types identified by the initial clustering). Cell number, n, for each cell type is listed in Supplementary Table 3.

Extended Data Fig. 4 Analysis of cell subtypes during mouse organogenesis.

a, t-SNE visualization of all cells (top plot, n = 2,026,641) and downsampled subset of high-quality cells (bottom plot, n = 50,000, UMI > 400), coloured by Louvain cluster IDs from Fig. 2a. b, t-SNE visualization of all endothelial cells (top plot, n = 35,878) and those from the downsampled subset (bottom plot, n = 1,173), coloured by Louvain cluster ID computed on the basis of the 35,878 endothelial cells. c, d, t-SNE visualization of the downsampled subset of 50,000 cells (c), and 1,173 endothelial cells (d), coloured by Louvain cluster ID computed on the basis of sampled cells only. The number of clusters and subclusters identified with the same parameters drops from 38 (a, bottom plot) to 27 (c) and 16 (b, bottom plot) to 12 (c), respectively. e, Histogram showing the distribution of subclusters with respect to cell number (median 1,869; range 51–65,894). f, Histogram showing the distribution of subclusters with respect to the number of contributing embryos (>5 cells to qualify as a contributor). g, Histogram showing the distribution of subclusters with respect to the ratio of cells derived from the most highly contributing embryo. h, Histogram showing the distribution of subclusters with respect to the ratio of doublet cells detected by Scrublet. i, Histogram showing the distribution of subclusters with respect to the number of marker genes (at least twofold (blue)- or fivefold (red)-higher expression when compared with the second-highest expressing cell subtype within the same main cluster; 5% FDR). Out of 655 subclusters, 644 (98%) have at least 1 such gene marker with a twofold difference, and 441 of 655 (67%) have at least 1 such marker with a fivefold difference. j, t-SNE visualization of subcluster-specific marker expression (for example, cell number n = 74,651): Calb1 (left), Nox3 (middle) and Tex14 (right) are gene markers for three endothelial subclusters. ‘High’ indicates cells with UMI count for Calb1 >0, Nox3 >0 or Tex14 >1. k, Cumulative histogram showing how many subtypes (out of a total of 572 non-doublet-artefact subtypes) can be distinguished from all other subtypes on the basis of 1 or several markers and >fourfold expression differences (see also Methods, Supplementary Table 5).

Extended Data Fig. 5 Cell-type correlation analysis between single-cell mouse atlases.

a, Cell-type correlation analysis (Methods) matched cell types between independently generated and annotated analyses of the adult mouse kidney (sci-RNA-seq component of sci-CAR19 (rows) versus Microwell-seq10 (columns)). All cell types identified by sci-RNA-seq are shown, but we only show Microwell-seq cell types that are top matches for one or more sci-RNA-seq cell types. Colours correspond to beta values, normalized by the maximum beta value per row. b, Left, we compared our subtypes against 130 fetal cell types annotated in the MCA10 with cell-type correlation analysis, matching 96 MCA-defined cell types (rows) to 58 subtypes in our mouse embryo atlas (columns). Colours correspond to beta values, normalized by the maximum beta value per row. All MCA cell types with maximum beta of matched cell type >0.01 are shown (rows; n = 96), as are mouse embryo atlas cell types that are top matches for one or more displayed MCA cell types (columns; n = 58). Right, zoom-in of a subset of matches shown on the left. Cell-type annotations are from MCA (rows) or our study (columns; major cell-type annotation and subcluster ID). c, Box plot showing the ratio of cells from E13.5 for subclusters with (subcluster number, n = 58) versus without (subcluster number, n = 514) a matched cell type in the MCA. Thick horizontal lines, medians; upper and lower box edges, first and third quartiles, respectively; whiskers, 1.5 times the interquartile range; circles, outliers. d, Left, we compared our subtypes against 265 cell types annotated by a recent mouse brain cell atlas (BCA)32 with cell-type correlation analysis, matching 48 BCA-defined cell types (rows) to 68 subtypes in our data (columns). Colours correspond to beta values, normalized by the maximum beta value per row. All mouse embryo cell types with maximum beta of matched cell type >0.01 are shown (column; n = 68), as are BCA cell types that are top matches for 1 or more displayed mouse embryo cell types (rows; n = 48). Right, zoom-in of a subset of matches shown on the left. Cell-type annotations are from BCA (rows) or our study (columns; major cell cluster and subcluster ID).

Extended Data Fig. 6 Analysis of mouse epithelium, endothelium and limb-AER cells.

a, b, Dot plot showing expression of one selected marker gene per epithelial (a) or endothelial (b) subtype. Doublet-derived subclusters (2/29 epithelial subtypes and 5/16 endothelial subtypes) are excluded from these plots, but are shown in Fig. 3a and in c, respectively. The size of the dot encodes the percentage of cells within a cell type, and its colour encodes the average expression level. c, t-SNE visualization and marker-based annotation of endothelial cell subtypes (n = 35,878). d, Heat map showing smoothed pseudotime-dependent differential gene expression (510 genes at FDR of 1%) in AER cells, generated by a spline fitting with a generalized linear model (assuming gene expression following the negative binomial distribution) and scaled as a percentage of maximum gene expression. Each row indicates a different gene, and these are split into subsets that are activated (top), repressed (middle) or exhibit transient dynamics (bottom) between E9.5 and E13.5. e, f, Plots showing the −log-transformed q value and Enrichr-based combined score of enriched Reactome terms (e) and transcription factors (f) for genes with expression that significantly decreases in AER development. The top enriched pathway terms (Reactome2016) for significantly decreasing genes include cell-cycle progression (‘mitotic cell cycle’, q = 0.0002, one-sided Fisher exact test with multiple comparisons adjusted) and glucose metabolism (metabolism of carbohydrates, q = 0.0002, one-sided Fisher exact test with multiple comparisons adjusted). The top enriched transcription factors with targets from decreasing genes include pluripotent factors such as Isl1 (q < 1 × 10−5), Pou5f1 (q = 0.002, one-sided Fisher exact test with multiple comparisons adjusted) and Nanog (q = 0.003, one-sided Fisher exact test with multiple comparisons adjusted).

Extended Data Fig. 7 Characterizing cellular trajectories during limb mesenchyme differentiation.

a, UMAP 3D visualization of limb mesenchymal cells coloured by development stage (cell number, n = 26,559; left and right represent views from two directions). b, Heat map showing top differentially expressed genes between different developmental stages for limb mesenchyme cells. c, Bar plot showing the −log10-transformed adjusted P value (one-sided Fisher exact test with multiple comparisons adjusted) of enriched transcription factors for significantly upregulated genes during limb mesenchyme development. d, t-SNE visualization of limb mesenchyme cells coloured by forelimb (Tbx5+; cell number, n = 2,085) and hindlimb (Pitx1+; cell number, n = 1,885). Cells with no expression or expression of both in Tbx5 and Pitx1 are not shown. e, h, i, k, Each panel illustrates a different marker gene. Colours indicate UMI counts that have been scaled for library size, log-transformed, and then mapped to Z-scores to enable comparison between genes. Cells with no expression of a given marker are excluded to prevent overplotting. e, Hindlimb marker Pitx1 and forelimb marker Tbx5. f, Scatter plot showing the normalized expression of Pitx1 and Tbx5 in limb mesenchyme cells. Only cells in which Pitx1 and/or Tbx5 were detected are shown. g, Volcano plot showing the differentially expressed genes (FDR of 5%, one-sided likelihood ratio test with multiple comparisons adjusted, coloured red) between forelimb (cell number, n = 2,085) and hindlimb (cell number, n = 1,885). Top differentially expressed genes are labelled. x axis, log2-transformed fold change between forelimb and hindlimb for each gene; y axis, −log10-transformed q value from differential gene expression test. h, Same visualization as e, coloured by normalized gene expression of proximal/chondrocyte (Sox6 and Sox9), distal (Hoxd13 and Tfap2b), anterior (Pax9 and Alx4) or posterior (Hand2 and Shh) markers. Only cells with the gene marker expressed are plotted. i, Same visualization as e. First row, proximal limb markers Sox6 (which also marks chondrocytes) and Sox9. Second row, distal limb markers Hoxd13 and Tfap2b. Third row, anterior limb markers65 Pax9 and Alx4. Fourth row, posterior limb markers Shh and Hand2. j, In situ hybridization images of Hoxd13 in E9.5 to E13.5 embryos (n = 5). k, Same visualization as e, coloured by normalized gene expression of Cpa2. Only cells with positive UMI counts are shown. Values are log-transformed, standardized UMI counts. The expression pattern of Cpa2 within this trajectory led us to predict that it is a distal marker of the developing limb mesenchyme, similiar to Hoxd13. l, In situ hybridization images of Cpa2 in E10.5 and E11.5 embryos (n = 5. Arrow, site of gene expression. m, Modules of spatially restricted genes in the limbs. A total of 1,783 genes were clustered via hierarchical clustering. The dendrogram was cut into eight modules using the cutree function in R, and the aggregate expression of genes in each module was computed. Colours indicate aggregate UMI counts for each module that have been scaled for library size, log-transformed and then mapped to Z-scores to enable comparison between modules. Cells with no expression of a given module are excluded to prevent overplotting.

Extended Data Fig. 8 Characterization of ten major developmental trajectories present during mouse organogenesis.

a, Heat map showing the proportion of cells from each of the 38 major cell types assigned to each of the 12 PAGA algorithm-identified groups. We merged 2 groups corresponding to sensory neurons (12 and 3) and another 2 groups corresponding to blood cells (6 and 7) as each pair was closely located in UMAP space upon visual inspection, yielding the 10 supergroups shown in a similar heat map in Fig. 4b. b, Same as Fig. 4a, but with colours corresponding to the 38 major cell clusters. c, Area plot showing the estimated proportion (top) and estimated absolute number (bottom) of cells per embryo derived from each of the ten major cell trajectories from E9.5 to E13.5. Although the estimated number of cells per embryo in each of these supergroups increases exponentially, their proportions remain relatively stable, with the exception of hepatocytes which expand their contribution by nearly tenfold during this developmental window (from 0.3% at E9.5 to 2.8% at E13.5). d, UMAP 3D visualization of epithelial subtrajectories (as in Fig. 4c), coloured as per the epithelial subtypes shown in Fig. 3a.

Extended Data Fig. 9 UMAP visualization of the ten major cell trajectories.

We iteratively reanalysed each of the ten major trajectories, nearly all of which further resolved into multiple subtrajectories. The 10 major cell trajectories are visualized with UMAP (as in Fig. 5) but coloured: as per the 38 major cell clusters (top left), subcluster ID (top right), developmental stage (bottom left) and pseudotime (bottom right). The lines correspond to the principal graph learned by Monocle 3. These images are also available at http://atlas.gs.washington.edu/mouse-rna as manipulatable 3D renderings.

Extended Data Fig. 10 UMAP visualization of the 56 subtrajectories, coloured by development stage.

We further iteratively reanalysed and visualized with UMAP each of the 56 subtrajectories. Although Monocle 3 did not have access to these labels, the subtrajectories are highly consistent with developmental time (that is, cells ordered from E9.5 to E13.5). The lines correspond to the principal graph learned by Monocle 3.

Extended Data Fig. 11 UMAP visualization of the 56 subtrajectories, coloured by inferred pseudotime.

To orient each subtrajectory (same projections as Extended Data Fig. 10), we identified one or several starting points as focal concentrations of E9.5 cells, and then computed developmental pseudotime for cells present along various paths. The lines correspond to the principal graph learned by Monocle 3.

Extended Data Fig. 12 Gene dynamics in the myogenic trajectory.

a, Genes that are differentially expressed between the Myf5 path and the Myod path highlighted in Fig. 6. Cells along each path were compared using Monocle’s differentialGeneTest function. Pseudotimes along each path were scaled from 0 to 100 independently. The full model formula was ‘~path ∗ sm.ns(Pseudotime, df=3)’, whereas the reduced model was ‘~sm.ns(Pseudotime, df=3)’. Differentially expressed genes (FDR <1%, one-sided likelihood ratio test with multiple comparisons adjusted) were clustered via Ward’s method and visualized as a heat map via the pheatmap package. b, Pseudotemporal kinetics for selected genes involved in Robo–Slit signalling. Red indicates cells on the Myod1 path, while blue corresponds to the Myf5 path. Standardized expression scores for each gene on the original myogenic trajectory are shown next to the expression curves for each. Only cells with detectable expression are rendered, to prevent overplotting. c, Modules of genes differentially expressed over the myogenic trajectory. A total of 2,908 genes were clustered via hierarchical clustering. The dendrogram was cut into 14 modules using the cutree function in R, and the aggregate expression of genes in each module was computed. Colours indicate aggregate UMI counts for each module that have been scaled for library size, log-transformed and then mapped to Z-scores to enable comparison between modules. Cells with no expression of a given module are excluded to prevent overplotting.

Supplementary information

Supplementary Information

This file contains Supplementary Note 1, Supplementary References and full legends for Supplementary Tables 1-11

Reporting Summary

Supplementary Tables

This file contains Supplementary Tables 1-11

Rights and permissions

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Cao, J., Spielmann, M., Qiu, X. et al. The single-cell transcriptional landscape of mammalian organogenesis. Nature 566, 496–502 (2019). https://doi.org/10.1038/s41586-019-0969-x

Download citation

  • Received:

  • Accepted:

  • Published:

  • Issue Date:

  • DOI: https://doi.org/10.1038/s41586-019-0969-x

This article is cited by

Search

Advanced search

Quick links

Nature Briefing

Sign up for the Nature Briefing newsletter — what matters in science, free to your inbox daily.

Get the most important science stories of the day, free in your inbox. Sign up for Nature Briefing