Individual organs and tissues form in human embryos during the first two months of pregnancy. Any errors during this crucial stage of human development can result in miscarriage or serious birth defects. Yet remarkably little is known about how this process works. What is known has been inferred from studies of how other animals develop, human stem cells grown in a laboratory, and babies born with genetic conditions that cause developmental problems.

Genes control the way that organs and tissues form, and are switched on or off in complex patterns during development to ensure that particular cells develop into one type of organ and not another. When genes are switched on, their DNA is copied into molecules called RNA. Many RNA molecules are used as templates to make proteins, which then perform critical roles in cell processes. One way to find out which genes are activated during development is to identify which RNAs are made by cells in the embryo.

Here, Gerrard, Berry et al. used a technique called RNA-sequencing to identify the RNAs that human embryos make while their organs and tissues form. The RNA came from many different tissues including the heart, limbs and the roof of the mouth. Gerrard, Berry et al. developed a new computational model that used the identity of the RNAs to decode the precise patterns of gene activity in the tissues. The model correctly identified many genes that were already known to cause developmental problems when faulty, and identified numerous others that are now predicted to cause developmental defects in humans.

Gerrard, Berry et al. also discovered over 6,000 RNAs in the human embryos that are unlikely to code for proteins. These “non-coding” RNAs may have other roles in cells, such as switching off genes, and many of them appear to be specific to human embryos. Together, these findings have uncovered new patterns of gene activity that drive development in human embryos and provide a resource for studying how organs and tissues form. Future challenges are to understand what controls these patterns of gene activity, and how the patterns change over time.

Embryogenesis encompasses the progression from fertilized zygote to blastocyst and through gastrulation to establish the three germ layers of ectoderm, mesoderm and endoderm, from which all organs and tissues subsequently arise during organogenesis. Remarkably little is known about this latter phase of assembling organs and tissues in human due to the restricted availability of human embryonic tissue and its tiny size. Previous transcriptomics post-implantation have sampled either the whole embryo by expression microarray (Fang et al., 2010), thus lacking organ-specific resolution and the vast majority of long non-coding (lnc) transcription; or included lnc expression by massively parallel short-read RNA sequencing (RNA-seq) but focussed on single sites such as limb bud (Cotney et al., 2013) or pancreas (Cebola et al., 2015). RNA-seq from NIH Roadmap and other studies during or after the end of the first trimester of pregnancy falls after the embryonic period (which ends at 56–58 days post-conception (Carnegie Stage 23)) and commonly reflects near terminal differentiation within heterogeneous fetal organs and tissues (Jaffe et al., 2015; Roadmap Epigenomics Consortium, 2015; Roost et al., 2015). As a consequence of these combined deficiencies, we set about compiling global transcriptomic data during the critical phase of human organogenesis, sampling each germ layer and including sites of mixed origin that are subject to major developmental disorders such as cleft palate and limb abnormalities (Figure 1a).

Figure 1 with 5 supplements see all

Download asset Open asset

Organs and tissues from fifteen human embryonic sites were sequenced in two sets of biological replicates (except pancreas and tongue) to generate 28 strand-specific RNA-seq datasets with 44–90 million uniquely mapped reads per replicate (Figure 1a; Supplementary file 1A, which contains information on embryonic stages). Global transcription rates across all organs and tissues were comparable over a high dynamic range; approximately 70% of protein-coding genes contained 100–10,000 mapped reads (Figure 1b; Supplementary file 2). We assessed whether our human embryo datasets identified earlier developmental processes than currently available fetal data (Roadmap Epigenomics Consortium, 2015). There were three-fold the number of differentially expressed genes in the fetal datasets but equivalent enrichment of gene ontology (GO) terms in the embryo, including many early developmental processes such as morphogenesis of an epithelial bud, anterior/posterior pattern specification and embryonic morphogenesis. These were in contrast to homeostatic processes enriched in the fetal dataset (Figure 1c; Supplementary file 1B–C).

Sampling gene expression across multiple sites allowed us to set about deciphering the precise transcriptomic codes responsible for the development of the different human embryonic organs and tissues. While ZNF and ZSCAN family members were broadly expressed discrete site-specific expression was more apparent for individual members of other transcription factor families (Figure 1—figure supplement 1) exemplified by the HOX gene clusters (Figure 1d). User-defined sets of up to five developmental transcription factors characteristic for a particular organ or tissue displayed very high levels of tissue specificity (Figure 1—figure supplement 2). However, while principal components (PC) analysis (PCA) or clustering grouped biological replicates, relationships between different organs and tissues other than the distinctiveness of brain and liver were not resolved (Figure 1—figure supplements 3–4). Non-negative matrix factorisation (NMF) also allows unbiased clustering of gene expression (Gaujoux and Seoighe, 2010). By setting the parameters such that representative genes were only extracted once, we identified eleven non-overlapping ‘metagenes’ from the complete expression dataset with clear tissue-specific signals for thyroid, liver, RPE, brain, heart and adrenal gland (Figure 1—figure supplement 5; Supplementary file 1D). We hypothesized that these new signals might allow benchmarking to assess the fidelity of in vitro differentiated stem cells, similar to a previous report (Roost et al., 2015). As an exemplar, we chose hepatocyte differentiation for which RNA-seq data are available including positive (primary adult hepatocytes) and negative (human embryonic fibroblasts) control data (Du et al., 2014). Clear enrichment for the stem cell-derived hepatocytes and the primary hepatocytes (but not the fibroblasts) was apparent in metagene 2, the cluster of 39 genes indicative of human embryonic liver. From this starting point, we wanted to move beyond the unique organ-specific signatures to study how patterns of gene expression co-varied across tissues. While relaxing NMF parameters would allow non-exclusive gene selection across metagenes, we also wanted to capture differences in gene expression between organs (e.g. aspects of what is not expressed as a contributing factor to an organ’s identity). Moreover, different embryonic organs are related according to developmental lineage. We reasoned that being able to apply a lineage structure would create natural assemblies of co-regulated genes (Figure 2a). Accordingly, we adapted a method from spatial ecology and phylogenetics (Jombart et al., 2010a, 2010b) to constrain PCA by imposing a hierarchical developmental lineage and termed this approach LgPCA. We also included RNA-seq from undifferentiated human embryonic stem cells (Roadmap Epigenomics Consortium, 2015) to represent pre-gastrulation human biology. Of the total thirty-one principal components (PCs) arising from LgPCA the first fifteen now identified patterns of gene expression across groups of related tissues in addition to unique organ-specific signatures (Figure 2b) while PCs 16–31 sampled heterogeneity within individual organs and tissues (Figure 2—figure supplement 1). In keeping with this transition PCs 1–4 ordered samples reminiscent of very early developmental events: pluripotency (extreme positive loadings in PC1; ‘PC1 high’), early brain formation (extreme negative loadings in PC2; ‘PC2 low’), foregut endoderm (PC4 low) and intermediate mesoderm (PC4 high). PCs 5–15 resolved the individual organs and tissues; for instance low PC5 loadings discriminated liver from the other foregut endoderm derivatives. Heatmaps illustrated the composite or tissue-specific signals emanating from the genes with the most extreme PC loadings which also underlay appropriate developmental gene ontology (GO) terms (Figure 2c–d and Supplementary file 1E–F).

