Comparative phenomics: a new approach to study heterochrony

Abstract

Understanding the links between development and evolution is one of the major challenges of biology. ‘Heterochronies’, evolutionary alterations in the timings of development are posited as a key mechanism of evolutionary change, but their quantification requires gross simplification of organismal development. Consequently, how changes in event timings influence development more broadly is poorly understood. Here, we measure organismal development as spectra of energy in pixel values of video, creating high-dimensional landscapes integrating development of all visible form and function. This approach we termed ‘Energy proxy traits’ (EPTs) is applied alongside previously identified heterochronies in three freshwater pulmonate molluscs (Lymnaea stagnalis, Radix balthica and Physella acuta). EPTs were calculated from time-lapse video of embryonic development to construct a continuous functional time series. High-dimensional transitions in phenotype aligned with major sequence heterochronies between species. Furthermore, differences in event timings between conspecifics were associated with changes in high-dimensional phenotypic space. We reveal EPTs as a powerful approach to considering the evolutionary importance of alterations to developmental event timings. Reimagining the phenotype as energy spectra enabled continuous quantification of developmental changes in high-dimensional phenotypic space, rather than measurement of timings of discrete events. This approach has the possibility to transform how we study heterochrony and development more generally.

1 Introduction

Heterochronies, changes in the timings of developmental processes between ancestors and their descendants, are proposed as an important mechanism of evolutionary change, and are frequently regarded as the main process linking development to evolution (Gould, 1977; Gould, 1982; McKinney, 1988; Smith, 2001; Smith, 2003; Keyte and Smith, 2014). Heterochrony research typically involves comparing the timings of development between extant, closely related taxa, and subsequently inferring evolutionary change by mapping these changes to a phylogeny (Smith, 2003). A relatively recent advancement has been the use of relative timings of developmental events, and changes in their sequence relative to one another (Smith, 2001; Jeffery et al., 2002; Bininda-Emonds et al., 2002). Historically, such studies have focussed on the rates of growth of morphological structures, primarily because of the reliance on morphology by palaeontologists studying heterochrony, who are often forced to use size as a proxy for age (Gould, 1977; McKinney, 1988; Spicer and Gaston, 1999; Smith, 2001; Smith, 2003). More recently there have been calls for more integrated approaches to the investigation of developmental event timings, via the inclusion of both functional and behavioural developmental characters (Spicer, 2006; Spicer and Rundle, 2006). Incorporation of functional developmental event timings into heterochrony research has enabled changes in functional events to be examined within the evolutionary context of heterochrony (Spicer and Rundle, 2006).

Measuring the timings of discrete developmental events enables direct comparison of developmental itineraries between species, but this comes at the expense of reducing complex and dynamic developmental processes down to a single point (Burggren, 2021). Furthermore, the selection of developmental events from a vast number of candidates is; i) reliant on their occurrence in the study species of interest, and therefore must be identified a priori (Walls et al., 2019), and ii) potentially introduces a significant element of chance to selection of developmental events related to biological processes of interest (Houle et al., 2010; Forsman, 2015). Consequently, our understanding of the links between development and evolution via heterochrony must be limited by the current frameworks and methodologies used to investigate it. What is required is a new, more objective, way of describing and analysing heterochrony and development generally drawing on as much biological information as possible.

The capacity to use bioimaging to continuously measure phenotypic change at the organism wide scale in developing animals offers new opportunities for interrogating the lines between heterochrony as a pattern, and its role as a process in macro-evolutionary change (Burggren, 2021). Bioimaging enables researchers to apply new computer-vision approaches to measuring phenotypic change, using methods with no manual equivalents. Energy proxy traits (EPTs) are a measure of fluctuations in pixel intensities quantified as a spectrum of energies across different temporal frequencies (Welch, 1967) and are proving a valuable approach to measuring complex phenotypes in developing embryos (Tills et al., 2018; Tills et al., 2021). EPTs, rather than selecting specific aspects of an organism’s morphology, physiology or behaviour, are indiscriminate measures of the phenotype applicable to different species and experimental designs that can be followed continuously during the course of development. As a method of quantifying features of developing embryos they therefore overcome the limitations associated with measuring discrete points in development time. EPTs are effective at capturing developmental transitions in embryos of an aquatic invertebrate, where traditional phenotypic measures are largely ineffective or non-transferable between stages of development. There is also evidence to suggest that they are indicative of energy turnover at the biochemical level (Tills et al., 2021). However, so far EPTs have not been used to compare species with different developmental itineraries, to investigate high-dimensional phenotypic change associated with evolutionary differences in the timings of development.

Consequently, our aim here is to present an alternative way of describing and analysing heterochrony, by investigating the extent to which evolutionary differences in the timings of developmental events (heterochronies) are associated with high-dimensional phenotypic change, using EPTs. We hypothesise that the onset of developmental events are associated with changes in EPT spectra, and that evolutionary differences in the timings of developmental events are reflected in time series of EPT spectra. To do this, we measured interspecific differences in EPTs in encapsulated embryos of species within a well-resolved phylogeny. Species of freshwater pulmonate snail exhibit sequence heterochronies at the familial level (Smirthwaite et al., 2007), thereby providing a tractable model system. Smirthwaite et al. (2007) investigated event timings in 13 species across three families and detected significant sequence heterochronies in embryos of the Lymnaeidae and Physidae. Embryos of the physid Physella acuta exhibit sequence heterochronies in the timings of muscular crawling and cardiac function, relative to embryos of the lymnaeids Lymnaea stagnalis and Radix balthica, specifically an earlier onset of cardiac function relative to muscular crawling in the two lymnaeids. Therefore, we used EPTs to characterise high-dimensional changes in observable phenotype throughout embryonic development in these species, alongside measuring the timings of major functional developmental events: i) the onset of ciliary driven rotation; ii) the onset of cardiac function; iii) attachment to the wall of the egg capsule and the transition to muscular crawling; iv) the onset of radula function (referred to throughout as rotation, heart, crawling and radula respectively) (Smirthwaite et al., 2007).

