Development modeling of Lucilia sericata (Diptera: Calliphoridae)

The relationship between insect development and temperature has been well established and has a wide range of uses, including the use of blow flies for postmortem (PMI) interval estimations in death investigations. To use insects in estimating PMI, we must be able to determine the insect age at the time of discovery and backtrack to time of oviposition. Unfortunately, existing development models of forensically important insects are only linear approximations and do not take into account the curvilinear properties experienced at extreme temperatures. A series of experiments were conducted with Lucilia sericata, a forensically important blow fly species, that met the requirements needed to create statistically valid development models. Experiments were conducted over 11 temperatures (7.5 to 32.5 °C, at 2.5 °C) with a 16:8 L:D cycle. Experimental units contained 20 eggs, 10 g beef liver, and 2.5 cm of pine shavings. Each life stage (egg to adult) had five sampling times. Each sampling time was replicated four times, for a total of 20 measurements per life stage. For each sampling time, the cups were pulled from the chambers and the stage of each maggot was documented morphologically through posterior spiracle slits and cephalopharyngeal skeletal development. Data were normally distributed with the later larval stages (L3f, L3m) having the most variation within and transitioning between stages. The biological minimum was between 7.5 °C and 10 °C, with little egg development and no egg emergence at 7.5 °C. Temperature-induced mortality was highest from 10.0 to 17.5 °C and 32.5 °C. The development data generated illustrates the advantages of large datasets in modeling Lucilia sericata development and the need for curvilinear models in describing development at environmental temperatures near the biological minima and maxima.


INTRODUCTION
In the past, insect development research has focused on agriculture pests and disease vectors, with temporal accuracy levels of days (or weeks) considered acceptable since the focus was on economic thresholds and vector prevention (Higley & Haskell, 2010). In the last 50 years, however, there has been a growing interest in the development of necrophagous insects. Unfortunately, unlike agricultural and medically important insects, whose biology had been studied down to eye color during pupal stages and mode of infection from digestive track to mouthparts, necrophagous flies had no such The solution is to explicitly limit the temperature range for which a degree-day model is valid to linear portions of the development curve (Nabity, Higley & Heng-Moss, 2007).
Another pertinent methodological point is how to replicate temperature treatments. By definition, an experimental unit is the thing to which a treatment is applied, and experimental units must have independent replication. Because temperatures are applied in growth chambers, the chamber is (by definition) the experimental unit (Snedecor & Cochran, 1989). Measurements of within growth chamber temperature variations and systematic differences between experienced and nominal chamber temperatures (Nabity, Higley & Heng-Moss, 2007) demonstrate that chamber replication and/or temperature measurements within chamber locations are necessary to provide proper estimates of experimental error and to avoid systematic measurement errors. Richards & Villet (2008) discuss how deficiencies in development data can reduce the accuracy of PMI estimations. In particular, sampling errors and models based on too few temperatures directly impact statistical validity and error rates. Besides the obvious reason of unknown or high error rates making it difficult to reach conclusions, known error rates are a requirement for Federal Rules of Evidence (Rule 702: Testimony by Expert Witnesses) under the Daubert standard. Judges can use the "the known or potential rate of error of the technique or theory when applied" as an assessment of reliability of the evidence being presented.
Many of the issues associated with current data sets likely relates to the sheer amount of resources required to establish a complete, statistically valid development model. Preparatory work averages between 15 and 18 h per temperature. These hours include cutting weighed liver, labeling and organizing experimental units, counting eggs, and putting all units together before they go in the chambers. After set up, the hours required for actual sampling can easily exceed 120 h (at an average of one hour per sampling time). A considerable time (not included in the hours above) is required for colony maintenance. Multiple colonies are required to for high egg production and to prevent time loss from loss of a colony.
Thus, although the experiments are time and labor intensive, we conducted a series of experiments that cover a broad range of temperatures, have large sample sizes, and have the consistent sampling times required to create a statistically valid development model.

METHODS AND MATERIALS Flies
Lucilia sericata were obtained from colonies maintained at the University of Nebraska-Lincoln (Lincoln, Nebraska). Colonies were established in October 2010, from fieldcollected insects from Morgantown, West Virginia. At research time, the colonies had achieved 100 generations without addition of new flies to reduce genetic variation within the colony. Adult flies were maintained in screen cages (46 cm × 46 cm × 46 cm) (Bioquip Products, California) in a rearing room at 27.5 • C (±3 • C), with a 16:8 (L:D) photoperiod. Multiple generations were maintained in a single cage, and ca. 1,000 adult flies were introduced every 1-2 weeks. Adults had access to granulated sugar and water ad libitum, and raw beef liver for protein and as an ovipositional substrate. After egg laying, eggs and liver were placed in an 89 ml plastic cup, which was surrounded by pine shavings in a 1.7 L plastic box. The pine shavings served as a pupation substrate. The 1.7 L box was placed in a I30-BLL Percival biological incubator (Percival Scientific, Inc., Perry, Iowa, USA) set at 26 • C (±1.5 • C). After eclosion, adults were released into the screened cages.

