The Fluctuating Asymmetry of the Butterﬂy Wing Pattern Does Not Change along an Industrial Pollution Gradient

: The rapid and selective responses to changes in habitat structure and climate have made butterﬂies valuable environmental indicators. In this study, we asked whether the decline in butterﬂy populations near the copper-nickel smelter in Monchegorsk in northwestern Russia is accompanied by phenotypic stress responses to toxic pollutants, expressed as a decrease in body size and an increase in ﬂuctuating asymmetry. We measured the concentrations of nickel and copper, forewing length, and ﬂuctuating asymmetry in two elements of wing patterns in Boloria euphrosyne , Plebejus idas , and Agriades optilete collected 1–65 km from Monchegorsk. Body metal concentrations increased toward the smelter, conﬁrming the local origin of the collected butterﬂies. The wings of butterﬂies from the most polluted sites were 5–8% shorter than those in unpolluted localities, suggesting adverse effects of pollution on butterﬂy ﬁtness due to larval feeding on contaminated plants. However, ﬂuctuating asymmetry averaged across two hindwing spots did not change systematically with pollution, thereby questioning the use of ﬂuctuating asymmetry as an indicator of habitat quality in butterﬂy conservation projects.


Introduction
Butterflies are valuable environmental indicators because they rapidly and selectively respond even to subtle changes in habitat structure and climate [1,2] (and references therein). Their colorful appearance and aesthetic appeal have attracted attention to these insects from both the general public and decision-makers [3] and have made butterflies focal objects of multiple national and international monitoring schemes [4,5], as well as of citizen science projects (e.g., www.iNaturalist.org (accessed on 20 March 2021)). Because of this attention, butterflies are one of the best studied groups of insects-the Global Biodiversity Information Facility at the beginning of 2019 included over 7.5 million complete and unique records of butterflies from around the globe [6].
Environmental pollution was recently identified as one of the five main drivers of biodiversity loss [7]. Along with direct (toxic) impacts on biota, pollution also changes habitat quality, which in turn affects biodiversity [8,9]. Therefore, pollution-induced changes in the diversity and abundance of butterflies are highly likely; however, the available evidence remains surprisingly scarce. Industrial pollution was hypothesized to be one reason for the declines in butterfly numbers in some regions [10,11], but mechanisms behind these declines have not been deciphered [12]. These mechanisms may be rather non-trivial; for example, the loss of thermophilous butterflies from the Netherlands was explained by microclimatic cooling due to the increased plant growth that arose due to a combined effect of increased nitrogen deposition and climate warming [13].
Counts of day-active moths and butterflies, conducted in 1991-1993 around the coppernickel smelter in Monchegorsk in northwestern Russia, demonstrated that butterflies were most abundant in slightly polluted forests but were almost absent in the industrial barrens (bleak open landscapes with small patches of vegetation surrounded by bare land) next to the smelter [14]. Although many insect populations decline with increases in pollution [12,15,16], thus making this pattern seem trivial, this decline may have emerged for various reasons. These include (but are not limited to) direct toxicity of pollutants for insects and/or multiple indirect effects, acting, e.g., through the disappearance of forest habitats and the resulting changes in microclimate.
In this study, we tested the hypothesis that the decline in butterfly populations near the smelter can be predicted based on phenotypic differences between the exposed and control populations. The concentrations of nickel and copper in plants growing near Monchegorsk in the 1990s were as high as 600-800 µg g −1 (i.e., 50 times greater than the regional background: [17]). We therefore predicted that the metal concentrations in butterflies would increase with increasing proximity to the polluter, as the larvae nearer to the polluter would have fed on more highly contaminated plants. The subsequently increased body concentrations of toxic metals are likely to cause physiological stress, which can be detected by phenotypic changes [18], including decreases in insect size [16]. Finally, pollution has been repeatedly reported to disturb the developmental stability of different organisms [19][20][21][22]. These observations led us to predict a greater fluctuating asymmetry (FA) of the wing pattern in butterflies from polluted sites than from unpolluted sites. We verified our predictions by measuring concentrations of nickel and copper, forewing length, and FA of the two wing spots in three butterfly species collected at different distances from the polluter.

