Environmental stress in Parnassius apollo reflected through wing geometric morphometrics in a historical collection with a possible connection to habitat degradation

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 post1960 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.


Introduction
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.

Material and methods
All 506 Parnassius apollo specimens used for this analysis belong to the historical collection  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 col-lector 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.

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).
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(Klingenberg , 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, ac- cording 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.

Statistics
Statistics were conducted separately for the fore and hindwing to rule out the possible misalignment of wings during preparation (Bookstein 1997;Alibert and Auffray 2002). Using the MorphoJ software (Klingenberg 2008) a Procrustes superimposition was conducted following 3 steps for fitting the landmarks: (1) scaling the landmarks to the same centroid size (hereinafter "size"); (2) shifting the landmarks to the same position and (3) rotating them to fit the orientation as closely as possible minimizing the sum of squared 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. distances of landmarks (Klingenberg 2015). Procrustes ANOVA in MorphoJ was used to determine the morphological variance found among individuals when comparing the left and right wing. The same method was used to compare the left and right wing within each individual to test if there is a significant deviation from bilateral symmetry of our selected landmarks (Klingenberg and McIntyre 1998). The deviation to bilateral symmetry across the temporal variable is further referenced as FA (fluctuating asymetry). Procrustes residuals were used for canonical variate analysis (CVA) to test for differences between years using MorphoJ and PAST software (version 3.20) (Hammer et al. 2001). A scatterplot produced within the MorphoJ environment was created to display the comparison of the morphological changes. Boxplots of the FA, CV1 and CV2 scores were generated in the R software (version 3.4.4), using RStudio (version 1.1.456) (RStudio Team 2015). To compare the variance of morphology between years and between artificially created pre-and post-1960 datasets a factorial analysis was conducted using the "aov" function in R. Also, a breakpoint factorial analysis using "segmented" package was conducted to test the observed change on the morphology (RStudio Team 2015).

Results
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).

Fluctuating asymmetry (FA) throughout time
Rejecting the null hypothesis using Procrustes ANOVA (p < 0.0001) for the fore and hindwing corresponds with statistically significant levels in wing asymmetry within each year measured. As there was no significant variance in asymetry throughout the years for the forewing (ANOVA, p = 0.417) as well as the hindwing (ANOVA, p = 0.0564), the asymmetry is statistically significant within each year, although with no fluctuation when comparing the time series.

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.
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

6
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).

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  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. same pattern of statistically significant (<0.0001; P = 0.05) changes in morphology after being pooled. 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.

-1958 1960 -1968 hindwing
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.

Weather and wing morphology correlation
Generally, the comparative statistics of temperature, precipitation and morphology of individuals and their FA showed no distinct correlation (r 2 < 0. 04). The strongest correlation r 2 = 0. 25 (forewing) and r 2 = 0.30 (hindwing) of standard deviation of the centroid wing size to the average monthly temperature was observed during February. The other positive correlation refers to January temperature vs. wing size (r 2 = 0.21 for forewing, r 2 = 0.175 for hindwing) and March temperature vs. wing size (r 2 = 0.18 for forewing, r 2 = 0.16 for hindwing).

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).
Similarly, the extreme weather event, in this case heavy rainfall in June 1958, does not show any correlation to wing morphology (R 2 < 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 2 = 0.12) was found for standard deviation

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).

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).
of monthly precipitation to the wing size, when the wings grew larger with increasing variation of monthly precipitation. (Figs 13, 14).

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.

Discussion
The Mann-Whitney's pairwise comparison and the Canonical variate analysis (CVA) revealed a significant change in wing morphology starting between 1958 and 1960, Figure 15A-C. Extensive construction works at the collection site, November 1959 (Kocián, 1959a, b, c).

A B C
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(Habel et al. , 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 morpho-logical 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.

Conclusion
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.