Speciation and the developmental alarm clock.

New species arise as the genomes of populations diverge. The developmental 'alarm clock' of speciation sounds off when sufficient divergence in genetic control of development leads hybrid individuals to infertility or inviability, the world awoken to the dawn of new species with intrinsic post-zygotic reproductive isolation. Some developmental stages will be more prone to hybrid dysfunction due to how molecular evolution interacts with the ontogenetic timing of gene expression. Considering the ontogeny of hybrid incompatibilities provides a profitable connection between 'evo-devo' and speciation genetics to better link macroevolutionary pattern, microevolutionary process, and molecular mechanisms. Here, we explore speciation alongside development, emphasizing their mutual dependence on genetic network features, fitness landscapes, and developmental system drift. We assess models for how ontogenetic timing of reproductive isolation can be predictable. Experiments and theory within this synthetic perspective can help identify new rules of speciation as well as rules in the molecular evolution of development.

Introduction

“Are certain developmental processes especially likely to be disrupted in hybrids? This question has been surprisingly neglected given that hybrid defects provide a rare window on those developmental processes and pathways that diverge rapidly between taxa.’ – Coyne and Orr, 2004, p.309.

Distinct species are those separate collections of genomes that, if you were to put them together in the same cells, the mixture would create a broken organism. Development in the hybrid individuals would go awry in such a way as to cause inviability, infertility, or a phenotypic mismatch to ecology or mating interactions that compromises further reproductive success. Given common descent from an ancestor, how could evolution produce such disastrous phenotypic consequences of genetic change? Darwin recognized this dilemma (Darwin 1859), in that natural selection will oppose changes that confer a net fitness cost. But evolution can circumvent this problem through interactions between multiple genetic factors, as intuited by W. Bateson, T. Dobzhansky and H.J. Muller and made explicit in the genetic mechanism of Dobzhansky-Muller incompatibility (DMI) for post-zygotic disruption in hybrids (Box 1A,B). Put simply: when evolution substitutes an allele at one locus, it makes no guarantee that this new derived genetic background will be compatible with allele substitutions occurring in other populations at other loci. Such inter-locus incompatibilities could involve two, or three, or many more interacting loci to create a DMI to genetically enforce species boundaries.

Visualizing, conceptualizing, and modeling Dobzhansky-Muller incompatibilities.

Dobzhansky-Muller Incompatibilities (DMIs) can be viewed as disrupted gene networks and as valleys (or holes) in fitness landscapes (Gavrilets, 2003; Fragata et al., 2019). Distinct biological species comprise separate groups of individuals that fail to successfully interbreed, though ‘good’ biological species may still yield hybrid F1 and F2 offspring, so long as they suffer clear fitness deficits (Coyne and Orr, 2004). Of special interest for development are genetically intrinsic post-zygotic species barriers that manifest as DMIs in hybrid individuals and do not depend on extrinsic circumstances, reflecting disruptive changes to developmental programs. DMIs between populations can evolve from a single common ancestral population through the independent substitution of two or more mutations distinguishing the descendant lineages. The substitutions may fix due to positive selection or genetic drift. This DM model also encapsulates the essence of developmental system drift (DSD): independent evolution in distinct lineages that causes divergence in genetic architectures while retaining within-lineage fitness (True and Haag, 2001).

Visualizations of genetic divergence and incompatibilities.

Species trees in (A) show the history of substitutions for loci in two genetic networks that differ in molecular evolution (case on left with e.g. fewer genes, stronger purifying selection, less adaptation, greater modularity, or lower pleiotropic effects) and in the potential for DMIs (solid purple lines = potential DMI for derived-derived substitutions between species, dashed lines = potential derived-ancestral DMIs; lowercase = ancestral alleles, uppercase = derived substitutions unique to one lineage). The number of potential DMIs (purple) scales faster than linear with number of substitutions (red and blue hashes); faster evolving genetic networks may be more likely to experience this ‘snowball effect’ of reproductive isolation (Orr, 1995). DSD arises when the outward phenotype remains constant despite molecular divergence between descendant species. Panel (B) shows how two loci (a and b) that diverge can potentially create a DMI upon formation of F1 hybrids between descendant species. Panel (C) illustrates with a Fisher’s geometric model visualization, for two traits with a shared genetic architecture, how adaptive evolution with respect to one trait (Trait 1) can generate DSD in another (Trait 2). Concentric circles represent lines of equal fitness; filled dots (black = ancestor, red and blue = descendant species) indicate genotypes (letters as in A) that evolve via three substitutions (arrows) toward the fitness optimum at the center. Note the DSD in Trait 2 due to no net phenotypic difference relative to the ancestor at the end of the adaptive walk for both species, despite underlying genetic changes. Panel (D) shows evolution along ridges of equal fitness in a fitness landscape comprised of a genetic architecture with many genes. Genotypic paths evolve independently in different species (ancestral black to derived red and blue species), similarly to DSD, such that hybrids between them (purple) occupy a portion of genotype space with low fitness (‘holes’).

In this way, the mutational substitutions that accumulate in one population are indifferent to the substitutions that accumulate in other populations, but only so long as the two populations do not share alleles with one another through gene flow. When the populations of incipient species do intermingle genetically, then developmental programs in hybrid individuals must confront whether alleles derived from such ‘out-of-sight, out-of-mind’ evolution interact in a way that allows ontogeny to proceed normally (Box 1A,B). Incompatible combinations of alleles in hybrid individuals at two or more loci can be thought of as negative epistasis or as a kind of non-linear genetic perturbation. Put this way, the genetics of speciation sound eerily similar to notions of cryptic genetic variation (Wade et al., 1997; Ledón-Rettig et al., 2014; Paaby and Rockman, 2014) and developmental system drift in the literature on the evolution of development (True and Haag, 2001; Pavlicev and Wagner, 2012; Tulchinsky et al., 2014a). The regulatory pathways that govern development are an important component of how and whether divergence will lead to genetic incompatibility and speciation (Porter and Johnson, 2002). Our aim here is to draw these connections more explicitly, to place speciation in a firmer developmental context (Sucena and Stern, 2000; Johnson and Porter, 2001), and to emphasize the relevance of the evolution of reproductive isolation to problems in developmental biology.

Distinguishing detectable versus undetectable differences in phenotype between species is key to thinking about divergence in the genes and genetic networks that underlie the developmental pathways that create an organism. Detectable differences in the development of traits between species, often due to adaptive divergence, clearly involve genetic changes. The converse, however, is not true: conservation of phenotype does not imply conservation of genetic controls. For example, stabilizing selection on expression levels leads to similar expression levels between species, despite evidence of widespread compensatory effects of molecular evolution to both cis- and trans-acting regulators of gene expression (Mack and Nachman, 2017). This idea of genetic change despite phenotypic stasis is known in evo-devo as developmental system drift (DSD), and in population genetics as evolution along a ridge in a fitness landscape (Box 1C,D). Adaptive evolution at the molecular level also can contribute to DSD for a given trait, particularly if the genetic changes affect fitness through pleiotropic effects with phenotypic consequences for only a subset of traits or from subsequent compensatory evolution (Box 1C).

