Consideration of microarray methodology
Previous studies have shown that using single-species microarrays in order to measure the expression levels of the mRNA of multiple species can produce biased estimates due to sequence divergence between the probe and hybridized mRNA if appropriate statistical thresholds are not employed . With this caveat in mind, we sought to minimize the effect of binding bias on our analysis by: (a) employing a minimum 1.5-fold change, expression threshold for expression differences to be considered significant; and (b) analysing only genes with expression information on all available replicate spots. Gilad et al.  noted that a 1.5-fold threshold difference in the expression level provided near 100% specificity (albeit at a cost in sensitivity) to accurately measuring significant expression differences between human and orangutan mRNA samples hybridized on human microarrays, which differ in nucleotide sequence by approximately 3% (D. melanogaster and D. simulans/D. sechellia differ ~3% ). While this has led to a reduction in the size of our dataset which is available for analysis, as well as a reduced sensitivity to significant changes in expression level, the increased specificity is necessary in order to minimize the possibility of spurious inferences. We further validated the results of our expression analysis by testing for a potential correlation between pairwise expression and sequence divergence among the three pure species and found no significant correlation (Spearman rank correlation test, P > 0.05) suggesting that our methods were appropriately conservative (see Methods).
However, it should be noted that the use of such conservative methods does impose a cost on the general applicability of the overall analysis, as our final dataset involves the analysis of ~20% of the genes spotted on the D. melanogaster 12Kv2 cDNA microarray. In order to ascertain the possibility that our final dataset was not representative of the genomic fraction probed by the array, we tested whether those genes analysed showed enrichment, or a paucity, of certain gene ontology (GO) terms as compared to the overall array using the FATIGO web tool . The distribution of GO terms among the genes analysed was not significantly different from that represented on the entire array (data not shown). We also tested for the possibility that genes used in the analysis were more conserved than the array-wide average, as would be expected if the reason that a significant portion of the genes were being rejected was due to binding artifacts generated by probe-sequence divergence. However, we found a statistically significant bias in favour of greater divergence among genes used in the analysis, when compared to those that were not used for all three species (see Methods). While this may suggest that the genes employed in the analysis are, on average, more rapidly diverging than the genes that were excluded from analysis, it does not suggest that our analysis is biased towards genes that show greater sequence conservation because of artifacts resulting from probe/RNA divergence.
Despite this, a certain degree of bias may be expected in our dataset given the limitation of our analysis to genes that were detectibly expressed in all stages among the species being compared. Applying such a restriction avoided the possibility that a gene would be considered not expressed at a particular developmental time-point simply because of technical errors associated with probe/RNA binding on the microarrays. However, it also carries the negative consequence of ignoring genes that are legitimately expressed at very low levels during certain periods of ontogeny. It is unlikely that much useful information was lost during such stages as microarrays are known to have low power to detect significant differences in comparisons between samples with very low absolute expression levels . Conversely, if genes that were excluded from our analysis due to their very low expression levels during particular stages showed significant between-species/hybrid differences in expression levels during other stages, such differences would be missed in our analysis. It is not clear if the value of inclusion of such genes would offset the potential biases introduced by including false-negatives (see above). Such considerations will be resolved as newer, more sensitive, techniques are employed to answer such questions (for example, RNA-seq). In addition, these techniques will also give us the capability of inferring meaningful evolutionary change in expression level among changes below the 1.5-fold expression threshold employed in this study. Nevertheless, a note of caution is warranted in studies that seek to compare loci with minute expression differences: substantial inter-strain variation has been observed in Drosophila  and, thus, it will be important to determine whether such differences are strain-specific or species-wide.
Within-species variation in expression levels over ontogeny
Our estimates of the percentage of genes that vary significantly in expression level over the course of ontogeny (~60%-80%) are similar to those reported in a previous analysis of gene expression over the entire course of ontogeny in D. melanogaster (~86%) . The reduction in the proportion of genes varying significantly in expression level between stages in our study is most likely a result of our having sampled only the latter portion of development. The significant increase in the proportion of developmentally modulated genes in D. sechellia is striking. However, we noted a systematic bias towards low between-replicate variance in D. sechellia compared to the other samples. Correcting for this bias in D. sechellia caused much of the elevated signal of developmental modulation to disappear. Despite this, the general pattern of an elevated number of developmentally modulated genes in D. sechellia remained significant in comparison with D. simulans which, again, suggests that D. sechellia may have undergone lineage specific divergence in terms of its developmental expression profiles. It is also interesting to note that the observed reduction in between replicate variance in D. sechellia is consistent with previous population genetics studies of nucleotide diversity in this species, which have found a significantly reduced level of within species polymorphism relative to D. simulans or D. melanogaster [20, 23] (see below).
Regulatory divergence in the context of ontogeny
When comparisons between stages were statistically significant, as was the case in the comparisons between D. melanogaster and D. simulans as well as D. simulans and D. sechellia, we found that expression level is more conserved between the pure species during earlier as compared to later stages (Figure 1). Thus, Von Baer's  classic observation that earlier stages of ontogeny are more conserved than later stages may apply at both the level of nucleotide sequence [14–16] and transcript abundance (see also Additional File 1 for an analysis of nucleotide divergence in the context of ontogeny using the expression data from the present study). Expression levels are not significantly more conserved during earlier developmental stages in our comparison between D. melanogaster and D. sechellia, which may be the result of two possibilities (Figure 1). First, our analysis may not have provided the sensitivity required to detect differences in the divergence patterns of expression over development in this comparison. Our observation of a reduced level of inter-replicate variability in D. sechellia may have led to our observation of a relatively uniform signal of significant divergence in expression level over ontogeny relative to D. melanogaster. Alternatively, it is possible that D. sechellia has been subject to substantial divergence in expression levels in early developmental stages, more of which are shared with its closer relative D. simulans, as compared to D. melanogaster. It is possible that, while we are able to detect significant difference among stages in comparisons with D. simulans, similar comparisons with D. melanogaster show a more uniform distribution of divergence over development.
Von Baer's pattern is not represented in comparisons among the hybrids and either of their parental species (Additional File 5). In the case of the comparison with D. simulans there is no significant difference in the proportion of genes differentially expressed among the four sampled developmental time-points, whereas, in the case of the comparison with D. sechellia, the intermediate time points show less divergence. Such an observation may indicate that the hybrids neither represent an intermediate phenotype between their parental species, nor are they more closely allied with one over the entire course of development. Rather, the hybrids may express phenotypic traits more consonant with one parent as compared to the other, which varies over the course of ontogeny (see below). This lack of a Von Baerian trajectory in expression level is also apparent among misexpressed genes in the hybrids, where the larval stage shows the highest proportion of misexpressed genes and the two sampled pupal time points the least (note that a gene is only considered 'misexpressed' if it is either under- or over-expressed in the male hybrids relative to both male parents). While this seems to contradict the evidence for a greater conservation of expression levels during earlier stages observed in comparisons of parental species, it is possible that comparisons of expression level among pure species and pure species and their hybrids are measuring two different aspects of regulatory divergence. Whereas interspecific comparisons reveal only the ultimate outcome or regulatory divergence (expression level), comparisons between hybrids and their parents allow us to observe the result of divergence at all levels of gene regulation, from cis-regulatory elements to the feedback/forward regulation occurring in the overall hybrid regulatory network. During the larval stage, we may be observing what True and Haag  have termed 'developmental systems drift': stabilizing selection appears to be the primary evolutionary force acting upon expression levels for the majority of genes , and while the underlying regulatory machinery may have diverged between D. sechellia and D. simulans, this divergence may be compensatory such that it does not manifest itself in terms of expression level differences between the parental species  (Figure 1). Rather it reveals itself only as improper regulation in the case of the hybrids (see below). The larval stage is characterized by a high rate of growth, which is associated with a rapid increase in transcription of total mRNA and translation of proteins , which may generate selective pressure in order to maintain uniformly high expression levels despite substantial divergence among underlying regulatory machinery. Such non-adaptive divergence of complex regulatory systems is facilitated in species with reduced effective population sizes, due to the increased probability of fixation of neutral or slightly deleterious mutations , which, as noted above, is known to be the case in D. sechellia. On the other hand, the complex organogenesis occurring during metamorphosis may involve greater integration among the regulatory circuits (transcription factors or cis-regulatory elements) than other stages, leading to its underlying machinery being more refractive to divergence causing misexpression in hybrids. Nevertheless, assuming that an organism's ultimate phenotype determined by gene expression levels (and subsequent post-transcriptional regulation), the conserved expression levels among parental species during earlier developmental stages is probably a compelling validation of Von Baer's third law.
Both our clustering analysis (Figure 2), as well as our analysis of the number of genes differentially expressed in pairwise comparisons, support phylogenetic expectations of between-species divergence in the case of the larval, late pupal, and adult stages (D. sechellia and D. simulans are more similar to one another). However, this is not the case during the early pupal stage, where D. melanogaster and D. simulans cluster together and show the fewest significantly differentially expressed genes (Figure 1). This may suggest that, as part of its adaptation to its host plant (Morinda citrifolia), D. sechellia may have been exposed to unique selective pressures that have altered particular aspects of its larval or early pupal development. This hypothesis is supported by previous studies that have examined developmental phenotypes in the D. melanogaster group [38–40] and have found evidence of altered developmental phenotypes in D. sechellia relative to other species, where the opportunity for selection may have been responsible for altering aspects of its ontogeny. An alternative explanation for our observation of D. sechellia's status as an outlier is that 'typical' gene expression levels in this species would most probably be observed when it is living on its native host, M. citrifolia . By raising our species/hybrids on standard cormeal/molasses/agar medium, we may have induced stress in D. sechellia, leading to a pronounced difference in expression patterns relative to the other species. We tested to see if genes that were uniquely differentially expressed at each stage in D. sechellia were over-represented in terms of certain GO categories relative to our entire dataset, although no differences were statistically significant (data not shown). The possibility does remain, however, that some proportion of D. sechellia's unique expression patterns results from stress (though these flies have been kept on standard cornmeal medium since the mid-1980s; W. Haerty personal communication), but the present analysis of the hybrids in particular would have been complicated were D. sechellia raised on medium containing M. citrifolia while D. simulans was not.
The developmental basis of hybrid misregulation
A number of studies involving genome-scale transcriptional profiling have revealed that a substantial proportion of the transcriptome (> 10%) is misexpressed in interspecific hybrids relative to their parental species . Work in the field of speciation has been quite equivocal about the relative number of genes involved in producing the sterile phenotype observed in hybrids. However, recent evidence has supported the notion that sterility in hybrids between closely related species such as those used in the present study is primarily the result of a small number of genes of relatively large effect (≥ 6) . Thus, it appears unlikely that the misexpression observed in our study stems from the incompatible divergence of cis-regulatory factors at such a large number of loci. More plausible is the hypothesis that the large scale-patterns of misexpression observed in interspecific hybrids result from divergence of a smaller number of loci having widespread effects in trans - perhaps corresponding to some of the loci of large effect revealed through introgression studies [43, 44]. However, our study does not support the suggestion that these trans acting factors are derived from the cascading effects of a smaller number of genes that are significantly misexpressed at earlier stages of development . In contrast, it would appear that there is considerable stage-specific autonomy of regulatory breakdown, with no obvious pattern of an increasing proportion of genes misexpressed during subsequent stages. MBGs are over-represented among under-expressed, misexpressed genes during three of the sampled developmental stages, indicating that common regulatory factors or selection pressures (for example, sexual selection on MBGs in the adult) may nevertheless underlie misexpression at multiple stages.
The highest proportion of genes is significantly misexpressed during the larval stage (Figure 3; Table 1). Two non-mutually exclusive hypotheses may account for an elevated degree of misexpression during this early stage. First, as mentioned above, the larval stage is characterized by a rapid increase in transcription of total mRNA and translation of proteins . Slight heterochronic changes in hybrid development (for example, later activation of transcriptional machinery) may manifest themselves as widespread under-expression or even over-expression of larval genes. Some evidence for this hypothesis is provided by previous observations which suggested that spermatogenesis may be delayed in third instar larval hybrids between D. simulans and D. mauritiana, which could lead to qualitative differences in expression pattern of genes involved in this process . Secondly, previous studies have suggested that D. sechellia has been subject to divergence in embryonic and larval ontogeny [38, 39]. As stabilizing selection appears to be the primary evolutionary force acting upon expression levels for the majority of genes , this divergence may not manifest itself in terms of expression level differences in the parental species (Figure 1). Rather, the underlying regulatory machinery may have diverged revealing itself as an elevated proportion of incompatibilities in hybrids (see above).
If we restrict our analysis to genes expressed in the hybrids that are significantly differentially expressed between males of D. simulans and D. sechellia, they are significantly more likely to be expressed at the D. simulans level in the hybrids during the larval and adult stages, whereas they are more likely to be expressed at D. sechellia levels during the two pupal stages (Figure 3, Table 1). One may assume that expression levels would generally show an overall dominance in hybrid males in the direction of the parent from which it inherits its X chromosome (in this case D. simulans), assuming, of course, that a significant number of regulatory loci are harboured on the X. Our results suggest that this is not the case during all stages and, while the ability of regulatory factors to interact is more conserved during the pupal time points, significant divergence has occurred between the two parental species that manifests itself dominantly with regards to D. sechellia. This hypothesis is supported by our observation that the degree to which genes vary in expression level between developmental transitions in hybrids is significantly more similar to the D. sechellia parent during the EP to LP transition (Additional File 6).
A classic study, testing the vulnerability of various stages of Drosophila ontogeny to induced mortality when exposed to X-rays, found that susceptibility is highest during pupation, which suggests that this stage is particularly sensitive to deleterious perturbation . Such an observation is particularly intriguing in the context of Raff's  'developmental hourglass' hypothesis, which argues that stages involved in organogenesis may be more resistant to evolutionary divergence than preceding or subsequent stages. The original hypothesis focused on embryogenesis. However, our observation of a significantly reduced number of misexpressed genes during the two sampled pupal time points suggests that the mechanisms underlying gene regulation during this stage may be more conserved (Table 1). Holometabolous insects such as Drosophila, undergo two rounds of extensive organogenesis (for example, embryogenesis and metamorphosis), and may have two periods of increased regulatory conservation. Interestingly, however, the decreased proportion of genes significantly misexpressed in the hybrids during the pupal stage does not coincide with a decreased proportion of genes significantly differentially expressed between the two parental species (Figure 1).