Low winter precipitation, but not warm autumns and springs, threatens mountain butterflies in middle-high mountains

Study area and species

Both the Jeseník Mts (highest summit: Praděd, 1,491 m, N 50°4.96′, E 17°13.85′, hereinafter “J”) and Krkonoše Mts (Sněžka, 1,603 m, N 50°44.16′, E 15°44.38′ E, hereinafter “K”) were formed by Variscan orogeny and consist of rolling ridges built from metamorphic (J, K) and crystallinite (K) rocks. The timberline, situated at ca 1,300 m, is formed by sparse Picea abies growths, replaced upslope by Pinus mugo shrubs (the latter non-native in J: Kasak et al., 2015). The summit elevations are covered by species-poor grasslands and heaths, dominated by Nardus stricta, Avenella flexuosa, Vaccinium myrtillus, and Calluna vulgaris (≈10 km² in J, ≈50 km² in K). Much richer and structurally diverse vegetation is found at valley headwalls near the timberline (Bures, 2018, Fig. 1). The (sub)alpine fauna of the mountains is impoverished due to isolation from other massifs overtopping the timberline (aerial distance J–K: 130 km; J–Mala Fatra Mts, Carpathians: 160 km; J–Schneeberg Mt, Eastern Alps: 270 km).

Erebia sudetica occurs in several mountain systems outside of the Alps: Mons du Cantal, France; the Carpathians (Rodna Mts, Retezat Mts, Godeanu-Tarcu Mts, Ciuicas Mts, Fagaras Mts, Cuvelier & Dinca, 2007); Jeseníky Mts, the Czech Republic. It occurs in two Western Alps districts (E. Isére, France; the Bernese Alps, Switzerland). In terms of taxonomy, it appears paraphyletic with E. melampus (Fuessly, 1775) (Haubrich & Schmitt, 2007; Pena et al., 2015) and J represents the northernmost and highly isolated locality of the melampus/sudetica species complex. In J, it forms spatially restricted colonies at subalpine tall-herb formations, plus at lower-elevated woodland openings (Konvicka et al., 2014).

Erebia epiphron is the most widely distributed European alpine Erebia. Besides J and K, it inhabits the mountains of Scotland and northern England, Cantabrian Mts in Spain, the Harz Mts in Germany (extinct), the Pyrenees, Massif Central and Vosges in France, Central Apennines in Italy, Alps, Carpathians, and the Dinaric mountains of the Balkans southward to northernmost Greece (Hinojosa et al., 2018, see Minter et al., 2020). The populations of particular mountain regions are genetically unique (Minter et al., 2020). Across this range, it prefers Nardus dominated grasslands on nutrient-poor bedrock (Ewing et al., 2020). It is restricted to elevations above the timberline in J (Kuras et al., 2003), whereas its non-native K population has also colonised cultural grasslands down to ≈1,100 m (Cizek et al., 2015).

Both species are in flight in July–August. During warm years, the adults emerge at the end of June. Larvae are oligophagous on grasses (Sonderegger, 2005). They feed during September and October, overwinter in grass tussock, and resume feeding after the snow melt in mid-May. Biennial development occurs in E. epiphron in the Alps (Sonderegger, 2005), whereas univoltine development was observed in E. melampus (Wipking & Mengelkoch, 1994), closely related to E. sudetica. In our study system, we assume that univoltine development prevails also in E. epiphron, as the species does not display biennial adult abundance fluctuations in J (Konvicka et al., 2016) and both species developed with a single overwintering in a rearing experiment (Kuras et al., 2001).

Adult numbers monitoring

The study is based on transect monitoring of relative abundance of adult butterflies (2009–2020, plus an earlier 5-year period 1995–1999 for one of the transects), combined with meteorological data from nearby weather stations. This time span is sufficient to detect butterfly population trends and to test the effects of climate variability on adult abundance changes (see Wilson & Fox, 2020). We established three transects in J (J1–J3, total length 4.2 km, altitude 1,325–1,492 m a.s.l.) and four in K (K1–K4, total length 7.1 km, altitude 1,120–1,510 m a.s.l.) (Fig. 1 and Supporting Information, Supplementary methods).

