Monitoring climate changes and habitat degradation in threatened species without negative impact to the populations can pose a considerable challenge. A rare chance to test the morphological response of wing shape and size to environmental factors on the mountain Apollo (Parnassius apollo) collected from 1938 to 1968 at a single location – Strečno mountain pass, N Slovakia presented itself in a historical collection. The canonical variate analysis showed a significant shift from a narrower to broader forewing, with more extremes in either extra broad or narrow forewings in the post- 1960 population. Analysis of existing data was conducted to determine the possible factors affecting this change. Generally, the comparative statistics of temperature and precipitation to morphology of individuals and their fluctuating asymmetry showed no, or weak, correlations. Two extreme weather events (ECEs), identified using the historical weather data, show no correlation of wing morphology to these events. Although no strong correlations can be drawn in case of the available weather data and morphology, the results of this study can be connected to strong anthropogenic effects of a large-scale road development project taking place in the vicinity of the collection site starting in November 1959 causing changes in the available habitat and therefore a shift in the wing morphology.

Anthropic impact, butterfly, geometric morphometry, historical collection, Parnassius Apollo, wing

Sudden or continual alterations to the environment caused by natural or anthropogenic processes have been of interest to the scientific community for some time (Stern et al. 1992). Monitoring these environmental changes has become an essential component of applied ecology, using suitable species (bioindicators) and their assemblages to analyze environmental stress and disturbance (Dubovský et al. 2010; Zvaríková et al. 2016). There is a large body of work concerning the negative effects in populations of arthropod fauna, including butterflies, to various environmental changes (Van Swaay and Martin 1999; Descimon et al. 2005; Wenzel et al. 2006; Nakonieczny et al. 2007; Van Dyck et al. 2009; da Rocha et al. 2010; Gibbs et al. 2011; Wallner et al. 2013). Due to the more extreme nature of these changes in mountainous regions the native populations face even higher risk of population reduction and potential extinction (Nogués-Bravo et al. 2007). The same risks must be expected for their host and nectar plants to which their life cycle is bound (Engler et al. 2011).

Among bioindicators, detectors may provide an important material to detect response based on their body morphology (Blake et al. 1994; Findlay and Houlahan 1997; Clarke et al. 2000; Spellerberg 2005; Zvaríková et al. 2016). Generally, all organisms exhibit a certain degree of morphological (phenotypic) plasticity (Ananthakrishnan and Whitman 2005; Hoffmann et al. 2005; Sukhodolskaya and Saveliev 2014), which describes the capacity of a genotype to display a range of phenotypes in response to changes in the environment (Garland and Kelly 2006; Whitman and Agrawal 2009) allowing organisms to maintain high fitness in the face of environmental heterogeneity (Pigliucci et al. 2006). The monitoring of these changes in morphology is heavily focused on the symmetries and asymmetries of individuals in attempts to describe the underlining factors of these changes (Ananthakrishnan and Whitman 2005; Zvaríková et al. 2016).

Asymmetry in individuals’ morphology is often used to illustrate isolation effects, reduction in genetic diversity, and lower fitness of populations or environmental stress (Valen 1962; Parsons 1990; Frankham 2005; Allendorf et al. 2013; Beasley et al. 2013; Woods et al. 1999). Within previous research, population genetic diversity was negatively correlated to the level of asymmetry on multiple occasions (Vrijenhoek and Lerman 1982; Blanco et al. 1990; Hoelzel et al. 2002). Higher levels of morphological asymmetry in birds corresponded with higher mortality during extreme climatic events (Brown and Brown 1998). In contrast to these claims are examples of naturally highly inbred populations where the benefits of outbreeding are not always maintained (Fountain et al. 2015). Fluctuating asymmetry (FA) of individuals’ morphology is often used as an indicator for population fitness, although the complex nature of the origin and expression of fluctuating asymmetry (FA) or lack thereof shows that a unification of these results may be impossible, and the problem therefore requires a case-by-case approach (Kaeppler 2012; Windig et al. 2000).

