Evolutionary fitness is determined by the match between an organism’s phenotype and its local environment. When mismatched, individuals may disperse to more suitable habitats. For flightless insects, however, the range of dispersal is typically limited. Numerous flightless species have, therefore, evolved a dispersal dimorphism, that is, some individuals in otherwise short-winged populations develop long wings. Wing development may be genetically or environmentally determined, but these two drivers have rarely been analysed together.

We studied the inheritance and density-dependent plasticity in the dispersal dimorphism of the meadow grasshopper Pseudochorthippus parallelus. Using a full-sib half-sib breeding design, we found that the development of long wings strongly depends on rearing density, with tactile stimulation being the most likely proximate cause. Additionally, we found heritable variation in the development of long wings, both in the propensity to produce long wings and in response to density (genotype-by-environment interactions). While at high and low densities, the environmental effect dominates, genetic variation takes effect mostly at intermediate densities.

Our results have implications for the phenotype-environment match and ultimately the evolution of individualized niches. Induced dimorphisms represent a form of niche conformance (or adaptive phenotypic plasticity) and both genetic and induced dispersal dimorphisms facilitate niche choice in allowing individuals to sample a greater range of environments. Our study shows that niche-related polymorphisms can evolve directly and indirectly via selection on the sensitivity threshold.

Parental generation

We set up a full-sib half-sib breeding design with meadow grasshoppers Pseudochorthippus parallelus. Parental individuals were caught from the field in the surroundings of Jena (50.94°N, 11.61°E) in June and July 2021. Only nymphal females were caught to ensure virginity. Males were mostly caught as nymphs, but a few adult males were also sampled. All parental individuals were short-winged, except for two long-winged males (0.8% of all parents). Parental individuals were housed in the laboratory in same-sex cages until maturity. A total of 172 mature females were then assigned to individual mating cages (22 x 16 x 16 cm), where they were maintained with freshly cut grass potted in small vials filled with water and a water tube for moisture. Grass pots were replaced three times per week. A small pot containing a 50:50 vermiculite-sand mixture was provided for egg deposition. A total of 68 males were mated to 1-4 females each by moving males every 2-3 days between mating cages. Males were kept in rotation and females were allowed to lay eggs until they died (or when breeding was terminated in early September 2021).

Sand pots were sieved once per week for the collection of egg pods. Egg pods are solid structures of 1-2 cm in length containing up to 12 eggs, and are typically buried in the sand. Egg pods were retrieved and placed in petri dishes lined with moist filter paper, while carefully documenting the cage of origin. All egg pods from the same cage and collection date were gathered on the same dish. A total of 1,560 egg pods were collected from 167 parental cages. In October, we replaced the filter paper with a 50:50 vermiculite-sand mixture and transferred petri dishes to refrigerators (approx. 5°C) for diapause. Egg pods were sprayed twice with a fungicide to prevent fungal growth, and sprayed with water weekly before diapause and biweekly during diapause.

Offspring generation

A single offspring generation was raised to adulthood in five cohorts (cohorts represented the same generation, but had to be raised in separate time spans due to constraints in housing capacities). Individuals that hatched from the same petri dish were released together into a single offspring cage (22 x 16 x 16 cm). As different petri dishes produced different numbers of offspring, the densities per cage were the result of the hatching success in each petri dish. However, some densities were experimentally manipulated. After the large hatching peak had passed, cages hosting more than 8 offspring were equally split into two cages (hereinafter referred as receiving cages as opposed to the source cage from which they originated). For logistic reasons, 32 cages (3% of the cages with more than 8 offspring) were not split and remained with densities ranging from 9 to 12 offspring per cage. Due to the large number of cages involved, it was operationally unfeasible to experimentally manipulate all offspring cage densities. Therefore, cages with up to 8 offspring were always left unsplit.

Some egg pods unexpectedly hatched before diapause on October 2021. Pre-diapause hatching came at a surprise and had not yet been documented for the meadow grasshopper. Offspring that hatched before diapause were raised as cohort 1. For this cohort, all egg pods collected from a parental cage on the same day were in a single petri dish. This happened because individual egg pods were isolated only after diapause to minimize the number of dishes to be tended during diapause. Therefore, for cohort 1, it is uncertain if hatchlings in an offspring cage came from the same or different egg pods of the same female. After diapause, each egg pod was placed in a Petri dish. For cohorts 2 to 5 (hatched in April, May, June, and July 2022, respectively), the offspring raised in one cage came from the same egg pod (with few exceptions). We emphasize that cohorts 1-5 represent a single offspring generation, and the separation into cohorts was done solely to produce a manageable number of active cages (at times, more than 400 offspring cages were maintained simultaneously).