When molecular evolutionary differences between species involve two or more loci, they may not interact properly in hybrid individuals that have copies of both genomes. This can happen regardless of whether those differences induce detectable phenotypic differences that distinguish the species and regardless of whether selection or genetic drift caused their fixation. Each mutational substitution that distinguishes species in protein coding sequence or gene regulatory control thus has some probability of contributing to the formation of a DMI in hybrid individuals. The more substitutions, the more chances for incompatibilities to create hybrid dysfunction at some point in development (Orr, 1995). In this way, inter-species hybrids can reveal genetic divergence in the control of even those developmental programs that yield seemingly equivalent phenotypic outputs. Hybrids can also reveal the incidence and role of different kinds of changes, such as cis- vs trans-acting regulators of gene expression (Box 2; Wittkopp et al., 2004; Mack et al., 2016).

Developmental disruption and regulatory inference in inter-species hybrids.

Analysis of inter-species hybrids can reveal genome-wide mechanisms of gene regulatory change.

(A) Comparison of allele-specific expression ratios in F1 hybrids to expression ratios for orthologous genes in parental individuals allows inference of expression changes due to local cis-acting regulatory differences, distant trans-acting changes, or distinct ways that both cis- and trans-regulatory divergence can jointly influence gene expression. (B) Another type of expression comparison between F1 hybrid individuals and parental species can characterize the dominance in allele-specific expression, with misexpression in hybrids reflecting disproportionately high (overdominant) or low (underdominant) gene expression. (C) The inference of misexpression in hybrids at the transcriptional level may be buffered at the translational level, leading to more severe misexpression of the transcriptome than proteome. Other approaches to using inter-species hybrids for deciphering the genetics of divergence in development include quantitative trait locus (QTL) mapping. Inter-species QTL analyses can incorporate screens with deletion libraries (Masly and Presgraves, 2007), recombinant inbred line (RIL) panels (Bedinger et al., 2011), near isogenic line (NIL) or introgression line panels (Guerrero et al., 2017), or multigeneration selection approaches such as X-QTL mapping (Ehrenreich et al., 2010). Panels A and B redrawn from data for Drosophila flies in McManus et al., 2010. Panel C data for Saccharomyces yeast from Wang et al., 2015.

Divergence in developmental genetic programs and intrinsic post-zygotic reproductive isolation in the speciation process are thus close conceptual brethren, despite their largely separate research traditions. But how do features of genetic networks and evolutionary forces intersect to create such profound developmental genetic divergence revealed as DMIs in hybrids? Are some stages of ontogeny predisposed to genetic architectures and evolutionary pressures to be more likely to yield dysfunctional development in hybrid organisms? Applying a timepiece metaphor: what gears and springs in the genetic clockwork will set the developmental alarm clock to ring at one time versus another during ontogeny to signal speciation? The more molecular evolutionary change, the more likely are genetic incompatibilities to manifest in inter-species hybrids. We must therefore set our goal to define the factors that make some genes evolve more quickly than others. The answers will help us to determine ways in which the molecular evolution that underpins development is predictable and, consequently, in what ways the genetics of speciation is predictable.

To address this goal, here we first consider how molecular evolution is influenced by the properties of genes and genetic networks, such as pleiotropy, modularity, robustness, and cis/trans regulation. These are the ‘gears and springs’ in the genetic architecture of development. We then explore how genetic architecture may sensitize some phases of ontogeny disproportionately to disruptive effects of misexpression in genetic networks. By integrating genetic architecture and ontogenetic timing, we arrive at distinct predictions for how molecular evolution and hybrid dysfunction manifest across developmental time. Finally, we summarize the literature on these issues for three case study systems (Caenorhabditis nematodes, Drosophila fruit flies, and anuran Bufo and Xenopus). In the present state of the field, we find that few general answers emerge from these factors considered in isolation, motivating deeper attention from theory and empirical study.

Speciation and development: divergence in genetic networks

As we explore how the evolution of development intersects with speciation, it is valuable to consider some key aspects of genetic architecture from the perspective of multi-gene networks that, in turn, control organismal development (Johnson and Porter, 2007; Palmer and Feldman, 2009). Here, we focus on how pleiotropy, network modularity, and robustness can influence the molecular evolution of coding sequences and non-coding regulatory elements. These links will help ground our expectations for incorporating ontogenetic time into our thinking to then consider predictions for molecular evolution and the production of incompatibilities in genetic networks of hybrid individuals.

Pleiotropic roles and effects

The mapping of genotype on phenotype and fitness leads us to predict that evolution will proceed along genetic lines of least resistance (Schluter, 1996). Genetic ‘resistance’ to evolutionary change is affected by the mutability and covariation of traits within and across developmental stages, in addition to natural selection. That is, the net rate of evolution integrates the likelihood and accumulation of mutational input (genetic variance) with the consequences of genetic variants for the ensemble of traits that comprise an organism (genetic covariances between traits) and fitness (natural selection). From the perspective of developmental biologists, this means that genetic hot-spots of evolutionary change over the long term ought to favor mutations with minimally pleiotropic effects (Carroll, 2005; Stern and Orgogozo, 2008), a form of developmental bias. As more molecular evolutionary change occurs, genetic incompatibilities become more likely to form from dysfunctional interactions between diverged sequences in inter-species hybrids. Our goal, then, is to enumerate the factors that facilitate molecular evolution to make some genes evolve more quickly than others. The structure of genes and gene regulatory networks provide clues to predict what these factors are (Garfield et al., 2013). These clues about molecular evolution can help us in thinking how to connect the temporal dynamics of developmental genetic networks to genetic incompatibilities between species.

The extent of pleiotropic effects induced by mutation to a gene is a major determinant of the likelihood of evolutionary change to that gene (Carroll, 2005; Stern and Orgogozo, 2008). This logic follows from the idea that most traits experience stabilizing selection on short timescales, and so the effects of mutation on most traits will be detrimental and offset any fitness benefit conferred to a particular trait by a mutation. On longer timescales, conserved traits are thought to track a moving but bounded fitness landscape (Estes and Arnold, 2007), which can facilitate compensatory evolution among the collection of interacting genes of the genetic network to result in sequence divergence and reproductive isolation between populations (Haag and Molla, 2005; Tulchinsky et al., 2014a; Tulchinsky et al., 2014b). Consequently, we might expect that more rapid evolution should occur in genes with fewer pleiotropic roles or for classes of changes that perturb fewer of a gene’s roles. The position of a gene within gene regulatory networks defines the pleiotropic roles it plays in development. Genes in more central and highly connected network positions can influence a greater fraction of the genome, and so changes to those genes have the potential to induce greater pleiotropic effects (Promislow, 2004; Hahn and Kern, 2005; He and Zhang, 2006; Kahali et al., 2009). However, the precise relationship between centrality and pleiotropy remains complex (Siegal et al., 2006). And, despite the idea that network ‘kernels’ may perdure through evolution (Erwin and Davidson, 2009), the molecular interactions among kernel components may nevertheless coevolve in ways that could generate DMIs (Haag, 2007; Tulchinsky et al., 2014a). Genes that are expressed at high levels also are more likely to have a large number of interaction partners and therefore to occupy central network positions (Bloom and Adami, 2003).

