Microgeographic differentiation in thermal performance curves between rural and urban populations of an aquatic insect

Abstract The rapidly increasing rate of urbanization has a major impact on the ecology and evolution of species. While increased temperatures are a key aspect of urbanization (“urban heat islands”), we have very limited knowledge whether this generates differentiation in thermal responses between rural and urban populations. In a common garden experiment, we compared the thermal performance curves (TPCs) for growth rate and mortality in larvae of the damselfly Coenagrion puella from three urban and three rural populations. TPCs for growth rate shifted vertically, consistent with the faster–slower theoretical model whereby the cold‐adapted rural larvae grew faster than the warm‐adapted urban larvae across temperatures. In line with costs of rapid growth, rural larvae showed lower survival than urban larvae across temperatures. The relatively lower temperatures hence expected shorter growing seasons in rural populations compared to the populations in the urban heat islands likely impose stronger time constraints to reach a certain developmental stage before winter, thereby selecting for faster growth rates. In addition, higher predation rates at higher temperature may have contributed to the growth rate differences between urban and rural ponds. A faster–slower differentiation in TPCs may be a widespread pattern along the urbanization gradient. The observed microgeographic differentiation in TPCs supports the view that urbanization may drive life‐history evolution. Moreover, because of the urban heat island effect, urban environments have the potential to aid in developing predictions on the impact of climate change on rural populations.

