Species and origin of the females used for the breeding experiment are listed in Table 1.

Larvae were kept at about 22 °C in plastic boxes of different sizes, i.e. 5 cm × 3 cm × 2.5 cm for L1 and L2 larvae and 8 cm × 4,.5 cm × 3.5 cm for L3 to L5 larvae. Various Poa species were used as food plant. All shedded head capsules were collected for further examination.

Widths of the head capsules were determined by means of an Olympus VMT-4 stereomicroscope equipped with an ocular micrometer at magnifications of 20× and 80× depending on the size of the head capsule. In some cases, digital images of head capsules in combination with a calibrator were used for measurements by means of IMAGEJ software.

Means and standard deviations of head capsule widths were calculated with EXCEL. GRAPHPAD PRISM was used for calculation of column statistics (t-test) and frequency distributions as well as for linear and nonlinear regression to produce the graphical presentations, i.e. charts, bar diagrams, box plots.

The original data were applied to non-linear regression using an exponential growth equation (y = a * e^bx) plotting the mean of head capsule widths against the instar number. The log-transformed values were analyzed using linear regression (y = ax + b). In addition, a second order polynomial function was applied.

The phylogenetic relationships of the studied Erebia species were analyzed using published sequence data of the mitochondrial cytochrome oxidase subunit I gene (barcodes). A phylogenetic tree was constructed with MEGA7, version 7.0.18 (Kumar et al. 2016) by the neighbor-joining method (Saitou and Nei 1987). For accession numbers of the used COI-5P sequences and for further details see Fig. 22.

Phylogenetic non-independence of growth parameters was checked using PAST4.17 (Hammer et al. 2001). A parsimony-based tree of COI sequences was constructed and phylogenetic generalized least squares (PGLS) was applied to calculate Pagel’s lambda for the growth parameters.

The number of moults and the growth characteristics of the larval head capsules of 28 species of the genus Erebia were examined. Usually, the different species of the genus Erebia undergo metamorphosis either via 4 or via 5 larval instars. As an exception, 6 instars were found in Erebia triaria De Prunner, 1798 (Table 3). In three cases, Erebia aethiops (Esper, [1777]), Erebia neoridas (Boisduval, [1828]) and Erebia medusa ([Denis and Schiffermüller], 1775), development either via 4 or via 5 instars has been observed for specimens from different geographical regions.

For analyses of head capsule growth several parameters were defined and calculated and checked for their suitability for interspecies comparisons. As a basis, means of head capsule widths in µm were determined for all available capsules of each instar of a species. Plotted values revealed that growth progression is not linear but can be described by nonlinear regression. A good fit can be achieved by application of an exponential growth function of the following simple form:

Hx=H0∗ekx (generally: y=a∗ebx) [1]

H_x: head capsule width (µm) of instar x; H₀: extrapolated value for H_x at x = 0; e: basis of natural logarithms (e = 2.718); x: larval instar number; k: growth constant defining the curvature.

After logarithmic transformation of the head capsule values, growth progression can be well described by linear regression resulting in the equation:

log⁡Hx=ax+S( generally: y=ax+b) [2]

a: slope of the line; S: extrapolated value for log H_x at x = 0.

The growth constant k of equation [1] and the slope a of equation [2] turned out to be suitable parameters to describe the size increment of the head capsule. An example of plots and of the equations is shown in Figs 1, 2 for Erebia gorge (Hübner, [1804]) with k = 0.345 and a = 0.155. Goodness of fit values for the two functions amounting to r² = 0.9831 and r² = 0.9932, respectively, indicate that the parameters are suitable to describe the head capsule growth. Furthermore, values of the r² coefficient approaching 1.0 show that the growth increment follows Dyar’s rule.

Another method used by Kim et al. (2016) to mathematically describe the head capsule growth is non-linear regression using a second order polynomial function of the form

y=ax2+bx+c [3]

The resultant graph is additionally shown in Fig. 1 for E. gorge. Goodness of fit value r² = 0.9885 is also very high (see also Table 2). However, comparative interspecies analysis using this method is aggravated by the presence of two varying parameters, a and b, which must be considered together for meaningful comparisons. Some examples of the calculated coefficients are presented in Table 2 which illustrates the difficulties for interpretation and comparative analyses in contrast to the growth factor k of the exponential function.

Dyar’s constant (r) was also calculated for head capsule widths of each pair of successive instars. To account for “variation of the constant” in a growth sequence the values were finally averaged and termed average per-moult growth rate (Minelli and Fusco 2013), which can be considered the averaged Dyar’s coefficient (r_m):

rm=(n−1)−1⋅∑x=2n(Hx−1/Hx) [4]

