Increased fluctuation in a butterfly metapopulation leads to diploid males and decline of a hyperparasitoid

Climate change can increase spatial synchrony of population dynamics, leading to large-scale fluctuation that destabilizes communities. High trophic level species such as parasitoids are disproportionally affected because they depend on unstable resources. Most parasitoid wasps have complementary sex determination, producing sterile males when inbred, which can theoretically lead to population extinction via the diploid male vortex (DMV). We examined this process empirically using a hyperparasitoid population inhabiting a spatially structured host population in a large fragmented landscape. Over four years of high host butterfly metapopulation fluctuation, diploid male production by the wasp increased, and effective population size declined precipitously. Our multitrophic spatially structured model shows that host population fluctuation can cause local extinctions of the hyperparasitoid because of the DMV. However, regionally it persists because spatial structure allows for efficient local genetic rescue via balancing selection for rare alleles carried by immigrants. This is, to our knowledge, the first empirically based study of the possibility of the DMV in a natural host–parasitoid system.


Introduction
The frequency of weather extremes is increasing under ongoing climate change [1]. One repercussion of extreme events is increased spatial synchrony of local populations [2], which can decrease stability regionally and alter biotic interactions and community structure [3]. Small populations are vulnerable to a cascading feedback between demography and loss of genetic diversity, leading to extinction via the extinction vortex [4]. The risk of extinction increases especially in species at high trophic levels, as resources they depend on become increasingly sparse and locally unstable [5,6]. As resources become scarce, declining populations may lose genetic diversity, further increasing the risk of extinction [7,8]. In insect communities, parasitoids occupy the higher trophic levels. They interact strongly with their hosts, and those with limited host ranges have smaller population sizes than do their hosts [9].
Parasitoid wasps lay eggs in or on other arthropods. The larvae then consume and eventually kill their hosts. As all other Hymenoptera, they are haplodiploid. Females are diploid, developing from fertilized eggs, and males are haploid, developing from unfertilized eggs. Inbreeding can be detrimental for haplodiploid species that have single locus complementary sex determination (sl-CSD). In these species sex is determined by a single locus, wherein haploid individuals develop as males and diploid individuals develop as females when heterozygous [10]. mating results in half the fertilized eggs being homozygous at the sex locus and developing into diploid males at the expense of female offspring (figure 1).
Diploid males are inviable [11] with a few exceptions [12,13]. Thus, production of diploid males is costly to the parents, and in the case of sterile diploid males, to their mates [14,15]. sl-CSD is the most prevalent and presumably the ancestral mode of sex determination in Hymenoptera [10,16,17]. Many parasitoid species produce diploid males under laboratory conditions (e.g. [12,18 -21]). They are also found in natural populations [20,22] and after introduction for biological control [23 -25]. Zayed & Packer [26] demonstrated theoretically that high production of diploid males in small and isolated populations of species with sl-CSD can elevate extinction risk through a 'diploid male vortex' (DMV), in which accelerating decline in population size triggers positive feedback between small population size and loss of alleles at the CSD locus. The DMV increases the risk of extinction in haplodiploid over diploid species living in similar ecological circumstances. Later theoretical studies challenged the assumptions made by Zayed & Packer [26], showing that DMVs require a set of stringent conditions if relevant behavioural, demographic, and ecological factors are included [27][28][29]. Particularly, a small amount of dispersal in spatially structured populations promotes the maintenance of CSD alleles through strong balancing selection for rare alleles [27,30]. In single-species population dynamic models incorporating CSD [19,25 -27,29], resource availability (carrying capacity) is kept constant over time, but large population fluctuations could cause genetic bottlenecks leading to loss of CSD alleles from the population [25]. The existing models have not incorporated both temporal host population fluctuations (fluctuation in resource availability) and spatial population structure. Such extensions are necessary to understand how populations constrained by CSD may fair with increased habitat fragmentation and fluctuations in environmental conditions that are anticipated under climate change.
Existing empirically based studies of diploid male production (DMP) do not come to a consensus on its significance for the persistence of natural populations [24,[31][32][33][34]. So far, evidence for DMVs occurring in the wild has not been demonstrated. To do so, sampling would have to be done over time. We studied DMP in a hyperparasitoid that is part of the insect community associated with the Glanville fritillary butterfly, Melitaea cinxia (Lepidoptera: Nymphalidae) in the Å land islands, southwest Finland [35], over a period of four years. The host butterfly lives as a classic metapopulation inhabiting 300-500 small meadows over an area of 3500 km 2 [36,37]. In recent years, spatial synchrony of the local butterfly population sizes has increased in Å land. This appears to be owing to increasing extreme and spatially synchronized weather, leading to spatial autocorrelation of population sizes. This synchrony is due to local rates of survival and reproduction and not because of homogenizing dispersal in the landscape [38]. The increasing butterfly population fluctuation must affect the community of parasitoids that depend on it, and indeed, we have observed a decline in the effective population size (N e ) of the hyperparasitoid [39]. Here, we address the population level consequences of DMP for the hyperparasitoid, using empirical data from large-scale field sampling over time and a spatially-structured multitrophic population model. Our study demonstrates the long-term impact of increased fluctuations in the population dynamics of the host butterfly (the second trophic level) on the population size of a hyperparasitoid (the fourth trophic level), mediated through both demographic processes and loss of genetic diversity at the CSD locus.

Material and methods (a) Study species
Mesochorus cf. stigmaticus (Hymenoptera: Ichneumonidae) is a hyperparasitoid wasp that lays eggs only into larval parasitoids developing inside M. cinxia caterpillars. The wasp is found throughout the butterfly metapopulation in Å land [39], where it uses almost exclusively the host Hyposoter horticola (Hymenoptera: Ichneumonidae) which is a specialized solitary egg-larval parasitoid of the butterfly [35,40]. A peculiarity of this system is that the primary parasitoid (H. horticola), under natural conditions consistently parasitizes about a third of the butterfly larvae in each gregarious nest (figure 2b). Along with its strong dispersal ability [41], this behaviour translates to a uniform rate of parasitism over the landscape independent of local or regional host density [40,42]. As a result of the lack of local density dependence of the primary parasitoid, the population dynamics of the butterfly can be considered as directly influencing the hyperparasitoid.

(b) Sample collection, rearing and genotyping
Melitaea cinxia larvae live gregariously in silk nests. Three larvae, some of which were naturally parasitized and hyperparasitized, were collected from each nest during the annual autumn survey of Å land [36] over four years (autumn 2008 -2011). The number of nests sampled ranged from 1224 to 4689 owing to the yearly fluctuation of butterfly population size. The collected larvae were reared in the laboratory until they became butterflies, parasitoids emerged from them, or they died. The number of adult hyperparasitoids reared differed between years (n ¼ 48 -361), depending on the host sample size, overwintering mortality, and the rate of hyperparasitism. Adult wasps were stored in 96% ethanol at 2208C for genetic analysis. The samples were used to estimate the rates of parasitism and hyperparasitism in the population, and to determine the fraction of hyperparasitoid males that were diploid. DNA was extracted from tissue of males using DNeasy isolation kit (Qiagen, Hilden, Germany), and genotyped using species-specific microsatellite loci [43]. The detailed genotyping methods are described in Nair et al. [39].

(i) Data analyses
Diploid males were identified based on genotyping results from 17-25 microsatellite markers [43]. An individual was identified as a diploid male if it lacked an ovipositor, and one or more microsatellite marker was heterozygous. The percentage of diploid male for each year was estimated as the fraction of males that were diploid in the field samples. The number of hyperparasitoid migrants per generation between survey areas was estimated in GENEPOP (ver. 4.2) [44] using neutral microsatellite data (total wasps ¼ 175) from Nair et al. [45].

(c) The model
We developed a discrete-time individual-based simulation model of spatially structured populations of a hyperparasitoid including explicit genetics for CSD. It simulates the population dynamics of the hyperparasitoid in response to fluctuation and spatial autocorrelation of local butterfly populations. The model is formulated using the known long-term butterfly dynamics [36,38] and parasitoid population dynamics, genetic structure and natural history [35,39,40]. Table 1 lists the state variables (a) and model parameters (b). More detail can be found in the electronic supplementary material. The model was implemented in MATLAB R2017a. The code is available from the Dryad Digital Repository: https://doi.org/10.5061/dryad.56qf11h [46].

(i) The spatial structure in the model
The study area was divided into 12 sub-regions, based on the butterfly survey areas (electronic supplementary material, figure S1; [36]) to represent realistic heterogeneity in host availability to the hyperparasitoid arising from spatial distribution of habitat patches for the butterfly. The hyperparasitoid is dispersive [39] but individual wasps (approximately 6 mm in size) are dispersal limited and their distribution in the landscape is limited by the heterogeneity of host availability [35]. We estimated the butterfly larval abundances annually using the autumn survey [36] data from 2003 to 2016. The populations in four sub-regions were designated as small (average number of larvae , 3000), four were medium (between 3000 and 6000), and four were large (greater than 10 000). To account for uncertainty in dispersal patterns of the hyperparasitoid

(ii) The butterfly larvae sub-model
The mean population size of butterfly larvae at time t in subregion i, B i,t , is drawn from a multivariate log-normal distribution with the means and the variance -covariance matrix of the associated multivariate normal distribution estimated from the annual autumn survey data. The time step is 1 year (generation) and, as until recently in the survey data, there is no temporal correlation between years [38]. The actual number of larvae in sub-region i, B i,t , was determined by drawing values from a Poisson distribution with mean equal to B i,t . Based on the survey data, the estimates for the correlation of nest counts among the sub-regions ranged from about 0.25 to 0.75, and the standard deviations increased 1.85-fold between 2003 -2009 and 2010 -2016. The spatial variation was due primarily to spatially uncorrelated weather [38]. We simulated variability in butterfly population fluctuations by varying the cross-correlation coefficient, r, and the factor M multiplying the standard deviations in calculating the variance-covariance matrix for population dynamics in the sub-regions.

(iii) The host parasitoid sub-model
The population size of the parasitoid before hyperparasitism at time t in sub-region i, N 0 i,t , is drawn from a binomial distribution with the butterfly larval population size as the sample size and probability of success equal to 1=3 [40,42]. We modelled the mean rate of hyperparasitism at time t in sub-region i, f i,t , as a saturating function of the ratio between hyperparasitoid female and parasitoid densities, and visually fit the function to empirical data (electronic supplementary material, E1 and figure S2). We used the data as a guide to estimate the parameters because the functional form is theoretically reasonable and available data are limited. To reflect yearly stochastic variation, the actual rate of hyperparasitism was drawn from a normal distribution with the mean equal to f i,t with standard deviation 0.1 to match the variation seen in the empirical data (electronic supplementary material, E2 and figure S2). The number of parasitized hosts at time t, N p i,t (p for parasitized) is drawn from a binomial distribution with probability of success equal to f i,t and the number of trials equal to the number of parasitoids.

(iv) The hyperparasitoid sub-model
In the model, each individual hyperparasitoid is represented by one sex locus. We assume diploid males to be sterile. An individual may disperse once with equal probability h (migration rate) to any directly neighbouring sub-region (two-dimensional stepping-stone dispersal [47,48]). We used three levels of h, 0.002, 0.014, and 0.027 (electronic supplementary material, table S1) to account for the uncertainty and variation of migration rate observed in the temporal genetic structure of the wasp [39]. After the dispersal step, the hyperparasitoid mates in the natal or destination sub-region. Mates are randomly paired and polyandry is assumed. The total number of offspring at time t in sub-region i is N p i,t . The parental genotypes (sl-CSD allele(s)) are randomly assigned to each offspring and passed on according to Mendelian inheritance for haplodiploidy. Each allele mutates into a new allele with a probability m of 10 À7 [25,49].

(v) Simulation experiments
We ran 30 replicates for each parameter set for 10 000 generations to ensure that transient dynamics disappeared and that sufficient numbers of local extinctions were observed to distinguish decline owing to demographics from that owing to CSD. We initialized the population with 10 alleles at the CSD locus, randomly distributed among hyperparasitoid individuals. We chose 10 to achieve [17,39,50] the 9-10% of diploid males present at the low level of fluctuation in butterfly population size observed between 2003 and 2009. We varied the multiplier on the level of fluctuation, M, to cover a parameter region from zero to three times the baseline (2009) fluctuation of the butterfly population sizes. The correlation parameter r was varied from 0 to 1. We summarized outputs from the last 5000 generations at the whole Å land scale and at the sub-regional scale. To isolate the consequences of DMP from demographic effects of host population dynamics, we compared the model outcomes under sl-CSD conditions with those from a hypothetical scenario in which all diploid offspring develop normally as female. We expressed the costs due to DMP by taking the difference between the two scenarios with respect to extinction rate, total population size, and the number of persisting sub-regions.

(b) Simulation experiments
The results of the simulation model illustrate that DMP is costly for the hyperparasitoid, over and above the direct demographic costs of fluctuating butterfly abundance at high fluctuation amplitudes and at low and intermediate migration rate Though locally sex alleles are lost, overall allelic diversity is maintained in the whole Å land even with relatively high host population fluctuation, and extinction of the hyperparasitoid is unlikely (electronic supplementary material, figure S4A-b). Once the amplitude of butterfly population fluctuation and autocorrelation surpass what has been observed in the natural system, the costs of DMP increase greatly (figures 4 and 5). Fewer sex alleles are maintained, which increases diploid male production. Overall population size declines, resulting in fewer local sub-regions persisting (electronic supplementary material, figure S4A). This pattern is most prominent at low hyperparasitoid migration rate.
The importance of balancing selection in maintaining sex alleles in the landscape is illustrated by comparing the number of alleles under sl-CSD (dark grey lines, figure 5c) with the allele number when all diploids are females (light grey lines, figure 5c). In the latter case, genetic drift owing to demographic stochasticity causes a large loss of allelic diversity without balancing selection, as host population size increasingly fluctuates.

Discussion
We present a scenario in which large fluctuations in population dynamics of a herbivore due to environmental change propagate up through the food chain, affecting a 4th-trophic-level hyperparasitoid, both demographically and genetically. In addition, we show the importance of spatial population structure and genetic rescue among local  populations for persistence under rapid environmental change. During a period of increased fluctuation of the host metapopulation size, we detected a precipitous increase in production of diploid males by the hyperparasitoid, up to approximately 28%, with a decline in population size that is evocative of a DMV. Our simulation model suggests that, despite signs of a DMV, the hyperparasitoid may not be facing regional extinction because of genetic rescue, which depends on spatial population structure and strong balancing selection resulting in global maintenance of high allelic diversity at the hyperparasitoid sex locus.

(a) The burden of single locus complementary sex determination and countermeasures
sl-CSD can be a constraint for haplodiploid species because inbreeding causes the production of inviable or sterile diploid males [10]. Many Hymenoptera minimize the cost of sl-CSD by avoiding inbreeding behaviourally via natal dispersal, protandry, or mate choice [51 -53]. Other parasitoid species have evolved multilocus (ml)-CSD to reduce DMP [18,23] or evolved means to produce fertile diploid males [12,13]. Some parasitoids have developed completely different sex determination mechanisms (e.g. [54 -56]). The relatively frequent evolution of countermeasures to DMP suggests strong selection against it in the evolutionary history of these species. The hyperparasitoid in this system is dispersive, as shown by its wide distribution in the landscape and very low spatial genetic structure relative to the host parasitoid and the host butterfly [39]. Being dispersive should reduce inbreeding and reduce diploid male load [27,29]. Thus, it may be that under normal conditions this wasp minimizes the costs of sl-CSD through dispersal.
We observed a relatively high percentage of diploid males (about 9%) despite a large N e before the wasp declined (figures 2 and 3). This suggests that the wasp has sl-CSD, rather than ml-CSD [25]. On islands species often have reduced genetic diversity owing to isolation and founder effects [57][58][59]. Thus, we expect Hymenoptera with CSD on islands to have relatively high DMP [20]. For example, captive (9.5%) and natural island populations (12%) of the parasitoid Venturia canescens have higher DMP than mainland populations (4.3%) [20].
In order to detect the cost of DMP and find evidence of DMVs, data should be collected over time to show an increase in DMP along with decreasing population size or N e . In one of the few studies encompassing more than a single season Weis et al. [25] followed up on a study by de Boer et al. [23] of an introduced population of the parasitoid Cotesia rubecula. They show that, despite a population bottleneck upon introduction, DMP decreased from 8-13% to 0-3% over 10 years. Additionally, balancing selection maintained high variation at the CSD locus. Similarly, Gloag et al. [30] saw the restoration of allelic diversity at the CSD loci during the invasion into Australia of the Asian honeybee. Unlike these studies, we found an increase in the percentage of diploid males in the hyperparasitoid with a drastic decline in N e over the period of four years, indicating the possibility of an ongoing DMV.

(b) Host population fluctuations and diploid male vortices in the hyperparasitoid
The metapopulation size of the host butterfly in Å land has historically been relatively constant, although local populations fluctuate strongly [36]. The hyperparasitoid has tolerated these local fluctuations well via dispersive foraging  that allows it to find hosts even in newly colonized local butterfly populations [35,39]. Over the last decade, the butterfly metapopulation has experienced increased spatial synchrony of local population sizes and fluctuation of population size that is attributable to increased weather extremes [38]. As part of this trend, 2010-2011 had the lowest recorded metapopulation size over the last 20 years (figure 2a). The primary parasitoid maintained a constant rate of parasitism (about 1/3 of the butterfly population) over the four years of the study (figure 2b). By contrast, the relative population size of the hyperparasitoid declined ( figure 2c). In addition, its absolute neutral allelic diversity declined, while spatial genetic structure increased with decreasing population size [39]. In fact, the hyperparasitoid population did not recover in spite of the dramatic recovery of the butterfly and primary parasitoid populations in 2011 -2012 (figure 2b). These observations suggest that some sex alleles were lost or became very rare during the 2010-2011 crash. The sharp decline in the hyperparasitoid along with exceptionally high frequencies of diploid males is in accordance with the theoretical predictions of the DMV [26].
Small populations are at severe risk of extinction simply because of demographic and environmental stochasticity [60] as well as genetic factors [8,61]. Our model results elucidate the cost of DMP in the hyperparasitoid separate from demographic extinction risks. To illustrate the potential effects of climate changes, we gradually increased spatial synchrony and variability of fluctuations in local butterfly population size in the model from completely asynchronous with minimal fluctuation to highly synchronous with large fluctuation (figure 5a). Under high fluctuation of the butterfly population, the hyperparasitoid experiences increased local extinctions associated with increased DMP, hence a DMV, in small-to medium-sized sub-regions of the landscape (figures 3 and 5). Regional extinction risk owing to a DMV increase only at fluctuation rates well beyond what has been observed in Å land ( figure 4). At the current level of fluctuation (years 2010-2017), the hyperparasitoid population is unlikely to be facing an imminent extinction from the entire Å land islands. In fact, the hyperparasitoid still persists now in 2018. On the other hand, at the sub-region scale the number of sex alleles declines, the fraction of diploid males increases, and total population sizes decrease even under the current conditions, which could explain why the hyperparasitoid population continued to decline despite the recovery of the butterfly.

(c) Genetic rescue and balancing selection
Balancing selection can reduce the proportion of diploid males in a population over time by maintaining high allelic diversity at CSD loci [25,30]. Dispersal of individuals between sub-regions can rescue local populations from extinction both demographically and genetically [62,63]. Dispersal together with balancing selection in a spatially structured population may be sufficient to prevent a DMV [27]. Genetic rescue is particularly effective when unique sex alleles are maintained in different parts of the landscape. In our model local hyperparasitoid populations collectively maintain all 10 sex alleles, and local populations in larger sub-regions are able to maintain a large fraction of them when conditions are not extreme (figure 5; electronic supplementary material, figure S4). Hence, our study generalizes the earlier finding by Hein et al. [27] to highly fluctuating host and parasitoid populations. The effectiveness of this process is illustrated by the sharp contrast between the maintenance of sex alleles and rapid loss and fixation of a neutral allele through genetic drift under the same conditions (figure 5c; electronic supplementary material, figure S4A-b, S4B-b). The reduction of population decline caused by inbreeding via maintenance of allelic diversity through balancing selection is not confined to Hymenoptera with CSD. For instance, it is also found in plant species that have the gametophytic self-incompatibility S-locus [64,65].
Zayed & Packer [26] showed that theoretically sl-CSD increases extinction rate of haplodiploid populations leading to a DMV. Later theoretical studies relaxed the assumptions that populations are small and isolated, but most of the models assume resource (host) availability to the parasitoid to be temporally constant [25 -27,29]. Bompard et al. [28] made a non-spatial host-parasitoid model that allowed for temporal fluctuation of host density arising from the explicit host -parasitoid interaction. They found DMP to stabilize rather than destabilize the parasitoid population by dampening endogenous cyclic fluctuations, which made DMVs less likely. Our model includes both temporal and spatial fluctuation of host availability caused by exogenous (i.e. environmental) fluctuations and spatial heterogeneity. In addition, we show that the occurrence of DMV-driven extinctions depends on the spatial scale: DMP elevates extinction risk in small local populations via a DMV, but the risk is substantially alleviated in a spatially structured population via genetic rescue and strong balancing selection. As many populations are spatially structured in nature [66], DMVs appear not to be as likely in natural populations as it was proposed by Zayed & Packer [26]. However, once a population enters a DMV initiated by climate change, habitat fragmentation, or pesticide use (in the case of bees), the trajectory towards extinction might proceed quickly. Monitoring for harbingers of a DMV, such as increasing proportions of diploid males [67], may be a good monitoring strategy, as larger environmental fluctuations are anticipated under climate change and there are many ecologically and economically important species of Hymenoptera.
Ethics. There are no ethical issues associated with species used in this research.