We monitored both species in J in 2009–2020 (12 years), plus E. epiphron, J1 only, in 1995–1999 (5 years). In K, E. epiphron was monitored in 2010–2020 (11 years) (Table S1). The transects were walked 2–4 times per week, weather permitting, between 1 July and 15 August, except for 2009, when the monitoring terminated earlier, and the extraordinarily warm year 2018, with the first random records from June 7. Monitoring started earlier (15 June) during 2019 and 2020 to cover the possibility of an earlier emergence of adults. However, we did not record any adults before the beginning of July during these 2 years with extended monitoring period. The walks were restricted to air temperature > 15 °C, clear or half-overcast sky, and wind < 11 km/h. Permissions for the field research were obtained from the administration of the Jeseníky Protected Landscape Area (permit no. 442/JS/10) and the administration of the Krkonoše National Park (permit no. KRNAP 03599/2013).

Abundance and phenology indices

We described the relative annual abundance using the population abundance index (PAI) (Kleckova, Vrba & Konvicka, 2015; Konvicka et al., 2016), a measure derived from the shapes of the population recruitment curves, defined by daily relative abundances related to day-of-year (Rothery & Roy, 2001) (Table S1). We fitted the curves (Figs. S1–S11) via the generalised additive models (GAMs) using mgcv (Wood, 2011) in R 3.4.4 (R Development Core Team, 2018), with cubic splines (k = 4) and quasi-Poisson errors. The number of daily records were standardised for unit transect length. To facilitate the convergence of the fitted curve to zero prior to and after the adult flight period, zero values were added 1 week before the first monitoring record and 2 weeks after the last adults observation on the transect. The areas under the GAM curves were computed using the rectangle method, which approximates a definite integral. Specifically, the date interval between the added zeroes for each species and transect was divided into subintervals of length 0.01 days (x), predicted abundance value (y) was generated for each subinterval by the relevant GAM, and the areas of the resulting rectangles (x * y) were summed across the entire date interval to obtain the approximate area under the curve, i.e., the population abundance index (PAI).

