Prevalence and population genetic analyses of parasites in invasive Vespa velutina and native Hymenoptera

Prevalence and population genetic analyses of parasites in invasive Vespa velutina and native Hymenoptera Carolina Bartolom´ e a,b,* , Damian Dasilva-Martíns a,b , Rosa Vali˜ nas a,b , Luís B. Gabín-García a,b , Anabela Nave c,d , Ana L. García-P´ erez e , Karine Monceau f , Denis Thi´ ery g , Alastair Christie h , Moon Bo Choi i , Beatriz Sobrino b,j , Jorge Amigo b,j , Xulio Maside a,b a Grupo de Medicina Xen´ omica, CIMUS, Universidade de Santiago de Compostela, 15782 Santiago de Compostela, Galiza, Spain b Instituto de Investigaci´ on Sanitaria de Santiago (IDIS), 15706 Santiago de Compostela, Galiza, Spain c Centro de Investigaç˜ ao e de Tecnologias Agro-Ambientais e Biol´ ogicas, CITAB, Universidade de Tr´ as-os-Montes e Alto Douro, UTAD, 5000-801 Vila Real, Portugal d Unidade Estrat´ egica de Sistemas Agr´ arios e Florestais e Sanidade Vegetal (UESAFSV) Instituto Nacional de Investigaç˜ ao Agr´ aria e Veterin´ aria, I.P. (INIAV), Av. da República, Quinta do Marquˆ es, 2780-157 Oeiras, Portugal e Department of Animal Health, NEIKER-Basque Institute for Agricultural Research and Development, Basque Research and Technology Alliance (BRTA), 48160 Derio, Bizkaia, Spain f UMR 7372, Centre d ′ ´ Etudes Biologiques de Chiz´ e, La Rochelle Universit´ e & CNRS, 79360 Villiers en Bois, France g INRAE, UMR1065 SAVE, 33140, Villenave d’Ornon, France h Government of Jersey, Natural Environment, Route de la Trinit´ e, Trinity, JE3 5JP Jersey, Channel Islands i Institute of Agricultural Science and Technology, Kyungpook National University, 41566 Daegu, Republic of Korea j Fundaci´ on Pública Galega de Medicina Xen´ omica, Servizo Galego de Saúde (SERGAS), 15706 Santiago de Compostela, Galiza, Spain ARTICLE INFO Keywords: Vespa velutina Invasive species Native Hymenoptera Nosematidae Trypanosomatidae Lipotrophidae ABSTRACT Invasive species pose a threat to the ecological balance of the ecosystems they invade by altering local hostpathogen dynamics. To investigate these relationships and their potential consequences, we examined the prevalence and genetic diversity patterns of Trypanosomatidae, Lipotrophidae, and Nosematidae in a collection of sympatric isolates of the invasive hornet Vespa velutina and local Hymenoptera from two recently colonized areas: Europe and South Korea. Data were gathered through PCR amplification and massive parallel sequencing, and analyses were conducted using population genetics tools. Parasite prevalences showed substantial variation depending on (i) the parasite family (Trypanosomatidae and Nosematidae were the most and less prevalent, respectively), (ii) location (e.g. Galicia displayed the highest pooled values), (iii) the season (highest in spring for Trypanosomatidae and Lipotrophidae), and (iv) the host. V. velutina exhibited significantly lower parasite occurrence than native Hymenoptera across all parasite families (consistent with the enemy release hypothesis), although this difference was less pronounced during the periods of heightened predatory activity, suggestive of trophic transmission. Parasite species displayed significant genetic differentiation between European and South Korean isolates, yet no differentiation was observed across hosts, suggesting that all Hymenoptera are exposed to a common local pathogen population. There was no indication that V. velutina acted as a carrier of foreign parasites to the invaded territories. 1. Introduction The introduction of invasive species has an prominent effect on the health of the ecosystems, especially if their populations are large and widely spread, both by causing losses in biodiversity and by increasing the risk of disease due to the disruption of local host-pathogen dynamics (Chinchio et al., 2020). Such disturbances may occur through the transmission of non-indigenous parasites to native communities, where they tend to cause more virulent infections (Daszak et al., 2000; Wei et al., 2022). Additionally, invasive species can act as novel hosts or * Corresponding author at: Grupo de Medicina Xen´ omica, CIMUS, Universidade de Santiago de Compostela, 15782 Santiago de Compostela, Galiza, Spain. E-mail addresses: [email protected] (C. Bartolom´ e), [email protected] (D. Dasilva-Martíns), [email protected] (A. Nave), agarcia@neiker. eus (A.L. García-P´ erez), [email protected] (K. Monceau), [email protected] (D. Thi´ ery), [email protected] (A. Christie), [email protected] (M.B. Choi), [email protected] (B. Sobrino), [email protected] (J. Amigo), [email protected] (X. Maside). Contents lists available at ScienceDirect Journal of Invertebrate Pathology journal homepage: www.elsevier.com/locate/jip https://doi.org/10.1016/j.jip.2024.108203 Received 22 May 2024; Received in revised form 5 September 2024; Accepted 18 September 2024 Journal of Invertebrate Pathology 207 (2024) 108203 Available online 21 September 2024 0022-2011/© 2024 The Author(s). Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND license ( http://creativecommons.org/licenses/bync-nd/4.0/ ). vectors for native pathogens, further amplifying the disease burden in local populations (Chinchio et al., 2020; Kelly et al., 2009). Both mechanisms have been documented in pollinators (Cilia et al., 2022; Nanetti et al., 2021; Pereira et al., 2021; Pislak Ocepek et al., 2021) linked to the share of resources, which provide an ideal interface for the transmission of pathogens (Graystock et al., 2015; Ke et al., 2022). Vespa velutina nigrithorax, also known as yellow-legged hornet or Asian hornet, is an invasive species native of South-East Asia that was accidentally introduced to South Korea (2003) and southwestern France (2004) and has now colonised many areas of Korea, Japan and Europe (France, Spain, Portugal, Belgium, the Netherlands, Luxembourg, Switzerland, Italy and Germany, with a few incursions in other regions like Hungary, Czech Republic, the UK or Ireland (Lioy et al., 2022; Inventaire National du Patrimoine Naturel). In August 2023, its presence was also reported in Georgia (USA) (https://agr.georgia.gov/yello w-legged-hornet). This hornet represents an important threat to insect populations, especially to honey bees, which constitute a large proportion of its diet (Rome et al., 2021). The predation pressure exerted by this species leads to changes in foraging habits and has a serious impact on pollination and beekeeping (Requier et al., 2019; Rojas-Nossa and Calvi˜ no-Cancela, 2020). Besides, its presence generates severe economic costs to beekeepers, agricultural producers and administrations (Barbet-Massin et al., 2020), as well as important safety concerns (Fe´ as, 2021). In addition, V. velutina poses a serious threat to the health of other sympatric species, since its abundance and feeding habits (involving predation and plant visiting (Monceau et al., 2014; Rojas-Nossa and Calvi˜ no-Cancela, 2020; Rome et al., 2021)) meet the conditions to potentially alter the prevalence and diversity of local parasite communities. However, despite this risk, the assessment of the relationship between the parasites harboured by this hornet and the native entomofauna has only been addressed in Galicia (SW-Europe) by applying a PCR-cloning and Sanger sequencing strategy (Gabín-García et al., 2021). The progress of the technology, namely the use of massive parallel sequencing (MPS), which enables the simultaneous identification of multiple pathogens in a large number of samples (Kulski, 2016), and the possibility of analysing specimens, both of V. velutina and native Hymenoptera, collected in different and distant invaded areas (specifically, Europe and South Korea), have prompted this work, whose main objective is to investigate the potential transmission of parasites between invasive and local hosts and its impact on these communities. 2. Material and Methods 2.1. Sample processing and DNA extraction One hundred and ten V. velutina specimens and 134 native Hymenoptera from different Atlantic European regions, as well as 240 hornet samples from South Korea (150 V. velutina and 90 local Vespidae), were obtained to provide initial data on parasite prevalence (Table 1 and Supplemental Table 1). These samples were collected either using traps (i.e., liquid baits, electric harps, entomological nets) or directly from nests, and submitted to CiMUS in 80 % Ethanol where they were stored at 4 ◦C until further processing. Specimens were dissected to remove their digestive tract, including midgut, Malpighian tubules, and fat body. To prevent any external contamination, prior to any manipulation, samples were washed individually with sterile distilled water for 1 min, transferred to a fresh 5 % sodium hypochlorite solution for 10–15 s, and finally rinsed three times with sterile distilled water. All dissection tools were disinfected with 5 % bleach for 20 min, rinsed twice in sterile distilled water and once in 80 % Ethanol, and air-dried prior to handling each specimen. The abdominal section of each specimen’s gut (i.e., ventriculus, ileum and rectum) was extracted by pulling from both ends of the abdomen with the aid of two forceps, transferred to a 1.5 μ L plastic tube and homogenized with 200 μ L of distilled water using a sterile disposable plastic pestle (VWR). DNA was extracted from these homogenates by applying a phenol: chloroform: isoamyl alcohol protocol. Briefly, 100 μ L of homogenate were digested overnight at 56 ◦C in 300 μ L of lysis buffer (50 mM Tris-HCl pH 8, 100 mM EDTA, 100 mM NaCl, 1 % SDS) and 5 μ L of proteinase K (~20 mg/mL, Thermo Scientific). Afterwards, lysates were treated with 400 μ L of phenol: chloroform: isoamyl alcohol (25:24:1; Sigma Aldrich) and the supernatants – obtained after a 6 min centrifugation at 12,000 rpm –transferred to fresh vials where the DNA was precipitated with 400 μ L of isopropanol. DNA pellets, resulting from a new 6 min centrifugation at 14,000 rpm, were washed with 80 % Ethanol, air-dried, and resuspended in 100 μ L of PCR grade water (Invitrogen). They were stored at −20 ◦C until further use. 2.2. PCR amplification The marker loci Actin, RPB1 (RNA polymerase II large subunit RPB1) and SSU (small-subunit ribosomal DNA) were used for the detection of Nosematidae, Trypanosomatidae and Lipotrophidae, respectively (Bartolom´ e et al., 2020), three major families of hymenopteran enteropathogens that are particularly prevalent in bees (Evans and Schwarz, 2011) and especially in honey bees, which are one the main preys of these invasive hornets. PCR amplifications were carried out using the Phusion High-Fidelity DNA Polymerase (Thermo Fisher Scientific). The reactions were performed in 25 μ L volumes containing 13.3 μ L of PCR grade water (Invitrogen), 5 μ L of 5X Buffer, 0.5 μ L of 10 mM dNTPs, 2.5 μ L of each primer pair (Supplemental Table 2), 0.3 μ L of Phusion DNA polymerase and 1 μ L of DNA template set to a concentration of 40 ng/ μ L. Cycling conditions consisted of an initial denaturalization at 98 ◦C for 30 s, followed by 45 cycles of 98 ◦C for 10 s, 60–66 ◦C for 30 s (see Supplemental Table 2 for priming details) and 72 ◦C for 10 s, and a final extension of 8 min at 72 ◦C. Negative controls were included in all reactions to monitor for any potential contamination. PCR results were checked by electrophoresis in 2.0 % agarose gels. Positive amplicons were purified with Agentcourt AMPure XP (Beckman Coulter) and quantified with Qubit 2.0 (Thermo Fisher Scientific). 2.3. PCR amplification, cloning and Sanger sequencing of new Nosematidae taxa Samples carrying unidentified Nosematidae taxa were subject to an additional PCR amplification to obtain sequences of the small subunit rRNA (SSU rRNA) locus. This was carried out using the conditions Table 1 Distribution of the samples analysed according to their origin. Species Portugal Galicia Basque Country France Jersey S. Korea Pooled Vespa velutina 25 25 25 25 10 150 260 Vespa crabro 5 16 15 16 −30 82 Vespa simillima − − − − − 30 30 Vespa analis − − − − − 30 30 Vespula germanica 2 9 10 10 10 −41 Vespula vulgaris 2 8 5 −5−20 Bombus spp. 7 10 4 − − − 21 All samples 41 68 59 51 25 240 484 C. Bartolom´ e et al. Journal of Invertebrate Pathology 207 (2024) 108203 2 described above and the primers detailed in Supplemental Table 2. Amplicons were gel extracted (NZYGelpure, NZYTech, Portugal) and cloned using CloneJET PCR Cloning Kit, Thermo Scientifc). Plasmid DNA was purified (NZYMiniprep, NZYTech, Portugal) and sequenced on an ABI 3730XL sequencing machine (GATC, Eurofins Genomics, Germany). 2.4. Massive parallel sequencing (MPS) 2.4.1. Library preparation and sequencing Among all PCR-positive samples, 182 were selected for MPS (78 from V. velutina and 104 from native Hymenoptera; Table 2) for an in-depth analysis of pathogen diversity. Sample selection was based on amplicon concentration, trying to maximize the representation of parasite families (52 samples were positive for more than one pathogen group; 26 were positive for Nosematidae, 155 for Trypanosomatidae and 44 for Lipotrophidae), host species and sampled regions. Purified amplicons from each sample were pooled at equimolar concentrations (50–100 ng/ μ L). Library preparation and sample indexing were carried out using the KAPA HyperPrep kit (Roche Sequencing Solutions Inc.) and the KAPA universal adapter together with the Unique Dual-Indexed primer mixes (Roche Sequencing Solutions Inc.), following the manufacturers protocols (https://rochesequencingstore. com/wp-content/uploads/2022/07/KAPA-HyperPrep-Kit-Technical- Data-Sheet.pdf). This allowed the indexation of up to 96 samples per run. Amplicons were used as input for library preparation directly from the end-repair and A-tailing step since their size (187 to 283 bp) did not require a fragmentation step. Library pools were normalized to a concentration of 4 nM and loaded at a concentration of 12 pM on an Illumina MiSeq instrument for 1 ×300 bp single-end sequencing (Flow cell Nano V2). 2.4.2. Sequence processing Raw reads were processed with fastp for quality filtering and sorted according to the sequences of the primers (which were discarded from further analyses). Once identified, they were divided into individual FASTA files providing the number and sequence of the haplotypes detected for each amplicon in a sample. Those containing less than five reads were excluded from the final dataset. Furthermore, only haplotypes present in at least two individuals were included in the analyses, except if identical sequences had already been deposited in a public database or if they were the only haplotype in a sample. These were submitted to GenBank, where they are available under accession numbers PP692852-PP693040, PP704714-PP706065, PP709277- PP709280, PP718138-PP718261 and PP725407-PP725408. 2.5. Evolutionary analyses The nucleotide diversity in Trypanosomatidae and Nosematidae species was estimated at synonymous and nonsynonymous sites −or just silent in the case of Lipotrophidae −using two statistics: π (the average number of nucleotide differences per site between two sequences; (Nei, 1987), which was calculated applying the Jukes and Cantor correction (Jukes and Cantor, 1969), and ϴ W (the number of segregating sites; Watterson, 1975). These parameters, as well as the Tajima’s D statistic (Tajima, 1989b), were obtained with DnaSP v6 (Rozas et al., 2017). The same software was used to quantify the synonymous divergence between species (K S ), which was estimated using the Nei-Gojobori model (Nei and Gojobori, 1986), applying the Jukes and Cantor correction (Jukes and Cantor, 1969). Phylogenetic analyses were conducted using MEGA v. 11 (Tamura et al., 2021). The best model for each locus was selected by applying the Akaike information criterion (AIC) and the reliability of the resulting tree topologies was tested by bootstrap support (1000 replicates). The evolutionary relationships among both Nosematidae and Trypanosomatidae were inferred using the Neighbour-Joining (NJ) method and evolutionary distances were estimated at all sites by applying the Tamura 3-parameter model. The rate of variation among sites was modelled with a gamma distribution (shape parameter =1). Primer sequences and positions containing gaps and/or missing data were eliminated (complete deletion option). The scrutiny of the genetic structure of the sampled populations was carried out using the analysis of the molecular variance (AMOVA) implemented in Arlequin 3.5.2.2 (Excoffier and Lischer, 2010). This method is similar to other approaches based on the variance of haplotypic frequencies, but it also combines information on the genetic distance between haplotypes. A hierarchical analysis of variance allows the partition of the total variance into covariance components corresponding to three levels of genetic structure: within-specimens, among specimens within groups and among groups of specimens. Their significance was checked by applying non-parametric permutation procedures (1000 permutations). The population structure was also investigated using the Snn statistic (Hudson, 2000) which estimates how often related sequences are found in the same population. It was calculated with the aid of DnaSP v6 (Rozas et al., 2017) and its significance assessed using permutation tests (1000 replicates). 2.6. Statistical analyses The hypergeometric test was used to determine whether the observed co-occurrence of parasite families adjusted to random expectations. Calculations were made with the aid of Hypergeometric Distribution Probability Calculator (Berman H.B., “Hypergeometric Distribution Probability Calculator”, Available at: https://stattrek. com/online-calculator/hypergeometric [accessed 20–5-2024]). Fit to the Poisson distribution was done by means of the Chi-square test of goodness of fit (Sokal and Rohlf, 1998). 3. Results 3.1. Parasite prevalence Two hundred fifty two out of the 484 specimens analysed by PCR (52.1 % of the dataset) produced PCR amplicons of the expected size for any of the three pathogen groups under screening. However, subsequent sequencing (see below) revealed that 16 of them yielded pollen, bacteria, or host-specific sequences, so they were discarded from any further Table 2 Number of samples used for massive parallel sequencing (MPS). Host species Portugal Galicia Basque Country France Jersey S. Korea Pooled Vespa velutina 9 14 3 10 2 40 78 Vespa crabro −5 9 3 −18 35 Vespa simillima − − − − − 17 17 Vespa analis − − − − − 15 15 Vespula germanica −5 9 3 1 −18 Vespula vulgaris −4 4 −2−10 Bombus spp. 2 5 2 − − − 9 All hosts 11 33 27 16 5 90 182 C. Bartolom´ e et al. Journal of Invertebrate Pathology 207 (2024) 108203 3 analyses, reducing the final number of positive samples to 236 (48.8 % of the dataset). Pairwise comparisons of the distributions of parasite families across samples revealed a significant excess of co-occurrence of the three pairwise combinations of families relative to random expectations (p ≪ 0.001 in all pairwise comparisons using a hypergeometric distribution; Material and Methods subsection 2.6). To compare the prevalences we used the relative parasite frequency (RPF), a parameter representing the number of parasite families present in each specimen relative to the maximum possible (N=3). Overall, parasite prevalences were consistently greater among European than S. Korean samples (RPF=31.4 % vs 15.0 %, respectively; Table 3 and Supplemental Tables 3 and 4). In Galicia (RPF=44.6 % ±6.03) and the Basque Country (RPF=33.3 ±6.14) this parameter was well above that of the pooled sample (23.3 % ±1.9). Trypanosomatidae was the most prevalent group (RPF=41.9 % ± 2.24; mean ±SE; from data in Supplemental Tables 3 and 4). It was particularly prominent in Galicia and the Basque Country, with regional pooled values above 62 %. Nosematidae stood out as the less frequent family with a pooled value of 9.1 % ±1.31, again with highest incidence in Galicia (26.5 % ±5.35). Parasite frequencies also varied across hosts (p <0.001 in a Chisquared test of independence). Overall, V. velutina’s RPF was about half that of native Hymenoptera (16.0 % ±2.28 vs. 31.7 % ±3.11, respectively; Table 3). This trend, although not always statistically significant, was observed across all regions except Portugal and Galicia. Given the likely influence of environmental factors on parasite prevalences and considering that the majority of specimens were collected during the summer (51.4 % compared to 19.4 % in spring and 29.1 % in autumn; see Supplemental Table 1), a seasonal analysis was conducted to assess potential differences in the relative fractions of positive samples in spring, summer, and autumn. The assessment of the frequencies in each of these periods revealed that the higher RPF observed in native Hymenoptera with respect to V. velutina extended to all pathogen families, although the differences between hosts faded during summer (Fig. 1). This analysis also showed that the occurrence of Trypanosomatidae and Lipotrophidae tended to decrease from spring to autumn (Fig. 1), while Nosematidae were more prevalent in summer. 3.2. Molecular identification of the parasite species by MPS One hundred and eighty-two specimens of the 252 PCR-positive samples were selected (Material and Methods, subsection 2.4.1) for an in-depth assessment of patterns of pathogen genetic diversity. Of those, 166 (91.2 %) produced parasite-related sequences (Supplemental Table 4). MPS yielded 3,246,955 raw reads which, after filtering, resulted in 602,431 high quality sequence reads for Trypanosomatidae, 106,720 for Lipotrophidae, and 163,073 for Nosematidae. These figures correspond to an average of 3886.7, 2475.8, and 4801.4 full-length sequences per sample for each of these parasite families, respectively. 3.2.1. Trypanosomatidae Crithidia mellificae and Lotmaria passim were identified as the predominant species within a sample of 155 Trypanosomatidae-positive specimens (77.4 % and 67.1 %, respectively), followed by C. bombi (27.7 %) and a recently identified taxon (Crithidia species 4, hereinafter abbreviated as C. sp. 4; 23.9 %) (Bartolom´ e et al., 2020; Gabín-García et al., 2021). C. acanthocephali and two other taxa (Crithidia species 2 and 5, hereinafter abbreviated as C. sp. 2 and C. sp. 5, respectively), were present at much lower frequencies (<2.0 %; Table 4 and Supplemental Table 4). Two of the samples exhibited recombinant sequences between L. passim and C. sp. 4 and were discarded from further analyses (Supplemental Table 4). Most specimens (67,3%) harboured more than one Trypanosomatidae species. The distribution of the number of different Trypanosomatidae taxa per sample did not deviate from random expectations assuming a Poisson distribution (Material and Methods, subsection 2.6). No pairwise associations between the most frequent species were detected (L. passim, C. mellificae, C. bombi and C. sp.4; p >0.05 in pairwise comparisons using the hypergeometric distribution). The relative frequencies of parasite species varied significantly across continents: C. mellificae was more frequent among Asian than European positive specimens (86.5 % and 69.1 %, respectively; p <0.001 in a twotailed Z-test). Conversely, L. passim and C. sp. 4 were more frequent among European than South Korean samples (86.4 % vs 45.9 % and 39.5 % vs 6.8 %, respectively; p <0.001 in two-tailed Z-tests). A phylogenetic tree constructed with previously published Trypanosomatidae sequences (Fig. 2) showed that the newly obtained sequences of C. sp. 4 clustered with others located within the R3 clade (Fig. 2 of Gabín-García et al., 2021) which displayed little genetic differentiation (K S between 0 and 1.1 %). In contrast, the sequences assigned to R1 (Fig. 2) exhibited considerable levels of synonymous pairwise divergence (Ks up to 13–15 %) which suggests that, although they were initially gathered into the same group under the name “Crithidia species 2” (GenBank accession numbers MW288790.1- MW288792.1), they likely represent two distinct taxa: Crithidia sp. 2 (accession numbers MW288792.1 and PP693038.1) and a new one called Crithidia sp. 5 (accession numbers MW288790.1, MW288791.1 and PP693040.1). However, it must be noted that this classification and the taxonomic relationships among the different taxa should be taken with care due to the short length of the sequences analysed. The most polymorphic members of the Trypanosomatidae family at synonymous sites were C. mellificae (pooled π S =4.4 %) and C. acanthocephali (pooled π S =5.3 %) (Supplemental Table 5). L. passim and C. mellificae were the only species showing consistent differences in polymorphism levels between Europe and South Korea. Specifically, the pooled π S of European samples were about two-fold those of South Korean ones for both species (2.9 % vs. 1.3 % for L. passim and 5.2 % vs. 3.1 % for C. mellificae, respectively). This pattern applied to both V. velutina and native hosts, e.g. π S for L. passim in European V. velutina and native Hymenoptera were 3.0 % and 2.9 %, respectively, compared to 1.5 % and 1.2 % in South Korea. Such differences were not detected at nonsynonymous positions, which were always less variable than synonymous sites (Supplemental Table 5). Table 3 Relative parasite frequency (RPF) across hosts and regions (average % ±standard error). Host species Portugal Galicia Basque Country France Jersey S. Korea Pooled V. velutina 18.7 ±7.79 57.3 ±9.89 12.0 ±6.50 20.0 ±8.00 6.7 ±7.89 9.3 ±2.38 16.0 ±2.28 V. crabro 13.3 ±15.20 41.7 ±12.33 46.7 ±12.88 35.4 ±11.96 −26.7 ±8.07 34.1 ±5.24 V. simillima − − − − − 27.8 ±8.18 27.8 ±8.18 V. analis − − − − − 18.9 ±7.15 18.9 ±7.15 V. germanica 33.3 ±33.33 40.7 ±16.38 53.3 ±15.78 36.7 ±15.24 6.7 ±7.89 34.1 ±7.41 V. vulgaris 33.3 ±33.33 20.8 ±14.36 46.7 ±22.31 −53.3 ±22.31 −36.7 ±10.78 Bombus spp. 23.8 ±16.10 40.0 ±15.49 50.0 ±25.00 − − − 36.5 ±10.51 All hosts 20.3 ±6.28 44.6 ±6.03 33.3 ±6.14 28.1 ±6.29 16.0 ±7.33 15.0 ±2.30 23.3 ±1.92 Native hosts 22.9 ±10.51 37.2 ±7.37 49.0 ±8.57 35.9 ±9.41 22.2 ±10.73 24.4 ±4.53 31.7 ±3.11 C. Bartolom´ e et al. Journal of Invertebrate Pathology 207 (2024) 108203 4 The AMOVA analysis also indicated significant genetic differentiation between the European and South Korean populations of Trypanosomatidae, accounting for up to 9 % of the total molecular variance (Table 5 and Supplemental Table 6). Furthermore, molecular variance for C. mellificae, L. passim and C. sp. 4 also displayed a significant level of geographical structuring among European locations, a pattern not observed in South Korea (Table 5). This finding was further supported by the Snn values calculated for each of these groups (Snn =0.60, 0.50 and 0.59 for C. mellificae, L. passim and C. sp.4, respectively, p <0.001 in all cases; Supplemental Table 6). In contrast, covariance across hosts explained a much smaller fraction of the observed variation, and was not statistically significant, albeit for C. mellificae in Europe, where it accounted for 6.2 % of the total variance, and for N. ceranae in the pooled sample (3.1 % of the total variance; Supplemental Table 7). 3.2.2. Lipotrophidae Apicystis bombi was present in 97.7 % of the Lipotrophidae positive samples (43 out of 44), 68.2 % of which (N=30) corresponded to native pollinators (Supplemental Table 4). A second Neogregarinorida species, Ng. sp., occurred in just two European specimens (Table 6). This taxon was first reported in Apis mellifera from central Iberian Peninsula, where it was detected in nearly 30 % of the sample (Bartolom´ e et al., 2020; Accession number MN031271.1). Neutral genetic variability within Lipotrophidae species was very low in both cases (0.7 % for A. bombi and 0.4 % for Ng. sp., respectively; Supplemental Table 8). No significant differences in π S values were observed across continents ( π S =0.7 % each), nor between V. velutina and the native hosts ( π S =0.7 % and 0.8 %, respectively). However, the AMOVA results and Snn estimates suggest that there is a moderate, but significant, differentiation (6.7 % of the total variance) between the European and Asian A. bombi populations (Table 5 and Supplemental Table 6). No local structuring was otherwise detected in Europe nor in S. Korea. Fig. 1. Relative frequency of positive samples (%) across seasons and hosts. Sp: spring, Su: summer; Au: autumn. Table 4 Trypanosomatidae species detected by MPS across regions and hosts. Region Host species N L. passim C. mellificae C. bombi C. acantocephali C. sp. 2 C. sp. 4 C. sp. 5 Portugal V. velutina 7 4 6 0 0 0 2 0 Bombus spp. 2 1 0 1 0 0 1 0 Galicia V. velutina 13 10 10 4 1 1 2 2 V. crabro 4 4 2 1 0 0 0 0 V. germanica 5 5 2 1 1 0 0 0 V. vulgaris 2 2 2 0 0 0 0 0 Bombus spp. 5 4 2 2 0 0 1 0 Basque Country V. velutina 3 3 3 1 0 0 0 0 V. crabro 9 9 8 1 0 0 5 0 V. germanica 9 9 9 4 0 0 6 0 V. vulgaris 4 4 4 0 0 0 4 0 Bombus spp. 2 1 2 0 0 0 0 0 France V. velutina 10 8 2 2 0 0 6 0 V. crabro 3 3 1 1 0 0 3 0 V. germanica 3 3 3 2 0 0 2 0 S. Korea V. velutina 28 11 26 8 1 0 2 0 V. crabro 16 7 13 6 0 0 1 0 V. analis 15 9 12 6 0 0 1 0 V. simillima 15 7 13 3 0 0 1 0 N, number of positive samples analysed. C. sp., Crithidia sp. C. Bartolom´ e et al. Journal of Invertebrate Pathology 207 (2024) 108203 5 3.2.3. Nosematidae Nosema ceranae (also known as Vairimorpha ceranae) was the most abundant Nosematidae species (84.6 % of the 26 Nosematidae-positive specimens; Table 7 and Supplemental Tables 3 and 4). N. thomsoni, N. bombi (also known as Vairimopha thomsoni and Vairimorpha bombi, respectively) and two other Nosema taxa (N. sp. and N. sp.1, see below) were also detected, although at much lower frequencies (7.7 %, each). Most Nosematidae were found in V. velutina (mainly from Galicia, where its prevalence was highest; see above), and in V. simillima from S. Korea (Supplemental Table 3). Four positive samples presented slightly divergent Actin haplotypes that likely correspond to two previously uncharacterized taxa, here named N. sp. and N. sp. 1 (Nosema species and Nosema species 1, respectively; Table 7). Given that the identification of Nosematidae taxa is usually based on the small subunit rRNA (SSU rRNA), sequences of this locus were obtained from the N. sp. and N. sp. 1 positive samples (Material and Methods subsection 2.3). A BLAST analysis showed that N. sp. had already been found in several Vespidae (V. velutina, V. crabro, Vespula and Polistes) (Gabín-García et al., 2021), Bombus (Li et al., 2012) and Lepidoptera (Tokarev et al., 2015) species, whereas N. sp. 1 was new to the public databases (this new taxon was also identified in V. velutina specimens from Galicia by means of PCR amplification, cloning and Fig. 2. Phylogenetic relationships of RPB1 sequences from Trypanosomatidae (using Trypanosoma cruzi as an outgroup). The evolutionary history was inferred using the Neighbour-Joining method and the distances were estimated by applying the Tamura 3-parameter model assuming a gamma distribution. Bootstrap values ≥50 % are shown. Sequences obtained in this study are highlighted in bold. The numbering of the clades (R1-R4) corresponds to that of Gabín-García et al. (2021). Sequences displaying an asterisk were previously classified as Crithidia species 2. Table 5 Variance components in parasite populations across geographical regions (AMOVA analyses). Regions Parasite Covariance Component EU-SK EU SK L. passim Among populations 4.33 *** 8.47 ** −0.79  Among specimens within populations 15.55 *** 14.34 *** 2.46  Within specimens 80.12 *** 77.19 *** 98.33        C. mellificae Among populations 5.08 *** 16.86 *** −2.40  Among specimens within populations 15.32 *** 2.44 * 22.04 *** Within specimens 79.59 *** 80.70 *** 80.36 ***       C. bombi Among populations 5.58 ** 7.88 7.1  Among specimens within populations 1.56 5.47 −15.25  Within specimens 92.85 86.65 108.15        C. sp. 4Among populations 8.78 * 4.29 * −11.64  Among specimens within populations 1.69 0.06 0.27  Within specimens 89.53 95.65 111.38        A. bombi Among populations 6.67 *** 4.20 10.67  Among specimens within populations −26.09 −36.70 −20.77  Within specimens 119.42 132.50 110.11        N. ceranae Among populations 11.55 *** 1.69 −0.47  Among specimens within populations −7.77 −13.53 −5.65  Within specimens 96.23 111.84 106.12  EU-SK: comparison between European and South Korean populations; EU: among European populations; SK: among South Korean populations; statistical significance was assessed using permutation tests (1000 replicates): *, p <0.05; **, p <0.01; ***, p <0.001. C. Bartolom´ e et al. Journal of Invertebrate Pathology 207 (2024) 108203 6 sequencing of the SSU rRNA locus; R. Vali˜ nas unpublished results). Phylogenetic trees constructed with sequences of these two loci revealed that N. sp. and N. sp. 1 were closely related to Nosema muscidifuracis (Becnel and Geden, 1994) (Fig. 3), although both topologies exhibited low bootstrap due to the short length of the sequences analysed. In terms of synonymous diversity ( π S ), N. bombi, N. sp. and N. thomsoni were monomorphic at the Actin locus whereas N. sp.1 and N. ceranae showed a pooled π S of 1.6 % and 1.2 %, respectively (Supplemental Table 9). Despite the similarity of these figures, an indepth analysis of polymorphism data revealed the existence of distinct patterns in these two species. While N. sp.1 displayed higher π S than π A (1.6 % and 0.9 %, respectively), N. ceranae exhibited similar overall values for both parameters (1.2 % and 1.4 %, respectively). In addition, in N. sp.1 nucleotide variants tended to be at intermediate frequencies, both at silent and at nonsynonymous sites (as suggested by π >θ and positive Tajima’s D values; Supplemental Table 9), whereas in N. ceranae mutations were mainly present at low frequencies —represented by π <θ estimates, and therefore by negative Tajima’s D values (Supplemental Table 9)-. Although the pooled π S for N. ceranae was 1.2 %, it should be noted that the rate of synonymous variability was two to four times lower in V. velutina than in any other Vespidae species (Supplemental Table 9), both in European and South Korean samples. Besides, when N. ceranae sequences obtained from native Hymenoptera were compared, those with an Asian origin displayed greater π S values than those from Europe, an observation that became even more evident when the sequences from all hosts were pooled across areas ( π S =1.5 % for South Korean vs. 0.5 % for European samples, respectively). In contrast, the levels of pairwise nonsynonymous diversity ( π A ) of N. ceranae did not exhibit significant differences neither across regions nor among hosts (Supplemental Table 9). The AMOVA analysis revealed that, although most of the variation of N. ceranae resided within specimens (Table 5), there was significant differentiation between European and South Korean populations – a finding also supported by the statistically significant values reached by the Snn estimator (Supplemental Table 6). 4. Discussion Four hundred and eighty-four samples collected both in the European Atlantic façade and South Korea were subject to a molecular analysis to investigate the patterns of parasite distribution in V. velutina and native Hymenoptera in these regions. First, the prevalences of three major groups of pathogens were assessed by means of PCR amplification. Despite regional variations, this analysis showed an overall pattern consistent with that described by Gabín-García et al. (2021), who reported that Trypanosomatidae were the most frequent enteropathogens in Vespidae (Vespa, Vespula, Polistes spp.) and Apidae (Bombus spp.) species in Galicia (NW-Iberian Peninsula), followed by Lipotrophidae and Nosematidae. The analysis of the distribution of these families across samples revealed that these parasite groups co-occurred significantly more often than expected by chance. This observation aligns with previous studies on pathogen populations of wild and managed pollinators (Bartolom´ e et al., 2020; Cilia et al., 2022; Evison et al., 2012; Gabín-García et al., 2021; Nanetti et al., 2021), and is likely due to shared use of habitats and resources (Graystock et al., 2015; Nanetti et al., 2021; Proesmans et al., 2021). However, despite these environmental interactions, the relative parasite prevalence varied across host species. Consistent with previous work (Gabín-García et al., 2021), V. velutina showed significantly lower parasite prevalence (measured as RPF) than native Hymenoptera, a difference that was observed across all parasite families in most regions (Fig. 1 and Table 3). The lesser occurrence of pathogens in introduced species has been attributed to the difficulty of local parasites to infect novel hosts (Torchin et al., 2003), as well as to other non-mutually exclusive factors such as: (i) Lower parasite load in the founder specimens of introduced populations (Torchin et al., 2003) —a single multimated V. velutina gyne, in the case of Europe (Arca et al., 2015). (ii) Greater immune competence against pathogens of the alien species —as reported in V. velutina gynes, but no other castes (Cappa et al., 2022). And/or (iii) larger levels of polyandria than those shown by other Hymenoptera (Arca et al., 2015), a characteristic that has been associated with greater within-colony genetic variation and lower parasite loads (Schmid-Hempel, 1999). However, this would be counter balanced by the observed reduction of genetic diversity of the invasive population as a consequence of the founder effect, at least in Europe (Arca et al., 2015; Quaresma et al., 2022), so its overall effects are difficult to ascertain. Table 6 Lipothrophidae species detected by MPS across regions and hosts. Region Host species N A. bombi Ng. sp. Portugal V. velutina 1 1 0 Galicia V. velutina 5 5 0 V. crabro 3 3 0 V. germanica 2 2 0 V. vulgaris 2 2 0 Bombus spp. 1 1 0 Basque Country V. velutina 1 1 0 V. crabro 4 4 0 V. germanica 4 4 0 V. vulgaris 1 1 0 Bombus spp. 2 2 1 France V. velutina 1 1 0 V. crabro 1 1 0 Jersey V. velutina 1 1 0 V. germanica 1 1 0 V. vulgaris 2 1 1 S. Korea V. velutina 5 5 0 V. crabro 3 3 0 V. analis 1 1 0 V. simillima 3 3 0 N, number of positive samples analysed. Ng. sp., Neogregarinorida species. Table 7 Nosematidae species detected by MPS across regions and hosts. Region Host species N N. ceranae N. bombi N. thomsoni N. sp. N. sp.1 Galicia V. velutina 8 8 0 1 1 0 V. crabro 1 1 0 0 0 0 V. germanica 1 0 0 1 1 0 Basque Country V. velutina 0 1 0 0 0 0 Bombus spp. 0 1 0 0 0 0 Jersey V. germanica 1 0 0 0 0 1 V. vulgaris 1 0 0 0 0 1 S. Korea V. velutina 3 3 0 0 0 0 V. crabro 2 2 0 0 0 0 V. analis 1 1 0 0 0 0 V. simillima 6 5 2 0 0 0 N, number of positive samples analysed; N. sp. and N. sp.1 (Nosema species and Nosema species 1, respectively). C. Bartolom´ e et al. Journal of Invertebrate Pathology 207 (2024) 108203 7 This whole scenario is consistent with the enemy release hypothesis, which suggests that reduced containment by local pathogens and predators could explain the reproductive success of some exotic invasive species, thereby posing a threat to the ecosystems (Keane and Crawley, 2002). At any rate, it should be noted that the nature of the interactions between the invasive hornets and the native pathogens remains largely unknown. The local host-parasite equilibria are likely dynamic and specific studies of host-pathogen interactions are needed to elucidate their potential role in V. velutina success in newly invaded areas. For instance, N. ceranae, currently the predominant Nosematidade species in Hymenoptera from Europe and S. Korea, shares a native range in Southeast Asia with V. velutina. It is therefore likely that they have interacted in the past and that their host-pathogen relationships have been fine-tuned over a long period of coexistence prior to their recent independent global expansions. A controlled infection assay using contaminated food demonstrated that N. ceranae spores from European honeybees could be efficiently introduced in the intestinal track of V. velutina specimens but produced no evidence of active infection (Maside et al, unpublished data). The fact that the parasite prevalences in V. velutina and native hosts become similar in summer (Fig. 1) may be caused by a combination of seasonal factors (Cilia et al., 2022; Graystock et al., 2020) as well as temporal changes in the feeding habits of V. velutina. These hornets feed on sap and nectar during their whole life cycle, choosing preferably flowers easily accessible with their short proboscis (Ueno, 2015). This morphological trait may limit their exposure to parasites (Pinilla-Gal- lego et al., 2022) during the periods of lower predatory activity, a situation that is likely reversed as their predation activity increases, promoting their contact with the pathogens carried by their prey. The maximum number of captures takes place between end of summer and early autumn, depending on the latitude (Choi, 2021; Monceau et al., 2015; Rojas-Nossa et al., 2022), which coincides with the period when the relative parasite frequency equals in both hosts. The identification of species by means of massive parallel sequencing revealed that the most abundant parasites in the dataset were those most commonly found in honeybees and bumblebees (Bartolom´ e et al., 2020; Pislak Ocepek et al., 2021), which goes in line with the findings of Gabín-García et al. (2021), although their abundance and distribution presented differences that will be discussed below. 4.1. Trypanosomatidae The most frequent Trypanosomatidae in this study were C. mellificae and L. passim, which contrasts with the lower prevalence of these species —as compared to C. bombi— detected by Gabín-García et al. (2021) in similar hosts. This observation is not particularly surprising as the distribution of pathogens in the environment is far from being homogeneous across geographical regions and seasons (Bartolom´ e et al., 2020; Graystock et al., 2020; Ivers and Jha, 2023). For instance, in this dataset, C. mellificae was more common in South Korean than in European samples, as opposed to L. passim or C. sp. 4, which displayed the reverse Fig. 3. Phylogenetic relationships of (a) SSU rRNA and (b) Actin sequences from Nosematidae (using Vittaforma corneae as an outgroup). The evolutionary history was inferred using the Neighbour-Joining method and the distances were estimated by applying the Tamura 3-parameter model assuming a gamma distribution. Bootstrap values ≥50 % are shown. Sequences obtained in this study are highlighted in bold. C. Bartolom´ e et al. Journal of Invertebrate Pathology 207 (2024) 108203 8 pattern. The use of MPS enabled the detection of a greater co-occurrence of Trypanosomatidae species compared to estimates obtained using a PCR- cloning and sequencing approach (67.1 % vs. 28.6 %, respectively) (Gabín-García et al., 2021). It also revealed higher levels of synonymous variation than previously reported (e.g. pooled π S for C. mellificae =4.4 % vs. 2.2 % or 2.4 % vs. 0.4 % for L. passim; Bartolom´ e et al., 2022). The MPS applied in this study, likely the most comprehensive and costeffective method for assessing parasite diversity (Kulski, 2016), also allowed the detection of Trypanosomatidae taxa previously described by Gabín-García et al. (2021) in V. velutina and other native hosts, which warrant further detailed investigation. The analysis of the populations of C. mellificae and L. passim revealed that both species exhibited geographic differentiation, not only in terms of variability (the pooled π S of European samples was about twice that of South Korean ones, regardless of the host) but also in what refers to their genetic structure. Although all Trypanosomatidae displayed a significant degree of geographic differentiation (Table 5), as predicted from the limited gene flow expected between European and Korean populations, this structuring was more evident in C. mellificae and L. passim. Most of the variation shown by these two species resided within European specimens, which also displayed significant haplotype subdivision among regions (Table 5). This stratification, together with the greater level of neutral diversity exhibited by European samples, seem to indicate that European populations of these species have a more ancient origin, greater population size, or both, than those from South Korea (C. mellificae likely being older than L. passim, as previously suggested by Bartolom´ e et al. (2022)). 4.2. Lipotrophidae Based on the MPS results, most Lipotrophidae-PCR positive samples corresponded to A. bombi, which was much more prevalent in European than in South Korean specimens (especially among native Hymenoptera). Despite this difference, the genetic variability of this parasite did not differ significantly between these two regions. However, it must be taken into account that the use of the SSU marker to estimate the levels of polymorphism within a species is not ideal, as the low diversity displayed by A. bombi could derive from the existence of processes of homologous recombination and/or gene conversion between paralogous ribosomal copies (both mechanisms drive to the homogenization of the sequences involved and may mask the true diversity of a species (Eickbush and Eickbush, 2007)). The best way to overcome this issue would be to investigate single copy markers, which nowadays are still unavailable. The AMOVA analysis of the haplotypes revealed differentiation between the populations from Europe and South Korea, although no significant regional subdivision was detected within each of these areas. This suggests that although the variability of A. bombi is low, at least in what concerns this genetic marker, it still enables to discriminate between populations, especially if these are distant and isolated. As in Maharramov et al. (2013), and regardless of the higher prevalence of this species in native Hymenoptera, there was no evidence of structuring by host suggesting that V. velutina specimens harbour the same haplotypes as the local entomofauna. 4.3. Nosematidae The relative prevalences of Nosematidae taxa obtained in this study varied with respect to those previously reported by Gabín-García et al. (2021) in the same region and hosts. Although the range of species detected was the same (with the exception of N. sp.1, which was newly found in this work), the most striking difference was that N. ceranae was the most frequent Nosematidae species in the present study, especially in Galicia where this species seems to be quite prevalent (Meana et al., 2017), which contrasts with the higher abundance of N. thomsoni described by Gabín-García et al. (about 10 % in V. velutina and V. crabro). This difference might reflect spatial and/or temporal variation in the prevalences of parasite species which, as mentioned above, is rather frequent. Despite the lack of diversity shown by N. bombi, N. sp. and N. thomsoni at the Actin locus, the overall level of neutral variation in N. ceranae and N. sp.1 was of the order of 1 %, as already reported for N. ceranae and N. apis at this and other loci using PCR and cloning (G´ omez-Moracho et al., 2015; G´ omez-Moracho et al., 2014; Maside et al., 2015; Roudel et al., 2013), which suggests that the filtering strategy applied to the raw MPS data was adequate and did not lead to an overestimation of the genetic variability. The reduced rates of synonymous polymorphism detected in N. ceranae isolates from V. velutina with respect to those obtained from native sympatric Hymenoptera are difficult to justify since the number of positive specimens sequenced was similar in the two groups (Table 7) and both shared the same local environment and were theoretically exposed to a similar array of parasite variants (Ke et al., 2022). This discordance could be explained by the differences observed in the gut microbial composition of these invasive hornets, which might play a role in the defence against some pathogens, and thus contribute to reduce their diversity (Cini et al., 2020; Do et al., 2023; Zhang et al., 2022). This would be in line with the significantly lower RPF found in V. velutina as compared with native Hymenoptera (see above). Consistent with the assumptions that N. ceranae comes from East Asia (Botías et al., 2012), and that the source populations should be more diverse than those derived from them, South Korean samples displayed higher levels of synonymous polymorphism ( π S ) than European ones, which lines up with previous reports describing more variability in Asian populations of N. ceranae than in those from areas recently colonised by this pathogen (G´ omez-Moracho et al., 2015; Ke et al., 2022; Maside et al., 2015; Valizadeh et al., 2022). In agreement with this scenario, the pooled π A was 2.4 times higher than π S in European isolates (1.2 % vs. 0.5 % for π A and π S , respectively), as a consequence of selection being less effective at eliminating weakly deleterious mutations in small (=derived) populations (Charlesworth et al., 1993), while South Korean ones exhibited similar values at both types of sites (1.3 % vs. 1.5 % for π A and ππ S , respectively). Another observation arising from the analysis of the patterns of mutation frequency of these samples was that these N. ceranae populations are experiencing a demographic expansion (G´ omez-Moracho et al., 2015), as suggested by the excess of rare variants detected at synonymous and nonsynonymous positions in both groups of isolates (Tajima, 1989a). In contrast, the frequency spectrum of the polymorphisms shown by N. sp.1 seems to indicate that, at least in Europe, this taxon underwent a population contraction that would also be reflected by the positive Tajima’s D values obtained both at silent and replacement sites (Tajima, 1989a). The population differentiation found between European and South Korean N. ceranae isolates may arise from different factors, in addition to the absence of gene flow between these isolated populations. One possibility is that the population/s that colonised Europe and South Korea had a different origin, which is rather likely if Asian populations are ancient and genetically structured. Another possibility, not mutually exclusive, could be that haplotypes that were at low frequencies in the ancestral population had drifted to higher frequencies after the colonisation bottleneck. In any case, to distinguish between these options, and clarify the demographic history of this pathogen in Asia, would require further population data from its native range. 5. Conclusions Overall, the main drivers of variation in parasite prevalences are the geographic origin of the samples and the nature of the hosts, depending on whether these are native or invasive. V. velutina exhibits lower parasite prevalence than the local entomofauna, potentially contributing to its successful spread into new territories. Another influential C. Bartolom´ e et al. Journal of Invertebrate Pathology 207 (2024) 108203 9