A reliable evaluation of insect morphology vs. environment interaction requires a long-term monitoring of a suitable model species with sufficient information on its autecology as well as environmental disturbance details (Webster and Sheets 2010). Using historical collections, we offer an analysis of the mountain Apollo, Parnassius apollo (Linnaeus, 1758), a xeromontainous papilionid butterfly with diminishing distribution across Europe. The species has been in long-term decline since the 1900s due to lack of land management and global climate change (Todisco et al. 2010). Current populations are now usually present in small, often isolated patches, unable to colonize new habitats (Fred et al. 2006). Although its physiological predispositions would allow the Apollo butterfly to migrate longer distances, its successful reproduction remains restricted to the distribution of its larval host plant (Sedum sp.) located on open rock formations (Brommer and Fred 1999). The survival of these small populations was hypothesized to be due to the ability to maintain a low long-term effective population, which may result from a strong historic bottleneck or a founder event (Habel et al. 2009). Due to the continuous decline in population this species has been listed in the IUCN Red List of Threatened Species, the appendix II in CITES, as well as the annex IV of Habitats Directive.

The isolation effect on the phenotype of Parnassius apollo was previously tested on western European populations, showing similar levels of asymmetry throughout the region (Habel et al. 2012), however higher levels of asymmetry were presented by Adamski and Witkowski (2002) in a post-bottleneck population when compared to the pre-bottleneck individuals. The populations of the Carpathian region (Eastern/Central European) were determined as genetically homogeneous (Todisco et al. 2010), although multiple subspecies based on morphology are considered (Hrubý 1964; Capdeville 1978; Kizek 1997).

Although the pressures on selection in butterfly wing morphology remain largely unknown, the strongest factors seem to be habitat, predators and sex-specific behavior (Le Roy et al. 2019). Contrasting habitats (Jugovic et al. 2018), restoration attempts (Sivakoff et al. 2016) or landscape structure (Merckx and Van Dyck 2006) show impact on the resulting wing morphology in butterflies. The material of mounted Parnassius apollo specimens by a single collector at a location in the vicinity of an important north Slovakian trade and travel route presents an ideal dataset for the analysis of ecological (temperature and precipitation) and environmental (stress and disturbance) factors along a 30-year long period in the mid-20^th century. A thorough analysis of the correlation strength of weather and extreme climatic events to the resulting changes in morphology was conducted. We hypothesize that the anthropic impact at the collection sites could induce a change in the morphology of the P. apollo wing observed through the geometric morphometry.