Coding sequence changes can affect protein activity everywhere they get expressed, and so ought to have greater potential for pleiotropic effects than regulatory changes (Wray et al., 2003; Carroll, 2005; Stern and Orgogozo, 2008). Alterations of cis-regulatory elements, in contrast, generally influence only a subset of a gene’s expression, and so will exert the fewest pleiotropic effects relative to coding sequence and trans-regulatory changes. Consequently, trans-regulatory evolution and coding sequence divergence ought generally to be slower for genes with high centrality and connectivity in genetic networks. We should expect such genes to have especially high ratios of cis:trans regulatory divergence, provided that those genes also have reasonably complex arrays of cis-regulatory elements controlling their expression that can facilitate divergence. Because trans-acting regulation tends to show greater condition-dependence than does cis-regulation (Smith and Kruglyak, 2008; Tirosh et al., 2009), however, pleiotropy-based predictions for their evolution may be tempered by the fact that varying conditions experienced by organisms could limit the net negative pleiotropic effects of trans-regulatory changes. Moreover, trans-regulatory evolution may often be coupled to subsequent cis-regulatory compensatory evolution that ameliorates negative pleiotropic effects of a population having fixed a beneficial trans-regulatory change. Such compensatory evolution involving both cis- and trans-acting regulatory evolution may be an especially important contributor to DMI formation (Landry et al., 2005; Ortíz-Barrientos et al., 2006; Takahasi et al., 2011; Mack et al., 2016).

Experiments show that trans-acting regulatory factors contribute disproportionately to segregating genetic variation for gene expression levels within species, whereas cis-acting regulatory differences disproportionately underlie fixed genetic differences in expression between species (Wittkopp and Kalay, 2012). The loci controlling DMIs also can be polymorphic or fixed (Cutter, 2012), and so to the extent that regulatory evolution is responsible for creating DMIs, we ought to expect polymorphic DMIs to often be controlled by trans-acting factors. Because trans-regulatory mutations tend to arise more readily due to larger mutational target size in the genome, and to be more recessive, than cis-acting mutations (Landry et al., 2007; Gruber et al., 2012), they are expected to be a common and persistent kind of polymorphism within populations.

Modularity of genetic network architecture

Modularity of genetic networks reduces the scope for genes to have highly pleiotropic roles (Wagner and Altenberg, 1996). By constraining the neighborhood of partner interactions, greater modularity limits the pleiotropic effects of changes to gene expression or function for any given gene. Like tissue-specific expression in the spatial modularity of gene networks (Larracuente et al., 2008; Packer et al., 2019), we can also consider temporal modularity for those portions of gene regulatory networks that show stage-specific expression. Modularity of genetic network structure will increase over ontogeny as cells differentiate and tissues establish autonomous regulatory programs (Packer et al., 2019). Some transient phases of development are thought to reverse this trend, however, as when tissues integrate during gastrulation or reorganize during metamorphosis (Chin-Sang and Chisholm, 2000). Consequently, if pleiotropy at the scale of the whole organism constrains evolution of individual genes, then we might expect faster molecular evolution for genes that experience highly modular genetic network structures due to spatially- or temporally restricted expression. This logic provides one way to frame the ‘early conservation’ and ‘hourglass’ models in the evolution of development (Figure 1).

Figure 1.

Predictions and hypotheses in the evolution of ontogeny and reproductive isolation.

(A–C) Models from population genetics and evo-devo suppose that some modes of natural selection may be more potent at particular life stages, as described in the main text. However, not all models make clear predictions about the ontogenetic dynamism of selection for all selection modes (purifying vs adaptive vs neutral) or for all stages of development. (D–E) Differential incidence of selection across development may be mediated by genetic architecture in terms of the pleiotropic effects of genetic changes and how that translates into the robustness of fitness-related traits. When fitness is disproportionately robust to changes to genes expressed at a given stage, then that stage will be more likely to accumulation cryptic genetic variation (CV) within species, divergence between species as developmental system drift (DSD), and to result in production of Dobzhansky-Muller incompatibilities (DMIs) in inter-species hybrid individuals.

Despite a trend of increasing genetic network modularity from single-celled zygote to adult, some genes will experience the modularity more than others at any given place or time in ontogeny. Genes with greater breadth of expression across tissues, and genes with greater temporal persistence of expression, will experience more scope to contribute to different genetic networks and their phenotypic outputs. Consequently, genes most likely to produce pleiotropic effects when genetically perturbed are those expressed at high levels across many tissues throughout a long span of developmental time (high expression breadth and low temporal expression bias; Larracuente et al., 2008; Coronado-Zamora et al., 2019). Even a modest role in any given tissue or at any single stage of development, however, may compound when integrated over space and time.

The appropriate weighting of the importance of spatial versus temporal expression breadth and modularity, however, is not entirely obvious. Therefore, empirically, it will be valuable to partition the ontogenetic trajectories of genes with similar expression breadths across tissues, or, complementarily, partition the spatial profiles of expression for genes with similar ontogenetic expression dynamics (Box 3). The transcriptome analysis of distinct cell lineages, as conducted for C. elegans precursors of endoderm, mesoderm and ectoderm (Hashimshony et al., 2015) or for single-cell transcriptomes throughout embryogenesis (Hashimshony et al., 2012; Tintori et al., 2016; Packer et al., 2019), provides one intriguing scheme for approaching this issue. Another approach could incorporate developmental time into transcriptome analysis of serial-sectioned samples to access four-dimensional expression profiles through developmental space and time (Ebbing et al., 2018).

Measuring molecular divergence across ontogeny.

Incorporating ontogeny into molecular evolution requires developmental time series of gene expression, which necessarily involves more experimental effort and sophistication of analysis than studies using a single developmental timepoint. This idea can be explored in several ways. First, one may quantify divergence in gene expression as a phenotypic output that closely maps to genotype. Studies in diverse organisms have compared transcriptome profiles over ontogeny with this approach (Table 1). Second, the contrast of gene expression levels between two parental species and their F1 hybrids provides a way to infer whether changes in expression result from cis- or trans-regulatory evolution (Box 1; Wittkopp et al., 2004; Signor and Nuzhdin, 2018). This approach has challenges and opportunities: environmental conditions also can influence transcriptome expression and the balance of cis:trans controls of the set of expressed genes (Tirosh et al., 2009), and it may be difficult to deconvolve maternal from zygotic trans-effects in F1 embryos. The technique also has not yet been applied in a developmentally dynamic way, making it ripe for future studies. Third, rates of coding sequence evolution (dN/dS or KA/KS) for genes expressed differentially across development provide a means of assessing selection on the encoded protein sequences used in genetic networks. A simple way of testing for trends in purifying selection and positive selection across ontogeny is to compare the average dN/dS value for the set of genes with peak expression at a given developmental timepoint (Cutter and Ward, 2005; Davis et al., 2005) or an expression-weighted mean dN/dS value for all genes expressed at a given timepoint (Quint et al., 2012; Liu and Robinson-Rechavi, 2018b). When a multi-species phylogeny is used, or if polymorphism data are incorporated in a McDonald-Kreitman testing framework, then positive selection may be distinguished from relaxed purifying selection (Liu and Robinson-Rechavi, 2018a). Genes also can be grouped according to patterns of co-expression to analyze for consistent differences in rates of molecular evolution among them or among parameters of mathematical functions that describe the temporal profiles of expression (Coronado-Zamora et al., 2019; Cutter et al., 2019). Spatial definition of expression can enhance such approaches using tissue- and sex-specificity (Kim et al., 2016), restricted cell-lineages (Hashimshony et al., 2015), or tomographic expression profiling (Ebbing et al., 2018). While mRNA expression is most commonly accessed, all these approaches can be extended to protein expression; protein levels appear to differ less between species than do transcriptomes (Leducq et al., 2012; Khan et al., 2013). Quantifying chromatin misregulation in a developmental time-course also could help discern the influence of epigenetic factors in hybrid dysfunction.