Incubators
Incubator information has been previously discussed in Lein (2013). Pertinent information has been revisited here. Incubators were customized model SMY04-1 DigiTherm ® CirKinetics Incubators (TriTech Research, Inc., Los Angeles, California, USA). The DigiTherm ® CirKinetics Incubator have microprocessor controlled temperature regulation, internal lighting, recirculating air system (to help maintain humidity), and use a thermoelectric heat pump (rather than coolant and condenser as is typical with larger incubators and growth chambers). Customizations included the addition of a data port, vertical lighting (so all shelves were illuminated), and an additional internal fan. The manufacturer's specifications indicate an operational range of 10-60 • C ± 0.1 • C. It is worth noting that a range of ±0.1 • C is an order of magnitude more precise than is possible in conventional growth chambers. Although growth chambers have been shown to display substantial differences between programmed temperatures and actual internal temperatures (Nabity, Higley & Heng-Moss, 2007), we tested the customized DigiTherm ® CirKinetics incubators in a replicated study and found internal temperatures on all shelves within incubators did not vary more than 0.4 • C from the programmed temperature (data not shown). Given this measured accuracy, the incubators could be used without the need for additional internal temperature measurements or risk of systematic error (as there was <0.4 • C internal variation from the programmed temperature).

Experimental Design
The study comprised eleven temperatures (7.5, 10, 12.5, 15, 17.5, 20, 22.5, 25, 27.5, 30, and 32.5 • C) with a light:dark cycle of 16:8. Twenty eggs (collected within 30 min of oviposition) were counted onto a moist black filter paper triangle and placed in direct contact with 10 g of beef liver in a 29.5 mL plastic cup. The cup was placed in a 7 cm × 7 cm × 10 cm plastic container that had 2.5 cm of wood shavings in the bottom. The container was then placed randomly in an incubator. There were 27 container locations within each chamber (9 locations per shelf). Containers were randomized by chamber using a random number generator in Excel (Microsoft Excel 2007). Each life stage (egg-1st stage, 1st-2nd stage, 2nd-3rd stage, 3rd-3rd migratory, 3rd migratory-pupation, pupation-adult) was calculated using Kamal's (1958) data, which was converted to accumulated degree hours (ADH) and divided equally into five sampling times (Table 1). Each sample was replicated four times, with a total of 20 samples per life stage. During each sample time, a container was pulled from each of the four incubators and the stage of each maggot was documented morphologically using the posterior spiracular slits and cephalopharyngeal skeleton.
During egg hatch, a larva was recorded as 1st stage if they had broken the egg chorion and were actively emerging. Pharate larvae (larvae that have undergone apolysis but not ecdysis) were recorded as the earlier stage (e.g., 3rd stage spiracular slits can be seen beneath the current spiracular slits would be recorded as 2nd stage), since they had not yet molted. Larvae were considered 3rd migratory when they stopped feeding, left the liver, and began burrowing through the pine shavings. Pupariation began when larvae had a shortened body length and no longer projected mouth hooks when put in the larval fixative KAAD (kerosene-acetic acid-dioxane). There were times when a larva appeared to be entering the puparium stage but would extend its body length and begin crawling if disturbed or placed in KAAD. These larvae were recorded as 3rd migratory. All life stages were preserved in 70% ethyl alcohol. Third and 3rd-migratory stages were fixed in KAAD for 48 h and transferred to 70% ethyl alcohol.