2 Materials and methods

2.1 Animal collection and maintenance

Adult snails L. stagnalis and P. acuta were collected using a sweep net (1 mm mesh) from Exeter canal and nearby streams, Devon, United Kingdom (50°41′57.8″N 3°30′43.7″W). Adult R. balthica were collected using the same method from drainage canals, Middle Furlong Rhynie Bridgwater, UK (51°11′23.9″N 2°52′47.9″W). Snails were immediately transferred to the laboratory in plastic containers containing water and pondweed within 24 h of collection. Upon arrival snails were divided between a number of plastic containers (volume = 4 L) each filled with continuously aerated artificial pond water (CaSO₄—120 mg L⁻¹, MgSO₄—245 mg L⁻¹, NaHCO₃—192 mg L⁻¹, KCl—8 mg L⁻¹), and maintained at T = 15°C. Snails were fed spinach and lettuce ad libitium. Snails were acclimated to laboratory conditions for a minimum of 1 week prior to experimentation under a 12 h light/12 h dark light regime, and with weekly water changes.

2.2 Embryo collection

Snails regularly deposited egg masses onto the walls and floor of rearing aquaria. These masses were carefully removed using a piece of thin laminate plastic. When viewed under low power magnification (10–40×) any eggs that had not developed past the 4-cell stage were removed. Embryos from a minimum of 3 egg masses were used for each species in order to account for any brood variation. Individual embryos were carefully removed from their egg masses selected haphazardly before being transferred into individual wells of microtitre plates containing artificial pond water (Nunc, Microwell, 96 wells, 350 µL per well).

2.3 Bio-imaging

Embryonic development from the 4-cell stage to hatching was recorded using an Open Video Microscope (OpenVIM), enabling long term repeated video imaging of aquatic embryos (Tills et al., 2018). Microtitre plates containing embryos were placed into incubation chambers (H101-K-Frame, Okolab™, Italy), and reared for the duration of their embryonic development at 20°C. Temperature was controlled by circulation of water through the chamber by a temperature bath (H101-CRYO-BL, Okolab™, Italy). Water in the incubation chambers was constantly but gently aerated using an air pump (OKO AP, Okolab™, Italy). To reduce evaporation within wells of microtitre plates, air was pre-humidified using a humidity module (Okolab™, Italy). Incubation chambers were mounted onto an aluminium frame, the position of which was controlled using a motorised XY stage (SCAN 130 × 85, Märzhäuser Wetzlar™, Germany). There were mortalities of 33.3%, 16.7% and 10.4% during various points in the development of embryos of L. stagnalis, R. balthica and P. acuta respectively. Mortalities were excluded from analyses given that these embryos did not undergo all developmental events used in this study.

Image sequences of individual embryos were acquired using an inverted lens at ×200 magnification (VH-720R, Keyence™, United Kingdom) attached to a Charged Couple Device digital camera (resolution: 2048 × 2048 pixels, Pike F421B, Allied Vision™, Germany). Dark field illumination was achieved using an LED ring light placed above the incubation chamber (LDR2-42-SW2, CCS, United Kingdom). Raw video data some of which was published in Tills et al., 2018; Tills et al., 2021), was used for R. balthica. Image sequences were acquired hourly for 30 s at between 30 frames s^-1 (R. balthica) and 48 frames s^-1 (P. acuta and L. stagnalis) for the duration of embryonic development, using the open source ImageJ plugin µManager (Edelstein et al., 2010). A resolution of 512 × 512 pixels was used for embryos of P. acuta and L. stagnalis, and 1,048 × 1,048 pixels for embryos of R. balthica.

2.4 Image analysis

Image sequences of developing embryos were analysed using an open source Python package Embryo Computer Vision (EmbryoCV) (16). Energy proxy traits (EPTs) were calculated for each 30 s video timepoint, for each embryo, using the method of Tills et al. (2021). In brief, mean-pixel values of the region containing the embryo in each frame were extracted as a time series. Signal decomposition using Welch’s method (1967) was then used to decompose mean pixel value signals into the temporal domain, thereby providing a spectra of energy values at different temporal frequencies.

For each embryo, the timings of onset of a number of key developmental events were quantified from video via manual analysis for each species. These were i) the onset of ciliary driven rotation, ii) the onset of cardiac function (determined by the first visible heart beat), iii) attachment to the wall of the egg capsule and transition to muscular crawling (the point at which there is a clear attachment of the foot, rather than the embryo simply ‘resting’ on the wall of the egg capsule), and iv) the onset of radula function, the last developmental stage in all three species used (Smirthwaite et al., 2007) (Table 1; Figure 1). Heterochronies in the timing of muscular crawling and cardiac function exist between P. acuta, L. stagnalis and R. balthica. Physella acuta exhibits a significantly earlier onset of attachment and crawling on the wall of the egg capsule, relative to the onset of cardiac function, and relative to the timing of this event in L. stagnalis and R. balthica (Smirthwaite et al., 2007).

FIGURE 1

To visualise developmental differences in EPTs between species, the sum of energy levels across all frequencies of the EPT spectra (hereafter referred to as total energy) were calculated, normalised within each individual (0–1), and expressed as a time series. The duration of embryonic development (i.e., the time taken to hatch) varies between individuals and species, therefore to generate a standardised developmental rate that would allow for comparison between species and individuals, and be invariant to differences in overall rates of development, the absolute timings from the 4-cell stage to hatching were converted to relative time (0–1).