Table 1.

Studies in the literature that characterize expression profiles across development (Supplementary file 1).

	Number of studies (%)		Number of species (%)
Developmental stage analyzed
Embryogenesis only	58	39.5%	45	42.5%
Non-embryo stages	37	25.2%	17	16.0%
All ontogeny	51	34.7%	44	41.5%
Type of molecular time series data
Transcriptome	128	87.1%	106	100.0%
Proteome	7	4.8%	5	4.7%
Other	14	9.5%	12	11.3%
Inter-species divergence feature
Gene expression	32	21.8%	55	51.9%
DNA sequence	5	3.4%	14	13.2%
Non-comparative	99	67.3%	42	39.6%
Taxonomic group
Angiosperm	16	10.9%
Arthropoda	26	17.7%
Chordata	55	37.4%
Nematoda	10	6.8%
Other group only	30	20.4%
Multiple groups	10	6.8%

Robustness to genetic change in gene regulatory networks

Selection on traits within a species often is stabilizing, favoring a particular phenotypic value and meaning that factors that lead to alternative phenotypic values will have lower fitness. A genetically ‘robust’ phenotype is insensitive to genetic perturbations, which permits an organism to produce the same phenotype and to maintain fitness despite mutational disruption to a gene or genetic network that contributes to the trait’s developmental program (Félix and Wagner, 2008). That is, a greater fraction of mutations to genes in more robust networks have effectively neutral effects (Ohta, 2011). A gene network’s robustness therefore ought to influence the mutations that can accumulate, which in turn affects rates of molecular evolutionary divergence and the manifestation of genetic incompatibilities in hybrids. Differences across gene networks in genetic robustness could arise from differential selection on robustness itself (adaptive robustness or canalization) or could simply arise as a byproduct of differences in their epistatic genetic architectures (Hermisson and Wagner, 2004). Genes and genetic networks that are more robust to genetic perturbation can be thought of as larger phenotypic capacitors (Paaby and Testa, 2018), letting more cryptic genetic divergence in developmental controls accumulate to be revealed upon the formation of inter-species hybrids. This effect is analogous to how specific environmental circumstances or genetic backgrounds can expose so-called cryptic genetic variation or conditionally-neutral variants within a species (Ledón-Rettig et al., 2014; Paaby and Rockman, 2014).

As a consequence, traits with greater genetic robustness ought to accumulate more cryptic genetic variation within a species that can then can contribute to greater developmental system drift between species as substitutions accrue (Félix and Wagner, 2008). This process will, in turn, lead separate populations to travel independently along a fitness ridge in the landscape of genotype space (Box 1D; Félix and Wagner, 2008; Fragata et al., 2019). More robust networks are more likely to contain more cryptic genetic variation, but the presence of cryptic genetic variation alone is insufficient to conclude whether or not selection directly favored greater robustness, and characterizing the robustness of genetic networks empirically can be difficult. The selective regime that shapes the fitness landscape will influence the speed at which developmental system drift can evolve for any given trait (Tulchinsky et al., 2014b). DSD is especially likely to result from directional selection on traits, perhaps facilitated by a role for cryptic genetic variation in adaptation to environmental shifts (Félix and Wagner, 2008), and from selection that drives molecular co-evolution among genetic elements (Tulchinsky et al., 2014b).

Phenotypes with a larger underlying genetic network are expected to be more robust to perturbation (Fragata et al., 2019), implying that they should more readily tolerate the accumulation of genetic changes. Genes holding more central positions in genetic networks also appear to be associated with greater robustness upon perturbation (Proulx et al., 2007), although this effect may depend more on expression level than connectivity per se (Siegal et al., 2006). This means that changes to any given gene in a network that is shared with other networks, and that confers a beneficial effect to that other network, are less likely to manifest negative pleiotropic effects for those genetic networks that are larger (Pavlicev and Wagner, 2012). Consequently, cryptic variation and DSD may be disproportionately prevalent in developmental programs comprised of larger genetic networks. Genetic networks with low modularity, by definition, have all constituent genes interacting, and networks with low modularity indeed are most robust to genetic perturbation (Tran and Kwon, 2013).

Recall, however, that we argued how such low modularity conditions maximizes the pleiotropic roles of genes, which should impede rather than tolerate molecular evolution. These effects thus appear to present a contradiction. The resolution may lie in assumptions about network size, such that the implications of low modularity for pleiotropy may be secondary to robustness for large genetic networks. Analysis of networks in high-resolution time series of expression showing distinct tissues with distinct transcriptomic profiles that vary in size may help in testing this idea (Packer et al., 2019). The resolution may also depend on whether we consider coding or regulatory sequence evolution, as coding sequence evolution does not appear to depend strongly on how connected the gene is in the genetic network (Jordan et al., 2003; Batada et al., 2006). The way that connections define the network output also can be crucial. Simple connectivity metrics may just poorly summarize functional and evolutionary properties of networks, making the contradiction more apparent than real (Siegal et al., 2006). For example, genetic networks with switch-like behavior may be more capable of accumulating cryptic variation, as appears to occur disproportionately among genes with early zygotic expression in sea urchin embryogenesis (Garfield et al., 2013).

Ontogenetic stages that have an especially high incidence of mutationally-robust phenotypes within species should thus be especially prone to experience developmental system drift. Such DSD can get revealed in the novel genomic environment of inter-species hybrids, taking the form of DMIs. A counter-argument to this idea is that mutationally-robust traits might also be predisposed to being robust to the genomic perturbation experienced in hybrid individuals, and therefore be less likely to show dysfunction. The logic of this counter-argument might be more pertinent at the very early stages of divergence between populations, becoming less and less applicable as substitutions accrue in the speciation process. It will be interesting for future research to resolve this question. In general, it remains unclear how network robustness within species versus between species might scale differently with the number of genetic changes. Measurements of mutational variance for traits shared across ontogenetic stages provide another possible way to test for stage-dependent changes in mutational robustness, for example, using mutation accumulation experiments that, so far, point to this possibility (Farhadifar et al., 2016; Zalts and Yanai, 2017). It will be interesting to see in future studies whether the genetic robustness of traits truly can predict the incidence of cryptic genetic variation, developmental system drift, and the propensity to reveal DMIs in hybrids.

Degeneracy via partially redundant genetic network pathways provides another way that phenotypes could be robust to genetic perturbation, potentially allowing DSD to accumulate without yielding DMIs in F1 hybrids. Indeed, genes associated with redundant networks evolve more quickly (Wagner, 2000; Kitami and Nadeau, 2002). Robustness mediated by distributed, rather than redundant (Félix and Wagner, 2008), genetic networks may thus more readily experience DSD in a way that would foster DMIs. Moreover, translational buffering of gene expression leads to the inference that protein production shows much less misexpression than do mRNA transcripts in F1 hybrids (Leducq et al., 2012; Khan et al., 2013; Artieri and Fraser, 2014; McManus et al., 2014). And, it is fitness as a phenotype that is the ultimate readout of robustness. Consequently, genetic networks important for development may be less disrupted in hybrids than transcriptome analyses might otherwise suggest. Overall, discerning the relative incidence of distinct mechanisms conferring phenotypic robustness will be important in defining whether different genetic networks and different stages of development will be more or less likely to contribute to post-zygotic reproductive isolation as divergence between species accumulates.

