Rewiring of interactions in a changing environment: nettle-feeding butterflies and their parasitoids

is elevated at sites where Ar. levana has been established for a longer period. In our study system, variations in butterfly species assemblages were associated in a predictable way with substantial variations in rates of parasitism. This relationship is likely to affect the dynamics of the butterfly host species, and potentially cascade to the larger number of species with which they interact. These results highlight the importance of indirect interactions and their potential to reorganise ecological communities, especially in the context of shifts in species distributions in a warmer world.


Introduction
Biotic interactions are key drivers shaping the composition and dynamics of ecological communities. Although these interactions occur locally, their impact on species assemblages can be detected across larger scales through their effects on population dynamics, species co-occurrence, community structures, distributions and abundances (Araújo and Luoto 2007, Meier et al. 2010, Wisz et al. 2013, Belmaker et al. 2015. Biotic interactions can be direct, such as predatorprey and insect-plant interactions, or indirect when the interactions between two species are mediated by other species, such as shared resources or enemies. Direct interactions are widely studied and documented, but the vast majority of the biotic interactions operate in complex networks where species are connected through both direct and indirect interactions Shaw 1974, 1986). Unfortunately, studies which examine biotic interactions and their effects on species assemblages are limited by the scarcity of empirical data, especially for indirect biotic interactions. Our understanding of indirect interactions is mainly derived from limited empirical data gathered under laboratory conditions or at relatively small spatial and temporal scales (Holt and Bonsall 2017). The difficulty of collecting detailed and community wide field data hampers our ability to analyse the stability of biotic interactions, how they vary in space and time, and how they impact populations and communities (Holt and Bonsall 2017).
Understanding the nature and outcome of biotic interactions is all the more important as they are likely to be altered by climate and land use changes, with further impacts on ecological communities (Tylianakis et al. 2008, Blois et al. 2013. Although the importance of biotic interactions is widely recognised in the literature on global environmental change, for example in the context of refining predictions of species' responses to environmental change (Wisz et al. 2013, Dormann et al. 2018, measuring how and understanding why these interactions alter as the environment alters remains a major challenge. This difficulty is partly due to their dynamic nature and the many ways that environmental change can affect them at different spatial scales (Wisz et al. 2013, Kissling and Schleuning 2015, Pellissier et al. 2017, Dormann et al. 2018. For example, changes in climate and land use can affect the distribution and demography of species, which in turn might alter the nature and strength of biotic interactions (Tylianakis et al. 2008, Early andKeith 2019). Likewise, differences in the sensitivity of interacting species to environmental change can alter their abundances and the spatial pattern in their co-occurrence, disrupt their temporal synchrony or lead to novel interactions in the case of invasive or range-expanding species (Blois et al. 2013).
Insects constitute more than half of the described biodiversity of Earth and underpin a wide range of ecosystem functions and services (Losey and Vaughan 2006). Identifying and understanding the forces that shape biotic interactions in this group, with consideration to changing environments, is therefore of high importance. Holt (1977) argued that indirect interactions mediated by shared natural enemies are strongly involved in structuring communities of herbivorous insects. These indirect interactions, where the population dynamics of species at the same trophic level can be linked via the action of shared natural enemies, is called apparent competition (Holt andLawton 1993, 1994). For example, apparent competition mediated by shared parasitoids was shown in leafhopper communities in California where the introduction of a new host species, combined with the strong preference for the native species, resulted in an overall increase in parasitoid pressure and decline of the native species (Settle and Wilson 1990). More generally, in natural systems, the impact of apparent competition can affect multiple species that share common enemies (Morris et al. 2004, Frost et al. 2016 and was shown to vary with the size of the community, the abundance of hosts and their phenology (Bonsall and Hassell 1997, Van Nouhuys and Hanski 2000, Morris et al. 2004, Blitzer and Welter 2011. Further studies also suggest that biotic interactions such as those between hosts and parasitoids will be influenced by climate warming (Jeffs and Lewis 2013). Microcosm experiments on Drosophila and their parasitoids have revealed that both the direct and indirect biotic interactions that determine the distribution and abundance of species are temperature-dependent (Davis et al. 1998). In a further example, significant shifts in the distribution of species following the establishment of a range-expanding species attest, on a larger scale, of changes in interspecific interactions in response to warming (Audusseau et al. 2017). Specifically, the distribution of the native butterfly species Aglais urticae and Aglais io in Sweden shifted following the arrival of a newly-established third butterfly species Araschnia levana which shares the same larval host-plant (Audusseau et al. 2017), with evidence that these shifts were the result of ongoing apparent competition, mediated by shared parasitoids.
Here we focus on a community of closely-related (Nymphalidae: Nymphalinae, Nymphalini) nettle-feeding butterflies (Ag. urticae, Ag. io, Ar. levana, Vanessa atalanta) and their shared parasitoids. We identified the complex of parasitoid species involved in the interactions among these butterfly species and examined how interactions varied over space and time in a field study spanning a 500 km latitudinal gradient in Sweden. This latitudinal gradient, which follows the range expansion of Ar. levana, is also an opportunity to investigate further the impact of this butterfly on the resident nettle-feeding butterflies. In that context, we investigated the signature of apparent competition within this community of nettle-feeding butterflies by studying the phenology of parasitism, its spatio-temporal structuring and the role of species co-occurrence. We also explored the potential impact of the arrival of Ar. levana in Sweden by studying how parasitism in the resident butterfly species covaries with the presence or absence of the newly-established butterfly species. In particular, we showed that the phenology and the spatial structure of parasitism rate vary with changes in butterfly host species assemblage, and how biotic interactions mediated by shared parasitoids might change when a novel expanding host enters the community. Our work thus provides novel insights on the relationships between herbivorous insects and their parasitoid complexes and how they vary latitudinally.

Study system
Aglais urticae, Aglais io, Araschnia levana and Vanessa atalanta are closely-related butterfly species from the same tribe (Nymphalini, family Nymphalidae, Wiemers et al. 2020). The larvae of all four species feed (practically exclusively) on nettle Urtica dioica, but the butterflies differ in their egg-laying behaviour, phenology and distribution (Eliasson et al. 2005).
Aglais urticae, Ag. io and Ar. levana are batch-laying species, with batches of 10-40 eggs for Ar. levana and of 200-300 eggs for Ag. urticae and Ag. io, while V. atalanta lays eggs singly (Ebert 1993). During the first three instars of their development, the larvae of Ag. urticae and Ag. io are gregarious and conspicuous as they feed near the apex of the nettle stem. At their fourth instar, larvae of Ag. urticae and Ag. io become solitary and feed over the entire plant and may hide in the foliage. Larvae of Ar. levana are also gregarious in the early instars and become solitary later in development. However, batches of this species are less conspicuous to the human eye, not only because the smaller batch size and larvae both result in less damage to the plant, but also because the larvae often feed from the lower surface of the leaf.
In Sweden, where we carried out our study, the four species have broadly overlapping phenologies, with adults flying from March to September. However, the periods during which larvae of each species are found in the field vary due to differences in voltinism (the number of generations per year) and yearly variations in weather conditions. Populations of Ag. urticae are bivoltine in Sweden. Larvae of this species are recorded from early May to the end of August with the first generation being found from early May and the second generation from late June. Aglais io is univoltine in Sweden and starts reproducing soon after Ag. urticae, with larvae observed from late May to early August. Araschnia levana is an obligate bivoltine species. In contrast to Ag. urticae and Ag. io, which overwinter as adults, Ar. levana hibernates in the pupal stage. Larvae from the first generation are found in the field in June; larvae from the second generation are found from the end of July to early September. Last, V. atalanta is a migratory butterfly in Sweden and its population depends on the migratory influx from the areas where the species is resident. It is univoltine in Sweden with larvae observed in the field from May to early September.
All species are distributed across Sweden, except for Ar. levana which is currently limited to the southern half of the country (Eliasson et al. 2005). Araschnia levana is a recent colonist with the first anecdotal observations reported in Skåne in 1982. It has progressively expanded from the county

Field sampling
We hand-collected larvae of the four study butterflies, Ag. urticae, Ag. io, Ar. levana and V. atalanta, over two years (2017)(2018) and fortnightly throughout the species' reproductive season (May-August). Sampling was conducted in 19 sites distributed along a 500 km latitudinal gradient from south Sweden to the Stockholm area (Fig. 1). The 13 sites located in the southern part of Sweden fall within the distribution range of all four butterfly species, while the six sites in the Stockholm area are north of Ar. levana's current range. As butterflies and parasitoids are influenced by their surrounding habitat (Shaw 2006, Van Halder et al. 2017), we selected sites within comparable landscapes (i.e. the proportion of forest and agricultural land within 1 km 2 ). Early in the first season (2017), we located a set of nettle patches at each site that we then visited for each sampling event in both years (see the Supporting information for details on the field sampling). Thus, the sampling effort was constant over time, but varied across sites as it was proportional to the number of nettle patches found at each site.

Larval sampling and monitoring
Although pupal parasitism is likely to cause high mortality in the species studied (Pyörnilä 1977, Shaw et al. 2009), solitary and concealed pupae are difficult to collect in sufficient numbers for reliable estimates of pupal mortality. We therefore focused on larval parasitism, sampling butterfly larvae at different stages. We stratified our sampling effort across larval stages to capture a representative sample of the parasitoid species diversity, as the temporal window of attack on butterfly larvae differs among parasitoid species and can be restricted to a few developmental stages. For example, ichneumonids of the genus Thyrateles attack very late larval or prepupal stages (Shaw unpubl.). In contrast, Cotesia vestalis, which can be an important opportunistic parasitoid of at least Ag. urticae, although not found in this survey, parasitises first instar larvae and emerges mainly from second instar larvae (Shaw unpubl.). Therefore, at each sampling occasion and for each butterfly species we aimed to collect seven second instar larvae per batch from a maximum of five batches, 20 fourth instar larvae per batch from a maximum of five batches, and up to 20 fifth instar larvae from different batches, where possible (see the Supporting information for details on the number of butterfly larvae of each species sampled per larval stage). Note that as we rarely reached these maxima, the number of larvae and butterfly nests collected is proportional to the abundance of the species at each site and sampling occasion. Larvae were reared in transparent plastic boxes (155 × 105 × 45 mm). Depending on the number of larvae collected per batch, the larvae were raised alone or in groups of up to five larvae from the same batch. We reared larvae under laboratory conditions (temperature 23°C, light regime 22L:2D) and fed them daily with Urtica dioica leaves collected from the same location as where the larvae had been sampled. This was because some of the Tachinidae parasitoids (Sturmia bella and Pales pavida, of those encountered) lay microtype eggs on nettle leaves and the butterfly larvae become parasitised only when they eat the infected leaves.
We recorded the date and stage from which the parasitoid emerged from the parasitised larvae and kept the parasitoids individually or per batch in plastic vials, under the same laboratory conditions as the butterfly larvae. We preserved freshly-dead adult parasitoids in 95% alcohol, before taxonomic identification. The parasitoid pupae that did not hatch by early September, as well as the pupae from the second generation of Ar. levana (which have an obligate diapause before adult emergence), were kept cold over the winter, until we broke their diapause around mid-April (see the Supporting information for details on the diapause conditions). The taxonomic determination of our samples relied on the experience and the good knowledge of all these species and their regular hosts held by the authors for the Ichneumonoidea (MRS) and the Tachinidae (CR), including an awareness of the wider taxonomic literature (see the Supporting information for details on the taxonomic identification).

Analyses
We performed all analyses in R ver. 3.5.1 (<www.r-project. org>). Parasitism rate was modelled as a binomial response with a two-vector variable equivalent to Bernouilli trials of 'success' and 'failure', where success is the number of parasitised larvae and failure is the number of non-parasitised larvae in each batch. For the analyses performed in a Bayesian framework, we used generalised linear and non-linear multivariate multilevel models. We fitted the model through Markov chain Monte Carlo (MCMC) sampling, using the Hamiltonian Monte Carlo algorithm implemented in Stan (Carpenter et al. 2017) and the R interface provided in the brms package (Bürkner 2017(Bürkner , 2018. We ran four chains for 10 000 iterations, with the first 4000 discarded as burn-in and a subsequent thinning of 2, and used the default noninformative priors. To test for significant differences in parasitism between groups, we compared the posterior probability distributions of the model parameters.

Effects of latitude and phenology on parasitism rates
We investigated variation in overall parasitism rates per butterfly species and county (Skåne, Kronoberg and Stockholm). We performed this analysis in a Bayesian framework. Parasitism was modelled assuming a binomial distribution and a logit link function. We tested for the effect of species, county, year, and the interaction between species and county on parasitism rate and included the week of sampling as a non-linear effect (with k up to 4) to control for phenological variations in parasitism rate for each species. We grouped sites by county (Skåne, Kronoberg and Stockholm, Fig. 1) to reflect the south-north progression of the establishment of Ar. levana and increase the power of our analysis along this gradient.

Effect of butterfly species assemblage on parasitism rates
We examined the effect of the butterfly assemblage on each species' parasitism rates. Specifically, we investigated how parasitism rates of each butterfly vary with the presence or absence of each of the other butterfly species, coded as a binary variable (0/1), and the total abundance of butterfly larvae. We also included in each model the non-linear effect of the sampling week (with k up to 4). The presence and absence of each butterfly and the abundance of butterfly larvae at a site were extracted from the data collected for each sampling occasion. The total abundance of butterfly larvae corresponds to the total number of larvae from all butterfly species collected per site and sampling week was zero-centred prior to inclusion in the models. We performed these analyses in a Bayesian framework. Parasitism was modelled assuming a zero-inflated binomial distribution with a logit link function. Araschnia levana experienced low levels of parasitism (24 out of 160 batches of Ar. levana sampled were parasitised), impeding analysing the impact of the butterfly assemblage on parasitism for this species.
Note that this analysis examines the effect of the butterfly assemblage on each species' parasitism rates, regardless of the parasitoids responsible for the parasitism rate recorded. Since the parasitoids responsible for the highest mortality were partially or entirely shared between the study butterflies (Table 2), this analysis explores how each butterfly species' cooccurrence was linked to the parasitism of each of the focal butterfly. In the Supporting information we provide a similar analysis but focused on the subset of parasitoids shared with Ar. levana. By limiting this analysis to parasitoids shared with Ar. levana, we study how the parasitism rates for each butterfly covary with the presence or absence of the newlyestablished species and the potential impact of Ar. levana on the resident butterflies through its influence on the population dynamics of shared parasitoids.

Parasitism rate and time since the establishment of Ar. levana
We explored the potential role of the time since Ar. levana has established on parasitism rate of the native species. The available observations of Ar. levana suggest that the species has first established in the southern part of the country and has since spread northward (Supporting information). If the establishment of Ar. levana has induced an increase in parasitism rate of the native species through apparent competition (as proposed by Audusseau et al. 2017), we would expect a decrease in parasitism rates with increasing latitude. In addition, the establishment of Ar. levana and its progression might not have strictly followed the latitudinal gradient as a result of the configuration of the landscape such as the presence of corridors or barriers affecting dispersal. We therefore also tested the effect on parasitism rate of the time since the first observation of Ar. levana within a 10 km buffer zone around each site with the hypothesis that parasitism would be higher in the earliest colonized sites.
For each species, we tested the effects of latitude and time since the first observation of Ar. levana in the 10 km the buffer zone around the sites, using generalised linear models and assuming a binomial distribution. We restricted these analyses to sites where Ar. levana has established (Skåne and Kronoberg). Time since colonisation was extracted from Artportalen (Swedish Species Observations System, <www. artportalen.se/>, 30/08/2019). Because the latitude and the time since the first observation of Ar. levana at a site are highly correlated, both reflecting the south-north gradient of progression of Ar. levana, we transformed the time elapsed since the first observation of Ar. levana into a 4-level ordinal variable. By dividing the variable into four quartiles, we grouped the sites by periods of establishment of the expanding species. The dates that divide the first, second and third quartiles, 16-05-2004, 22-07-2006 and 02-08-2007, follow a sigmoid curve as would be expected from the establishment and expansion processes. Latitude was zero-centred before it was included in the model.

General patterns of incidence of butterfly species and parasitoid attack
Over the two sampling seasons, we sampled 6777 butterfly larvae across the 19 sites (Aglais io = 2259, Ag. urticae = 2254, Araschnia levana = 1583, Vanessa atalanta = 681). The three resident butterfly species (Ag. io, Ag. urticae, V. atalanta) occurred in each region. However, Ag. io was absent at three sites (Odensjö, Åsvägen, 31), and Ag. urticae was absent from site 31. As expected, Ar. levana was found at all sites in the two southern counties but not at the latitude of the Stockholm area.
Of the 6777 collected larvae, 1508 were parasitised and produced parasitoids from three families: Tachinidae (Diptera), Ichneumonidae (Hymenoptera) and Braconidae (Hymenoptera  Fig. 2). Phryxe vulgaris and M. subcompleta were present Table 1. Distribution of larva death according to counties and sampling sites, by parasitoid family and species, and due to unknown causes, which cover infection by virus, bacteria or fungi. The sites are ordered from high to low latitude. Note that five larvae were simultaneously parasitised by two species, which lead to the discrepancy between the total by family and the grand total. Larval death by:  at most of the sampling sites and across the three counties (Table 1, Fig. 2). We recorded other parasitoid species in low numbers, providing limited information about their latitudinal distribution. Still, T. haereticus (n = 21) was restricted to the two northern counties and C. vanessae (n = 30) to the two southern counties (Table 1, Fig. 2). The parasitoid complex varied among the butterfly hosts. Vanessa atalanta was the host of the most parasitoid species, including representatives of all three families (Table 2). Aglais urticae was also found to be parasitised by a wide range of species from the three families (Table 2). Aglais io and Ar. levana were not parasitised by braconids and Ar. levana was almost exclusively parasitised by S. bella, except on two occasions by Pho. confusa ( Table 2). The three most abundant parasitoids, Pe. tibialis, Pho. confusa and S. bella, were shared between the resident butterfly hosts and Ar. levana, except for Pe. tibialis which was never observed parasitising Ar. levana. We also recorded cases where the cause of larval death was unknown. Although such mortality can sometimes be related to parasitoids that failed to complete their development within the body of their host, either due to a late attack of the parasitoid or to the immune response of their host (HA, pers. obs.), we also recorded cases of mortality due to viral, bacterial or fungal infection. The overall percentage of dead larvae due to unknown causes varied from 4.2% for Ar. levana to 19.8% for Ag. io ( Table 2). The high mortality of Ag. io is not surprising as this species is relatively sensitive to laboratory rearing conditions, especially during the early instars (Audusseau unpubl.).

Effects of latitude and phenology on parasitism rates
Parasitism was responsible for high mortality, particularly in Ag. urticae and Ag. io (Fig. 3a, Supporting information) and showed a gradual decrease towards higher latitudes, from Skåne to Stockholm (Fig. 3a). Over the two field seasons, 40.2% of Ag. urticae and 37.0% of Ag. io larvae collected in Skåne were parasitised. These rates decreased to 20.4% and 17.4% in Stockholm County for Ag. urticae and Ag. io, respectively. Aglais urticae showed higher parasitism rates than Ag. io, although this effect was driven mainly by the difference observed in the Stockholm area (Fig. 3a, Supporting information). Across counties, Ag. urticae and Ag. io were parasitised at a significantly higher frequency than V. atalanta and Ar. levana (Fig. 3a). Over the two field campaigns, parasitism rate in V. atalanta was highest in Skåne, of 39.9%, than in the counties of Kronoberg and Stockholm, where it was of 12.0% and 13.1%, respectively. Araschnia levana was very weakly parasitised, with parasitism rates of 4.1% in Skåne and 3.9% in Kronoberg.
The overall parasitism rate was significantly lower in 2017 compared to 2018 (estimated difference between 2017 and 2018 = −0.33, 95% credible interval (CI) = [−0.49, −0.17], Supporting information). Furthermore, parasitism rates were seasonal and specific to each butterfly species (Fig. 3b, Supporting information), as a result of differences in their phenology and the phenology of their parasitoids. Parasitism rate in Ag. urticae followed a bimodal distribution that reflects the bivoltine life cycle of the species in Sweden. In contrast, the parasitism rate in Ag. io Table 2. Numbers of dead larvae per butterfly species caused by parasitism or due to unknown causes, covering infections by virus, bacteria or fungi. The table also summarises the contribution of each parasitoid species to the total parasitism found per butterfly species and intermediate summaries show parasitoids contribution by family. The percentages of larvae dead due to unknown causes are related to the total amount of larvae of each sampled species. The single rearing of Microgaster subcompleta from Aglais urticae is exceptional for this parasitoid (MRS pers. obs.), and Ag. urticae should be regarded as outside its host range (Shaw 1994   and V. atalanta followed a unimodal pattern with a peak at the end of July. We observed a similar unimodal pattern of parasitism in Ar. levana, but the low parasitism rate observed in this species makes it difficult to derive reliable estimates of its phenological variations.

Effect of butterfly species assemblage on parasitism rates
Parasitism rate by species varied with the butterfly species assemblage at the time of collection. That is, the rates in Skåne. Non-overlapping credible intervals correspond to significant differences in parasitism rate between groups. Note that we have adjusted for week 4.4 as at this time, differences in parasitism between species reflect the overall differences observed the season. The phenology of parasitism is illustrated in Skåne but follows the same pattern in the other two counties, modulated by a variation in the intercept. The red line on (b) indicates week 4.4, being the time of the reproductive season for which the marginal means shown in (a) were extracted for Skåne. We only plotted the estimated variation in parasitism rate over the time window for which each species was sampled in the field.  . Contrasting effects of butterfly species assemblage, taken as the presence/absence of the other species, including Araschnia levana, on parasitism rate of (a) Aglais urticae, (b) Aglais io and (c) Vanessa atalanta. Estimation of marginal means of parasitism rates (%) are given at representative values (week = 4.74) and parasitism rates of each of the focal species are ordered on the x-axis according to the number of species which co-occur (also visually recognizable with the gray scale). The first bar on each plot corresponds to parasitism rate of the focal species found alone (mean ± 95% credible interval) and the letter stands for the identity of the focal species with A for Ag. urticae, B for Ag. io, C for V. atalanta. The following bars correspond to parasitism rate of the focal species (mean ± 95% credible interval) when cooccurring with other nettle-feeding butterflies, with +A when the species co-occur with Ag. urticae, +B with Ag. io, +C with V. atalanta and +D with Ar. levana. Non overlapping credible intervals correspond to significant differences in parasitism rate between groups. of parasitism are species-specific and vary with the number and species identity of co-occurring larvae, as well as the total abundance of larvae at a given time (Fig. 4, Supporting  information).

Parasitism rate and time since the establishment of Ar. levana
The time period since the first observation of Ar. levana significantly explained variations in parasitism rate in Ag. urticae and Ag. io (LR χ 2 (3) = 35.15, p < 0.001 for Ag. urticae and LR χ 2 (3) = 15.88, p = 0.001 for Ag. io, Fig. 5, Supporting information), which showed higher parasitism rates in the earliest colonised sites. In addition to the effect of time, parasitism observed in Ag. io also decreased with latitude (LR χ 2 (1) = 5.22, p = 0.022, Supporting information). Parasitism in V. atalanta was not explained by differences in the time period since the first observation of Ar. levana but decreased with latitude (LR χ 2 (1) = 24.61, p < 0.001, Supporting information).

Discussion
Our results highlight the influence of butterfly assemblages on the parasitism of nettle-feeding butterflies. We showed that parasitism was responsible for high mortality rates in two of the native species, Aglais urticae and Ag. io. In comparison, parasitism caused relatively lower mortality in Vanessa atalanta and Araschnia levana. Although the commonest parasitoid of V. atalanta, Microgaster subcompleta, was not found in the other species (Table 2), the parasitoid complex was largely shared among the nettle-feeding butterflies, but Ar. levana, the newcomer in Sweden, was almost exclusively parasitised by the tachinid Sturmia bella. We observed that parasitism was influenced by the butterfly assemblage and that this effect was specific to each butterfly species. In addition, we found that at sites where Ar. levana has established for a longer time period, parasitism rates of the native species were significantly higher.
The low parasitism rate in V. atalanta and Ar. levana relative to the other butterfly species might be the result of morphological, physiological, behavioural, and immunological differences, compared to the other study species. Vanessa atalanta larvae are solitary, which may complicate the search for host larvae by their parasitoids, in comparison to the other species that lay batches of eggs Dyer 2002, Hawkins 2005). However, V. atalanta larvae also live concealed in folded leaves, a shelter-building behaviour that has been shown to concentrate chemical and visual signals that facilitate the localisation of individual larvae by parasitoids Gentry 1999, Sugiura 2007). Nevertheless, the low parasitism measured in V. atalanta is difficult to explain, considering that this species was host for the largest diversity of parasitoids and that it has been documented to be highly parasitised in other parts of its range (Rice 2012). Variation in V. atalanta parasitism rates across its range may be related to its migratory behaviour, conditions at overwintering sites, and synchrony between the butterfly and its parasitoids. For example, its migration  -2004, 22-07-2006 and 02-08-2007, respectively. phenology is dependent on the weather conditions locally (Brattström et al. 2008), which is therefore unlikely to be synchronised with the phenology of the parasitoids at the reproductive sites. The pattern is different for Ar. levana, which is resident in Sweden and has been found to be weakly parasitised in other parts of its distribution (Wagner et al. 2011). Arachsnia levana larvae show a pronounced dropping behaviour, which in other species has been shown to be effective against parasitoids which lose track of the chemical and sensory cue of their hosts (Gross andPrice 1988, Fitzpatrick et al. 1994). Alternatively, lower parasitism in Ar. levana could be a result of its recent establishment in Sweden. The enemy release hypothesis (Jeffries andLawton 1984, Keane andCrawley 2002) predicts that, in a new area, species experience a period when they escape their natural enemies until interactions with the local parasitoid complex are established (Menéndez et al. 2008). In Sweden, Ar. levana was first reported in the 1980s but probably became established more recently, as there are very few reports of the species before 2000 (Supporting information). Considering the relatively short time that was available for recruitment of local parasitoids (Cornell and Hawkins 1993), we cannot rule out the possibility that lower levels of parasitism observed in Ar. levana are partly a consequence of its recent establishment and that populations have escaped from parasitism during its expansion phase. This hypothesis is strengthened by the fact that Söderlind (2009) reported no parasitism in Ar. levana in south Sweden, while our data reveal that the species is now parasitised by local parasitoid populations. Future monitoring of parasitism load in Ar. levana populations in Sweden and across the wider distribution range of the species would be necessary to disentangle the relative importance of these two possibilities.
Butterfly species assemblage was significantly associated with parasitism in Ag. urticae and Ag. io but not in V. atalanta. The differences in egg-laying behaviour mentioned above, where parasitoids prefer gregarious species when present, is again one potential explanation. However, V. atalanta was mostly parasitised by M.. subcompleta, a parasitoid associated almost exclusively with this species in our study, limiting the influence of indirect interactions mediated by parasitoids shared with the other butterflies. Aglais io seemed to benefit from the co-occurrence of Ag. urticae, by showing reduced parasitism, while parasitism in Ag. urticae was higher when it co-occurred with Ag. io, suggesting that these two species can interact indirectly through their shared parasitoids. Most interestingly, parasitism in Ag. io was higher when co-occurring with Ar. levana and we observed a similar impact of the co-occurrence of Ar. levana on Ag. urticae when restricting our analysis to only shared parasitoids (Supporting information). Our data did not allow us to quantify the proportion of parasitoids that had been recruited from the newly-established butterfly and thus its indirect impact on the parasitism of the resident butterflies as proposed by Müller et al. (1999). Nonetheless, this latter result, and the association between the decrease in the parasitism of the native species and the progression gradient of Ar. levana, where the native species showed lower parasitism rate in sites recently colonised by Ar. levana, are consistent with previous work conducted on this study system which hypothesised that the arrival of Ar. levana would influence the dynamics and spatial distribution of the resident butterflies through apparent competition (Audusseau et al. 2017). Furthermore, parasitism in Ag. urticae increased with the total abundance of larvae, a phenomenon that might partly be associated with the arrival of the novel host. Differences in the phenology of parasitism between hosts also suggest that Ar. levana could provide a refuge for parasitoids at a time when the native hosts (Ag. urticae and Ag. io) are rare. Thus, species co-occurrence at a site over the season, rather than at a sampling event, may influence their level of parasitism. Although our study focused on larval parasitoids (for reasons previously mentioned), it is important to note that pupal parasitoids are also known to be shared among our study butterflies and to cause high mortality (Pyörnilä 1977, Shaw et al. 2009). In particular, the restricted host range of the pteromalid Pteromalus puparum, which includes the butterflies of our study (Shaw et al. 2009), and the size of its brood, make this species a strong candidate for driving apparent competition in our study community. All these elements constitute additional evidence in favor of an increase in the realised host repertoires of the parasitoid species (which gradually integrate Ar. levana as a host), to explain the high parasitism rates of the native species in the southern counties.
From our study, we cannot rule out the possibility of other spatial confounds causing the differences in parasitism rates across counties. Among the factors not considered here, the differences in parasitoid population dynamics, habitat quality, parasitoids community or the variation in phenological synchrony between the butterflies and their parasitoids (Audusseau et al. 2020), may all contribute to the observed decrease in parasitism rates along the latitudinal gradient. For example, the occurrence of other hosts for the parasitoids in the landscape, but not included in our study, may influence the population dynamics of parasitoids and mediate apparent competition (Davis 1991, Gaston et al. 2005. Parasitoids are also responding to the conditions of their habitat (Shaw 2006), which may vary between counties, despite our effort to select sites within comparable landscapes. The latitudinal decrease in parasitism could also be associated with a latitudinal trend in weather conditions as temperature affects insect-parasitoid interactions (Thomas and Blanford 2003). While in some systems parasitoid activity can increase with temperature (Mann et al. 1990), which could lead to a higher activity period and oviposition rate in the parasitoids at lower and warmer latitude, the literature does not provide consistent evidence of such a pattern (Hawkins 2005). Further, the recorded differences in microclimatic conditions across sites do not align with the latitudinal pattern observed for parasitism (Supplementary material Appendix 1 in Audusseau et al. 2020). However, we observed latitudinal differences in the parasitoid community. For example, Sturmia bella, one of the most abundant parasitoid species in our sample, was only found in the two southern counties. While no evidence testifies of its role, S. bella has also recently established in the UK and its arrival coincided with the decline of Ag. urticae (Gripenberg et al. 2011).
Although rarely conducted, it is only through concerted surveys, conducted across spatial scales and over years, on a set of co-occurring species, that we can hope to assess the full impact of environmental change on populations and communities. Here, our systematic sampling enabled us to study how species assemblage, variation in abundance, and species phenology influence local biotic interactions. We showed that parasitoid pressure has a substantial effect on the mortality rate of nettle-feeding butterflies in Sweden. This suggests that the responses to environmental change, such as the arrival and establishment of a new species, may potentially have cascading effects on resident host species, even if they do not interact directly. In that respect, we showed that although weakly parasitised, the occurrence of Ar. levana in the landscape can increase the abundance of available host butterflies and result in increased parasitism at the community level. The complexity and dynamic nature of biotic interactions in trophic networks demonstrated in our study are all the more important in common species as their ubiquity and abundance often make them connect with a large number of species through trophic interactions (Gaston 2010).