Study Region and Study Sites
The town of Monchegorsk (67 • 56 N, 32 • 55 E), the site of a large copper-nickel smelter, is located 150 km south of the northern tree line. A century ago, this area was covered by impenetrable forests of Norway spruce (Picea abies (L.) Karst.) and Scots pine (Pinus sylvestris L.).
The mean temperature at Monchegorsk is −13.8 • C in January and 14.1 • C in July, and the mean annual precipitation is 561 mm; the frost-free period ranges from 50 to 100 days. The smelter started its operations in 1939, and the peak annual values of its emissions amounted to 278,000 metric tonnes (t) of sulphur dioxide in 1983 and 13,150 t of non-ferrous metals in 1987; the current annual emissions are close to 40,000 t of sulphur dioxide and 1000 t of metals [9,23]. Along with these main pollutants, the smelter emits dozens of other potentially toxic substances [24].
The pollution has transformed over 250 km 2 of the previously forested areas around the smelter into industrial barrens. Both Scots pine and Norway spruce are practically absent in barren habitats, and low-stature (0.2-3 m tall) mature mountain birches (Betula pubescens var. pumila (L.) Govaerts) growing 5-15 m apart dominate the landscape. In the intermediate zone, the top canopy of sparse forests is formed by Norway spruce that show visible signs of damage (upper parts of crowns are dead, needle longevity is low) and mountain birches, whereas field-layer vegetation is sparse. Mountain birches are common also in Norway spruce forests that are visibly unaffected by pollution, where Empetrum nigrum ssp. hermaphroditum (Hagerup) Böcher, Vaccinium myrtillus L., and V. vitis-idaea L. dominate in the dense field layer. In dry pine forests, the field layer includes Calluna vulgaris (L.) Hull, whereas V. uliginosum L. and Rhododendron tomentosum Harmaja (=Ledum palustre L.) are common in wet microsites [9,23].

Study Objects
We explored three of the most common butterfly species in the study region. Boloria euphrosyne (L.), an orange butterfly with black spots, has a wingspan of 38-46 mm and is on the wing during the first half of the summer. Its larvae feed on Viola sp., R. tomentosum, and V. uliginosum. The second species, Plebejus idas (L.), has a wingspan of 17-28 mm and flies during the second half of the summer. The males of this species have iridescent blue wings, whereas the females have brown wings with orange spots. The larvae of P. idas feed on C. vulgaris, V. uliginosum, E. nigrum, and on various Fabaceae species (most of which are absent from our study sites). The third species, Agriades optilete (Knoch), has a wingspan of 23-29 mm and flies simultaneously with P. idas; it also resembles P. idas in wing color and pattern, although the females of A. optilete may be almost as blue as the males. The larvae of A. optilete feed on V. uliginosum, V. vitis-idaea, V. myrtillus, V. oxycoccos L., E. nigrum, and Andromeda polifolia L. All these butterfly species hibernate as larvae.
We collected butterflies using standard insect collecting nets in 17 sites located 1-65 km from the smelter ( Figure 1, Table 1) during the periodic assessment of insect diversity in the Monchegorsk pollution gradient from 1989-2001. Each year, we visited each site three to 10 times. During each visit, we collected all moths and butterflies seen during two hours. All collected butterflies (Supplementary Materials, Data S1) were pinned, and most of them were spread. We aimed at measuring 10-15 butterflies for each species-by-site combination; when the number of collected butterflies exceeded this desired value, we randomly selected individuals for this study from the available specimens. After measurements, the larger part of the collected butterflies was donated to the Zoological Museum, University of Helsinki, Finland. We collected butterflies using standard insect collecting nets in 17 sites located 1-65 km from the smelter ( Figure 1, Table 1) during the periodic assessment of insect diversity in the Monchegorsk pollution gradient from 1989-2001. Each year, we visited each site three to 10 times. During each visit, we collected all moths and butterflies seen during two hours. All collected butterflies (Supplementary Materials, Data S1) were pinned, and most of them were spread. We aimed at measuring 10-15 butterflies for each species-by-site combination; when the number of collected butterflies exceeded this desired value, we randomly selected individuals for this study from the available specimens. After measurements, the larger part of the collected butterflies was donated to the Zoological Museum, University of Helsinki, Finland.  Table 1.  [9], and unpublished; ¶ BWC, secondary birch-and willow-dominated community; DF, slightly damaged spruce forest; IB, industrial barren; SDF, severely damaged spruce forest; UF, undamaged spruce forest.