Evolvability of distinct genetic components

Mode of selection

Traits and sequences with greater propensity to diverge are said to be more evolvable (Wagner and Zhang, 2011). One way to detect evolvability is when traits and sequences diverge as a result of natural selection, because adaptive divergence in phenotypes between species makes it easy to conclude that there must have been changes to the underlying genetic networks. Such selection can involve multiple changes on an adaptive path, for example, if first a large-effect trans-regulatory mutation gets fixed and is followed by subsequent compensatory cis-regulatory substitutions that ameliorate suboptimal pleiotropic effects of the initial trans-acting substitution (Box 1C; Goncalves et al., 2012); similar logic can also apply to coding sequence changes (Clark et al., 2009). Moreover, positive selection and multi-locus antagonistic coevolution will drive molecular evolution that is much more rapid than will genetic drift. Repeated adaptive evolution of orthologous genes in disparate lineages is implicated as a general feature of some types of phenotypic change, including aggregate effects of multiple mutations affecting the same locus (Stern and Orgogozo, 2008; Martin and Orgogozo, 2013). The genes contributing to such adaptive divergence may thus be especially likely to contribute to post-zygotic DMIs (Presgraves, 2010b), in addition to pre-mating, gametic, or ecological reproductive isolation barriers.

Alternately, some genes may evolve faster than others simply because they experience weaker purifying selection or greater mutational target size. The evolutionary accumulation of changes in this way is rate-limited by the input of mutation and the speed of genetic drift, and so will be faster in species with small population sizes.

Haldane’s rule as developmentally predictable evolution

One developmentally predictable rule in evolution says something about sex: Haldane’s rule, in which the heterogametic sex is more likely to suffer inviability or sterility in inter-species hybrids (Haldane, 1922; Delph and Demuth, 2016). Dominance theory provides one explanation for this pattern: incompatibility loci linked to sex chromosomes will reveal themselves in F1 hybrids disproportionately for the sex that has only one copy of a given sex chromosome (Turelli and Orr, 1995). This tells us that genomic compartmentalization of how traits are genetically encoded also is important for the implications of molecular evolution; when genes contributing to a given developmental process are biased in genomic location, it may confer greater or lesser tendency to yield DMIs in hybrids.

When males are heterogametic (X/Y or X/O sex-determination), other factors may also contribute to Haldane’s rule (Delph and Demuth, 2016). Genes that control male-biased traits may evolve especially rapidly due to sexual selection or mutational biases (‘faster male’ theory), or gene regulatory networks controlling male traits may be unusually sensitive to genetic perturbation in hybrids (‘fragile males’), or genes linked the X-chromosome may evolve especially fast (‘faster X’) (Charlesworth et al., 1987; Wu and Davis, 1993). Consequently, X-linked genes may show distinct patterns of misexpression in hybrids (Moehring et al., 2007; Turner et al., 2014; Civetta, 2016; Sanchez-Ramirez et al., 2020). Developmental pathways related to spermatogenesis seem especially sensitive to disruption in hybrids, perhaps being predisposed to disproportionate compensatory cis-trans regulatory coevolution (Mack et al., 2016; Mack and Nachman, 2017). Sex-limited genetic networks also may differ for males versus females in size, location of genomic encoding, or predominant mechanism of regulatory control, and so influence the relative accumulation of cryptic genetic variation and developmental system drift.

Coding vs regulatory evolution

Both coding sequences and regulatory sequences diverge between species, despite the fact that most selection on each of them is expected to be purifying (Casillas et al., 2007). The rate of evolution for a coding sequence and its cis-regulatory regions, however, correlate only weakly (Castillo-Davis, 2004; Liao and Zhang, 2006; Tirosh and Barkai, 2008). This observation implies that the strength and mode of selection affecting mutations to protein structure tells us little about the strength and mode of selection on mutations affecting expression, and vice versa. Consequently, molecular evolution associated with ontogenetic timing may show distinct patterns for coding and regulatory sequences, and so also yield different implications for when DMIs manifest over ontogeny.

The type of regulatory change may be key, however, in understanding the relation between regulatory and coding sequence divergence (i.e. evolvability), as well as to the likelihood of contributing to a DMI. For example, genes showing evidence of trans-regulatory divergence appear to evolve slower in their coding sequences than do genes with cis-regulatory divergence (Goncalves et al., 2012). And, DMIs due to misregulation generally involve non-compensatory evolution of both cis- and trans-regulators of expression (Ortíz-Barrientos et al., 2006; Mack and Nachman, 2017). Specifically, regulatory divergence involving reinforcing cis+trans changes (thought to result from directional selection within species) may actually be less likely to create DMIs than changes with opposing cis-trans effects (thought to result from stabilizing selection within species) (Mack et al., 2016) (but see Tulchinsky et al., 2014b). Consequently, positive selection affecting regulatory versus coding sequences may differ in their propensity to yield DMIs.

Might the functional role of genes also represent an axis of predictability to factors that instigate DMIs in inter-species hybrids? The relatively few known ‘speciation genes’ in animals give little clue to whether particular molecular functions may be predisposed to involvement in inter-species incompatibilities (Blackman, 2016). At a coarse level, because DMIs require interaction, we should expect sequences that affect interactions between DNA, RNA, and proteins to be more prevalent among DMI loci than, say, enzymes that interact predominantly with metabolites. From the perspective of gene regulatory networks, transcription factors experience faster sequence evolution than other genes in the genome (Gilad et al., 2006; Haerty et al., 2008). piRNA genes also turnover rapidly and are implicated in hybrid dysfunction (Assis and Kondrashov, 2009; Bagijn et al., 2012; Kelleher et al., 2012), along with other classes of regulatory endogenous small RNAs (Li et al., 2016). Genes that act cell non-autonomously, as for secreted proteins or diffusible signaling molecules, may tend to evolve slower if they also tend to exhibit greater expression breadth and lower modularity, unless they are predisposed to co-evolutionary dynamics. Reproductive isolation between species need not only involve classic regulators of development, however, as attested by incompatibilities that commonly seem to involve chromosome segregation and cyto-nuclear interactions (Blackman, 2016; Lima et al., 2019). Despite the fundamental role of the mitochondrial genome in energy metabolism, mitochondrial genes experience adaptive molecular evolution (Bazin et al., 2006) and mito-nuclear incompatibilities can provide important reproductive barriers between species (Hill, 2015; Lima et al., 2019). To the extent that some developmental stages may be more sensitive to the disruption of chromosome segregation or mitochondrial function, such as phases of heightened cell division or metabolism, such developmentally tangential genetic pathways might nevertheless contribute to ontogenetic patterns of hybrid dysfunction.

Gene duplication and gene origination

