Mathematical modeling of the interaction between yolk utilization and fish growth in zebrafish, Danio rerio

Optimal embryonic development plays a major role in the health and longevity of an individual past the developmental stage (Wyness et al., 2013). It has long been hypothesized and since supported that nutritional and toxicological perturbation during development is associated with cardiovascular, metabolic and chronic disease later in life, known as the ‘developmental origins of health and disease’ (DOHaD) (Barker, 2004; McMillen et al., 2008). Additionally, the association between low birthweight and increased rates of coronary heart disease have been shown in various extensively replicated studies (Barker, 2004). However, most human studies rely solely on maternal nutrition as a proxy to estimate embryonic and fetal nutrition. Furthermore, the majority of these studies have been cross-sectional and do not look for specific critical windows of nutritional susceptibility during pregnancy. With both nutritional uptake and overall growth being risk factors for eventual health, it is necessary to understand not only the behavior of the processes during development but also their interactions.

Human embryonic nutrition is primarily histiotrophic until placentation during the fetal period. There is an initial maternal contribution of yolk granules in the matured oocyte, and then additional nutrition is acquired via uptake of nutrients into the yolk sac by processes such as receptor-mediated endocytosis (Burton et al., 2001; Burton et al., 2002). This histiotrophic nutrition is vital until the onset of maternal intraplacental circulation near the end of the first trimester. Because the first trimester is a dynamic period (encompassing implantation, gastrulation, embryogenesis and organogenesis), a reliable nutrient supply is vital to support these energy-intensive processes. However, our understanding of the dynamics of embryonic nutrition and their temporal impacts on growth and DOHaD is limited due to the inaccessibility of the intrauterine environment. Studies have had to rely on using maternal nutrition as a proxy for embryonic nutrition, but the appropriateness of this strategy is unknown. There is a need to understand the kinematics of embryonic nutrition, and their susceptibility to external cues such as maternal nutrition, environmental exposures and stress.

Zebrafish (Danio rerio) are an excellent model for human development and embryonic nutrition (Link and Megason, 2008). They share >70% genetic homology with humans, and are vertebrates yielding structural similarity (Howe et al., 2013). Because zebrafish eggs are laid and fertilized externally, embryonic growth can be easily visualized, unlike in mammalian models. They share a protruding yolk sac structure with human embryos, unlike rodent models which have inverted yolk sacs (reviewed by Sant and Timme-Laragy, 2018). Unlike humans, the zebrafish yolk is purely lecithotrophic, meaning that a finite amount of yolk is supplied prior to external fertilization (Kunz-Ramsay, 2013). The composition of the yolk is well-characterized, composed mostly of lipids, such as cholesterol and triglycerides, and bulk proteins, such as vitellogenins (Anderson et al., 2011; Fraher et al., 2016; Miyares et al., 2014; Quinlivan and Farber, 2017). Therefore, the zebrafish is a unique, well-controlled model in which to study the fundamental contributions of yolk nutrition to growth and development.

From zygote until about 5 days post-fertilization (dpf), zebrafish development is reliant on yolk sac metabolism for its energy. Once feeding begins to take place, the yolk plays a secondary role in nutrient uptake until full yolk absorption at 7 dpf (Wilson, 2012). Previous studies have investigated yolk sac volume perturbations at particular time points, such as 5 or 7 dpf, as a model for the amount of nutrients used (Kalasekar et al., 2015). It is believed these time points are sufficient for reasonable conclusions; however, this study aims to test this belief with a continuous time course investigation. The addition of multiple discrete time points, in comparison to just day 5 or 7, makes it possible to analyze continuous time through ordinary differential equation modeling. Not only is it important to understand how nutrient absorption occurs over time, but also how the fish is growing over time and the interaction of nutrient absorption and growth during this dynamic period of development. To address this, we used mathematical models to quantify yolk absorption from time-lapse fluorescence microscopy images of zebrafish.