For phenology, we used the GAM-fitted population recruitment curves (Figs. S1–S11) and the function predict (library mgcv) to infer the days when 20% and 80% of individuals were in flight (an example script is available in Figshare at https://doi.org/10.6084/m9.figshare.14642472). The 20% day is herein the flight period onset, and 20–80% day interval length is the flight period duration (Table S1). This diminishes the uncertainty which would arise from using the dates of the first and last observed individual to describe the flight period duration.

Climate variables

Monthly weather data were provided by the Czech Hydrometeorological Institute (https://www.chmi.cz/historicka-data/pocasi/mesicni-data). For J, 1995–1997 records were from the Mt Praděd meteorological station (1,491 m, N 50.0831° E 17.2310°) adjoining the transect J3, and 2009–2020 records from Mt Šerák station (1,328 m, N 50.1874° E 17.1082°), located 15 km NW, but climatically equivalent to the study sites. For K (2010–2020), we used records from the Luční Bouda station (1,410 m, N 50.7344° E 15.6973°) for the summit transects K3 and K4, and from the Pec pod Sněžkou station (820 m, N 50.6918° E 15.7287°) for the lower-altitude transects K1 and K2.

Adult abundance and phenology may be affected by weather experienced by the immatures (Turlure et al., 2010), including extremes such as ground frosts. For the pre-hibernation larvae (September–October), we considered total autumn precipitation (P_totAut), as well as monthly average (T_avgAut), maximum (T_maxAut), minimum (T_minAut), and ground minimum (T_groundAut) temperatures, averaged across the 2 months (Table S1). For the overwintering (November–April), when the mountains are usually snow-covered, we used the variables as above (indexed P_totWin, T_avgWin, T_maxWin, T_minWin, T_groundWin) and the snow cover duration (in days, Snow_days). For the spring larval period (May–June), we used the same predictors as for autumn (indexed P_totSpr, T_avgSpr, T_maxSpr, T_minSpr, T_groundSpr). The same applied for summer (July–August), when the weather directly affects the adults (P_totSum, T_avgSum, T_maxSum, T_minSum, T_groundSum).

Whereas most of the climate variables were available for all years in J (1995–1997 + 2009–2020, herein J_all) the monthly number of snow days (Snow_days) and ground minimum temperature (T_ground) were available only for the J recent period (2009–2020), J_rec. The ground minimum temperature, measured at 5 cm above the ground, provides a proxy for days with ground frost, possibly critical for larval survival (Vrba, Konvicka & Nedved, 2012). T_ground and Snow_days were not available for K.

We calculated the Pearson’s correlation coefficients for all climate variables for J_all, J_rec, K lower altitudes (transects K1, K2), and K higher altitudes (K3, K4), and visualised the matrices using the package corrplot (Dray, 2008) (Figs. S12–S15). We also computed the Principal Component Analysis (PCA) in vegan for R (Oksanen et al., 2019) (Fig. S16). These analyses revealed multiple correlations among the climate variables, but without a clear general pattern.

Statistical analyses

We conducted a separate analysis for each mountain range and for the two periods in J: J_all and J_rec. We ran the analyses separately for the two species in J_rec. We first tested whether there was a temporal trend in the abundance and phenology of the butterflies, and in the climate variables. We used generalised linear models (GLM) in R to test the linear effect of year on the response variables, with transect identity as a fixed effect. We also tested for possible biennial abundance fluctuations in population abundance index of E. sudetica and E. epiphron using a GLM with odd vs even year as a predictor, including the transect identity and linear effect of the year as covariates. We similarly tested whether PAI depended on Onset and Duration, and whether the two phenological variables were correlated.

For the effect of climate predictors on PAI, Onset and Duration, we tested the effects of conditions experienced by autumn, overwintering, and spring larvae, and the effect of conditions experienced by adults during their flight period on PAI and Duration. Separate analyses were performed for four datasets: E. sudetica, E. epiphron K, E. epiphron J_rec, and E. epiphron J_all to account for differences in the availability of climate variables between J and K, and between J_rec and J_all. Analysing E. epiphron data separately for J and K allowed detection of possible inter-population differences; separate analyses of data from J_rec and J_all also facilitated the comparison of trends in E. epiphron and in E. sudetica, because the latter was monitored only during the recent period (2009–2020).

We used an information theoretic approach, which is more appropriate for comparing large numbers of models with different predictors than a frequentist hypothesis testing (Symonds & Moussalli, 2011). Understanding that climate variables are always causally interdependent and numerically intercorrelated, and given the rather low degrees of freedom (3/4 transects for J/K $\times$ 11/12/15 monitoring years for K/J_rec/J_all), we did not fit models with multiple predictors and their interactions, as this would risk overfitting and unreliable estimates (Harrell, 2015). Instead, we evaluated individual predictors separately.

We used GLM with population abundance index (PAI, log-transformed) per year and transect as the response variable and transect as a fixed factor. The predictors were standardised to zero mean and unit variance, which enabled relative comparisons of their effects. The residuals were checked against the predicted values using simulateResiduals in DHARMa for R (Hartig, 2016).

We first constructed, for each of the four data sets, a null model not containing predictors: log(PAI)~1 + Transect. Its AIC was used to evaluate subsequently fitted models according to ΔAIC between the null and a fitted model(s). We considered ΔAIC < −6 a strong and < −2 a moderate indication of an effect of the respective predictor (Richards, Whittingham & Stephens, 2011). Because population abundances may reflect dependence on previous years, we tested for the presence of a temporal autocorrelation using the corAR1 function (autocorrelation structure of order 1). For all data sets, the null models without autocorrelation had smaller AICs than those with autocorrelation, and we did not consider autocorrelation further.

We used an analogous procedure for flight period Onset and Duration. We did not consider autocorrelation with the previous year, as there is no reason to expect dependency of flight timing on the previous year abundance.

Correlations among variables

We found a strong negative correlation between the flight period Onset and Duration in E. sudetica (r = −0.70, 95% confidence interval = [−0.84, −0.48]), E. epiphron, J (J_all: r = −0.73, [−0.84, −0.53]; J_rec: r = −0.71, [−0.84, −0.50]) and K (r = −0.64, [−0.79, −0.41]). Thus, the earlier onset resulted in longer duration of the flight period. In E. sudetica and E. epiphron Jeseniky Mts., the relative abundance (PAI) significantly increased with Duration but was not affected by Onset (Table 1, Fig. S17). Some of the climate variables were correlated, in particular different temperature variables within individual parts of the year (e.g., T_avgWin and T_minWin) (Supporting information, Figs. S12–S16).

Temporal trends

We detected several temporal trends in climate variables over the entire study period in J_all (i.e., including the 1990s data), but most of them were not apparent over the shorter recent period (J_rec and K). Specifically, most measures of temperature increased in J_all, winter precipitation increased in J_rec, and summer precipitation decreased in K. The trends in these variables were similar in both mountains, but often weak in terms of statistical significance (Table S2).

Despite large inter-annual variation, adult abundance significantly increased in E. sudetica J and E. epiphron K, but not in E. epiphron J (Table 2, Fig. 2). There was no temporal trend in Onset (Fig. S18) and in most cases in Duration (Fig. S19, Table 2) with the exception of E. epiphron J_all, for which the flight period Duration increased since the 1990s (Fig. S19). E. epiphron K displayed an evidence of biennial fluctuations, with population abundance index ≈ 1.8 times higher in odd years than in even years (GLM, the model including the odd/even year had ΔAIC = −4.7 compared to the null model), but no such pattern applied for either of the species in J.

Climate variables and abundance

In E. sudetica, warmer autumns (T_minAut, T_groundAut), higher winter precipitation (P_totWin) and higher spring temperatures (T_minSpr, T_maxSpr, T_avgSpr) increased PAI in the following year. In contrast, increasing spring precipitation (P_totSpr) and increasing minimum summer temperatures (T_minSum) decreased population abundance index in the following year (Table 3).

For E. epiphron J, adult population abundance index increased with autumn temperatures (T_minAut, T_groundAut), winter precipitation (P_totWin) and spring temperatures (T_avgSpr, T_maxSpr, T_minSpr), but decreased with winter temperatures in (T_avgWin) (Table 3) and higher summer minimum ground temperatures (T_groundSum). The negative effect of winter temperatures (T_avgWin, T_minWin, T_groundWin) was more distinct for J_rec than for J_all (Table 3). For J_rec, T_minWin and T_minSum displayed moderate effects, not detected for J_all.

For E. epiphron K, the predictors increasing subsequent population abundance index consistently with E. epiphron J were higher autumn temperatures (T_minAut), higher winter precipitation (P_totWin), and higher spring temperatures (T_avgSpr) (Fig. 3). The negative effect of warm winters, prominent in J, was not detected in K. Another difference was higher precipitation in summer (P_totSum), reducing adult abundances in J (ΔAIC < −2) but not in K. The major trends of warm springs increasing adult abundance, and winters with little precipitation decreasing it, were consistent for the two mountains (Table 3).

Climate variables and phenology

Adult phenology of E. sudetica (Table S3 and S4, Figs. 4 and 5) was affected by climate in all larval phases. Higher autumn ground temperatures (T_groundAut) advanced Onset, i.e., accelerated larval development. Higher spring precipitation (P_totSpr) postponed Onset and shortened Duration. Warm springs (T_avgSpr, T_maxSpr, T_minSpr) advanced Onset and prolonged Duration. Summer weather had no effect on Duration (Table S3).

In E. epiphron J, warmer autumns (T_minAut, T_groundAut) and springs (T_avgSpr) advanced Onset (Tables S3, Fig. 4) and prolonged Duration (Table S4). Higher spring precipitation (P_totSpr) had opposite effects. In summer, flight period Duration increased with T_avgSum.

In E. epiphron K, the factors affecting phenology partly differed from J, but there were commonalities (Table S3, Figs. 4, 5). As in J, rainy springs (P_totSpr) delayed Onset and shortened Duration, whereas warm springs (T_avgSpr, T_minSpr) advanced Onset. Contrary to J, warm autumns (T_maxAut) delayed and rainy autumns (P_totAut) advanced Onset (Table S3). In addition, warm winters (T_avgWin, T_minWin) and summers (T_maxSum) shortened Duration, which was not detected in J (Table S4).