Bergmann's rule is maintained during a rapid range expansion in a damselfly

Climate‐induced range shifts result in the movement of a sample of genotypes from source populations to new regions. The phenotypic consequences of those shifts depend upon the sample characteristics of the dispersive genotypes, which may act to either constrain or promote phenotypic divergence, and the degree to which plasticity influences the genotype–environment interaction. We sampled populations of the damselfly Erythromma viridulum from northern Europe to quantify the phenotypic (latitude–body size relationship based on seven morphological traits) and genetic (variation at microsatellite loci) patterns that occur during a range expansion itself. We find a weak spatial genetic structure that is indicative of high gene flow during a rapid range expansion. Despite the potentially homogenizing effect of high gene flow, however, there is extensive phenotypic variation among samples along the invasion route that manifests as a strong, positive correlation between latitude and body size consistent with Bergmann's rule. This positive correlation cannot be explained by variation in the length of larval development (voltinism). While the adaptive significance of latitudinal variation in body size remains obscure, geographical patterns in body size in odonates are apparently underpinned by phenotypic plasticity and this permits a response to one or more environmental correlates of latitude during a range expansion.


Introduction
A major challenge in the study of global ecological change is to predict ultimate consequences from proximate biotic responses to environmental change. One predictive framework arises from static patterns of clinal variation in body size (and associated life-history traits), a common phenomenon whereby individuals become either larger (Bergmann, 1847;Atkinson & Sibly, 1997;Blackburn et al., 1999) or smaller (Park, 1949;Mousseau, 1997) with increasing latitude, presumably a response to variation in the length of the growth season (Chown & Gaston, 1999) and/or temperature (Atkinson, 1994). For example, most ectotherms exhibit an increase in body size at lower developmental temperatures (the temperature size rule; Atkinson, 1994), albeit with the underlying mechanisms poorly understood (Atkinson, 1994;Forster et al., 2012). Moreover, under counter-gradient variation (Conover & Schultz, 1995), individuals at high latitudes attempt to compensate for the shorter season by growing faster than low latitude conspecifics; perfect compensation for seasonal time constraints will result in a uniform phenotype, but under-or overcompensation of growth rate will be manifest as apparent converse Bergmann or Bergmann clines, respectively. While there is some indication that larger terrestrial arthropods almost uniformly exhibit converse Bergmann clines (Blanckenhorn & Demont, 2004), several reviews of latitude-size relationships in insects have yielded few consistent patterns (Chown & Gaston, 2010;Shelomi, 2012).
A second component of the predictive framework involves the dynamics of species' ranges. Morphological changes in colonists during a range expansion have been associated with selection for increased dispersal ability (Parmesan, 2006;Hill et al., 2011). As many high latitude species are apparently attempting to maintain their fundamental niches, the direction of range expansion correlates with latitude (Hill et al., 2011). The extent to which dispersive phenotypes complement or are antagonistic with the expected latitudinal variation in body size during the process of a range expansion is unclear. For example, larger arthropod species with development times that are relatively long compared with the growth season are susceptible to seasonal time constraints (Chown & Gaston, 1999) and tend to decrease in size with latitude (Blanckenhorn & Demont, 2004). However, the relationship between body size and dispersal ability within arthropod taxa is inconsistent, being either positive (Lepidoptera; Stevens et al., 2012), negative (Odonata: Anisoptera;McCauley, 2010) or variable (Odonata: Zygoptera; Conrad et al., 2002). Some additional complexity to the pattern of phenotypic expression that occurs during a range expansion may be associated with the numbers of individuals that comprise the expanding front. High levels of gene flow can limit the capacity for local adaptation (Lenormand, 2002), while colonization by relatively few individuals allows founder effects and drift (stochastic processes) to predominate. As the phenotypic response to a heterogeneous environment can determine a species' capacity to shift its range (Gunnarsson et al., 2012), it is important to understand the extent to which trait expression differs from the background level of population differentiation due to neutral, demographic processes (McKay & Latta, 2002;Brommer, 2011).
Odonata (dragonflies and damselflies) are excellent model organisms to study the processes that occur during range expansions as many species are altering their ranges in response to environmental warming (Hickling et al., 2005;, 2010. Moreover, body size is a key life-history parameter for odonates, for example affecting clutch size (Cordero, 1991), fitness (Michiels & Dhondt, 1991;but cf. Thompson & Fincke, 2002) and trophic ecology (Thompson, 1978). Studies that quantify latitudinal gradients in odonate morphology are limited. Hassall (2013) demonstrated that temperature and time stress generate U-shaped relationships between size and latitude (i.e. converse Bergmann clines and Bergmann clines at lower and higher latitudes, respectively) in Calopteryx maculata, supporting Johansson's (2003) research on Enallagma cyathigerum but raising doubts about whether the higher latitude response was driven solely by changes in voltinism. At a more regional scale, Hassall et al. (2009) found the strength of latitudinal variation in body size to differ between zygopteran species in the UK (ca. 51-56.6°N). Northern individuals of Coenagrion puella were heavier than their southern counterparts (with no significant correlation between wing length and latitude), and while a similar trend was apparent for Pyrrhosoma nymphula, neither body mass nor wing length significantly correlated with latitude. The contrast may be explained by relative position within a species' range, with C. puella collected from its northern range margin where there are both semivoltine and univoltine populations, while samples of P. nymphula were at the middle of its latitudinal range and are semivoltine.  noted a difference in range expansion, with C. puella a more recent colonizer than P. nymphula. Odonates thus exhibit complex patterns of geographical variation in body size and can exhibit a weak 'Bergmann cline' at higher latitudes. Previous work has implicated a genetic component to phenotypic variation in E. cyathigerum (De Block et al., 2008); however, whether clinal variation in body size is affected by conditions of high gene flow, for example during a rapid range expansion, remains to be determined.
The zygopteran Erythromma viridulum (Charpentier) is a coenagrionid damselfly that inhabits ponds, lakes and ditches. E. viridulum is bivoltine at the southern edge of its range in Greece (Galletti & Pavesi, 1983) and semivoltine in the UK (Keat, 2007), but no data exist for other European populations. E. viridulum has made an extensive range expansion into north-western Europe in the last 30 years (Ketelaar, 2002) and is the first recorded example of a migrant damselfly establishing colonies in the British Isles (Dewick & Gerussi, 2000). Until the 1980s, this species was absent from northern France and was recorded infrequently in the Netherlands (Askew, 1988), but by the mid-1990s, E. viridulum had colonized northern continental Europe. Large numbers of E. viridulum were observed at a few sites in south-east UK in 1999 (Dewick & Gerussi, 2000) and in southern UK in 2000 (Cham, 2001). This species subsequently expanded its British range at roughly 30 km yr À1 (Fig. 1). Concomitant with these observations, genetic studies have revealed high levels of gene flow among UK populations (Watts et al., 2010). The early records and sampling of this invasion permitted a uniquely detailed investigation of how high levels of gene flow impacts phenotypic variation at the forefront of an expanding range margin.
We quantify morphological variation along the north-west edge of E. viridulum's range to identify whether a latitude-body size relationship is maintained during a rapid range expansion. It might be expected that the recently established populations at the leading edge of the range expansion would exhibit similar morphology due to constraints imposed by high levels of gene flow. We examine the ability of latitude, temperature, location and the level of genetic divergence among populations to explain the pattern of morphological variation and find that variation in body size is primarily influenced by latitude and is not constrained during a rapid range expansion.

Sample collection and processing
Between 2002 and 2006, 777 adult male E. viridulum were captured from a total of 21 sites: 12 in southern UK and 9 (from 4 countries: Germany, Belgium, the Netherlands and France) in northern continental Europe ( Fig. 1; Table 1). While it is impossible to discriminate immigrants from residents at any site, we assume that most, if not all, of the individuals that were sampled had emerged from the water body at which they were caught. Despite the recent immigration of E. viridulum from the continent, coenagrionid damselflies rarely disperse more than 1-2 km when they reside within a suitable habitat matrix (Conrad et al., 1999;Watts et al., 2004). If immigrants were present in the samples, they likely constituted a small proportion of the individuals, and thus, the patterns here are assumed to represent a developmental response to the local environment. Seven phenotypic traits were measured on all specimens to the nearest 0.1 mm using either digital dial callipers or a light microscope and eyepiece micrometre: head width, wing length, wing width, thorax length, thorax width, tibia length and abdomen length. A subsample of 713 specimens from the 12 largest samples were genotyped at 10 microsatellite loci (Keat et al., 2005) with genotyping procedures provided in Watts et al. (2010), as incorporating data from the other sites with relatively few genotypes (mean sample size n = 9) would introduce noise. Temperature data were extracted from the 30 arc seconds (ca. 1 km) resolution WorldClim data set (Hijmans et al., 2005) to provide estimates of mean annual temperature (Temp ANN ) and mean winter temperature (Temp WIN , the mean temperature of December, January and February) for every site.

Data analyses
As all phenotypic traits were significantly correlated (r > 0.5, P < 0.001 in all cases), a principal components analysis (PCA) was used to summarize the overall variation in body size. The first principal component (PC1) explained 92% of the variation and was negatively correlated with all traits: PC1 was 'absolute' transformed (-PC1) to give a positive measure of overall body size, hereafter 'SIZE'. We tested for a correlation between latitude and mean population SIZE in (i) the full sample (n = 21); (ii) the UK sites only (n = 12); and (iii) the European sites only (n = 9). We then identified the subset of variables that best predicted variation in body size using generalized linear mixed effects models (GLMMs), with the individual as the unit of replication, constructed in the nlme package (Pinheiro et al., 2013) in R (R Development Core Team, 2013); SIZE was the response variable examined against four predictor variables, latitude, Temp ANN , Temp WIN , and location and with site always included as a random effect. While latitude as a predictor has no biological relevance per se, odonate body size is influenced by photoperiodic time stress which varies systematically with latitude (Hassall, 2013). The eight potential models (including the null model with a floating intercept) were compared using AIC in the AICcmodavg (Mazerolle, 2012) package in R (R Development Core Team, 2013) to select the most parsimonious model. We also examined for a latitude-location interaction (i.e. SIZE~latitude + location + latitude 9 location) within the GLMM framework to test for a difference in the latitude-size pattern between the UK and continental Europe. We evaluated the potential for a confounding effect of multicollinearity among predictors using variance inflation factors (VIFs) calculated for the models' fixed effects.
To examine the effect of genetic divergence among populations upon phenotypic variation, we transformed the  Table 1. pattern of phenotypic variation to the measure 'P ST ', a measure of phenotypic differentiation among pairs of populations that is analogous to F ST , a measure of genetic distance among samples using methods described by McKay & Latta (2002). In its original derivation, phenotypic divergence 'Q ST ' requires an estimate of additive genetic variance that should be obtained from common garden experiments (McKay & Latta, 2002); because such extensive larval breeding was not possible, the measure of phenotypic divergence P ST includes nonadditive effects (Brommer, 2011). Components of trait variance are partitioned within and among populations, with the proportion of the total variance among populations acting as a measure of the differentiation in the trait(s) between populations. Where P ST is greater than F ST (for neutral loci), the trait has diversified more than would be expected purely on the basis of gene flow and genetic drift implying a selective response to the environment.
Genetic differentiation between all pairs of samples (F ST ) was calculated using FSTAT v.2.9.3 (Goudet, 2001). To quantify the effect of spatial structuring upon genetic and phenotypic variation, the matrices of pairwise F ST and P ST were correlated against a matrix of the linear distances separating sample pairs (ln km; Rousset, 1997) using Mantel tests. A third Mantel test was used to test for a correlation between the matrix of F ST and the corresponding matrix of P ST to determine whether there is any relationship between the level of genetic and phenotypic divergence among samples. The significance of each Mantel test was calculated by permutation procedure (n = 999 permutations).

Results
E. viridulum exhibited substantial variation in body size among sites, with individual wing length varying from 15.9 to 20.6 mm and abdomen length varying from 20.9 to 26.9 mm. Body size (SIZE) was positively associated with latitude in the data set as a whole (Pearson's correlation, r = 0.717, P < 0.001, n = 21) continental European populations (r = 0.664, P = 0.051, n = 9) and in the UK (r = 0.671, P = 0.017, n = 12) (Fig. 2). The best-fit model describing variation in body size involved latitude, location and the latitude 9 location interaction, although the explanatory power was similar to the same model without that interaction (ΔAICc = 1.63, Table 2). After centring the latitude predictor, VIFs for the predictors in the top model of SIZE~latitude + location + latitude 9 location were 3.117, 1.139 and 2.919, respectivelywell below the guideline value of 10 where multicollinearity is thought to interfere with models (Kutner et al., 2005) and providing confidence in the parameterization of our models (which are presented with uncentred predictors for ease of interpretation). Body size was marginally greater in the UK than in the older populations in continental Europe (hence the selection of the location term in the two best models; Tables 2 and 3), although there is no conspicuous difference to the latitude-size relationship between the populations from the UK and continental Europe that would be indicative of a fundamental location effect (Fig. 2). Indeed, while the location and latitude 9 location interaction terms contributed to improved fit in the two top models (Table 2), neither were statistically significant ( Table 3) which suggests that any location effects are minor. There was no evidence for a sawtooth cline in the relationship between size and latitude that might indicate a change in voltinism across the study area (Fig. 2). The level of genetic differentiation among all pairs of populations was low, with average F ST = 0.032 AE 0.003 (SE) and all pairwise values of F ST < 0.09. Distance separating the sample sites had no significant effect on the amount of genetic differentiation between sites (r = 0.148, P = 0.202, Fig. 3a), which is consistent with the recent range expansion. Mean P ST = 0.264 AE 0.033 (SE) was an order of magnitude greater than F ST , and accordingly, pairwise estimates of phenotypic variation (P ST ) among locations generally were greater than the corresponding estimates of F ST (a pattern consistent with other studies, Brommer, 2011). We did not find a significant relationship between P ST and the distance separating populations (r = À0.186, P = 0.904, Fig. 3b) or between genetic (F ST ) and phenotypic (P ST ) differentiation among sample locations (r = 0.036, P = 0.372, Fig. 3c).

Discussion
Body size is an important biological trait, influencing behaviour, physiology and fitness. Like many animals, E. viridulum body size varies with latitude, with a general increase in size at high latitudes, an apparent Bergmann cline. We demonstrate that this latitudinal variation in body size occurs during this species' extensive northwards range expansion and under conditions of high gene flow that might have been expected to constrain morphology. In fact, emerging adults of the newly founded populations in the UK almost immediately (i.e. within a few generations) fit into the 'appropriate place' on a latitudinal cline in body size. The lack of evidence for a shift in voltinism in our sample populations implies that this pattern may occur within a group of populations exhibiting the same mode of voltinism.  While few consistent latitude-body size patterns have been identified in insect taxa (Blanckenhorn & Demont, 2004;Chown & Gaston, 2010;Kivel€ a et al., 2011;Shelomi, 2012), predicting the latitude-body size response of odonates will be affected by the semi-aquatic life history. For example, Forster et al. (2012) showed a trend towards a negative relationship between temperature and body size in large aquatic species, but a more positive (and more variable) relationship in large terrestrial species. Structural growth occurs only during the aquatic (larval) period in Odonata. This may explain why the clear increase in body size of adult E. viridulum, and other odonates (Johansson, 2003;Hassall, 2013), at higher latitudes contrasts with the typical pattern exhibited by large terrestrial insect taxa (Blanckenhorn & Demont, 2004) but is consistent with the broader pattern described by Forster et al. (2012). In the case of E. viridulum, the apparent speed of response appears to suggest that phenotypic plasticity during larval growth drives the latitude-size relationship. However, does the broader pattern in Odonata and other aquatic insects imply greater selection on the aquatic phase? Even so, the body-size response of odonates is somewhat counter-intuitive as small size at high latitudes is an expected response to a shorter larval growth season (Mousseau, 1997;Blanckenhorn & Demont, 2004). As the impact of seasonal time constraints upon growth will be more severe in more northerly populations, one mechanism by which individuals may complete development and emerge as viable adults is to extend the duration of the larval phase (Johansson, 2003;Cabanita & Atkinson, 2006). For organisms with a flexible larval developmental period, an interaction between shifts in voltinism and growth season is predicted to manifest as a sawtooth cline in body size whereby size decreases with latitude until the growth season is so short that there is an obligate extension to the larval period and concomitant increase in size (Roff, 1980). A sawtooth cline was not apparent from Johansson's (2003) study of E. cyathigerum, most likely because the large distances (ca. 2000 km) separating high latitude samples (from Belgium, Sweden and Norway) coincided with shifts in voltinism. The finer spatial resolution of our sampling (populations separated by up to ca. 400 km) addresses this issue, but still provides no evidence for a sawtooth cline in body size. It is hard to argue that larval developmental period in our sample of continental European populations was different to the UK populations, which are semivoltine (Keat, 2007). At least for E. viridulum, a two-year developmental period may relax the effects of seasonal constraints compared with an attempt to complete development within a single growth season. With this in mind, it is notable that univoltine adult E. cyathigerum decreased in size with an increase in latitude between 42°N and 57°N (Johansson, 2003). The larval stage is an important component in determining the rate of range expansion (e.g. Crozier, 2004), and variation in latitudinal body-size relationships among arthropod taxa may reflect the ability of the immature phase to tolerate (cold) winters and thus develop over successive growth seasons.
Latitude encompasses a range of abiotic components that affect growth, particularly temperature and photoperiod, and odonate larvae certainly are capable of marked physiological responses to these factors (De Block & Stoks, 2003;De Block et al., 2008). Like Johansson (2003), we observe a strong effect of latitude, but in contrast to that study, we found no effect of temperature. Indeed, there is a contrast between laboratory experiments that have shown body size in ectotherms to increase at lower developmental temperatures (larger size at lower temperatures, Atkinson, 1994) and the complex latitude-body size relationships exhibited by empirical systems (Chown & Gaston, 2010), presumably due to additional drivers of body size. Temperature determines the length of the growing season for many insects through thermal thresholds for development, although where range expansions maintain a thermal niche the effect on body size may not be important. Where photoperiod determines the length of the growing season, movement of individuals into a new latitudinal zone with shifting climate space may result in shorter growing seasons at a given temperature regime. Decreasing photoperiods reduce odonate size in laboratory experiments when temperature is kept constant (Johansson & Rowe, 1999), and with photoperiod and temperature, working synergistically to reduce body size when the latter is increased (De Block & Stoks, 2003). Such an effect of photoperiod (reduced growth season) is not observed here, implying for example some adaptive changes to any new photoperiod regimes (Harada et al., 2005), a response to some other correlate(s) of latitude, or over-compensation in growth rate. The latitude-size relationship is essentially similar for the older populations from continental Europe and in the recent UK colonists (Fig. 2) and presumably is adaptive. A key issue is whether any phenotypic response relies upon standing genetic variation or is plastic. In contrast to the effect of latitude per se (Table 1), the pattern of phenotypic variation exhibited by E. viridulum is decoupled from the distance separating populations (Fig. 3b) and the amount of genetic differentiation between areas (Fig. 3c). The implication is that an environmental correlate of latitude elicits a strong selective pressure on body size and that E. viridulum is capable of mounting an appropriate phenotypic response despite the potential homogenizing effect of high gene flow. Of course, gene flow may increase the efficacy of selection under some circumstances (Holt & Gomulkiewicz, 1997). However, interpreting the action of selection using estimates of P ST is complicated when the additive genetic component to the trait is unknown (Brommer, 2011). While many species are assumed to have a substantial additive genetic component to their body size (Brommer, 2011), this does not appear to be the case for odonates; Lowe et al. (2009) failed to find significant heritability for body size in Coenagrion puella. More relevant would be the role of phenotypic plasticity in facilitating colonization of novel habitats (Knop & Reusser, 2012). For many odonates, larval development rate appears to be quite plastic, with individuals capable of altering their developmental period in response to local conditions De Block et al., 2008;Watts & Thompson, 2012); delayed larval development has been observed in other aquatic insect taxa (Cabanita & Atkinson, 2006). Conversely, a genetic component to latitudinal variation in body mass in E. cyathigerum, consistent with the size variation described by Johansson (2003), was manifest through variation in growth rates among areas and an interaction with development time (De Block et al., 2008).
Under conditions of high gene flow during a rapid range expansion, odonates are capable of sufficient developmental phenotypic plasticity that ultimately determines adult body size to allow an appropriate response to a new environment ( Fig. 2; see also Watts et al., 2010). Such plasticity could also facilitate subsequent improved adaptation through selection at colonized areas. Studies of fitness components in odonates generally have emphasized the importance of stabilizing selection (Fincke, 1982;Stoks, 2000) for male body size, and yet the apparent optimal body size varies considerably with latitude. It would be interesting to quantify the fitness consequences of variation in body size throughout a species' range. Such large-scale patterns in the extent of phenotypic variation, the underlying drivers in the field, and their adaptive and ecological consequences remain to be uncovered.