Mathematical modeling is an excellent method that has been used successfully to study many developmental processes (Fletcher et al., 2017; George et al., 2015; George and Lubkin, 2018; Vasieva et al., 2013). It has significantly improved our understanding of zebrafish development (Sharpe, 2017; Wertheim and Roose, 2017; Yamaguchi et al., 2007; Zhu et al., 2017). The majority of the mathematical models that describe zebrafish development have been used to investigate the movement of pigment cells, pigment patterns and formation of zebrafish stripes (McGuirl et al., 2020; Volkening and Sandstede, 2018) and the mathematical models that have been developed to study the growth dynamics of zebrafish focus on embryonic (0-3 dpf) or early larval periods (Zhu et al., 2017). Augustine et al. (2011) described fish length over time throughout the entire juvenile stage of the fish (until about 80 dpf) using ‘dynamic energy budget’ (DEB) theory, which quantifies the uptake and use of food. DEB theory requires excess data to complete analysis and is reliant on stage transitions (e.g. larval to juvenile) to produce growth predictions, making the model invalid for the sole investigation of the embryonic stage. Investigation of yolk and length interaction during the embryonic stage has yet to be completed using continuous time mathematical models. Ramlan et al. (2017) examined fish length and yolk diameter at only 1, 2 and 3 dpf in order to evaluate the effects of ethanol exposure in an experimental group. No mathematical models were used in the study by Ramlan et al. (2017) to describe yolk and fish length independently or associatively.

In this study, we use mathematical techniques, such as differential equations modeling, image analyses, curve fittings and parameter estimation, to quantify the rate of yolk absorption over time and its effect on early fish development from 1 to 7 dpf. Findings from this study establish a nonlinear functional relationship between nutrient absorption and early fish growth. We found that the graph of the size of the fish yolk over time is a logistic decay curve. That is, it decays rapidly for a period of time and then levels out. This implies that the rate of yolk utilization in zebrafish is nonlinear.

Understanding the time at which yolk usage has its greatest impact on fish growth is of paramount importance to the DOHaD paradigm, providing a critical window of nutritional sensitivity and susceptibility during which nutritional deviance could yield the highest likelihood of adverse health effects. The goal of this study is to resolve the temporal relationship between embryonic nutrition and growth in a vertebrate model (Danio rerio), and to characterize the prioritization of yolk nutrition towards overall growth throughout the developmental period.

The amount of nutrients used by the zebrafish can be inferred from the amount of yolk used at a particular time point. Using the yolk area as the measure for embryonic nutritional uptake, we determined a mathematical model that explains the yolk change during 1-7 dpf. Embryos of adult wild-type zebrafish (AB strain) were used in this study. Microscopy images from 14 zebrafish were obtained for each day of the experimentation. The yolk area in mm² and fish length in mm was computed (Fig. 1), and each data point was normalized according to the maximum value in each set. The experimental data for yolk area and fish length are shown in Tables S1 and S2, respectively. The normalized data for yolk size and fish length are shown in Tables S3 and S4, respectively. Normalized data were used for the modeling.

The raw experimental data of yolk area plotted over time foreshadows a decreasing sigmoidal behavior. Here, we examined multiple potential mathematical models in order to conclude with a function that is indubitably the best representation of the experimental data and the inherent biological properties. The biologically relevant model restrictions that we considered were yolk depletion at 7 dpf, consistent depletion of yolk over time with no regeneration of yolk, and the initial yolk size taken at 1 dpf. Therefore, we aimed to determine a function Y(t) that describes the yolk area over time with mathematical properties ⁠,

dY/dt<0 for t∈[1, 7], and Y(1)=Y_ini, where Y_ini is the initial yolk size. The three mathematical models used in this investigation are: the exponential model, Hill function and the logistic model.

As with human and other vertebrate development, successful initial fish growth is an important factor in determining fish longevity and health at the adult stage. A continuous time investigation of fish length over the first 7 days of development was conducted to allow for proper analysis of nutritional effects on this developmental process and the potential adverse health effects. The fish length experimental data depict an increasing trend with a decreased growth rate near the end of the developmental stage. Therefore, we wished to determine a function, L(t), that describes the fish length during days 1-7 with L(1)=L_ini, where L_ini is the initial fish length, and dL/dt>0 for t∈[1, 7] dpf. In this section, we investigated two different mathematical models in order to confidently describe the experimental data in the best way possible.

