Parasitism is a key factor in the population dynamics of many herbivorous insects, although its impact on host populations varies widely, for instance, along latitudinal and altitudinal gradients. Understanding the sources of geographical variation in host-parasitoid interactions is crucial for reliably predicting the future success of the interacting species under a context of global change.

Here, we examine larval parasitism in the butterfly Aglais urticae in south-west Europe, where it is a mountain specialist. Larval nests were sampled over two years along altitudinal gradients in three Iberian mountain ranges, including the Sierra Nevada, home to its southernmost European population. Additional data on nettle condition and adult butterflies were obtained in the study areas.

These data sources were used to investigate whether or not differences in parasitism rates are related to the geographical position and phenology of the host, and to the availability of the host plants.

Phenological differences in the host populations between regions were related to the severity of summer drought and the corresponding differences in host plant availability. At the trailing-edge of its distribution, the butterfly’s breeding season was restricted to the end of winter and spring, while in its northern Iberian range the season was prolonged until mid-summer. Although parasitism was an important source of mortality in all regions, parasitism rates and parasitoid richness were highest in the north and lowest in the south. Moreover, within a region, there was a notable increase in parasitism rates over time, which probably led to selection against an additional late-summer host generation in northern regions. Conversely, the shorter breeding season in Sierra Nevada resulted in a loss of synchrony between the host and one important late-season parasitoid, Sturmia bella, which may partly explain the high density of this butterfly species at the trailing-edge of its range.

Our results support the key role of host phenology in accounting for differences in parasitism rates between populations. They also provide insights into how climate through host plant availability affects host phenology and, ultimately, the impact of parasitism on host populations.

Study system

We studied the complex of larval parasitoids of Aglais urticae in three regions encompassing the whole of its latitudinal range in the Iberian Peninsula: the Pyrenees in the north of the Iberian Peninsula, the Sierra de Guadarrama in central Spain, and Sierra Nevada in southern Spain. Sampling sites were established along an altitudinal gradient in each region that covered most of the altitudinal range in which this species breeds. Nine sites were sampled for parasitoids in the Pyrenees at 1,127–2,560 m a.s.l., seven sites in the Sierra de Guadarrama at 1,150–2,004 m a.s.l., and six sites in Sierra Nevada at 975–2,532 m a.s.l. As part of a larger project aimed at investigating various aspects of the ecology of A. urticae, additional sites were surveyed in each region. The information gathered at these additional sites was used in this work to improve knowledge of the phenology of this butterfly.

Field sampling

To study the phenology of A. urticae adults, we used 500-m transects on which butterflies were counted every two weeks from March to September (a total of 15 sampling visits), following the standard methodology of the Pollard walks (Pollard & Yates, 1993). Butterflies were classified either as overwintered or freshly emerged based on wing colouration (i.e. dull or brightly coloured, respectively). Transects were walked at 15, 24 and 14 sites in 2016, and at 16, 24 and 20 sites in 2017, in the Pyrenees, Sierra de Guadarrama and Sierra Nevada, respectively.

To study the phenology of A. urticae larvae, we counted all larval nests found in four (exceptionally, just two and three at two sites) focal U. dioica patches in a subsample of the sites used for adult counts in each region. The focal patches were randomly selected along the butterfly transects and, if not available, at other accessible sites that were as close as possible to the transect route. The focal patches were visited every two weeks from March to September, whenever possible during the same visits as for the adult transect counts.

To study larval parasitism, larval nests detected at focal U. dioica patches (see above) were marked and, if larvae were in the third or later instars, 20 individuals were collected to assess parasitism. Because the total number of larvae per nest was sometimes less than 20, the overall average number of larvae per sample (± SD) was 16.2 ± 6.9. Moreover, given that the opportunistic parasitoid, Cotesia vestalis, is known to parasitize first instar larvae of the small tortoiseshell and to emerge mainly from the second instar (Audusseau et al., 2021), in 2017 we also collected 8 samples of five second instar larvae in all three regions (3 samples in the Pyrenees, 2 in the Sierra de Guadarrama and 3 in Sierra Nevada). We did not assess pupal parasitism, even if it may be important (Pyörnilä,1977; Shaw et al., 2009) because pupae are difficult to locate in the field, thereby precluding any reliable estimates of mortality.

