Comparative Transcriptomics Reveals Distinct Patterns of Gene Expression Conservation through Vertebrate Embryogenesis

Abstract

Despite life's diversity, studies of variation often remind us of our shared evolutionary past. Abundant genome sequencing and analyses of gene regulatory networks illustrate that genes and entire pathways are conserved, reused, and elaborated in the evolution of diversity. Predating these discoveries, 19th-century embryologists observed that though morphology at birth varies tremendously, certain stages of vertebrate embryogenesis appear remarkably similar across vertebrates. In the mid to late 20th century, anatomical variability of early and late-stage embryos and conservation of mid-stages embryos (the "phylotypic" stage) was named the hourglass model of diversification. This model has found mixed support in recent analyses comparing gene expression across species possibly owing to differences in species, embryonic stages, and gene sets compared. We compare 186 microarray and RNA-seq data sets covering embryogenesis in six vertebrate species. We use an unbiased clustering approach to group stages of embryogenesis by transcriptomic similarity and ask whether gene expression similarity of clustered embryonic stages deviates from a null expectation. We characterize expression conservation patterns of each gene at each evolutionary node after correcting for phylogenetic nonindependence. We find significant enrichment of genes exhibiting early conservation, hourglass, late conservation patterns in both microarray and RNA-seq data sets. Enrichment of genes showing patterned conservation through embryogenesis indicates diversification of embryogenesis may be temporally constrained. However, the circumstances under which each pattern emerges remain unknown and require both broad evolutionary sampling and systematic examination of embryogenesis across species.

Keywords: developmental hourglass; diversification; evo-devo; phylotypic stage.

Fig. 1

Anatomical and gene expression similarity predicted under different models of conservation through embryogenesis (A). Expression conservation was assessed using 186 publicly available microarray and RNA-seq data sets through embryogenesis across a phylogeny of six vertebrate species (B). Divergence times at each node are shown in millions of years.

Fig. 2

The number of studies testing the developmental hourglass and alternative hypotheses has increased considerably since the start of the new millennium (A, gray line: r² = 0.5, F_{(1, 19)} = 21, P = 0.0002), driven in large part by an increase in the number of gene expression studies (A, black line: r² = 0.27, F_{(1, 19)} = 8.3, P = 0.01). Despite this increased research effort, whether variation in embryogenesis follows an hourglass pattern has remained unresolved (B, orange points: r² = 0.008, F_{(1, 19)} = 1.2, P = 0.29). Neither divergence time of the species compared (C), nor number of species included in any given study (D) affect whether an hourglass pattern is observed. Quantitative literature analysis of studies examining early embryonic development across species was carried out according to the PRISMA flow diagram (supplementary fig. 1, Supplementary Material online; Moher et al. 2009). Detailed methods are provided in supplementary material (supplementary tables S1 and S2 and fig. S1, Supplementary Material online).

Fig. 3.

k-means clustering of the microarray (A) and RNA-seq (B) data sets analyzed. Reduction of within cluster variance increases as the number of cluster (k) increases. Gains asymptote at approximately k = 5 (dashed line).

Fig. 4

Spearman rank correlations were used to group stages into five clusters. Shown are all pairwise correlations of stages for all species and both gene expression profiling technologies. Grouping of stages is in indicated color.

Fig. 5

Overlap of embryonic stage cluster and major developmental events. Variation in sampling and heterochrony among species result in differences of developmental events captured by the available data across species and platforms. For three species, gene expression data sets were obtained from both microarray (A) and RNA-seq (B) platforms.

Fig. 6

Spearman rank correlations for pairwise comparisons of species at each cluster of embryogenesis for microarray (A) and RNA-seq data (B). Gene expression correlations (as a measure of conservation) vary through embryogenesis for both microarray and RNA-seq data and show support for a developmental hourglass. Colored boxes indicate observed correlations; gray boxes indicate rank correlations after permutation analysis randomizing stage association with cluster. Asterisks indicate that the observed median correlation differs significantly from the null expectation at P < 0.01. Note that only the RNA-seq data showed a temporal pattern (consistent with the early conservation hypothesis) that differed from the null expectation (i.e., no temporal pattern present).

Fig. 7

Enrichment of temporally patterned gene expression conservation across vertebrates. The expected number of genes for a given temporal expression pattern is determined by the proportion of trajectories of each conservation pattern (A.2: microarray; B.2: RNA-seq). The proportion of trajectories differed across conservation patterns and between platforms (A.2 and B.2). The enrichment/depletion values from permutation analysis (1,000 iterations) are shown as violin plots (A.1 and B.1). The number of genes associated with all conservation patterns were enriched compared with random expectation, except for the inverse hourglass (RNA-seq only), whereas the number of genes with no relationship between expression conservation across embryogenesis (A.1: microarray; B.1: RNA-seq) was significantly less than expected by chance. For all plots, colors indicate the expression conservation pattern. Enrichment at other evolutionary nodes (tetrapods, amniotes, and anurans) are provided in supplementary figure S4, Supplementary Material online.

Fig. 8

A total of 255 one-to-one orthologs shared between microarray and RNA-seq data sets largely differed in expression conservation pattern (A). Most genes with shared expression conservation patterns were “no relationship” genes followed by hourglass or inverse hourglass genes (B).