Phenotypic integration plasticity in Daphnia magna: an integral facet of G × E interactions

Phenotypic integration can be defined as the network of multivariate relationships among behavioural, physiological and morphological traits that describe the organism. Phenotypic integration plasticity refers to the change in patterns of phenotypic integration across environments or ontogeny. Because studies of phenotypic plasticity have predominantly focussed on single traits, a G × E interaction is typically perceived as differences in the magnitude of trait expression across two or more environments. However, many plastic responses involve coordinated responses in multiple traits, raising the possibility that relative differences in trait expression in different environments are an important, but often overlooked, source of G × E interaction. Here, we use phenotypic change vectors to statistically compare the multivariate life‐history plasticity of six Daphnia magna clones collected from four disparate European populations. Differences in the magnitude of plastic responses were statistically distinguishable for two of the six clones studied. However, differences in phenotypic integration plasticity were statistically distinguishable for all six of the clones studied, suggesting that phenotypic integration plasticity is an important component of G × E interactions that may be missed unless appropriate multivariate analyses are used.


Introduction
The concept of phenotypic plasticity has been framed in terms of properties specific to individual characters rather than the whole organism (Bradshaw, 1965). This, combined with the fact that single-trait studies and models are often much simpler and easier to obtain than multivariate ones, has led to an understanding of phenotypic plasticity predominantly focussed on single traits (Schlichting, 1989a;Pigliucci, 2004;Relyea, 2004). However, organisms are integrated entities characterized by numerous correlated traits that are linked by genetic, developmental, physiological or functional associations (Schlichting, 1989a;Roff, 2002;Pigliucci & Preston, 2004). Hence, there are potentially a number of problems with univariate studies. First, if phenotypic plasticity in different traits is not correlated, then the relative plasticity of different genotypes may depend entirely upon which trait is being studied (Schlichting, 1989a). Second, it seems unlikely that the processes proposed to explain the limitations of a plastic response operate at the level of individual traits (Steiner & Van Buskirk, 2008). Finally, univariate measures of plasticity do not capture differences in plastic responses that involve coordinated change in multiple traits (Schlichting, 1989a,b;Schlichting & Piglucci, 1998).
A significant G 9 E interaction in a multivariate analysis of variance (MANOVA) can be used to test for significant genetic difference in multivariate phenotypic plasticity in much the same way as it is currently used in univariate analyses (Langerhans & DeWitt, 2002). However, identifying the most plastic genotype is complicated by the fact that when we are dealing with multivariate phenotypes, a significant G 9 E interaction can arise in two ways (Collyer & Adams, 2007;Adams & Collyer, 2009). Genotypes might differ in the magnitude of phenotypic change, with different genotypes showing the same response to an environmental change but to a greater or lesser degree. Alternatively, genotypes might differ in the nature of their phenotypic change, with different genotypes varying in their response to the same environmental change through differences in their phenotypic integration plasticity (Schlichting, 1989a,b). Examples of adaptive phenotypic integration plasticity include shade avoidance in plants (Schmitt et al., 1999;Weinig, 2000;Sultan, 2003;Pigliucci, 2004) or predator avoidance in animals (Lively, 1986;Boersma et al., 1998;Donohue et al., 2000;Van Buskirk, 2002;Relyea, 2004;Dennis et al., 2011) which both involve coordinated shifts in the expression of many different traits. However, phenotypic integration plasticity may also be nonadaptive (Schlichting, 1989a,b;Schlichting & Piglucci, 1998;Pigliucci & Preston, 2004).
Because studies of phenotypic plasticity are typically univariate (Schlichting, 1989a;Relyea, 2004), it is currently not clear to what extent multivariate G 9 E interactions are a product of differences in the magnitude of phenotypic change, differences in phenotypic integration plasticity or both. For example, Collyer and Adams (2007) found that phenotypic integration plasticity explained character displacement in the head shape of sympatric and allopatric populations of two Plethodon salamander species, whereas differences in the degree of sexual size dimorphism in White Sands pupfish in different populations were explained by differences in the magnitude of phenotypic response. In contrast, Chun et al. (2007) found that differences in the multivariate plasticity of native and invasive populations of the plant purple loosestrife, Lythrum salicaria varied according to the stressor.
Understanding what causes multivariate G 9 E's is important for a number of reasons. First, if phenotypic integration plasticity is a major component of a multivariate G 9 E interaction, univariate studies that are not able to detect it are clearly not suitable for comparing plastic responses between different genotypes. Second, acknowledging that a G 9 E interaction may be the product of differences in the magnitude of phenotypic change and/or differences in phenotypic integration plasticity may facilitate a better understanding of the mechanisms underpinning plastic responses and the selection pressures operating on them (Chun et al., 2007). Finally, quantifying the extent that phenotypic integration plasticity underpins G 9 E interactions may help to explain how the environment can often induce changes in the strength or direction of genetic correlations between traits (Pigliucci & Preston, 2004;Townley & Ezard, 2013).
The 'two-state multivariate approach' (Collyer & Adams, 2007;Adams & Collyer, 2009) is a general statistical framework for comparing multivariate reaction norms, summarizing differences in the magnitude of phenotypic change as differences in the length of phenotypic change vectors, and differences in phenotypic integration plasticity as differences in the angle between phenotypic change vectors in multivariate trait space. In this study, we use traditional univariate and multivariate analyses (MANOVA, PCA) and the two-state multivariate approach to compare the multivariate lifehistory plasticity in response to resource availability for six Daphnia magna clones collected from four different European populations. We test whether Daphnia magna demonstrate multivariate G 9 E interactions, and if so, whether these derive from differences in the magnitude of phenotypic change, differences in phenotypic integration plasticity or both.