Larvae were reared indoors in transparent plastic containers (155 x 105 x 45 mm) in groups of up to five individuals, all from the same sample. To avoid possible contamination, larvae were always reared with nettle leaves collected from their original nettle patch; if not available, nettles were harvested from sites where A. urticae and its closest congener, A. io, were absent, since some common parasitoids (e.g. the tachinids Sturmia bella and Pales pavida) lay microtype eggs on nettle leaves that can infect caterpillars if they eat these leaves. When a larva or pupa (in the case of larva-pupal parasitoids) produced a parasitoid, we recorded the stage at which the host was killed and kept the parasitoid individually in a vial until the adult emerged. Adults were preserved in pure ethanol for identification (Ichneumonoidea by MRS, Tachinidae by DH). Although hatching success was generally poor for most tachinids, careful inspection of puparia allowed for correct identification in almost all cases.

To examine host plant availability, we recorded the growing condition (i.e. quality level and height) of nettles over the season at a subsample of sites used for larval nest counts in each region. At each visit, two stems were randomly selected from each nettle patch. Their height was measured (in cm) and they were given a categorical value from 1 (worst quality) to 4 (best quality) in which (1) corresponds to already dry or withered plants, with senescent leaves; (2) to flowering plants and plants with green but not fresh leaves; (3) to old plants in which regrowth leaves were beginning to become visible (a common situation at the end of summer after rain or after herbivory); and (4) to vigorous plants, with fresh leaves. Category ranking was based on previous work showing how nettles in these various phenological stages differently affect larval growth rates, pupal and adult weights in the small tortoiseshell (Pullin, 1987) and the map butterfly, Araschnia levana (Mevi-Schütz & Erhardt, 2005).

Analyses

Host phenology

A combination of the standardized adult and larval count data was used to define the phenology of the species. GAM models were fitted to the adult (overwintering and fresh butterflies separately) and larval nest counts, which allowed us to extract the Julian day corresponding to each peak of abundance in a given region and season. GAM models were built using the package mgcv in R (Wood, 2011). In these models we used pooled data from 2016 and 2017 to increase the sample size and to improve the overall phenological picture in each of the study regions.

To investigate the potential altitudinal delay in larval phenology, we regressed separately the timing of larval nest appearance against site altitude for each year and region. The timing of larval nest appearance was summarised as the weighted mean appearance date.

Host plant phenology

We tested for differences in the phenology of nettles between regions using GAMM models, in which either nettle condition or height were the response variables and altitude, region, year, visit number (i.e. the timing of the season, used as the smoothing term) and the interactions of region with both visit number and altitude were the predictors. In these models, each individual stem was used as a data point and 'nettle patch' was entered as a random factor.

Impact of parasitism on host populations

To test whether or not the number of parasitoid species was comparable between regions (because sample sizes differed greatly between regions, see below), we computed the most common nonparametric estimators of species richness for each region and year separately (based on all recorded parasitoid species and genera) using the SpadeR package in R (Chao & Chiu, 2016).

To assess the impact of parasitoids on host populations, we calculated the parasitism rate for each larval nest as the number of larvae killed as a result of parasitism, after discounting those that died for unknown reasons (i.e. our calculations were always based on effective larval samples). To avoid biases resulting from low sample sizes, the parasitism rate was calculated for effective samples of ≥ 5 larvae. We obtained very similar results (not shown) when models of parasitism rate were built following a more restrictive criterion of effective samples of ≥ 10 larva.

The parasitism rate was modelled with generalized linear models (GLMs) using a binomial distribution and logit link function, with region, Julian day (date of nest collection) and altitude as predictors. However, because altitude and Julian day were highly correlated (r=0.76), only models with just one of these two variables were retained in the end. Models were built separately for 2016 and 2017 because the sampling sites in Andorra differed slightly between the two years. All possible models were built using the package lme4 in R (Bates et al., 2015); the best models were selected using the MuMIn package in R (Barton, 2015), with model selection being based on the Akaike Information Criterion (AIC). Models that differed by < 2 points from the lowest AIC (∆AIC < 2) were considered the top-ranked models (statistically equivalent to the best model of the set).

Phenological overlap between the host and its main parasitoids

The overlap (i.e. temporal co-occurrence) between the host and its two main parasitoids, Pelatachina tibialis and Sturmia bella, was estimated using the Overlap Parasitoid-Host index (OPH), as described by Audusseau et al. (2020), which is bouinbded from zero to one. This index was calculated for all possible combinations of site and year in each region. The maximum value of 1 is obtained in the hypothetical case when all individuals recorded in a given season, both of the parasitoid and the host, are concentrated in the same sampling event k. The opposite situation occurs (index value equaling to zero) when in all available samples one of the interacting species is always missing.

To understand which factors explain the degree of overlap between the parasitoid and the host, we used GLM models with a binomial distribution and a logit-link function. For each parasitoid species and year we applied a model in which OPH was the dependent variable and site altitude and region were the variable predictors.

All models were run in R (R Core Team, 2018).