DOI: 10.5061/dryad.18931zd4t
title: Phenotypic plasticity, heritability, and genotype-by-environment interactions in an insect dispersal polymorphism

description: 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.
description: <em>Parental
 generation</em> We set up a full-sib half-sib
 breeding design with meadow grasshoppers <em>Pseudochorthippus
 parallelus</em>. 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. <em>Offspring
 generation</em> 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 <em>receiving
 cages</em> as opposed to the <em>source cage</em> 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
 <em>ad libitum</em> 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 <em>hatchling density</em> 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 <em>nymphal
 density</em>. 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 <em>embryo
 density</em> 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. <em>Offspring
 phenotyping</em> 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. 
 <em>General
 statistical analyses</em> 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
 χ<sup>2</sup> 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<em><sub>i</sub></em> represents an indicator of
 wing length of offspring <em>i</em>
 (y<em><sub>i</sub></em> = 1 for long-winged
 and y<em><sub>i</sub></em> = 0 for short-winged
 individuals) that is modelled by a Binomial distribution with probability
 <em>p</em><sub>i</sub> for offspring
 <em>i</em> (note that parameter n = 1 for binary data).
 Probability <em>p</em><sub>i</sub> is modelled on
 the logit scale where α represents the intercept,
 x<em><sub>i</sub></em> represents the density
 experienced by individual
 <em><strong>i</strong></em> (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<em><sub>i</sub></em> represents the cage
 identity (permanent environment) random effect of individual
 <em><strong>i</strong></em> with a
 population-level variance , a<em><sub>i</sub></em>
 represents the additive genetic effect (individual identity linked to the
 two-generation pedigree) with a population-level additive-genetic variance
  and e<em><sub>i</sub></em> represents the
 residual effect for individual
 <em><strong>i</strong></em> 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<sup>th</sup> 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 π<em>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 <em>P</em>), 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 <em>V</em><sub>A,M</sub> = <em>V</em><sub>A,I</sub> + <em>V</em><sub>A,S</sub></em>V<em><sub>x</sub> + µ²V<sub>A,S</sub> + 2 µ C<sub>A,IS</sub>, where <em>V</em><sub>A,I</sub>, <em>V</em><sub>A,S</sub> and C<sub>A,IS</sub>** are** as above (for the random-slope model), and µ and <em>V</em><sub>x</sub> the mean and variance of nymphal densities, respectively. The variance exclusively explained by random-slopes is <em>V</em><sub>s</sub> = <em>V</em><sub>A,S</sub></em>V<em><sub>x</sub> + µ²</em>V*<sub>A,S</sub>. Both, <em>V</em><sub>A,M</sub> and <em>V</em><sub>S</sub> can be set in relation to the total phenotypic variance.
description: # Phenotypic plasticity, heritability, and genotype-by-environment
 interactions in an insect dispersal polymorphism
 [https://doi.org/10.5061/dryad.18931zd4t](https://doi.org/10.5061/dryad.18931zd4t) ## Description of the data and file structure We used a full-sib half-sib breeding design to assess the heritability and the phenotypic plasticity of the wing dimorphism of the meadow grasshopper *Pseudochorthippus parallelus*. We found that the development of long wings strongly depends on rearing density. Additionally, we found significant heritable variation in the development of long wings. Code version: R version 4.2.2 **Overview of files:** \- **WingLength_ParentsData_240405.txt** = data relative to the generation used as parents in our full-sib half-sib breeding design. Column names: CageID = a single identifying code for each cage housing a mating pair FemaleID = a single identifying code for each female (Uxxx: from generation U, xxx being the single identifying code with xxx &lt; 500 for females) MaleID = a single identifying code for each male (Uxxx: from generation U, xxx being the single identifying code with xxx &gt; 500 for males) MatingWingLengthFemale = wing length morph of the mating female (Short: short-winged female, NA: do not apply) MatingWingLengthMale = wing length morph of the mating male (Short: short-winged male, Long: long-winged male, NA: do not apply) \- **WingLength_OffspringData_240405.txt** = data relative to the offspring generation. Each line corresponds to a single offspring born from the mating pairs of the full-sib half-sib breeding design. Column names: IndID = a single identifying code for each offspring (VPAxxx, xxx being the single identifying code) Sex = sex of the offspring Cohort = the cohort from which the offspring has been raised (cohort 1 hatched before diapause in October 2021, and cohorts 2 to 5 hatched in April, May, June, and July 2022, respectively) WingLength = wing length morph of the offspring (Short: short-winged individual, Long: long-winged individual, NA) WingIntermediate = individuals with ambiguous wing morph (x) as opposed to clearly distinguishable wing morph (.) InconsistentWingLength = individuals for which the re-scoring of wing morph was inconsistent with the first scoring (x) as opposed to consistent re-scoring (.) CageIDOrg = a single identifying code for each cage originally housing the offspring at hatching (referred as source cage in the manuscript; VPxxxx, xxxx being the single identifying code) CageID = a single identifying code for each cage housing the offspring after splitting the cages with more than 8 offspring (referred as receiving cage in the manuscript; VPxxxx, xxxx being the single identifying code) nIndTotalIn = total number of offspring raised in the source cage of the focal individual (corresponds to Embryo density if all offspring came from a single egg pod, and to Hatching density if offspring came from several egg pods [see nEggPods]) nEggPods = number of egg pods from which offspring from the same source cage hatched (from 0 to 3 egg pods, NA when unknown, x when coming from the 1st cohort [unknown as well]) GrSize = total number of offspring raised in the receiving cage of the focal individual (corresponds to Nymphal density) nNymphsOut = number of offspring that left the source cage to form the receiving cage InFrom = this column tells FROM which SOURCE cage the focal individual originated when it has been moved into a receiving cage (individuals that stayed in their source cage are labelled with ".") OutTo = this column tells TO which RECEIVING cage the focal individual has been moved into (individuals that stayed in their source cage are labelled with ".") FemaleID = the single identifying code corresponding to the mother (Uxxx: from generation U, xxx being the single identifying code with xxx &lt; 500 for females) MaleID = the single identifying code corresponding to the father (Uxxx: from generation U, xxx being the single identifying code with xxx &gt; 500 for males) WingLengthFemale = wing length morph of the mother (Short: short-winged female, NA) WingLengthMale = wing length morph of the father (Short: short-winged male, Long: long-winged male, NA) 
 
 
 \- **WingLength_241220.Rmd** = Rmarkdown script with the complete
 analyses and figures 
 
 
 \- **WingLength_241220.html** = script in publishable format with the
 complete analyses and figures outputs 
 
 
 \- **animalModel_*xxxxx*.RData** = backups of the outputs of the animal
 models run in the script, '***xxxxx***' describes the type of
 animal model that was run (these files can be directly loaded from the
 Rmarkdown script) 
 
 ## Code/Software Software: R Code version: R version 4.2.2