2.5 Dimensionality reduction

EPT spectra are high-dimensional representations of observable movement, and therefore to compare interspecific, and development stage specific differences in combinatorial signals across frequency spectra, dimensionality reduction in the form of principal component analysis (PCA) was used. Temporal frequency data were binned to 0.1 Hz increments (0.03—6.0 Hz, 60 frequency bins). Frequency data were restricted to 6.0 Hz as the upper limit at which biologically meaningful signals would be expected. Mean energy within each frequency bin for each normalised time point (0–1) was calculated and log transformed (McCoy et al., 2023). PCA (prcomp, package ‘stats’, v4.0.3) was applied to these data and the resulting eigenvectors were used to determine whether the onset of developmental events were associated with changes in high dimensional phenotypic space.

2.6 Statistical analysis

All data were analysed in R v4.0.3 (R Core Team, 2020). To determine whether time series of total energy data were different between species, a repeated measures analysis of variance (ANOVA, p < 0.05) was applied to total energy data. Additionally, Bayesian structural time series (BSTS) model using the R package CausalImpact (Brodersen et al., 2014) was applied to determine whether the onset of developmental events were associated with changes in total energy time series. To determine whether the onset of developmental events was associated with frequency specific differences in EPTs, a Kruskal–Wallis test (p < 0.05) was applied to mean energy values within each of the 60 frequency bands of 5 time point values before and after the onset of each developmental event. Multiple testing correction was applied using the Bonferroni method.

3 Results

3.1 Interspecific differences in developmental EPTs

Before we discuss the differences in EPT profiles between species, we will first describe the general trends in EPT profiles observed in these embryos. First, the timing of onset of rotation was associated with pronounced increases in total energy in the first quartile of development (Table 2) in each species, and reach a peak in total energy (Figure 2), followed by reductions in total energy prior to the onset of the heart and crawling approximately midway through development (Figure 2; Table 2). Total energy rapidly declines prior to the onset of crawling, before gradually reducing in the final quartile of development in all species (Figure 2).

FIGURE 2

Clear differences occurred in EPT profiles between species (repeated measures ANOVA, F_{2, 199} = 54.52, p < 0.0001, Figure 2). Additionally to this, the onset of functional developmental events were also associated with pronounced differences in EPTs within each species (Table 2). The reduction in total energy levels prior to the onset of crawling appeared to occur earlier in Physella acuta relative to Lymnaea stagnalis and Radix balthica. Additionally, this period of total energy decline preceding the onset of crawling was compressed in P. acuta relative to L. stagnalis and R. balthica. In L. stagnalis and R. balthica, embryos remain free swimming during the development of the heart and other functions. Before crawling, embryos will intermittently rest with increasing frequency on the wall of the egg capsule, which is evident via the gradual reduction in total energy levels before crawling. However, in P. acuta this free swimming stage is absent which likely explains the comparatively compressed period of total energy decline prior to crawling. Finally, the generally later onset of each developmental event in R. balthica was associated with a shift of the entire time series later into relative developmental time (Figure 2).

3.2 Dimensionality reduction for EPT spectra differentiation

To determine whether the onset of developmental events was associated with changes in EPTs, PCA was applied to multivariate EPT time series. Data collected from individuals for each species were first averaged (mean) by normalised time point (0–1). Following this, levels of energy within each frequency band were normalised (0–1) across relative developmental time.

Principle components 1 and 2 incorporated 97.2%, 95.2% and 93.5% of the variance in the EPT spectrum for L. stagnalis (PC1: 85.5%, PC2: 11.7%) (Figure 3A), R. balthica (PC1: 75.2%, PC2: 20%) (Figure 3B), and P. acuta (PC1: 69.5%, PC2: 24%) (Figure 3C) respectively. The onset of some developmental events were associated with distinct separation of points, notably the onset of muscular crawling and radula function in L. stagnalis and R. balthica. In L. stagnalis we observed clear separation of points along the axis of PC1 between the onset of cardiac function and muscular crawling. Variation along the axis of PC1 was predominantly driven by frequencies ranging from 1.7 to 3.7 Hz (Figure 3A). Additionally, frequencies up to 0.9 Hz were predominantly driving variation along the axis of PC2, suggesting that differences in embryos between the onset of crawling and radula function were driven predominantly by frequencies within this range (Figure 3A) (Supplementary Information S2).

FIGURE 3

3.3 High-dimensional phenotypic change associated with the onset of developmental events

EPT are a spectra of energy values across different temporal frequencies. Therefore, to determine whether the onset of the developmental events used in this study were associated with composition differences of EPT spectra across temporal frequency bands, pairwise comparisons of energy before and after the onset of each developmental event were carried out for binned frequencies. In all species, the onset of ciliary driven rotation was associated with a significant increase in the levels of energy in all temporal frequency bands (Kruskal–Wallis, p < 0.0001, Figure 4; Supplementary Information S3). The onset of cardiac function was associated with a significant increase in energy up to 0.7 Hz and 1.6–1.8 Hz in L. stagnalis. In R. balthica and P. acuta the onset of cardiac function coincided with a significant increase in energy levels within the range of 1.2–1.6 Hz and 1.6–4.0 Hz respectively. (Kruskal–Wallis, p < 0.0001, Figure 4). The onset of muscular crawling on the wall of the egg capsule was associated with a significant reduction in energy within all temporal frequency bands for P. acuta (Kruskal–Wallis, p < 0.0001, Figure 4). Finally, following the onset of radula function, there were significant reductions in the levels of energy up to 0.5 Hz in P. acuta (Kruskal–Wallis, p < 0.0001, Figure 4).

FIGURE 4

4 Discussion