Measurements
We measured forewing length with a ruler (to the nearest 0.5 mm) as the distance between the base of the costal wing margin and the wing apex. This characteristic was only used to quantify the size of butterflies, because the locomotory traits are more buffered against developmental perturbations than non-locomotory traits are [25] and because precise measurements of wing size would require detaching the wings from the butterfly bodies. The FA was assessed in two spots on the underside of the hind wings ( Figure 2). The size of the selected spots was measured (to the nearest 0.05 mm) on both the left and the right wings using a stereomicroscope micrometer. In B. euphrosyne, we measured the distance between the points at which the external margins of the selected spots approached the adjacent veins ( Figure 2a), whereas in P. idas and in A. optilete, we measured the maximum distance between the opposite spot margins (Figure 2b,c). In 20 butterflies of each species, the measurements were repeated 12 months after the first measurement (Data S2). All measurements were performed by V.Z., blindly with respect to the collection site.  [9], and unpublished; ¶ BWC, secondary birch-and willow-dominated community; DF, slightly damaged spruce forest; IB, industrial barren; SDF, severely damaged spruce forest; UF, undamaged spruce forest.

Measurements
We measured forewing length with a ruler (to the nearest 0.5 mm) as the distance between the base of the costal wing margin and the wing apex. This characteristic was only used to quantify the size of butterflies, because the locomotory traits are more buffered against developmental perturbations than non-locomotory traits are [25] and because precise measurements of wing size would require detaching the wings from the butterfly bodies. The FA was assessed in two spots on the underside of the hind wings ( Figure 2). The size of the selected spots was measured (to the nearest 0.05 mm) on both the left and the right wings using a stereomicroscope micrometer. In B. euphrosyne, we measured the distance between the points at which the external margins of the selected spots approached the adjacent veins (Figure 2a), whereas in P. idas and in A. optilete, we measured the maximum distance between the opposite spot margins (Figure 2b,c). In 20 butterflies of each species, the measurements were repeated 12 months after the first measurement (Data S2). All measurements were performed by V.Z., blindly with respect to the collection site.

Chemical Analyses
Concentrations of nickel and copper were measured in P. idas and A. optilete specifically collected for this purpose 1-170 km from Monchegorsk in 2003-2004. The third species, B. euphrosyne, was found too infrequently in the heavily polluted sites to obtain a representative sample of this species. The butterflies were handled and packed in paper envelopes with clean forceps to avoid cross-contamination. The randomly selected butterflies were digested in a 'Kjeldatherm' (Gerhardt, Germany) block digestion system with species, B. euphrosyne, was found too infrequently in the heavily polluted sites to obtain a representative sample of this species. The butterflies were handled and packed in paper envelopes with clean forceps to avoid cross-contamination. The randomly selected butterflies were digested in a 'Kjeldatherm' (Gerhardt, Germany) block digestion system with analytical-grade nitric acid (Merck). The residue was diluted to 5 mL with Milli-Q deionized water and filtered through paper (pore size 1.5-2 µm). Concentrations of Ni and Cu were measured by means of graphite furnace atomic absorption spectrometry using an AAnalyst 800 atomic absorption spectrometer (Perkin Elmer, Waltham, MA, USA) equipped with a THGA graphite furnace with Zeeman-effect background correlation and autosampler AS 800. The analyses were conducted at the Institute of North Industrial Ecology Problems, Apatity, Russia.