Nymphs were maintained with ad libitum access to freshly cut grass potted in small vials filled with water and a water tube for moisture. Dead nymphs were recorded and removed. In total, 6,598 nymphs hatched, 1,683 (25%) nymphs died, 401 (6%) nymphs were lost for other reasons and 4,436 (67%) reached the imago stage about four weeks after hatching. Average full-sib family size of successfully developing offspring was 30.9 ± 20.9 individuals (mean ± SD, range 1-95) with 7.9 ± 4.2 (range 1-18) cages per family and 4.6 ± 1.9 (range 1-10) nymphal density classes represented per family (see below and Table 1 for offspring sample sizes).

Offspring densities were quantified as the number of offspring per cage. We refer to the number of hatchlings that were released into a cage as hatchling density and to the number of adults that developed in a cage (after early nymphal deaths and splitting of cages with more than 8 individuals) as nymphal density. Most nymphal mortality occurred within the first few days after hatching, such that hatchling density roughly represents density within the first week of nymphal development, while nymphal density mainly represents the density during weeks 2-4 of nymphal development (in which little mortality occurred). Nymphal density was analysed by including both source (unsplit) cages and receiving cages (with densities after splitting). Hatchling density did not include receiving cages (9% of all cages) and included source cages with their densities before splitting. Hatchling density, in principle, also represents the number of embryos that successfully emerged from a single egg pod. However, offspring from cohort 1 and 4.5% of the cages from cohorts 2 to 5 came from more than one egg pod and for those cases the number of developing embryos per egg pod was unknown. Since embryos may influence each other during development, we analysed embryo density as the number of hatchlings among the subset of cages that derived from a single egg pod.