Gene duplication is a powerful factor in the evolution of phenotypic novelty (Kaessmann, 2010). Moreover, regulatory or structural subfunctionalization in the evolution of gene copies following duplication could lead to dysfunction in hybrids and so contribute to a DMI (Mack and Nachman, 2017). Similarly, a DMI could result from ‘divergent resolution’ in the loss of alternate copies in different species for functionally equivalent gene duplicates (Lynch and Force, 2000). The de novo origin of new genes also can instigate species differences in gene network structure (Neme and Tautz, 2014), as could divergent use of alternative splice forms of a gene (Ortíz-Barrientos et al., 2006). Thus, the phenotypic evolvability promoted by gene duplication and gene origination also confers the potential to induce reproductive isolation as gene copies evolve through distinct trajectories in different evolutionary lineages.

Ontogenetic timing of genetic networks and molecular evolution

Gene expression profiles and genetic network architecture are dynamic over the course of ontogeny. This dynamism suggests that molecular evolution could differ for the distinct subsets of genes associated with different stages of development (Box 3). Importantly, this molecular divergence may or may not correspond to divergence in organismal phenotypes (True and Haag, 2001), despite the common emphasis on how changes to gene regulatory networks alter the development of phenotypes (Erwin and Davidson, 2009). Stages especially prone to rapid molecular evolution, of coding sequences or regulatory elements, may also translate into a greater incidence of DMIs. Ontogenetic predispositions toward sequence evolution may therefore confer predictable ontogenetic detection of reproductive isolation between species, letting us know when to listen most carefully for the developmental alarm clock of speciation. In our search of the literature, we identified 147 studies involving 106 species that quantified transcriptome, proteome, or related molecular data across multiple developmental stages (Table 1, Supplementary file 1). We found 40% of studies to focus on embryogenesis only. Over 60% of the studies included a comparative analysis of expression divergence between species (52%) or of DNA sequence divergence (13%). Integrating such studies with developmental time courses of hybrid dysfunction would provide a powerful collection of empirical tests for how incompatibilities in genetic networks arise through ontogeny. In the meantime, ideas from both development and population genetics lead to several partially-distinct predictions for the evolution of genes expressed differentially across ontogeny, which we enumerate below (Figure 1). We conclude that a conceptual gap is the lack of a modeling framework that integrates these disparate perspectives on the role of ontogenetic timing in molecular evolution.

Early conservation model

The early conservation model (‘von Baer’s third law’) derives from the fact that time is unidirectional, so changes early in development may cascade catastrophically as cell division and differentiation proceeds. Consequently, genes expressed early in development would experience stronger purifying selection, with slower molecular and phenotypic evolution early in ontogeny (Kalinka and Tomancak, 2012). Some perspectives on gene regulatory network structure follow the spirit of this view, as well (Erwin and Davidson, 2009). This idea would also be consistent with genes and genetic networks expressed in early development having lower robustness, modularity, and evolvability, and with greater pleiotropy when perturbed (Figure 1). The early conservation model is mostly applied to embryogenesis, with some empirical support from vertebrates based on gene expression divergence (Roux and Robinson-Rechavi, 2008; Irie and Kuratani, 2014). In principle, however, the logic of this model applies to the entirety of ontogeny. If the likelihood of DMIs scales with rate of evolution, then we ought to expect post-zygotic reproductive isolation to be more likely to manifest later in development and that species at earlier stages of divergence would manifest defects later in development.

A caveat about inferring the developmental timing of hybrid dysfunction is that early-acting incompatibilities may preclude detection of DMIs that would otherwise be revealed later (see uniform chance model, below). Disproportionate observation of early-acting developmental defects in hybrids thus may not imply that molecular evolution disproportionately accrues for genes in early-acting developmental programs. This ‘pull of the early’ represents a general challenge in characterizing the profile across ontogeny of disrupted genetic networks that confer post-zygotic reproductive isolation. There are at least four ways to potentially address this issue empirically: (i) assess distinct species comparisons from different phylogenetic depths, (ii) exploit partial penetrance of F1 hybrid dysfunction (Bundus et al., 2015), (iii) use early-acting hybrid-rescue genotypes to test for late-acting hybrid dysfunction (cf. Hmr in Drosophila [Hutter et al., 1990]), or (iv) focus on misregulation of gene expression over ontogeny for species with relatively weak post-zygotic isolation (i.e. without catastrophic effects in early life stages).

Uniform chance model

We suggest that the ‘uniform chance’ model represents a null model for the manifestation of hybrid incompatibilities over the course of ontogeny. If molecular evolution is unbiased with respect to the timing of expression of genes that have diverged between species, then every point in development can be considered to have an equal chance of expressing a DMI to yield hybrid disruption of a tissue type or termination of development (Figure 1). This could arise from each stage adapting continuously to the ‘habitat’ it experiences distinctively from other stages, whether inside an egg (or uterus), juvenile environmental circumstances, or post-metamorphosis adulthood. Consequently, species pairs with greater overall divergence would be more likely to have hybrids that terminate development earlier in ontogeny; the probability that the earliest developmental stage avoids the effects of DMIs declines exponentially (p with n DMIs each with probability p of not occurring in the earliest stage). This null model therefore predicts a negative relationship of the genetic distance between species and the terminal stage to which hybrids develop, and predicts that later stages will have more tissue types exhibiting dysfunction in hybrids. Note that these predictions overlap with several other models that we describe, but do not depend on differential selection pressures, pleiotropy, or network features for genes expressed at different times in development. This model, however, predicts no association of molecular evolutionary rates for genes (or the mode of selection) with the developmental timing of their expression.

Mutation accumulation model of aging

The mutation accumulation model of aging and senescence predicts that genes and genetic networks expressed earlier in development will evolve more slowly (Promislow and Tatar, 1998; Partridge, 2001). Here, the logic and developmental timing, however, is different from the early conservation model: the reproductive value of individuals declines following the onset of reproductive maturity, meaning that purifying selection is weaker against deleterious mutations that affect adult phenotypes relative to embryonic and juvenile phenotypes (Medawar, 1952; Flatt and Schmidt, 2009). This perspective presumes similarly strong purifying selection for all genes expressed prior to maturity, and makes no specific prediction about positive selection across ontogeny (Figure 1). Not explicitly formulated by the mutation accumulation model, however, is whether the strong early-life selection would indirectly favor mutational robustness to embryonic and juvenile developmental programs. Such robustness could facilitate greater molecular evolution of early-expressed genes than otherwise anticipated, which also could facilitate the accumulation of developmental system drift.

Antagonistic pleiotropy model of aging

Gain in fitness from traits that increase survival or reproduction early in life will experience disproportionate selection pressure. Consequently, positive selection may more fiercely favor beneficial mutations to genes expressed early in life, irrespective of negative pleiotropic consequences on fitness-related traits expressed later on. This idea is the essence of the antagonistic pleiotropy model for the evolution of aging and senescence (Williams, 1957; Partridge, 2001; Flatt and Schmidt, 2009). In terms of molecular evolution, it ought to yield a signature of excess positive selection on genes first expressed at pre-reproductive stages (Figure 1). However, subsequent selection may lead to compensatory evolution for genes expressed late in life to ameliorate negative pleiotropic effects. This sequential evolutionary series of events could produce developmental system drift, especially for late-life traits. Regulatory divergence between species often reflects compensatory changes in cis and trans, with trans-acting changes being more likely to exert pleiotropic effects and arising at a higher rate (Gruber et al., 2012). It will be interesting to determine whether the developmental timing of cis versus trans changes are biased over ontogeny in a way that would be consistent with the antagonistic pleiotropy model, that is, disproportionate cis-regulatory changes that modulate late-acting expression.