Figure 2 with 1 supplement see all

Download asset Open asset

Identifying the master regulators that differentially orchestrate organogenesis across the body has not previously been possible directly in human embryos. We undertook this in two different ways, both based on studying the 1000 genes with the most extreme loadings in PCs 1–15 that identified gene co-expression patterns across tissues and within individual organs (Figure 2b; Supplementary file 1E). We interrogated these gene sets for regulatory networks based on the enrichment of transcription factor motifs (Janky et al., 2014). Numerous well known master regulators were recovered alongside previously unappreciated factors for either broad tissue groups (e.g. foregut endoderm derivatives) or individual organs (Figure 3a). As proof-of-principle, this also included proven regulators of human pluripotency, NANOG, OCT4 and MYC, at an extreme of PC1. Remarkably, in several instances approaching half of the 1000 genes with the most extreme PC loadings imputed co-regulation by a single transcription factor, such as HNF4A in the liver or SRF in the heart. Alongside NR5A1, the data predicted RUNX and BAD as novel regulators of human adrenal and gonadal development (Figure 3a). As a second approach to study gene regulation, we extracted the transcription factors (typically <5%) from amongst the 1000 most extreme genes in PCs 1–15 (Supplementary file 1G). We searched the Mouse Genome Informatics database (MGI) and in 309/594 instances there was a relevant mouse mutation phenotype supporting the notion that the transcription factors identified by LgPCA are key regulators of human organogenesis. At the lowest extreme of PC5 (liver) the twenty-two transcription factors contained all three of those required for reprogramming fibroblasts directly to hepatocytes (Huang et al., 2014). This suggests novel fate programming roles for transcription factors at the extremes of other PCs (including new potential regulators of pluripotency amongst the sixteen factors containing zinc fingers in PC1). In keeping with these regulatory roles, the extreme PC loadings in the LgPCA data also prioritized those transcription factors responsible for major congenital disorders (Supplementary file 1G). Because LgPCA is not limited to individual organs this included a novel ability to predict multisystem abnormalities such as the combined heart and limb defects of Holt-Oram syndrome (OMIM 142900, TBX5, PC13 low) or the palate and limb abnormalities associated with mutations in TP63 (OMIM 603543, PC3 high).

Figure 3

Download asset Open asset

Mutations in genes encoding transcription factors are over-represented causes of congenital disorders, most likely due to their critical function during organogenesis and inadequacy when haploinsufficient. The enrichment of transcription factors with specific disease-associations at the extremes of the LgPCA implicates the co-enriched genes as leading candidates for unanswered clinical syndromes. To test this model we identified some of the earliest chromosomal mapping or patient deletion data for the known disease-associated transcription factors from Supplementary file 1G. 53 disorders were suitable for assessment with an average critical region of 13.7 Mb each containing an average of 111 protein-coding genes (Supplementary file 1H). Strikingly, in 37 instances (73%) LgPCA uniquely selected the correct transcription factor and in 48 instances (91%) narrowed the field down to three or fewer transcription factors. When applied to 13 syndromes (mostly deletion disorders) where the causative gene remains unresolved clear predictions of causality emerge, for instance in cleft palate (DLX5, DLX6, LHX8 and FOXF2) or cerebellar disorders (ZIC1 and ZIC4) (Supplementary file 1H). Frequently, there is an appropriate mutant mouse phenotype such as CASZ1 in cardiac malformations, part of Chr1p36 deletion syndrome, or SOX10 in the 46,XX disorder of sex differentiation (DSD) linked to duplication on Chr22 (Figure 3b and Supplementary file 1H).