Data Analysis
The FA values were calculated as follows: FA = 2 × abs(WL − WR)/(WL + WR), where WL and WR refer to the spot size measured on the left and right wings of the same butterfly. The use of this index, labelled FA2 by Palmer and Strobeck [26], is justified by the significant positive correlation between the absolute difference in trait measurements between the left and right wings and an average trait size (r = 0.27, n = 1452 measurements, p < 0.0001).
We explored the data on 60 butterflies, which were measured twice, for the presence of FA and directional asymmetry (DA) relative to the measurement error by means of mixed-model ANOVA with the butterfly wing (right or left) considered as a fixed factor and the individual butterfly as a random factor. The reproducibility of the measurements was quantified by the index ME5 = (MSi − MSm)/(MSi + MSm), where MSi and MSm are the interaction and error mean squares from a side × individual ANOVA for two measurements of each trait in each individual [26]. This index expresses FA variation as a proportion of the total variation between sides, which includes variations due to both FA and measurement error. When mixed-model ANOVA revealed a significant difference between the right and left sides (i.e., the existence of DA), we compared the DA value with the FA4a index (FA4a = 0.798 √ var (R − L)), as suggested by Palmer and Strobeck [26]. Sources of variation in metal concentrations, wing length, and FA (averaged between two measured spots) were explored by means of a mixed-model ANCOVA (SAS GLIMMIX procedure, type 3 tests: [27]). We considered butterfly species and sex as fixed effects, the log-transformed distance from the smelter (a proxy of pollution load: Kozlov et al. [9]) as a covariate, and the site as a random effect. Following the practice commonly accepted in observational studies [28][29][30], we preferred to use the distance from the polluter rather than the concentration of one of the main pollutants (nickel or copper) in plant foliage or in a litter. This was done to avoid misinterpretation of our results, because we do not know which of the pollution-related factors (which all change with the distance from the smelter in a coordinated manner) may have affected the butterflies. In our study area, concentrations of metal pollutants in plants and in soil strongly correlate with logtransformed distance from the polluter [9,23]. To facilitate accurate F-tests of the fixed effects, we adjusted the standard errors and denominator degrees of freedom as described by Kenward and Roger [31]. The significance of a random factor was evaluated by testing the likelihood ratio against the Chi-squared distribution [32]. The associations between the study variables were explored by calculating the Pearson product-moment correlation coefficients (SAS CORR procedure [27]).

Results
The concentrations of nickel and copper did not differ between P. idas and A. optilete, and the concentrations increased by factors of 12 and 4, respectively, with increasing proximity to the smelter (Table 2, Figure 3).  The forewing length in the three butterfly species similarly decreased toward the smelter (Table 3), being 5%-8% smaller in the most polluted sites than in the unpolluted sites ( Figure 4). However, the species-specific correlations between forewing length and distance from the smelter were significant only in B. euphrosyne and P. idas (Figure 4a,b). Table 3. Sources of variation in forewing length and fluctuating asymmetry of butterflies (SAS GLIMMIX procedure, type 3 tests). The side × individual interaction was significant for both measured spots in all butterfly species (Table 4), thereby confirming the existence of FA in spot size and our ability to quantify it using measurements of the given accuracy. The DA was significant in the first spot in B. euphrosyne and P. idas, but not in A. optilete (Table 4). However, the DA in both B. euphrosyne and P. idas was smaller than the FA4a index, suggesting that DA's con- The forewing length in the three butterfly species similarly decreased toward the smelter (Table 3), being 5-8% smaller in the most polluted sites than in the unpolluted sites ( Figure 4). However, the species-specific correlations between forewing length and distance from the smelter were significant only in B. euphrosyne and P. idas (Figure 4a,b). Table 3. Sources of variation in forewing length and fluctuating asymmetry of butterflies (SAS GLIMMIX procedure, type 3 tests). the distance from the smelter (Table 3) or by nickel concentrations in birch leaves (results not shown). The spot-specific correlations with distance were generally non-significant (Figure 5a-e), with the exception of the second spot in A. optilete, which showed a decreasing FA with increasing proximity to the smelter (Figure 5f).    The side × individual interaction was significant for both measured spots in all butterfly species (Table 4), thereby confirming the existence of FA in spot size and our ability to quantify it using measurements of the given accuracy. The DA was significant in the first spot in B. euphrosyne and P. idas, but not in A. optilete (Table 4). However, the DA in both B. euphrosyne and P. idas was smaller than the FA4a index, suggesting that DA's contribution to the total variation in size of the first spot was small and can therefore be neglected.

Source of Variation
The FA of the two spots did not correlate either with each other (B. euphrosyne: r = 0.05, n = 84, p = 0.65; P. idas: r = 0.08, n = 312, p = 0.15; A. optilete: r = 0.10, n = 298, p = 0.09) or with forewing length (r = −0.10 . . . 0.17, n = 86 . . . 314, p = 0.08 . . . 0.96). The FA averaged across the two spots varied among the study sites, but this variation was not explained by the distance from the smelter (Table 3) or by nickel concentrations in birch leaves (results not shown). The spot-specific correlations with distance were generally non-significant (Figure 5a-e), with the exception of the second spot in A. optilete, which showed a decreasing FA with increasing proximity to the smelter (Figure 5f).