Our aim was to present an alternative way of describing and analysing heterochrony, by determining the extent to which evolutionary differences in the timings of developmental events (heterochronies) were associated with high-dimensional phenotypic change, using energy proxy traits (EPTs). We hypothesised that the onset of developmental events used in this study are associated with changes in EPT spectra, and that evolutionary differences in the timings of developmental events would be associated with changes in time series of EPT spectra between species. Causal impact analysis and PCA revealed that the onset of developmental events were associated with changes in EPT spectra. Furthermore, from our results we conclude that evolutionary differences in the timings of developmental events are associated with pronounced changes in high-dimensional phenotypic space. EPTs measure complex phenotypes as a spectrum of energy and application of these to the complete embryonic development of three closely related snail species Lymnaea stagnalis, Radix balthica and Physella acuta revealed interspecific differences in these spectra, tightly associated with previously described sequence heterochronies. Furthermore, investigation of EPT time series revealed functional changes before and after the onset of developmental events, but also striking differences between species for the same developmental events. Combinatorial analysis of EPTs using dimensionality reduction revealed distinct separation of points between developmental events, highlighting transitions in the functional phenotype of embryos in high-dimensional space. This suggests that developmental events, while a useful approach to comparative studies of development, may bely complex differences in functional phenotypes that could perhaps themselves be the object of natural selection.

4.1 Interspecific differences in EPTs coincide with pre-established sequence heterochronies

Time series of EPTs reflected interspecific differences in developmental event timings. Firstly, differences in the relative timing of onset of rotation was evident from time series of total energy between species. Ciliary driven rotation occurred earliest in P. acuta and this is evident as an earlier increase in total energy in relative developmental time. Following the onset of ciliary driven rotation, levels of total energy gradually increased in all species. Following this, embryos of all species showed clear reductions in total energy leading up to the transition to muscular crawling on the wall of the egg capsule, the timing of which reflected developmental differences between these species. In P. acuta, the onset of muscular crawling is brought forward significantly relative to the timing of cardiac function, and the onset of this event in L. stagnalis and R. balthica (Smirthwaite et al., 2007; Supplementary Information S1). From time series of EPT data, we can see a marked reduction in total energy occurring at approximately 37% of relative developmental time in P. acuta whereas this pronounced reduction in energy occurs at approximately 41% and 52% of relative developmental time in L. stagnalis and R. balthica respectively (Figure 2).

The observed changes in energy both before and after the timings of discrete developmental events used in this study showcase potential shortcomings in current approaches used to quantify the timings of development. Burggren (2021) observed that current research often focusses on development as a series of discrete events, when development should be seen as a continuum. Experimental protocols should gather data at multiple points of the developmental continuum, thereby allowing phenotypic measurements to be put into the context of the entirety of development. Here, we demonstrated considerable phenotypic change both before and after the onset of discrete developmental events, including the onset of ciliary driven rotation, and attachment to the wall of the egg capsule. The gradual increases and decreases in total energy likely reflect the continuous nature with which major multi-faceted developmental transitions take place in these embryos.

Firstly, the gradual increases in total energy with time co-occurs with increases in rates of rotation as development progresses (Figure 2). Without generalising to the other two species tested, previous research showed that embryos of L. stagnalis exhibited greater rates of rotation at the hippo stage (equivalent to approx. 40% developmental time) relative to the veliger stage (approx. 30% developmental time), which may contribute to explaining the gradual increases in total energy following the onset of rotation. Byrne et al. (2009) found that in embryos of L. stagnalis rates of rotation at the hippo stage were approximately 1.3–1.6 times greater than that of the veliger stage, and suggested that such differences in rotation rate between these two stages may be due to increases in anatomical complexity at the hippo stage.

Secondly, we also observed gradual decreases in total energy prior to attachment to the wall of the egg capsule and the onset of muscular crawling. Before this developmental event, embryos of R. balthica and L. stagnalis remain free swimming in the egg capsule (Meshcheryakov, 1990; Smirthwaite et al., 2007). During this period, embryos will increase the frequency at which they intermittently ‘rest’ their head on the wall of the egg capsule, before firmly attaching with the foot and commencing muscular crawling, which likely explains the gradual reductions in energy (and therefore to some extent overall rates of movement) proceeding the onset of this developmental event. When viewed as a discrete developmental event, the considerable phenotypic changes the developing embryo undergoes prior to its onset becomes masked, and arguable the event itself becomes increasingly arbitrary.

Additionally to this, there were features of both total energy time series and time series PCA that cannot be explained by the onset of, or changes in the timings of developmental events used in this study. For example, in each of the total energy time series, there is a distinct peak in total energy occurring at approximately 35%, 50% and 33% relative developmental time for L. stagnalis, R. balthica and P. acuta respectively (Figure 2). This peak occurred prior to the onset of cardiac function in L. stagnalis and R. balthica and muscular crawling in P. acuta, yet did not coincide with the onset of any of the developmental events used in this study. Embryos of L. stagnalis and R. balthica exhibit an intermittent resting behaviour prior to crawling. This peak represents the onset of this intermittent resting behaviour, prior to which rates of rotation during free swimming are greatest (Byrne et al., 2009). This indicates that the transitions to muscular crawling may better be regarded as a gradual reduction in rates of movement prior to attachment of the foot, rather than as a discrete event in developmental time, and highlights a potential limitation of using the timings of discrete developmental events for assessing evolutionary changes in early development.