Estimated parameters, shown in Table S5, for the yolk model (Eqns 1-3) were obtained by fitting each model to individual fish experimental data. The graphs for the individual fish yolk samples are shown in Fig. 2 (for logistic model, Eqn 3), Fig. S1 (for exponential model, Eqn 1) and Fig. S2 (for Hill function, Eqn 2). Plots of the mean yolk data±s.e.m. are shown along with the mean best-fit curve in Fig. S3. After computing the errors for each model describing yolk area over time (Fig. S3), we determined the logistic decay model to be the most representative of the data (exponential model: R²=0.86; Hill function: R²=0.98; logistic model: R²=0.99). We note that even though the logistic model is the best fit for the data, the Hill function is also an excellent fit with R²=0.98. This is because both the Hill function and logistic model have sigmoidal curves and therefore can describe the sigmoidal curve for the yolk data.

The logistic model has a mean decay rate, c, of 1.5 time⁻¹ and mean dimensionless carrying capacity, k₁, of 0.92. This model shows that nutrition is being used at a faster rate during the initial stages of development but less nutrition is available near the end, as expected (Fig. 2). An interesting finding is that the rate of yolk utilization (i.e. the rate of change in yolk with respect to time) is nonlinear, inferred from the slope of the graph of the logistic model. Additionally, the logistic decay model yields insight into how fast the yolk is being absorbed.

Parameter values for the models for fish length (Eqns 4 and 5) were estimated by fitting the models to the data for each individual fish sample. The graphs for the individual fish length plotted against dpf for individual fish samples are shown in Fig. 3 (for logistic model, Eqn 5) and in Fig. S4 (for Hill function, Eqn 4).

Parameter estimates found using the least squares method are shown in Table S6 and plots showing each best fit model and experimental data±s.e.m. to fit the data are shown in Fig. S5. Both models taken into consideration visually fit the model well (Fig. S5), but the logistic growth model (Eqn 5) yields the least error (Hill function: R²=0.98; logistic model: R²=0.99). We made note of the logistic behavior of the fish length to be investigated in subsequent sections. Our findings showed that there is rapid growth at the early stages of development but the rate slows down near the end of the developmental (Fig. 3) time range we considered in the study. We aimed to further investigate the link between the rapid early growth and the rapid yolk absorption taking place between 1 and 4 dpf.

Based on the preceding results, we conclude that yolk size and fish length vary with time in a logistic pattern. The equations best describing their change over time are given by Eqns 3 and 5. Knowing exactly how yolk size and fish length vary independently of each other will pave the way for studying their interactions. Here, we develop a system of two differential equations that describe the effects one developmental process has on the other.

Parameter estimation for the coupled model was performed using the least squares regression method as described previously. The model to fit the data and the experimental data±s.e.m. is shown in Fig. 4 with parameter estimates for each fish shown in Table S7.

Initial parameter estimates were inferred from the computed individual function parameters in Eqns 1-3. The Runge–Kutta-4 method was used to solve the differential equations to ensure accuracy and efficiency. The two-dimensional ordinary differential equation model (Eqn 11) proved to be a novel and accurate representation of experimental data (R²=0.99, Fig. 4).

The system of differential equations (Eqn 11) predicts that yolk depletion that fuels fish length is logistic, whereas the yolk depletion that fuels other organs and/or decay within the body is exponential. In Eqn 11, α represents the rate at which nutrients are being used elsewhere in the body. The average value of α was 0.058 day⁻¹ computed from data in Table S7. In Eqn 11, the parameter β describes rate at which the interaction between yolk and length occurs. Similarly, the average value of β was 1.55 day⁻¹, computed from data in Table S7. The average values for γ, K₁ and K₂ were 1.47 day⁻¹, 0.891 and 0.998, respectively. The variability of parameter values across samples for the various models was determined by the percentage of sample parameters that fall within one standard deviation from the mean. This was done for each parameter in the logistic yolk over time model, the logistic length over time model and the ODE system model. The models are robust, and the variability of each parameter value was minimal across samples (Table S8).

The system of differential equations (Eqn 11) was then used to investigate the impact of yolk usage on daily growth.