Non-coding transcription has emerged as a critical regulator of cell and developmental biology (Goff and Rinn, 2015). A dedicated programme operating during human organogenesis seemed likely as 81 out of the 1571 genes enriched in embryogenesis compared to the fetal datasets were annotated long intergenic non-coding (LINC) transcripts (Supplementary file 1B). To look beyond this we assembled strand-specific transcripts not recognized by current genome annotation [GENCODE 18 (Harrow et al., 2012)] and systematically named them individually according to recommended criteria (Mattick and Rinn, 2015). 6251 unique loci accounted for in excess of 9 Mb of novel polyadenylated transcription from the human genome (Figure 4a and Supplementary file 1I). The vast majority of transcripts fulfilled criteria as lnc RNAs by assessment of coding potential (CPAT score <0.2) (Figure 4b), length over 200 base pairs (bp) and an absence of reads spanning splice junctions to currently annotated genes (Mattick and Rinn, 2015). These lncRNAs were classified as either bidirectional, antisense or overlapping, or by exclusion intergenic, according to orientation and position in relation to the annotated genome (Mattick and Rinn, 2015). Transcripts were most commonly 500–1,500 bp but could extend to over 600 Kb (Figure 4c) and showed high tissue-specificity with the median Tau value (Yanai et al., 2005) of 0.86, much higher than for protein-coding genes (0.63) but consistent with previously annotated non-coding genes (0.89). We investigated the association between this novel human embryonic transcriptome and the annotated genome. Reduced physical distance to expressed annotated genes markedly increased the likelihood of novel transcript co-expression (Figure 4e), although the best correlations were by no means always with the closest gene (Figure 4f–g). The median distance to the closest annotated gene was 7.7 Kb (Figure 4—figure supplement 1) while on average the best correlation was at 188 Kb (random prediction was 476 Kb). Over half (3634) of the lnc transcripts were classified positionally as LINC RNAs. While LINC RNAs can harbour important regulatory function, how to forecast their relationship(s) with the protein-coding genome and prioritize the investigation of thousands of new transcripts is immensely challenging (Goff and Rinn, 2015). As a first step, the multi-tissue nature of our dataset allowed intricate correlative patterns to be deciphered implying putative relationships; for instance over a 2 Mb window and across numerous genes on chromosome 7 between HE-LINC-C7T121 and TBX20, which encodes a developmental cardiac transcription factor mutated in a wide range of congenital heart disease (Figure 4h).

Figure 4 with 1 supplement see all

Download asset Open asset

Taken together, this study reports the first comprehensive transcriptomic atlas during human organogenesis to complement parallel initiatives from later development and adulthood (Jaffe et al., 2015; Roadmap Epigenomics Consortium, 2015; Roost et al., 2015). Subjecting transcription from many sites to a method of analysis that incorporated developmental lineage deciphered novel genetic signatures, predicted causality in many human developmental disorders and associated novel non-coding transcription with expression from the surrounding protein-coding genome. At present, the data arise from a relatively narrow window of embryonic development but set the stage for future longitudinal studies for individual organs over time. The tiny amounts and scarcity of human embryonic tissue also necessitated aspects of pooling across different Carnegie stages for some sites but it is striking that this had no impact on ascertaining organ and tissue-specific transcriptomic signatures by LgPCA. The integrated data are expected to be particularly valuable to stem cell researchers examining the fidelity of PSC differentiation in vitro or searching for transcription factors for direct reprogramming of chosen cell lineages. Finally, the discovery of a major new programme of non-coding transcription adds a fresh layer of detail on the spatiotemporal regulation of the human genome.

1. 1. Dave T Gerrard
   2. Andrew A Berry
   3. Rachel E Jennings
   4. Karen Piper Hanley
   5. Nicoletta Bobola
   6. Neil A Hanley

   (2016) An integrative transcriptomic atlas of organogenesis in human embryos

   Publicly available at EBI ArrayExpress (accession no: E-MTAB-3928).

   http://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-3928/

1. 1. Santiago A Rodriguez-Segui
   2. Andrew A Berry
   3. Rachel E Jennings
   4. Neil A Hanley
   5. Jorge Ferrer

   (2015) TEAD and YAP regulate the enhancer network of human embryonic pancreatic progenitors

   Publicly available at EBI ArrayExpress (accession no: E-MTAB-3061).

   http://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-3061/

1. Software

   1. Alexa A
   2. Rahnenfuhrer J

   (2010) Bioconductor, version topGO: topGO: Enrichment Analysis for Gene Ontology

   Bioconductor.

   http://bioconductor.org/packages/topGO/
2. Software

   1. Bolstad B

   (2007) Bioconductor, version preprocessCore: A Collection of Pre-Processing Functions

   Bioconductor.

   http://bioconductor.org/packages/preprocessCore/
3. 1. Brunet J-P
   2. Tamayo P
   3. Golub TR
   4. Mesirov JP

   (2004) Metagenes and molecular pattern discovery using matrix factorization

   PNAS 101:4164–4169.

   https://doi.org/10.1073/pnas.0308531101

   • Google Scholar
4. 1. Cebola I
   2. Rodríguez-Seguí SA
   3. Cho CH
   4. Bessa J
   5. Rovira M
   6. Luengo M
   7. Chhatriwala M
   8. Berry A
   9. Ponsa-Cobas J
   10. Maestro MA
   11. Jennings RE
   12. Pasquali L
   13. Morán I
   14. Castro N
   15. Hanley NA
   16. Gomez-Skarmeta JL
   17. Vallier L
   18. Ferrer J

   (2015) TEAD and YAP regulate the enhancer network of human embryonic pancreatic progenitors

   Nature Cell Biology 17:615–626.

   https://doi.org/10.1038/ncb3160

   • Google Scholar
5. 1. Cotney J
   2. Leng J
   3. Yin J
   4. Reilly SK
   5. DeMare LE
   6. Emera D
   7. Ayoub AE
   8. Rakic P
   9. Noonan JP

   (2013) The evolution of lineage-specific regulatory activities in the human embryonic limb

   Cell 154:185–196.

   https://doi.org/10.1016/j.cell.2013.05.056

   • Google Scholar
6. 1. Dillies MA
   2. Rau A
   3. Aubert J
   4. Hennequet-Antier C
   5. Jeanmougin M
   6. Servant N
   7. Keime C
   8. Marot G
   9. Castel D
   10. Estelle J
   11. Guernec G
   12. Jagla B
   13. Jouneau L
   14. Laloë D
   15. Le Gall C
   16. Schaëffer B
   17. Le Crom S
   18. Guedj M
   19. Jaffrézic F
   20. French StatOmique Consortium

   (2013) A comprehensive evaluation of normalization methods for Illumina high-throughput RNA sequencing data analysis

   Briefings in Bioinformatics 14:671–683.

   https://doi.org/10.1093/bib/bbs046

   • Google Scholar