Yet, surprisingly few studies looked at trait differentiation in response to urban heat islands, and we lack evidence whether these differences are genetic (reviewed in Chown & Duffy, 2015;Diamond, Dunn, Frank, Haddad, & Martin, 2015; but see . Considering that the temperature difference between urban and rural areas fall frequently within the range of the expected temperature increases by 2,100 under IPCC (2013) scenarios, adaptation to urban environments can inform on the impact of climate change on organisms (Chown & Duffy, 2015;Youngsteadt, Dale, Terando, Dunn, & Frank, 2015;Youngsteadt, Ernst, Dunn, & Frank, 2017).
Thermal performance curves (TPCs), the continuous reaction norms of organismal performance traits (e.g., growth rate) in response to temperature, are useful tools for the study of differentiation in thermal responses (Sinclair et al., 2016;Stinchcombe & Kirkpatrick, 2012).
The spatial differentiation of TPCs to compensate differences in local temperatures may take three not mutually exclusive adaptive patterns (Angilletta, 2009;Yamahira & Conover, 2002): (i) The TPCs may be shifted horizontally ("hotter-colder" model), with the thermal optimum being higher in warm-adapted genotypes than in cold-adapted genotypes as predicted by theory (Gilchrist, 1995;Lynch & Gabriel, 1987). (ii) The TPCs may be shifted vertically ("faster-slower" model), with cold-adapted genotypes having a higher performance than warm-adapted genotypes at all temperatures (for continuous thermal spatial gradients referred to as countergradient variation) (Conover, Duffy, & Hice, 2009;Conover & Schultz, 1995). Under countergradient variation, genetic and environmental effects oppose each other to produce similar performance in the natural populations along a thermal gradient. (iii) Based on the suggested trade-off between maximal performance and thermal breadth, the TPCs may differ in their width ("generalist-specialist" model), with the performance breadth being narrower in genotypes with higher maximum performance than in genotypes with lower maximum performance. It is theoretically predicted that genotypes adapted to high compared to low temperature variations would have broader TPCs (Gilchrist, 1995;Lynch & Gabriel, 1987), yet studies have provided mixed empirical support (Angilletta, 2009). While "faster-slower" differentiation is not predicted by optimality models of thermal evolution (Gilchrist, 1995;Lynch & Gabriel, 1987), it is often documented for growth rate where it is believed to be driven by time constraints as the increased growth rates in colder environments compensate for the shorter growing season (Conover et al., 2009).
We here focused on geographic differentiation in TPCs for growth rate, a commonly used performance trait (Angilletta, 2009), between replicated urban and rural populations of an aquatic insect using a common garden rearing experiment with a range of temperatures.
While urbanization gradients correspond to temperature gradients (Parris, 2016; for the here studied urbanization gradient: De Ridder, Maiheu, Wouters, & Lipzig, 2015;, we focused on the extreme urban and rural populations of the gradient. Note that this approach is a relevant and powerful setting to test for the type of spatial differentiation in TPCs (horizontal vs vertical shift). To better interpret the TPC pattern for growth rate and possible costs, we also reconstructed TPCs for three other key life-history traits: larval survival, egg development time, and hatchling body size. As study species, we chose the damselfly Coenagrion puella, which is very abundant in both rural and urban areas in Europe (Goertzen & Suhling, 2013).

| Study species and populations
Adult C. puella reproduce in early summer and eggs hatch ca. 3 weeks later. Larval development takes ca. 10 months (Lowe, Harvey, Watts, & Thompson, 2009). This species is univoltine in central Europe (Corbet, Suhling, & Soendgerath, 2006). We studied three rural populations (Bierbeek, Bornem, and Houwaart) and three urban populations (Leuven, Mechelen, and Oudenaarde), all situated within a 45 km radius in Flanders, Belgium ( Figure 1, Appendix S1). All ponds are shallow water bodies with abundant aquatic vegetation. The selection of urban and rural ponds was carried out following a two-step procedure F I G U R E 1 Location of the six study ponds. All ponds were located in Flanders, Belgium (inset) using geographic information system (GIS). First, three urban and three rural 3 × 3 km plots were selected based on the percentage of built-up area: >15% for urban plots and <3% for rural plots. Second, we selected in each plot a pond in a subplot of 200 × 200 m with the same urbanization level. This ensured that both the direct environment (subplot) and the broader surroundings (plot) reflected the same urbanization level. This sampling design was applied in several recent studies where effects of urbanization on organisms are investigated Kaiser, Merckx, & Van Dyck, 2016;Piano et al., 2017;Tüzün et al., 2015Tüzün et al., , 2017. We collected eggs from 10 mated females from each of six ponds in July 2013.

| Experimental setup
We set up a full factorial common garden experiment with five rearing temperatures crossed with two levels of urbanization (urban and rural, each represented by three populations). Animals were reared from the egg stage at one of five constant temperatures: 16, 20, 24, 28, and 30°C. Mean water temperatures experienced by the study species are typically in the range 16-20°C during the early and intermediate periods of development in summer and early fall (Nilsson-Örtman, Stoks, De Block, & Johansson, 2013). Eggs and larvae may experience higher water temperatures during summer (R. Stoks, unpublished data).
Simulations with the lake model Flake (Lake Model Flake 2016) confirmed the occurrence of the here used extreme temperatures in the study region during summer. Although we did not use the more realistic daily fluctuating temperature regime, the ranking of trait values among mean temperature treatments has been shown to be largely unaffected by fluctuating vs. constant temperature regimes (e.g., Fischer, Kölzow, Höltje, & Karl, 2011). We assigned 30 individuals (three larvae per female) per population to each temperature (total of 900 individuals).
Eggs and larvae were kept in dechlorinated tap water and placed in incubators at a photoperiod of 14:10-h light/dark (reflecting the late summer -early fall photoperiod at the study region) at one of the five rearing temperatures. Larvae were reared individually in 200-ml plastic cups and were fed Artemia nauplii 5 days a week (mean ± SE: 212 ± 67 nauplii per feeding portion, n = 12 feeding portions), corresponding to high food levels.
To estimate growth rate, we measured the head width of each larva on days 0 (newly hatched), 30, and 50 using a digital camera attached to a binocular microscope. Head width is an often used measure to estimate size and growth in damselfly larvae (for an example in the study species : Mikolajewski, Brodin, Johansson, & Joop, 2005). The repeated measurements of head width allowed testing for ontogenetic changes in the thermal growth curves .
As the thermal sensitivity of growth rates in the study species strongly depends on the ontogenetic stage Van Doorslaer & Stoks, 2005), we calculated growth rate separately for the periods between days 0-30 and days 30-50 (from here on referred to as first and second period, respectively). This period spans the important growth period of C. puella larvae during late summer and early fall. To calculate growth rate, we used the formula [ln(final head size)ln(initial head size)]/number of days. In animals experiencing seasonal time constraints, such as the here studied damselfly C. puella (Lowe et al., 2009;Mikolajewski, De Block, & Stoks, 2015), growth rate is a relevant performance trait. Moreover, rapid larval growth is important to reach a size advantage in cannibalistic interactions in this study species (Rolff, 1999), especially immediately following hatching when the densities are high.
We calculated larval survival as the ratio of larvae that survived up to day 50. Details on egg size, egg development time, and larval size at hatching are reported and discussed in Appendix S2. The experiment was terminated at day 50, when larvae were old enough to weigh them without causing damage (see Appendix S3 for details).

| Statistical analyses
Unless stated otherwise, all analyses were conducted with R version 3.2.2 for Windows (R Core Team 2015). We used the package "lme4" (Bates, Maechler, Bolker, & Walker, 2015) for mixed-effects models, and the package "car" to compute Wald χ 2 statistic and p-values for fixed effects (Fox & Weisberg, 2011). Significant interactions were further analyzed by comparing least-square means using contrast analysis.
To assess the effects of urbanization level and rearing temperature on growth rate and survival, we used separate (generalized) linear mixed-effects models. We tested for the effects of urbanization level (urban and rural) and temperature (both linear and quadratic term) by including these terms, and their interactions, as fixed effects. As an exception, growth period (first and second period) was included as an additional fixed effect to the growth rate model. We included hatchling size as a covariate to the growth rate model. The following random effects were added where appropriate: population, nested within urbanization level (accounting that animals from the same pond are not independent replicates, thereby avoiding pseudoreplication), female identity of the offspring (accounting for among-brood variation), individual identity of larvae (there were two growth rate estimates, one per period, for each larva; i.e., repeated-measures design). We provide a detailed summary of the model structures in Appendix S4.
To decompose the variation in TPCs into contributions of "hotter-colder", "faster-slower" and "generalist-specialist" models, we used the Template Mode of Variation (TMV; Izem & Kingsolver, 2005) using the code by Izem and Kingsolver (2005) implemented for Matlab (v.8.6.0). This method uses a polynomial function for the decomposition, where each direction of variation is represented by changes in curve-specific parameters (i.e., height, width and optimum temperature). Based on the detected urban-rural differentiation in TPCs using linear mixed-effects models (see Results), we applied this method on growth rates during the second period.
During the first period, increasing temperatures resulted in an increase in growth rate up to 28°C where a plateau was reached (Figure 2a).
During the second period, the TPC had a concave downward shape with decreasing growth rates at temperatures above 24°C (Figure 2b).
Growth rates were higher in rural larvae than in urban larvae, but only in the second period (Urbanization level × Growth period: χ 2 = 6.72, df = 1, p = .01; contrast test, first period: p = .97; second period: p = .02, Figure 2b). The effect of temperature on growth rate was similar for rural and urban larvae, and this was consistent across growth periods (all interactions with temperature p > .4).
The TMV analysis explained 64.65% of the total variation in TPCs.
In line with the results from the linear mixed-effects model, the TMV analysis of the growth rate during the second period revealed that the vertical shift ("faster-slower" model) explained the majority of the variation (58.82%) in TPCs of urban and rural larvae, whereas the horizontal shift of the TPCs ("hotter-colder" model) and the "generalistspecialist" model accounted only for 3.76% and 2.07% of the variation, respectively. Additional outputs of the TMV analysis are reported in Appendix S5.

| DISCUSSION
We found solid support for microgeographic differentiation of the TPCs for growth rate between urban and rural populations. Average In line with this, a recent study reported a 4.03°C difference in mean summer maximum water temperature in a subset of two urban and two rural ponds used in this study . Despite the higher water temperatures in urban ponds, the TPCs for growth rate were not shifted horizontally toward higher optima in urban damselfly populations. Furthermore, although daily temperature fluctuations are more pronounced in urban than in rural ponds in the study region (J. Engelen and K. Brans, unpublished data), there was no differentiation in the width of the TPCs. Instead, TPCs were shifted vertically with rural larvae having higher growth rates across rearing temperatures during the second period. The observation that the vertical shift in growth TPCs was consistent across the three sampled populations per urbanization level suggests this pattern to be general (see Appendix S6 and Appendix S7 for more detail). While the here chosen design of replicated populations studied at the extreme ends of the urbanization gradient is powerful to detect microgeographic shifts in TPCs, our setup does not allow to infer that these shifts occur gradually across the urbanization gradient. As we will argue that time constraints as- There was a pronounced ontogenetic shift in the TPC for growth rate. While during the first period, growth rates were faster and the TPC had a steep slope and a high thermal optimum, there was a shift in the second period toward lower growth rates and a downward concave-shaped TPC with a lower thermal optimum. Such ontogenetic shifts in TPC for growth have been associated with the decreasing temperatures larvae are exposed to during early and intermediate growth periods, and have been reported before in Coenagrion damselfly species, including the study species C. puella Van Doorslaer & Stoks, 2005). The high thermal optimum in the first period (28°C) is not often reached in the study ponds. This fits the pattern that thermal optima are typically higher than modal environmental temperatures (Angilletta, Huey, & Frazier, 2010). This has been explained by the typically asymmetric shape of TPCs resulting in drastically lower performance above the thermal optimum than below the thermal optimum (Deutsch et al., 2008). Even if the optimal temperatures are met infrequently, it may still be beneficial to optimize growth at higher temperatures as even short exposures to higher temperatures may allow significant gains (Kingsolver, 2000). Selection for high growth rates and reaching a larger size may be stronger in early larval stages because of the initial high larval densities (leading to increased cannibalism; Hopper, Crowley, & Kielman, 1996) and higher temperatures (leading to increased competition and intraguild predation; Brown, Gillooly, Allen, Savage, & West, 2004). This, together with the identified trade-off between early and late growth rates in Coenagrion damselfly larvae , may have caused the higher growth rates in the first compared to the second period. Alternatively, the lower growth rates in the second period could be an artifact driven by size-dependent growth rate (Tammaru & Esperk, 2007), where the lower growth rates of the larvae during the second period were due to their on average larger size. This scenario, however, is unlikely to drive the observed pattern as we controlled for the larval size (by including either initial size or time-varying sizes as covariate; see Appendix S3) when analyzing growth rate. Note that given size is fixed per developmental stage, this also corrects for any differences in developmental stage between temperature treatments.
Although we followed the growth rate TPCs only until day 50, this captures the major part of the prewinter growth period. Moreover, for the study species, the shape of the TPCs between days 30 and 50 reflects the TPCs until the end of the prewinter growing season (based on Van Doorslaer & Stoks, 2005;Nilsson-Örtman, Rogell, Stoks, & Johansson, 2015).
A key finding of our study was the vertical shift in TPC for growth rate during the second period with higher growth rates in rural larvae compared to urban larvae. This vertical shift in TPCs is consistent with a pattern of countergradient variation (Conover & Schultz, 1995;Conover et al., 2009). We did not find any evidence for the theoretically predicted higher thermal optima ("hotter-colder" model, Lynch & Gabriel, 1987;Gilchrist, 1995) of urban damselfly larvae that would suggest thermal adaptation. The majority of evidence for thermal adaptation comes from interspecific studies (e.g., Angilletta et al., 2010;Frazier, Huey, & Berrigan, 2006), whereas intraspecific studies usually favor the "faster-slower" model (e.g., Izem & Kingsolver, 2005;Richter-Boix et al., 2015;Yamahira & Conover, 2002). It has been suggested that thermal adaptation ("hotter-colder" model) requires radical changes in the genetic structure (e.g., mutations that would allow for an increased growth via horizontal shifts in thermokinetics of enzyme function), which requires time of the magnitude similar to what is needed for the divergence of species (Yamahira & Conover, 2002;Yamahira, Kawajiri, Takeshi, & Irie, 2007). Furthermore, in spite of the higher daily temperature fluctuations in urban ponds (J. Engelen and K. Brans, unpublished data), urban populations did not have wider TPCs ("generalist-specialist" model; Lynch & Gabriel, 1987;Gilchrist, 1995). This matches the general pattern that thermal specialization is only encountered infrequently (Angilletta, 2009; but see Latimer, Wilson, & Chenoweth, 2011;Richter-Boix et al., 2015).
The majority of studies reporting vertical shifts in TPCs for growth rates show this pattern across large geographic gradients with little gene flow among populations (Conover et al., 2009), whereas only a handful of studies demonstrated this pattern at a microgeographic scale where gene flow can be high (but see Blanckenhorn, 1991;Skelly, 2004;Richter-Boix, Teplitsky, Rogell, & Laurila, 2010;Richter-Boix et al., 2015). Also in the here studied C. puella, gene flow is high at the studied scale in Flanders (as indicated by the very low F STvalues of ca. 0.02 for this study species, see Figure 2 in Johansson et al., 2013). This indicates that selection for higher growth rates is especially strong to counteract the homogenizing force of gene flow, thereby supporting accumulating evidence that local adaptation may occur at small spatial scales in the presence of gene flow (reviewed in Richardson et al., 2014).
Given that we obtained the TPC differentiation in a common garden experiment suggests that the pattern reflects genetic adaptation, rather than environmental differences. Although we used first generation animals, maternal effects seem unlikely to have played a major role because rural and urban populations did not differ in egg size and hatchling size (Appendix S2), while differences in egg size are the most common way how maternal effects are transferred (Mousseau & Dingle, 1991). More importantly, larval growth rates did not differ between rural and urban populations during the first period. It seems unlikely that any nongenetic effects transferred by the mother would have caused differences in growth rate between population types that would not be present during the first growth period. Indeed, if anything, maternal effects tend to decay throughout ontogeny (e.g., Lindholm, Hunt, & Brooks, 2006). This fits the pattern based on quantitative genetic rearing experiments that maternal effects on growth rate are absent or small in damselfly larvae (Shama et al., 2011;Sniegula, Golab, Drobniak, & Johansson, 2016). Nevertheless, we cannot fully exclude mechanisms acting via maternal effects such as transgenerational plasticity (e.g., Richter-Boix, Orizaola, & Laurila, 2014).
The observation that the vertical shift in TPCs only was apparent in the second period may be explained by the stronger time constraints, the major selective force underlying a vertical shift in TPCs for growth rates (Conover et al., 2009). Indeed, during the second period, temperatures encountered in the field are lower while larvae face stronger time constraints to reach a certain developmental stage before the onset of winter. In support of a stronger selection to accelerate growth rates in rural populations, growing seasons are reported to be considerably shorter in rural than in urban areas (Yang et al., 2013;Somers et al., 2013; for the study region: J. Engelen and K. Brans, unpublished data). Higher growth rates under time constraints have been theoretically predicted (Abrams, Leimar, Nylin, & Wiklund, 1996) and empirically shown in Coenagrion damselfly larvae, including the study species (Mikolajewski et al., 2015). However, systematic differences between urban and rural ponds other than temperature might have contributed to the differentiation in TPCs. Notably, differences in predation pressure have been identified as a key factor driving the shift in TPCs for growth rate (Richter-Boix, Quintela, Kierczak, Franc, & Laurila, 2013). Predators of damselfly larvae did, however, not differ in densities between the here studied urban and rural ponds (see Appendix S8). In agreement with our finding, aquatic macro-invertebrate community structures are not driven by urbanization in the study region (C. Souffreau, personal communication), and a recent study reported no difference in species richness of aquatic macro-invertebrates from urban and nonurban ponds in the United Kingdom (Hill et al., 2017).
Nevertheless, even when predator densities did not differ between urban and rural ponds, predators may impose stronger selection for reduced growth rates in the urban ponds because of the typical higher predation rates at warmer temperatures (e.g., De Block et al., 2013;Sentis, Morisson, & Boukal, 2015). Larvae of the study species indeed reduce growth rates in the presence of predators (Mikolajewski et al., 2005). The higher temperatures in the urban populations may therefore be driving the lower growth rates in urban populations both by causing lower time constraints and by generating higher predation rates. Furthermore, differences in pollution levels between pond types might have contributed to the observed shift in TPCs. Indeed, tadpoles exposed to a herbicide had on average upwards shifted ("fasterslower" differentiation) and narrower TPCs ("generalist-specialist" differentiation) for swimming speed (Katzenberger et al., 2014). Yet, the expected higher contamination levels in urban waterbodies (e.g., Gilliom, 2007), together with the absence of a generalist-specialist trade-off in our study, make this scenario unlikely.
Despite the widespread evidence of urban heat islands (Gaston et al., 2010), surprisingly few studies documented differentiation in TPCs between rural and urban populations, and as these all directly measured traits on field-collected animals, it is unknown to what extent these differentiations are genetic rather than environmentally driven (reviewed in Chown & Duffy, 2015). In the only study considering continuous TPCs, which is essential to discriminate between the models (Richter-Boix et al., 2015), two different patterns of thermal differentiation were observed when rearing four species of soil fungi directly isolated from one urban population and from one rural population (McLean, Angilletta, & Williams, 2005). Two fungus species followed the "hotter-colder" model, with isolates from the urban population growing faster at 26°C, but slower at 18°C compared to isolates from the rural population. In the two other fungus species, urban isolates grew as fast or faster at all temperatures than the rural isolates ("faster-slower" model). This is the opposite vertical shift in TPCs of that we here documented. The authors suggested this is because high temperatures in urban areas inhibit fungal growth, resulting in restricted growing seasons for urban populations (McLean et al., 2005), while the higher temperatures in the urban damselfly populations, instead, were beneficial and increased the length of the growing season.
In contrast, two studies report higher heat tolerance in urban compared to rural populations of field-collected leafcutter ants (Angilletta et al., 2007) and common garden-reared acorn ants (Diamond, Chick, Perez, Strickler, & Martin, 2017), yet none of them tested differential thermal responses in life history. Finally,  have shown genetic adaptation to urbanization in terms of a higher heat tolerance in urban populations of the water flea Daphnia magna.
Studies on vertical shifts in TPCs typically assume life-history trade-offs where the faster growing genotypes pay costs that preclude them from occupying the entire thermal gradient (Conover & Schultz, 1995;Conover et al., 2009). Costs of rapid growth may be manifold (Dmitriew, 2011), and in damselfly larvae include decreased investment in energy storage (Stoks, De Block, & McPeek, 2006), reduced immune function (De Block & Stoks, 2008b), and increased oxidative stress (De Block & Stoks, 2008a). These costs may explain the overall lower survival in rural compared to urban populations. Similarly, Hong and Shurin (2015) showed that northern populations of the tidepool copepod Tigriopus californicus grew faster than southern populations, yet suffered a reduced survival. However, we cannot exclude the possibility that adaptations to urbanization-related stressors (e.g., pollution) might have resulted in an overall higher resilience of urban populations, reflected in their higher survival.
Urbanization is a major driver of microevolutionary change (Alberti, 2015;Alberti et al., 2017). Urban environments have the potential to provide unique information and aid in developing predictions on the impact of climate change on organisms (Chown & Duffy, 2015;Youngsteadt et al., 2015Youngsteadt et al., , 2017, and the use of cities as natural experiments has recently been promoted. Our finding that individuals from urban and rural populations consistently differed in growth rate and survival is relevant for predicting climate change impact, as the summer temperature difference between urban and rural ponds falls within the range of the expected temperature increases by 2100 under several IPCC (2013) scenarios. The urban populations can therefore be used as proxies to understand and predict the impact of global warming in rural populations under gradual evolution (Stoks, Geerts, & De Meester, 2014). Aside from contributing to the limited evidence for vertical shifts in TPCs at a microgeographic scale (Richter-Boix et al., 2015;Skelly, 2004), this study is the first using replicated populations in a common garden experiment from the egg stage to report evidence for urbanization-associated countergradient variation in a key performance trait.