An irregular hourglass pattern describes the tempo of phenotypic development in placental mammal evolution

Organismal development is defined by progressive transformations that ultimately give rise to distinct tissues and organs. Thus, temporal shifts in ontogeny often reflect key phenotypic differences in phylogeny. Classical theory predicts that interspecific morphological divergence originates towards the end of embryonic or fetal life stages, i.e. the early conservation model. By contrast, the hourglass model predicts interspecific variation early and late in prenatal ontogeny, though with a phylogenetically similar mid-developmental period. This phylotypic period, however, remains challenging to define within large clades such as mammals. Thus, molecular and morphological tests on a mammalian hourglass have not been entirely congruent. Here, we report an hourglass-like pattern for mammalian developmental evolution. By comparing published data on the timing of 74 homologous characters across 51 placental species, we demonstrated that variation in the timing of development decreased late in embryogenesis––when organ formation is highly active. Evolutionary rates of characters related to this timeframe were lowest, coinciding with a phylotypic period that persisted well beyond the pharyngula ‘stage’. The trajectory culminated with elevated variation in a handful of fetal and perinatal characters, yielding an irregular hourglass pattern. Our study invites further quantification of ontogeny across diverse amniotes and thus challenges current ideas on the universality of developmental patterns.

. The phylopars function also simultaneously estimated maximum likelihood ancestral state states for each character [17].
The rank-ordered ancestral character sequence was highly concordant with a character sequence ranked by means for relative timing in extant species (using phylogenetically and nonphylogenetically imputed data sets) (Fig. S8). Data sets imputed with missMDA and with Rphylopars were highly congruent (Fig. S8). Subsequent analyses on interspecific character variation were performed using the phylogenetically imputed data set and the rank-ordered ancestral character (developmental) sequence, though character variation was also plotted against relative timing for comparison (Fig. S9). Reconstructed ancestral states ultimately depend on the availability of trait data, as well as the taxa of choice [19]. To explore whether the chosen taxa and characters introduced biases in ancestral state reconstructions, a taxonomically and phenotypically restricted data set (N = 44 species; N = 66 characters), though with a lower missing data rate of 15%, was analyzed and compared to the analyses based on 51 species and 74 characters (21% missing data rate). Reconstructed ancestral means for the relative timing of characters and their sequence ranks were highly similar between restricted and full data sets, indicating that no major biases were introduced by performing ancestral state reconstructions on the phenotypically broader and more taxonomically inclusive data set of 51 species (Fig. S10).
To test for phylogenetic signal in our data set, we computed the multivariate version of Blomberg's K statistic (Kmult; [20]) using the geomorph R package [21]. As our data set exhibited phylogenetic signal (Kmult = 0.72; P < 0.0001, permutations = 10,000), we explored and compared interspecific differences in the relative timing of characters by performing a phylogenetic PCA in Phytools [22]. Graphical inspection of the phylogenetic PCA was performed via ggplot2 using the custom ggphylomorpho function (https://github.com/wabarr/ggphylomorpho). Principal components (PCs) 1 and 2 explained 80% of variation in the data set (Figs. S11-12; Table S1). PCs 1 and 2 were reconstructed for internal nodes and mapped onto the phylogeny of placental mammals using the contMap function of Phytools [23]. With the exception of Fig. S13, graphs were created using ggplot2, Phytools or base R functions.
Evolutionary rates (BM rate parameter: σ 2 ) for each character were estimated, using Phytools, by computing phylogenetic independent contrasts on log-transformed data (Table S2). To ascertain the extent to which character evolutionary rates differed among each other, we compared the likelihood of a model in which characters featured distinct evolutionary rates to the likelihood of a model in which all characters shared a common rate [15]. This test assumed a BM model of evolution and was performed using the likelihood-based methodology (including R code) developed by D.C. Adams [see 15]. Illustrated specimens exemplifying the mid-developmental period during which interspecific variation decreases and rates of evolution were lowest, see Table S2, are depicted in Fig.  S13. Rates of evolution for all characters and their 95% confidence intervals are depicted in Figs. S14-15.    Figure S8. Diagnostic plots displaying correspondence (Spearman's r = 0.99) between the rank-based character sequence in extant placental mammals versus the reconstructed ancestral sequence (A). Means for relative timing of characters were also strongly concordant with reconstructed ancestral values (B; r = 0.99). Character sequences based on non-phylogenetically imputed versus phylogenetically imputed data sets were tightly correlated (C;

Supplementary Figures and Tables
Spearman's r = 0.99). Coefficients of variation based on non-phylogenetically imputed versus phylogenetically imputed data sets were also congruent (D; r = 0.97), though the relationship was moderate for extreme values representing fetal and perinatal (including birth) characters (see Table S2 below). Coeff. variation (phylo. imputed) Coeff. variation (non-phylo. imputed) D Figure S9. The coefficient of variation is plotted against relative timing scores; based on a character data set that included phylogenetically imputed missing values for 51 placental mammal species. Locally estimated scatterplot smoothing was applied to the data to visualize the mean trend (blue dashed line).      Figure S13. Representative embryos displaying characters associated with the embryo-to-fetus transition in placental mammals (see Tables S1-2). Images are modified from [2][3][4][5].
Figure S14. Evolutionary rates for developmental characters in placental mammals. Characters are in ascending order according to evolutionary rates (range bars = 95% confidence intervals).