Analysis
In the literature, the time in a stage is typically reported as a single number. Because variation exists among individuals in their development times, it is essential to use an appropriate indicator reflecting when individuals transition from one stage to another. With normally distributed variation this estimator is a mean, but with alternative distributions other measures are more appropriate. Consequently, the distribution of individuals during stage transitions must be determined. Details on stage transition modeling and its importance in properly determining development rates are reported elsewhere (Roe, 2014).
Here, we summarize the procedures used in determining time in stage. First, we modeled stage transitions using TableCurve 2d, version 5.01 (SYSTAT Software Inc, San Jose, California, http://www.sigmaplot.com/products/tablecurve2d/ Tablecurve2d.php), and Prism, version 6.02 (GraphPad Software, Inc., La Jolla, California, http://www.graphpad.com/scientific-software/prism/). We fit one of four Gaussian functions (specifically, a regressed proportion (percentage) in stage versus time, at each temperature tested. Then, stage durations were determined, probit models were constructed, and the 50% transition point was determined.
For all regression analyses, the data were examined closely to determine their propriety for inclusion in analysis. In a few instances, individuals were sampled with extraordinarily extended durations. These were treated as outliers and excluded from analysis. Details on all data used are included in Appendix S1.

DEGREE DAYS
Degree day requirements were calculated with a combination of regression analyses and iterative analyses to ensure the resulting degree day models reflected only the linear portion of the insect development curve. The outline of these procedures is: 1. Determine the stage transitions by fitting Gaussian curves to the proportion of insects entering the new stage vs. time for each temperature (curves were calculated for L1, L2, L2f, L3m, P, and A). Only data for the first portion of each curve (0%-100%) was included in the regression which reflects the stage transition.
2. Calculate the 50% transition point from the Gaussian curve for each stage and temperature combination.
3. With data from 2, determine time in stage by subtraction between 50% transition points.
4. Express development times in days (rather than hours, as data was initially determined) and calculate 1/days for each time to transition and stage duration.
5. By linear regression, estimate the relationship between development rate (1/days to transition or stage) vs. temperature) to determine the slope and x-intercept. Each resulting regression was runs tested to identify non-linearity, and where nonlinearity was indicated, points were excluded from the regression until any non-linearity was eliminated. (Runs testing is a procedure by which iterative calculations are used to distinguish linear from non-linear points in a regress (Mutulsky, 1995;GraphPad Software, Inc., 2014)). Primarily, non-linearity was associated with low and high temperatures (as expected) and indicated in development graphs. The regression of 1/days vs. temperature is conventional in degree day determination, but the use of runs testing to identify non-linear points in the regression has not been. To the best of our knowledge this approach was first used in Nabity, Higley & Heng-Moss (2007) to ensure that assumptions underlying degree day analysis were met.
6. From the resulting linear regressions, the x-intercept represents the developmental minimum and 1/slope represents the accumulated degree days required for an event (stage transition or stage duration) (Arnold, 1959). Although this point usually represents the end of most degree day determinations, it is still possible at this point to have included data in the linear regressions that are not properly part of the linear portion of the development curve. Consequently, we did additional calculations and corrections to determine the validity of our degree day models.
Using regression results we calculated degree day accumulations for each experimentally determined combination of temperature and time of transition or stage duration. We then did a linear regression of these data and evaluated the resulting lines for linearity and slope. To meet the core assumption of degree day models, a regression of degree day accumulations must be linear and have no slope. Where our results did not meet these requirements, we removed points (again, at high and low temperatures), and recalculated both the 1/days regression and the accumulated degree days regressions (steps 5-7). We repeated this process until we arrived at linear relationships meeting all degree day assumptions, and noted the range of temperatures for which the resulting equation was valid.

