Climate shapes the spatiotemporal variation in color morph diversity and composition across the distribution range of Chrysomela lapponica leaf beetle

Color polymorphism offers rich opportunities for studying the eco‐evolutionary mechanisms that drive the adaptations of local populations to heterogeneous and changing environments. We explored the color morph diversity and composition in a Chrysomela lapponica leaf beetle across its entire distribution range to test the hypothesis that environmental and climatic variables shape spatiotemporal variation in the phenotypic structure of a polymorphic species. We obtained information on 13 617 specimens of this beetle from museums, private collections, and websites. These specimens (collected from 1830–2020) originated from 959 localities spanning 33° latitude, 178° longitude, and 4200 m altitude. We classified the beetles into five color morphs and searched for environmental factors that could explain the variation in the level of polymorphism (quantified by the Shannon diversity index) and in the relative frequencies of individual color morphs. The highest level of polymorphism was found at high latitudes and altitudes. The color morphs differed in their climatic requirements; composition of colour morphs was independent of the geographic distance that separated populations but changed with collection year, longitude, mean July temperature and between‐year temperature fluctuations. The proportion of melanic beetles, in line with the thermal melanism hypothesis, increased with increasing latitude and altitude and decreased with increasing climate seasonality. Melanic morph frequencies also declined during the past century, but only at high latitudes and altitudes where recent climate warming was especially strong. The observed patterns suggest that color polymorphism is especially advantageous for populations inhabiting unpredictable environments, presumably due to the different climatic requirements of coexisting color morphs.