Furthermore, measurement of this change in energy as a continuum also reveals evolutionary differences in developmental phenotype prior to the onset of attachment and muscular crawling, rather than just a difference in the timing as a discrete event in developmental time. The period of total energy decline in P. acuta prior to the onset of muscular crawling is compressed relative to L. stagnalis and R. balthica (Figure 2). In P. acuta, there is an absence of a free swimming stage, with embryos attaching to the egg capsule and commencing muscular crawling following the trochophore stage (Smirthwaite, 2007), thereby producing the relatively condensed period of total energy decline observed in the total energy time series (Figure 2). Consequently, by viewing development as a continuum of phenotypic change as we have done here, we were able to identify evolutionary differences in phenotype at the whole organism scale before, during and after the onset of discrete developmental events.

4.2 A high-dimensional phenotyping approach to physiological heterochrony

Development is an inherently complex phenomenon characterised by massive functional and spatial variation, variation that may act as a source of evolutionary change and innovation (West-Eberhard, 2003). There has been a historical emphasis on heterochrony as the main mechanism linking development to evolution (Gould, 1982; Smith, 2001; Smith, 2003). Whilst measurement of the timings of developmental events or stages in investigations of heterochrony facilitates direct interspecific comparison (Walls et al., 2019) quantification of the timings of developmental events which are in themselves discrete points in time, fails to capture any notion of development as a continuous and dynamic process. Furthermore, quantifying phenotypic change at the scale of the whole organism alongside the timings of discrete developmental events may provide a means of investigating the evolutionary significance of variation in developmental event timings, by linking the phenotype, which is the ultimate object of selection, to the timings of developmental events. Most biologists currently confine phenotyping efforts to a small number of observable traits, given the often overwhelming complexity of organismal biology and the challenge of quantifying it in a discrete and reproducible way. However, selection typically does not act on single traits, rather on multiple traits simultaneously (Lande and Arnold, 1983; Phillips and Arnold, 1999). When presented with high-dimensional datasets, dimensionality reduction allows for the visualisation of structure and cumulative drivers in high-dimensional phenotypic space. In the current study, we showed considerable separation of points associated with different points in relative developmental time from principal component analysis (PCA) (Figure 3). Additionally, PCA of time series data revealed distinct clustering of points following the onset of various developmental events, but also revealed considerable variation between these developmental events (Figure 3). Results from PCA of time series analysis further demonstrates how developmental events are only a snapshot in developmental time, and that focussing on these discrete events occludes phenotypic change during the periods of development surrounding them. Given what we know about how selection operates on multiple traits simultaneously, expanding the scale at which we capture phenotypic information during periods when phenotypic complexity is at its greatest may provide greater insight into the developmental mechanisms driving evolutionary change (Blows, 2007).

Despite this, application of high-dimensional phenotyping approaches are rarely extended to species and systems that may provide effective models for investigating the evolutionary significance of changes in the timings of development. In areas of research where high-dimensional phenotyping approaches are most developed, e.g., in plant and medical phenomics, there are already well-established and standardised phenotyping approaches applicable to the model species of interest (Houle et al., 2010; Furbank and Tester, 2011; Alexandrov et al., 2016; Tardeiu et al., 2017). A potential shortcoming in the field of comparative phenomics, and across evolutionary developmental biology more generally is an absence of phenotyping approaches that are readily transferable to non-model species of interest (Tills et al., 2018). One of the cornerstones of comparative developmental physiology is the selection of species that are best suited to answer a particular biological question (Krogh, 1929), yet approaches to high-dimensional organismal phenotyping remain constrained to model animals of interest, particularly the zebrafish Danio rerio (Pelkowski et al., 2011; Kalueff et al., 2016; Peravali et al., 2018), nematode worm Caenorhabditis elegans (Yemini et al., 2013; Cornaglia et al., 2015; Olmedo et al., 2015), and the fruit fly Drosophila melanogaster (Chung et al., 2010; Levario et al., 2016). Here, EPTs were effective at characterising high-dimensional functional change in embryos of three species of freshwater gastropod, despite significant differences in their patterns of development. We suggest EPTs provide an effective and transferable phenotyping approach to quantifying phenotypic change in early life stages of a range of non-model species.

We suggest that evolutionary differences in the timings of developmental events are associated with differences in high-dimensional phenotypic space. As well as detecting interspecific differences in the timings of development between embryos of each species, we also found differences in the timings of EPTs between conspecifics. Previous research has shown that considerable standing variation exists in developmental event timings between conspecifics, which may provide the raw variation on which heterochronies are formed (Reilly et al., 1997; Mabee et al., 2000; Schmidt et al., 2004; de Jong et al., 2009; Kawajiri et al., 2009; Rundle et al., 2011; Tills et al., 2011; Rager et al., 2014). In addition to the observed interspecific differences in EPTs, individual level EPT data indicate considerable variation in EPTs between conspecifics. Comparison of developmental event timings alongside EPT data indicates that this variation corresponds with intraspecific differences in the timings of developmental events (Supplementary Informations S4-S6). The consequences of such variation are not currently understood but provide interesting lines of enquiry with which to investigate the performance and fitness implications of variation in the timings of developmental event between conspecifics. From previous research, EPTs appear to be related to some components of organismal fitness. Tills et al. (2021) showed that higher levels of total energy (reported in Figure 2) were associated with a faster growth rate in R. balthica, and posited that EPTs may in themselves be indicative of biochemical energy turnover, given that rates of biochemical energy turnover (metabolism) have been positively correlated with growth rates in a number of studies (Metcalfe et al., 1998; Nylin et al., 1998; Yamamoto et al., 1998) (although this is not always the case as such relationships are often context dependent, particularly between wild and laboratory reared specimens: Álvarez et al., 2005; Burton et al., 2011). Consequently, differences in the timings of developmental events may ultimately manifest as differences in biochemical energy turnover in developing embryos. Rates of energy expenditure are a significant object of selection (Bartheld et al., 2015), however, our understanding of the extent to which EPTs are linked with biochemical energetic turnover is currently poorly understood, and so is our understanding of how variation in EPTs ultimately influence components of organismal performance and fitness post-hatch. Further research is now needed to understand: i) the extent to which EPTs are related to biochemical energy turnover in developing embryos; and ii) the consequences of ontogenetic variation in EPTs for aspects of organismal performance and fitness.