Hourglass model

The hourglass model is a pattern in search of a mechanism that derives from classic phenotypic observations of a ‘phylotypic stage’ in mid-embryogenesis, a point of greatest morphological similarity across organisms. Profiles of gene expression divergence between species, a molecular phenotype, also are often interpreted to be consistent with this pattern (Raff, 1996; Kalinka and Tomancak, 2012). In terms of DNA sequence evolution, it is less obvious why genes expressed in mid-embryogenesis would exhibit slowest rates of sequence evolution, whether due to stronger purifying selection or less frequent adaptive divergence (Figure 1). Ontogenetic trends in genetic network architecture may provide some intuition, as it is proposed that the phylotypic stage may coincide with developmental periods of intense integration of distinct cell lineages, often around gastrulation (Levin et al., 2012). Consequently, mid-embryogenesis may experience a shift toward fewer and larger genetic networks (low modularity) that leads to greater pleiotropic effects when perturbed. Changes in network structure related to switch-like or threshold traits also could contribute (Garfield et al., 2013). If such a genetic network architecture in mid-embryogenesis also confers greater within-species genetic robustness, then it may be especially prone to DSD and the production of DMIs and hybrid dysfunction in crosses between species. If, instead, it represents a point of low robustness to genetic perturbation, then DSD from genes expressed at earlier stages may induce a tipping point of dysfunction that manifests in hybrids during such a phase of integration.

The potential for distinctive properties in the ‘waist of the hourglass’ has drawn most scrutiny by researchers, but an alternative hypothesis might suppose that it is the earlier points in development that are unusual and require special explanation. This possibility has recently received greater theoretical and empirical attention (Demuth and Wade, 2007; Dapper and Wade, 2016; Zalts and Yanai, 2017; Coronado-Zamora et al., 2019). For example, theory predicts more rapid accumulation of mutations to genes acting in the earliest genetic networks that derive from maternal resources, leading to their unusually fast molecular evolution with the potential to drive DSD and DMIs (Demuth and Wade, 2007). Faster-than-anticipated molecular evolution at the very earliest stages of development might also arise from especially high robustness of genetic networks involving maternally deposited gene products, or from co-evolutionary dynamics associated with parent-offspring or other genomic conflicts (Brandvain and Haig, 2005; Crespi and Nosil, 2013). In such scenarios, the point in time that coincides with ‘waist of the hourglass’ might simply represent the onset of conditions consistent with the early conservation model.

Sex-biased selection

Selection on genes with sex-limited expression can lead to their more rapid evolution. This faster evolution can arise in two ways: sexual selection/conflict and weaker purifying selection. Sexual selection and sexual conflict that occurs within and between the sexes for reproductive adults can promote rapid phenotypic and molecular evolution via positive selection, especially in reproduction- and gamete-related traits (Swanson and Vacquier, 2002; Mank, 2017; Rowe et al., 2018). Thus, the timing of expression for genes that are part of the developmental programs that build such traits in late juvenile and adult stages would be expected to be most strongly impacted. This ontogenetic timing is similar to the mutation accumulation model, except that it should reflect adaptive evolution rather than weaker purifying selection and that it deals primarily with the subset of traits and genes associated with sexual interactions between individuals (Figure 1). Genes controlling gamete traits may be especially prone to rapid molecular evolution (Swanson and Vacquier, 2002), with the effects on hybrids potentially reflected in the earliest-stage zygotes. Speciation genetics research has incorporated this idea of rapid evolution of sex-limited genes into the ‘faster male’ theory for Haldane’s rule and into the explanation for why hybrid sterility often appears to arise earlier in the speciation process than does hybrid inviability (Coyne, 1989; Wu and Davis, 1993; Ortíz-Barrientos et al., 2006; Presgraves, 2010a).

In addition to the possibility of more prevalent positive selection on genes with sex-limited expression, they may also experience weaker purifying selection (Dapper and Wade, 2016). Weaker purifying selection means greater accumulation of divergence between species. This weaker efficacy of selection arises because half of the individuals in the population, one sex or the other, do not express the gene and so mask the effects of new detrimental mutations. This logic applies to both sex-limited gene expression for adult traits, as well as for maternally provisioned transcripts delivered to eggs (Demuth and Wade, 2007; Dapper and Wade, 2016). Rapid evolution of maternal-provisioning genes may lead to incompatibilities when they interact with zygotically-expressed genes in inter-species hybrids, potentially elevating the incidence of mid-embryonic hybrid dysfunction independently of any temporal changes in genetic network modularity or robustness. Genetic incompatibilities involving uniparentally inherited genetic factors also may lead to predictable asymmetries in reproductive isolation between species pairs (Turelli and Moyle, 2007), including roles for genes encoded on sex chromosomes and mitochondrial genomes (Bolnick et al., 2008; Barreto et al., 2018; Cutter, 2018).

Toward predictable rules in the evolution of development in speciation

To decipher how predictable molecular evolution and post-zygotic reproductive isolation might be, it is important to consider the ontogenetic context of organisms to identify trends in the tempo of gene expression, cell division, and differentiation (Table 1). The kind, number, and molecular evolutionary implications of genetic interactions also likely are sensitive to classes of embryogenesis programs (e.g., syncytial with late cellularization as in Drosophila and other insects, totipotent versus highly cell-autonomous cell lineage development as in C. elegans, nourishment by minimal vs very large yolk vs maternal tissues as seen in many insects, birds, and mammals) (Church et al., 2019). While here we emphasize animal systems, similar consideration of the evolutionary genetics of development and speciation for plant systems may also reveal valuable insights (Rieseberg and Blackman, 2010; Bedinger et al., 2011; Baack et al., 2015). Here, we distill in abbreviated form what is known for a few concrete motivating example systems (Caenorhabditis nematodes, Drosophila insects, Bufo toads) to frame these issues about how ontogenetic features link to molecular evolution and hybrid dysfunction. Surprisingly, few broad conclusions can be drawn from even these well-studied systems, and general principles await concerted research efforts. Nevertheless, these study systems present clear promise for future research to disentangle the constellation of causal contributing factors that link microevolutionary mechanisms to developmental programs and macroevolutionary patterns.

C. elegans ontogenetic profiles of gene expression, molecular evolution, and hybrid dysfunction