These densities could be confounded with maternal fecundity. However, most of the fecundity of female grasshoppers is determined by the total number of egg pods produced rather than the number of fertilized eggs per egg pod (which is also, perhaps mostly, affected by male fertility (Butlin et al., 1987). The total number of offspring per female was correlated with offspring density per cage in our experiment, though the correlation was not strong (r = 0.33). While we cannot rule out an effect of maternal fecundity on density, we consider it more likely that any effects operated via embryo number (within egg pods) or nymphal density.

Offspring phenotyping

After the final moult, mature offspring were collected from their offspring cages. We recorded offspring cage, sex and wing morph. Short-winged meadow grasshoppers have their hind wings rudimentarily developed. Fore wings of short-winged females reach the middle of the abdomen (Figure 1). Short-winged males have better developed forewings that can reach the tip of the abdomen. Fore and hind wings exceed the abdomen tip in long-winged males while almost reaching the abdomen tip in long-winged females. In some intermediate cases (2.5%), we checked the length of the hind wings and assigned to long-winged individuals the ones with well-developed hind wings.

General statistical analyses

Before delving into the animal model analyses (see below), we estimated sex-biases in the probability of developing short or long wings. This was done be χ² tests on pooled data. To further evaluate sex-differences in the wing length development across the environmental gradient, we used correlation analysis for the proportion of short-winged individuals for each density class on pooled data (females vs. males). Since there was no evidence of any sex-specific effect even in these analyses on pooled data (see results), we did not include offspring sex in the animal model analysis.

Furthermore, we evaluated the different measures of density as described above (nymphal density, hatchling density, embryo density) and also proportional survival per cage before fitting the animal model. This was done on pooled data in a generalized linear model with binomial error structure, logit link and density as the only predictor.

Animal model analysis

We fitted animal models (Kruuk, 2004) to decompose the variance into environmental and genetic sources of variation using the R package brms version 2.20.4 (Bürkner, 2017). Our basic animal model is represented by the following phenotypic equation and population-level variances:

 y_i ~ B(1,p_i)

logit(p_i) = α + βx_i + c_i + a_i + e_i

c_i ~ N(0, V_C)

a_i ~ N(0, V_A)

e_i ~ N(0, V_R)

where y_i represents an indicator of wing length of offspring i (y_i = 1 for long-winged and y_i = 0 for short-winged individuals) that is modelled by a Binomial distribution with probability p_i for offspring i (note that parameter n = 1 for binary data). Probability p_i is modelled on the logit scale where α represents the intercept, x_i represents the density experienced by individual i (centred to zero at a density of 6.5, which represents the midpoint between densities 1 and 12), β represents the population-level slope of the density effect, c_i represents the cage identity (permanent environment) random effect of individual i with a population-level variance , a_i represents the additive genetic effect (individual identity linked to the two-generation pedigree) with a population-level additive-genetic variance  and e_i represents the residual effect for individual i with a population-level residual variance . The additive genetic relatedness matrix used to estimate  consisted of values 0 for unrelated individuals, 0.5 for full-siblings and 0.25 for half-siblings. We refer to the variance  as the pedigree effect here (sometimes also referred to as the animal effect).

Furthermore, we fitted a random-slope animal model that estimates the additive-genetic variance of response to nymphal density (as well as the random intercept-slope correlation). The random-slope animal model is represented by the following phenotypic equation and population-level variances:

y_i ~B(1, p_i)

logit(p_i) = α + (β + a_i,2)x_i + c_i + a_i,1 + e_i

c_i ~ N(0, V_C)

(a_i,1   a_i,2) ~ MVN( (0   0), variance-covairance matrix of V_A,I and V_A,S)

e_i ~ N(0, V_R)

where a_i,1  now represents the additive genetic effect with the population-level additive-genetic random intercept variance V_A,I , a_i,2   represents the additive genetic effect for the response to density with the population-level additive-genetic random slope variance V_A,S , and C_A,IS  represents the population-level additive genetic intercept-slope covariance. Other terms are the same as above.

We used default minimally informative priors for all parameters estimated by the model. The (non-)informativeness of the priors was checked with posterior predictive checks. For the fixed effect slope β the default prior is an improper flat prior. For random effects variances , ,  and  default priors are half student-t priors on the respective standard deviation with 3 degrees of freedom and a scale parameter adjusted to the standard deviation of the response on the link scale (http://paul-buerkner.github.io/brms/reference/set_prior.html). Finally, for the additive genetic variance-covariance matrix there is a prior for the Cholesky factor of the correlation matrix with η = 1. We ran two chains of the Bayesian model with a warmup of 50,000 iterations followed by 120,000 iterations for sampling from the posterior distributions keeping every 70^th sample. Convergence was checked from trace plots, by  criterion (all R_hat  ≤ 1.008) and posterior effective sample size (all ESS ≥ 1,774, out of 2,000 posterior samples across the two chains). All analyses were done in R 4.1.3 (R Core Team, 2022).

Evaluation of the statistical significance of variance components can be difficult in Bayesian models, since variances are bound at zero and posterior medians and credibility intervals are often larger than zero even in the absence of a statistically significant effect (Pick et a. 2024).  Since the statistical significance of the G x E effect (estimated in the form of random-slope variation) is of particular interest in our analysis, we used simulations to evaluate the posterior distribution in comparison to a null-model reference scenario. This was done by randomizing the vector of densities and fitting this randomized vector as an additional fixed effect and random-slope term – instead of the random slope term for the actual density data. All other terms in the model remained the same, that is, we fitted a population-level slope for actual density data and individual identity linked to the pedigree as a random (intercept) effect. We compared the posterior distribution for the random slope variance for the real-data model to the posterior distribution for the randomized random-slope term for 19 iterations of randomization.

Variance components were calculated following (Anonymous, 2010) with the distribution-specific variance of the binomial distribution approximated by π2/3 and the residual variance fixed to 1. These variances give the link-scale variances. We can also calculate the data-level variances (following  Anonymous, 2010), Table 1, using the average proportion of long-winged individuals for P), which, reassuringly, gave highly similar variance ratios. For random-slope models, however, variance components depend on the value of the covariate and we used a recently proposed method for conditional variance decomposition  Anonymous, 2022a). The total (marginalized) additive genetic variance (arising from random intercept and random-slope additive genetic variation) is given by V_A,M = V_A,I + V_A,SV_x + µ²V_A,S + 2 µ C_A,IS, where V_A,I, V_A,S and C_A,IS** are** as above (for the random-slope model), and µ and V_x the mean and variance of nymphal densities, respectively. The variance exclusively explained by random-slopes is V_s = V_A,SV_x + µ²V*_A,S. Both, V_A,M and V_S can be set in relation to the total phenotypic variance.