Understanding the time at which yolk usage has its largest impact on fish growth is of paramount importance because nutritional deviance at this key time point could yield the highest likelihood of adverse health effects. The first step in determining the impact on yolk usage is to look at the change in yolk over a 24-h period and compare it to the change in fish length over the same period using experimental data (Fig. 5).

Our results show that there is a strong association between yolk usage and length change between 4 and 5 dpf (P=0.017) but no other time frame showed statistical significance (Fig. 5). To support this, we found the percentage change for each fish individually and determined that the same time frame is statistically significant (P=0.04) with no other time points demonstrating statistical significance. This discrete time analysis shows that yolk metabolism between 4 and 5 dpf plays a large role in successful fish growth but fails to give a precise time at which the rate of yolk usage is at its greatest.

In order to further study the rate of yolk utilization over time, we plotted |dY/dt| from Eqn 11 over time in days. This represents the rate of yolk utilization over time rather than the yolk area over time as displayed previously. After substituting the values for the parameters from the data fitting, we obtained the functions Y(t) and L(t) by solving the differential Eqn 11. These were used to determine the magnitude of the rate of yolk utilization |dY(t)/dt| from Eqn 11 over time in days. This was done for every sample and the mean values are plotted in Fig. 6A. The overall rate of yolk utilization, |dY/dt|, is nonlinear with maximum utilization occurring at 3.4 dpf.

In Fig. 6B, the two curves on the graph represent the rate of yolk utilization based on both terms in the differential equation evaluated separately. Thus, the rate of yolk utilization that drives the increase in fish length is dY₁/dt=−βL(t)Y(t)(1−Y(t)/K₁) whereas the rate of yolk utilization for other processes including yolk decay is dY₂/dt=− αY(t), where dY/dt=dY₁/dt+dY₂/dt. We plotted computed values for |dY₁/dt| and |dY₂/dt| against dpf (Fig. 6B). The plots were generated in a similar way as described above for Fig. 6A, where Y(t) and L(t) are functions of time. The rate of change in yolk necessary for fish growth |dY₁/dt| is nonlinear with maximum value at 3.5 dpf. The rate of yolk decay, |dY₂/dt|, is nonlinear in time with maximum value at 1 dpf. In addition, |dY₂/dt| is larger in value than |dY₁/dt| until approximately 2 dpf (Fig. 6B) after which |dY₂/dt| becomes larger. This implies that between 1 and 2 dpf the fraction of yolk depleted to fuel the increase in fish length was smaller than that used for other processes in the body.

Through this investigation, we were able to determine the exact time at which yolk utilization reaches its maximum for nutrition fueling the increase in fish length versus other processes, including yolk decay. In both cases, we observed that the maximum rate of utilization, (|dY/dt|, |dY₁/dt|), occurs between 3 and 4 dpf (3.4 dpf and 3.5 dpf, respectively). The precise time at which maximum utilization occurs (3.5 dpf) assists in the identification of utilization sensitivity, which can be informative when determining sensitive time points for further analysis. This finding from the system of differential equations complements the results from the raw data statistical analysis approach (Fig. 5), in which the yolk usage between 4 and 5 dpf was found to be highly correlated with fish growth during 4-5 dpf. We see that once the maximum utilization rate occurs, the impact on fish growth is increased with a slight delay. These findings suggest there may be a delay between yolk usage and fish growth, which is to be discussed in a later section.

The model has two ultimate potential equilibrium points, one at zero and another at its carrying capacity (⁠⁠). Algebraic stability analysis reveals that the zero equilibrium is stable whereas the maximum yolk equilibrium is unstable (see supplementary Materials and Methods). The change in length over time is given in Eqn 5 with its dynamical form given in supplementary Materials and Methods. In this case, the model has similar equilibrium points as the yolk model at zero and its carrying capacity ⁠. Whereas the yolk utilization over time reaches the zero equilibrium, algebraic stability analysis reveals stability of the maximum fish length, ⁠. This is consistent with our expectation, because we do not expect the fish length to increase or decrease after yolk depletion (i.e. in the absence of outer sources of nutrients).

