Exploring Cold Hardiness within a Butterfly Clade: Supercooling Ability and Polyol Profiles in European Satyrinae

Simple Summary In insects distributed in temperate and cold zones, cold hardiness during overwintering crucially affects the distribution, including range shifts due to climate change. Our previous work on the genus Erebia, a cold-adapted and species-rich group of the sub-family Satyrinae (Nymphalidae), disclosed unexpected diversity of cold hardiness strategies, with closely related species surviving or not surviving freezing of larval body fluids. Asking whether this diversity is peculiar to this genus, or may be common in the Satyrinae clade, we investigated supercooling ability, contents of sugars and polyols in overwintering larvae tissues, and evolutionary signal of these traits of eight European Satyrinae species (from seven genera) and compared them with the Erebia representatives investigated earlier. We show that cold hardiness strategies are indeed diverse in the group and that high mountain and continental steppe species employ similar cryoprotection mechanisms, differing from those employed by species of more mesic environments. Abstract The cold hardiness of overwintering stages affects the distribution of temperate and cold-zone insects. Studies on Erebia, a species-rich cold-zone butterfly genus, detected unexpected diversity of cold hardiness traits. We expanded our investigation to eight Satyrinae species of seven genera. We assessed Autumn and Winter supercooling points (SCPs) and concentrations of putatively cryoprotective sugars and polyols via gas chromatography–mass spectrometry. Aphantopus hyperantus and Hipparchia semele survived freezing of body fluids; Coenonympha arcania, C. gardetta, and Melanargia galathea died prior to freezing; Maniola jurtina, Chazara briseis, and Minois dryas displayed a mixed response. SCP varied from −22 to −9 °C among species. Total sugar and polyol concentrations (TSPC) varied sixfold (2 to 12 μg × mg−1) and eightfold including the Erebia spp. results. SCP and TSPC did not correlate. Alpine Erebia spp. contained high trehalose, threitol, and erythritol; C. briseis and C. gardetta contained high ribitol and trehalose; lowland species contained high saccharose, maltose, fructose, and sorbitol. SCP, TSPC, and glycerol concentrations were affected by phylogeny. Species of mountains or steppes tend to be freeze-avoidant, overwinter as young larvae, and contain high concentrations of trehalose, while those of mesic environments tend to be freeze-tolerant, overwinter as later instars, and rely on compounds such as maltose, saccharose, and fructose.


Introduction
The recent debates on the effects of warming climate on biotic communities raise interest in the cold-adapted insects inhabiting narrow cold-climate zones, such as those in high mountains [1][2][3]. Due to warming climate, cold-adapted species may be impaired 1.
What are the supercooling abilities and concentrations of putatively cryoprotectant compounds in the studied species, do these measures of cold hardiness correlate with each other, and where do the studied species stand relative to the high-elevation Erebia spp.? 2.
How do the above characteristics of cold hardiness change from late autumn to winter? 3.
Which of the sugars and polyols, based on correlations between their concentration and cold hardiness measures, play a cryoprotective role in the studied species? 4.
Is there a relationship among the Satyrinae cold hardiness characteristics and their vertical distribution (the latter standing for the harshness of conditions experienced by overwintering larvae)? 5.
Is there a phylogenetic signal of the identity of cryoprotectant compounds used?

Study Species
The eight species we investigated (Table 1) all develop on grasses and typically form a single generation per annum. Table 1. Overview of the Satyrinae butterflies analyzed, including Erebia butterflies from earlier studies [14,15]. Localities and dates of sampling, their respective elevation, and numbers (n) of larvae used to assess supercooling point (SCP), lower lethal temperature (LLt), and total sugar and polyol concentration (TSPC). * The two numbers separated by "/" denote the numbers of larvae used for Autumn/Winter treatment.