Discussion
Insects living in polluted habitats often accumulate heavy metals, including nickel and copper [12]. The concentrations of nickel in butterflies were of the same range as in other plant-feeding insects collected from the same sites [33,34], whereas the concentrations of copper amounted to one-fifth of those in other insects. Importantly, the metal concentrations in the butterflies were much smaller than those reported in plants [17,35]. This result is in line with the earlier finding that larvae of Eriocrania semipurpurella excrete 90%-95% of the nickel and 50%-80% of the copper consumed from contaminated birch leaves [33]. The metals are excreted with feces [33], larval exuviae, and pupal shells [36,37]. More generally, our findings confirm an opinion by Laskowski [38] that the accumulation of heavy metals in food webs (biomagnification) is an exception rather than the general rule.

Discussion
Insects living in polluted habitats often accumulate heavy metals, including nickel and copper [12]. The concentrations of nickel in butterflies were of the same range as in other plant-feeding insects collected from the same sites [33,34], whereas the concentrations of copper amounted to one-fifth of those in other insects. Importantly, the metal concentrations in the butterflies were much smaller than those reported in plants [17,35]. This result is in line with the earlier finding that larvae of Eriocrania semipurpurella excrete 90-95% of the nickel and 50-80% of the copper consumed from contaminated birch leaves [33]. The metals are excreted with feces [33], larval exuviae, and pupal shells [36,37]. More generally, our findings confirm an opinion by Laskowski [38] that the accumulation of heavy metals in food webs (biomagnification) is an exception rather than the general rule.
The increased concentrations of nickel and copper in butterflies collected near Monchegorsk indicate their local origin, because their larvae have consumed polluted food. The concentrations of nickel and copper in leaves of different plants growing in the most polluted sites near the Monchegorck smelter reached 500-800 µg g −1 [17,35].
The earlier experiments demonstrated that 50-200 µg g −1 of copper added to the larval diet decreased the pupal weight of the oriental leafworm moth, Spodoptera litura (F.), to 80% of the control value [39], whereas 100 µg g −1 of nickel did not affect wing size in the cabbage white butterfly, Pieris rapae (L.) [40]. At the same time, the wing length of the autumnal moth Epirrita autumnata (Bkh.) did not change along the Monchegorsk pollution gradient [41], and the wing length of the brassy tortrix Eulia ministrana (L.) increased by 10% [30]. Thus, although metal toxicity is a likely reason for the recorded decrease in the size of butterflies accompanying the increase in pollution, experimental data are required to uncover the immediate reason behind the observed effect.
The abundance of butterflies near the smelter was reduced to 10% relative to unpolluted and slightly polluted sites [14]. Similarly, the biomass of dwarf shrubs, which include food plants of B. euphrosyne, P. idas, and A. optilete, in the most polluted sites was reduced to 0-20% of unpolluted controls [23]. Thus, the abundance of butterflies decreased proportionally to the decline in larval food resources. Nevertheless, the densities of butterflies in polluted sites were so small (0.3 butterflies seen during an hour [14]), that we do not consider larval starvation to be a possible reason for the small size of the butterflies. Instead, we suggest that this reduction in body size results from the toxicity of metal pollutants accumulated in larval host plants [42,43], combined with the physiological costs of metal detoxification and excretion. Metal toxicity could increase larval mortality [43], whereas a decrease in body size reduces the fecundity of female butterflies [44,45].
The decline in the butterflies in a Dutch nature reserve during the 1990s was attributed to the adverse effects of heavy metals on the nectar plants [46]. However, the proportion of flowering patches of V. myrtillus and V. vitis-idaea near Monchegorsk was higher, and these patches had more flowers compared to the unpolluted localities [47]. This means that nectar availability to adult butterflies decreases with an increase in pollution at a smaller rate than food availability to the larvae. Therefore, we doubt that a shortage of nectar could be the primary reason for the decline in butterflies near the Monchegorsk smelter.
Based on the increased body content of metals and the significant reduction in wing length, we conclude that butterflies living near Monchegorsk suffer from adverse environmental conditions. Following the arguments of Waddington [48], these adverse conditions are likely to disturb developmental stability, resulting in an increase in FA. Due to its seeming simplicity, this concept gave rise to a wealth of studies that explored the impacts of environmental stress on FA in different organisms [20,49,50].
Contrary to optimistic expectations [51][52][53], the accumulated evidence does not confirm that FA consistently increases in living beings facing unfavourable conditions during their development [22,54]. The studies of insects detected both positive associations between pollution and FA [55] and absences of correlation between these variables [30,56]. We estimate that, across organisms and traits, the support for this hypothesis was found in no more than half of the examined data sets [57]. The present study further confirms that life in extreme environments does not necessarily result in an increase in FA.
Only about a half of the previous studies measured FA in more than one trait per individual, and only a small fraction of these studies statistically combined information across traits [50] or explored among-trait correlations [25]. The single-trait approach reflects an implicit assumption of evolutionary models that any trait showing FA is suitable for the analysis of environmental impacts on developmental stability. However, we found that the FA values in the two wing spots were not related to each other. This finding, which is in line with earlier studies of FA in multiple traits of butterflies and moths [30,58,59], seemingly contradicts the outcomes of meta-analysis, which confirmed the existence of the organism-wide response in FA [25]. However, the among-trait correlation in FA is so small (an average effect size of 0.05 [25]) that it can hardly be detected in an individual study, keeping in mind the obvious constraints of sample size.
Although the FA of different traits is practically uncorrelated, we are not aware of any theoretical model that predicts whether the FA in a particular trait of a given species will change in response to a specific environmental factor [57]. Nevertheless, the popularity of the theory linking environmental stress with a decrease in developmental stability and an increase in FA is so high that researchers exploring FA are forced to provide explanations for each 'negative' result, erecting hypotheses as to why their results do not fit the theoretical predictions (e.g., [60,61]). At the same time, the confirmatory evidence has long been published without deep examination of the research methodology used to arrive at these 'positive' results.
Two decades ago, when discussing the criticism directed toward the theory outlined above, Gangestad and Thornhill [62] (p. 414) wrote: " . . . bodies of evidence could have turned out quite differently and, hence, falsified prevailing notions about developmental imprecision and asymmetry. That that they did not but, rather, fit nicely with theory is a strange coincidence if those notions are entirely mistaken". We now know that this coincidence may have resulted from confirmation bias, defined as the tendency of humans to seek out evidence in a manner that confirms their beliefs and hypotheses. We experimentally demonstrated that when scientists expected to find high FA in some samples, the results of their measurements confirmed their expectations [63]. False discoveries of this type can be avoided by blinding [64], which we regard as obligatory in future studies of FA.
We suggest that authors and reviewers always ensure that (1) samples are collected either randomly or blindly with respect to the expected result; (2) the measurements of FA are conducted blindly (i.e., the measurer is not aware of the object's origin or of the hypothesis tested); and (3) at least a portion of the objects are remeasured and the FA is tested against the measurement error. We recommend that studies violating these criteria not be published or, if published, their outcomes not be included in future meta-analyses exploring the sources of variation in the FA of living beings. We also encourage the authors to publish their 'negative' and inconclusive results when they were obtained using adequate methodology, and to open up their file drawers, as Lane et al. [65] did. This is the only way to make the publication portfolio more representative of the actual findings.
In conclusion, the increased concentrations of nickel and copper in butterflies collected from polluted sites suggest that these butterflies are of local origin. The decline in abundance of butterflies with an increase in pollution, as observed near the Monchegorsk smelter, is likely driven by a combination of the pollution-induced decrease in biomass of the larval host plants and the toxicity of the metal pollutants accumulated in these plants. However, we did not detect any increase in the FA of the wing spots in the butterflies persisting in the polluted sites. This finding adds to a growing body of 'negative' results (in terms of the hypothesis predicting an increase in FA in response to environmental stress) and questions the use of FA as an indicator of habitat quality in butterfly conservation projects.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/sym13040626/s1, Data S1: Results of the first measurement, Data S2: Results of the second measurement.