For E. gorge, with 5 larval instars, the sequential values are H₂/H₁ = 0.634; H₃/H₂ = 0.698; H₄/H₃ = 0.723; H₅/H₄ = 0.737. The values are increasing with higher instar number showing that in this case Dyar’s constant is not strictly a constant. To account for deviations of Dyar’s constant from constancy the standard deviation (SD) of the averaged data pairs (H_x+1/H_x) is calculated. For interspecies comparisons, however, absolute values of SD are less suitable and are therefore transformed into a percentage of the mean. In the case of E. gorge, the mean and SD are 0.698 ± 0.046 and the transformed SD amounts to 6.59% [(0.046 * 100)/0.698)].

Another method to account for deviations from Brooks-Dyar’s ratios was proposed by Hawes (2020), who introduced the following equation [5] to evaluate the geometric property of growth progression (GP):

GP=1/ns⋅∑[(b2−ac)s1+(b2−ac)s2+(b2−ac)s3+……] [5]

corresponding to

GP=(n−2)−1⋅∑x=1n−2(Hx+1)2−(Hx⋅Hx+2) [6]

n = 4, n = 5 etc. for the number of larval instars. H: head capsule width (µm).

The optimal fit, when the growth constant is absolutely constant, results in GP = 0 which is achieved and when the terms (H_x+1)² – (H_x ∙ H_x+2) = 0, i.e., when (H_x+1)² = (H_x ∙ H_x+2) or b² = ac for equation [5] of Hawes (2020), respectively, for all instars. However, a GP of 0 may also occur when growth deviates from Brooks-Dyar’s rule. This is the case when individual terms derive from zero but are positive (b² > ac) or negative (b² < ac) and cancel each other in the sum. It is shown in Fig. 3 that GP values and percentage-transformed SDs of averaged Dyar’s r are not correlated. Thus, GP values may be misleading in interpretation.

Parameters describing head capsule growth such as averaged Dyar’s r (r_m, see Methods) and the growth constant of the exponential growth function k vary over a large range for different Erebia species. Dyar’s r ranges between 0.63 and 0.79 and the growth constant k ranges between 0.23 and 0.44 (Figs 4, 5).

The observed variation in growth parameters appears to be related to the number of moults in different Erebia species. The frequency distribution of the growth constant k reveals 3 distinct clusters comprised of species which develop via 4, 5 or 6 larval instars. A fourth cluster includes species with 4 or 5 instars, which cannot be clearly separated (Fig. 6). The detailed composition of this cluster is depicted in Fig. 7 in which samples are arranged according to increasing k-values and are additionally distinguished by their instar number. Samples with 4 or 5 instars tend to be grouped close to their respective clusters of the frequency distribution. Means and standard deviations of the growth constant k for species with 4, 5, or 6 instars amount to 0.419 ± 0.022 (n = 13), 0.341 ± 0.031 (n = 19), and 0.232 (n = 1), respectively (Tables 3, 4). Corresponding values of the slope (a) derived from the logarithmic functions are 0.181 ± 0.010 (n = 13), 0.151 ± 0.012 (n = 19) and 0.108 (n = 1) (see also Tables 3, 4 below).

The frequency distributions of k-values for species developing via 4 or 5 instars reveal two distinct clusters within each group that differ significantly from each other with p < 0.0001. Means and standard deviations of the growth constant k for species of group 4A, 4B, 5A and 5B are 0.435 ± 0.006, 0.397 ± 0.011, 0.381 ± 0.004 and 0.323 ± 0.016, respectively.

The determined growth parameters are summarized for the different Erebia species in Tables 3, 4 in which the species are arranged according to the clusters defined above. For comparison data for some Satyrinae species belonging to tribes other than Erebiini are shown in Table 5. In addition, Table 6 presents calculated parameters of published data including results from three North American Erebia species and from two non-Satyrinae species.

The averaged growth parameters for the recognized groups (4, 5, 6) and subgroups (4A, 4B, 5A, 5B) as listed in Tables 3, 4 are presented as box plots in Figs 10–12 which demonstrate that the groups are convincingly distinguishable by the growth parameters k and a of the mathematical functions. They differ significantly (P < 0.0001), except for the subgroup pair 4B/5A with lower significance (P = 0.0052 and P = 0.0121 for the growth parameters k and a, respectively; Table 7). Dyar’s r also differs significantly between all groups except for the subgroup pair 4B/5A. A box plot for the calculated start values H₀ of the exponential growth function for the different clusters is presented in Fig. 13. A corresponding significance analysis does not reveal any relevant differences between the groups or subgroups (Table 7).