7. 1. Du Y
   2. Wang J
   3. Jia J
   4. Song N
   5. Xiang C
   6. Xu J
   7. Hou Z
   8. Su X
   9. Liu B
   10. Jiang T
   11. Zhao D
   12. Sun Y
   13. Shu J
   14. Guo Q
   15. Yin M
   16. Sun D
   17. Lu S
   18. Shi Y
   19. Deng H

   (2014) Human hepatocytes with drug metabolic function induced from fibroblasts by lineage reprogramming

   Cell Stem Cell 14:394–403.

   https://doi.org/10.1016/j.stem.2014.01.008

   • Google Scholar
8. 1. Fang H
   2. Yang Y
   3. Li C
   4. Fu S
   5. Yang Z
   6. Jin G
   7. Wang K
   8. Zhang J
   9. Jin Y

   (2010) Transcriptome analysis of early organogenesis in human embryos

   Developmental Cell 19:174–184.

   https://doi.org/10.1016/j.devcel.2010.06.014

   • Google Scholar
9. 1. Gaujoux R
   2. Seoighe C

   (2010) A flexible R package for nonnegative matrix factorization

   BMC Bioinformatics 11:367.

   https://doi.org/10.1186/1471-2105-11-367

   • Google Scholar
10. 1. Goff L
    2. Rinn J

    (2015) Linking RNA Biology to lncRNAs

    Genome Research 25:1456–1465.

    https://doi.org/10.1101/gr.191122.115

    • Google Scholar
11. 1. Harrow J
    2. Frankish A
    3. Gonzalez JM
    4. Tapanari E
    5. Diekhans M
    6. Kokocinski F
    7. Aken BL
    8. Barrell D
    9. Zadissa A
    10. Searle S
    11. Barnes I
    12. Bignell A
    13. Boychenko V
    14. Hunt T
    15. Kay M
    16. Mukherjee G
    17. Rajan J
    18. Despacio-Reyes G
    19. Saunders G
    20. Steward C
    21. Harte R
    22. Lin M
    23. Howald C
    24. Tanzer A
    25. Derrien T
    26. Chrast J
    27. Walters N
    28. Balasubramanian S
    29. Pei B
    30. Tress M
    31. Rodriguez JM
    32. Ezkurdia I
    33. van Baren J
    34. Brent M
    35. Haussler D
    36. Kellis M
    37. Valencia A
    38. Reymond A
    39. Gerstein M
    40. Guigó R
    41. Hubbard TJ

    (2012) GENCODE: the reference human genome annotation for The ENCODE project

    Genome Research 22:1760–1774.

    https://doi.org/10.1101/gr.135350.111

    • Google Scholar
12. 1. Huang P
    2. Zhang L
    3. Gao Y
    4. He Z
    5. Yao D
    6. Wu Z
    7. Cen J
    8. Chen X
    9. Liu C
    10. Hu Y
    11. Lai D
    12. Hu Z
    13. Chen L
    14. Zhang Y
    15. Cheng X
    16. Ma X
    17. Pan G
    18. Wang X
    19. Hui L

    (2014) Direct reprogramming of human fibroblasts to functional and expandable hepatocytes

    Cell Stem Cell 14:370–384.

    https://doi.org/10.1016/j.stem.2014.01.003

    • Google Scholar
13. 1. Jaffe AE
    2. Shin J
    3. Collado-Torres L
    4. Leek JT
    5. Tao R
    6. Li C
    7. Gao Y
    8. Jia Y
    9. Maher BJ
    10. Hyde TM
    11. Kleinman JE
    12. Weinberger DR

    (2015) Developmental regulation of human cortex transcription and its clinical relevance at single base resolution

    Nature Neuroscience 18:154–161.

    https://doi.org/10.1038/nn.3898

    • Google Scholar
14. 1. Janky R
    2. Verfaillie A
    3. Imrichová H
    4. Van de Sande B
    5. Standaert L
    6. Christiaens V
    7. Hulselmans G
    8. Herten K
    9. Naval Sanchez M
    10. Potier D
    11. Svetlichnyy D
    12. Kalender Atak Z
    13. Fiers M
    14. Marine JC
    15. Aerts S
    16. Janky Rekin’ s

    (2014) iRegulon: from a gene list to a gene regulatory network using large motif and track collections

    PLoS Computational Biology 10:e1003731.

    https://doi.org/10.1371/journal.pcbi.1003731

    • Google Scholar
15. 1. Jennings RE
    2. Berry AA
    3. Kirkwood-Wilson R
    4. Roberts NA
    5. Hearn T
    6. Salisbury RJ
    7. Blaylock J
    8. Piper Hanley K
    9. Hanley NA
    10. Hanley N

    (2013) Development of the human pancreas from foregut to endocrine commitment

    Diabetes 62:3514–3522.

    https://doi.org/10.2337/db12-1479

    • Google Scholar
16. 1. Jennings RE
    2. Berry AA
    3. Strutt JP
    4. Gerrard DT
    5. Hanley NA

    (2015) Human pancreas development

    Development 142:3126–3163.

    https://doi.org/10.1242/dev.120063

    • Google Scholar
17. 1. Jombart T
    2. Balloux F
    3. Dray S

    (2010a) Adephylo: new tools for investigating the phylogenetic signal in biological traits

    Bioinformatics 26:1907–1909.

    https://doi.org/10.1093/bioinformatics/btq292

    • Google Scholar
18. 1. Jombart T
    2. Pavoine S
    3. Devillard S
    4. Pontier D

    (2010b) Putting phylogeny into the analysis of biological traits: a methodological approach

    Journal of Theoretical Biology 264:693–701.

    https://doi.org/10.1016/j.jtbi.2010.03.038

    • Google Scholar
19. 1. Mattick JS
    2. Rinn JL

    (2015) Discovery and annotation of long noncoding RNAs

    Nature Structural & Molecular Biology 22:5–7.

    https://doi.org/10.1038/nsmb.2942

    • Google Scholar