All 506 Parnassius apollo specimens used for this analysis belong to the historical collection (1938–1968) of the Slovak National Museum (Bratislava, Slovakia) and were mounted by Ján Zelný. Since the whole collection has been the work of a single collector it presents an ideal study opportunity, reducing the variable of multiple collectors or preparation techniques used. The material refers to a single location at Strečno mountain pass, N Slovakia (49°10'07.6"N, 18°52'51.6"E) (Fig. 1), at the altitude of 450 m a. s. l.

Figure 1.

Location of Strečno mountain pass in European context and a close-up of the collection site on a historic map from 1955.

Digitization and morphometric analyses

Each individual was photographed from the dorsal side using a digital camera (Canon 60D, 50 mm lens) under standardized light conditions, with the added size standard (5cm with 1mm increments) fixed at the height of the wing. Geometric morphometric analysis was based on 17 landmarks (10 landmarks on forewing and 7 on hindwing), situated terminally or on vena intersection (Fig. 2) to assure the repeatability of the landmarks on the highest numbers of individuals (Habel et al. 2012).

Figure 2.

Left fore and hindwing of P. apollo used for the morphometric analysis, including the morphometric landmarks Forewing: 1 – R4 and R5 intersection, 2 – R4 terminally, 3 –R5 terminally, 4 –M2 terminally, 5 –M3 terminally, 6 –Cu1 terminally, 7 –Cu2 terminally, 8 –A2 terminally, 9 – Discal cell and Cu1 intersection, 10 – Discal cell and Cu2 intersection. Hindwing: 11 –Rs terminally, 12 – M1 terminally, 13 – M2 terminally, 14 – M3 terminally, 15 – Cu1 terminally, 16 – Discal cell and M1 intersection, 17 – Discal cell and Rs intersection.

A database (.TPS file) of all digitized specimens was created using tpsUtil64 (version 1.70) software (http://life.bio.sunysb.edu/morph/soft-utility.html, accessed 23.11.2018) and imported to tpsDIG2×32 (version 2.26) software (http://life.bio.sunysb.edu/morph/soft-dataacq.html, accessed 23.11.2018) where all landmarks were set. TPS file containing the landmarks for each specimen was imported to MorphoJ software (version 1.06d) (Klingenberg 2008, 2011).

To address the potential influence of digitization or setting of the landmarks on our results, individuals were tested for digitization and measurement error. To test for digitization error, 50 specimens were randomly selected, photographed twice, measured and the results compared. The measurement error was addressed by measuring the whole dataset twice by the same person and compared (Palmer and Strobeck 1986; Klingenberg and McIntyre 1998).

The weather data was obtained from the archives of the Slovak Hydrometeorological Institute (SHMI) for the locations closest to the collection site. The precipitation data was recorded at location Vrútky, Slovakia (49°06'42.7"N, 18°55'26.2"E), 7 km south from the collection site and Varín, Slovakia (49°12'05.4"N, 18°52'12.4"E) 3 km north of the collection site. The closest site collecting temperature data was maintained by the Slovak Army at the Žilina Airport, Slovakia (49°13'59.1"N, 18°36'47.4"E) located 20 km west of the collection site. All data on temperature prior to 1944 was probably lost, therefore comparative statistics between morphometric data and temperature was conducted only for the dataset 1944–1968. The data of monthly precipitation and temperature, and yearly precipitation and temperature, were correlated to symmetric and asymmetric components, centroid size and Procrustes FA scores of Procrustes coordinates exported from MorphoJ and analyzed using “corrr” package in R (RStudio Team 2015). The extreme climatic events (ECEs) were selected according to the ability of the collected weather data to display these events. In the case of extreme temperatures, the highest monthly and yearly temperature did not exceed the averages for the location and we therefore did not consider high temperatures from 1944 – 1968 to display any ECEs. In the case of low temperatures, February 1956, according to the SHMI, saw the lowest monthly temperature ever recorded for the location during the measured period 1944–2019. The heavy rainfall and flooding in 1958 can also be supported by the collected data. Based on this information we selected the low monthly February 1956 temperatures and 1958 heavy rainfalls and flooding as the two most important ECEs during our research period.

Data on the anthropic impact and road construction timeline was pieced together from articles and archives by The News Agency of the Slovak Republic (TASR), using mostly contemporary photographs with descriptions.

All 506 Parnassius apollo specimens (Tab. 1), including 342 males and 164 females, used for this analysis were checked for strength of measurement and digitization error, with negligible results when compared to the morphological changes. In accordance with no significant size, shape or fluctuating asymmetry differences of hind or forewings between sexes (Tab. 2), the specimens were pooled together for further analyses.

Testing for digitization error by photographing 50 individuals under the same light conditions and digitization setup, no significant difference in landmark placement was detected (ANOVA, digitization error of forewing p = 0.76, hindwing p = 0.56). The same result of non-significant difference was recorded for the measurement error, remeasuring all 506 individuals twice (ANOVA, measurement error of forewing p = 0.59, hindwing p = 0.45).

Wing shape change patterns

The Mann-Whitney’s pairwise comparison, to test the significance of the differences in wing shape between each year of the time series (Appendix 1), and the Canonical variate analysis (CVA) show a significant change in the wing morphology starting with 1960 for the forewing and 1961 for the hindwing, continuing until the end of our dataset. Canonical variate 1 (CV1) and 2 (CV2) (Figs 3–6) describes 81% (for the forewing) and 76% (for the hindwing) of the total variance in wing morphology.

Figures 3, 4.

Mountain Apollo forewing shape changes displayed on the canonical variate 1 and 2. The boxplots cover the upper and lower 25%, separated by a line showing median inside the box, whiskers cover remaining upper and lower 25%, the white circles represent outliers and black stars are extreme values.

Figures 5, 6.

Mountain Apollo hindwing shape changes displayed on the canonical variate 1 and 2. The boxplots cover the upper and lower 25%, separated by a line showing median inside the box, whiskers cover remaining upper and lower 25%, the white circles represent outliers and black stars are extreme values.

The shift observed starting in 1960 was tested using 2 factor analysis of variance where we tested the artificial pre/post 1960 datasets and the natural variation between years using CV1 (Tab. 3). To measure the shift in morphology observed in Figs 3–6 a breakpoint regression analysis was tested on CV 1 of the fore- and hindwing (Figs 7, 8).

Table 3.

Download as

CSV

XLSX

Canonical variate 1 (CV1) used in the analysis of variance for comparison of natural variance in wing morphology between years and the shift observed around 1960.

		Df	Sum sq	Mean sq	F value	Pr(>F)
CV1 forewing	pre vs post 1960	1	381.9	381.9	381.922	< 2 x 10-16	***
	Years	19	60.7	3.2	3.194	4.81 x 10-6	***
	Residuals	989	989	1
CV1 hindwing	pre vs post 1960	1	117.6	117.6	117.601	< 2 x 10-16	***
	Years	19	64.8	3.41	3.412	1.14 x 10-6	***
	Residuals	993	993	1

Figures 7, 8.

Breakpoint regression analysis of the CV1 for the fore- and hind wing displaying the estimated breakpoint of the forewing CV1 in 1960 (se = 0.0058, Pr(>F) = 6.42×10^-8) and for the hindwing in 1960 (se = 0.0051, Pr(>F) = 0.00511).

Wing shape dynamics pre/post 1960

Due to the results of the Mann-Whitney’s pairwise comparison (Appendix 1) and the CVA, observing a significant change in the wing morphology after 1960 (1961) we pooled the dataset into two groups (1938 – 1958 with 307 individuals and 1960 – 1968 with 199 individuals) to mitigate the possibility of a bias due to low numbers of individuals collected in some of the years. The scatter plots (Figs 9, 10) display the same pattern of statistically significant (<0.0001; P = 0.05) changes in morphology after being pooled.

Figure 9.

Scatter plot of the P. apollo forewing shape change along the first two canonical variates as axes. CV1 and CV2 cumulatively displaying 81.12% of the variability with 80% confidence ellipses.

Figure 10.

Scatter plot of the P. apollo hindwing shape change along the first two canonical variates as axes. CV1 and CV2 cumulatively displaying 76.39% of the variability with 80% confidence ellipses.

The scatter plots of CV1 and CV2 are displaying a shift from a narrower to broader forewing, with more extremes in either extra broad or narrow forewings in individuals from the 1960 – 1968 group (Fig. 9). The apex of the forewing is getting narrower and the landmarks M2 (4) and Cu2 (7) are moving closer along the time series.

For the hindwing (Fig. 10) the CV2 refers to a shift of the wing from broad to narrow along the Y axis. There is no general direction of hindwing landmark dynamics for CV1 and the movement on the hindwing is strongest at the M3 (14) and Cu1 (15), with an inwards landmark movement along the time series.

Extreme climatic events

Using historical weather data (temperature and precipitation) and photographic documentation of ECEs, we identified two extreme weather events (ECE) and compared them with the morphometric data (symmetric and asymmetric components, centroid size and Procrustes FA scores). The first ECE (a deep drop in monthly temperature during February 1956: -11.5 °C) corresponds with no deviation from the average wing size despite the extremely cold weather, as the wing size parameter lies within the range of the dataset (Figs 11, 12).

Figures 11, 12.

The correlation of the Centroid size standard deviation of the P. apollo fore and hindwing to the average temperature in February (the extreme weather event highlighted with a green circle).

Similarly, the extreme weather event, in this case heavy rainfall in June 1958, does not show any correlation to wing morphology (R² < 0.02), when comparing the daily, monthly, or annual precipitation total sum, average, min or max values at the study site. A weak positive correlation (R² = 0.12) was found for standard deviation of monthly precipitation to the wing size, when the wings grew larger with increasing variation of monthly precipitation. (Figs 13, 14).

Figures 13, 14.

The correlation of monthly sd of precipitation to Centroid size of the P. apollo fore and hindwing (the extreme weather event highlighted with a green circle).

Environmental impact

Since the changes in wing morphology show only weak correlations to temperature or precipitation, there is a strong suggestion of high anthropic effect. A large-scale road development project in November 1959 with heavy modification or complete destruction of some of the cliff faces was conducted in the immediate proximity of the collection site when constructing a new concrete road at Strečno castle, the Domašín meander and Žilina (Fig. 15A–C). The construction works were underway until at least 1965 and in 1964 heavy traffic was reported for the location.

Figure 15A–C.

Extensive construction works at the collection site, November 1959 (Kocián, 1959a, b, c).

The Mann-Whitney’s pairwise comparison and the Canonical variate analysis (CVA) revealed a significant change in wing morphology starting between 1958 and 1960, continuing to the end of the dataset. This is supported by the analysis of variance and the breakpoint regression analysis. The fluctuating asymmetry was statistically significant in all tested years, however with no significant changes over the period of 30 years, proving previous results and the hypothesis that P. apollo populations are able to maintain low long-term effective populations despite the high but constant amounts of asymmetry (Habel et al. 2009, 2012).

In principle, any external factor can produce plastic responses in organisms with the most prominent factors being temperature, photoperiod and humidity, where a single population can display annual or even seasonal variation (Blanckenhorn 2009). However, the analysis of the interaction of temperature and precipitation with morphology showed only weak or no distinct correlations. The first ECE (a deep drop in monthly temperature during February 1956) corresponds with no deviation from the average wing size despite the extremely cold weather, as the wing size parameter lies within the range of the dataset. February may be of high importance in the ontogeny of P. apollo as the eggs usually hatch at the end of February and the beginning of March (Žltková and Havranová 2017). Lower temperatures generally show positive effects on the overwintering eggs, although the overwintering larvae may be threatened by sudden temperature increases (Radchuk et al. 2013). Interestingly, Yu et al. (2012) found higher winter temperatures, especially in February, to have the highest negative correlation coefficient for the numbers of collected mountain Apollo specimens. Although our data does not show such a relationship, the wing size does correlate with the average monthly temperatures in February.

The second ECE with high short-term precipitation had little to no impact on the mountain Apollo populations at the study site; however we hypothesize that the abnormal rainfall may have had a negative impact on adults, pupas as well as host plants (McDermott Long et al. 2017). Although we do not reject the possibility of the ECEs’ impact on wing morphology, the available meteorological data suggest only weak or no correlations to changes in pre/post 1960 datasets, therefore a different effect had to be considered when evaluating impact on wing morphology.

A road development in a complicated terrain of an incised meander of the Strečno mountain pass was carried out from 1959 to 1965. The populations of mountain Apollo were present along the cliff faces and the surrounding meadows, creating ideal conditions with all the essential properties present for this species. We hypothesize that the removal of large parts of cliff faces which were used for the leveling of the road and the use of the open meadows for storage of material or by movement of the construction machinery changed the habitat and reduced the number of habitable patches that the population can occupy and limited the resources at the location (Blake et al. 1994; Brommer and Fred 1999; Fred et al. 2006). If we consider that even positive attempts to restore the habitats can change the morphology of butterflies by indirectly altering the host plant presence or quality (Sivakoff et al. 2016), the road work must be considered a major factor. Direct and indirect alteration or complete destruction of the habitats during the construction phase had most likely an impact on the morphology of the mountain Apollo butterflies, although the complex and often subtle morphological changes do not allow us to hypothesize further regarding the possible outcomes to the foraging behavior or overall fitness of the affected populations.

A more subtle, but lasting effect of the newly built road was a reported increase in the numbers of vehicles at the location. Since the morphological changes of the wings were observable until the end of our dataset in 1968 we here hypothesize that the newly built road could have continually affected the mountain Apollo and their host and nectar plants (Bengtsson et al. 1989; Trombulak and Frissell 2000; Muñoz et al. 2015), until the butterfly completely disappeared from this location in the 1990s.

Due to its conservation status and the vulnerability of current populations of mountain Apollos, the historic collections often provide the only and possibly the last opportunity to analyze large numbers of these individuals. The populations collected by Ján Zelný at Strečno mountain pass were gathered in such numbers due to the attempts of discovering a new subspecies based on the wing color patterns, where the method of collection was to gather the largest number of individuals during each visit. In most cases the selective nature of the amateur collectors could include a bias by selecting the largest or subjectively most beautiful specimens. Due to the nonselective nature of the collection method used by Zelný and the fact that he was the only person collecting and mounting the specimens we concluded that even if he had a specific preference in collection of butterflies the bias is the same in all collected populations. The statistical analysis of environmental effects on the wing morphology of the 506 individuals of P. apollo from a single location, collected periodically almost at the same time in June and July, over 30 years creates a unique look into a historical population of Parnassius apollo no longer present at the location.

The asymmetry of wings did not significantly change over time, nor could it be correlated to the analyzed environmental factors. Only weak or no statistical correlation with the meteorological data within the analyzed timeseries was detected. There is a strong suggestion of anthropogenic impact due to road construction on the changes in wing morphology at the studied site.

We would like to extend our gratitude to the staff of the Slovak National Museum in Bratislava – Natural Science Museum and the archives of the Slovak Hydrometeorological Institute for providing the data and invaluable help during our research.We are also grateful to Ivana Cocherová, Júlia Haruštiaková, Vladimír Janský, Ľubomír Rajter, and Michal Vlačiha for their help during the digitization of the samples. The research was funded by VEGA 1/0286/20.