RESULTS AND DISCUSSION
All calculations can be found in Appendix S1. We observed substantial variation in stage transition times (Roe, 2014) and stage durations comparable to that reported by Tarone & Foran (2006). The largest variation was observed during the L3m and pupation stages, regardless of temperature, with the variation largest at 10.0 • C through 17.5 • C (L3m and pupation) and 32.5 • C (L3m) (Fig. 1). These stages are also the longest life stages (by proportion) ( Table 2).
With minimal egg development and no egg eclosion at 7.5 • C, no data were reported. There is evidence, however, that the biological developmental minimum for L. sericata is between 7.5 • C and 10.0 • C, since there were individuals that successfully emerged as adults at 10.0 • C. Although the adults at 10.0 • C and 12.5 • C were normal-sized, the total number of individuals that survived into adulthood was very small compared to the other temperatures. High mortality rates and reduced developmental rates have been reported at 35.0 • C (Ash & Greenberg, 1975), indicating suboptimum temperatures on both ends of the spectrum can impact survivorship and growth. Since extreme temperatures lead to extreme biological variation, there is a disruption to normal gene expression, which can alter the hormones and proteins needed in molting and maintenance. The biological variation at these temperatures may be an inherent variation in L. sericata that allows the species to survive in suboptimal conditions and also successfully maintain populations throughout the world. Based on temperatures measured and methodology, the previous study most directly comparable with our data was that of Kamal (1958). In the early life stages (E to L2), the percent in stage observed here was within 2.1% of that reported by Kamal (1958) (Table 4). In the later life stages there was a difference of 10.5% (L3f) to 16.2% (L3m) ( Table 4). Although we might expect greater variation at later life stages, differences in development times might also reflect differences in methods. Kamal (1958) used constant lighting during his experiments, and Nabity, Higley & Heng-Moss (2007) showed significant delays in developmental under constant light compared to 16:8 L:D cycles. However, the development times are comparable, indicating that geographic variation in development may be less than the inherent variation in development of L. sericata.
In comparison with other previous work (besides Kamal), differences among studies seem likely to be associated with methodology (see also Tarone & Foran, 2006). We say this because sampling times used to determine stage duration were not consistent among previous studies and in some instances may have been inappropriate (Table 5). For example, sampling extensively at the start of a stage would shift the developmental distribution to the left, skewing the mean. Alternately, sampling the largest individuals (a common practice in forensic entomology casework) changes the age cohort data by using the maxima and treating them as normal data points, not outliers (Richards & Villet, 2008).
Our data indicate that stage transitions are integral to generating realistic, accurate development data (Roe, 2014). Transitions occur over a period of hours to days and it is difficult to discuss development without knowing the transition periods between stages, since the vast majority of time in stage is spent as a mixed-age population. Most error occurs in two areas: temperature and stage/age calculations. Being aware of transition times and having consistently spaced sampling times and temperatures lets us attach error margins to the data, reducing the error in the age/stage calculations. By reducing error in one of two areas, we accomplish two things: we can focus attention on the error associated with temperature (fluctuating versus constant, unknown versus close by, etc.) and we can attach known error rates to PMI estimates.
Currently, the data sets available make attaching error rates difficult or impossible, depending on the data presentation. While these data let us generate error rates, they are only accurate for the linear portion of the life stages (or the temperature range that the ADD are valid). The assumption of linearity was met for all life stages, but not all temperatures (Fig. 2). Not surprisingly, L3m has the shortest linear temperature range from 17.5 to 30.0 • C. L1 also has a shortened range from 12.5 to 32.5 • C (Table 3).
In conventional uses of degree days (e.g., Arnold, 1959), using multiple methods to ensure only linear development data are used in determining degree-day models is not undertaken. Presumably this omission has occurred because it is well recognized that degree-day models use assumptions of linearity to describe what is known to be a curvilinear relationship, so approaches for improving accuracy have focused on curvilinear model development (Wagner et al., 1984;Higley & Haskell, 2010) rather than on improving linear degree day accuracy. Additionally, most conventional uses of degree days with insects involve modeling population level phenomena, where other sources of error (particularly in temperature data) and resolution (of days) are such that more precision in how degree-day models are developed may not be warranted. In contrast, with forensic use of degree-day models, the potential inaccuracy associated with including non-linear data in the calculation model introduces systematic error that could easily be significant in using degree days for estimating postmortem intervals. As a practical matter, systematic Table 3 Linear regression results. Linear regression results (from Graph Pad Prism) of Lucilia sericata, with excluded (non-linear) points indicated by an empty cell. Accumulated Degree Days (ADD) were indicated by the slope of the regression line, and the developmental minimum was indicated by the x-intercept value (Arnold, 1959). For comparison, regression-based ADD were compared to mean ADD calculated across temperatures. The range of the linear regression indicates the temperature limits at which the assumption of linearity between temperature and development is valid.

Temperature
Transition ADD by 1/Day Stage ADD by 1/Days  error will result in under or overestimates in the PMI, depending on temperatures errors not of hours but of days. Looking past the specifics reported here on L. sericata, the most far-reaching implication of this work is the recognition that existing data on the development of forensically important insects may not be comprehensive enough for precision in PMI estimates. Irrespective of the type of developmental modeling used, whether linear (degree day) or curvilinear, the ability to make estimates of insect development at a forensically necessary  15,17,19,20,21,22,25,28,30,34 Minimum Beef liver E, L1, L2, L3f, L3m, P resolution of hours or 1-2 days, requires better estimates of stage transitions and data sets than currently exist. Additionally, the high levels of variation we observed in development further illustrate the crucial need for proper replication in developmental work, if we are going to be able to make scientifically valid statements on variation. Consequently, despite the time, difficulty, and expense of comprehensive developmental studies, our results indicate such data for all forensically important blow flies are essential to meet the need for accurate PMI estimates based on insect development.