20. 1. Quinlan AR
    2. Hall IM

    (2010) BEDTools: a flexible suite of utilities for comparing genomic features

    Bioinformatics 26:841–842.

    https://doi.org/10.1093/bioinformatics/btq033

    • Google Scholar
21. 1. Roadmap Epigenomics Consortium
    2. Kundaje A
    3. Meuleman W
    4. Ernst J
    5. Bilenky M
    6. Yen A
    7. Heravi-Moussavi A
    8. Kheradpour P
    9. Zhang Z
    10. Wang J
    11. Ziller MJ
    12. Amin V
    13. Whitaker JW
    14. Schultz MD
    15. Ward LD
    16. Sarkar A
    17. Quon G
    18. Sandstrom RS
    19. Eaton ML
    20. Wu YC
    21. Pfenning AR
    22. Wang X
    23. Claussnitzer M
    24. Liu Y
    25. Coarfa C
    26. Harris RA
    27. Shoresh N
    28. Epstein CB
    29. Gjoneska E
    30. Leung D
    31. Xie W
    32. Hawkins RD
    33. Lister R
    34. Hong C
    35. Gascard P
    36. Mungall AJ
    37. Moore R
    38. Chuah E
    39. Tam A
    40. Canfield TK
    41. Hansen RS
    42. Kaul R
    43. Sabo PJ
    44. Bansal MS
    45. Carles A
    46. Dixon JR
    47. Farh KH
    48. Feizi S
    49. Karlic R
    50. Kim AR
    51. Kulkarni A
    52. Li D
    53. Lowdon R
    54. Elliott G
    55. Mercer TR
    56. Neph SJ
    57. Onuchic V
    58. Polak P
    59. Rajagopal N
    60. Ray P
    61. Sallari RC
    62. Siebenthall KT
    63. Sinnott-Armstrong NA
    64. Stevens M
    65. Thurman RE
    66. Wu J
    67. Zhang B
    68. Zhou X
    69. Beaudet AE
    70. Boyer LA
    71. De Jager PL
    72. Farnham PJ
    73. Fisher SJ
    74. Haussler D
    75. Jones SJ
    76. Li W
    77. Marra MA
    78. McManus MT
    79. Sunyaev S
    80. Thomson JA
    81. Tlsty TD
    82. Tsai LH
    83. Wang W
    84. Waterland RA
    85. Zhang MQ
    86. Chadwick LH
    87. Bernstein BE
    88. Costello JF
    89. Ecker JR
    90. Hirst M
    91. Meissner A
    92. Milosavljevic A
    93. Ren B
    94. Stamatoyannopoulos JA
    95. Wang T
    96. Kellis M

    (2015) Integrative analysis of 111 reference human epigenomes

    Nature 518:317–330.

    https://doi.org/10.1038/nature14248

    • Google Scholar
22. 1. Robinson MD
    2. McCarthy DJ
    3. Smyth GK

    (2010) edgeR: a Bioconductor package for differential expression analysis of digital gene expression data

    Bioinformatics 26:139–140.

    https://doi.org/10.1093/bioinformatics/btp616

    • Google Scholar
23. 1. Roost MS
    2. van Iperen L
    3. Ariyurek Y
    4. Buermans HP
    5. Arindrarto W
    6. Devalla HD
    7. Passier R
    8. Mummery CL
    9. Carlotti F
    10. de Koning EJ
    11. van Zwet EW
    12. Goeman JJ
    13. Chuva de Sousa Lopes SM
    14. Iperen Liesbeth van

    (2015) KeyGenes, a Tool to Probe Tissue Differentiation Using a Human Fetal Transcriptional Atlas

    Stem Cell Reports 4:1112–1124.

    https://doi.org/10.1016/j.stemcr.2015.05.002

    • Google Scholar
24. 1. Shannon P
    2. Markiel A
    3. Ozier O
    4. Baliga NS
    5. Wang JT
    6. Ramage D
    7. Amin N
    8. Schwikowski B
    9. Ideker T

    (2003) Cytoscape: a software environment for integrated models of biomolecular interaction networks

    Genome Research 13:2498–2504.

    https://doi.org/10.1101/gr.1239303

    • Google Scholar
25. 1. Trapnell C
    2. Roberts A
    3. Goff L
    4. Pertea G
    5. Kim D
    6. Kelley DR
    7. Pimentel H
    8. Salzberg SL
    9. Rinn JL
    10. Pachter L

    (2012) Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks

    Nature Protocols 7:562–578.

    https://doi.org/10.1038/nprot.2012.016

    • Google Scholar
26. 1. Trapnell C
    2. Williams BA
    3. Pertea G
    4. Mortazavi A
    5. Kwan G
    6. van Baren MJ
    7. Salzberg SL
    8. Wold BJ
    9. Pachter L

    (2010) Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation

    Nature Biotechnology 28:511–515.

    https://doi.org/10.1038/nbt.1621

    • Google Scholar
27. 1. Wang L
    2. Park HJ
    3. Dasari S
    4. Wang S
    5. Kocher JP
    6. Li W

    (2013) CPAT: Coding-Potential Assessment Tool using an alignment-free logistic regression model

    Nucleic Acids Research 41:e74.

    https://doi.org/10.1093/nar/gkt006

    • Google Scholar
28. 1. Yanai I
    2. Benjamin H
    3. Shmoish M
    4. Chalifa-Caspi V
    5. Shklar M
    6. Ophir R
    7. Bar-Even A
    8. Horn-Saban S
    9. Safran M
    10. Domany E
    11. Lancet D
    12. Shmueli O

    (2005) Genome-wide midrange transcription profiles reveal expression level relationships in human tissue specification

    Bioinformatics 21:650–659.

    https://doi.org/10.1093/bioinformatics/bti042

    • Google Scholar

Thank you for submitting your article "An integrative transcriptomic atlas of organogenesis in human embryos" for consideration by eLife. Your article has been reviewed by three peer reviewers, and the evaluation has been overseen by Janet Rossant as the Senior Editor and Reviewing Editor. One of the three reviewers has agreed to reveal their identity: Majlinda Lako.

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission. Given the nature of this contribution, we would like to consider it in the category of a Tools and Resources paper.

