Daniela Polic1*
Yeşerin Yıldırım1†
Roger Vila2
Paulo Ricardo Ribeiro Cardoso2
Markus Franzén1
Anders Forsman1- 1Department of Biology and Environmental Science, Linnaeus University, Kalmar, Sweden
- 2Institut de Biologia Evolutiva (CSIC), Universitat Pompeu Fabra, Barcelona, Spain
Understanding how large-scale environmental variability may shape the distribution of phenotypic variation remains central to evolutionary biology. Across-species comparisons of trait variation alongside environmental gradients may offer valuable insights into how different species may respond to similar selective pressures. We conducted a comparative morphological study (>32° latitude and >47° longitude) on three closely related Argynnini butterfly species, Speyeria aglaja, Fabriciana adippe, and F. niobe. We measured wing size and coloration to assess (1) whether they respond similarly or differently to environmental factors (longitude, latitude, altitude, temperature, precipitation, solar radiation, wind speed); (2) if these factors correspond with those associated with the species’ genetic structure based on a previous study; and (3) whether correlations between phenotypic traits within individuals are species-specific. We found common and species-specific associations of climatic (precipitation, wind speed) and geographic (longitude, altitude) factors with the composite phenotypic variation. Wing size was associated with different variables in the studied species, while melanisation mainly increased in cooler regions in all species, suggesting that the need for temperature regulation is a strong selective pressure on melanisation. Wing size was associated with the genetic structure in all species, highlighting the functional importance of this trait. The environmental drivers associated with the phenotypic variation in S. aglaja and F. adippe were largely the same as those associated with their genetic structure, hinting at a genetic underpinning of the observed morphological variation due to local adaption. We report some distinct intraspecific trait correlations in S. aglaja and F. adippe, indicative of independent phenotypic integration. These phenotypes seem to be associated with protection against predators and thermal regulation in the respective habitats of both species, suggesting that similar selective pressures have resulted in the evolution of different trait combinations. Some of the inter-specific differences could be related to diverging niche breadths and dispersal capacities, exemplifying that the evolution of trait integration and spatial phenotypic differentiation may differ between closely related species with overlapping distribution ranges. Our findings highlight the importance of comparative assessments of variation, and demonstrate that the relative effects of drivers of variability may vary between sister species.
Introduction
There is a seemingly boundless variation in nature among species, populations, and individuals, and to find and explain patterns of biodiversity has been of long-standing interest in evolutionary biology (Darwin, 1859; Wiens, 2015; Rangel et al., 2018). By selecting for characteristics that increase ecological and evolutionary success of individuals, populations, and species, the large-scale variation in environmental conditions can influence the geographic distribution of trait variation across species (Violle et al., 2014). Large-scale comparisons of species’ variability in phenotypic characters along environmental gradients thus offer an opportunity to understand the selection pressures involved and to explore whether species with different life-histories and legacies show shared or unique responses. If species respond similarly to prevalent environmental conditions, phenotypic patterns might be spatially analogous in sympatric congeneric species as a result of convergent evolution or shared developmental plasticity responses (Losos, 2011; Stern, 2013; Forsman, 2015). Conversely, different spatial patterns and associations of phenotypes with environmental variables in distinct species might either point to species-specific adaptations or plasticity responses owing to differences in ecology and genetic make-up, or reflect important roles of neutral evolution associated with stochastic effects such as genetic drift and founder events (Lande, 1976, 1980; Huynen, 1996; Bahrndorff et al., 2006; Tinnert and Forsman, 2017; Vargas et al., 2017; Qu et al., 2021). Further, distinct phenotypic traits within a species do not develop and evolve independently of each other (Forsman, 2015). This is in part because correlational selection favors individuals with certain combinations of trait values over individuals with other trait value combinations, and this promotes the evolution of genetic correlations and phenotypic integration, aka complex phenotypes, with the optimal solutions depending on environmental conditions and selective regimes (Brodie, 1992; Forsman and Appelqvist, 1998; Forsman et al., 2002; Hagman and Forsman, 2003; McKinnon and Pierotti, 2010; Forsman, 2018; Dingemanse et al., 2020). Comparing variation in different phenotypic dimensions across hierarchical levels of organization, and their associations with environmental factors may thus allow for evaluating the interacting roles of external conditions and intrinsic species characteristics for population differentiation and evolutionary modifications of phenotypic integration.
Butterflies are well suited for assessing associations of phenotypes with environmental and geographic factors within and across species. Their wings consist of arrays of colored scales, come in numerous patterns, colors, shapes, and sizes, and have many important functions with ample opportunity for evolutionary modifications. They are important features used for flight, and variation in their shapes and color patterns are involved in inter- and intraspecific visual communication, temperature regulation, protection from UV-radiation, and predator avoidance (Betts and Wootton, 1988; Scoble, 1992; Brunton and Majerus, 1995; Tsai et al., 2020). There is a tremendous variation in wing phenotypes among species, among populations within species, between sexes, among individuals within populations (e.g., polymorphism), and even within individuals (e.g., between fore and hind wings, between dorsal and ventral sides).
Because of the diverse functional roles of butterfly wings, different selective pressures may have contributed to the immense variety in lepidopteran wing morphology. For example, the variation in the degree of melanism might depend on the ambient temperature and solar radiation, as melanisation is important for thermoregulation and protection against UV-light (Clench, 1966; Matute and Harris, 2013; Heidrich et al., 2018). Melanisation can also increase pathogen resistance in insects, such that darker individuals might be more effective in encapsulating parasitoids, and more resistant against viral and fungal infection (Wilson et al., 2001; True, 2003; Mikkola and Rantala, 2010). Green color scales may reduce predation by providing camouflage while perching on leaves or grass, and aid intraspecific communication, as in Callophrys rubi (Michielsen et al., 2010). Structural iridescent scales might act as protection against high UV-radiation, constitute a water repellent, and strengthen the wing structure (Potyrailo et al., 2007; Zheng et al., 2007; Fang et al., 2015), and they are also involved in signaling and communicating during species recognition and mate choice (Brunton and Majerus, 1995; Brunton, 1998; Kemp, 2008). Variation in wing size is also influenced by both biotic and abiotic factors. For example, butterfly populations in areas with longer vegetation periods might have larger wings resulting from a prolonged period of larval feeding and growth (Chown and Gaston, 2010; Zhang et al., 2019). Additional factors believed to influence the size of butterflies are temperature, solar insulation, and the presence of predators (Zeuss et al., 2017; Le Roy et al., 2019). Body size—in butterflies traditionally accounted for by wing size—has previously been associated with variation in flight capacity in insects (Nieminen et al., 1999; Öckinger et al., 2010; Kuussaari et al., 2014), and habitat fragmentation, connectivity, and other environmental stressors acting on mobility might therefore contribute to the diversification of wing length in butterflies. Most previous investigations have attempted to evaluate the contribution and responses to these potential sources of variation using spatial and temporal comparisons of single species. While such approaches can generate important insights, multispecies comparisons provide added value by allowing for the assessment of reproducibility and generality of research findings (Voelkl et al., 2020; Polic et al., 2022).
Fritillary butterflies (i.e., Nymphalidae butterflies with mostly checkered wing patterns) in the tribe Argynnini offer suitable model systems for studying and comparing associations of wing morphology with climatic and geographic factors within and among closely related species that are exposed to similar environmental conditions. Speyeria aglaja, Fabriciana adippe, and F. niobe are morphologically similar and have partly overlapping distribution ranges, yet they differ in their degree of ecological generalism, habitat preferences, host plant use, and mobility patterns (Forster and Wohlfahrt, 1955; Higgins and Riley, 1970; Komonen et al., 2004; Eliasson et al., 2005; Fric et al., 2005; Öckinger et al., 2006; Polic et al., 2020). A previously published large-scale comparative study of genetic structure in these three Argynnini species (Polic et al., 2022) reports on species-specific associations of geographical and climatic factors with genetic variation. We hypothesize that the same combinations of these environmental variables might also be linked to the large-scale morphological trait distributions. Additionally, individuals from regions that constituted far diverged genetic groups in the phylogeographic comparative study might also be morphologically distinct from other studied regions.
Here, we quantified, characterized, and compared phenotypic (viz. morphological) variation in these three Argynnini species based on data collected across a large spatial scale that spanned more than 32 degrees latitude and more than 47 degrees longitude (Figure 1), with the overarching aim to identify and compare eco-evolutionary drivers of spatial differentiation and phenotypic integration. To achieve this, we set out to answer the following specific questions. (i) Do different species respond similarly to environmental conditions, such that the spatial variations in wing color patterns and wing size, and their associations with geographic and climatic variables are parallel in the three studied species? Or is the spatial variation in phenotypic traits seen in one species independent of the spatial variation seen in another species? (ii) Is the variation in different phenotypic dimensions (wing size, degree of melanisation, and the amount of white and green color on the wings) independent or correlated within individuals (i.e., between fore and hind wings, or between dorsal and ventral sides of the wings), and are the correlations between different traits species-specific or shared between species? Finally, (iii) Do results from these large-scale morphological analyses parallel findings from a recent comparative phylogeographic study on essentially the same individuals (Polic et al., 2022), such that the same geographic and climatic variables are associated with genetic and phenotypic distribution of variation in the respective species?
Figure 1. Map of sampling localities for S. aglaja (turquoise), F. adippe (yellow), and F. niobe (red). The left parts of the butterflies show the dorsal wings and the right parts show the ventral wings.
Materials and methods
Study species
Speyeria aglaja is a rather widespread, and often abundant species, inhabiting wet and dry grasslands, woodland meadows, coastal areas, and scrublands (Eliasson et al., 2005; Ellis et al., 2010). The species’ distribution area comprises Northern Africa, Europe, Russia, and Southeast Asia, reaching as far north as Southern Lapland (Fox et al., 2011; Tolman and Lewington, 2012). The two Fabriciana species do not occur in the Maghreb and are not distributed as far north as S. aglaja, possibly due to a higher tolerance of colder temperatures in the latter species (Zimmermann et al., 2009; Ellis et al., 2010). Fabriciana adippe reaches as far north as central Sweden (Warren, 1995; Franzén and Johannesson, 2007; Fox et al., 2011; Tolman and Lewington, 2012), while F. niobe’s northern range limit is southern Sweden (van Swaay and Warren, 1999; Eliasson et al., 2005; Reinhardt and Bolz, 2011). Fabriciana adippe is relatively constrained to open areas within woodland landscapes, and is facing steep declines in some parts of northern Europe, such as the UK (Warren, 1995; Zimmermann et al., 2009; van Swaay et al., 2010; Fox et al., 2011; Ellis et al., 2019), while F. niobe prefers nutrient-poor warm micro-habitats within grasslands with a high abundance of Viola spp. from sea-level up to sub-alpine altitudes (Salz and Fartmann, 2009; Spitzer et al., 2009; Salz and Fartmann, 2017). During the last decades, F. niobe has strongly declined throughout central Europe, and is now considered endangered in some countries, e.g., Germany (van Swaay and Warren, 1999; Reinhardt and Bolz, 2011). Although the overall physical appearance in these Argynnini species is similar, there are some differences in wing patterns both between and within species. Fabriciana adippe, and F. niobe, for example, exhibit an analogous intraspecific color polymorphism, both species having one morph with silvery reflective dots on the ventral side of the hindwings (F. adippe f. adippe, and F. niobe f. niobe), and another morph where these dots are yellow (F. adippe f. cleodoxa, and F. niobe f. eris) (Eliasson et al., 2005; Tolman and Lewington, 2012). In both species, the silvery morphs are dominating in the northern parts of the distribution area, whereas the yellow morphs gradually replace the nominate forms in the southern ranges, especially in the Southern Balkans (Coutsis and Ghavalás, 2001; Tolman and Lewington, 2012). Contrastingly, S. aglaja is a monomorphic species, morphologically similar to the silvery morphs in the other two species. The dorsal wing sides in all three species are orange with dark (brown/black) spots and they display green scales on the underside of the hindwings (Eliasson et al., 2005; Tolman and Lewington, 2012). A large-scale comparative phylogeography study revealed that the genetic structure in these three species was linked to distinct combinations of climatic and geographic factors, indicative of species-specific responses to shared environments (Polic et al., 2022).
Study area and data collection
We selected 174 Argynnini individuals, viz. 53 S. aglaja, 85 F. adippe (including both morphs F. adippe f. adippe, and F. adippe f. cleodoxa), and 36 F. niobe (including both morphs F. niobe f. niobe, and F. niobe f. eris) from butterfly collections at the Institut de Biologia Evolutiva (IBE), CSIC-Universitat Pompeu Fabra in Barcelona and the Department of Biology and Environmental Science, Linnaeus University in Kalmar, covering a wide range of locations within their European (and Moroccan for S. aglaja, as the other two species do not inhabit the Maghreb) distribution areas (33.542–65.856° latitude, and − 8.227–38.875° longitude, Figure 1). Sampling collection took place between 2008 and 2017 and coordinates for each locality were recorded directly in the field. After capture, the wings were either directly cut from the dead specimens and dried, or the whole specimens were dried and wings were cut just before taking pictures for morphological analyses of coloration and size (see below).
Quantification of wing size and color patterns
The butterfly wings were analysed according to coloration and size. For each individual, we placed the ventral side of the left forewing and hindwing and the dorsal side of the right forewing and hindwing (4 wing parts) on a standard white paper sheet, together with a color card and a ruler (Figure 2). Pictures were then taken with a Nikon D3000 camera with an 18–55 mm lens without flash using the following camera settings: focal length 35 mm, ISO 400, shutter speed 1/200 s, aperture f/4.5. Each picture was taken in the exact same position with 34 cm distance between the camera lens and the wings. Light was provided by a lamp with a Compact Fluorescent Bulb (FPL 18 W 5400 K) in a fixed position. The same person took all photographs in the IBE laboratory in Barcelona. The pictures were analysed according to wing coloration and wing size using the image analysis software GIMP 2 (The GIMP Development Team, 2019). For each respective wing and side, the pixels containing a specific color were selected using the select by color tool with a threshold of 15. Wing colors were assigned into separate categories of interest, “dark brown” (#4a3822), “black” (#0a0907), “white” (#f2f4e7), “yellow” (#c5a035), “light green” (#7b7f20), and “dark green” (#565d17), and the number of pixels covered by each respective color was measured. By measuring the total amount of colored pixels in each wing, the percentage of each color in the wing and respective side was calculated. The degree of melanism in each individual was defined as the total percentage of black and dark brown pixels of both dorsal and ventral side of both fore and hind wing. Dorsal and ventral melanism were defined as the total percentages of black and dark brown pixels of dorsal and ventral sides, respectively, of both fore-and hindwing. Similarly, we calculated forewing and hindwing melanism. We added the percentages of light green and dark green in the underside of the hindwing and named the variable green. White hereafter refers to the percentage of white in the underside of the hindwing (silvery white spots only occur in the underside of the hindwings in S. aglaja, F. adippe f. adippe, and F. niobe f. niobe). Forewing size was measured as the linear distance from the base to the apex of the forewing, and hindwing size from the hindwing base to the tip of the median vein M3 (Comstock and Needham, 1898) by using the ruler tool in GIMP 2.
Figure 2. Placement of wings for size and color measurements, here in F. adippe f. adippe. Forewing size was measured as the straight-line distance from the base to the apex of the forewing and hindwing size from the hindwing base to the tip of the median vein M3.
Geographic locations and climatic conditions
For each sampling location, altitude was extracted from raster data attained from U.S. Geological Survey Earth explorer using QGIS 2.18.25 (QGIS.org, 2022). In addition, we obtained climate data (mean annual temperature, mean summer temperature, mean annual solar radiation, mean summer solar radiation, mean annual precipitation, mean summer precipitation, and mean summer wind speed) from WorldClim—Global Climate Data (Fick and Hijmans, 2017) for each sampling point. The R package raster 3.5 (Hijmans et al., 2021) was utilized to import these data as raster layers, and the function extract to connect the climatic values to each sampling location. We chose these variables as they have previously been reported to influence wing patterns in Lepidoptera (Sullivan and Miller, 2007; Nygren et al., 2008; Heidrich et al., 2018; Zaman et al., 2019). Please see the section “Data analysis” for specifications of which environmental variables were used for each analysis.
Data analysis
To justify the use of linear regression models, we visualized our data by plotting all possible combinations of phenotypic values with geographic and environmental variables. Since the associations shown in the plots followed rather linear relationships, we are confident that all variables should enter the equations linearly (Supplementary Figures S1–S3).
Multivariate analyses of overall phenotypic variation and its association with environmental factors
In addition to the univariate statistical approach outlined below, in which the different phenotypic traits were analysed one by one to test for differences between species and evaluate associations with the geographical and climatic environmental conditions, we also used a multivariate statistical approach. To this end, we used a principal component analysis (PCA, a dimension reducing approach) to describe the total phenotypic variation in all the original phenotypic traits in two dimensions (principal component axes, PC1 and PC2). This was first done for all three species pooled, and the results visualized in a biplot, to evaluate whether composite phenotypic variation overlapped or differed between the three species. To also formally evaluate whether the composite phenotypic variation differed markedly according to species or sex we performed a multivariate analysis of variance (MANOVA) using PC1 and PC2 as dependent variables, and species, sex, and their interaction as explanatory factors. Next, we performed a separate PCA for each species, and used the obtained PC-scores for each species as input to evaluate associations of overall phenotypic variation with environmental factors. To this end, we created biplots according to redundancy analysis (RDA), a constrained version of PCA, of phenotypic data to visualize associations of geographic and climatic variables with the phenotypic structuring in each studied Argynnini species.
Analysing spatial variation in wing size
We performed a three-way ANOVA to test for differences in wing size between species, sex, and wing types (forewing, hindwing). We also grouped individuals in each species into geographic regions, representing different European (and Moroccan for S. aglaja) climatic, ecological, and biogeographic zones, and, in some cases, glacial refugia during the Last Glacial Maximum (LGM, 26,500–20,000 years ago). We selected these regions to detect potential divergences in wing phenotypes, which are not only associated with environmental or geographic factors following a continuous gradient, but rather caused by large-scale landscape characteristics and/or the biogeographic history of these areas. These regions were also used in the comparative large-scale phylogeography study on essentially the same individuals (Polic et al., 2022). Depending on where the individuals in each species were collected, these regions included Morocco (only for S. aglaja, as the other two species do not occur in the Maghreb), Iberia, Central Europe, Eastern Europe, Sicily, Southern Italy, Southern Balkans, and Scandinavia. For each species, we performed 2-way ANOVA to model potential associations of wing type (forewing, hindwing), and geographic region, and used a post-hoc Tukey test to assess which regions differed significantly from each other. Then, we used linear mixed regression models, separately for each species, to explore if wing size was associated with geographic location (longitude, latitude, altitude), mean annual temperature, mean annual solar radiation, mean summer precipitation, mean summer wind speed, sex, or color morph (only for the polymorphic species F. adippe and F. niobe). We used mean annual temperature and mean annual solar radiation, as these variables are associated with the length of the vegetation period in temperate climates (Piao et al., 2006), which correlates with the duration of the larval food stage, thereby potentially affecting adult size in butterflies. Mean summer precipitation and mean summer wind speed were used, as these variables might affect mobility patterns during the flight period (i.e., summer), which might lead to adaptations in wing architecture and size. Wing type (forewing or hindwing) was included as a dummy variable, and individual identity was entered as a random factor to account for independence of repeated wing measurements of each individual. We also tested for interaction effects of sex, wing type and geographic and environmental variables, but kept only the significant ones (p < 0.05) to avoid problems associated with over-parametrization of the statistical models.
Analysing spatial variation in color patterns
We used a four-way ANOVA to test for differences in the degree of melanism between species, sex, wing types, and wing sides. Further, we performed ANOVA separately for each species to test for associations of melanism and geographic regions, equivalently to our models for wing size. Thereafter, we performed mixed linear models, separately for each species, to test for associations of the degree of melanism with geographic location (longitude, latitude, altitude), mean summer temperature, mean summer solar radiation, mean summer precipitation, mean summer wind speed, sex, and color morph (only for F. adippe and F. niobe). We used mean summer temperature and mean summer solar radiation, since the degree of melanism in the wings might be influenced by these variables during the butterflies’ flight period (i.e., summer), when darker wings might be beneficial for temperature regulation and protection from UV radiation. We included wing type (forewing/hindwing) and wing side (dorsal/ventral) as separate dummy variables, and the individual identity was included as a random factor to account for repeated measurements per individual. We also tested for interaction effects of sex, wing side, and wing type with climatic and geographic variables, but we kept and report only those interactions that were significant (p < 0.05).
We performed ANOVA and linear regressions to model associations of the amount of green on the ventral side of the hindwings (since green is only present there) with species identity, sex, the above-mentioned geographic regions, and climatic and geographic variables, respectively. Similarly, we analysed associations of the amount of white on the ventral side of the hindwings (white is only present in this part of the wings) with these variables. In the polymorphic species F. adippe and F. niobe, white spots only occur in the adippe and niobe morphs, respectively. Because of the small sample size in F. niobe f. niobe (8 individuals), this analysis was only performed for F. adippe f. adippe and S. aglaja.
Analysing phenotypic integration
To study phenotypic integration and evaluate whether evolutionary modifications and the expressions of different phenotypic dimensions are independent or not, we performed pairwise correlation analyses separately for each species. We evaluated all possible pairwise associations (viz. total melanism, dorsal melanism, ventral melanism, forewing melanism, hindwing melanism, forewing length, hindwing length, green in the ventral hindwing, and white in the ventral hindwing).
To determine whether phenotypic integration between wing size and the studied species of wing coloration were similar or different in the three studied species, we used the Mantel test (Mantel, 1967) to quantify the correspondence between the pairwise correlation matrices obtained for each species. To assess all possible comparisons, we performed three separate Mantel tests. With this approach, we evaluated whether pairs of phenotypic traits that were highly positively or weakly negatively correlated in one species were also highly positively or weakly negatively correlated in the other species. The statistical significance of the Mantel correlations was assessed using 999 permutations.
Associations of genetic structure with morphological variables
To assess associations of genetic structure with morphological variables (total melanism, forewing length, hindwing length, green, and white), we conducted distance-based redundancy analyses (db-RDA), a constrained version of a PCA (Legendre and Anderson, 1999), for each species separately. The genomic data was derived via double digest restriction-site associated DNA sequencing (ddRADseq). Details about DNA extraction and library preparation can be found in Polic et al. (2022). We generated Manhattan dissimilarity matrices for each species using pairwise distances between all individuals based on the proportion of shared alleles, which was used as input for the db-RDA. To build db-RDA models, we used the capscale function in the R package vegan 2.5 (Oksanen et al., 2019) to test which morphological variables contribute the most to the PCoA ordination. Using the function anova.cca with 999 permutations, we determined the overall statistical significance of the model and the explanatory variables. We used Holm’s correction for adjusting the p-values for multiple testing (Holm, 1979). We used longitude and latitude as a condition in the model to account for spatial autocorrelation between individuals, through which we partialled out these terms from the test prior to analysing the other predictor variables.
Model selection
For all regression and RDA models, we applied a backward elimination process, viz. we started with the most complex model and removed one variable at a time, starting with the least significant one. Non-significant interactions were removed before non-significant main effects. If non-significant main effects were included in a significant interaction, the main effect was kept in the model. The final most informative model included the highest number of significant predictor variables. The Akaike Information Criterion (AIC) was used for model comparisons for univariate analyses, and adjusted R2 was used for RDA model selection. To check for multi-collinearity, we calculated the variance inflation factor (VIF) for our explanatory variables, which were, depending on the analysis, mean annual temperature, mean summer temperature, mean annual solar radiation, mean summer solar radiation, mean annual precipitation, mean summer precipitation, mean summer wind speed, altitude, latitude, longitude, color morph (only for F. adippe and F. niobe), total melanism, forewing length, hindwing length, green, and white. The presence of multi-collinearity indicates that two or more predictor variables are highly correlated, and thus, the impact of one predictor on the response variable might be inaccurately estimated and not informative (James et al., 2013; Bruce and Bruce, 2017). Variables with a VIF > 10 were therefore excluded from the final model as this implies a potentially problematic level of collinearity (James et al., 2013). When we analysed all three species together, and when F. adippe and F. niobe were analysed separately, the variable green (underside of hindwing) was square root transformed to achieve normal distribution. The same was done for the variable altitude for S. aglaja and F. adippe, and for the variable melanism in all regression analyses. All statistical analyses were performed in R 3.6.1 (R Core Team, 2020).
Results
Associations of overall multivariate phenotypic variation with geography and climatic conditions
Principal component analysis on all three species pooled illustrated that the overall phenotypic variation based on the studied morphological variables was similar in the studied species (Figure 3). The results from multivariate analysis of variance suggest that composite phenotypic variation (as estimated using PC1 and PC2 in combination) did not vary to an important degree according to species or sex (MANOVA: effect of species, Pillai’s Trace = 0.04, F4,306 = 1.78; p = 0.13; effect of sex: Pillai’s Trace = 0.01, F2,152 = 0.92, p = 0.40; effect of species*sex: Pillai’s Trace = 0.07, F4,306 = 2.79, p = 0.03). Performing RDA based on PCA scores for each species separately to link the overall phenotypic variation to geographic and climatic variables revealed some similarities and some differences between the species (Figure 4). For S. aglaja, the final selected model included longitude and mean summer wind speed (adjusted R2 = 0.08), of which longitude was significantly associated with the overall phenotypic variation (F1,44 = 4.8, p < 0.01). For F. adippe, the final selected model included longitude, altitude, mean annual temperature, and mean summer precipitation (adjusted R2 = 0.22), of which longitude (F1,72 = 8.2, p < 0.001), altitude (F1,72 = 9.6, p < 0.001), and mean summer precipitation (F1,72 = 4.9, p < 0.01) were significantly associated with the general phenotypic variation. For F. niobe, the final selected model included mean annual temperature, mean summer precipitation, and mean summer wind speed (adjusted R2 = 0.19), of which mean summer precipitation (F1,30 = 4.8, p = 0.01), and mean summer wind speed (F1,30 = 3.5, p = 0.04) were significantly associated with the composite phenotypic structure.
Figure 3. Biplot illustrating the composite phenotypic differentiation between species resulting from PCA, indicated by dots colored according to species identity. The diamonds show the mean values of PC1 and PC2 scores and the error bars show standard deviations for each species.
Figure 4. Biplots according to RDA based on the PC-scores resulting from PCA for each species separately. Colors represent geographic regions where butterflies were collected and arrows indicate the direction of association of climatic and geographic variables with the overall phenotypic variation.
Wing size in relation to species, sex, geography, and climatic conditions
Wing size differed significantly between the studied species, with F. adippe being largest, S. aglaja intermediate and F. niobe smallest, but the relative size difference between fore-and hindwings did not vary between species. Females were significantly larger than males, which was most pronounced in F. niobe, and in the hindwings (ANOVA, main effect of species identity, F2, 329 = 21.08, p < 0.001; main effect of sex, F1, 329 = 71.18, p < 0.001; main effect of wing type, F1, 329 = 621.71, p < 0.001; interaction effect of species identity and sex, F2, 329 = 4.46, p = 0.01; interaction effect of sex and wing type, F1, 329 = 4.11, p = 0.04; interaction effect of species identity and wing type, F2, 329 = 1.1, p = 0.33, Table 1; Figure 5). According to results from ANOVA and linear mixed models for each species, the large-scale spatial variation in wing size and its association with environmental factors was largely species-specific (details in the Supplementary material, and the statistical results and their visualization in Tables 2, 3; Supplementary Table S1, Supplementary Figure S4).
Table 1. Mean values of phenotypic variables included in the analyses with standard deviations in brackets for S. aglaja, F. adippe, and F. niobe.
Figure 5. (A) Boxplot showing the significant differences in wing size between species, sexes and wing types. There was a significant interaction effect of sex with species identity and wing type, but the interaction of species identity and wing type was not significant. (B) Boxplot showing the significant 3-way interaction of species identity, wing type, and wing side on melanism (square root-transformed). Sex did not have an effect on the degree of melanism.
Table 2. Associations of wing size with geographic regions and wing type for S. aglaja, F. adippe, and F. niobe based on results from two-way ANOVA.
Table 3. Overview of geographic and climatic factors associated with phenotypic variables and interaction effects in the three studied fritillary species according to linear (mixed) regression models.
Dark coloration (melanism) in relation to species, sex, geography, and climatic conditions
The degree of wing melanisation varied in a complex interactive manner among species, and according to wing type (between the fore and hind wings) and wing side (between the dorsal and ventral sides), as evidenced by a significant three-way interaction between species, wing type and wing side. Sex was excluded from the analysis, as it did not have a significant effect (main or interaction) on the degree of melanism. In all three species, the dorsal wing side was darker than the ventral side. Dorsal hindwings were darker than dorsal forewings in all species, while ventral forewings were darker than ventral hindwings in S. aglaja and F. adippe, but not in F. niobe (3-way ANOVA, main effect of species identity, F2, 655 = 2.78, p = 0.06; main effect of wing type, F1, 655 = 6.18, p = 0.01; main effect of wing side, F1, 655 = 1369.44, p < 0.001; interaction effect of species identity and wing type, F2, 655 = 5.5, p < 0.01; interaction effect of species identity and wing side, F2, 655 = 14.26, p < 0.001; interaction effect of wing type and wing side, F1, 655 = 122.34, p < 0.001; interaction effect of species identity, wing type, and wing side, F2, 655 = 7.95, p < 0.001; Table 1; Figure 5). According to results from ANOVA and linear mixed models for each species, the large-scale spatial variation in dark coloration was mainly associated with aspects of the thermal environment (details in the Supplementary material, and the statistical results and their visualization in Tables 3, 4; Supplementary Table S2, Supplementary Figure S5).
Table 4. Associations of the degree of melanism with geographic regions, wing type and wing side for S. aglaja, F. adippe, and F. niobe based on results from 3-way ANOVA.
Amount of white on the ventral hindwings
The amount of white in ventral hindwings (only present there) differed significantly between F. adippe f. adippe, F. niobe f. niobe, and S. aglaja, while sex did not affect this trait (ANOVA, effect of species identity, F2, 111 = 7.16, p < 0.01; effect of sex, F1, 111 = 2.46, p = 0.12), i.e., in both F. adippe f. adippe (p < 0.01) and S. aglaja (p = 0.04) the amount of white was higher than in F. niobe f. niobe. Results from ANOVA and linear regression analyses for each species revealed that the associations of white wing coloration with geography and climatic conditions were species-specific (details in the Supplementary material, and the statistical results in Table 3 and Supplementary Table S3).
Amount of green on the ventral hindwings
The amount of green in the ventral hindwings (only present there) differed significantly between the three species, while sex did not have an effect on this trait (ANOVA, effect of species identity, F2, 160 = 27.11, p < 0.001; effect of sex, F1, 160 = 1.69, p = 0.2). Speyeria aglaja showed a significantly higher amount of green than F. adippe (p < 0.001) and F. niobe (p < 0.001), but there was no difference between F. adippe and F. niobe. ANOVA and regression analyses for each species suggested that the associations of green wing coloration with geography and climatic conditions were species-specific (details in the Supplementary material, and the statistical results in Table 3 and Supplementary Table S4).
Pairwise correlations between traits—phenotypic integration
The analyses of pairwise correlations between traits indicate that phenotypic integration was both general and species-specific. In all three species, dorsal, ventral, forewing, and hindwing melanism were positively correlated, although some of these associations were not significant in S. aglaja and F. adippe (Figure 6, for correlation coefficients please see Supplementary Table S5). Forewing and hindwing length correlated positively and significantly with each other in all species. In S. aglaja, the amount of green in the ventral hindwing was significantly and positively associated with both forewing and hindwing length, while it was significantly and negatively associated with ventral melanism and the amount of white in the ventral hindwing (Figure 6). Dorsal melanism was significantly and positively correlated with hindwing length, while there was a significant and positive correlation of forewing melanism with forewing and hindwing length in S. aglaja. In F. adippe, the amount of green in the ventral hindwing was significantly and negatively associated with ventral melanism, forewing melanism and total melanism (Figure 6). Forewing length was significantly and positively associated with the amount of white in the ventral hindwing in F. adippe. While there were some apparent differences in phenotypic integration between species as described above, results from Mantel correlation analyses based on comparisons of the matrices of estimated pairwise trait correlations for each species suggested significant and positive correlations between the three matrices (Mantel test, S. aglaja – F. adippe, r = 0.58, p < 0.01, S. aglaja – F. niobe, r = 0.47, p = 0.01, F. adippe – F. niobe, r = 0.45, p = 0.03).
Figure 6. Correlation matrices showing pairwise correlations in all analysed phenotypic variables for S. aglaja, F. adippe, and F. niobe. The colors show the correlation coefficients with blue illustrating positive correlations and red negative correlations. The larger and darker the circles, the stronger the correlations. The significance codes in the circles can be interpreted as follows: ***p < 0.001, **p < 0.01, *p < 0.05.
Associations of genetic structure with morphological variables
In S. aglaja and F. adippe (but essentially not in F. niobe), the phenotypic trait variation was largely associated with similar environmental variables that were previously found to be associated with their genetic structure (Polic et al., 2022). Wing size and melanism were associated with the genetic structure in all three species or in both Fabriciana species, respectively, hinting at a strong genetic signal and emphasizing the functional importance of these traits, e.g., for flight performance and temperature regulation. For S. aglaja, the final selected model included total melanism, forewing length, hindwing length, white, and green (adjusted R2 = 0.02). The genetic structure was significantly associated with forewing length (F1,30 = 1.42, p = 0.01). For F. adippe, the final selected model included color morph, total melanism, forewing length, hindwing length, and green (adjusted R2 = 0.02). The genetic structure was significantly associated with color morph (F1,60 = 1.27, p = 0.01), total melanism (F1,60 = 1.19, p = 0.05), forewing length (F1,60 = 1.31, p = 0.01), hindwing length (F1,60 = 1.38, p < 0.01), and green (F1,60 = 1.17, p = 0.05). For F. niobe, the final selected model included color morph, total melanism, forewing length, and hindwing length (adjusted R2 = 0.03). The genetic structure was significantly associated with total melanism (F1,22 = 1.25, p = 0.04), and forewing length (F1,22 = 1.34, p = 0.02).
Discussion
In this study, we analysed how phenotypic variation (wing size, wing melanisation, and the proportion of green and white coloration in the ventral hindwings) varied between species and how it was associated with climatic and geographic factors to identify and compare eco-evolutionary drivers of large-scale spatial differentiation and phenotypic integration in three Argynnini species, viz. S. aglaja, F. adippe and F. niobe. Notably, the results of our multivariate analyses show that the three species are phenotypically largely overlapping and nearly indistinguishable based on the morphological variables studied here. This overlap is in sharp contrast to the results of a recent comparison of genetic structure (based on analyses of ddRADseq data for the same set of individuals as used in the present study) which showed that these three species are genetically clearly distinct (cf Figure 3 in this paper with Figure 3F in Polic et al., 2022). While our analyses of drivers of phenotypic variation revealed some parallel trends between species (e.g., the composite phenotypic variation was associated with longitude in S. aglaja and F. adippe, and with precipitation in F. adippe and F. niobe), most associations of wing size and different aspects of wing coloration with climatic and geographic factors were species-specific, showcasing that even closely related species with similar distribution ranges may respond differently to common environments. We hypothesized that the combinations of environmental variables that shape the large-scale morphological trait distributions in these three species might be similar to those previously reported to be associated with the genetic structure (Polic et al., 2022). For S. aglaja and F. adippe, this hypothesis was largely qualitatively supported, while the phenotypic trait variation in F. niobe was mainly associated with different factors than the genetic structure. Wing size and melanism stood out as morphological variables that were associated with the genetic structure in all three or the two Fabriciana species, respectively. Further, we report some distinct intraspecific trait correlations for S. aglaja and F. adippe, suggesting that phenotypic integration is in part species-specific.
Large-scale patterns of spatial variation in wing size were mostly incongruent among species
Assessing and comparing large-scale associations of phenotypic characters with climatic and geographic variables revealed mostly divergent patterns between the three studied fritillary species. For example, only altitude was associated with wing size in more than one species, while the associations of wing size with the other factors were species-specific. Moroccan S. aglaja and Iberian F. adippe individuals were on average larger than from other regions, which coincides with the findings of a large-scale comparative phylogeographic study on the essentially same set of individuals, in which Moroccan S. aglaja and Iberian F. adippe samples appeared to be genetically distinct from other European sampling locations (Polic et al., 2022). That these regions constituted important glacial refugia during the LGM, when populations retracted to warmer areas during the ice ages, is still apparent in the phylogeography of many species (e.g., Cheddadi et al., 2006; Ferrero et al., 2011). Today, southern regions within our study area have longer vegetation periods and therefore offer more time for larval development, which may translate into larger adult butterflies. Further, the Strait of Gibraltar and the Pyrenees, respectively, might act as barriers to dispersal for both species, which might increase differences in wing size and genetic structure between individuals from these regions and the rest of the study area.
At higher altitudes, F. adippe individuals had larger wings, while this association was stronger in the forewings of S. aglaja whose hindwing length did not increase substantially with altitude. In high altitude habitats with low temperatures, larger wings with higher surface areas might be beneficial for thermoregulation, as the cooling rates are reduced in larger individuals (Shelly and Ludwig, 1985; Heinrich, 1993). Hindwings are not necessary for flight in lepidopterans, although maneuverability is substantially decreased when hindwings are missing or significantly reduced (Jantzen and Eisner, 2008). The varying associations of wing length with altitude in S. aglaja and F. adippe might reflect differences in ecological specialization and mobility patterns. For the highly specialized F. adippe, a high flight capacity may be required to locate potentially suitable new habitat patches, while this might not be as crucial for the relative generalist S. aglaja that can occupy a broad range of habitats. In line with this argument, a small-scaled comparative study on these two fritillary species found higher dispersal distances and lower population densities in F. adippe than in S. aglaja (Polic et al., 2020). The spatial variation in wing size was positively correlated with temperature in F. adippe, but not in the other two species. These differences might again reflect distinct ecological niche breadths and dispersal capacities, such that the specialist F. adippe might rather invest in large wings in order to locate potentially suitable new habitat patches as opposed to the relative generalist S. aglaja. The sampling strategy for F. niobe (less samples than in the other two species and a smaller latitudinal range) might have masked associations of wing size and most environmental factors. Our results suggest that the evolution of wing size in these Argynnini butterflies is influenced in part by selection for temperature regulation in combination with environmental constraints on the time (mediated by the duration of the vegetation period) available for development and growth, and that these different drivers are modulated by geographic and climatic conditions. That this has contributed to species-specific evolutionary modifications of wing size in closely related sympatric species is possibly owing to differences in ecological niche breadth and dispersal strategies as mentioned above.
Variation in dark coloration (melanisation) was mainly linked to aspects of the thermal environment across species
The degree of melanism was associated with similar environmental factors in all three fritillary species, showing that this trait is likely influenced by similar selective pressures across species. Dorsal wings were darker than ventral wings in all species, probably reflecting that the dorsal side is exposed during sun basking and melanised scales increase the solar heat absorption (Watt, 1968; Kingsolver, 1985, 1987). Melanisation mainly increased in cooler regions (e.g., with lower solar radiation in S. aglaja and lower temperature in F. niobe), which is in line with several studies on different ectothermic taxa (ants, butterflies, moths, dragonflies, beetles, and reptiles) that found darker individuals in colder environments, presumably owing to more efficient solar absorption with increasing melanism (Clusella Trullas et al., 2007; Bishop et al., 2016; Schweiger and Beierkuhnlein, 2016; Pinkert et al., 2017; Heidrich et al., 2018; Rosa and Saastamoinen, 2020; Kozlov et al., 2022). In F. adippe and F. niobe, melanisation increased with precipitation. We propose that because precipitation is often negatively correlated with temperature in European summers (Lhotka and Kyselý, 2022), exhibiting darker wings might be beneficial in areas with high precipitation. Together, these results strengthen the conclusion that the thermal milieu is a crucial selective driver for the degree of melanism in insects, and ectotherms generally. In contrast to these previous studies, however, we also found some opposing trends, i.e., melanisation increased at lower latitudes (S. aglaja) and altitudes (F. adippe), which might highlight the important anti-pathogenic properties of melanised scales. Melanin protects against pathogens and UV radiation, which would select for darker phenotypes at warmer and low latitude localities, where these pressures intensify (Mackintosh, 2001; Burtt and Ichida, 2004; Caro, 2005). We argue that the selective pressures on wing melanisation related to temperature regulation and protection from pathogens and solar radiation interact in a complex manner, and that the relative importance of these two mechanisms strongly depend on environmental factors. Dark coloration and immune defense are developmentally and genetically associated, partly because they derive from differentiation of the same cells at early embryonic stages, and darker individuals are therefore better protected against some pathogens (True, 2003; McKinnon and Pierotti, 2010). The need for temperature regulation might be a stronger driving force for melanisation than its anti-pathogenic function, as has been suggested for moths and grasshoppers (Heidrich et al., 2018; Srygley and Jaronski, 2022). This might become detrimental to insects in a changing climate: if ectotherms suppress melanisation in warmer climates to absorb less sunlight, their innate immune response can be lowered and they might be more vulnerable to, e.g., fungal attacks (Srygley and Jaronski, 2022).
Latitudinal trend in white coloration
The proportion of white scales increased in northern latitudes in S. aglaja, a pattern that parallels the geographic color morph frequency distribution in F. adippe. Fabriciana adippe f. adippe has silvery-white spots on the ventral hindwings (similar to S. aglaja) and occurs more frequently in northern Europe, while the yellow form F. adippe f. cleodoxa becomes progressively more frequent toward the south, and is dominating in the Balkans and Greece (Tolman and Lewington, 2012). The silvery-reflective patterns might represent an adaptation to forest habitats, where they may complicate the chase and capture of butterflies by visually oriented predators through the reflection of light in the prevalent complex light conditions (Olofsson et al., 2010; Kjernsmo et al., 2020). This is in line with results from a study on neotropical nymphalid butterflies, which found that species with iridescent wing patterns were more likely to occur in forest habitats (Douglas et al., 2007). The evolution of polarized light signaling in forest dwelling species may be favored due to the complex light environment and the lack of polarization in most forest backgrounds (Endler, 1993; Sweeney et al., 2003). In Europe, the forest cover increases toward higher latitudes (Schuck et al., 2002), which might render butterflies with more iridescent white scales better suited for northern European landscapes, as showcased in both S. aglaja and F. adippe.
Phenotypic integration in S. aglaja and F. adippe
The results for our analyses of trait correlations showed that fore- and hindwing length, as well as different aspects of melanisation were correlated in all studied species. This might suggest that the different dimensions of size and melanism, respectively, are inherited together in these species, e.g., via pleiotropic effects or linked genes, or influenced by similar environmental cues as a result of plasticity (Andersson, 2001; Schlichting and Wund, 2014; Fruciano et al., 2016). We detected some signs of species-specific phenotypic integration in S. aglaja and F. adippe, but essentially not in F. niobe (Figure 6). The question arises which ecological and selective pressures are involved in the evolution of these trait correlations in S. aglaja and F. adippe. Phenotypic trait integration may evolve when certain combinations of characteristics are advantageous in a specific environment, or it may result from developmental constraints or trade-offs (Forsman and Appelqvist, 1998; McKinnon and Pierotti, 2010; Dingemanse et al., 2020). In S. aglaja, but not in F. adippe or F. niobe, wing length and melanism were positively correlated, which possibly makes S. aglaja better adapted to colder habitats compared to the other two fritillary species, as both wing size and darkness can be beneficial for thermoregulation (Heinrich, 1993; Pinkert et al., 2017; Heidrich et al., 2018). Speyeria aglaja can breed at lower temperatures than the other two studied species, and its distribution range reaches further north (Warren, 1995; Ellis et al., 2010; Fox et al., 2011; Tolman and Lewington, 2012). The observed trait correlation might contribute to S. aglaja’s ecological generalism (it occupies a broad range of habitats), while the studied Fabriciana species are restricted to warm microclimates (Salz and Fartmann, 2009; Ellis et al., 2019).
That the proportion of green scales increased with wing size in S. aglaja might represent an adaptation to open grassland habitats where larger individuals might be more easily detected by predators, while the green scales might provide camouflage during perching in grassland, as in the butterfly C. rubi (Michielsen et al., 2010). An adaptation to grassland habitats might also have driven a trait correlation in F. adippe, namely that paler individuals also had more green pigments on the wings. Individuals in shady woodland habitats might have a higher need to efficiently absorb solar light, which would make darker individuals better suited for such habitats. Lighter and greener individuals might thus be better adapted in open grassy areas, which also corresponds with the findings from a study on two F. adippe color morphs showing a parallel trait value combination of the yellow spotted cleodoxa morph, which is likely better adapted to grassland habitats as opposed to the silvery spotted adippe morph (Polic, 2023). That wing size increased with the proportion of white-reflective scales in F. adippe might constitute an adaptation to forest habitats. Individuals with larger wings can accumulate solar heat faster, as shown for example for the forest dwelling butterfly Pararge aegeria (Berwaerts et al., 2001). The silvery scales in F. adippe are likely beneficial for camouflage in complex light conditions as in forest dominated areas (Olofsson et al., 2010; Kjernsmo et al., 2020), and have been suggested to function as intra-specific signals during mate choice (Wilts et al., 2013). This is in line with findings from the above-mentioned study (Polic, 2023), which revealed that the adippe morph was also larger than the cleodoxa morph, possibly resulting from adaption to distinct microhabitats in this species.
Associations of phenotypic variation with genetic structure, and some common environmental factors linked with both phenotypic and genetic variation
In all three species, forewing length was significantly associated with the genetic structure. That this association was observed across species and that it was the only one in S. aglaja, suggests a strong genetic signal (which might have masked other associations in S. aglaja) and highlights the importance of this trait, e.g., for flight capacity and temperature regulation. Similarly, the degree of melanism was associated with the genetic structure in both Fabriciana species, which might emphasize this trait’s relevance in this genus, such as for temperature regulation, UV protection, and pathogen resistance. In S. aglaja and F. adippe, there was a broad overlap of environmental variables associated with the large-scale variation in phenotypic traits and the genetic structure as assessed on the same set of individuals in a phylogeographic study (Polic et al., 2022). In S. aglaja, some variables were associated with the genetic structure, the degree of melanism (latitude, solar radiation, precipitation), and the amount of white (latitude). This might hint at adaptive modifications with respect to melanisation and white coloration along these geographic and climatic factors, although a direct association of these phenotypic traits with the genetic structure could not be verified. Alternatively, plastic responses in accordance to these environmental drivers might be involved in the trait variation distribution. In F. adippe, almost all variables that were associated with the genetic structure were also correlated with different aspects of the phenotypic appearance, viz. wing size (altitude, temperature), the degree of melanism (altitude, precipitation), and green coloration (solar radiation), suggesting that these environmental factors pose selective pressures on these phenotypic traits. These phenotypic dimensions were also associated with the genetic structure in F. adippe, hinting at adaptive genetic modifications of these traits in congruence with spatial environmental variation. Further, color morph was significantly associated with the genetic structure, corroborating results from a small-scaled study on color polymorphism in F. adippe (Polic, 2023).
Conclusions and future directions
This study revealed some parallel (e.g., similar responses to selective pressures related to the thermal environment regarding melanisation) and diverging (e.g., species-specific responses to environmental factors regarding size) large-scale patterns of variation in the three studied fritillary species. Some of the observed differences could be related to diverging ecological niche breadths and dispersal capacities, thereby highlighting the importance of comparative assessments of variation, as the relative effects of drivers of variability may vary between species.
The species-specific combinations of climatic and geographic drivers associated with the phenotypic variation in S. aglaja and F. adippe were widely the same variables that shape their large-scale genetic structure (Polic et al., 2022), which might suggest a genetic underpinning of the observed trait distribution in these species as a result of adaption to prevalent environmental conditions. This could be further tested by using a combination of rearing experiments, and phenotypic and genomic analyses. Such studies could help understand whether phenotypic trait variation stems from evolutionary genetic modifications or from plasticity responses.
Associations of wing size and environmental conditions were species-specific, posing the questions if different selective drivers shape this trait depending on the species’ biology and if responses to environmental conditions are modulated by species-specific developmental constraints. The above-mentioned approaches could further the identification of the specific mechanisms involved, and help clarify if there are trade-offs between size and other traits depending on environmental conditions. That melanism was mainly correlated with factors shaping the thermal milieu in all three species is in line with the thermal melanism hypothesis (Clusella Trullas et al., 2007; Pinkert et al., 2017). Since the need for temperature regulation seems to impose a substantial selective pressure on melanisation, the trade-off between thermoregulatory and anti-pathogenic properties of melanised scales might become disadvantageous on a warming planet. Host-pathogen interactions have changed due to global warming, in part because of shifts in species range distributions (Harvell et al., 2002; Traill et al., 2010), and in part because the replication and the transmission of viruses, bacteria, and fungi are temperature dependent (Altizer et al., 2013). The increased pathogenic pressure and the suppression of melanisation in warmer climates might thus pose a threat to ectotherms (Caro, 2005; Srygley and Jaronski, 2022). Monitoring changes in ectotherm infections over time and along spatial temperature gradients could offer important insights for making informed decisions for management and conservation plans.
Finally, our results suggest that although phenotypic integration was overall similar in the different species, S. aglaja and F. adippe exhibit some distinct trait correlations with respect to color and size of the wings, possibly driven by a combination of correlational selective pressures related to susceptibility to visually oriented predators and thermoregulatory requirements in their respective habitats. Again, rearing experiments combined with genomic assessments could help identify the underlying mechanisms involved in shaping these complex phenotypes, including an assessment of the relative roles of genetic correlations and plasticity integration.
Data availability statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.
Author contributions
DP, AF, MF, RV, and YY conceived the study. PR photographed the butterflies. PR and DP collected data on wing size and wing coloration. DP analysed the data and wrote the first draft of the manuscript, under supervision by AF. All authors contributed to the article and approved the submitted version.
Funding
Funding was provided by Linnaeus University to AF, by Birgit och Hellmuth Hertz foundation to MF, and by project PID2019-107078GB-I00/MCIN/AEI/10.13039/501100011033 to RV.
Acknowledgments
The authors thank José Antonio García-Alamá, Miguel López Munguira, Helena Romo, Camille Pitteloud, Dan Leština, Enrique García-Barros, Félix Javier González, Gheorghe Ardeleanu, Javier Pérez López, Joan Carles Hinojosa, Jordi Jubany, Juan Hernández-Roldán, Jurgen Couckuyt, Leonardo Dapporto, Michel Raymond Tarrier, Marin Goia, Paula Escuer, Raluca Vodă, Santi Viader, Sergio Montagud, Sylvain Cuvelier, and Vlad Dincă for their contribution with sample collection. The authors thank two reviewers for their comments on an earlier version of the manuscript.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2023.1087859/full#supplementary-material
References
Altizer, S., Ostfeld, R. S., Johnson, P. T., Kutz, S., and Harvell, C. D. (2013). Climate change and infectious diseases: from evidence to a predictive framework. Science 341, 514–519. doi: 10.1126/science.1239401
Andersson, L. (2001). Genetic dissection of phenotypic diversity in farm animals. Nat. Rev. Genet. 2, 130–138. doi: 10.1038/35052563
Bahrndorff, S., Ward, J., Pettigrove, V., and Hoffmann, A. A. (2006). A microcosm test of adaptation and species specific responses to polluted sediments applicable to indigenous chironomids (Diptera). Environ. Pollut. 139, 550–560. doi: 10.1016/j.envpol.2005.05.024
Berwaerts, K., Van Dyck, H., Vints, E., and Matthysen, E. (2001). Effect of manipulated wing characteristics and basking posture on thermal properties of the butterfly Pararge aegeria (L.). J. Zool. 255, 261–267. doi: 10.1017/S0952836901001327
Betts, C. R., and Wootton, R. J. (1988). Wing shape and flight behaviour in butterflies (Lepidoptera: Papilionoidea and Hesperioidea): a preliminary analysis. J. Exp. Biol. 138, 271–288. doi: 10.1242/jeb.138.1.271
Bishop, T. R., Robertson, M. P., Gibb, H., van Rensburg, B. J., Braschler, B., Chown, S. L., et al. (2016). Ant assemblages have darker and larger members in cold environments. Glob. Ecol. Biogeogr. 25, 1489–1499. doi: 10.1111/geb.12516
Brodie, E. D. (1992). Correlational selection for color pattern and antipredator behavior in the garter snake Thamnophis ordinoides. Evolution 46, 1284–1298. doi: 10.1111/j.1558-5646.1992.tb01124.x
Bruce, P., and Bruce, A. (2017). Practical statistics for data scientists: 50 essential concepts. Sebastopol, California: O'Reilly Media, Inc.
Brunton, C. F. (1998). The evolution of ultraviolet patterns in European Colias butterflies (Lepidoptera, Pieridae): a phylogeny using mitochondrial DNA. Heredity 80, 611–616. doi: 10.1046/j.1365-2540.1998.00336.x
Brunton, C. F., and Majerus, M. E. N. (1995). Ultraviolet colours in butterflies: intra-or inter-specific communication? Proc. R. Soc. B Biol. Sci. 260, 199–204. doi: 10.1098/rspb.1995.0080
Burtt, E. H. Jr., and Ichida, J. M. (2004). Gloger's rule, feather-degrading bacteria, and color variation among song sparrows. Condor 106, 681–686. doi: 10.1093/condor/106.3.681
Caro, T. (2005). The adaptive significance of coloration in mammals. Bio Sci 55, 125–136. doi: 10.1641/0006-3568(2005)055[0125:TASOCI]2.0.CO;2
Cheddadi, R., Vendramin, G. G., Litt, T., François, L., Kageyama, M., Lorentz, S., et al. (2006). Imprints of glacial refugia in the modern genetic diversity of Pinus sylvestris. Glob. Ecol. Biogeogr. 15, 271–282. doi: 10.1111/j.1466-822X.2006.00226.x
Chown, S. L., and Gaston, K. J. (2010). Body size variation in insects: a macroecological perspective. Biol. Rev. 85, 139–169. doi: 10.1111/j.1469-185X.2009.00097.x
Clench, H. K. (1966). Behavioral thermoregulation in butterflies. Ecology 47, 1021–1034. doi: 10.2307/1935649
Clusella Trullas, S., van Wyk, J. H., and Spotila, J. R. (2007). Thermal melanism in ectotherms. J. Therm. Biol. 32, 235–245. doi: 10.1016/j.jtherbio.2007.01.013
Comstock, J. H., and Needham, J. G. (1898). The wings of insects. Am. Nat. 32, 81–89. doi: 10.1086/276780
Coutsis, J. G., and Ghavalás, N. (2001). The skippers and butterflies of the Greek part of the Rodόpi massif (Lepidoptera: Hesperioidea & Papilionoidea). Phegea 29, 143–158.
Dingemanse, N. J., Barber, I., and Dochtermann, N. A. (2020). Non-consumptive effects of predation: does perceived risk strengthen the genetic integration of behaviour and morphology in stickleback? Ecol. Lett. 23, 107–118. doi: 10.1111/ele.13413
Douglas, J. M., Cronin, T. W., Chiou, T.-H., and Dominy, N. J. (2007). Light habitats and the role of polarized iridescence in the sensory ecology of neotropical nymphalid butterflies (Lepidoptera: Nymphalidae). J. Exp. Biol. 210, 788–799. doi: 10.1242/jeb.02713
Eliasson, C., Ryrholm, N., and Gärdenfors, U.. (2005). Nationalnyckeln till Sveriges flora och fauna. Fjärilar. Dagfjärilar, Hesperiidae–Nymphalidae. Artdatabanken SLU, Uppsala.
Ellis, S., Wainwright, D., Dennis, E. B., Bourn, N. A. D., Bulman, C. R., Hobson, R., et al. (2019). Are habitat changes driving the decline of the UK’s most threatened butterfly: the high Brown fritillary Argynnis adippe (Lepidoptera: Nymphalidae)? J. Insect Conserv. 23, 351–367. doi: 10.1007/s10841-019-00134-0
Ellis, S., Warren, M., Brereton, T., Toynton, P., and Chandler, D.. (2010). Fact sheet-dark green fritillary. Butterfly Conservation, Dorset. Available at: https://butterfly-conservation.org/sites/default/files/dark-green-fritillary-regional-priotrity-species-factsheet.pdf.
Endler, J. A. (1993). The color of light in forests and its implications. Ecol. Monogr. 63, 1–27. doi: 10.2307/2937121
Fang, Y., Sun, G., Bi, Y., and Zhi, H. (2015). Multiple-dimensional micro/nano structural models for hydrophobicity of butterfly wing surfaces and coupling mechanism. Sci. Bull. 60, 256–263. doi: 10.1007/s11434-014-0653-3
Ferrero, M. E., Blanco-Aguiar, J. A., Lougheed, S. C., Sánchez-Barbudo, I., De Nova, P. J. G., Villafuerte, R., et al. (2011). Phylogeography and genetic structure of the red-legged partridge (Alectoris rufa): more evidence for refugia within the Iberian glacial refugium. Mol. Ecol. 20, 2628–2642. doi: 10.1111/j.1365-294X.2011.05111.x
Fick, S. E., and Hijmans, R. J. (2017). Worldclim 2: new 1-km spatial resolution climate surfaces for global land areas. Int. J. Climatol. 37, 4302–4315. doi: 10.1002/joc.5086
Forsman, A. (2015). Rethinking phenotypic plasticity and its consequences for individuals, populations and species. Heredity 115, 276–284. doi: 10.1038/hdy.2014.92
Forsman, A. (2018). On the role of sex differences for evolution in heterogeneous and changing fitness landscapes: insights from pygmy grasshoppers. Phil Trans R Soc B 373:20170429. doi: 10.1098/rstb.2017.0429
Forsman, A., and Appelqvist, S. (1998). Visual predators impose correlational selection on prey color pattern and behavior. Behav. Ecol. 9, 409–413. doi: 10.1093/beheco/9.4.409
Forsman, A., Ringblom, K., Civantos, E., and Ahnesjö, J. (2002). Coevolution of color pattern and thermoregulatory behavior in polymorphic pygmy grasshoppers Tetrix undulata. Evolution 56, 349–360. doi: 10.1111/j.0014-3820.2002.tb01345.x
Forster, W., and Wohlfahrt, T. A. (1955). Die Schmetterlinge Mitteleuropas. 2. Tagfalter: Diurna (Rhopalocera und Hesperiidae) Stuttgart, Germany: Franckh.
Fox, R., Warren, M. S., Brereton, T. M., Roy, D. B., and Robinson, A. (2011). A new red list of British butterflies. Insect Conserv. Divers. 4, 159–172. doi: 10.1111/j.1752-4598.2010.00117.x
Franzén, M., and Johannesson, M. (2007). Predicting extinction risk of butterflies and moths (Macrolepidoptera) from distribution patterns and species characteristics. J. Insect Conserv. 11, 367–390. doi: 10.1007/s10841-006-9053-6
Fric, Z., Klimova, M., Hula, V., and Konvicka, M. (2005). Catterpillars of Argynnis aglaja (LINNAEUS, 1758) feeding on Bistorta major. Atalanta 36, 119–121.
Fruciano, C., Franchini, P., Kovacova, V., Elmer, K. R., Henning, F., and Meyer, A. (2016). Genetic linkage of distinct adaptive traits in sympatrically speciating crater lake cichlid fish. Nat. Commun. 7:12736. doi: 10.1038/ncomms12736
Hagman, M., and Forsman, A. (2003). Correlated evolution of conspicuous coloration and body size in poison frogs (Dendrobatidae). Evolution 57, 2904–2910. doi: 10.1111/j.0014-3820.2003.tb01531.x
Harvell, C. D., Mitchell, C. E., Ward, J. R., Altizer, S., Dobson, A. P., Ostfeld, R. S., et al. (2002). Climate warming and disease risks for terrestrial and marine biota. Science 296, 2158–2162. doi: 10.1126/science.1063699
Heidrich, L., Friess, N., Fiedler, K., Brändle, M., Hausmann, A., Brandl, R., et al. (2018). The dark side of Lepidoptera: colour lightness of geometrid moths decreases with increasing latitude. Glob. Ecol. Biogeogr. 27, 407–416. doi: 10.1111/geb.12703
Heinrich, B. (1993). Heat managers. Book reviews: the hot-blooded insects. Strategies and mechanisms of thermoregulation. Science 260, 1155–1156. doi: 10.1126/science.260.5111.1155
Higgins, L. G., and Riley, N. D.. (1970). A field guide to the butterflies of Britain and Europe. UK: Collins. London
Hijmans, R. J., van Etten, J., Sumner, M., Cheng, J., Baston, D., Bevan, A., et al. (2021). Raster: Geographic data analysis and modeling. R package version 3.5–11. Available at: http://CRAN.R-project.org/package=raster.
Huynen, M. A. (1996). Exploring phenotype space through neutral evolution. J. Mol. Evol. 43, 165–169. doi: 10.1007/BF02338823
James, G., Witten, D., Hastie, T., and Tibshirani, R. (2013). An introduction to statistical learning New York, USA: Springer.
Jantzen, B., and Eisner, T. (2008). Hindwings are unnecessary for flight but essential for execution of normal evasive flight in Lepidoptera. Proc. Natl. Acad. Sci. 105, 16636–16640. doi: 10.1073/pnas.0807223105
Kemp, D. J. (2008). Female mating biases for bright ultraviolet iridescence in the butterfly Eurema hecabe (Pieridae). Behav. Ecol. 19, 1–8. doi: 10.1093/beheco/arm094
Kingsolver, J. G. (1985). Thermoregulatory significance of wing melanization in Pieris butterflies (Lepidoptera: Pieridae): physics, posture, and pattern. Oecologia 66, 546–553. doi: 10.1007/BF00379348
Kingsolver, J. G. (1987). Evolution and coadaptation of thermoregulatory behavior and wing pigmentation pattern in pierid butterflies. Evolution 41, 472–490. doi: 10.2307/2409250
Kjernsmo, K., Whitney, H. M., Scott-Samuel, N. E., Hall, J. R., Knowles, H., Talas, L., et al. (2020). Iridescence as Camouflage. Curr. Biol. 30, 551–555.e3. doi: 10.1016/j.cub.2019.12.013
Komonen, A., Grapputo, A., Kaitala, V., Kotiaho, J. S., and Päivinen, J. (2004). The role of niche breadth, resource availability and range position on the life history of butterflies. Oikos 105, 41–54. doi: 10.1111/j.0030-1299.2004.12958.x
Kozlov, M. V., Oudendijk, Z., Forsman, A., Lanta, V., Barclay, M. V. L., Gusarov, V. I., et al. (2022). Climate shapes the spatio-temporal variation in colour morph diversity and composition across the distribution range of Chrysomela lapponica leaf beetle. Insect Sci. 29, 942–955. doi: 10.1111/1744-7917.12966
Kuussaari, M., Saarinen, M., Korpela, E. L., Poyry, J., and Hyvonen, T. (2014). Higher mobility of butterflies than moths connected to habitat suitability and body size in a release experiment. Ecol. Evol. 4, 3800–3811. doi: 10.1002/ece3.1187
Lande, R. (1976). Natural selection and random genetic drift in phenotypic evolution. Evolution 30, 314–334. doi: 10.2307/2407703
Lande, R. (1980). Genetic variation and phenotypic evolution during allopatric speciation. Am. Nat. 116, 463–479. doi: 10.1086/283642
Le Roy, C., Debat, V., and Llaurens, V. (2019). Adaptive evolution of butterfly wing shape: from morphology to behaviour. Biol. Rev. 94, 1261–1281. doi: 10.1111/brv.12500
Legendre, P., and Anderson, M. J. (1999). Distance-based redundancy analysis: testing multispecies responses in multifactorial ecological experiments. Ecol. Monogr. 69, 1–24. doi: 10.1890/0012-9615(1999)069[0001:DBRATM]2.0.CO;2
Lhotka, O., and Kyselý, J. (2022). Precipitation–temperature relationships over Europe in CORDEX regional climate models. Int. J. Climatol. 42, 4868–4880. doi: 10.1002/joc.7508
Losos, J. B. (2011). Convergence, adaptation, and constraint. Evolution: international journal of organic. Evolution 65, 1827–1840. doi: 10.1111/j.1558-5646.2011.01289.x
Mackintosh, J. A. (2001). The antimicrobial properties of melanocytes, melanosomes and melanin and the evolution of black skin. J. Theor. Biol. 211, 101–113. doi: 10.1006/jtbi.2001.2331
Mantel, N. (1967). The detection of disease clustering and a generalized regression approach. Cancer Res. 27, 209–220.
Matute, D. R., and Harris, A. (2013). The influence of abdominal pigmentation on desiccation and ultraviolet resistance in two species of Drosophila. Evolution 67, 2451–2460. doi: 10.1111/evo.12122
McKinnon, J. S., and Pierotti, M. E. (2010). Colour polymorphism and correlated characters: genetic mechanisms and evolution. Mol. Ecol. 19, 5101–5125. doi: 10.1111/j.1365-294X.2010.04846.x
Michielsen, K., De Raedt, H., and Stavenga, D. G. (2010). Reflectivity of the gyroid biophotonic crystals in the ventral wing scales of the green hairstreak butterfly, Callophrys rubi. J. R. Soc. Interface 7, 765–771. doi: 10.1098/rsif.2009.0352
Mikkola, K., and Rantala, M. J. (2010). Immune defence, a possible nonvisual selective factor behind the industrial melanism of moths (Lepidoptera). Biol. J. Linn. Soc. 99, 831–838. doi: 10.1111/j.1095-8312.2010.01398.x
Nieminen, M., Rita, H., and Uuvana, P. (1999). Body size and migration rate in moths. Ecography 22, 697–707. doi: 10.1111/j.1600-0587.1999.tb00519.x
Nygren, G. H., Bergström, A., and Nylin, S. (2008). Latitudinal body size clines in the butterfly Polyommatus icarus are shaped by gene-environment interactions. J. Insect Sci. 8, 1–13. doi: 10.1673/031.008.4701
Öckinger, E., Hammarstedt, O., Nilsson, S. G., and Smith, H. G. (2006). The relationship between local extinctions of grassland butterflies and increased soil nitrogen levels. Biol. Conserv. 128, 564–573. doi: 10.1016/j.biocon.2005.10.024
Öckinger, E., Schweiger, O., Crist, T. O., Debinski, D. M., Krauss, J., Kuussaari, M., et al. (2010). Life-history traits predict species responses to habitat area and isolation: a cross-continental synthesis. Ecol. Lett. 13, 969–979. doi: 10.1111/j.1461-0248.2010.01487.x
Oksanen, J., Blanchet, F. G., Friendly, M., Kindt, R., Legendre, P., McGlinn, D., et al. (2019). Vegan: community ecology package. R package version 2.5-5.
Olofsson, M., Vallin, A., Jakobsson, S., and Wiklund, C. (2010). Marginal eyespots on butterfly wings deflect bird attacks under low light intensities with UV wavelengths. PLoS One 5:e10798. doi: 10.1371/journal.pone.0010798
Piao, S., Fang, J., Zhou, L., Ciais, P., and Zhu, B. (2006). Variations in satellite-derived phenology in China's temperate vegetation. Glob. Chang. Biol. 12, 672–685. doi: 10.1111/j.1365-2486.2006.01123.x
Pinkert, S., Brandl, R., and Zeuss, D. (2017). Colour lightness of dragonfly assemblages across North America and Europe. Ecography 40, 1110–1117. doi: 10.1111/ecog.02578
Polic, D. (2023). Ecology and evolution of genetic and phenotypic variation in fritillary butterflies. Doctoral dissertation, Linnaeus University Press.
Polic, D., Tamario, C., Franzén, M., Betzholtz, P. E., Yıldırım, Y., and Forsman, A. (2020). Movements and occurrence in two closely related fritillary species. Ecol. Entomol. 46, 428–439. doi: 10.1111/een.12987
Polic, D., Yıldırım, Y., Lee, K. M., Franzén, M., Mutanen, M., Vila, R., et al. (2022). Linking large-scale genetic structure of three Argynnini butterfly species to geography and environment. Mol. Ecol. 31, 4381–4401. doi: 10.1111/mec.16594
Potyrailo, R. A., Ghiradella, H., Vertiatchikh, A., Dovidenko, K., Cournoyer, J. R., and Olson, E. (2007). Morpho butterfly wing scales demonstrate highly selective vapour response. Nat. Photonics 1, 123–128. doi: 10.1038/nphoton.2007.2
QGIS.org (2022). QGIS geographic information system QGIS Association Available at: http://www.qgis.org.
Qu, Y., Chen, C., Chen, X., Hao, Y., She, H., Wang, M., et al. (2021). The evolution of ancestral and species-specific adaptations in snowfinches at the Qinghai–Tibet plateau. Proc. Natl. Acad. Sci. 118:e2012398118:118. doi: 10.1073/pnas.2012398118
R Core Team . (2020). R: A language and environment for statistical computing. R Foundation for Statistical Computing., Vienna. Available at: https://www.R-project.org/.
Rangel, T. F., Edwards, N. R., Holden, P. B., Diniz-Filho, J. A. F., Gosling, W. D., Coelho, M. T. P., et al. (2018). Modeling the ecology and evolution of biodiversity: biogeographical cradles, museums, and graves. Science 361:eaar5452. doi: 10.1126/science.aar5452
Reinhardt, R., and Bolz, R. (2011). Rote Liste und Gesamtartenliste der Tagfalter (Rhopalocera) (Lepidoptera: Papilionoidea et Hesperioidea) Deutschlands. Natursch. Biol. Vielfalt 70, 165–194.
Rosa, E., and Saastamoinen, M. (2020). Beyond thermal melanism: association of wing melanization with fitness and flight behaviour in a butterfly. Anim. Behav. 167, 275–288. doi: 10.1016/j.anbehav.2020.07.015
Salz, A., and Fartmann, T. (2009). Coastal dunes as important strongholds for the survival of the rare Niobe fritillary (Argynnis niobe). J. Insect Conserv. 13, 643–654. doi: 10.1007/s10841-009-9214-5
Salz, A., and Fartmann, T. (2017). Larval habitat preferences of a threatened butterfly species in heavy-metal grasslands. J. Insect Conserv. 21, 129–136. doi: 10.1007/s10841-017-9961-7
Schlichting, C. D., and Wund, M. A. (2014). Phenotypic plasticity and epigenetic marking: an assessment of evidence for genetic accommodation. Evolution 68, 656–672. doi: 10.1111/evo.12348
Schuck, A., Van Brusselen, J., Päivinen, R., Häme, T., Kennedy, P., and Folving, S. (2002). Compilation of a calibrated European forest map derived from NOAA-AVHRR data. European Forest Institute. EFI Internal Report 13:44.
Schweiger, A. H., and Beierkuhnlein, C. (2016). Size dependency in colour patterns of Western Palearctic carabids. Ecography 39, 846–857. doi: 10.1111/ecog.01570
Scoble, M. J. (1992). The Lepidoptera. Form, function and diversity Oxford, UK: Oxford University Press.
Shelly, T. E., and Ludwig, D. (1985). Thermoregulatory behavior of the butterfly Calisto nubila (Satyridae) in a Puerto Rican forest. Oikos 44, 229–233. doi: 10.2307/3544693
Spitzer, L., Beneš, J., and Konvicka, M. (2009). Oviposition of the Niobe fritillary (Argynnis niobe (Linnaeus, 1758)) at sub-mountain conditions in the Czech Carpathians (Lepidoptera, Nymphalidae). Nachrichten des Entomologischen Vereins Apollo 30, 165–168.
Srygley, R. B., and Jaronski, S. T. (2022). Increasing temperature reduces cuticular melanism and immunity to fungal infection in a migratory insect. Ecol. Entomol. 47, 109–113. doi: 10.1111/een.13088
Stern, D. L. (2013). The genetic causes of convergent evolution. Nat. Rev. Genet. 14, 751–764. doi: 10.1038/nrg3483
Sullivan, J. B., and Miller, W. E. (2007). Intraspecific body size variation in macrolepidoptera as related to altitude of capture site and seasonal generation. J. Lepidopterist's Soc. 61, 72–77.
Sweeney, A., Jiggins, C., and Johnsen, S. (2003). Polarized light as a butterfly mating signal. Nature 423, 31–32. doi: 10.1038/423031a
The GIMP Development Team . (2019). GIMP. Available at: https://www.gimp.org.
Tinnert, J., and Forsman, A. (2017). The role of dispersal for genetic and pehnotypic variation: insights from comparisons of sympatric pygmy grasshoppers. Biol. J. Linn. Soc. 122, 84–97. doi: 10.1093/biolinnean/blx055
Tolman, T., and Lewington, R.. (2012). Schmetterlinge Europas und Nordwestafrikas. 2. Franckh-Kosmos Verlag, Stuttgart.
Traill, L. W., Lim, M. L. M., Sodhi, N. S., and Bradshaw, C. J. A. (2010). Mechanisms driving change: altered species interactions and ecosystem function through global warming. J. Anim. Ecol. 79, 937–947. doi: 10.1111/j.1365-2656.2010.01695.x
True, J. R. (2003). Insect melanism: the molecules matter. Trends Ecol. Evol. 18, 640–647. doi: 10.1016/j.tree.2003.09.006
Tsai, C.-C., Childers, R. A., Nan Shi, N., Ren, C., Pelaez, J. N., Bernard, G. D., et al. (2020). Physical and behavioral adaptations to prevent overheating of the living wings of butterflies. Nat. Commun. 11:551. doi: 10.1038/s41467-020-14408-8
van Swaay, C., and Warren, M.. (1999). Red data book of European butterflies (Rhopalocera). Nature and Environment, No. 99, Council of Europe Publishing, Strasbourg.
van Swaay, C., Wynhoff, I., Verovnik, R., Wiemers, M., López Munguira, M., Maes, D., et al. (2010). Argynnis adippe. The IUCN red list of threatened species 2010. e.T174203A7028071. Downloaded on 06 February 2020.
Vargas, C. A., Lagos, N. A., Lardies, M. A., Duarte, C., Manríquez, P. H., Aguilera, V. M., et al. (2017). Species-specific responses to ocean acidification should account for local adaptation and adaptive plasticity. Nat. Ecol. Evol. 1:0084. doi: 10.1038/s41559-017-0084
Violle, C., Reich, P. B., Pacala, S. W., Enquist, B. J., and Kattge, J. (2014). The emergence and promise of functional biogeography. Proc. Natl. Acad. Sci. 111, 13690–13696. doi: 10.1073/pnas.1415442111
Voelkl, B., Altman, N. S., Forsman, A., Forstmeier, W., Gurevitch, J., Jaric, I., et al. (2020). Reproducibility of animal research in light of biological variation. Nat. Rev. Neurosci. 21, 384–393. doi: 10.1038/s41583-020-0313-3
Warren, M. S. (1995). “Managing local microclimates for the high brown fritillary, Argynnis adippe” in Ecology and conservation of butterflies. ed. A. S. Pullin (Dordrecht: Springer), 198–210.
Watt, W. B. (1968). Adaptive significance of pigment polymorphisms in Colias butterflies. I. Variation of melanin pigment in relation to thermoregulation. Evolution 22, 437–458. doi: 10.1111/j.1558-5646.1968.tb03985.x
Wiens, J. J. (2015). Faster diversification on land than sea helps explain global biodiversity patterns among habitats and animal phyla. Ecol. Lett. 18, 1234–1241. doi: 10.1111/ele.12503
Wilson, K., Cotter, S. C., Reeson, A. F., and Pell, J. K. (2001). Melanism and disease resistance in insects. Ecol. Lett. 4, 637–649. doi: 10.1046/j.1461-0248.2001.00279.x
Wilts, B. D., Pirih, P., Arikawa, K., and Stavenga, D. G. (2013). Shiny wing scales cause spec(tac)ular camouflage of the angled sunbeam butterfly, Curetis acuta. Biol. J. Linn. Soc. 109, 279–289. doi: 10.1111/bij.12070
Zaman, K., Hubert, M. K., and Schoville, S. D. (2019). Testing the role of ecological selection on colour pattern variation in the butterfly Parnassius clodius. Mol. Ecol. 28, 5086–5102. doi: 10.1111/mec.15279
Zeuss, D., Brunzel, S., and Brandl, R. (2017). Environmental drivers of voltinism and body size in insect assemblages across Europe. Glob. Ecol. Biogeogr. 26, 154–165. doi: 10.1111/geb.12525
Zhang, C., Settele, J., Sun, W., Wiemers, M., Zhang, Y., and Schweiger, O. (2019). Resource availability drives trait composition of butterfly assemblages. Oecologia 190, 913–926. doi: 10.1007/s00442-019-04454-5
Zheng, Y., Gao, X., and Jiang, L. (2007). Directional adhesion of superhydrophobic butterfly wings. Soft Matter 3, 178–182. doi: 10.1039/B612667G
Keywords: evolution, geography, melanism, phenotypic integration, thermal environment, wing coloration, wing size
Citation: Polic D, Yıldırım Y, Vila R, Ribeiro Cardoso PR, Franzén M and Forsman A (2023) Large-scale spatial variation and phenotypic integration in three Argynnini species inform about functions and evolutionary drivers of butterfly wings. Front. Ecol. Evol. 11:1087859. doi: 10.3389/fevo.2023.1087859
Edited by:
Daniel de Paiva Silva, Goiano Federal Institute (IFGOIANO), BrazilReviewed by:
Simeão Moraes, State University of Campinas, BrazilFlorian Menzel, Johannes Gutenberg University Mainz, Germany
Copyright © 2023 Polic, Yıldırım, Vila, Ribeiro Cardoso, Franzén and Forsman. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Daniela Polic, ZGFuaWVsYS5wb2xpY0BsbnUuc2U=
†Present address: Yeşerin Yıldırım, Department of Ecology and Genetics, Uppsala University, Uppsala, Sweden