Origin and maintenance of clones
The clones used in this study were collected from various disparate European populations. Clones B5 and B7 originated from Weston Park, Sheffield (53°38 0 20″N 1°49 0 07″W). Clones DKN 1-3 and DKN 1-6 came from Kniphagen, Ostholstein, Germany (54°10 0 36 0 86″N 10°48 0 24 0 06″E), and clone H01 came from Bogarzo-to, Kikungsagi-nemzeti park, Hungary (46°48″N 19°08″E). Finally, clone Grosse originated from Kreis Pl€ on, Schleswig-Holstein, Germany (54°19 0 29″N 10°37 0 39″E). All clones were maintained as laboratory stock cultures in a temperature-controlled incubator at 21°C AE 1°on a 14 : 10 L/D cycle. Prior to experimentation, all clones were conditioned for at least 3 generations in order to control for maternal effects. Three individuals from each clone line were put into separate 175 mL jars containing 150 mL of hard water ASTM enriched with an organic extract (Marinure) (see Baird et al., 1989 for details). All individuals were fed 200 000 cells of batch-cultured Chlorella vulgaris algae each day (High food), and the media in each jar was changed every other day. On producing the 3rd clutch, a single offspring was collected and used to set up the subsequent generation.

Measuring life histories
For each clone, thirty-third clutch offspring were randomly allocated into either high (200 000 cells mL À1 ) or low (40 000 cells mL À1 ) food treatments. All individuals were photographed as neonates using a Canon EOS 350D digital camera connected to a Leica MZ6 dissecting microscope at 9 2.5 magnification. All individuals were fed once a day and their developmental stage was recorded. Individuals that had moulted were photographed as described above. Body size was measured from the ventral base of the tail spine to the anterior edge of the carapace using the ImageJ 1.28u image analysis package (Rasband, 1997(Rasband, -2008. Individuals were recorded as being mature once eggs were observed in the brood pouch. For each clutch that females produced, the number of offspring was counted and five neonates photographed in order to determine the mean body size of the clutch. This procedure was continued until all individuals had produced their 3rd clutch. The seven life-history traits measured for each individual therefore included neonate size (mm), size at maturity (mm), age at maturity (days), size upon dropping the 3rd clutch (mm), age upon dropping the 3rd clutch (days), mean fecundity per clutch (numbers) and mean neonate size (mm).

Statistical analyses
Univariate G 9 E (clone 9 food) interactions were tested for using ANOVA. Subsequently, a difference in the multivariate plasticity exhibited by different clones in response to food availability (G 9 E) was tested for using a multivariate analysis of variance (MANOVA). Differences between different clones responses were visualized by projecting multivariate clone 9 food means onto the first two principal component axes of the multivariate set of trait values and clone 9 food effects (Chun et al., 2007). For pairwise comparisons between clones, a two-state multivariate analysis was used to describe the multivariate reaction norms as phenotypic change vectors, following the methodology proposed by Collyer and Adams (2007). These phenotypic change vectors are simply the difference in multivariate mean vectors: D Y ¼ Y ij À Y ik , for treatment group i in environments j and k. The magnitude of the vector is calculated as the Euclidean distance of a phenotypic change vector: where T represents a vector transpose. To determine whether the phenotypic change vector of one clone is greater than another, the test statistic |D E1 À D E2 | is calculated (see Adams & Rohlf, 2000). To test whether the clones differ in the angle of their phenotypic change vectors, their vector correlation is calculated as the inner product of the two vectors scaled to unit length: (Cheverud & Leamy, 1985;Schluter, 1996). The arccosine of this value is the angle, h, between vectors, which describes the similarity of their direction (i.e. h°i s equivalent to parallel vectors). The significance of the two test statistics can then be evaluated by comparing them to distributions created from random pairs of vectors created using a permutation procedure (see

Principal component analysis (PCA)
PC1 and PC2 together explained 77% of the total variation in the data set (Table 1; Fig. 2a). PC1 explained more than 45% of the total variation and included clonal differences in PC1 scores irrespective of environment and the plastic response of clones with respect to food availability (see Fig. 2b). Individuals reared in a high food environment had positive PC1 scores that were associated with more offspring and larger body size after dropping the third clutch of offspring (see Table 1, Fig. 2). Individuals reared in a low food environment had negative PC1 scores that were associated with delayed age at maturity, delayed age after dropping the third clutch of offspring and a larger mean offspring size (see Table 1, Fig. 2). In contrast, PC2 explained 32% of the total variation and varied between clones irrespective of which food environment they were in. Thus, clones showed considerable variation in the mean size of offspring produced, neonate length and size at maturity that was unrelated to a plastic response to food environment, although four of the six clones showed a tendency to produce larger offspring in low food environments (see Fig. 1).

A comparison of phenotypic change vectors
A comparison of phenotypic change vectors for each clone revealed that the magnitude of phenotypic change (i.e. the length of the phenotypic change vector) was indistinguishable for clones DKN16, DKN13, B5 and B7, but significantly longer for clone 'H01' and even longer for clone 'Grosse' (Table 2). However, the direction of the phenotypic change vectors was significantly different for all clone comparisons ( Table 2), demonstrating that differences in phenotypic integration plasticity underpin the multivariate G 9 E interaction observed in this study.

Discussion
G 9 E interactions have traditionally been viewed as differences in the magnitude of the plastic responses of different genotypes (Bradshaw, 1965;Schlichting & Piglucci, 1998;Nussey et al., 2007). However, we show here that when a multivariate approach to phenotypic plasticity is adopted, phenotypic integration plasticity is an important component of a G 9 E interaction that, in this case at least, may even overshadow differences in the magnitude of phenotypic change. This is an important finding given that phenotypic integration plasticity would probably not have even been detectable in a univariate study and may have been underestimated by previous techniques used such as PCA, which sometime fails to distinguish differences in plasticity from differences in trait means. PCA also provides no statistical inferences with respect to differences between clone responses. The two-state multivariate approach avoids this problem by accounting for fixed effects prior to analysing G 9 E interaction effects (Collyer & Adams, 2007). For these reasons, it is a useful tool for comparing multivariate plasticity phenotypes.
Incorporating phenotypic integration plasticity into studies of phenotypic plasticity is important for a number of reasons. First, being able to reliably quantify phenotypic plasticity is an essential precursor to any study attempting to measure a cost of plasticity (Relyea, 2002;Callahan et al., 2008). As most studies to date have been univariate and have therefore not incorporated phenotypic integration plasticity (but see Steiner & Van Buskirk, 2008), it may not be surprising that evidence for the existence of costs or limits of plasticity is currently at best equivocal (DeWitt, 1998;Scheiner & Berrigan, 1998;Donohue et al., 2000;Steiner & Van Buskirk, 2008;Callahan et al., 2008). Second, acknowledging phenotypic integration plasticity is imperative in shaping our understanding of how phenotypic plasticity evolves (Schlichting, 1989a,b;Schlichting & Piglucci, 1998;Ernande et al., 2004;Pigliucci & Preston, 2004;Tonsor & Scheiner, 2007;Gianoli & Palacio-Lopez, 2009). In any environment, selection should favour a suite of traits that result in a functional phenotype (Arnold, 1983;Schlichting, 1986Schlichting, , 1989aSchlichting & Piglucci, 1998;Pigliucci & Preston, 2004;Draghi & Whitlock, 2012). Because the selection pressures operating on individuals vary in different environments, adaptive phenotypic plasticity may often involve changes in the relative plasticity of several traits rather than the absolute plasticity of a single trait (Lively, 1986;Boersma et al., 1998;Schmitt et al., 1999;Donohue et al., 2000;Weinig, 2000;Van Buskirk, 2002;Sultan, 2003;Badyaev, 2004;Pigliucci, 2004;Relyea, 2004;Draghi & Whitlock, 2012). Variation in the responses generated by individual genotypes in a population, and their adaptive significance, can only be quantified and understood using a methodology capable of capturing and testing for differences in the magnitude and the nature of the response. In this study, a MA-NOVA confirmed the existence of a multivariate G 9 E interaction. However, the two-state multivariate approach was essential for quantifying and interpreting the differences in the multivariate plastic responses of the six different clones. Understanding how the clones generate such differences will require much more work examining how genetic and nongenetic variations (such as maternal effects or epigenetic factors) translate into differences in the developmental processes responsible for generating integrated, functional phenotypes in different environments.
Our results support previous studies of D. magna that demonstrate that a reduction in food availability generally results in delayed maturation and a reduced fecundity (reviewed in Pietrzak et al., 2010;Nogueira et al., 2004). However, we were also able to statistically compare the clonal variation in multivariate plastic responses. Clones Grosse and H01 were the most plastic with respect to the magnitude of the response, which was statistically indistinguishable in the other four clones we studied. But, all the clones were statistically distinguishable with respect to the nature of their plastic responses. Clones DKN13, DKN16, B5 and H01 had phenotypic change trajectories that differed in angle by less than 30 h . These clones matured at similar sizes in high and low food environments, but produced fewer larger offspring in low food environments and higher numbers of smaller offspring in high food environments. In contrast, clones Grosse and B7 demonstrated high plasticity in size at maturity but no plasticity in offspring size and typically smaller clutches. Consequently, their phenotypic change trajectories differed from the other 4 clones by more than 30 h . Understanding and quantifying different plastic responses in this way may be useful for understanding the mechanisms underpinning plastic responses and the selection pressures operating on them (Chun et al., 2007), but also for understanding the consequences of plastic responses. For example, the fact that some clones transmitted food perturbation to their offspring more than others may have implications for understanding the relative importance of genetic and nongenetic inheritance in different populations (Hallsson et al., 2012).
The reduction in fecundity in Grosse could be interpreted as a cost of plasticity, because Grosse also produced the largest plastic response. But this would not explain why clone B7 also showed a reduced fecundity, whereas clone H01 did not. Instead, we suggest that clone Grosse probably has smaller clutch sizes because it matures at much smaller body sizes, whereas clone B7 has smaller clutch sizes because it always invests in very large offspring. These differences were independent of plastic responses to a large extent and were instead explained by clonal differences in neonate size that are already known to be an important driver of phenotypic variation in Daphnia magna (reviewed in Barata & Baird, 1998).
Incorporating phenotypic integration plasticity into studies of phenotypic plasticity also has implications for understanding phenotypic evolution at a more general level. An increasing number of studies have now demonstrated that P-matrices and G-matrices can be altered by the external environment (reviewed in Schlichting & Piglucci, 1998;Murren, 2002;Sgr o & Hoffmann, 2004;Pigliucci, 2006;Kruuk et al., 2008;Townley & Ezard, 2013). As a result, the trait covariances and the constraints that they impose on evolutionary trajectories may change on short timescales potentially altering patterns of adaptation and community dynamics (Tonsor & Scheiner, 2007). In environments that are frequently encountered, selection is expected to optimize multivariate plastic responses resulting in the evolution of adaptive phenotypic integration plasticity as previously demonstrated by plastic responses to different predators (Relyea, 2004), shade avoidance (Schmitt et al., 2003;Bradshaw, 2006), and thigmotropic (Pigliucci, 2004)-, hydrophilic (Bradshaw, 1965)-and competition-mediated (Dudley, 2004) responses in plants. However, the generation of novel patterns of trait integration in response to new or rare environments might also play an important role in allowing populations to adapt to novel environments (Badyaev, 2009;Chevin et al., 2010;Chevin & Lande, 2011;Moczek et al., 2011). The emerging field of eco-devo (Dusheck, 2002;Sultan, 2007;Gilbert & Epel, 2009) proposes that genetic variation in reaction norms (i.e. G 9 E interactions) represents the raw material for microevolutionary responses to environmental change, irrespective of whether the variation is random or adaptive (Sultan, 2007). The results of this study suggest that phenotypic integration plasticity is an important component of any G 9 E interaction that may be missed unless appropriate multivariate analyses are used.

Supporting information
Additional Supporting Information may be found in the online version of this article: Table S1 Univariate ANOVAs comparing the effect of clone, food treatment and their interaction for the traits analysed as part of this study of multivariate plasticity.