Introduction
More than a century ago, Gerould (1911, p. 257) wrote: "There is perhaps no phenomenon of greater general interest to students of organic evolution than polymorphism." This statement remains valid today. Why are some species or geographic populations polymorphic, whereas others are not? Which factors influence the contributions of individual morphs to the composition of polymorphic populations, and how does this composition affect the fates of these populations both ecologically and over evolutionary time scales? Variation in the diversity and composition of color morphs among populations and over time is likely influenced by the combined effects of natural selection, stochastic genetic processes, and the consequences that interindividual variation may have for the ecological success of populations. The possible benefits that may accrue to polymorphic populations and species first captivated evolutionary biologists more than half a century ago (e.g., Fisher, 1930;Dobzhansky, 1951). Thereafter, interest in these questions temporarily declined, but recently scientific attention to the consequences of intraspecific trait variation has shown a resurgence (reviewed in Forsman, 2016;Forsman & Wennersten, 2016).
The idea of adaptive polymorphism was reinvigorated by Forsman et al. (2008), who developed a verbal model and predicted that populations consisting of two or more alternative color morphs that occupy different niches will utilize a greater diversity of resources, enjoy enhanced stability and persistence, have a higher capacity for establishment and range expansion, and be less vulnerable to environmental change and local extinctions compared to less variable populations. These predictions emphasize the importance of studies of geographical variation in polymorphism, which may help not only reveal the contributions of fundamental eco-evolutionary processes to the population divergence that can potentially lead to speciation (Forsman et al., 2008;McLean & Stuart-Fox, 2014), but also provide insights that benefit the planning of conservation measures (Forsman, 2016). The latter aspect is especially critical nowadays when habitat degradation, fragmentation, and climate change are causing global biodiversity redistributions and losses (Butchart et al., 2010;Pecl et al., 2017), including rapid declines in insect populations (Wagner, 2020).
Studying spatiotemporal variation in color morph diversity and composition, and its associations with environmental variables, is interesting for two reasons as this has potential to uncover both the causes and the consequences of polymorphisms. Evidence is mounting that selection imposed by spatially variable and temporally changing environmental conditions plays an important role in shaping the phenotypic and genetic variation seen among and within natural populations (e.g., Halsson & Björklund, 2012;Yildirim et al., 2018). At the same time, results from experimental studies and phylogenybased comparative analyses largely support the idea that populations and species with greater variability are better at coping with heterogeneous, novel and changing environments (Forsman, 2016;Forsman & Wennersten, 2016). However, attempts to evaluate the predictions that color polymorphisms likely are spatiotemporally dynamic, changing in relative morph frequencies, shifting from a polymorphic to a monomorphic state, or gradually developing into geographic variation if divergent selection favors alternative morphs in unlike environments, are rare (Forsman et al., 2008). Testing these predictions requires observational data that ideally cover the entire distributional range of a species.
Studies based on museum collections that represent different environments and time periods have a unique potential to deepen our understanding of how eco-evolutionary processes shape biodiversity and how polymorphic conditions might influence the ecological success and resilience of populations and species. However, large-scale and long-term spatiotemporal variation in morph frequencies is adequately documented for only a few polymorphic species (McLean & Stuart-Fox, 2014;Forsman, 2016), thus hampering generalizations across taxa and environments. In this paper, we explore spatiotemporal variations in color polymorphism of a Palaearctic leaf beetle Chrysomela lapponica Linnaeus, 1758 ( Fig. 1) to uncover the causes and consequences of phenotypic and genetic variability between individuals. The accumulated knowledge about C. lapponica, gained by previous investigations (Mikhailov, 2001;Zvereva et al., 2002Zvereva et al., , 2019Fatouros et al., 2006;Machkour-M'Rabet et al., 2008;Mardulyn et al., 2011;Geiselhardt Zverev et al., 2017Zverev et al., , 2018Doktorovová et al., 2019), identifies this species as a prospective model for uncovering mechanisms that shape the distribution and dynamics of color polymorphism.
Climate has multiple direct and indirect impacts on living beings (Bale et al., 2002;Buckley et al., 2012) and is frequently hypothesized to influence the spatial variations observed in the occurrence and frequency of color morphs in ectotherms (Broennimann et al., 2014;McLean & Stuart-Fox, 2014). Based on the theory of adaptive polymorphism (Forsman et al., 2008), we predict that environments with more seasonal and fluctuating temperatures should favor and select for phenotypically more diverse (color polymorphic) populations compared with environments with stable climates. This is because polymorphic populations cope better (compared with monomorphic populations) with temporally variable environments, such as in areas with greater seasonal changes where temperatures fluctuate strongly. We also hypothesize that color morphs of C. lapponica differ in their climatic requirements, and that these differences allow this species to occupy diverse environments.
Perhaps the best studied example of differential responses of color morphs to climate is variation in the level of melanization. The thermal melanism hypothesis states that, at low temperatures, dark individuals have an advantage over light individuals because they warm up more quickly at any given level of solar radiation (reviewed by Clusella-Trullas et al., 2007). Based on this hypothesis, we predict that the proportion of dark morphs in C. lapponica populations will increase in cold climates, with increases in both latitude and altitude. We also expect that the frequency of dark beetles will decrease with collection year due to climate warming. The photoprotection hypothesis (Law et al., 2020) yields similar prediction, that is, an increase in melanization in alpine regions, due to the protection provided by melanin against UV-B radiation (Majerus, 1998).
The stochastic processes (e.g., gene flow and genetic drift) and historical events (e.g., founder effects) can also affect geographic variation in polymorphic species (Reillo & Wise, 1988;Yildirim et al., 2018). These processes may contribute to increased differentiation in color morph composition among populations and to decreased neutral (but not necessarily functional) variability within polymorphic populations due to the organisms' limited dispersal, even if the environment is completely homogeneous (Soininen et al., 2007). We therefore examined whether the signatures of these processes can be found in the current distribution of the color polymorphism in C. lapponica populations.
The ultimate goal of the present study was to test the hypothesis that climate shapes the geographic variation in the phenotypic structure of a polymorphic species in such a way that the highest diversity is observed in cold and variable environments. To achieve this goal, we explored the spatial and temporal variation in the withinpopulation polymorphism level and in the frequencies of different color morphs of the leaf beetle C. lapponica across its entire distribution range, and we identified the environmental conditions explaining this variation. We answered the following questions: (1) How is the spatial variation in the phenotypic diversity and in the relative frequencies of color morphs of C. lapponica associated with geographic location of the population (latitude, longitude and altitude) and with climate variables (mean values of temperature and precipitation, their seasonal cycle amplitudes and between-year fluctuations)? (2) Do dissimilarities in color morph composition increase with increasing geographic distance between populations?
(3) How do phenotypic diversity and color morph frequencies vary over time (decades to centuries), and are these temporal changes consistent with recent climate change? We answered these questions by analyzing information on 13 617 specimens of C. lapponica collected from 1830 to 2020 across the continent-wide distribution range of this species, spanning 33°latitude, 178°longitude, and 4200 m altitude.

Study species
The leaf beetle C. lapponica is 5-9 mm in length, and it has variable, generally aposematic, color patterns ( Fig. 1) signaling its chemical defense. The analyses of individual progenies of C. lapponica females from crossing experiments conducted in different populations (Zverev et al., 2018;Kozlov, Mikhailov, Zverev, & Zvereva, unpublished) suggest that elytra coloration in C. lapponica, as in its close relatives (Brown, 1956), is at least partly heritable. Almost any combinations of color morphs occur in natural populations (Data S1 and S2), and hence taxonomists have never assigned subspecific rank to C. lapponica populations that differ in color morph frequencies, although some of these populations demonstrate high genetic differentiation (Machkour-M'Rabet et al., 2008;Mardulyn et al., 2011) and even reproductive isolation (Fatouros et al., 2006).
Chrysomela lapponica is univoltine across its entire range and feeds on both Salicaceae and Betulaceae, depending on the population. Adults emerge from the soil after hibernation, copulate in spring soon after leaf flush, and oviposit on leaves of their host plants. Larvae grow throughout summer and pupate on leaves. New generation adults emerge in the second half of the summer and feed for several weeks before burying in soil for overwintering (Mikhailov, 2001;Zvereva et al., 2002;Fatouros et al., 2006).

Assessment of color morph diversity
Photographs or specimens of C. lapponica and their label data were obtained from multiple museums, personal collections, publications and web pages (Data S1). We excluded the specimens that had no collection locality information. Coordinates of the localities and collection year, when missing from the labels, were obtained from Google Earth, collector's diaries and published expedition reports. Accuracy of the coordinates was estimated by the diameter of a circle that likely includes the exact collection locality (Data S1). The data with the lowest accuracy (error range >160 km; 293 individuals) were included only in the summary statistics of morph frequencies.
Each individual was attributed to one of five color morphs ( Fig. 1) based on percentage of orange, black and metallic color coverage on elytra (Table S1). This classi-fication was performed by Z.O. and controlled by E.L.Z. and M.V.K. The beetles from the personal collection of O.Y.K. were classified by the collection owner.
For our analyses, we aggregated the individuals by locality (i.e., the site with unique coordinates) and by collection year. In 58% of localities, all the individuals were collected during a single year. In six localities, where the difference between the last and the first collection year exceeded 20 years, the individuals were aggregated by samples in such a way that collection period did not exceed 20 years. For the aggregated samples, the collection year was averaged across individuals.
Substantial fraction of our samples consisted of a single individual (Fig. S1) and therefore had zero diversity. The proportion of samples which include only one morph of Chrysomela lapponica (diversity of which is zero) decreased rapidly ( Fig. S2) when we excluded small samples from the analysis. Based on the latter pattern, we selected the minimum sample size of ten beetles for the analyses of phenotypic diversity. This sample size gives a reasonable balance between the number of available samples (163) and the accuracy of the obtained diversity estimates. These 163 samples (large samples hereafter; Data S2) jointly included 11 577 beetles, that is, 85% of all studied individuals (Data S1).
For each large sample, we calculated (1) the proportion of each color morph; (2) the proportion of melanic (dark patterned and black) beetles among nonmetallic (orange, light patterned, dark patterned and black) beetles; and (3) the Shannon diversity index based on natural logarithms of proportions of the five color morphs. The metallic beetles were excluded from the analysis (2) because their coloration differs greatly from coloration of the remaining four morphs, which show gradual variation in melanization.

Climate data
The mean temperatures of January and July, mean annual temperature (MAT), and mean annual precipitation (MAP) from 1981-2019, extracted from NASAPOWER archive (power.larc.nasa.gov) with R package "nasapower" (Sparks, 2018), were used to calculate for each locality: (1) the mean values of all these four characteristics of climate; (2) the average difference in temperature and precipitation between July and January, as a measure of seasonality (seasonal cycle amplitude); and (3) the temperature and precipitation differences between maximum and minimum values observed in July between 1981 and 2019, as a measure of annual climatic variability.

Data analysis
We separately conducted two forward stepwise regression analyses to identify the models that best explain the variation in level of color polymorphism (quantified using Shannon diversity index) and in the proportion of melanic individuals. These analyses included both linear and quadratic forms of the ten explanatory variables (latitude, longitude, altitude, collection year, mean temperature of July, MAP, MAT, seasonality of temperature, and annual variability in temperature and precipitation). We excluded mean temperature of January (which was used for calculation of seasonality) from these analyses to avoid a singularity problem.
To aid interpretation of the obtained results, we regressed the level of color polymorphism and the proportion of melanic beetles (in separate analyses) to each of ten explanatory variables one by one. In these analyses, we disregarded geographic factors, because Mantel test (see below) did not reveal spatial autocorrelation. We excluded mountain populations (altitude >800 m) from the analyses of latitudinal and longitudinal patterns, and excluded European populations from the analysis of altitudinal pattern, because only 4 of 88 European samples were collected above 800 m. These restrictions did not affect either the direction or the statistical significance of the detected patterns, but increased the proportion of variation explained by the models by a factor of 1.5 to 2 (see caption of Fig. 3). The analysis of temporal changes in proportions of melanic beetles was limited to samples that were collected after 1850 and that included at least one melanic specimen. We selected between models based on the lowest Akaike information criterion (AIC) value. These analyses were performed in R (R Core Team, 2020) and SAS (SAS Institute, 2009).
We explored how variation in the frequencies of color morphs was associated with geographic and environmental factors using a distance-based redundancy analysis (dbRDA). The case scores were first obtained from Principal Coordinates Analysis performed on a distance matrix, and these were further constrained by environmental predictors using RDA with forward selection procedure (Legendre & Anderson, 1999). The statistical significance of the predictors was assessed using Monte-Carlo test (499 permutations). In addition to this multivariate analysis, the frequencies of color morphs were compared between Europe and Asia and between lowlands and mountains by chi-square tests. We also evaluated the associations of color morph frequencies with geographical variables (latitude, longitude and altitude) using mixed model ANCOVA (SAS GLIMMIX procedure, type 3 tests) with multinomial distribution and cumlogit link function (SAS Institute, 2009), and by calculating Spearman's rank correlation coefficients for frequencies of individual color morphs followed by sequential Bonferroni's correction.
We used the Mantel statistic calculated by the R package "vegan" (Oksanen et al., 2019) based on Spearman's rank correlation coefficient (999 permutations) to evaluate whether the similarity in color morph frequencies between populations depends on the geographical distance that separates populations. We calculated pairwise similarities (SI) in morph frequencies among 163 large samples as follows: , and p m are proportions of orange, light patterned, dark patterned, black and metallic beetles in the compared populations (i and j). This index varies from 0, when the compared samples do not share any morph, to 1, when morph frequencies in the compared samples are identical (Zhivotovsky, 1979). Geographical distances between localities were calculated by the "rgeos" and "rgdal" R packages (Bivand, 2020). The climatic requirements (i.e., climatic niches, in terms of mean temperatures of January and July, MAT and MAP) of individual morphs were estimated by averaging climatic characteristics across localities, from which these morphs were recorded (Fig. S3), and were compared among morphs by ANOVA followed with Tukey's HSD test.

Data overview and overall distribution range of Chrysomela lapponica
We obtained data on 13 617 individuals, which were collected from 1830-2020 (median year 1998) in 31 countries (Table S2) across the Palaearctic (Fig. 2): from Spain (8°W) throughout the central and northern Eurasia to Chukotka (170°E). The latitudinal range of C. lapponica in Europe extends from Spain (42°N) to Norway (71°N), and in Asia it extends from North Korea (38°N) to Yakutia in Russia (68°N). This species was collected from the sea level in many localities to 4200 m above sea level in the mountains of Kyrgyzstan. It was found in a variety of habitats, from Arctic shrubby tundra to broadleaved forests at the latitudinal scale and from lowlands to alpine habitats at the vertical scale. Mean annual temperature in the collecting localities ranged from +12.5°C (Pontevedra, Spain) to −13.9°C (Verkhoyansk, Russia), and total annual precipitation ranged from About one-third of specimens (4716) originated from 38 museums, 3 private collections, and 13 web projects; the remaining specimens were collected by the coauthors for morphological, biogeographical, ecological, and/or environmental research (Data S1). The vast majority of individuals belonged to light patterned (55.9%) and metallic (27.9%) morphs (Table S1).
All beetles were grouped into 963 samples, 163 of which contained at least 10 individuals. The 889 samples from museum collections included 1-127 individuals (median value 2), and 74 ecological samples included 1-1494 individuals (median value 50).

Variation in the level of the color polymorphism
The forward stepwise regression analysis, in which the 10 explanatory variables (listed in Table S3) were evaluated together, showed that variation in the level of color polymorphism was significantly associated with altitude, longitude, July temperature, MAP, betweenyear temperature fluctuations and collection year (Tables 1, S4). When the associations of polymorphism with the explanatory variables were analyzed one by one (Fig. 3, Table S3), the color morph diversity increased with increasing latitude (Fig. 3A), altitude (Fig. 3C), and between-year fluctuation in July temperatures (Fig. 3H), decreased with increasing July temperature (Fig. 3E) and showed hump-shaped changes with longitude (Fig. 3B).

Variation in color morph composition
The distribution ranges of color morphs overlap to a great extent (Fig. S3). Nevertheless, climatic requirements significantly differed among color morphs in terms of all four climatic variables (Fig. 4).
At the continental scale, the proportion of orange and light patterned beetles in both lowland and mountain habitats in Asia was higher than in Europe ( Fig. 5; χ 2 = 1996.8, df = 1, P < 0.0001), whereas dark patterned and black beetles mostly originated from northern Europe and central Asia (Fig. S3). The metallic beetles were most common in central Eurasia. Across both continents, the proportion of melanic beetles was higher in mountains (>800 m above sea level) than at lower latitudes ( Fig. 5; χ 2 = 1460.1, df = 1, P < 0.0001).
The multinomial model that included frequencies of all five color morphs did not converge. However, when morphs were combined into three groups: nonmelanic (orange and light patterned), melanic (dark patterned and black), and metallic, the model revealed that phenotypic composition of our samples changed with latitude (F 2, 500.1 = 21.42, P < 0.0001), longitude (F 2, 590.6 = 4.44, P = 0.01), and altitude (F 2, 542.7 = 12.10, Table 1 Characteristics of the best-fit models explaining variation in the Shannon diversity index (SDI) and in the proportion of melanic beetles (PMB): the outcomes of the forward step-wise regression analysis. MAP, mean annual precipitation; Year, collection year; Seasonality, the difference between mean temperatures of July and January.
P < 0.0001). The correlation analyses, in which the different morphs were analyzed one by one, demonstrated that the proportion of dark patterned beetles increased with increases in latitude and altitude; the proportion of metallic morphs decreased with increasing latitude; and the proportion of light patterned beetles increased from West to East (Table S5).
The dbRDA analysis showed that the most influential drivers of color morph composition were mean July temperature, collection year, longitude, and altitude, which together accounted for 87.4% of the total variation explained by the model (Table S6). The variation in frequencies of black, orange, dark, and light patterned beetles was primarily associated with collection year and variation in July temperature, while variation in occurrence of metallic beetles was mostly associated with July temperature (Fig. 6). The similarity in color morph frequencies among 163 large samples was independent of the geographic distance separating respective localities (Mantel test: r S = 0.02, P = 0.31).

Variation in the proportion of melanic beetles
The forward stepwise regression analysis, in which the 10 predefined explanatory variables were analyzed together, suggests that the proportion of dark morphs decreased with the increase in seasonality (Table S4). For the relationships between the proportion of dark morphs and all explanatory variables in separate analyses consult Table S7.
The temporal changes in the proportions of melanic beetles depended on summer temperatures (interaction term: F 1, 97 = 6.92, P = 0.01). This proportion declined over time in populations inhabiting localities with colder climates (mean July temperature below 14°C) at high latitudes and altitudes (Fig. 7A), whereas no temporal shifts were observed in localities with warmer climates (Fig. 7B).

Geographic variation in within-population polymorphism
We documented high variation among populations in terms of the expression of color polymorphism, which ranged from monomorphic populations (metallic: Belarus and the middle Urals; light patterned: Central Europe) to almost equal representation of all five morphs (the northern Urals). Spatial variation in the expression of color polymorphism, both systematic (clinal) and random, has previously been reported in other insect species (Sánchez-Guillén et al., 2011;Noriyuki & Osawa, 2015;Strickland et al., 2019). However, the majority of these studies considered only a small number of populations and covered only part of the distribution range of the study species. In this respect, our study is unique in its geographical coverage, because it includes samples from 959 populations across the entire continent-wide distribution range of C. lapponica (Text S1), thereby allowing to analyze polymorphism pattern of a species. The stepwise regression and dbRDA analyses consistently indicate that the among-population mosaic of color polymorphism in this leaf beetle is primarily shaped by climatic factors.
The greater phenotypic and genetic variation between individuals has been hypothesized to make polymorphic populations less vulnerable to environmental changes compared to monomorphic populations (Forsman et al., 2008;Forsman & Wennersten, 2016). In agreement with Fig. 3 The relationships of color polymorphism diversity (Shannon diversity index) of Chrysomela lapponica populations with geographic position of the locality (A, latitude; B, longitude; C, altitude), collection year (D) and climate (E, mean temperature of July; F, mean annual temperature; G, mean annual precipitation; and H, between-year fluctuations in July temperature). The analyses of the diversity relationships with latitude and longitude are based on 123 samples from localities situated <800 m above sea level; analysis of the effect of altitude is based on 75 samples from Asia; other analyses involved all 163 samples. The analyses of all 163 samples demonstrated apparently the same relationships of diversity with latitude (R 2 = 0.10, P < 0.001), longitude (R 2 = 0.13, P < 0.001) and altitude (R 2 = 0.03, P = 0.02); nevertheless, the removal of some samples from these analyses is justified by the unbalanced distribution of our samples between continents and altitudes, as explained in Materials and Methods. this hypothesis, we found an increase in the diversity of color morphs with distinctive climatic requirements (Fig. 4) in response to both the summer climate (i.e., with a decrease in June temperatures) and to the between-year variability of these temperatures (Table 1). This suggests that polymorphic populations of C. lapponica have a better ability to establish in and to cope with habitats that have temporally variable weather conditions, where tem- peratures fluctuate in an unpredictable manner. This is because each morph has advantages under different conditions, thereby allowing the population as a whole to persist. Oscillating selection may also contribute to the maintenance of genetic and phenotypic polymorphism in seasonal environments with fluctuating temperatures (Bertram & Masel, 2019).
In addition to climatic drivers, the results of our bestfit model suggest that a significant part of the variation in polymorphism level is explained by the longitude and altitude of the collecting localities. In particular, the phenotypic diversity of C. lapponica showed a hump-shaped longitudinal pattern and attained its maximum in western and central Asia (Fig. 3B), in approximately the same region where the larvae of this species exhibit the ancestral state of their defensive secretion system (Geiselhard et al., 2015;Zverev et al., 2017). This cooccurrence can be seen as evidence of the origin of this species in western Asia; therefore, the irregularities observed in the largescale distribution of color morphs may have resulted from climatic changes during the Pleistocene, in combination with random genetic drift. However, we cannot exclude the possibility that other as yet unknown factors that systematically change with longitude and altitude have also substantially contributed to the observed spatial variation of color polymorphism in C. lapponica.
Despite the systematic geographical variation, the similarity among C. lapponica populations in color morph frequencies appeared to be independent of the distance between the localities. This finding conforms to the results from analyses of the genetic similarity of nine European populations of C. lapponica (Machkour-M'Rabet et al., 2008) and points to the conclusion that stochastic and historical factors are paramount in shaping the current diversity of C. lapponica populations. These historical factors, which include range fragmentation caused by glaciations and sequential founder events during postglacial expansions, are no longer operating, but their genetic signatures may still be present in contemporary populations (Eckert et al., 2008). An alternative explanation for the lack of an association between geographic distance and color morphs composition in a population is that important environmental conditions and selective pressures vary in a nonsystematic manner Fig. 6 dbRDA ordination based on similarities among large samples of Chrysomela lapponica, with five leaf beetle morphs (black font) projected in reduced ordination space relative to significant explanatory variables (red font). T_July_mean, mean July temperature; T_July_diff, between-year fluctuations in July temperature. over small spatial scales, so that local adaptations of populations appear uncorrelated with the geographic distances that separate populations.

Geographic variation in color morph frequencies
Geographic variation in morph frequencies between populations may be driven by variations in biotic and abiotic factors (McLean et al., 2014). The increase in frequency of dark individuals in colder climates (at both high latitudes and altitudes, where ectothermic animals need to simultaneously develop adaptations to low temperatures and to high levels of solar radiation: Körner, 2007;De Frenne et al., 2013) is in line with several other studies and in agreement with the thermal melanism hypothesis (Clusella-Trullas et al., 2007;and references therein). At the same time, the consistency between Fig. 7 The relationship of the proportion of melanic beetles in Chrysomela lapponica populations with the collection year in cold (A, mean July temperature below 14°C; n = 47 samples) and warm (B, mean July temperature 14°C and above; n = 54 samples), based on the analysis of samples that include at least one melanic individual. changes in C. lapponica melanization with increases in both elevation and latitude suggests that these changes, like in geometrid moths (Heindrich et al., 2018), are unlikely explained by the photoprotection hypothesis, because UV decreases toward the poles but increases with elevation. More likely, latitudinal and altitudinal trends in frequencies of melanic individuals are both driven by temperature, which decreases with both latitude and altitude. However, protective role of melanin-based dark coloration against ultraviolet radiation can partly contribute to the observed altitudinal pattern.
Low temperatures may increase the frequencies of dark morphs in C. lapponica populations due to the mating advantages of dark patterned males over light patterned males and due to developmental plasticity (Zverev et al., 2018). This temperature effect is so strong that it maintains high frequencies of melanic beetles in cold climates, despite their lower fecundity and higher sensitivity to variations in host plant quality (Zvereva et al., 2002;Zverev et al., 2018).
The positive correlation across populations between the frequency of dark morphs of C. lapponica and precipitation is similar to the pattern observed in some other polymorphic species (Harris et al., 2012;Broenniman et al., 2014). This finding adds to the limited evidence for an increase in insect melanization due to higher eumelanin deposition in habitats that are more humid-an observation that is tentatively associated with better protection against more abundant parasites and pathogens (Delhey, 2019). The latter could be achieved in several ways, in particular because pleiotropic effects lead to better anti-inflammatory, antipyretic, and antioxidative responses in darker animals (Ducrest et al., 2008).
The high frequency of the metallic morph at lower latitudes and elevations is opposite to the distribution observed for the dark morph, and dbRDA ordination confirmed that the metallic and non-metallic beetles are associated with different explanatory variables. We suggest that the iridescence of the elytra (i.e., the high surface reflectance) protects metallic beetles from overheating (Vulinec, 1997;Mikhailov, 2008), thereby providing them with advantages in warmer climates. In addition, spatial variation in morph frequencies can also be shaped by selective pressure from predators (McLean et al., 2014;Matthews et al., 2018), which can vary geographically. The metallic morph of C. lapponica differs considerably from other morphs in terms of the strength and memorability of its aposematic signal, thereby resulting in better survival of metallic morphs exposed to bird predation (Doktorovová et al., 2019). This survival benefit of anti-predatory defenses (including warning coloration) is greater at higher predator density. Bird predation often decreases from low to high latitudes (Matthews et al., 2018;Zvereva et al., 2020); therefore, the increase in frequency of metallic beetles toward the south may also be related to the higher pressure from bird predation at lower latitudes.

Temporal changes in polymorphism and morph frequencies
A decline in the proportion of dark beetles during the past decades has been reported for populations of C. lapponica inhabiting the Kola Peninsula in northern European Russia (Zvereva et al., 2019). Our current study demonstrates that this pattern is typical for all C. lapponica populations inhabiting cold regions, in which the proportion of dark beetles decreased during the past 100 years from 40% to 20% (Fig. 7A), reaching the level observed in warm regions before the 1950s (Fig. 7B). Similar declines in the proportion of melanic morphs have been observed in some other insect species (Brakefield & de Jong, 2011;Clusella-Trullas & Nielsen, 2020;Scriber, 2020). The reasons for this decline are not always clear (Scriber, 2020), but the decline in C. lapponica likely occurs because dark-patterned beetles lose their mating advantages due to climate warming (Zverev et al., 2018;Zvereva et al., 2019). This suggestion is supported by the discovery of this declining trend only in regions with cold climates (i.e., at high latitudes and altitudes), where recent increases in mean temperatures have been especially strong (Walther et al., 2002;Zvereva et al., 2016). A similar pattern was observed in the butterfly Colias meadii, which showed a decrease in wing melanization from 1953 to 2013 in northern regions but not in southern regions (MacLean et al., 2016). This finding supports our suggestion that the rates of temporal changes in leaf beetle coloration are associated with the rates of climate change.
A temporal decrease in the level of color polymorphism was detected when a suite of explanatory variables was analyzed simultaneously (Table 1). The decrease in color polymorphism with climate warming is in line with pattern observed in geographic gradients, that is, decrease in color polymorphism from high to low latitudes and altitudes. This decrease in phenotypic diversity with climate warming can make populations more vulnerable to future climates, because climate change brings about more extreme weather conditions (Seneviratne et al., 2012), including greater among-year variability.

Conclusions
Understanding the geographical distribution of the phenotypic and genetic variation within and between populations is crucial for assessing species vulnerability to projected environmental changes. Indeed, growing evidence suggests that polymorphism increases the ability of a species to persist in the face of increasing disturbance. Despite the contribution of stochastic and historical factors to shaping the current diversity of C. lapponica at large spatial scales, the level of polymorphism was highest in populations inhabiting the cold and variable environments associated with high latitudes and altitudes. This pattern likely emerged because differences in the climatic requirements of coexisting color morphs provide survival benefits for the populations living in these conditions, while also enabling them to cope with continuing warming and with an increasing frequency of extreme weather events. The spatiotemporal variation in the frequencies of melanic morphs was in line with the thermal melanism hypothesis. The loss of phenotypic diversity, including the decline of melanic beetles in regions with cold climates, adds a new dimension to global biodiversity loss and may indicate an increased vulner-ability of persisting populations to future environmental disturbances.

Acknowledgments
Collection of the substantial part of the data and the completion of the study were supported by the Academy of Finland (projects 122133, 122144, 122180, 127047, 203156, 208016, 214653, 268124, 276671, 311929, and 316182

Data availability
The data supporting the results are included in this article as Supporting Information Data S1 and Data S2.

Supporting Information
Additional supporting information may be found online in the Supporting Information section at the end of the article.  .  Table S1. Characteristics of color morphs of Chrysomela lapponica. Table S2. Distribution of recorded beetles by countries. Table S3. The comparisons between linear and quadratic regression models explaining variation in Shannon diversity index. Table S4. Performance of models explaining variation in the Shannon diversity index and in the proportion of melanic beetles: the outcomes of the forward step-wise regression analysis. Table S5. Geographic patterns in proportions of color morphs in large samples of Chrysomela lapponica. Table S6. Results of db-RDA with forward selection for frequencies of beetle color morphs. Table S7. The comparisons between linear and quadratic regression models explaining variation in the proportion of melanic beetles.
Data S1. Data on studied individuals of Chrysomela lapponica.
Data S2. Data on samples of Chrysomela lapponica used in the analysis of patterns in morphological diversity.