Species (Phenomena Studied)
Origin (CZ-Czechia, AT-Austria) Elevation n (SCP) * n (LLt) n (TSCP) * Aphantopus hyperantus (Linnaeus, 1758): A Palearctic species, distributed from the Pyrenean peninsula to temperate China and Korea. In Europe, it reaches South Scandinavia and vertically occurs from the lowlands to the mountain zone (Czech Republic: <1200 m a.s.l.). Its larvae retain feeding activity during mild winters.
Chazara briseis (Linnaeus, 1764): Ponto-Mediterranean species, distributed from NW Africa to Central Asia, reaching Eastern Germany to the North. North of the Alps, its habitat is open-turf dry grasslands. The current distribution is severely fragmented [28,29]. In captive rearing, we observed larval feeding even during mild winters.
Coenonympha arcania (Linnaeus, 1761): European species, distributed from Western Europe to Asia Minor, the Caucasus and Ural Mts. Inhabitant of open woodlands, forest mantles and edges, requiring mosaics of grasslands and shrubs. A partial second generation appears in warm years [30].
Coenonympha gardetta (de Prunner, 1798): A sister taxon of the above [31]. It is a subalpine species restricted to the Alps, where it inhabits biotopes near the timberline.
Hipparchia semele (Linnaeus, 1758): European species with a prominently oceanic distribution [32]. Its range follows coastal areas from the Baltic countries through the British Isles to the Mediterranean region and Southern Russia. Populations in Central Europe are declining [33].
Maniola jurtina (Linnaeus, 1758): A west-Palearctic species, distributed from Northwestern Africa and the British Isles to Western Siberia and Iran. A generalist occurring from the lowlands to the mountain zone. Its European distribution reaches South Scandinavia latitudinally and ca 1000 m a.s.l. vertically. Widespread and common.
Melanargia galathea (Linnaeus, 1758): A west-Palearctic species, distributed from NW Africa to the Urals, Asia Minor and Transcaucasia; it is absent in more northerly areas (e.g., Scandinavia).
Minois dryas (Scopoli, 1763): A Euro-Siberian species, distributed from the Pyrenees across Europe (except southern peninsulas) and as far east as Korea and Japan. The distribution in Central Europe is discontinuous and includes both xeric and humid grasslands [34].

Captive Rearing
The pre-hibernating larvae of A. hyperantus, C. arcania, C. gardetta, M. jurtina, and M. galathea were obtained from wild-caught females. Larvae of C. briseis, H. semele, and M. dryas originated from an ex situ conservation rearing ( [29], Table 1). In both cases, the females oviposited in outdoor cages (wire frame 50 cm × 50 cm × 100 cm covered by nylon mesh) located in a half-shaded garden corner inČeské Budějovice, the Czech Republic (48 • 58 N, 14 • 28 E, 400 m a.s.l.) (the wild-caught species) and in a similar facility in Barchov, Czech Republic (50 • 12 N, 15 • 34 E, 250 m a.s.l.) (the captive-reared species). A maximum of five females were placed into one cage, the cages contained potted grass used as host plants at the localities. All the females readily accepted feeding by saccharose solution.
The larvae were taken from these rearing facilities immediately prior to experimental treatments, which were Autumn (November 2018) and Winter (January 2019). The Autumn treatment targeted the effects of the first frosts, the Winter treatment aimed at the highest cold hardiness level in the coldest month.
Identical rearing but different acclimation conditions were used in the earlier Erebia studies (cf. Table 1). Vrba, Konvicka, and Nedved [14] acclimated the larvae to a constant 5 • C, and Vrba et al. [15] used treatments simulating early autumn, late autumn, and winter. Here, we consider only late autumn and winter results (identical with Autumn and Winter here).

Supercooling Point and Polyol Profiles
The larvae were weighed and used for supercooling point measurement. SCP was measured individually using a PICO recorder with hand-made type K thermocouples [35], attached to the body of the experimental caterpillars enclosed in a syringe [36]. Larvae were gradually cooled above liquid nitrogen and the cooling rate was controlled at 1 • C/min. After an exotherm appeared on the recorder, the larva was kept in the cooling device until its body temperature decreased again to the crystallization temperature and then removed and warmed up at room temperature. Observation of movement of warmed-up caterpillars enabled us to determine if they were freeze-tolerant. Supplementary Table S1 contains the original data.