Through our results, we aimed to lay a foundation for further pharmacological modification studies investigating yolk utilization and fish growth rates. However, it is unknown whether our model (Eqn 11) is still applicable under conditions of stress and environmental modulation. To test this, we also exposed embryos daily to a pollutant, perfluorooctanesulfonic acid (PFOS, 32 µM). PFOS has been previously shown to disrupt yolk utilization (Jantzen et al., 2016; Sant et al., 2017, 2018), with no significant impact on fish length at the selected concentration (see supplementary Materials and Methods). Therefore, PFOS was an ideal candidate to test the robustness of our mathematical model because of the known outcome and because we can test whether our model may examine independently controlled variables. We performed a nonlinear least squares regression analysis to determine the parameter values (Table S9) for which Eqn 11 best fits the mean sample data. There was no significant increase in heterogeneity (variance) due to exposures (Levene's test; P>0.05) and the system of differential equations (Eqn 11) is still representative of the data (R²=0.99) (Fig. S7). PFOS samples displayed an increase in yolk utilization as shown by the best fit values for the rate parameters α and β (Table S9, Fig. S7). α displayed an increase of 84.47% and β increased by 51.28%. Therefore, the mathematical model (Eqn 11) can capture how pharmacological modifications affect the per capita yolk utilization rates that fuel the fish growth (β) and also the decay due to yolk use elsewhere in the body (α). This is interesting, because it demonstrates that exposure of zebrafish to the pollutant PFOS leads to the yolk being used up at an increased rate compared with control. The average rate parameter γ (i.e. per capita growth rate of zebrafish) displayed a 1.17% decrease in PFOS-exposed samples (Table S9). This little to no change in γ between controls and PFOS samples demonstrates the models’ ability to capture the similarities in fish length growth behavior between groups (Fig. S7). These results support our previous crude findings that yolk utilization, and not fish length, are impacted by this concentration of PFOS (Sant et al., 2017; Sant et al., 2018). The extension of the model (Eqn 11) to a ‘stressed’ or non-control environment demonstrates that even in the presence of toxicological modifications, logistic behavior is expected and that our model is robust enough to withstand environmental perturbation.

Embryonic development is a dynamic process, requiring a coordinated and energy-intensive continuum of events crucial for success. During this period, embryogenesis and later organogenesis lay a pattern that serves as a foundation for health throughout life.

The DOHaD paradigm has repeatedly associated the relationship between developmental nutrition, in utero environment, and prenatal growth with chronic disease outcomes. Here, we assess the relationship between embryonic nutrition and growth using the zebrafish model in order to better characterize these relationships temporally and identify key windows of metabolic development.

The impact of yolk utilization over time is a question unanswered by previous literature. Our findings demonstrate the continual time interactions between yolk utilization and fish growth through a system of two differential equations. The mathematical model is novel and it has been validated by experimental observations, making it a powerful tool for replication and future simulations. Our models show a strong association between the rates of change in yolk and fish length between 4 and 5 dpf (P=0.017), highlighting that yolk metabolism between 4 and 5 dpf plays a crucial role in successful fish growth. Therefore, disruption of yolk metabolism during this period may have the greatest impact on larval growth, and this is a critically sensitive metabolic window for DOHaD in zebrafish studies.

Zebrafish length over time has been investigated in various toxicology and biology studies, in which a mathematical model may or may not be reported to describe the experimental data (Collery et al., 2014; Dang et al., 2018; Ho and Burggren, 2012; Imrie and Sadler, 2010). Additionally, the time frame for fish length measurements varies greatly among studies with the majority of measurements being taken cross-sectionally during the embryonic stage or much later in the larval or juvenile stages (Augustine et al., 2011). As a result, the embryonic time frame lacks fish length mathematical modeling and studies investigating later time points possess few models to describe the experimental data, although one study is of particular interest (Zhu et al., 2017). Zebrafish length and body mass was assessed at 180 days post-hatching, where hatching takes place between 2 and 3 dpf (Kimmel et al., 1995). The experimental data for the study was fit to a four-parametric logistic model and it was concluded that the fish length over time behaves logistically (Zhu et al., 2017). With the embryonic stage also behaving logistically, as described extensively throughout this study, we see there may be a correlation between early and juvenile fish length development of the zebrafish. Early development studies that investigate length over time lack the necessary time points for the developed mathematical model of this study to be applied (Dang et al., 2018; Schäfers et al., 2007). One study in particular reported paired yolk and length data at 1, 2 and 3 dpf but the lack of further time points deem the data impractical for a continuous time analysis over the embryonic and larval stages (Ramlan et al., 2017). The need for a time course model of fish length throughout these embryonic and larval stages is apparent and the model proposed in this study lays a foundation for a growing body of literature surrounding this topic.

Zebrafish are an important model for developmental biology and toxicological research. Because zebrafish are lecithotrophic, the amount and rate of yolk utilization can be an important indicator of adaptive mechanisms due to environmental stressors, such as toxicant exposures or environmental modulation (temperature, pH, etc.). Yolk use is a common endpoint used in toxicological testing, and therefore we needed to ensure that our model was robust enough to withstand these stressors. Here, we showed that our model could indeed account for toxicant-induced stress (Fig. S7), and that sample heterogeneity, including differences in sample means between groups due to treatment, could be readily accounted for. Therefore, we expect this model to be applicable to diverse studies for use across disciplines, including pharmacology and toxicology.

Zebrafish have been shown to be a high-throughput, simplified model in which to study the relationship between nutrition and growth in vertebrates during the embryonic period (previously reviewed by Sant and Timme-Laragy, 2018). The majority of research relating developmental nutrition with growth outcomes has been largely associative due to the inaccessibility of the intrauterine environment in mammals. Unlike humans, who rely on maternally supplied yolk granules during early embryogenesis and histiotrophic uptake of nutrients following implantation, zebrafish are a clarified nutritional model. Because zebrafish have a finite amount of yolk available until exogenous feeding, the use of this supply over time can serve as a proxy for total embryonic nutrition.

Using our mathematical model, we were able to determine the time at which yolk metabolism fueling fish growth is at its peak. At approximately 84 hpf, the yolk utilization rate reached its maximum, with a 12 h delay until its impact on fish length growth at 96 hpf. Although this window of sensitivity requires additional validation across other vertebrate models using comparative biology (Carnegie stages), this finding is especially interesting biologically because this is a peak period of endodermal maturation in the zebrafish. During this time window, structural organogenesis wanes, but physiological functions commence. Organs and tissues that are vital for nutrient uptake, distribution and metabolism, such as the liver, become more physiologically and metabolically active between 4 and 5 dpf (reviewed by Goessling and Sadler, 2015; Wilkins and Pack, 2013). Although we cannot make any conclusions from this study about chronic disease and the relationship between exogenous feeding and growth after the developmental period, this work further supports that 4-5 dpf may be a critical metabolic window for growth and development in zebrafish. Furthermore, this window is also likely then to be a window of susceptibility (within the DOHaD paradigm) during which nutritional, toxicological or pharmacological exposures and deficiencies may produce deleterious effects that could persist throughout life.

Augustine et al. (2011) calculated the energy of the zebrafish using DEB theory. This method uses data on observed fish lengths, weights, reproduction and survival to estimate the cumulative energy invested in development. This cumulative maturity level is measured in millijoules and was found for the entirety of the juvenile life stage for the fish (day 1-218). For 1-7 dpf, maturity was estimated to be a linear process, which differs from the logistic behavior seen here. In our experiment, the energy invested in development was taken from the yolk area and we assume use the sole use of this deposit as energy during development. This differs from the study of Augustine et al. (2011), in which they assume a smooth transition from yolk usage to food consumption without a perturbation of the maturity level. The simplicity of the model described in this study is valuable for its accessibility for usage even with small data sets. The accessibility of the model will also allow for widespread usage of the model to investigate nutritional deviation (perturbation or supplementation) during development using measurements that are time and cost efficient (microscopy-based measurements). Further research and data collection should be executed to extend the proposed model to include the complex processes the DEB model considers while keeping the necessary data for conclusions to a minimum.

We acknowledge some limitations of this study. The mathematical model assumes that yolk metabolism occurs largely to fuel fish length, thereby using fish length as a measure of overall fish growth. Future studies obtaining growth data from various other organs could aid in the extension of the mathematical model to better incorporate growth in other areas of the body. Additionally, using 2D area as a proxy for total 3D yolk limits the exact evaluation of yolk present at a particular time. We acknowledge that other 3D methods are available, such as Oil Red O staining, that would help with the understanding of lipid-rich yolk, but Oil Red O staining would not allow for longitudinal assessment over time (Kalasekar et al., 2015). Note that when we use 2D area measurements, we measure yolk size in the x and y plane, but its size in the z plane is not taken into account (Fig. S8A,B). In the early stages of development, particularly at 2-3 dpf, the fish yolk may easily be approximated by a combination of an ellipsoid and a cylindrical tube (Fig. S8A,B). The tube at this stage may add a disproportionate amount to the area when using 2D measurements as a proxy for total 3D yolk. This is because the ellipsoidal portion holds much more mass and this may not be captured in the 2D plane (Fig. S8B). We estimated the yolk volume from sagittal and coronal images of zebrafish (Fig. S8A,B). We found that the yolk volume and yolk area are positively correlated at 2-3 dpf when the yolk has a significant depth not captured by a 2D image (see supplementary Materials and Methods and Fig. S8C). From this, we state that 2D modeling of the yolk, although not perfect, is an acceptable proxy for yolk volume in this study. Statistical analysis reveals a strong association (R²>0.85, P<0.01).

We understand that our model is only applicable to embryos and early larvae prior to the onset of exogenous feeding (during 100% yolk-supplied, lecithotrophic nutrition). Beyond that dynamic period, yolk is no longer available nor a viable predictor of overall growth. Additionally, parameter estimates made by the model are based on a relatively small data set and are setting a foundation that has not been laid before. Yolk utilization has previously been thought to be a labile process that is sensitive to perturbation, but the relationship between yolk utilization and time has not been well characterized (reviewed by Sant and Timme-Laragy, 2018). Here, we conclude that yolk utilization occurs logistically. Without a foundation for the logistic behavior, we identified the best mathematical models by using a combination of data from laboratory experiments, parameter estimations and least squares error analysis, rather than building on findings from existing literature. It is also of note the model has been made for a wild-type (AB) zebrafish line. Other strains of zebrafish may behave differently, and a larger data set is needed to overcome this limitation. Lastly, the model is able to make predictions on the first 7 days of development, without the added variable of exogenous feeding. Growth occurring from food intake may be incorporated for future extensions of the mathematical model.

In this study, we characterize the relationships between embryonic growth and nutrition in a controlled model. Unlike mammalian models, we are able to exclude exogenous nutrient sources in zebrafish. Although this is a simplified model that allows us to reduce the number of convoluted or extraneous variables, we are unable to account for any adaptive responses which may occur as a result. ‘Catch-up growth’ and adaptive changes in feeding behaviors, during which rapid growth is observed followed by a period of impaired or stunted growth, are well-documented pediatric observations (Boersma and Wit, 1997). Catch-up growth has been widely studied, and has been associated with increased risk for chronic disease later in life (Eriksson et al., 1999; Huxley et al., 2000; Ong et al., 2000). Because this adaptive response is a potential mechanism of the DOHaD paradigm, future studies should examine flux following the onset of exogenous diets (after 5-7 dpf). Because zebrafish do not rely on lactation and instead consume purely exogenous diets, they could be a unique tool to examine this relationship.

Here, we examined the relationship between nutrition and growth under ideal laboratory conditions for the development of zebrafish. Although this reductionist approach allowed us to construct an ideal model, this does not account for genetic or environmental diversity, which may contribute to physiological perturbations. Deviations from these ideal conditions, such as nutritional modulation (under- and over-nutrition, yolk composition and rate of use) or exposures to relevant pharmacological or toxicological agents may hinder or exacerbate the rate or timing of the relationships proposed in this study. We have previously shown that exposures to toxicants and pharmacological agents can impair embryonic uptake and availability of nutrients in rodent models, and that these exposures are associated with hindered growth (Harris et al., 2015; Sant et al., 2016a; Sant et al., 2016b; Sant et al., 2013). In future work, our laboratory will assess how these factors impact these relationships in order to show the robustness and sensitivity of our model.

In conclusion, this study applies mathematical techniques, such as differential equations modeling, image analyses, curve fittings and parameter estimation, to quantify the rate of yolk absorption over time and its effect on early fish development. It identifies a nonlinear functional relationship between nutrient absorption and early fish growth, and a critical window of yolk utilization at approximately 84 hpf. We highlight a susceptible time point during development at which exposures and nutrition may have deleterious and potentially lasting effects.

Animal care and use was performed in strict accordance with protocols approved by the San Diego State University Institutional Animal Care and Use Committee (PHS Assurance Number 16-00430). Adult wild-type zebrafish (AB strain) were obtained from existing populations at the University of California San Diego and maintained as a colony at San Diego State University. Zebrafish were maintained in an automated zebrafish habitat (Aquaneering), at 28.5°C and on a 12 h light:12 h dark cycle. Water quality conditions were monitored throughout the study, and safe levels of chemicals including total chlorine, ammonia, nitrites and nitrates were maintained. Adult zebrafish were fed a Gemma Micro 300 powdered diet (Skretting) daily, according to manufacturer's recommendations.

Embryos for each experiment were randomized from six breeding tanks containing 20 females and ten males in order to reduce clutch effects. Embryos were collected and confirmed for fertilization prior to 1 hpf. Embryos were rinsed and transferred to 100 mm polystyrene Petri dishes containing 0.3× Danieau's medium [17 mM NaCl, 2 mM KCl, 0.12 mM MgSO₄, 1.8 mM Ca(NO₃)₂, 1.5 mM HEPES, pH 7.6], and placed into a light-cycling incubator (28.5°C) overnight. At 24 hpf, embryos were manually dechorionated using watchmaker's forceps and individually transferred into wells of a polystyrene 24-well plate containing 1 ml of 0.3× Danieau's medium supplemented with 0.01% dimethyl sulfoxide (DMSO; Thermo Fisher Scientific, BP231-100). DMSO was utilized to mimic typical vehicle control conditions for zebrafish studies, and the 0.01% v/v concentration has been previously shown to be non-toxic and have minimal impact on growth and development (Hallare et al., 2006; Maes et al., 2012; Sant et al., 2016c). Media was refreshed daily to prevent hypoxia throughout the study. To test the applicability of our model to an environmental stressor, we additionally exposed embryos to a contaminant, heptadecafluorooctanesulfonic acid (PFOS) (Sigma-Aldrich, 33607). Additional experimental details related to exposures are provided in the supplementary Materials and Methods.

All microscopy was performed using a Nikon Ti-2 inverted light microscope using Nikon NIS Elements Advanced software. Prior to imaging, embryos and larvae were rinsed thoroughly in 0.3× Danieau's and briefly anesthetized in 0.168 mg/ml Tricaine-S MS-222 solution. Embryos and larvae were rinsed then transferred to drops of 3% methylcellulose for imaging. Images were acquired at 20× and 40× magnification. At the conclusion of imaging, embryos and larvae were thoroughly rinsed in fresh 0.3× Danieau's medium and transferred back to their wells with fresh media. We used 2D imaging in order to decrease handling and anesthesia time for each fish and to support high-throughput microscopy for obtaining data at successive time points.

Visualization of the developmental processes in vivo allows for accurate analyses of nutritional uptake and fish growth. Microscopic images were used to determine yolk area and fish length in zebrafish each day during experimentation beginning at 24 hpf (n=14). Measurements were taken using the NIS Elements Advanced Research software to extract yolk area in mm² and fish length in mm. Fish length is defined as total anterior-posterior length. The experimental data for yolk area and fish length are shown in Tables S1 and S2, respectively. For analysis, each data point was normalized according to the maximum value in each set. The normalized data for yolk size and fish length are shown in Tables S3 and S4, respectively.

Data are presented in mean fit figures as mean±s.e.m. Statistical tests for mean differences were performed using an independent two-sample t-test. For linear regression analysis, testing for association between two variables, independent ANOVAs were used to test whether the slope of the regression line is significantly different from zero using R software. In all cases, a 95% confidence level was utilized (α=0.05) to deduce statistical significance. All data used in this study are normally distributed (Kolmogorov–Smirnov test, α=0.05) and consistent with historical laboratory data.