Comparative analysis of larval growth in Lepidoptera reveals instar-level constraints

1. Juvenile growth trajectories evolve via the interplay of selective pressures on age and size at maturity, and developmental constraints. In insects, the moulting cycle is a major constraint on larval growth

completely impossible (Roff & Fairbairn, 2007). As a result, trait values may not reach optima that would be expected on the grounds of ecologically based selection pressures. For example, it has been argued that viviparity would be optimal for birds under certain conditions, but it has never evolved in birds due to physiological constraints (Blackburn & Evans, 1986;Dunbrack & Ramsay, 1989). In phylogenetic, among-species comparisons, constraints become apparent as conservatism of trait values (see Schwenk & Wagner, 2004;Tammaru, Esperk, Ivanov, & Teder, 2010;Tammaru, Vellau, Esperk, & Teder, 2015), yet conservatism may also result from strong and consistent stabilizing selection acting on a clade.
Age and size at maturity are central life-history traits that are, perhaps universally, under strong ecological selection (Roff, 2002;Stearns, 1992). Yet, juvenile growth trajectories-the proximate determinants of age and size at maturity-are shaped by developmental and physiological mechanisms that may act as evolutionary constraints. Hence, conservative elements of juvenile growth trajectories identified by phylogenetic comparative analyses could reveal constraints on age and size at maturity.
However, research has focused on just a few model species. As a consequence, the missing comparative perspective precludes generalizations about potential constraints on insect growth trajectory evolution.
In insects, ontogenetic growth is restricted to the larval stage, and regular replacement of the exoskeleton in a process called moulting makes the growth trajectory discontinuous. Body mass increases between the moults (i.e. during a larval instar), whereas the dimensions of the strongly sclerotized (i.e. inextensible) parts of the exoskeleton only increase at moults (Chapman, 1998). The number of moults and the size increment of each instar determine the peak mass of the larva, and thereby, largely, the size of an insect at maturity. Age at maturity is directly affected by the number of moults because physiological preparation for moulting and the moulting process itself take time; across all moults, this time forms a considerable part of the duration of the larval stage. Therefore, the time taken to complete a moult and the number of instars have a major influence on the minimum possible age at maturity. Developmental and physiological mechanisms associated with the moulting process have, thus, the potential to act as important constraints on insect life-history evolution (Davidowitz et al., 2016;Nijhout et al., 2006Nijhout et al., , 2010Tammaru, 1998;Tammaru et al., 2010Tammaru et al., , 2015. The insect moulting cycle appears to be conservative both within and among species. The number of larval instars is usually invariant within a species Esperk, Tammaru, Nylin, & Teder, 2007), but varies between insect species (Chapman, 1998;. However, in some species, the number of instars is sexually dimorphic , or it may vary in response to various environmental factors so that larvae growing under poor conditions add instars to reach a size threshold for a successful metamorphosis Greenberg & Ar, 1996;Grunert, Clarke, Ahuja, Eswaran, & Nijhout, 2015).
Within-instar growth patterns also appear relatively conservative as divergent selection on adult body size may only slightly change relative size increments in any particular instar Meister, Hämäläinen, Valdma, Martverk, & Tammaru, 2018). Moreover, the relative size increment often remains constant across larval instars, which is an empirically broadly supported pattern known as Dyar's Rule (Berg & Merritt, 2009;Dyar, 1890;Hutchinson, McNamara, Houston, & Vollrath, 1997). However, there are deviations from the rule (Hutchinson et al., 1997;Kivelä, Friberg, et al., 2016;Meister et al., 2018), so quantification of conservative instar-level patterns in insect growth, as well as deviations from such patterns, may allow the identification of evolutionary constraints, which is needed for advancing our understanding of life-history evolution.
The key to locating instar-level constraints on larval growth in insects is in identifying the reasons for the moulting process itself.
Ultimately, it is evident that an organism with an exoskeleton cannot grow indefinitely without moulting. However, it appears that moulting occurs before the theoretical instar growth limit is reached (cf. Kamimura & Kiuchi, 1998;Nijhout et al., 2006). Most hypotheses concerning moulting are related to allometry. Surface area to volume ratio inevitably changes with the growth of the organism, so various surfaces need to be periodically renewed to maintain physiological efficiency and functionality. Growth efficiency may decrease during an instar because of hypoallometric growth of the midgut (Esperk & Tammaru, 2004;Grunert et al., 2015), a fixed sized or hypoallometrically growing tracheal respiratory system (Callier & Nijhout, 2011;Greenberg & Ar, 1996;Greenlee & Harrison, 2004, 2005, fixed-sized structures used in food acquisition (Hutchinson et al., 1997), or because of the restricted capacity of the cuticle to expand (Bennet-Clark, 1963). Owing to allometry, the extent to which the area of a surface becomes limiting depends on the absolute size of the larva (cf. Kühsel, Brückner, Schmelzle, Heethoff, & Blüthgen, 2017). As these patterns are based on fundamental physical principles, we would expect a high degree of similarity in these relationships across species and instars.
To test the null hypothesis of isometry (i.e. Dyar's Rule; invariance of relative growth patterns across instars), we recorded two key parameters of larval growth trajectories-instar-specific relative mass increments (RMI) and end-of-instar growth deceleration (GD; Figure 1)-from 30 species of Lepidoptera. Our analysis identifies deviations from isometry in the measured traits. We proceed by evaluating phylogenetic conservatism of the examined traits, with the aim to quantify the degree to which these traits may constrain adaptive evolution of age and size at maturity. Finally, we discuss the data in the light of different hypotheses proposed to explain the ultimate causes and proximate triggers of the moulting process.

| Growth trajectory data acquisition
We obtained detailed growth trajectory data from 30 lepidopteran species from 11 families (Table 1). We used earlier-published data on five species (Leimar, 1996;Esperk et al., 2013;Kivelä, Friberg, Wiklund, Leimar, & Gotthard, 2015; see Table 1), but most of the data were collected specifically for the purposes of the present study. Butterflies and moths were captured in Oulu (Finland), Tartu County (Estonia) and Stockholm (Sweden). Captured females were allowed to oviposit on natural host plants in the laboratory (see Table 1 for the host plants used). Neonate larvae were weighed and individually placed in plastic containers with ad libitum natural host plant. A split-brood design was used for larvae reared at the University of Tartu, with broods being split between 16 and 24°C constant temperatures (day length 12 hr) whenever possible (see Table 1). Larvae reared at the University of Oulu experienced a day length of 16 hr and an associated thermal cycle between 18°C (photophase) and 10°C (scotophase). At Stockholm University, extra Pieris napi larvae, in addition to those included in the previously published dataset (Kivelä et al., 2015), were reared at a constant 16°C under a 17 hr day length. All larvae were weighed one to three times per day until pupation, depending on temperature; weighing was more frequent at higher temperatures. Moulting events were recorded both by direct observations and indirectly from recording the growth trajectory (larvae lose mass during moulting; see below for details). Sphingids were not weighed or handled during moulting because handling at this stage may easily be fatal for these species. Once the pupal moult was completed and the pupal cuticle had hardened, pupal mass was measured and sex was determined based on genital slits in the pupal cuticle. We combined all the datasets to test whether the analysed traits depend on temperature.

| Analysed traits
We focused on relative instar-specific mass increments (RMI) and end-of-instar GD. We calculated RMI as the ratio of the peak mass of the focal instar to the initial mass of the same instar. For quantification of GD, we used a modification of the approach used by Esperk and Tammaru (2004) in which GD is estimated by comparing the instar-specific peak mass that would have been attained without GD and the peak mass of the instar actually attained. We extrapolated growth until the time point where mass accumulation ceases in the focal instar by assuming exponential growth instead of the allometric power 2/3 relationship used by Esperk and Tammaru (2004). Then, we calculated the ratio of extrapolated instar-specific peak mass to the observed instar-specific peak mass (see below for details). A higher ratio indicates a stronger deceleration of growth at the end of the instar.
Next, we quantified ontogenetic change in RMIs and GD as regression slopes capturing across-instar changes in the instar-specific characteristics. Ontogenetic change in RMIs was analysed by fitting a linear regression model of ln-transformed RMI values on ln-transformed instar-specific initial masses separately for each individual ( Figure S1). Logarithmic transformations were used to linearize the relationship between the traits and also because this is conventional in analyses of allometric scaling of trait values to body size (see e.g. Kilmer & Rodríguez, 2017). The slopes of these RMI regressions (hereafter the relative mass increment slope, RMIS) capture ontogenetic across-instar change (or lack of it) in the relative amount of growth within an instar. F I G U R E 1 Illustration of an insect larval growth trajectory describing how relative instar-specific mass increment (RMI) and end-of-instar growth deceleration (GD) are estimated. Top panel shows the growth trajectory of an individual Xanthorhoe fluctuata (Lepidoptera: Geometridae) larva chosen arbitrarily for illustrative purposes. The four larval instars are indicated with different colours. The bottom-centre panel shows the same growth trajectory with ln-transformed mass, clearly visualizing moults and associated mass loss. The grey triangles project the differences (in ln-scale) between instar-specific peak mass and initial mass to the bar plot on the left. Each horizontal bar represents the difference between instar-specific ln[peak mass] and ln[initial mass], which is equal to ln[RMI] (i.e. RMI = [instar peak mass]/[instar initial mass] = exp(ln[instar peak mass] -ln[instar initial mass])). The yellow lines overlaying instar-specific growth trajectories are fitted regressions for instar-specific data, but with the latter part of the instar-specific growth trajectory excluded from regression analysis (see text and Figure S2 for details). Consequently, the regression lines extrapolate growth in the latter part of each instar as if growth continued exponentially until the end of the instar. The difference between extrapolated instar-specific peak mass and the observed instar-specific peak mass is illustrated by the vertical yellow lines. These differences are equal to ln(GD) (i.e. GD = [extrapolated instar peak mass]/[observed instar peak mass] = exp(ln[extrapolated instar peak mass] -ln[observed instar peak mass])). The yellow triangles project instar-specific GD values to the bar plot on the right (the bars have been scaled up to make reading easier), the horizontal bars illustrating instar-specific values of GD. The vertical dashed line in the right panel indicates zero (corresponds to one in the arithmetic scale), below which growth in fact accelerates towards the end of an instar, which is the case in the first instar in this example For GD, we tested the robustness of the results to the definition of the exponential growth phase at the beginning of the instar. For this purpose, we alternatively used the first 40%, 50%, 60%, 70% or 80% of within-instar mass measurements (rounded to an integer number of observations) from the growth phase of the instar (i.e. premoult mass loss was excluded) to describe exponential growth. A regression model was then used to estimate the hypothetical peak mass of the focal instar that would be attained if the exponential growth characteristic of the phase of 'free' growth (Grunert et al., 2015;Kivelä, Friberg, et al., 2016;Nijhout et al., 2006; could be maintained through the instar ( Figure S2). We This is because low thresholds (40% and 50%) cause a lot of noise in the data due to uncertainty of the regression used in GD derivation, and 60% and 80% thresholds resulted in analytical problems (non-converging models and unrealistic estimates of phylogenetic signal; see Appendix A in Supporting Information).
Because both RMI and GD were quantified as ratios of masses and we use body mass as a measure of body size, the null hypothesis of isometry would mean constancy of RMI and GD across instars.
Therefore, the isometric scaling exponent (i.e. value of RMIS or GDS) would be zero for both traits, if larvae grew isometrically.
We had to omit some instars of some individuals because moulting events could not always be unambiguously identified (7.2% of instars completed by the studied individuals). This is because, at high temperatures, some individuals moulted so rapidly that the actual moulting process could not be observed. Occasionally, moult-associated loss of mass (larvae cannot feed some time before and during moulting) could not be identified in the growth trajectory either. Therefore, we only analysed individuals from which the beginning and the end of a minimum of three instars could be unambiguously determined.  Kivelä et al. (2015); P. icarus data originally from Leimar (1996).

| DNA extraction and sequencing
i Temperature was 18°C in the daytime and 10°C at night (i.e. thermoperiod of 16 hr 18°C:8 hr 10°C).

TA B L E 1 (Continued)
the phylogeny, the sampling was supplemented with 10 additional closely related species (Table S1). These 10 species were pruned from the final phylogenetic tree used for subsequent analyses. net/Molec ular.htm. This set of markers has proven to be efficient in Lepidoptera phylogenetics  and has also been used in many studies focusing on the taxa included in this study (e.g. Heikkilä, Kaila, Mutanen, Pena, & Wahlberg, 2012;Sihvonen et al., 2011;Zahiri et al., 2010). Altogether, the gene regions used comprised a total of 6,286 base pairs. All sequences for each taxon were manually aligned and edited using BioEdit (Hall, 1999). All DNA sequences are available at the US National Center for Biotechnology Information (NCBI) GenBank (accessions numbers given in Table S2).

| Phylogeny derivation
An ultrametric phylogenetic tree based on the combined genes dataset was inferred using Bayesian inference (BI) through the CIPRES Science Gateway (Miller, Pfeiffer, & Schwartz, 2010). The BI analysis was run in BEAST 1.8.0 (Drummond & Rambaut, 2007).
Base frequencies were estimated, the number of gamma categories was set to six and a randomly generated initial tree was used.
A lognormal relaxed clock was employed and the tree prior was set to speciation: Yule Process. Parameters were estimated using two independent runs of 50 million generations each, and convergence was checked using the program TRACER 1.6 (Rambaut, Suchard,

| Statistical analyses
To analyse among-species variation in RMIS and GDS, we calcu-  The model included an exponential variance function ('varExp'; Pinheiro et al., 2017) to take increasing residual variance with increasing fitted value into account.

TA B L E 2
Parameter estimates of (nonphylogenetic) linear mixed-effects models for analysing variation in RMIS and GDS in relation to temperature. Species was set as a random effect in these models as several individuals per species were included in the analysis, data coming from more than one rearing temperature for some species (see Table 1). The factor 'temperature' was nested within species In a PGLS framework, we first analysed the relationship between RMIS (response) and neonate larval mass (predictor; ln-transformed).
For this, we used the functions 'gls' (package nlme; Pinheiro et al., 2017) and 'corPagel' (package ape; Paradis, Claude, & Strimmer, 2004), the latter of which first estimates (via maximum likelihood) the value of λ (i.e. phylogenetic signal) and then automatically applies the λ value in modelling correlations among residuals, thereby accounting for any phylogenetic non-independence of the observations. Neonate larval mass was set as a predictor because it captures variation in surface area to volume (A:V) ratio in first-instar larvae (high A:V ratio in small larvae), with this initial A:V ratio being potentially important to physiological functions affecting RMIS. Residual variance was much lower in butterflies (Lycaenidae, Nymphalidae and Pieridae; note that Iphiclides feisthamelii [Papilionidae] was not included in this analysis due to missing data on neonate larval mass) than in moths (all other families represented in this study). The division between butterflies and moths represents the deepest phylogenetic split in our sample of species (note that this division is not the case in Lepidoptera as a whole-butterflies nest within moths ). The phylogenetic independence of these lineages in our sample ( Figure   S3; see also Figure S4) allowed us to model residual variance separately for butterflies and moths by adding a 'varIdent' variance function (Pinheiro et al., 2017) with 'major group' (butterfly/moth) as a variance covariate to the model. The model with the variance function was better than the one without it (ΔAICc = 6.73).
Because initial investigation of the raw data indicated variation among taxonomic families ( Figure S5 We analysed variation in GDS similarly as explained above for RMIS, except that peak larval mass (ln-transformed) was used as an explanatory variable in the models instead of neonate larval mass.
This was to test whether GDS depends on the final body size of a species, with peak larval mass being a good index of species-specific body size. The PGLS models for GDS were not heteroscedastic, so no variance functions were added to the models. For consistency, we tested family-level variation similarly as for RMIS (see Figure S3 for the phylogenetic trees used in the analyses), although investigation of the GD estimates ( Figure S6) did not hint at such variation. Finally, we also tested whether GDS values differ from zero by fitting an   (Table S3; see also Figure   S7). In Sphingidae, RMIS was positive, and in the studied species of Notodontidae and Nymphalidae, scaling exponents did not differ from zero (Table S3). In most species with negative RMIS, it was in general slightly more negative when final instar was excluded from RMIS derivation (  There was, however, variation in RMIS among families (Table 3).

| Relative mass increment regression
RMIS was significantly negative in the reference family Endromidae, which was represented by a single species (E. versicolora) in our sample. In Notodontidae and Sphingidae, RMIS deviated significantly from the value in Endromidae, with values being positive in both Notodontidae and Sphingidae (Table 3). RMIS values in the rest of the families did not differ significantly from the value in Endromidae (Table 3; see also   There was no statistical evidence of an association between GDS and ln(peak larval mass) (coefficient for ln There was no evidence of biologically meaningful GDS variation among families (Table S5). Yet, the moderate phylogenetic signal in analysis of GDS in relation to peak larval mass (see PGLS results above) appears not to be a consequence of the data simply being noisy. This is because GDS varied among species when analysing data consisting of individual-specific GDS values with a linear model in which only species was set as a factor (F 29,393 = 3.05, p < 0.0001, R 2 = 0.18). As further evidence of GDS being a biologically meaningful index, GDS values were consistently positive, the only exception being Endromis versicolora (Figure 3; Table S4; see also Figure S8).

| Growth deceleration regression
Across species, GDS values differed statistically from zero as the intercept-only PGLS model estimated GDS to be significantly positive (intercept = 0.0715, SE = 0.0112, t = 6.40, p < 0.0001), and the confidence intervals of most species-specific GDS estimates did not encompass zero (Table S4) Parameter estimates of a phylogenetic generalized least squares model explaining variation in ontogenetic change in instar-specific relative mass increments (RMIS) among families. The phylogenetic tree used in this analysis was modified so that amongfamily correlations were removed, resulting in a polytomy of family lineages (see text and Figure S3 for full details) to avoid problems due to the inclusion of the 'family' factor that partially explains the topology of the phylogenetic tree in the model   (Table S4).

| D ISCUSS I ON
Our comparative analysis of larval growth trajectories in Lepidoptera reveals both evolving and constrained growth trajectory characteristics. We focused on relative instar-specific mass increments (RMIs) and end-of-instar growth deceleration (GD), both of which typically changed across instars independently of thermal regime.
Across-instar change in RMI (measured as RMI slope across instars, i.e. RMIS; see above) varies among taxonomic families, indicating that the mechanism determining RMIS is flexible enough to change at the macroevolutionary scale, yet it appears constrained within families.
Alternatively, the within-family conservativeness of RMIS arises because species belonging to the same family experience similar selection pressures, but this explanation seems unlikely due to the ecological diversity and consequently divergent selection pressures between species within lepidopteran families. In most of the studied families, RMIS was negative, which means that RMI decreases across instars. The opposite pattern, increasing RMI across instars (i.e. positive RMIS), was found in Sphingidae. Conversely, ontogenetic change in GD, measured as GD slope across instars (i.e. GDS; see above), was very similar in all studied families. GDS was generally positive and independent of body size, which means that GD becomes stronger with increasing body size across instars during ontogeny, and this ontogenetic change was similar across species with different body sizes.
It is also worth noting that the allometric exponent of GD (i.e. GDS)

| Constrained evolution and plastic variation in growth trajectories
The evolution of both RMIS and GDS appears constrained. There is variation in RMIS among taxonomic families and, besides that, subtle variation in the way RMI changes across instars in families with negative RMIS, as indicated by an increasing RMI in the last instar in Drepanidae and Geometridae. Ontogenetic change in RMI is consistent within families, suggesting that either some physiological or developmental constraints slow down within-family evolution (see Roff & Fairbairn, 2007) or that species belonging to the same family experience very similar selection pressures on RMIS.
On the other hand, the low variance in GDS among family lineages suggests that GD is constrained by fundamental developmental or physiological mechanisms so that only a subset of theoretically possible phenotypes can be attained (Roff & Fairbairn, 2007). A shared mechanistic basis of the trait would make it evolutionarily conservative across taxa (see Schwenk & Wagner, 2004;Tammaru et al., 2010Tammaru et al., , 2015. Phenotypic plasticity of growth and development (e.g. Esperk et al., 2013;Frazier, Woods, & Harrison, 2001;Gotthard, 2008;Greenberg & Ar, 1996;Välimäki, Kivelä, Mäenpää, & Tammaru, 2013) affects the instar-level characteristics investigated in this study (Grunert et al., 2015;Kivelä et al., 2018;Tammaru et al., 2015). We pooled data from several larval rearing experiments where the rearing conditions were not identical (Table 1), so natural variation in RMI and GD values due to phenotypic plasticity is included in the data.
Nevertheless, we focused on the ontogenetic change of RMI and GD across instars (i.e. RMIS and GDS) rather than the instar-specific trait values as such. This across-instar change appears quite insensitive to environmental conditions, as indicated by only a few statistically significant temperature effects on RMIS and GDS (

| Potential mechanisms underlying the results
Although several mechanisms have been proposed to explain the need to moult for maintaining a high growth rate, there is limited support for many of them. Out of the ultimate hypotheses for moulting, the role that the potentially hypoallometric growth of the midgut has remains speculative with insufficient data (Callier & Nijhout, 2012;Esperk & Tammaru, 2004;Grunert et al., 2015). The size of feeding structures (e.g. mandibles) does not explain growth rate in lepidopteran larvae, indicating that food acquisition rate does not constrain growth (Esperk & Tammaru, 2004). In fact, absolute growth rate increases within an instar despite constant mandible size, and moulting to the next instar does not result in such a growth benefit as would be expected if food acquisition rate limited growth (Esperk & Tammaru, 2004); therefore, the need to renew feeding structures is an unlikely reason for moulting. Moreover, a restricted capacity for extending (i.e. unfolding) the cuticle (Bennet-Clark, 1963) seems to be an unlikely trigger for moulting because juvenile-hormone-treated final-instar Bombyx mori and Manduca sexta larvae grow unnaturally large during the instar (see Kamimura & Kiuchi, 1998;Nijhout et al., 2006), suggesting that normal growth does not take cuticle expansion capacity close to the limit.
The ODIM mechanism would shape the growth trajectory and predicts decreasing RMI across instars, which should become pronounced in the last instar(s) (Kivelä, Friberg, et al., 2016). Within the framework of the ODIM hypothesis, we might also predict a positive GDS. This is because GD should be proportional to the mass of tissues that are short of oxygen at the end of an instar, and the mass of oxygen-deficient tissue will increase with increasing total body size across instars. However, there is no empirical evidence that oxygen limitation would become stronger in larger instars (Lundquist et al., 2018; see also Greenlee & Harrison, 2005).
Finally, Callier and Nijhout (2014) hypothesized that compromised availability of metabolic intermediates that act as precursors for growth may result in GD, but it remains unclear how this hypothesized mechanism would shape growth trajectories. Hence, the prediction concerning GD remains vague and may be based on too simplistic assumptions.
The present results were consistent with predictions of the ODIM hypothesis in many respects. RMIS was negative in 8 out of 10 families investigated (Table 3; Figure 2) and GDS was positive in all families except Endromidae (Figure 3; Tables S4 and S5).
However, RMI increased in the last instar despite overall negative RMIS in Drepanidae and Geometridae ( Figures S5 and S7), contradicting the prediction of the ODIM hypothesis. In the Sphingidae, RMIS was positive, which also contradicts the prediction. There are two potential explanations that would reconcile the deviating RMI results within the ODIM framework: hyperallometric tracheal growth across instars and within-instar tracheal growth.
Hyperallometric tracheal growth would increasingly postpone the point at which oxygen demand meets supply across instars, resulting in an across-instar increase in RMI. There is a paucity of data on across-instar tracheal growth, but it appears to be isometric in Manduca sexta (Callier & Nijhout, 2011;Lundquist et al., 2018).
Within-instar tracheal growth would similarly postpone the point at which oxygen demand surpasses supply, and would thus allow a larva to grow more during the instar. Interestingly, in the sphingid M. sexta, tracheal size increases within the final instar (Callier & Nijhout, 2011;Helm & Davidowitz, 2013) but not within earlier instars (Callier & Nijhout, 2011;Lundquist et al., 2018), which could partially explain the presently observed positive RMIS in Sphingidae. Whether the same mechanism could also explain the last-instar increase in RMI in Drepanidae and Geometridae remains unknown. Comparative data on within-and across-instar tracheal growth would be needed for a rigorous assessment of the ODIM hypothesis and mechanisms underlying the observed patterns.
If the balance between oxygen supply and demand has a role in moult induction, then we also need to consider potential cutaneous respiration (see Chapman, 1998) in very small larvae. This is because small larvae have a high surface area to volume ratio, potentially facilitating efficient gas exchange through the cuticle, which would allow more growth within an instar than would be expected considering the tracheal system alone. There is some indirect evidence that cutaneous respiration may occur in first-instar larvae in some species (Kivelä, Friberg, et al., 2016). The growth benefit attainable via cutaneous respiration should decrease rapidly with increasing body size and is expected to be found only in the very first instar(s) in species with very small neonate larvae.
A high RMI in the first instar in relation to the subsequent instars and negative RMIS in Noctuidae, Drepanidae, Endromidae and Geometridae (Table 3; Figure S5) are consistent with a prediction based on the cutaneous respiration hypothesis. However, Endromis versicolora (the only studied endromid species) has very large neonate larvae, suggesting that something else other than cutaneous respiration is responsible for the steep decrease in RMIs across instars in this species. Also, the neonate larvae of the sphingids, Cerura vinula (Notodontidae) and Pararge xiphia (Nymphalidae) are so large that it is unlikely that they would get a significant growth benefit from cutaneous respiration and, consistently, growth trajectories of these species do not show signs of cutaneous respiration ( Figure S7). Our results thus suggest cutaneous respiration in very small larvae, but rigorous physiological studies are needed to demonstrate it.

| CON CLUS IONS
Evolutionary life-history studies cannot ignore developmental and physiological mechanisms that shape growth trajectories because these mechanisms form the proximate basis of traits such as age and size at maturity. Research on the model organisms M. sexta and Drosophila melanogaster has highlighted the evolutionary significance of some proximate mechanisms that operate during the final larval instar in insects (Callier & Nijhout, 2013;Davidowitz et al., 2016;Nijhout et al., 2006Nijhout et al., , 2010Nijhout et al., , 2014Shingleton, 2011). Our work adds to this by showing how comparisons across instars and species can contribute to the understanding of the constraints on the evolution of age and size at maturity. We also demonstrate that there is some variation in growth patterns among evolutionary lineages, which means that caution is required when generalizing results obtained in model organisms.

ACK N OWLED G EM ENTS
We thank Jon Harrison and three anonymous reviewers for constructive comments on an earlier version of the manuscript. We are grate-

DATA AVA I L A B I L I T Y S TAT E M E N T
Phenotypic data deposited in the Dryad Digital Repository http://doi. org/10.5061/dryad.p2ngf 1vmz (Kivelä et al., 2020). DNA sequence data archived in GenBank (accession numbers given in Table S2).