There is a good correlation between Dyar’s r and the growth constant k (r² = 0.9674) depicted in Fig. 14 in which the different subgroups are also distinguished. All the differences in growth parameters between the groups are highly significant (P < 0.0001). This is also true for Dyar’s r of the groups with 4 and 5 instars.

The head capsule size of a larval stage varies among individuals. However, the ranges of variation between successive larval stages rarely overlap. Therefore, in most cases, it is possible to assign a larva to a defined stage based on head size. This is exemplified in Figs 15, 16, which show the size distribution of larval head capsules as box plots for Erebia gorge and Erebia aethiops. Similar results were obtained for the other Erebia species studied.

For a few species, variability in growth parameters was examined for samples from different localities or populations, including cases in which the same species can develop via 4 or 5 larval instars. In Figs 17, 18 the growth functions of head capsule widths for Erebia aethiops from Hautes Alpes, Isère, and Bavaria are shown. Growth characteristics of larva from Isère and Bavaria, which develop via 4 instars are similar and clearly differ from the Hautes Alpes specimens (5 instars). Regardless of the developmental type, however, the head capsule widths of the final instars (L4 or L5), were nearly identical (Figs 17, 18). As the initial head capsule widths in L1 are similar, the ratios of the head capsule widths of the final and the first instar are also identical: namely 3.62 (Hautes Alpes), 3.62 (Isère) and 3.61 (Bavaria). Thus, the same final result is achieved regardless of developmental type. Similar results were obtained for Erebia medusa (not shown) and also for Erebia neoridas (Figs 19, 20) with H_t/H₀ ratios of 3.83 and 3.54, respectively.

For most of the studied species their life cycle is known. Respective data were well summarized by Lafranchis (2000) and by Sonderegger (2005) for the French and the Swiss Erebia species, respectively. The arrangement of Sonderegger (2005) was adopted here and related to the development groups defined above. The arrangement considers larval developmental time (1 or 2 years) as well as the stage of hibernation, and has been coded using numerical values according to Table 8. The assigned values were then plotted against the growth constant k of the exponential function (Fig. 21).

The plot shows that species with 4 larval instars (4A and 4B) develop within a 1-year life cycle with one exception only (Erebia lefebvrei), while species with 5 instars of group 5B may develop over one or two years. All five species of group 5A demand a 2-year life cycle.

To test the hypothesis that the various growth strategies are associated with phylogenetically based species clusters of the Erebia radiation several analyses were performed. First, a neighbor-joining tree and a maximum parsimony tree were constructed for the species examined in this study. They are based on COI DNA-sequence sections of the 5'-region with a length of 658 base pairs accessible from GenBank (National Center for Biotechnology Information). For Erebia melancholica and the outgroup species Callerebia polyphemus (Oberthür, [1876]), slightly shorter sequence sections were available only.

Fig. 22 shows the resultant neighbor joining tree in which the species have been marked by colored circles and diamonds according to the growth types described in Tables 3–6. Additionally, four species clusters were marked with A – D, including the species of the pronoe-group (A) as defined by Warren (1936) and the species of the tyndarus-group (C), see Lorkovic (1957, 1958). They are considered monophyletic for the following reasons. They consistently appear in trees based on sequence data of further mitochondrial and nuclear genes, namely ND1, ND5, Wingless, GAPDH, RpS5, and 16SrRNA as well as on restriction data of ND5 (Roos, unpublished). The 4 clusters are also identified in the recently published multi-gene phylogeny of European butterflies (Wiemers et al. 2020). In addition, their monophyly is supported by imaginal and preimaginal morphological characters (Roos, unpublished).

Obviously, there is no association between growth strategies and phylogenetically defined species clusters, in particular pointed up by the monophyletic melas-group (A) and tyndarus-group (C).

In a second approach, phylogenetic non-independence of all growth parameters of Tables 3–6 was checked by PGLS (phylogenetic generalized least squares). A parsimony-based tree was used with Proterebia afra as outgroup. Pagel’s lambda was calculated for all growth parameters. The low lambda-values, most are zero, show that growth parameters are independent of the phylogenetic background (Table 9).

As is evident from Table 10 the ratio of head capsule widths (H) from successive instars (Dyar’s r) of a species is not constant in most cases and thus does not follow Dyar’s rule. The deviation from constancy can be expressed by the calculated standard deviation (SD) of averaged Dyar’s r which is further transformed into % deviation for better inter-species comparability. For some data sets the deviation is below 2% such as for E. medusa, E. cassioides and E. melancholica so that constancy of the size increment, i.e. consistency with Dyar’s rule, can be assumed. High deviations from the rule, however, were observed for E. hispania, E. alberganus, developing via 4 instars, and E. melas, E. neoridas, E. aethiops developing via 5 instars.