https://orcid.org/0000-0002-9840-1050
Abstract
Understanding the interplay between genetic differentiation, ancestral plasticity, and the evolution of plasticity during adaptation to environmental variation is critical to predict populations’ responses to environmental change. However, the role of plasticity in rapid adaptation in nature remains poorly understood. We here use the invasion of the horned beetle Onthophagus taurus in the United States during the last half century to study the contribution of ancestral plasticity and post-invasion evolution of plastic responses in rapid population differentiation. We document latitudinal variation in life history and morphology, including genetic compensation in development time and body size, likely adaptive responses to seasonal constraints in the North. However, clinal variation in development time and size was strongly dependent on rearing temperature, suggesting that population differentiation in plasticity played a critical role in successful adaptation on ecological timescales. Clinal variation in wing shape was independent of ancestral plasticity, but correlated with derived plasticity, consistent with evolutionary interdependence. In contrast, clinal variation in tibia shape aligned poorly with thermal plasticity. Overall, this study suggests that post-invasion evolution of plasticity contributed to range expansions and concurrent adaptation to novel climatic conditions.
Natural selection on standing genetic variation and novel mutations affecting mean phenotype expression represents a major avenue toward adaptation to environmental change (Dobzhansky 1982; Barrett and Schluter 2008). However, organisms may also exhibit plastic responses to environmental variation, thereby allowing phenotypic adjustments to environmental conditions in the absence of changes in a population's genetic composition (Schlichting and Pigliucci 1998; West-Eberhard 2003; Whitman and Ananthakrishnan 2009). Such plastic, within-generation responses to environmental conditions are abundant and may impact subsequent evolutionary changes in a variety of ways (West-Eberhard 2003; Crispo 2008; Hendry 2015). For example, while plasticity can be maladaptive or neutral, it may also facilitate later adaptation by maintaining a population's fitness until novel, beneficial mutations emerge (e.g., Corl et al. 2018), or alternatively, may hamper selection from removing maladaptive alleles by buffering their phenotypic effects and shielding them from selection (e.g., Huey et al. 2003). In addition, plasticity itself can evolve, impacting the phenotypic variation visible to selection with the potential to influence direction and magnitude of adaptive responses (West-Eberhard 2003; Moczek et al. 2011; Hendry 2015).
The role of plasticity and the evolution of population differentiation in plasticity during local adaptation are particularly relevant to our understanding of how organisms react to global climate change (Merilä 2012; Kingsolver and Buckley 2017; Garnas 2018; Oostra et al. 2018). Thermal plasticity is especially widespread (Angilletta 2009; Sinclair et al. 2012), facilitates immediate, short-term responses to variation in weather, and has been proposed to not only precede, but also affect subsequent genetic changes in populations beyond a given individual's lifetime as a response to changing climates (Kelly 2019). Yet, our understanding of how ancestral plasticity influences subsequent evolution in natural populations remain scarce, especially on “ecological timescales” (Fox et al. 2019). This is because, in order to study how ancestral plasticity relates to genetic change, plasticity must first be quantified in the ancestral population, and then be contrasted to genetic differentiation after evolution has taken place. As ancestral plasticity and evolutionary change are not simultaneously accessible, it has been exceedingly difficult to investigate if and how ancestral plasticity may bias evolutionary change (but see, e.g., Waddington 1952; Suzuki and Nijhout 2006; Shaw et al. 2007; Wund et al. 2008; Levis et al. 2018; Casasa et al. 2020).
Studying species that successfully and rapidly invaded new habitats that are climatically different from the ancestral range offers an exceptional opportunity to study how organisms cope with and adapt to rapidly changing environments (Johnston and Selander 1964; Gilchrist et al. 2001; Kingsolver et al. 2007; Geng et al. 2016), and, provided that the ancestral source population (or a proxy thereof) is known and still accessible, to investigate the role of ancestral plasticity therein (Moczek 2007; Casasa and Moczek 2018). A species that presents itself as particularly useful in this context is the dung beetle Onthophagus taurus. Native to the Mediterranean region and parts of Central Europe (see hatched area in Fig. 1), this species was accidentally introduced to the south-eastern United States in the 1970s, successfully established exotic populations, and progressively expanded its range northward (Hoebeke and Beucke 1997) and was recorded close to the Canadian border by 2011 (Rounds and Floate 2012). That is, O. taurus overcame a distance of ∼1700 km within just four decades (∼80 to 100 generations), and now inhabits a climatic niche that is colder and more humid compared to its ancestral distribution range (DaSilva et al. 2016). Leveraging the colonization and invasion history of O. taurus, we here investigate genetic differentiation and thermal plasticity along an environmental gradient in the exotic range, and contrast patterns of rapidly evolved population differentiation to thermal plasticity still present in the ancestral Mediterranean range.
(A) Sampling locations in the exotic and native range of Onthophagus taurus (also see Table S1). The native distribution range in Europe, based on DaSilva et al. (2016), is highlighted by the hatched area. Landmarks used to obtain morphometric measurements for fore tibiae (B) and hind wings (C). Semi-landmarks used to quantify wing shape are shown with smaller symbols.
Means of climatic adaptation are manifold, but life history responses to season length and temperature, and morphological as well as physiological adaptations linked to thermal regimes are documented most often (Blanckenhorn and Demont 2004; Sinclair et al. 2012; Roy et al. 2018; Flatt 2020). In large insects with long development times and one (or few) generations per year, body size and development time typically decrease with increasing latitude—an adaptive response to the continuous shortening of the reproductive season with latitude (Blanckenhorn and Demont 2004; Chown and Gaston 2010). However, in addition to changes in life history and physiology, the evolution of functional morphological traits allowing individuals to interact with their environment may also contribute to local adaptation. Buffering potentially detrimental environmental variation via behavioural mechanisms, the evolution of traits related to niche construction and (micro)habitat choice might represent an alternative way to respond to changing environments (Odling-Smee et al. 2013; Saltz and Nuzhdin 2014; Clark et al. 2020). As behavior typically shows high levels of plasticity and context-dependency, such traits may contribute to local adaptation during range expansions. Being a medium-sized coleopteran, O. taurus regulates its body temperature mainly by behavioral means, such as microhabitat choice (May 1979). In the adult, this is mostly driven by flight and digging underground. Because wing shape and size play a role in flight efficiency, they may also relate to local adaptation to varying climates (Stalker 1980; Frazier et al. 2008). In insects, it has for instance been argued that an increase in relative wing size enables more efficient dispersal at cool temperatures. This is because larger wings reduce the weight strained per unit wing area (“wing loading,” ρw), and thereby reduce power requirements for flight and increase lift production (see, e.g., Dudley 2002; Gilchrist and Huey 2004). Corresponding plastic as well as evolutionary variation in relative wing size and/or shape found in several intra- as well as interspecific comparisons might therefore relate to thermal adaptation (Stalker and Carson 1949; Azevedo et al. 1998; Rohner et al. 2018). In contrast to the adult stage, eggs, larvae, and pupae are confined to a subterranean brood chamber with limited or no capacity to evade unfavorable temperatures. To shield developing offspring from environmental variation, mothers evolved adaptive maternal behavioral plasticity in brood burying depth (Macagno et al. 2016, 2018). As tibia shape plays a key role in digging efficiency (Linz et al. 2019), population differentiation in foreleg morphology has been linked to temperature-dependent maternal care in terms of thermal niche construction (Macagno et al. 2016, 2018). Both tibia and wing size and shape may therefore be relevant for thermoregulatory capacity and thus contribute to local adaptation to novel climates.
Taking advantage of the well-documented natural history of O. taurus, we sought to assess the potential contribution of ancestral plasticity and the evolution of thermal plasticity during rapid range expansion and concurrent adaptation to novel climates. By rearing one population originating from the ancestral range and four exotic populations along a latitudinal cline at two temperatures under laboratory conditions, we investigate variation in major life history components, and use geometric morphometric tools to assess potential functional variation in the shapes of wings and tibiae. We expect Northern populations to (i) develop and grow relatively faster due to seasonal time constraints and, as a consequence, be smaller in final adult size (countergradient variation), (ii) survive better when reared at low temperatures, (iii) have disproportionately large wings (i.e., lower wing loading) and hence better dispersal capacity in cool environments. In addition, if plasticity in wing and tibia shape is an adaptive response to thermal regimes, we expect latitudinal variation in the exotic range to mirror ancestral thermal plasticity. We find evidence for rapid evolution of plasticity in major life history traits and morphology, and discuss the role of developmental plasticity during rapid range expansions/invasions and adaptation to climate change.
Materials and Methods
COMMON GARDEN REARING
To test for local adaptation among exotic populations, we first collected from four wild O. taurus populations along a latitudinal cline within the United States covering 1550 km (Fig. 1; Table S1). To contrast differentiation in the exotic range to ancestral patterns of plasticity, we also collected a wild Mediterranean population (Italy) to serve as a proxy for thermal plasticity within the native range. After having spent at least 2 weeks under standard laboratory conditions in mixed-sex colonies (to ensure sexual maturity), individual wild-caught females were transferred into rectangular oviposition containers (27cm × 8cm × 8cm) that were filled with a sterilized sand-soil mixture and topped off with 200 g defrosted cow dung (this rearing method has been applied previously in this species: Beckers et al. 2015; Macagno et al. 2016, 2018) and kept at 24°C. Reproductively active females dig vertical tunnels (typically 10–30 cm deep) immediately underneath the dung pat and, pulling dung form the surface, construct several compact spheres out of dung in which a single egg is laid. After 5 days, these so-called “brood balls” were sifted from the soil. Because morphological and life history traits are strongly dependent on larval nutrition and maternal investment in this species (Moczek 1998), we reared the F1 offspring in standardized, artificial brood balls as described previously (Shafiei et al. 2001). In brief, we opened all natural brood balls and transferred eggs into separate wells of a standard 12-well tissue culture plate provisioned with previously frozen, thoroughly homogenized cow dung (for details, see Shafiei et al. 2001; Moczek and Nijhout 2002b). The brood of each female was then evenly allocated to two temperature treatments that mimic local soil temperatures at a depth of 20 cm in the breeding season of the most southern (Florida; 27°C) and the most northern (Michigan; 19°C) population (based on data provided by Syngenta, the National Resources Conservation Service, and the National Oceanic and Atmospheric Administration). We used F1 offspring as opposed to later generations because European populations do not produce a second filial generation under laboratory conditions (Casasa and Moczek 2018). In order to keep the number of generations spent under laboratory conditions constant across all populations, we therefore only used offspring of wild-caught individuals.
F1 offspring reared in 12-well plates were inspected every 24 h to measure the duration of the three larval instars (L1–L3) and the pupal stage. To quantify larval growth rate, larvae were carefully removed from their artificial brood ball using featherweight forceps, and individually weighed on their second and fourth day into their third instar using a Mettler Toledo (AL54 Ohio, USA, d = 0.1 mg) scale. This is the instar in which larvae grow most, and this time window captures larvae during their nearly exponential growth phase and therefore best approximates the maximal intrinsic growth rate (Moczek and Nijhout 2002a). Following Tammaru and Esperk (2007), instantaneous relative larval growth rates (RGR) were computed as the ratio of larval mass at day 4 divided by mass at day 2. This ratio was then log-transformed and divided by the number of hours in between the two measurements (up to the nearest half hour). Upon eclosion and complete hardening, adult beetles were sacrificed and stored in 70% ethanol until dissection. Pronotum width, a standard measure of size in onthophagids (c.f., Emlen 1994; Moczek 2003), was measured as an estimate of overall body size.
To quantify tibia shape, we removed the protibia and photographed it using a digital camera (Scion, Frederick, MD, USA) mounted on a Leica MZ-16 stereomicroscope (Bannockburn, IL, USA). Hindwings were removed with micro scissors, embedded in glycerol, mounted on a glass slide, and photographed as described above. For tibiae, we acquired the coordinates for nine full landmarks (Fig. 1) using tpsDig2 (Rohlf 2009). For wings, we used 11 full landmarks and 24 semi-landmarks (Fig. 1). Landmark coordinates were subjected to a full Procrustes superimposition using the function gpagen() of the R-package geomorph version 3.1.1 (Adams et al. 2020). The position of semi-landmarks was optimized by minimizing bending energy. As an estimate of wing loading (ρw), we divided pronotum width3 by wing centroid size2 (i.e., body size per unit wing area; c.f., Starmer and Wolf 1989; Rohner et al. 2019). This index is proportional to the average pressure exerted on the surrounding air by the wings in non-accelerating flight (Dudley 2002).
STATISTICAL ANALYSIS
The effect of temperature, latitude, sex, and their interactions on life history, wing loading, and relative tibia size in the invasive range were analyzed with linear mixed models using the functions lmer() and lmerTest() as implemented in the R-packages lme4 version 1.1-21 (Bates et al. 2015) and lmerTest version 3.1-1 (Kuznetsova et al. 2017), respectively. Temperature, latitude, and sex were treated as fixed effects (the former two were mean-centered). To statistically account for the experimental design, we used the identity of the mother nested within population, as well as the 12-well plate an individual was reared in, as random intercepts. Non-significant interactions were removed. Survival probability was analyzed using a generalized linear mixed model with a binomial error distribution (using the glmer() function of the lme4 package). Because we could not assess the sex in all dead individuals, we could not test for sex-specific mortality. Similar models were used to test for thermal plasticity and sex differences in the Italian population, although without a latitude fixed effect and without nesting mothers within population (because only one population could be studied).
GEOMETRIC MORPHOMETRIC ANALYSIS
We averaged wing and tibia shape by family and tested for the effects of latitude, temperature, sex, and log centroid size, and all their interactions in the exotic range using a Procrustes ANOVA (as implemented in the R-package geomorph using the function procD.lm()). Non-significant interaction terms were removed. Similar models were run for the Italian population, although without fitting a latitude effect.
To compare latitudinal variation in tibia and wing shape in the exotic range with the plastic response to temperature in the ancestral range, we calculated correlations among the respective shape deformation vectors. To this end, we applied a Bayesian multivariate general linear mixed effects model utilizing Markov Chain Monte Carlo sampling (R-package MCMCglmm: Hadfield 2010). We first estimated the simultaneous main effects of latitude, temperature, and size on shape in populations in the exotic range. Because raw Procrustes-transformed coordinates are often prone to show high collinearity, we used their Principal Components (PCs) for further analysis. PCs are orthogonal and hence cause no computational issues related to multicollinearity. Because Procrustes superimposition results in a deficiency of four ranks, and semi landmarks contribute only about one dimension to a PC-space, we only fitted the first 14 PCs for tibia and 42 PCs for wing morphology, respectively. MCMCglmms were fitted separately for each sex using the identity of the mother within populations as random effects. The off-diagonal elements of the covariance matrix were set to zero (using the idh() function of MCMCglmm) given the orthogonal structure of the PCs at the individual level. Uninformative priors based on population identity were used for the residual and random covariance matrices (R, G1: v = 0.002). Models were run for 250,000 iterations using a thinning interval of 100, with the first 50,000 iterations being discarded (burn-in), resulting in 2000 uncorrelated posterior estimates stored for further analysis.
That is, the shape data (Y) was projected onto a vector in the direction of the plastic or latitudinal response (; as derived from a multivariate regression on family means). The regression score can then be used to visualize the strength and shape of the overall relationship between latitude or temperature and shape.
Results
CLINAL DIFFERENTIATION IN LIFE HISTORY AND DISPERSAL TRAITS
We sought to investigate the relationship between ancestral plasticity and subsequent population differentiation during the rapid colonization of North America by O. taurus. We found that the invasion of the United States is accompanied by rapid clinal differentiation across populations. Egg-to-adult development time decreased markedly with increasing latitudes, yet much more strongly so when O. taurus was reared at low temperatures (temperature × latitude interaction, Table 1A, Fig. 2). Because all individuals were reared under standardized environmental conditions, this clinal variation is most likely driven by genetic differentiation, indicating rapid evolution of plasticity that compensates for the drastic prolongation of development time driven by lower temperatures in the ancestral population (countergradient variation, also called “genetic compensation” c.f. Grether (2005); a common phenomenon in ectotherms: Blanckenhorn and Demont 2004; Conover et al. 2009; Sinclair et al. 2012). Analyzing development time separately for the three larval instars and the pupal stage revealed that this is (mostly) driven by the duration of the third instar while of little or no effect in the other stages (Fig. S1, Table S2). Mirroring the patterns of variation in development time, body size also decreased with latitude. This decline was also more pronounced at low temperatures (temperature × latitude interaction, Table 1).
Mixed models testing for sex differences, latitudinal clines, and thermal plasticity in exotic (a) and native (b) populations of O. taurus. Maternal identity and the 12-well plate an individual was reared in were used as random effect. When testing for latitudinal clines (a), maternal identity was further nested within populations. Non-significant interactions were removed. Three-way interactions were never significant
| a) Clinal variation in exotic range | ||||||||||
|---|---|---|---|---|---|---|---|---|---|---|
| Survival | Growth rate | Body size | Development time | Wing loading | ||||||
| P | P | P | P | P | ||||||
| Intercept | 103.06 | <0.001 | 2726.51 | <0.001 | 24577.55 | <0.001 | 29469.61 | <0.001 | 5264.59 | <0.001 |
| Latitude | 0.96 | 0.328 | 0.02 | 0.888 | 1.95 | 0.163 | 11.85 | 0.001 | 3.14 | 0.076 |
| Temperature | 3.86 | 0.049 | 7.47 | 0.006 | 4.08 | 0.043 | 1914.67 | <0.001 | 1.85 | 0.173 |
| Sex | 5.13 | 0.024 | 0.14 | 0.707 | 1.24 | 0.266 | 10.82 | 0.001 | ||
| Latitude × temperature | 7.63 | 0.006 | 9.42 | 0.002 | 7.87 | 0.005 | ||||
| Temperature × sex | 9.14 | 0.003 | 12.94 | <0.001 | ||||||
| b) Thermal plasticity in ancestral range | ||||||||||
| Survival | Growth rate | Body size | Development time | Wing loading | ||||||
| P | P | P | P | P | ||||||
| Intercept | 3.92 | 0.048 | 586.16 | <0.001 | 22724.2 | <0.001 | 9273.33 | <0.001 | 5078.3 | <0.001 |
| Temperature | 0.89 | 0.347 | 0.34 | 0.56 | 12.85 | <0.001 | 853.84 | <0.001 | 8.18 | 0.004 |
| Sex | 0.43 | 0.514 | 1.78 | 0.183 | 0.05 | 0.827 | 0.62 | 0.431 | ||
| a) Clinal variation in exotic range | ||||||||||
|---|---|---|---|---|---|---|---|---|---|---|
| Survival | Growth rate | Body size | Development time | Wing loading | ||||||
| P | P | P | P | P | ||||||
| Intercept | 103.06 | <0.001 | 2726.51 | <0.001 | 24577.55 | <0.001 | 29469.61 | <0.001 | 5264.59 | <0.001 |
| Latitude | 0.96 | 0.328 | 0.02 | 0.888 | 1.95 | 0.163 | 11.85 | 0.001 | 3.14 | 0.076 |
| Temperature | 3.86 | 0.049 | 7.47 | 0.006 | 4.08 | 0.043 | 1914.67 | <0.001 | 1.85 | 0.173 |
| Sex | 5.13 | 0.024 | 0.14 | 0.707 | 1.24 | 0.266 | 10.82 | 0.001 | ||
| Latitude × temperature | 7.63 | 0.006 | 9.42 | 0.002 | 7.87 | 0.005 | ||||
| Temperature × sex | 9.14 | 0.003 | 12.94 | <0.001 | ||||||
| b) Thermal plasticity in ancestral range | ||||||||||
| Survival | Growth rate | Body size | Development time | Wing loading | ||||||
| P | P | P | P | P | ||||||
| Intercept | 3.92 | 0.048 | 586.16 | <0.001 | 22724.2 | <0.001 | 9273.33 | <0.001 | 5078.3 | <0.001 |
| Temperature | 0.89 | 0.347 | 0.34 | 0.56 | 12.85 | <0.001 | 853.84 | <0.001 | 8.18 | 0.004 |
| Sex | 0.43 | 0.514 | 1.78 | 0.183 | 0.05 | 0.827 | 0.62 | 0.431 | ||
Mixed models testing for sex differences, latitudinal clines, and thermal plasticity in exotic (a) and native (b) populations of O. taurus. Maternal identity and the 12-well plate an individual was reared in were used as random effect. When testing for latitudinal clines (a), maternal identity was further nested within populations. Non-significant interactions were removed. Three-way interactions were never significant
| a) Clinal variation in exotic range | ||||||||||
|---|---|---|---|---|---|---|---|---|---|---|
| Survival | Growth rate | Body size | Development time | Wing loading | ||||||
| P | P | P | P | P | ||||||
| Intercept | 103.06 | <0.001 | 2726.51 | <0.001 | 24577.55 | <0.001 | 29469.61 | <0.001 | 5264.59 | <0.001 |
| Latitude | 0.96 | 0.328 | 0.02 | 0.888 | 1.95 | 0.163 | 11.85 | 0.001 | 3.14 | 0.076 |
| Temperature | 3.86 | 0.049 | 7.47 | 0.006 | 4.08 | 0.043 | 1914.67 | <0.001 | 1.85 | 0.173 |
| Sex | 5.13 | 0.024 | 0.14 | 0.707 | 1.24 | 0.266 | 10.82 | 0.001 | ||
| Latitude × temperature | 7.63 | 0.006 | 9.42 | 0.002 | 7.87 | 0.005 | ||||
| Temperature × sex | 9.14 | 0.003 | 12.94 | <0.001 | ||||||
| b) Thermal plasticity in ancestral range | ||||||||||
| Survival | Growth rate | Body size | Development time | Wing loading | ||||||
| P | P | P | P | P | ||||||
| Intercept | 3.92 | 0.048 | 586.16 | <0.001 | 22724.2 | <0.001 | 9273.33 | <0.001 | 5078.3 | <0.001 |
| Temperature | 0.89 | 0.347 | 0.34 | 0.56 | 12.85 | <0.001 | 853.84 | <0.001 | 8.18 | 0.004 |
| Sex | 0.43 | 0.514 | 1.78 | 0.183 | 0.05 | 0.827 | 0.62 | 0.431 | ||
| a) Clinal variation in exotic range | ||||||||||
|---|---|---|---|---|---|---|---|---|---|---|
| Survival | Growth rate | Body size | Development time | Wing loading | ||||||
| P | P | P | P | P | ||||||
| Intercept | 103.06 | <0.001 | 2726.51 | <0.001 | 24577.55 | <0.001 | 29469.61 | <0.001 | 5264.59 | <0.001 |
| Latitude | 0.96 | 0.328 | 0.02 | 0.888 | 1.95 | 0.163 | 11.85 | 0.001 | 3.14 | 0.076 |
| Temperature | 3.86 | 0.049 | 7.47 | 0.006 | 4.08 | 0.043 | 1914.67 | <0.001 | 1.85 | 0.173 |
| Sex | 5.13 | 0.024 | 0.14 | 0.707 | 1.24 | 0.266 | 10.82 | 0.001 | ||
| Latitude × temperature | 7.63 | 0.006 | 9.42 | 0.002 | 7.87 | 0.005 | ||||
| Temperature × sex | 9.14 | 0.003 | 12.94 | <0.001 | ||||||
| b) Thermal plasticity in ancestral range | ||||||||||
| Survival | Growth rate | Body size | Development time | Wing loading | ||||||
| P | P | P | P | P | ||||||
| Intercept | 3.92 | 0.048 | 586.16 | <0.001 | 22724.2 | <0.001 | 9273.33 | <0.001 | 5078.3 | <0.001 |
| Temperature | 0.89 | 0.347 | 0.34 | 0.56 | 12.85 | <0.001 | 853.84 | <0.001 | 8.18 | 0.004 |
| Sex | 0.43 | 0.514 | 1.78 | 0.183 | 0.05 | 0.827 | 0.62 | 0.431 | ||
Clinal variation and thermal plasticity for life history, size, and wing loading. (Jittered) family means and population means (±95 CIs) are given. Thermal plasticity in the ancestral range is indicated with black diamonds. North American populations are sorted by latitude (round symbols, colored from blue to green).
In contrast, juvenile survival did not vary with latitude. Animals developing at lower temperatures survived better on average, but this effect was not present in the ancestral population (Table 1). Growth rate did not show clinal variation either. Males grew faster than females and temperature had a positive effect on growth rate regardless of sex (Table 1). However, variation among individuals was considerable and thermal plasticity for growth rate was not significant in the ancestral population (although temperature still tended to affect growth rate in a similar manner, see Fig. 2). This indicates that other components of larval growth, such as the duration of the exponential growth phase or the weight lost during the prepupal stage must cause genetic differentiation (see, e.g., Rohner et al. 2017).
Wing loading was lower in animals reared at cool temperatures (Table 1), as is expected if, during development, individuals plastically adjust wing size to enhance flight capacity in the cold (Gilchrist and Huey 2004). This plastic response showed clinal variation in which the two southern populations showed only modest thermal plasticity, while the two northern populations (temperature × latitude interaction, Table 1, Fig. 2) showed much stronger responses. However, because wing loading shows static allometry, clinal variation, and plasticity in body size could introduce confounding effects. When accounting for body size by adding pronotum width as covariate, clinal variation in thermal plasticity indeed became non-significant (temperature × latitude interaction: Χ2(1) = 0.03, P = 0.858). The latitude effect, however, persisted (latitude main effect: Χ2(1) = 10.35, P = 0.001). Wing loading was also lower in females compared to males, but this effect was weaker among individuals reared at cool temperatures (sex × temperature interaction, Table 1). This sexual dimorphism and its thermal plasticity were not statistically significant in the ancestral population, although patterns point in similar directions (Table 1).
It is important to note that we here investigate phenotypic variation in offspring of wild-caught individuals. Even though the parental generation was kept under standardized laboratory conditions for several weeks before rearing began, and even though our larval rearing setup was designed to exclude potential influences of the most important maternal effects (brood ball mass, dung quality, rearing temperature), we still cannot fully exclude other nongenetic mechanisms that could contribute to the clinal variation observed here (due to limitations imposed by the beetles’ natural history (Casasa and Moczek 2018; see Discussion).
POPULATION DIFFERENTIATION IN FUNCTIONAL MORPHOLOGY
Next, we sought to investigate whether populations had also diverged in morphological traits related to flight and digging, specifically hind wings and foretibia, and the degree to which possible divergences along latitude were mirrored in plastic responses to rearing temperatures. We found that northern populations developed less round and more elongated wings (Fig. 3, Table 2). This clinal variation was mirrored by thermal plasticity in the exotic range as animals of both sexes exposed to cooler temperatures also developed longer and more slender wings (Figs. 3 and 4A). Corresponding vector correlations were moderate (r♂ = 0.37 [0.02, 0.56]; r♀ = 0.35 [0.01,0.57]) but significantly larger than what would be expected by chance under a uniform distribution. Similar plastic effects are evident in the ancestral population, yet despite high correlations between the ancestral and derived thermal plasticity (r♂ = 0.72 [0.44, 0.86], r♀ = 0.70 [0.67, 0.86]), there was no association between the former and clinal differentiation (Figs. 3 and 4A). Thermal plasticity, allometry, and clinal variation were very similar between the sexes, as evidenced by high correlations between sex-specific shape changes (Fig. 4B).
Deformations in male wing shape associated with the ancestral plastic response to temperature in the species’ native range (A), the derived thermal plasticity in the exotic range (B), as well as the latitudinal variation in the exotic range (C). Scatterplots show Drake and Klingenberg's regression score (s; see Materials and Methods). Deformation grids show shape deformations associated with hotter (red) and cooler (blue) environments. The change with latitude is magnified twofold, while thermal plasticity is magnified fourfold. Northern populations developed less roundish and more elongated wings. This clinal variation was mirrored by thermal plasticity in the exotic as well as the native range as animals of both sexes exposed to cooler temperatures also developed longer and more slender wings.
Procrustes ANOVAs (type III SS) testing for clinal variation, thermal plasticity, allometry, and sexual dimorphism in protibia and wing shape. Analyses are based on family means. Non-significant interactions were removed. Three-way interactions were never significant
| a) Exotic range | ||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Tibia shape | Wing shape | |||||||||||
| Df | SS | MS | F | Z | P | Df | SS | MS | F | Z | P | |
| Centroid size | 1 | 2.31 × 10−03 | 2.31 × 10−03 | 3.58 | 1.95 | 0.036 | 1 | 4.90 × 10−04 | 4.90 × 10−04 | 2.36 | 2.06 | 0.017 |
| Sex | 1 | 1.79 × 10−03 | 1.79 × 10−03 | 2.77 | 1.65 | 0.058 | 1 | 7.09 × 10−04 | 7.09 × 10−04 | 3.41 | 2.78 | 0.003 |
| Temperature | 1 | 1.11 × 10−02 | 1.11 × 10−02 | 17.11 | 3.77 | 0.001 | 1 | 1.58 × 10−03 | 1.58 × 10−03 | 7.62 | 4.26 | 0.001 |
| Latitude | 1 | 9.90 × 10−03 | 9.90 × 10−03 | 15.31 | 3.52 | 0.001 | 1 | 1.03 × 10−03 | 1.03 × 10−03 | 4.96 | 3.52 | 0.001 |
| Centroid size × sex | 1 | 4.21 × 10−03 | 4.21 × 10−03 | 6.52 | 2.62 | 0.004 | ||||||
| Sex × temperature | 1 | 4.26 × 10−03 | 4.26 × 10−03 | 6.59 | 2.58 | 0.005 | ||||||
| Centroid size × temperature | 1 | 7.27 × 10−04 | 7.27 × 10−04 | 3.50 | 2.83 | 0.003 | ||||||
| b) Native range | ||||||||||||
| Tibia shape | Wing shape | |||||||||||
| Df | SS | MS | F | Z | P | Df | SS | MS | F | Z | P | |
| Centroid size | 1 | 1.16 × 10−03 | 1.16 × 10−03 | 1.30 | 0.80 | 0.222 | 1 | 8.89 × 10−04 | 8.89 × 10−04 | 3.43 | 2.59 | 0.005 |
| Sex | 1 | 9.16 × 10−02 | 9.16 × 10−02 | 102.47 | 5.90 | 0.001 | 1 | 2.20 × 10−04 | 2.20 × 10−04 | 0.85 | −0.10 | 0.555 |
| Temperature | 1 | 4.01 × 10−03 | 4.01 × 10−03 | 4.48 | 2.28 | 0.017 | 1 | 1.66 × 10−03 | 1.66 × 10−03 | 6.38 | 3.79 | 0.001 |
| a) Exotic range | ||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Tibia shape | Wing shape | |||||||||||
| Df | SS | MS | F | Z | P | Df | SS | MS | F | Z | P | |
| Centroid size | 1 | 2.31 × 10−03 | 2.31 × 10−03 | 3.58 | 1.95 | 0.036 | 1 | 4.90 × 10−04 | 4.90 × 10−04 | 2.36 | 2.06 | 0.017 |
| Sex | 1 | 1.79 × 10−03 | 1.79 × 10−03 | 2.77 | 1.65 | 0.058 | 1 | 7.09 × 10−04 | 7.09 × 10−04 | 3.41 | 2.78 | 0.003 |
| Temperature | 1 | 1.11 × 10−02 | 1.11 × 10−02 | 17.11 | 3.77 | 0.001 | 1 | 1.58 × 10−03 | 1.58 × 10−03 | 7.62 | 4.26 | 0.001 |
| Latitude | 1 | 9.90 × 10−03 | 9.90 × 10−03 | 15.31 | 3.52 | 0.001 | 1 | 1.03 × 10−03 | 1.03 × 10−03 | 4.96 | 3.52 | 0.001 |
| Centroid size × sex | 1 | 4.21 × 10−03 | 4.21 × 10−03 | 6.52 | 2.62 | 0.004 | ||||||
| Sex × temperature | 1 | 4.26 × 10−03 | 4.26 × 10−03 | 6.59 | 2.58 | 0.005 | ||||||
| Centroid size × temperature | 1 | 7.27 × 10−04 | 7.27 × 10−04 | 3.50 | 2.83 | 0.003 | ||||||
| b) Native range | ||||||||||||
| Tibia shape | Wing shape | |||||||||||
| Df | SS | MS | F | Z | P | Df | SS | MS | F | Z | P | |
| Centroid size | 1 | 1.16 × 10−03 | 1.16 × 10−03 | 1.30 | 0.80 | 0.222 | 1 | 8.89 × 10−04 | 8.89 × 10−04 | 3.43 | 2.59 | 0.005 |
| Sex | 1 | 9.16 × 10−02 | 9.16 × 10−02 | 102.47 | 5.90 | 0.001 | 1 | 2.20 × 10−04 | 2.20 × 10−04 | 0.85 | −0.10 | 0.555 |
| Temperature | 1 | 4.01 × 10−03 | 4.01 × 10−03 | 4.48 | 2.28 | 0.017 | 1 | 1.66 × 10−03 | 1.66 × 10−03 | 6.38 | 3.79 | 0.001 |
Procrustes ANOVAs (type III SS) testing for clinal variation, thermal plasticity, allometry, and sexual dimorphism in protibia and wing shape. Analyses are based on family means. Non-significant interactions were removed. Three-way interactions were never significant
| a) Exotic range | ||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Tibia shape | Wing shape | |||||||||||
| Df | SS | MS | F | Z | P | Df | SS | MS | F | Z | P | |
| Centroid size | 1 | 2.31 × 10−03 | 2.31 × 10−03 | 3.58 | 1.95 | 0.036 | 1 | 4.90 × 10−04 | 4.90 × 10−04 | 2.36 | 2.06 | 0.017 |
| Sex | 1 | 1.79 × 10−03 | 1.79 × 10−03 | 2.77 | 1.65 | 0.058 | 1 | 7.09 × 10−04 | 7.09 × 10−04 | 3.41 | 2.78 | 0.003 |
| Temperature | 1 | 1.11 × 10−02 | 1.11 × 10−02 | 17.11 | 3.77 | 0.001 | 1 | 1.58 × 10−03 | 1.58 × 10−03 | 7.62 | 4.26 | 0.001 |
| Latitude | 1 | 9.90 × 10−03 | 9.90 × 10−03 | 15.31 | 3.52 | 0.001 | 1 | 1.03 × 10−03 | 1.03 × 10−03 | 4.96 | 3.52 | 0.001 |
| Centroid size × sex | 1 | 4.21 × 10−03 | 4.21 × 10−03 | 6.52 | 2.62 | 0.004 | ||||||
| Sex × temperature | 1 | 4.26 × 10−03 | 4.26 × 10−03 | 6.59 | 2.58 | 0.005 | ||||||
| Centroid size × temperature | 1 | 7.27 × 10−04 | 7.27 × 10−04 | 3.50 | 2.83 | 0.003 | ||||||
| b) Native range | ||||||||||||
| Tibia shape | Wing shape | |||||||||||
| Df | SS | MS | F | Z | P | Df | SS | MS | F | Z | P | |
| Centroid size | 1 | 1.16 × 10−03 | 1.16 × 10−03 | 1.30 | 0.80 | 0.222 | 1 | 8.89 × 10−04 | 8.89 × 10−04 | 3.43 | 2.59 | 0.005 |
| Sex | 1 | 9.16 × 10−02 | 9.16 × 10−02 | 102.47 | 5.90 | 0.001 | 1 | 2.20 × 10−04 | 2.20 × 10−04 | 0.85 | −0.10 | 0.555 |
| Temperature | 1 | 4.01 × 10−03 | 4.01 × 10−03 | 4.48 | 2.28 | 0.017 | 1 | 1.66 × 10−03 | 1.66 × 10−03 | 6.38 | 3.79 | 0.001 |
| a) Exotic range | ||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Tibia shape | Wing shape | |||||||||||
| Df | SS | MS | F | Z | P | Df | SS | MS | F | Z | P | |
| Centroid size | 1 | 2.31 × 10−03 | 2.31 × 10−03 | 3.58 | 1.95 | 0.036 | 1 | 4.90 × 10−04 | 4.90 × 10−04 | 2.36 | 2.06 | 0.017 |
| Sex | 1 | 1.79 × 10−03 | 1.79 × 10−03 | 2.77 | 1.65 | 0.058 | 1 | 7.09 × 10−04 | 7.09 × 10−04 | 3.41 | 2.78 | 0.003 |
| Temperature | 1 | 1.11 × 10−02 | 1.11 × 10−02 | 17.11 | 3.77 | 0.001 | 1 | 1.58 × 10−03 | 1.58 × 10−03 | 7.62 | 4.26 | 0.001 |
| Latitude | 1 | 9.90 × 10−03 | 9.90 × 10−03 | 15.31 | 3.52 | 0.001 | 1 | 1.03 × 10−03 | 1.03 × 10−03 | 4.96 | 3.52 | 0.001 |
| Centroid size × sex | 1 | 4.21 × 10−03 | 4.21 × 10−03 | 6.52 | 2.62 | 0.004 | ||||||
| Sex × temperature | 1 | 4.26 × 10−03 | 4.26 × 10−03 | 6.59 | 2.58 | 0.005 | ||||||
| Centroid size × temperature | 1 | 7.27 × 10−04 | 7.27 × 10−04 | 3.50 | 2.83 | 0.003 | ||||||
| b) Native range | ||||||||||||
| Tibia shape | Wing shape | |||||||||||
| Df | SS | MS | F | Z | P | Df | SS | MS | F | Z | P | |
| Centroid size | 1 | 1.16 × 10−03 | 1.16 × 10−03 | 1.30 | 0.80 | 0.222 | 1 | 8.89 × 10−04 | 8.89 × 10−04 | 3.43 | 2.59 | 0.005 |
| Sex | 1 | 9.16 × 10−02 | 9.16 × 10−02 | 102.47 | 5.90 | 0.001 | 1 | 2.20 × 10−04 | 2.20 × 10−04 | 0.85 | −0.10 | 0.555 |
| Temperature | 1 | 4.01 × 10−03 | 4.01 × 10−03 | 4.48 | 2.28 | 0.017 | 1 | 1.66 × 10−03 | 1.66 × 10−03 | 6.38 | 3.79 | 0.001 |
Correlation between shape change caused by thermal plasticity and latitudinal clinal differentiation (A), and similarity of the effects of thermal plasticity, size, and latitude between sexes (B). (A) The posterior distributions of vector correlations between thermal plasticity in the exotic range (=derived plasticity) and latitudinal population differentiation in tibia and wing shape indicate a moderate overlap in both sexes. That is, the shape changes caused by developmental temperature overlap somewhat with clinal variation. In contrast, ancestral patterns of plasticity do not align with latitudinal variation for wing shape but do so for tibia shape. However, only the correlation between clinal variation and derived thermal plasticity was greater in magnitude than what could be expected by chance (asterisks indicate significance based on the distribution of vector correlations drawn at random from a uniform distribution with the proper dimensionality (see Materials and Methods); **P < 0.01, ***P < 0.001). (B) The effects of thermal plasticity and clinal variation correlate strongly between sexes for tibia and wing shape, and significantly more so than what can be expected by chance. The covariation between size and shape however differs strongly between sexes for tibiae, while allometry is similar for wings.
Northern populations also evolved stockier fore tibiae with a relatively short base, whereas southern populations exhibited more arched and slender tibiae, with tibial teeth located more toward the distal region of the tibia (Fig. 5). Similar shape deformations were associated with thermal plasticity of both sexes in the ancestral as well as the exotic range (Figs. 4A and 5), with correlations ranging from r = 0.21 to r = 0.28. The distribution of posterior values derived from MCMCglmms indicate weak to moderate correlations between thermal plasticity and latitude effects for tibia shape. However, the probability that a vector with the same dimensionality drawn at random from a uniform distribution has a vector correlation equal to or larger than the observed mean posterior value was larger than the critical threshold of 0.05. Correlations of this magnitude therefore fall in the range that could be expected by chance. Females had much broader tibiae with stronger protrusions compared to males and showed distinctly different allometric variation (Fig. 5B, Table 2A). This contrasts to the effects of thermal plasticity and latitude that aligned closely between the sexes (Fig. 4B; despite a significant sex-by-temperature interaction, see Table 2A). Combined, these findings suggest rapid evolution of clinal variation and plasticity in the exotic range. However, ancestral plasticity does not seem to be a major predictor of subsequent genetic change. Instead, latitude-by-temperature interactions for development time and body size, as well as the alignment between the derived thermal plasticity in wing shape and the newly evolved wing cline suggest that plasticity has the potential for rapid evolution during range expansions.
Deformations in male tibia shape associated with the ancestral plastic response to temperature in the species’ native range (A), the derived thermal plasticity in the exotic range (B), as well as the latitudinal variation in the exotic range (C). Scatterplots show Drake and Klingenberg's regression score (see Materials and Methods). Deformation grids show shape deformations associated with hotter (red) and cooler (blue) environments. The change with latitude is magnified twofold, while thermal plasticity is magnified fourfold. Northern populations evolved stockier fore tibiae with a relatively short base, whereas southern populations exhibited more arched and slender tibiae, with tibial teeth located more toward the distal region of the tibia. Similar shape deformations were associated with thermal plasticity of both sexes in the ancestral as well as the exotic range.
Discussion
In this study, we sought to assess the role of latitudinal population differentiation, phenotypic plasticity, and population differentiation in plasticity in rapid climatic range expansion during the invasion of O. taurus in eastern North America over the past half century. Studying ancestral and derived populations led to four salient results. First, we find that rapid range expansion in the exotic range coincides with the establishment of clinal variation in development time and body size, likely adaptations to seasonal time constraints in cooler environments (Blanckenhorn and Demont 2004; Conover et al. 2009; see below). Second, northern populations have evolved disproportionately large wings, consistent with the hypothesis that cool temperatures favor the evolution of increased dispersal capacity at cool temperatures via a decrease in wing loading (Stalker 1980; Frazier et al. 2008). Third, clinal variation in development time and body size is to a large extent driven by population differentiation in plasticity along the cline, suggesting that the evolution of plasticity might have played a critical role in successful invasion and adaptation on ecological timescales. Fourth, O. taurus also evolved clinal variation in wing and tibia shape. While the wing shape cline significantly aligns with thermal plasticity, the association between thermal plasticity and latitudinal cline was weaker for tibia shape. Below, we discuss the implications of our results in the context of the adaptive value of clinal variation and the role of plasticity evolution in facilitating the ability of populations to withstand and rapidly adapt to novel climatic conditions.
The role of plasticity during rapid biological invasions in natural populations remains poorly understood. In this study, we documented rapid evolution of clinal variation and/or thermal plasticity for development time, body size, wing loading, and wing and tibia shape. Clinal variation in development time is likely adaptive as northern populations have to complete their reproductive cycle within a much shorter time window compared to their southern counterparts (the average low temperature exceeds 10°C from April to October in Florida, while the same is true only in June, July, and August in Michigan; www.usclimatedata.com). Because the lower mean temperatures in the north further prolong development times and northern populations presumably only have one generation per year, seasonal constraints are likely a major selective driver and elicit the evolution of countergradient variation (i.e., genetic compensation sensu Grether 2005) in development time (Blanckenhorn and Demont 2004; Conover et al. 2009). Decreasing body size with latitude then results as a side effect of selection on development time due to their close genetic and developmental correlation. However, in this study, clinal variation depended strongly on rearing temperature. This indicates that selection did not bring about the evolution of decreased mean development time and body size per se but instead led to the evolution of thermal plasticity. In O. taurus, the evolution of the plastic response, rather than the evolution of mean phenotypes, therefore, compensates for the prolonged development time in the north, thereby rescuing a phenotype critical to population establishment and persistence.
Clinal variation was not restricted to major life history traits but was also apparent in relative wing size. Variation in wing size and corresponding changes in wing loading have previously been linked to the power output of flight and in particular the capacity for take-off at low temperatures (Stalker 1980; Dudley 2002; Fraimout et al. 2018). If individuals with lower wing loading are better able to disperse at cooler temperatures, this may provide a competitive advantage not only in terms of increased foraging efficiency and quicker escape behavior, but also in terms of increased mobility that allows evading cool microhabitats on the fly. As most insects regulate their body temperature mainly by flight and microhabitat choice (May 1979; Dillon et al. 2009), clinal variation and thermal plasticity in relative wing size may play a critical role in adaptation to cooler climates (Gilchrist and Huey 2004). Mirroring patterns of local adaptation in other insects (Stalker and Carson 1949; Azevedo et al. 1998; Rohner et al. 2018, Rohner et al. 2019), northern O. taurus populations evolved reduced wing loading, predicted to yield enhanced lift at cool temperatures. Clinal variation was again more pronounced at cool rearing temperatures, suggesting latitudinal differentiation in thermal plasticity. However, this was confounded by correlated variation in overall size, and whether such plastic responses are indeed a direct effect remains unclear. Nevertheless, the evolution of wing-body size scaling relationships, and possibly their dependence on environmental variation, may be a critical contributor to local adaptation, especially to cooler habitats.
We also found clinal variation and thermal plasticity in tibia shape in that animals from low latitudes and those reared at warm temperatures developed more slender tibiae (Fig. 4B). Although the cline and plasticity at least superficially have similar phenotypic effects, the vector correlation was not significantly stronger than what could be expected by chance. Previous work documented rapid divergence in tibial shape between Italian O. taurus and its syntopic sister species, O. illyricus, as well as between exotic eastern United States and Western Australian O. taurus populations. In both cases more slender tibiae correlated with more shallow brood provisioning, whereas broader tibiae were associated with deeper digging (Macagno et al. 2016, Macagno et al. 2018). Subsequent work further identified differential plasticity among exotic populations in burial depth in response to ambient temperature (Macagno et al. 2018) and functionally related tibial morphology to burying efficiency (Linz et al. 2019). Thus, clinal variation in tibial shape and shape plasticity documented here may reflect adaptive responses to divergent, temperature-dependent differences in selection pressures on burial depth.
Wing shape also showed plastic as well as clinal variation. Thermal plasticity in wing shape in the native range, did not correlate with the newly evolved cline, suggesting that ancestral plasticity did not directly affect clinal differentiation. However, the correlation between derived plasticity in wing shape and the cline are larger than expected by chance. This suggests that thermal plasticity in wing morphology evolved in the exotic range and interacted with subsequent population differentiation, thereby possibly mediating local adaptation. Analogous patterns have been found for the wing shape of Drosophila melanogaster (Pitchers et al. 2013) and the sepsid fly Sepsis punctum (Rohner et al. 2019), where plastic shape changes are related to genetic differentiation along climatic gradients. However, these patterns are not ubiquitous as thermal plasticity was unrelated to population differentiation in ecologically similar species (as for instance in the sepsid Sepsis fulgens, Rohner et al. 2019, or the yellow dung fly Scathophaga stercoraria, Schäfer et al. 2018). The mechanisms underlying such patterns, for instance the potential contribution of genetic accommodation (sensu West-Eberhard 2003), warrant further scrutiny.
However, it is worth pointing out that genetic variation in dispersal ability and traits related to dispersal, such as wing loading and wing shape, may assume a particular role especially during range expansions (Phillips et al. 2010; Ochocki and Miller 2017). Thus, clinal variation as detected in this study may not be a result of selection for thermoregulatory behavior in local environments alone, but also, or perhaps predominantly, a result of selection for dispersal capacity itself. This is because the invasion front is expected to harbor individuals that are disproportionately good dispersers. Nevertheless, independent of whether variation in dispersal capacity is driven by migration, selection on dispersal in general, or thermoregulatory behavior in the cold, (temperature-dependent) dispersal capacity is likely to have contributed to range expansion. Still, future work is needed to investigate the role of drift and repeated bottlenecks throughout the range expansion, as well as the potential contribution of population structure and ongoing introgression. These and other approaches will open up promising avenues to uncover the population genetic causes and consequences alongside the genetic underpinnings of rapidly evolving clinal variation, such as the potential role of genomic inversions (Kapun et al. 2016), allele surfing (Moreau et al. 2011; Gralka et al. 2016), or release of cryptic genetic variation (e.g., Ledón-Rettig et al. 2010; Rohner et al. 2013) during local adaptation.
It is important to note that, by studying phenotypic variation of F1 individuals, the clinal variation observed subsumes genetic differentiation among families and populations, as well as parental and other extra-genetic sources of heritable variation (e.g., Bonduriansky and Day 2020). In fact, parental effects are well known to affect life history and morphology in O. taurus, yet all routes of extra-genetic inheritance documented so far relate to the quality and size of the brood ball, as well as brood ball burial depth (Moczek 1998; Hunt and Simmons 2000; Hunt and Simmons 2002; Beckers et al. 2015; Macagno et al. 2018). By rearing all individuals in highly standardized artificial broodballs from the earliest stages of juvenile development, such parental effects are accounted for. The clinal variation observed is therefore likely to be driven mostly by quantitative genetic differentiation. Nevertheless, we cannot completely exclude that other sources, such as epigenetic modifications, maternal manipulation of egg contents, or the transgenerational transmission of microbial endosymbionts, contribute to population differentiation and/or thermal plasticity.
The rate of change in climatic conditions that North American O. taurus populations experienced during their recent range expansion is very likely much greater than climatic changes driven by global environmental change currently occurring at any given location. Thus, the mere observation that a species can thrive despite such rapid climatic change, is good news in general. However, O. taurus is well-known for its rapid population differentiation in diverse fitness-related traits (Moczek 2003; Beckers et al. 2015; Casasa and Moczek 2018). Whether this is the case in other species remains unclear, in particular because not all O. taurus populations have exhibited the same rapid range expansion. For instance, Australian O. taurus purposefully introduced around the same time as their North American counterparts (Bornemissza 1976) have diverged from the ancestral Mediterranean population in diverse morphological, physiological, and life history traits to a similar degree as their eastern North American counterparts (although typically opposite in direction; Moczek 2003; Beckers et al. 2015; Casasa and Moczek 2015). Yet climatic niche differentiation remained rather modest compared to the lineage that invaded eastern North America (DaSilva et al. 2016). Whether this disparity is driven by demography, differences in the composition of the founder population, or other biotic and abiotic interactions remain unknown but opens up interesting avenues of future comparative research on the causes of differential invasion success in the same species.
Conclusions
The role of plasticity and its genetic differentiation in rapid adaptation is complex (Crispo 2008; Conover et al. 2009; Hendry 2015). Growing evidence suggests that the evolution of plasticity may be especially significant during rapid range expansions and concurrent adaptation to new climatic conditions (e.g., Kelly 2019). Our data show that Onthophagus taurus rapidly evolved clinal population differentiation and suggest that post-introduction evolution of plasticity contributed to a significant degree to the successful invasion of North America by this species. In particular, our work highlights plasticity-mediated clinal differentiation in development time and morphological traits related to thermoregulatory capacity as potentially critical contributors to invasion success. More generally, our work illustrates how invasive species can be utilized to assess the interplay between development and environment in shaping evolutionary trajectories.
AUTHOR CONTRIBUTIONS
P.T.R. and A.P.M. designed the study, P.T.R. collected and analyzed all data, and P.T.R. and A.P.M. wrote the manuscript.
ACKNOWLEDGMENTS
We would like to thank Levi Burdine, Kayla Copper, Jordan Crabtree, and Anna Macagno for taking care of laboratory cultures. We are also indebted to Erik Parker and Maura Bocci for collecting O. taurus in the field. In addition, we thank Ian Dworkin and two anonymous reviewers as well as the Editors for helpful and constructive comments on earlier versions of this manuscript. P.T.R. was supported by an Early Postdoc.Mobility fellowship by the Swiss National Science Foundation (grant no. P2ZHP3_184003). Additional support was provided by National Science Foundation grants IOS 1256689 and 1901680 to A.P.M. as well as grant 61369 from the John Templeton Foundation. The opinions, interpretations, conclusions, and recommendations are those of the authors and are not necessarily endorsed by the National Science Foundation, or the John Templeton Foundation.
DATA ARCHIVING
Dryad doi: https://doi.org/10.5061/dryad.x69p8czfv
LITERATURE CITED
Associate Editor: T. Flatt
Handling Editor: T. Chapman