5 Summary and conclusion

Evolutionary biologists have long sought to establish mechanistic links between development and evolution. Current research frames heterochrony, alterations in the timings of development as the main mechanism by which development leads to evolutionary change. However the current focus of heterochrony as the timings of discrete points in development may hinder quantification of phenotypic change associated with differences in the timings of these events. Here, through the application of a novel spectral phenotyping approach (EPTs), we captured a continuous functional time series of the embryonic development of three freshwater snails. Analysis of these time series provided evidence that evolutionary differences in the timings of development are associated with high-dimensional phenotypic change. Additionally to this, we provide preliminary evidence that intraspecific differences in EPTs coincide with differences in developmental event timings between conspecifics. We suggest that EPTs may provide an alternative approach to investigating the evolutionary significance of variation in the timings of development by allowing for the continuous quantification of phenotypic change at the scale of whole organism, associated with intraspecific and evolutionary differences in the timings of development. This study has the possibility of transforming the way we study heterochrony and development more generally.

References

• 1

  AlexandrovV.BrunnerD.MenalledL. B.KudwaA.Watson-JohnsonJ.MazzellaM.et al (2016). Large-scale phenome analysis defines a behavioral signature for Huntington’s disease genotype in mice. Nat. Biotechnol.34, 838–844. 10.1038/nbt.3587

  • CrossRef
  • Google Scholar
  • View reference in article

• 2

  ÁlvarezD.NiciezaA. G. (2005). Is metabolic rate a reliable predictor of growth and survival of brown trout (Salmo trutta) in the wild?Can. J. Fish. Aquatic Sci.62, 643–649. 10.1139/f04-223

  • CrossRef
  • Google Scholar
  • View reference in article

• 3

  BartheldJ. L.Gaitán-EspitiaJ. D.ArtachoP.Salgado-LuarteC.GianoliE.NespoloR. F. (2015). Energy expenditure and body size are targets of natural selection across a wide geographic range, in a terrestrial invertebrate. Funct. Ecol.29, 1463–1474. 10.1111/1365-2435.12451

  • CrossRef
  • Google Scholar
  • View reference in article

• 4

  Bininda-EmondsO. R. P.JefferyJ. E.CoatesM. I.RichardsonM. K. (2002). From Haeckel to event pairing: the evolution of development sequence. Theory Biosci.121, 297–320. 10.1111/j.1525-142X.2007.00143.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 5

  BlowsM. W. (2007). A tale of two matrices: multivariate approaches in evolutionary biology. J. Evol. Biol.20, 1–8. 10.1111/j.1420-9101.2006.01164.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 6

  BrodersenK. H.GallusserF.KoehlerJ.RemyN.ScottS. L. (2014). Inferring causal impact using Bayesian structural time-series models. Ann. Appl. Statistics9, 247–274. 10.1214/14-aoas788

  • CrossRef
  • Google Scholar

• 7

  BurggrenW. W. (2021). Developmental physiology: grand challenges. Front. Physiology12, 706061. 10.3389/fphys.2021.706061

  • CrossRef
  • Google Scholar

• 8

  BurtonT.KillenS. S.ArmstrongJ. D.MetcalfeN. B. (2011). What causes intraspecific variation in resting metabolic rate and what are its ecological consequences?Proc. R. Soc. B Biol. Sci.278, 3465–3473. 10.1098/rspb.2011.1778

  • CrossRef
  • Google Scholar
  • View reference in article

• 9

  ByrneR. A.RundleS. D.SmirthwaiteJ. J.SpicerJ. I. (2009). Embryonic rotational behaviour in the pond snail Lymnaea stagnalis: influences of environmental oxygen and development stage. Zoology112, 471–477. 10.1016/j.zool.2009.03.001

  • CrossRef
  • Google Scholar
  • View reference in article

• 10

  ChungK.KimY.GongE.ShvartsmanS.LuH. (2010). “High-throughput end-on imaging of drosophila embryo for quantitative analysis of morphogens and signalling,” in Proceedings of the 14th International conference on Miniaturized Systems for Chemistry and Life Sciences, Groningen, The Netherlands, October 2010, 3–7.

  • Google Scholar
  • View reference in article

• 11

  CornagliaM.MouchiroudL.MaretteA.NarasimhanS.LehnertT.JovaisaiteV.et al (2015). “Multi-dimensional imaging and phenotyping of C. elegans embryos via an automated microfluidic device,” in Proceedings of the MicroTAS 2015 - 19th International Conference on Miniaturized Systems for Chemistry and Life Sciences, Gyeongju, Korea, October 2015, 25–29.

  • Google Scholar

• 12

  de JongI. M. L.ColbertM. W.WitteF.RichardsonM. K. (2009). Polymorphism in developmental timing: intraspecific heterochrony in a Lake Victoria cichlid. Evol. Dev.11, 625–635. 10.1111/j.1525-142X.2009.00370.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 13

  EdelsteinA. D.AmodajN.HooverK.ValeR.StuurmanN. (2010). Computer control of microscopes using μmanager. Curr. Protoc. Mol. Biol.14.20.1–14.20.17. 10.1002/0471142727.mb1420s92

  • CrossRef
  • Google Scholar
  • View reference in article

• 14

  ForsmanA. (2015). Rethinking phenotypic plasticity and its consequences for individuals, populations and species. Heredity115, 267–284. 10.1038/hdy.2014.92

  • CrossRef
  • Google Scholar
  • View reference in article