Taking C. elegans development as a point of reference, embryonic cell numbers grow approximately exponentially until ~520 cells upon which cell count increases relatively slowly to the 959 somatic and ~2000 germline cells that comprise the adult hermaphrodite animal (Giurumescu et al., 2012; Figure 2A). Cell lineages derived from five blastomeres exhibit distinct transcriptome profiles in ontogenetic timecourses (Hashimshony et al., 2015), with single-cell transcriptome sequencing of embryos in a time series showing even finer resolution (Packer et al., 2019). The number of genes expressed increases over most of embryogenesis, with relatively similar numbers of genes expressed at stages post-hatching (Boeck et al., 2016; Figure 2B). The change in identity of expressed genes, however, is greatest early in embryogenesis (primarily due to down-regulation) and late in embryogenesis (primarily due to up-regulation) (Boeck et al., 2016; Figure 2C). Gene connectivity peaks in early embryogenesis, declining until a spike in adulthood (Liu and Robinson-Rechavi, 2018a). Cellular development in embryos is conserved across species, being nearly indistinguishable at least to the 350-cell stage (Zhao et al., 2008; Memar et al., 2019). In hybrid crosses of most species pairs, however, embryos arrest around gastrulation (Baird and Yen, 2000; Baird and Seibert, 2013), near the time that expression divergence appears minimized across Caenorhabditis (Levin et al., 2012) and that expression-weighted coding sequence divergence is lowest (Cutter et al., 2019). Coding sequences with fastest molecular evolution show peak expression very early in embryogenesis or toward adulthood (Cutter and Ward, 2005; Cutter et al., 2019; Box 3). Hybrids of C. briggsae and C. nigoni show high embryonic inviability, but those genetically identical individuals that hatch successfully exhibit little larval mortality (Bundus et al., 2015). Gene misexpression is widespread in adults of both sexes for these hybrids (Sanchez-Ramirez et al., 2020) and misregulation of small-RNAs in hybrid genetic backgrounds is implicated in spermatogenic dysfunction (Li et al., 2016). Together, these observations suggest the possibility that developmental system drift may accrue more readily and generate hybrid incompatibility disproportionately at stages showing greatest selective constraint on gene expression levels and coding sequences.

Figure 2.

Developmental dynamics in C. elegans nematodes and in Bufo and Xenopus amphibians.

(A) C. elegans cell counts grow exponentially in early embryogenesis, before slowing later. Redrawn with data from Giurumescu et al., 2012 and wormatlas.org. (B) Gene expression changes dynamically over ontogeny in terms of number of genes expressed and the incidence of up- and down-regulated genes over time. Redrawn from Boeck et al., 2016. (C) Ontogenetic timing in the accumulation of reproductive isolation with genetic divergence for Bufo toads. Hybrid individuals between more closely-related species develop to later stages than do hybrids from distantly-related species pairs. Redrawn from Malone and Fontenot, 2008. (D–E) Xenopus gene expression level differentiation decreases over developmental time (yolk comprises ~½ of embryo volume; maternal-zygotic transition at stage ~8; yolk consumption begins around gastrulation). Redrawn from Yanai et al., 2011.

Drosophila divergence in sequences and transcriptomes over ontogeny

Coding sequences expressed across D. melanogaster development show slowest evolution in mid- to late-embryogenesis (Davis et al., 2005; Coronado-Zamora et al., 2019). The faster sequence evolution of maternally deposited transcripts in early embryogenesis is due disproportionately to non-adaptive divergence whereas it is adaptive divergence that is disproportionately implicated in faster sequence evolution of genes expressed post-embryonically (Coronado-Zamora et al., 2019; Figure 2E). Gene connectivity similarly shows a peak in early embryogenesis and in adulthood (Liu and Robinson-Rechavi, 2018a). Transcriptome divergence across species, however, is low throughout most of embryogenesis and highest after reproductive maturity; within embryogenesis, however, expression divergence is greatest at the earliest stages (Kalinka et al., 2010; Liu and Robinson-Rechavi, 2018a). Sterility tends to arise before inviability in hybrids of a given phylogenetic distance in Drosophila (Coyne and Orr, 1997), and hybrid misexpression among spermatogenesis genes is much more pronounced in adults than in late-stage larvae (Moehring et al., 2007). Genes with greater tissue-specificity show faster coding sequence evolution (Larracuente et al., 2008). Hybrid misexpression less often involves genes implicated in transcriptional regulation than expected (Moehring et al., 2007), and defects of chromosome condensation in mitosis, nucleoporins, and transcriptional regulators of selfish elements have been implicated in the inviability of hybrid larvae (Orr et al., 1997; Barbash et al., 2003; Tang and Presgraves, 2009; Satyaki et al., 2014). It will be interesting for future work to assess whether those embryonic stages that show greatest conservation in expression and coding sequence evolution might show disproportionate hybrid dysfunction, as hinted from experiments in Caenorhabditis. The syncytial nature of Drosophila embryos through the ~6000 cell stage, however, may lead to distinct genetic network architecture and expectations in the timing of selective pressures on gene products deposited maternally versus expressed zygotically across embryogenesis.

Bufo ontogenetic profiles of hybrid dysfunction

The developmental timing of hybrid inviability in frogs and toads enjoys much richer literature than for many other systems. Hybrid developmental data exist for several genera, including Hyla (Fouquette, 1960; Mecham, 1960; Mecham, 1965; Kuramoto, 1984; Kawamura et al., 1990), Pseudacris (Mecham, 1965), Rana (Kuramoto, 1974; Frost and Platz, 1983; Sumida et al., 2003), and most notably Bufo (Blair, 1972; Malone and Fontenot, 2008). The developmental stages of hybrid inviability are described for over 600 inter-species Bufo cross combinations, with late-stage hybrid dysfunction being most prevalent for species pairs that are less genetically divergent (Figure 2). By contrast, genome and transcriptome information is rare, due to the large genome sizes of most frogs and toads. One key exception among anurans is for the model species Xenopus laevis and X. tropicalis: comparative transcriptome analysis across their embryonic development showed changes in expression profiles of many genes despite a background of conservation in expression for most genes (Yanai et al., 2011). Interestingly, earlier stages of embryogenesis showed greater divergence in gene expression (Figure 2). Post-embryonic development was not analyzed, however, making it unclear what profile of gene expression divergence describes all of ontogeny. It will be valuable in future work to link directly the ontogeny of hybrid dysfunction to the degree of divergence in transcriptomes and sequences for genes expressed differentially across development.

Conclusions

A general understanding of genetic mechanisms in the speciation process integrates the ontogenetic timing of gene expression and gene action in developmental programs. Predictable ontogenetic trends in the molecular evolution of proteins and their regulation will introduce predictability into how genetic networks diverge and how Dobzhansky-Muller incompatibilities manifest dysfunctional phenotypes in hybrids. In this underexploited way, the genetics of post-zygotic isolation in the speciation process dovetails with research programs in evolutionary developmental genetics. The diversity of theoretical perspectives that contribute predictions to ontogenetic patterns of molecular evolution are, however, at best, incompletely integrated with conceptions about genetic architectures (pleiotropy, modularity, robustness). We still have only a rudimentary understanding of how the genetic clockwork might set the ontogenetic timing of this developmental alarm clock in the dawn of new species. Consequently, it is a challenge to extract consensus on patterns and predictions about sequence evolution and DMI incidence over the course of ontogeny.

Key features of an ontogenetic view of molecular evolution include the pleiotropic roles of genes and the pleiotropic effects of genetic perturbation, as well as the dynamism of genetic network modularity over development, to influence the evolvability of genes and the robustness of phenotypic outputs. A growing body of empirical literature is documenting the dynamics of gene expression and molecular evolution over developmental time. What is missing in the nascent state of the field is an integrated set of theoretical expectations for how these features can produce emergent trends of genome evolution and of dysfunctional genetic networks in the divergent set of genomes of hybrid individuals. Empirical tests of existing theoretical predictions, from models that focus primarily on only a subset of ontogeny, will benefit from exploring the relative influence of adaptive molecular evolution and purifying selection on genes that vary in expression over all ontogeny. Answering the neglected question -- ‘Are certain developmental processes especially likely to be disrupted in hybrids?’ -- offers the promise to identify new rules of speciation as well as rules in the molecular evolution of development.