High hymenopteran parasitoid infestation rates in Czech populations of the Euphydryas aurinia butterfly inferred using a new molecular marker

in Czech populations


Introduction
Hymenopteran parasitoids are one of the most diverse groups of animals in terrestrial ecosystems and play a key role in the natural regulation of their host populations (La Salle and Gauld 1991;Forbes et al. 2018).The impact of parasitoids on their hosts can vary depending on their ecological specialisation, but in general they are known to cause significant levels of mortality in their hosts (Hawkins 1994).High levels of parasitism may also pose a potential threat to many threatened butterfly species, especially to various specialists in fragmented landscapes (cf.Anton et al. 2007).Species of the genus Euphydryas Scudder, 1872 (Nymphalidae: Melitaeini) represent a suitable system for studying host-parasitoid interactions, as egg clutches and webs with gregarious caterpillars can easily be detected in the field (Stamp 1981;Hula et al. 2004;Johansson et al. 2019).Previous studies using rearing have shown that the level of parasitism varied annually and depended mainly on weather conditions, which play a key role in the synchronisation between larval development and the emergence of parasitoids (Porter 1979(Porter , 1983)).
The Marsh Fritillary, Euphydryas aurinia (Rottemburg, 1775) is an EU-protected butterfly, declining in many European countries (van Swaay et al. 2010), including the Czech Republic (Hejda et al. 2017).It belongs to a genetically polymorphic group of closely related taxa, the "E.aurinia complex" (cf.Korb et al. 2016), with a wide Palearctic distribution and regional habitat and host plant specificity (e.g., Munguira et al. 1997;Singer et al. 2002;Junker et al. 2010;Korb et al. 2016).In Central and Western Europe, its main habitats are oligotrophic grasslands, and the most frequently used host plant is Succisa pratensis Moench (Dipsacaceae) (Warren 1994;Anthes et al. 2003;Konvička et al. 2003;Meister et al. 2015).The butterfly is monovoltine, with flight period in late spring/early summer when the mated females oviposit on leaf rosettes of the host plant.The larvae feed gregariously in silken webs on the plants until overwintering.They enter hibernation with the host plants' senescence in mid-September and resume feeding solitarily in April.In the Czech Republic, the distribution is restricted to the western part of the country (Fig. 1), where it forms three distinct metapopulation clusters inhabiting ≈90 separate oligotrophic meadow patches interconnected by dispersal (Zimmermann et al. 2011;Junker et al. 2021;Tájek et al. 2023).This system is monitored annually by counting larval nests (cf.Ojanen et al. 2013) and displays remarkable within-site and inter-annual dynamics with booms and bursts (John et al. in rev.).
Diverse methods were used so far to study the parasitoids of E. aurinia, and related butterflies, ranging from field counts of hymenopteran cocoons (Ford and Ford 1930;Porter 1983), captive rearing (Eliasson and Shaw 2003), field experiments with captive-reared material (Stamp 1981) to population genetic studies targeting parasitoid adults (Lei and Hanski 1997;Van Nouhuys and Lei 2004).However, the question pivotal to the butterfly population dynamics and conservation, that of infestation rates relative to population cycle and state of the butterfly colonies, seems to be little explored.This is probably due to the work requirements for rearing both butterflies and parasitoids (cf.Klapwijk and Lewis 2014), combined with destructivity of such methods for field populations.To quantify the parasitism rates in the Czech Republic populations of E. aurinia, we developed a molecular method, allowing rapid and lowcost detection of Hymenoptera parasitoids' incidence.
DNA-based methods are increasingly used in studies of parasitoid-hosts interactions (Zhu et al. 2019;Jeffs et al. 2021).The protocols so far developed for Lepidoptera/Hymenoptera systems mainly focused on COI locus, commonly referred as barcode (Folmer et al. 1994;Hebert et al. 2004), with the hope that the solid barcoding databases will assist species` identification (Toro-Delgado et al. 2022).Particularly good results were obtained via a reversal approach when the adult parasitoids were screened for host DNA shortly after their emergence (Rougerie et al. 2011) or in a species-poor natural system (high Arctic: Wirta et al. 2014).However, use of COI-based primers may be unreliable without subsequent sequencing, because deeply phylogenetically conserved bases are few and too far between in COI to place a group specific primer.Therefore, it was necessary to find a novel primer or primer pair which would amplify Hymenoptera but not Lepidoptera in the mixed samples containing the DNA of known lepidopteran host and unidentified hymenopteran parasitoids.We found such a potential primer in the nuclear region encoding the 28S ribosomal DNA.The novel primer, together with a primer published by Larsen (1992), targets part of the 28S gene and aims to amplify only Hymenoptera.
In this paper, we quantify Hymenoptera parasitoids infestation rates in a selection of the Czech Republic populations of Euphydryas aurinia and relate the infestation rates to the stage of the butterfly population cycle.Additionally, we document utility of our primers' combination for rapid Hymenoptera infestation assessment in butterflies.

Material and methods
We sampled E. aurinia caterpillars in western Bohemia (Fig. 1) in late August and early September.We sampled two caterpillars per larval web (105 webs in total) from 13 sites in 2019 and four per web (90 webs in total) from 9 sites in 2020.
While sampling the caterpillars, we recorded the following: Julian date, to account for infestation changes during larval period; longest and shortest dimension of the larval web (cm); sward height, i.e., visually estimated height of surrounding vegetation in 2.5 metre radius circles around each larval web sampled; host plant density, expressed as the number of Succisa flowerheads in the circle; and webs density, expressed as the number of larval webs in the circle.
The material was stored in 96% ethanol, the DNA was extracted using the Tissue Genomic DNA Mini Kit (GenAid Biotech, Taiwan) following the manufacturer's protocol.
We targeted a part of the 28S D1 region.We used primer 28SD1F (GGG-GAGGAAAAGAAACTAAC; Larsen et al. 1992) in combination with a new primer HymR157 (TGGCCCCATTCAAGATGG) with a resulting product of 164-167 bp.For the primer design, we assembled a library of target sequences (Hymenoptera parasitoids) and non-target sequences (Lepidoptera) from sequences available in GenBank (primarily PopSet 300390962, Heraty et al. 2011).We aligned the sequences in GENEIOUS PRIME 2020.2.4 (https://www.geneious.com)software and used AMPLICON software (Jarman 2003) to identify sections with concentrat-ed nucleotides that consistently differ between the target and non-target groups.We then manually screened these promising sections and used general rules of thumb to design candidate primers.We aimed for at least three differences between target and non-target groups in the first five positions at the 3' end of the primer for reliable specificity.We then tested the candidate primer in an in-silico PCR in GENEIOUS and optimized melting temperature in the primer pairs by extending or shortening the primers.
To test the utility of the primer used, we also carried out control reactions which contained DNA extracted from various adult hymenopteran parasitoids and butterflies (Fig. 2) to assure that we amplified only the potential parasitoids and not the butterfly.Some of the obtained PCR products (n=4) of parasitoids from positive E. aurinia samples were sequenced in SEQme (Czech Republic) to confirm hymenopteran origin.We checked the identity of the obtained sequences by using BLAST (nBLAST algorithm) (https://blast.ncbi.nlm.nih.gov,Altschul et al. 1990).
Positive results were recalculated to infestation levels per larval web (1/0 factor, infested or not) and per site (% of larvae sampled).To relate the per-web infestation level to larval web properties, we carried out logistic regressions (binomial error distribution; in R 4.2.3,R Core Team 2018) with infestation (1/0) as the dependent variable; and the longest and shortest web ground projections, surrounding vegetation height, Succisa density, and number of larval webs as predictors.We used the information theory approach, comparing the fitted regression Akaike information criteria (AIC) with AIC of the null model, y~+1, and considered models with ΔAIC > ≈2.0 as fitting the data.
To relate the per site percentual infestation to larval counts at the sites, we used data from annual monitoring of the sites, ongoing since 2001 (Hula et al. 2004).Over this time, larval counts were obtained for roughly ¾ of site x year combinations (John et al. in rev.,Suppl. material 1).

Results
Out of 210 (year 2019) and 358 (year 2020) E. aurinia caterpillars assayed for hymenopteran DNA, we obtained 70 and 144 positive results, respectively; i.e., the total infestation rates were 33.3% and 40.2% per individuum.On a per-site basis, this translates to mean±SD / median / range 38.5±29.89/ 40 / 0-100 per cent in 2019, and 40.1±26.51/ 50 / 0-77.5 per cent in 2020.The sequences of the positive samples were 124-126 bp long.The most similar sequences in GenBank according to nblast algorithm are those of Cotesia glomerata (Linnaeus, 1758), with the query identity 95.2-96.03%.The next similar sequences did not even reach a match of 93%.
According to the logistic regressions, none of the recorded properties of larval webs were related to infestation of the web (Table 1).
At the level of individual colonies, the infestation rates were highly variable (Fig. 3a).The per-site infestation levels did not correlate with larval web counts from The fitted single-term models are compared with the null model(s) following the information theory approach.None of the models was significant, as ΔAIC (null -fitted model) were always < 2.0. a) We fitted two null models, because measurements of E. aurinia larval webs were not available for 34 webs, which were disintegrated at the time of sampling the caterpillars.the same year (Spearman's r s = 0.11, t (17df ) = 0.46, p = 0.65) or with web counts in the subsequent year (r s = 0.31, t (16df ) = 1.32, p = 0.21), but did correlate positively with the larval webs' counts from the previous year (i.e., 2018 web counts for 2019 sampling, and 2019 web counts for 2020 sampling: r s = 0.46, t (16df ) = 2.07, p = 0.054) (Fig. 3b).
Because the absolute values of web counts highly varied among the sites (Fig. 3a), depending, e.g., on the site areas, we also recalculated web counts into percentage fractions of a ten-year (2011-2021) maximum for the given site, and recalculated the correlations, with no significant results (identical year percentage web counts: r s = -0.23,t (17df ) = -0.99,p = 0.37; previous year percentage web counts: r s = 0.31, t (16df ) = 1.31, p = 0.21, subsequent year percentage web counts: r s = 0.10, t (16df ) = 0.40, p = 0.69).

Discussion
The novel combination of primers 28SD1F (Larsen 1992) and HymR157 allowed us labour-efficient detection of high infestation rate in the declining E. aurinia butterfly by Hymenoptera parasitoids.Sequencing a selection of the obtained PCR products suggested that some of the parasitoids belong to the genus Cotesia Cameron, 1891.This is supported by the fact that other hymenopteran parasitoids known from our region, Ichneumon spp.and Pteromalus spp., attack pre-pupation larvae and pupae, respectively, whereas we worked with pre-diapause larvae.The gene 28S D1 is currently little represented in nucleotide databases for genus Cotesia, which, together with the unsettled taxonomy of Cotesia wasps infecting Melitaeinae butterflies (cf.Kankare et al. 2005a, b), precludes specific identification at this moment.This highlights the need to continue building comprehensive reference libraries for species identification.Such libraries are currently well developed for the COI gene, but supplementing it with a marker in another locus could improve discrimination power in some complex cases.Still, our approach allowed quantifying hymenopteran infestation rates in the declining butterfly, revealing that per-colony infestation rate is affected by larval webs density in the previous season.
Although the role of parasitoids on population fluctuations of E. aurinia, and related species, had been proposed almost a century ago (Ford and Ford 1930), relatively few authors quantified the natural infestation rates, with widely varying results.Infestation rates of <5% were reported, e.g., for the American congeneric species E. editha (Boisduval, 1852) (Singer and Erlich 1979) and E. chalcedona (Doubleday, 1847) (Lincoln et al. 1982).Similar results were obtained for E. aurinia from Sweden, with rates 2.6% (Eliasson and Shaw 2003), and Spain, where the rates were 2.4-5.1% (Stefanescu et al. 2009).The latter study in fact covered a newly recognised species, E. beckeri (Lederer, 1853), feeding on Lonicera spp.The rates detected by us, 33.3% and 40.2% in two consecutive years, are more comparable to the situations reported for American E. phaeton (Drury, 1773) (up to ≈10% in larvae prior to diapause) (Stamp 1981), E. maturna (Linnaeus, 1758) in the Czech Republic (69%) (Dolek et al. 2006) and Sweden (32%) (Eliasson and Shaw 2003).For E. aurinia, values higher than ours (≈90%) were reported by Ford and Ford (1930) from Britain during peaks of cyclic fluctuation of the butterfly population, whereas Klapwijk and Lewis (2014) reported a high range of infestation rates within individual webs (4-83%) from Britain.This all points to a high variation among sites, seasons, parasitoids' generations, and Euphydryas species in hymenopteran infestation rates.In detailed studies of the closely related Melitaeinae model species, Melitaea cinxia (Linnaeus, 1758), this variation was attributed to spatial positions of butterfly colonies, competition among parasitoids and hyperparasitism, and annual variation in weather (e.g., Lei and Hanski 1997;Van Nouhuys and Lei 2004).Arguably, some of the reported variation may also be due to the diversity of methods applied by various authors, ranging from field counts of infested caterpillars (Lei and Hanski 1997), through captive rearing (Stamp 1981;Eliasson and Shaw 2003;Klapwijk and Lewis 2014), to molecular methods as used here.Possibly, the molecular detection reveals higher infestation rates than rearing, because some of the infested larvae may die prior to the parasitoid emergence.Causes of this mortality are then interpreted as "unknown" (e.g., in Eliasson and Shaw 2003).Klapwijk and Lewis (2014) observed that the probability of caterpillar web infestation increased in webs isolated from other E. aurinia larval webs.We did not detect any relationship to the webs or surrounding vegetation parameters, but these results may be biased, because -for conservation concerns -we sampled the caterpillars solely from colonies containing a high number of larval webs during the sampling.The mean±SD webs' number of sampled colonies were 116±97.3and 60±66.6,whereas the numbers across all colonies were 53±57.8(n = 44) and 19±38.0(n = 70) in 2019 and 2020, respectively (John et al. in rev.).Conservation concerns also prevented quantifying the proportion of infestation per larval web, which would require killing all the caterpillars (cf.Klapwijk and Lewis 2014).
It also should be noted that our approach did not distinguish parasitoids from hyperparasitoids; i.e., the insects developing within parasitoids and thus killing them (Nair et al. 2016).A hyperparasitoid, however, can infest only a larva already infested by a parasitoid, and hence hyperparasitoids presence does not affect our findings on hymenopteran infestation rates.
With all the limitations, we found that the infestation rate per site positively correlated with per-site caterpillar webs' numbers of the previous year.This is fully expected if the parasitoids need a rich resource supply (i.e., high host density) to multiply in a butterfly colony, depleting the hosts' numbers in the process (Ford and Ford 1930;Frazer 1954;Porter 1981Porter , 1983)).A time delay in parasitoids infestation, and higher likelihood of infestation of larger and more connected host colonies, were found by Lei and Hanski (1997) in the metapopulation system of M. cinxia and its parasitoids.Although the inter-annual abundance changes of Melitaeini butterflies' colonies are likely influenced by numerous other factors, including variation in weather (Brunbjerg et al. 2017) or site vegetation management (Johansson et al. 2019;Tájek et al. 2023), natural enemies' pressure certainly plays a significant role.

Figure 1 .
Figure 1.The distribution of Euphydryas aurinia in the Czech Republic, historic records included, based on Beneš et al. (2002), with actualisations.The inset in upper right corner shows the position of the country in Europe.

Figure 2 .
Figure 2. Electrophoresis gels used to assess whether the primers used can discriminate lepidopteran hosts and hymenopteran parasitoids a various adult hymenopteran parasitoids (PCRs are positive) and positive (P) and negative (N) samples of E. aurinia b four species of adult butterflies; PCRs are negative.The adult specimens of Hymenoptera and Lepidoptera were identified by M. Rindoš, M. Konvička, and Z. Faltýnek Fric.

Figure 3 .
Figure 3. Per-site hymenopteran parasitoids infestation rates in colonies of the butterfly Euphydryas aurinia in two consecutive years (2019-20), with information of caterpillar web counts in the respective colonies in 2018-2021 (above), and illustration of the relationship between hymenopteran parasitoids infestation rates and E. aurinia caterpillar web counts in the previous year (below).

Table 1 .
Logistic regressions relating field-measured properties of larval webs to infestation of the sampled larvae.