• 15

  FurbankR. T.TesterM. (2011). Phenomics - technologies to relieve the phenotyping bottleneck. Trends Plant Sci.16, 635–644. 10.1016/j.tplants.2011.09.005

  • CrossRef
  • Google Scholar
  • View reference in article

• 16

  GouldS. J. (1982). “Change in developmental timing as a mechanism of macroevolution,” in Evolution and development. Editor BonnierJ. T. (Berlin: Springer-Verlag), 333–346.

  • Google Scholar

• 17

  GouldS. J. (1977). Ontogeny and phylogeny. Cambridge, MA: Harvard University Press.

  • Google Scholar

• 18

  HouleD.GovindarajuD. R.OmholtS. (2010). Phenomics: the next challenge. Nat. Rev. Genet.11, 855–866. 10.1038/nrg2897

  • CrossRef
  • Google Scholar
  • View reference in article

• 19

  JefferyJ. E.RichardsonM. K.CoatesM. I.Bininda-EmondsO. R. (2002). Analyzing developmental sequences within a phylogenetic framework. Syst. Biol.51, 478–491. 10.1080/10635150290069904

  • CrossRef
  • Google Scholar
  • View reference in article

• 20

  KalueffA. V.EchevarriaD. J.HomechaudhuriS.StewartA. M.CollierA. D.KaluyevaA. A.et alInternational Zebrafish Neuroscience Research Consortium (ZNRC) (2016). Zebrafish neurobehavioral phenomics for aquatic neuropharmacology and toxicology research. Aquat. Toxicol.170, 297–309. 10.1016/j.aquatox.2015.08.007

  • CrossRef
  • Google Scholar
  • View reference in article

• 21

  KawajiriM.KokitaT.YamahiraK. (2009). Heterochronic differences in fin development between latitudinal populations of the medaka Oryzias latipes (Actinopterygii: adrianichthyidae). Biol. J. Linn. Soc.97, 571–580. 10.1111/j.1095-8312.2009.01233.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 22

  KeyteA. L.SmithK. K. (2014). Heterochrony and developmental timing mechanisms: changing ontogenies in evolution. Seminars Cell Dev. Biol.34, 99–107. 10.1016/j.semcdb.2014.06.015

  • CrossRef
  • Google Scholar
  • View reference in article

• 23

  KroghA. (1929). The progress of physiology. Am. J. Physiology90, 243–251. 10.1152/ajplegacy.1929.90.2.243

  • CrossRef
  • Google Scholar
  • View reference in article

• 24

  LandeR.ArnoldS. J. (1983). The measurement of selection on correlated characters. Evolution37, 1210–1226. 10.1111/j.1558-5646.1983.tb00236.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 25

  LevarioT. J.ZhaoC.RouseT.ShvartsmanS. Y.LuH. (2016). An integrated platform for large-scale data collection and precise perturbation of live Drosophila embryos. Sci. Rep.6, 21366. 10.1038/srep21366

  • CrossRef
  • Google Scholar
  • View reference in article

• 26

  MabeeP. M.OlmsteadK. L.CubbageC. C. (2000). An experimental study of intraspecific variation, developmental timing and heterochrony in fishes. Evolution54, 2091–2106. 10.1111/j.0014-3820.2000.tb01252.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 27

  McCoyJ. C. S.SpicerJ. I.RundleS. D.TillsO. (2023). A phenomics approach reveals interspecific differences in integrated developmental responses to chronic elevated temperatures. J. Exp. Biol.226, jeb245612. 10.1242/jeb.245612

  • CrossRef
  • Google Scholar
  • View reference in article

• 28

  McKinneyM. L. (1988). Heterochrony in evolution. New York: Plenum Press.

  • Google Scholar
  • View reference in article

• 29

  MeshcheryakovV. N. (1990). “The common pond snail Lymnaea stagnalis,” in Animal species for developmental studies. Editors DettlaffT. A.VassetzkyS. G. (New York: Consultants Bureau), 69–132.

  • Google Scholar
  • View reference in article

• 30

  MetcalfeN. B. (1998). The interaction between behaviour and physiology in determining life history patterns in Atlantic salmon (Salmo salar). Can. J. Fish. Aquatic Sci.55, 93–103. 10.1139/d98-005

  • CrossRef
  • Google Scholar
  • View reference in article

• 31

  NylinS.GotthardK. (1998). Plasticity in life-history traits. Annu. Rev. Entomology43, 63–83. 10.1146/annurev.ento.43.1.63

  • CrossRef
  • Google Scholar
  • View reference in article

• 32

  OlmedoM.GeibelM.Artal-SanzM.MerrowM. (2015). A high-throughput method for the analysis of larval developmental phenotypes in Caenorhabditis elegans. Genetics201, 443–448. 10.1534/genetics.115.179242

  • CrossRef
  • Google Scholar
  • View reference in article

• 33

  PelkowskiS. D.KapoorM.RichendrferH. A.WangX.ColwillR. M.CretonR. (2011). A novel high-throughput imaging system for automated analyses of avoidance behavior in zebrafish larvae. Behav. Brain Res.223, 135–144. 10.1016/j.bbr.2011.04.033

  • CrossRef
  • Google Scholar
  • View reference in article

• 34

  PeravaliR.GehrigJ.GiselbrechtS.LütjohannD. S.HadzhievY.MüllerF.et al (2018). Automated feature detection and imaging for high-resolution screening of zebrafish embryos. BioTechniques50, 319–324. 10.2144/000113669

  • CrossRef
  • Google Scholar
  • View reference in article