This paper provides transcriptional profiles from different tissues during early human embryogenesis and as such represents a significant and novel contribution. To date most of our knowledge on human embryogenesis is based on extrapolation from model animals. These new data will provide a resource to the community enabling more understanding of tissue differentiation and organogenesis. It will also help with identification of new disease genes and understanding developmental defects.

Essential revisions:

1) There are some important issues raised by two of the reviewers with regard to the statistical analysis of the data. These concerns are provided in detail below and must be addressed in the revised manuscript.

2) The linkage of the transcriptional data to developmental disorders is not well validated. There is considerable mouse developmental and phenotype data available in databases that could be mined to enhance the value of your findings and help determine whether the genes identified are truly key regulators of tissue development.

3) The claim that novel transcripts are likely lncRNAs was questioned and needs to be analyzed more carefully.

Statistical issues:

There are concerns of the methodology that would require further clarification:

For the NMF, while the implementation of the methods is likely to be appropriate, there may be a potential flaw in the interpretation of what the analysis may achieved. Regarding "non-overlapping" metagenes, by default, genes can be represented in multiple metagenes, based on the Brunet algorithm in the NMF R package (also the default algorithm). The Methods state that this algorithm has been used, so following on from that, the assertion that "non-overlapping metagenes" could be extracted from the complete dataset and the "potential of coordinated deployment of overlapping genes" is ignored may be incorrect. Furthermore, the critique that "NMF failed to discriminate transcriptional signatures for a number or organs or tissues, or discern the relationships between them… " may also apply to the outcome of lgPCA. Therefore, this is not a compelling argument for discarding the NMF.

The analysis could have been taken further to examine the enrichment of all the metagenes, similarly to the one for liver (Figure 1—figure supplement 3B) for all the specific tissue metagenes (“clear tissue-specific signals for thyroid, liver, RPE, brain, heart and adrenal gland”), and a similar approach could be taken for the downstream functional analysis of lgPCA.

For lgPCA, while the interpretation of the results would be correct, there is confusion in the nomenclature of loading and scores, which are fundamental concepts in PCA. PC scores refer to the eigenvalues (magnitude of variance) which determine the separation between the samples (as can be gleaned from Figure 2B). PC loadings are eigenvectors and describe the contribution of the variables (in this case genes) that causes the separation of the samples. Therefore, PC scores should rightly be PCA loadings (main text, third and fourth paragraphs). This confusion should be attended to. Additionally, it would be helpful to mention at least once in the main text that when referring to "PC1 low" for example, this is meant to be low/negative PC1 scores.

Both NMF and lgPCA are doing similar things, i.e. clustering samples based on gene expression to find which genes drive sample separation. The differences are that NMF can reveal the similarity between gene expression patterns for samples, which are clustered with the same metagenes, in this case kidney and testis to metagene 3, (Figure 2—figure supplement 3B). LgPCA can discern the differences in gene expression between certain samples (for example Brain and liver in PC2 of Figure 2B). If it was intended to include the NMF analysis in this study, a point should be made as the rationale of which samples are to be clustered together and why, before the results may be compared to those from the lgPCA. Alternatively, the NMF analysis can be omitted since the lgPCA method is sufficient to accomplish the data analysis.

The analysis presented in Figure 1C compares the authors' data with that presented in the fetal datasets from the NIH Roadmap. The authors claim substantial up-regulation of sets of genes in their dataset but this is not quantified statistically (instead an arbitrary cut-off of >2-fold enrichment was employed). While there is relatively little data from which to perform a rigorous analysis it should be possible to use appropriate models (gaining power across genes such as GSEA) to find enriched categories in a more robust manner than presented here (i.e., this would allow appropriate correction for multiple testing).

Essential revisions:

1) There are some important issues raised by two of the reviewers with regard to the statistical analysis of the data. These concerns are provided in detail below and must be addressed in the revised manuscript.

Thank you. We have addressed these issues in turn in the relevant later section.

2) The linkage of the transcriptional data to developmental disorders is not well validated. There is considerable mouse developmental and phenotype data available in databases that could be mined to enhance the value of your findings and help determine whether the genes identified are truly key regulators of tissue development.

We agree that this is important. The signposting in the original submission was weak. We have addressed this by re-writing the main text to reference clearly Supplementary file 1G (was F), which details 594 transcription factors at the extremes of principal components in the LgPCA. For each of these genes we mined the OMIM database for human developmental disorders (Column G of Supplementary file 1G). To link to mouse developmental biology we mined the Mouse Genome Informatics (MGI) database for each of the 594 genes providing each MGI ID in Column I and consistent mouse phenotypic features (identifiable in 309 instances) in Column J. In addition, each gene was mined in Pubmed; where additional information from mouse models was particularly insightful over and above the MGI data PMIDs are given in Column K. These data reinforce the message that our identified transcription factors are indeed developmental regulators. In addition, it is relevant that the LgPCA led to consistent GO terms (e.g. Figure 2D, ‘metanephros development’ for the embryonic kidney signature). These GO terms and the genes underlying them relate to regulatory roles discovered during mouse development.

3) The claim that novel transcripts are likely lncRNAs was questioned and needs to be analyzed more carefully.

We have debated this fast-moving area amongst ourselves throughout the course of the study. We considered it best to draw on up-to-date criteria and used the commentary by Mattick & Rinn in Nature Structural & Molecular Biology last year.