Statistical Analyses
Two-way ANOVA with species, treatment, and species x treatment interaction effects were used to compare SCPs, and the total concentrations of sugars and polyols (TSPCs) across the eight non-Erebia Satyrinae. An additional nine (SCPs) and five (TSPCs) Erebia species were assayed under somehow different acclimation conditions (Table 1); we compared only means and standard deviations of the earlier results.
Pearson's correlations were used to investigate relations between SCP, the elevation of origin (Table 1), TSPCs, and concentrations of individual sugars and polyols. The correlations were computed for both treatments, for Autumn and Winter separately, and after adding the Erebia results.
The composition of sugars and polyols were analyzed using multivariate statistics, the canonical correspondence analysis (CCA) in CANOCO [38]. CCA ordinates samples according to their composition and constrains the ordination according to predictors of interest. The numeric response variables were contents of individual compounds (logtransformed), whereas the factorial predictors were individual species. We ran the analyses with centering by species and samples and tested their significances using 999 Monte Carlo permutations.

Phylogeny of Sugar and Polyol Profiles
We used a phylogenetic tree of the studied species (Appendix A) to detect a possible phylogenetic signal in the sugar and polyol concentrations and the supercooling point of the larvae.
We computed, separately for Autumn and Winter, Blomberg's K [39] and Pagel's λ in R package "phytools" [40] for the major compounds (glycerol, fructose, glucose, sucrose, and trehalose), SCP and TSPC. These statistics compare the observed signal in a trait (a log-transformed in case of concentrations) to the signal under a Brownian motion model of trait evolution on a phylogeny. Blomberg's K is based on mean square errors, while Pagel's λ transforms phylogeny to fit the trait data. In both statistics, the values ≈ 0 correspond to a random or convergent evolution, while values ≈ 1 (both statistics) or >1 (Blomberg's K) indicate a phylogeny dependency. If a phylogenetic signal was present, we reconstructed the ancestral states on individual nodes of the tree using Felsenstein's phylogenetic independent contrasts (PICs) using the ace function in R "ape" package [41].

Sugar and Polyol Profiles versus Cold Hardiness
The fourth-corner analysis [42] relates three data tables, one with identity of the samples (here, butterfly species), one with the samples' properties (i.e., sugar and polyol concentrations), and one with species traits. The species traits potentially related to overwintering obtained from this study were: cold hardiness strategy (three states: freeze-avoidant, freeze-tolerant, mixed), SCP, TCSP (both continuous numeric: means across both treat-ments used for simplicity), overwintering larval instar (three states: 1-freshly hatched, 2-medium instars, 3-grown-up larvae before pupation), and elevation of sampling localities (continuous numeric; Table 1). In CANOCO, this analysis proceeded in four steps: (1) CCA constraining sugar and polyol composition~species|treatment; (2) the matrix of phylogenetic distances (from Appendix A) imported and subject to principal coordinate analysis (PcoA); (3) PcoA scores used to constrain CCA axes from the first step (control for phylogeny); (4) the resulting axes constrained by functional traits via redundancy analysis (RDA), a multivariate variant of linear regression. We used the forward-selection procedure to select the best fitting combination of traits.

Figure 1.
Overview of supercooling points, SCP (A), and total sugars and polyols concentrations, TSPC (B), found for the eight Satyrinae species assayed for this study, and the eight Erebia spp. assayed in [14,15]. Filled circles stand for Autumn, and empty diamonds for Winter treatments. The letters accompanying the marks denote significant differences among species and treatments, as revealed for the eight newly assayed Satyrinae species by Tukey's HSD test (Supplementary Table  S2). Color codes: freeze-avoidant species are blue, freeze-tolerant species red, mixed-strategy species violet. , found for the eight Satyrinae species assayed for this study, and the eight Erebia spp. assayed in [14,15]. Filled circles stand for Autumn, and empty diamonds for Winter treatments. The letters accompanying the marks denote significant differences among species and treatments, as revealed for the eight newly assayed Satyrinae species by Tukey's HSD test (Supplementary  Table S2). Color codes: freeze-avoidant species are blue, freeze-tolerant species red, mixed-strategy species violet.

Sugar and Polyol Profiles
The eight non-Erebia species differed in sugar and polyol profiles ( Table 2, Supplementary Table S4). In Autumn (Figure 3A), the alpine C. gardetta and steppe C. briseis contained high sorbitol, ribose, trehalose (mainly the former), and glucose (the latter). They contrasted from M. galathea with high maltose, fructose, glycerol, and saccharose. The remaining species were intermediate, with the steppe M. dryas containing high glucose, and the mesic grassland species C. arcania, A. hyperantus, and M. jurtina containing high ribitol, mannitol, and saccharose. After adding the earlier Erebia results ( Figure 3B), all Erebia spp. formed a loose group, associated with high trehalose (highest in the alpine E. cassioides), ribitol, threitol, and erythritol (the latter two unique for Erebia spp.). C. gardetta and C. briseis assumed intermediate positions, displaying association with sorbitol and glucose. The dry grassland species M. dryas and H. semele also inclined towards sorbitol and glucose, whereas the remaining non-Erebia mesic grassland species displayed high concentrations of saccharose, fructose, and myo-inositol. In Winter, glucose and myo-inositol increased in H. semele, glucose increased in M. galathea ( Figure 3C), and trehalose and glycerol increased in E. pronoe and E. aethiops ( Figure 3D).   Figure S5 for the treatments analyzed together.

Phylogenetic Signal and Ancestral State Reconstruction
The pruned tree topology corresponded to previous results on Satyrinae butterflies [12,43]. From a common ancestor, the two Coenonympha spp. branched off first, and the second division was between (Erebia + (A. hyperantus + M. jurtina)) and M. galathea + (H. The main patterns were retained for both treatments analyzed together, and the treatment effect was considered covariate (Table 2, Supplementary Figure S1). The mesic grassland/lowland species were still associated with high glucose, maltose, ribose, sorbitol, and fructose, whereas the alpine species contained high trehalose, threitol, erythritol, and ribitol. The secondary distinction between steppe (high glucose, ribose, and sorbitol) and mesic grassland (high saccharose, mannitol, and fructose) species was also retained. Treating species identities as covariates showed that maltose, saccharose, and ribitol were higher in Autumn, whereas arabinitol, erythritol, sorbitol, ribitol, and trehalose were higher in Winter.

Phylogenetic Signal and Ancestral State Reconstruction
The pruned tree topology corresponded to previous results on Satyrinae butterflies [12,43]. From a common ancestor, the two Coenonympha spp. branched off first, and the second division was between (Erebia + (A. hyperantus + M. jurtina)) and M. galathea + (H. semele + (M. dryas + C. briseis))). Blomberg's K values were significant, implying a phylogenetic dependency for Winter glycerol and Winter SCP, and marginally significant for winter TSPC (Table 3). Pagel's λ was significant only for winter glycerol, but the value was also close to 1 for Winter SCP. Table 3. Results of Blomberg's K and Pagel's λ testing the phylogenetic signal in major cryoprotectants, supercooling points (SCP) and total sugar and polyol concentrations (TSPC). See Appendix A for description of inference of the phylogenetic tree, and Figure 4 for mapping the traits onto the tree.

Discussion
Expanding the interest of overwintering butterfly larvae cold hardiness from the cold-adapted genus Erebia to a broader sample of European univoltine Satyrinae reveals that the high diversity of strategies and mechanisms detected earlier [14,15] is not restricted to the single genus. Among the eight non-Erebia species, two survived freezing of their body fluids, four were killed by the freezing, and three displayed a mixed strategy. The supercooling point values varied among all species by >10 °C, and concentrations of sugars and polyols differed fourfold (sixfold in Autumn) and eightfold when considering the previously assayed Erebia spp. There was neither a straightforward relationship between per-species SCP and TSPC, nor a clear correlation between concentrations of individual compounds and SCP. Instead, we detected phylogenetic signals in Winter glycerol, Winter SCP, and Winter TSPC. High-elevation species with high TSPC tend to be freeze-

Discussion
Expanding the interest of overwintering butterfly larvae cold hardiness from the coldadapted genus Erebia to a broader sample of European univoltine Satyrinae reveals that the high diversity of strategies and mechanisms detected earlier [14,15] is not restricted to the single genus. Among the eight non-Erebia species, two survived freezing of their body fluids, four were killed by the freezing, and three displayed a mixed strategy. The supercooling point values varied among all species by >10 • C, and concentrations of sugars and polyols differed fourfold (sixfold in Autumn) and eightfold when considering the previously assayed Erebia spp. There was neither a straightforward relationship between per-species SCP and TSPC, nor a clear correlation between concentrations of individual compounds and SCP. Instead, we detected phylogenetic signals in Winter glycerol, Winter SCP, and Winter TSPC. High-elevation species with high TSPC tend to be freeze-avoidant and overwinter in earlier larval instars. During acclimation, they accumulate trehalose and glycerol, as well as threitol and arabinitol, the latter two biosynthesized from glucose via the pentose phosphate pathway [44]. Low-elevation species, on the other hand, tend to overwinter in later instars and with elevated levels of glucose, saccharose, fructose, and maltose, i.e., compounds directly involved in primary metabolism.
The difference between high-and low-altitude species is unlikely due to larval food, as all the species develop on closely related plants (which can be interchangeable in captive rearing [10]). The overwintering in early larval instars in high-altitude species (and the steppe species C. briseis) is probably due time constraint for pre-diapause development, due to short season in the mountains [7] and late-season reproduction in C. briseis [28].
A sample of eight species (sixteen with the earlier Erebia results) assayed may seem high, given that studies of insect cold hardiness rarely sample multiple representatives of well-defined clades (but see [45][46][47][48][49][50]). It is admittedly poor compared to the radiation of temperate Satyrinae [12,51,52], but still allows cursory inference regarding ecological and geographical correlates.
Comparisons of insect antifreeze strategies across distant [53] and related [45,54] species indicate that ability to survive the freezing of body prevails in regions with oceanic climates (subantarctic islands: [45]; coastal mountains: [20]), which experience short, unpredictable, but frequent frosts. Continental conditions with temporarily predictable and undisrupted freezing periods favor freeze avoidance. In our sample, the connection between oceanicity and freeze tolerance applies for H. semele [32], whose coastal populations are prospering at present, whereas inland populations are declining [33]; however, this is also true for A. hyperantus, a drought-sensitive species [55,56]. Erebia medusa, also displaying this strategy, is peculiar among European congenerics by a broad distribution spanning from lowlands to mountains and from damp to xeric habitats [57]. A mixed strategy, previously reported for E. aethiops [15], which shares wide vertical distribution with E. medusa [58], was detected here in three additional species. Minois dryas inhabits a wide Euro-Siberian range with various climatic conditions; moreover, in Central Europe it inhabits both dry and humid habitats [34]. Chazara briseis is a specialist of continental steppes [28,29], whereas M. jurtina is the most common European grassland butterfly, persisting even in intensively exploited landscapes [59,60]. Other lepidopteran examples displaying mixed responses to freezing of body fluids are Papilio zelicaon (Lucas, 1858) (Lep.: Papilionidae), a North American species distributed from the humid coast to arid inland [16], and Pieris rapae (Linnaeus, 1758) (Lep.: Pieridae), a multivoltine generalist distributed from subtropical to arctic regions [61]. The mixed strategy thus appears beneficial in variable conditions, climatic or otherwise, among seasons or even within an individual lifespan. Climate variation among seasons is certainly the case of continental steppes, with prominent inter-seasonal variation in snow cover and temperature [62]; recall that C. briseis larvae may feed during mild winters. From all species assayed, C. briseis and M. jurtina also displayed the largest variation in SCP between seasonal treatments, indicating flexible reactions to autumn conditions. E. aethiops displays local adaptations in adult thermoregulation [63]; other adaptations, variable among individuals, might exist in larval cold hardiness.
The conflicting information implies that individual compounds fulfill different roles among taxa and during phases of overwintering [78]. What matters eco-evolutionarily is the functional outcome, expressed here as the value of SCP. Our results suggest that the SCP reflects a phylogenetic signal (or convergent evolution) towards low values in alpine Erebia spp., and in species of continental steppe environments. A special case was the two Coenonympha representatives, the subalpine C. gardetta and lowland C. arcania, both displaying low SCP but also low concentrations of putative cryoprotectants. Possibly, larvae of this genus, phylogenetically most distant to the remaining species [31], employ other cryoprotective agents, such as amino acids and antifreeze proteins. The case of the subtropical orchard pest A. ceratoniae indicates that a species experiencing rapid and unpredictable changes of winter weather can be supremely protected against freezing [68]. The low SCP in C. arcania could be a pre-adaptation of its ancestor for colonization of high elevations by its descendants (cf. [31]).
The evolutionary signals in glycerol concentration, TSPC, and SCP were apparent only in the Winter treatment. The conditions experienced by pre-hibernation larvae likely vary with species and habitats, and diverse external signals may launch cryoprotection [71]. Phasing of diapause is crucial for cold hardiness development [79]. Constitutive (Autumn) cryoprotectants can be replaced by inducible (Winter) ones [77]. Although we found some Autumn vs. Winter changes in the cryoprotectant profiles, there was no clear pattern allowing to distinguish between constitutive and inducible compounds, or between the role of polyols vs. sugars.
High trehalose and several polyols occurred in alpine species (Erebia spp., C. gardetta), in which they correlated with elevation, but also in the steppe C. briseis. Across biotic realms, trehalose enhances cryoprotection in extreme environments [80,81]. However, its high content did not decrease SCP among our species. It more likely protects cells against desiccation [82], which may be critical both at dry steppes (C. briseis) and wind-exposed alpine cliffs (E. pluto, E. cassioides) (cf. [83]). On the other hand, mono-and disaccharides were more typical for low-elevation species. Among them, inhabitants of dry grasslands (M. dryas, C. briseis, H. semele) contained high concentrations of sorbitol, glucose, and ribose, whereas mesic habitat species (M. galathea, M. jurtina, C. arcania, A. hyperantus) contained high saccharose, maltose, and mannitol.
The biochemical pathways synthesizing and assimilating individual polyols are closely related [84]. Individuals modify both cold hardiness and polyol profiles during the season, and variation among populations exists [67]. Given the diversity of cold hardiness strategies among temperate Satyrinae, contrasting with their uniform feeding on grasses and developmental modes (overwintering as larvae), it is tempting to speculate whether adaptation to varying climatic niches could have propelled the evolution of their remarkable species and habitat use diversity.    Institutional Review Board Statement: Not applicable. The study did not require ethical approval as it was not involving humans or other vertebrates.

Informed Consent Statement: Not applicable.
Data Availability Statement: The data on SPC and survival of overwintering larvae are provided in Supplementary Table S1, and the mean (±SD) concentrations of sugars and polyols in Supplementary  Table S4. More detailed data can be provided on request by the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Reconstruction of Phylogenetic Tree Used for Phylogeny of Sugar and Polyol Profiles in Satyrinae Butterflies
We used database sequences to reconstruct a Bayesian tree, which was further used for detection of a possible phylogenetic signal in cryoprotectant concentrations and supercooling points of Satyrinae butterfly larvae.
We downloaded sequences for 46 individuals from GenBank (Table A1, https://www. ncbi.nlm.nih.gov/genbank/ accessed on 10 September 2019), covering the 13 studied species in several specimens, related species to avoid long-branch attractions, and related Nymphalid butterflies as outgroups. We used a set of five genes available for this group (COI, EF-1 alpha, GAPDH, RpS5, and wingless; EF1alpha was mostly missing for the genus Erebia), resulting in a concatenated alignment of 4112 bp. A Bayesian ultrametric tree was constructed in Beast 1.8.2 [85], with the partition scheme and substitution models estimated by PartitionFinder2 [86], under speciation: birth-death incomplete sampling method, with 50,000,000 generations and three independent runs. The convergence of the runs was checked in Tracer v.1.7 [87]. Results were summarized onto a single target tree. The final tree was pruned to contain only a single specimen of each of the 13 target species. The PartitionFinder2 divided the concatenated alignment into three subsets and selected the GTR + I + G model for the first subset and GTR + G model for the second and third subsets. The tree topologies ( Figures A1 and A2) resembled the previous analyses from which the majority of the data used here originated [12,43].