• 35

  PhillipsP. C.ArnoldS. J. (1999). Hierarchical comparison of genetic variance-covariance matrices. I. using the flury hierarchy. Evolution53, 1506–1515. 10.1111/j.1558-5646.1999.tb05414.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 36

  R Core Team (2020). R: a language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing. Available at: https://www.R-project.org/.

  • Google Scholar
  • View reference in article

• 37

  RagerL.HautierL.ForasiepiA.GoswamiA.Sánchez-VillagraM. R. (2014). Timing of cranial suture closure in placental mammals: phylogenetic patterns, intraspecific variation, and comparison with marsupials. J. Morphol.275, 125–140. 10.1002/jmor.20203

  • CrossRef
  • Google Scholar
  • View reference in article

• 38

  ReillyS. M.WileyE. O.MeinhardtD. J. (1997). An integrative approach to heterochrony: the distinction between interspecific and intraspecific phenomena. Biol. J. Linn. Soc.60, 119–143. 10.1111/j.1095-8312.1997.tb01487.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 39

  RundleS. D.SmirthwaiteJ. J.ColbertM. W.SpicerJ. I. (2011). Predator cues alter the timing of developmental events in gastropod embryos. Biol. Lett.7, 285–287. 10.1098/rsbl.2010.0658

  • CrossRef
  • Google Scholar
  • View reference in article

• 40

  SchmidtK.StarckJ. M. (2004). Developmental variability during early embryonic development of zebra fish, Danio rerio. J. Exp. Zoology302, 446–457. 10.1002/jez.b.21010

  • CrossRef
  • Google Scholar
  • View reference in article

• 41

  SmirthwaiteJ. J. (2007). Developmental sequences in freshwater gastropods. Ph.D. thesis. Plymouth, England: University of Plymouth.

  • Google Scholar
  • View reference in article

• 42

  SmirthwaiteJ. J.RundleS. D.Bininda‐EmondsO. R. P.SpicerJ. I. (2007). An integrative approach identifies developmental sequence heterochronies in freshwater basommatophoran snails. Evol. Dev.9, 122–130. 10.1111/j.1525-142X.2007.00143.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 43

  SmithK. K. (2001). Heterochrony revisited: the evolution of developmental sequences. Biol. J. Linn. Soc.73, 169–186. 10.1111/j.1095-8312.2001.tb01355.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 44

  SmithK. K. (2003). Time’s arrow: heterochrony and the evolution of development. Int. J. Dev. Biol.47, 613–621.

  • Google Scholar
  • View reference in article

• 45

  SpicerJ. I. (2006). “A physiological approach to heterochrony,” in Comparative developmental physiology: contributions, tools, and trends. Editors WarburtonS. J.BurggrenW. W.PelsterB.ReiberC. L.SpicerJ. I. (Oxford: Oxford University Press), 191–202.

  • Google Scholar
  • View reference in article

• 46

  SpicerJ. I.GastonK. J. (1999). Physiological diversity and its ecological implications. Oxford, United Kingdom: Blackwell Science Limited.

  • Google Scholar
  • View reference in article

• 47

  SpicerJ. I.RundleS. D. (2006). Out of place and out of time - towards a more integrated approach to heterochrony. Anim. Biol.56, 487–502. 10.1163/157075606778967810

  • CrossRef
  • Google Scholar
  • View reference in article

• 48

  TardieuF.Cabrera-BosquetL.PridmoreT.BennettM. (2017). Plant phenomics, from sensors to knowledge. Curr. Biol.27, R770–R783. 10.1016/j.cub.2017.05.055

  • CrossRef
  • Google Scholar
  • View reference in article

• 49

  TillsO.RundleS. D.SalingerM.HaunT.PfenningerM.SpicerJ. I. (2011). A genetic basis for intraspecific differences in developmental timing?Evol. Dev.13, 542–548. 10.1111/j.1525-142X.2011.00510.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 50

  TillsO.SpicerJ. I.GrimmerA.MariniS.JieV. W.TullyE.et al (2018). A high-throughput and open-source platform for embryo phenomics. PLoS Biol.16, e3000074. 10.1371/journal.pbio.3000074

  • CrossRef
  • Google Scholar
  • View reference in article

• 51

  TillsO.SpicerJ. I.IbbiniZ.RundleS. D. (2021). Spectral phenotyping of embryonic development reveals integrative thermodynamic responses. BMC Bioinforma.22, 232. 10.1186/s12859-021-04152-1

  • CrossRef
  • Google Scholar
  • View reference in article

• 52

  WallsR. L.CooperL.ElserJ.GandolfoM. A.MungallC. J.SmithB.et al (2019). The plant ontology facilitates comparisons of plant development stages across species. Front. Plant Sci.10, 631. 10.3389/fpls.2019.00631

  • CrossRef
  • Google Scholar
  • View reference in article

• 53

  WelchP. (1967). The use of fast Fourier transform for the estimation of power spectra: a method based on time averaging over short, modified periodograms. IEEE Trans. Audio Electroacoustics15, 70–73. 10.1109/tau.1967.1161901

  • CrossRef
  • Google Scholar
  • View reference in article

• 54

  West-EberhardM. J. (2003). Developmental plasticity and evolution. Oxford: Oxford University Press.

  • Google Scholar
  • View reference in article

• 55

  YamamotoT.UedaH.HigashiS. (1998). Correlation among dominance status, metabolic rate and otolith size in masu salmon. J. Fish Biol.52, 281–290. 10.1111/j.1095-8649.1998.tb00799.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 56

  YeminiE.JucikasT.GrundyL. J.BrownA. E. X.SchaferW. R. (2013). A database of Caenorhabditis elegans behavioral phenotypes. Nat. Methods10, 877–879. 10.1038/nmeth.2560

  • CrossRef
  • Google Scholar
  • View reference in article