An important point is scale: we have discovered a huge amount of unannotated transcription and catalogued it as systematically as possible. This approach and publicising its existence is important because its scope (>6,000 transcripts) is beyond the capacity for functional investigation by any one lab. Thus, we view our job as to make the data open access for others to interrogate. This point was made by Mattick & Rinn: ‘Although these transcripts may be collectively referred to as long ‘noncoding’ RNAs (i.e. lncRNAs), this is a necessarily vague catch all term that is useful only as a moniker until they are better characterized’. We recognise, agree with and have adopted this initial starting point of using the ‘catch all term’, lncRNA. In addition, we have gone to great length to build on their suggested nomenclature (Figure 4A and Supplementary file 1I) rather than pursue the easier option of simple ‘in-house’ numbering (prevalent in other manuscripts) which lacks reference to genomic structure and location and is less helpful to the community (‘Given the enormous numbers of long apparently noncoding transcripts, it is important to have a logical and flexible structure for cataloguing and parsing newly found transcripts, in a way that will foster rather than impede their functional and mechanistic exploration.’).

With this in mind: ‘The broad term lncRNA refers to a transcript >200 nt in length that does not appear to contain a protein-coding sequence.’ The 200bp threshold is regarded as a practical lower limit for standard RNA-seq thus excluding smaller RNAs, such as tRNAs, snRNAs and MIRs. We assessed coding potential using CPAT (Figure 4B; Wang et al., NAR, cited in the manuscript), a method recommended by Mattick & Rinn. Our cut-off was <0.2, stricter than the 0.364 limit (determined by two-graph user operator characteristic) from the CPAT User Guide to increase our confidence that these transcripts are truly lncRNAs. Of note, over 90% of our novel transcripts had this very low coding potential.

We recognize that histone modifications can be useful for increasing the confidence with which transcriptional start sites of lncRNAs can be ascribed. We have started these studies (e.g. H3K4me3 ChIP-seq) for a number of organs as a part of a separate analysis of genome regulation (e.g. organ-specific enhancer usage during human embryogenesis). While too preliminary for inclusion here and well beyond the 2-month requested turnaround for this revision, it is reassuring that our early data on H3K4me3 marks do indeed correlate with the start of our lncRNAs.

Taken together, we feel that ‘lncRNA’ is the best descriptor for our novel transcripts filtered for length >200 bp and CPAT scores <0.2.

Statistical issues:

We note that the end-point of the Referees’ discussion on the NMF data is that we could simply remove it (which is a very helpful option—thank you!). We would like to reserve that opportunity; however, for now, we would like to try and preserve the NMF subject to the clarifications below for the following reasons. In previous advice we were encouraged to include an alternative approach that specifically targets organ-specific signatures. This seems reasonable. In addition, we feel that our manuscript will interest a broad audience not immediately familiar with the wide range of possible computational analyses. Therefore, it seems fair and reasonable to publicise the potential of NMF alongside our own approach of LgPCA. Clearly, if our clarifications and revisions are not convincing, the NMF analysis can be removed.

There are concerns of the methodology that would require further clarification:

For the NMF, while the implementation of the methods is likely to be appropriate, there may be a potential flaw in the interpretation of what the analysis may achieved. Regarding "non-overlapping" metagenes, by default, genes can be represented in multiple metagenes, based on the Brunet algorithm in the NMF R package (also the default algorithm). The Methods state that this algorithm has been used, so following on from that, the assertion that "non-overlapping metagenes" could be extracted from the complete dataset and the "potential of coordinated deployment of overlapping genes" is ignored may be incorrect.

Thank you for this point. Our original submission was not clear, some of the text needed moderation and we failed to reference Supplementary file 1D. We applied the algorithm specifically to decipher tissue-defining gene-sets. Therefore, we used only those genes with ‘basis contribution’ >80% across metagenes meaning that each of the genes featured in only one metagene. This is depicted in Figure 1—figure supplement 5 and now referenced properly as Supplementary file 1D. In not describing this adequately, we gave the false perception that this was the limit of the NMF approach when, as the reviewers correctly highlight, we could have reduced the threshold to allow overlapping sets of genes. We have amended to text to discuss this point and clarify our rationale.

Furthermore, the critique that "NMF failed to discriminate transcriptional signatures for a number or organs or tissues, or discern the relationships between them… " may also apply to the outcome of lgPCA. Therefore, this is not a compelling argument for discarding the NMF.

The original text was poorly chosen and we have tempered this section and the final paragraph. As a general aside, we feel that signatures were more comprehensively discerned in the LgPCA, in part by the opportunity to study profiles across PCs (e.g. adrenal gland across PC8 and PC9 as discussed in the text and as in Supplementary file 1F–G). Moreover, we wanted to explore a methodology that permitted a lineage to be imparted upon the global dataset to model human embryogenesis. We have tried to get these points across in revised text.

The analysis could have been taken further to examine the enrichment of all the metagenes, similarly to the one for liver (Figure 1—figure supplement 3B) for all the specific tissue metagenes (“clear tissue-specific signals for thyroid, liver, RPE, brain, heart and adrenal gland”), and a similar approach could be taken for the downstream functional analysis of lgPCA.

The inclusion of this figure panel and these data was in response to a helpful suggestion from a stem cell expert who commented on the importance of our data in allowing researchers to benchmark in vitro differentiated PSCs. In the example shown, we have taken the NMF metagenes (Figure 1—figure supplement 5A) and tested for gene enrichment in data from Du et al. (2014, Cell Stem Cell) for human PSCs differentiated to hepatocytes with primary human adult hepatocytes (positive control) and human embryonic fibroblasts (negative control). The hepatocytes (both native and PSC-derived) from Du et al. showed appropriate signals in Metagene 2 (liver) but not the other metagenes. (It is interesting to note some signal for metagene 1 (thyroid – another anterior endoderm derivative) in the in vitro PSC differentiated hepatocyte-like cells.) This analysis required robust differentiation protocols and suitable RNA-seq data from differentiated PSCs; hence why we chose the exemplar of liver were data are available from Du and colleagues and hence why we feel it would be premature to extend the analyses to other NMF metagene or LgPCA organ signatures (lack of in vitro protocols/RNA-seq data).

For lgPCA, while the interpretation of the results would be correct, there is confusion in the nomenclature of loading and scores, which are fundamental concepts in PCA. PC scores refer to the eigenvalues (magnitude of variance) which determine the separation between the samples (as can be gleaned from Figure 2B). PC loadings are eigenvectors and describe the contribution of the variables (in this case genes) that causes the separation of the samples. Therefore, PC scores should rightly be PCA loadings (main text, third and fourth paragraphs). This confusion should be attended to.

Thank you. We are very grateful for these important points. We have corrected this oversight.

Additionally, it would be helpful to mention at least once in the main text that when referring to "PC1 low" for example, this is meant to be low/negative PC1 scores.

We have added additional text.

Both NMF and lgPCA are doing similar things, i.e. clustering samples based on gene expression to find which genes drive sample separation. The differences are that NMF can reveal the similarity between gene expression patterns for samples, which are clustered with the same metagenes, in this case kidney and testis to metagene 3, (Figure 2—figure supplement 3B). LgPCA can discern the differences in gene expression between certain samples (for example Brain and liver in PC2 of Figure 2B). If it was intended to include the NMF analysis in this study, a point should be made as the rationale of which samples are to be clustered together and why, before the results may be compared to those from the lgPCA. Alternatively, the NMF analysis can be omitted since the lgPCA method is sufficient to accomplish the data analysis.

We have adapted the text. NMF can indeed generate metagenes across tissues (e.g. kidney and testis as pointed out) as can LgPCA (e.g. limbs and palate in PC3). As commented on, LgPCA can also discern differences in gene co-expression between organs (e.g. tongue and heart in PC15). To note, NMF does not start from a position of pre-determined sample clustering. Hoping to resolve this issue and recognising that the original text was unhelpful on NMF versus LgPCA, we have made revisions. As above, if these measures are inadequate, we can always remove the NMF.

The analysis presented in Figure 1C compares the authors' data with that presented in the fetal datasets from the NIH Roadmap. The authors claim substantial up-regulation of sets of genes in their dataset but this is not quantified statistically (instead an arbitrary cut-off of >2-fold enrichment was employed). While there is relatively little data from which to perform a rigorous analysis it should be possible to use appropriate models (gaining power across genes such as GSEA) to find enriched categories in a more robust manner than presented here (i.e., this would allow appropriate correction for multiple testing).

Thank you for this valid point. We are grateful for the Referees noting the relative paucity of data for this analysis. Nevertheless, we feel that it is a useful early component of the paper to show that there are indeed differences between human embryonic datasets and those from later development from NIH Roadmap. We have undertaken alternative analyses based on statistical cut-offs and with methodology that facilitates multiple testing. These data are now included in the manuscript in place of the previous Figure 1C.

Author details

1. Dave T Gerrard

   Division of Diabetes, Endocrinology & Gastroenterology, School of Medical Sciences, Faculty of Biology, Medicine and Health, Manchester Academic Health Science Centre, University of Manchester, Manchester, United Kingdom

   Contribution

   DTG, Wrote the manuscript, Conducted the bioinformatics analyses, Devised the study and planned experiments, Conception and design, Analysis and interpretation of data, Drafting or revising the article

   Contributed equally with

   Andrew A Berry

   Competing interests

   The authors declare that no competing interests exist.

   "This ORCID iD identifies the author of this article:" 0000-0001-6890-7213

2. Andrew A Berry

   Division of Diabetes, Endocrinology & Gastroenterology, School of Medical Sciences, Faculty of Biology, Medicine and Health, Manchester Academic Health Science Centre, University of Manchester, Manchester, United Kingdom

   Contribution

   AAB, Collected, dissected and prepared material, Devised the study and planned experiments, Conception and design, Acquisition of data, Drafting or revising the article

   Contributed equally with

   Dave T Gerrard

   Competing interests

   The authors declare that no competing interests exist.

3. Rachel E Jennings

   1. Division of Diabetes, Endocrinology & Gastroenterology, School of Medical Sciences, Faculty of Biology, Medicine and Health, Manchester Academic Health Science Centre, University of Manchester, Manchester, United Kingdom
   2. Endocrinology Department, Central Manchester University Hospitals NHS Foundation Trust, Manchester, United Kingdom

   Contribution

   REJ, Collected, dissected and prepared material, Devised the study and planned experiments, Conception and design, Acquisition of data, Drafting or revising the article

   Competing interests

   The authors declare that no competing interests exist.

4. Karen Piper Hanley

   Division of Diabetes, Endocrinology & Gastroenterology, School of Medical Sciences, Faculty of Biology, Medicine and Health, Manchester Academic Health Science Centre, University of Manchester, Manchester, United Kingdom

   Contribution

   KPH, Involved in study design, Conception and design, Drafting or revising the article

   Competing interests

   The authors declare that no competing interests exist.

5. Nicoletta Bobola

   Division of Dentistry, School of Medical Sciences, Faculty of Biology, Medicine and Health, Manchester Academic Health Science Centre, University of Manchester, Manchester, United Kingdom

   Contribution

   NB, Wrote the manuscript, Conducted the analysis of developmental disorders, Devised the study and planned experiments, Conception and design, Analysis and interpretation of data, Drafting or revising the article

   Competing interests

   The authors declare that no competing interests exist.

6. Neil A Hanley

   1. Division of Diabetes, Endocrinology & Gastroenterology, School of Medical Sciences, Faculty of Biology, Medicine and Health, Manchester Academic Health Science Centre, University of Manchester, Manchester, United Kingdom
   2. Endocrinology Department, Central Manchester University Hospitals NHS Foundation Trust, Manchester, United Kingdom

   Contribution

   NAH, Wrote the manuscript, Guarantor, Devised the study and planned experiments, Conception and design

   For correspondence

   neil.hanley@manchester.ac.uk

   Competing interests

   The authors declare that no competing interests exist.

   "This ORCID iD identifies the author of this article:" 0000-0003-3234-4038

• 8,021

  Page views

• 1,230

  Downloads

• 45

  Citations

Article citation count generated by polling the highest count across the following sources: Crossref, Scopus, PubMed Central.