The honey bees Apis mellifera and Apis cerana are neither protected nor endangered species. Ethical approval is not required for this study.

We designed a two-by-two factorial experiment to study parasite virulence and proliferation in two honey bee species. We used two parasite sources P_Acer (Parasite spores purified from the honey bee Apis cerana) and P_Amel (Parasite spores purified from the honey bee Apis mellifera), as well as two host species H_Acer (Host honey bee Apis cerana) and H_Amel (Host honey bees Apis mellifera) ( Supplementary Table S1 ). We combined bees from different hives for the inoculation, and uninfected bees (H_Acer and H_Amel) as controls for this multi-host-parasite experiment. The honey bee colonies are maintained in the experimental apiary at Jiangxi Agricultural University.

Three hundred honey bees of each A. mellifera and A. cerana were captured near the hive entrance using an insect net. The midgut was dissected and homogenized to isolate N. ceranae spores that were further purified using the Percoll gradient (Chen et al., 2013). The spores were counted under the light microscope using a hemacytometer. The sealed brood frames from three hives of H_Amel and H_Acer were kept in an incubator to collect newly emerged bees (35°C, 75% humidity).

The newly emerged honey bee workers (H_Acer and H_Amel) were individually inoculated with 2 µL of sugar solution with 10⁵ spores. Additional freshly emerged honey bees (< 24 h after emerging) were each fed 2 µL of sugar solution without spores as the uninfected control group. One hundred fifty honey bees were inoculated in each treatment group, and the cohorts were divided into 3 rearing cups (50 bees per cup) in an incubator (35°C, 75% humidity) (3 replicates, Supplementary Table S1 ) to investigate the general bee response and parasite proliferation. During the experiment, sucrose (50% w/w) was provided ad libitum as the only food. In each rearing cup, three bees were collected at 24 h intervals from 1 to 5 dpi (day post-inoculation) for RNA-seq. The remaining bees were dissected to count spores at 14 dpi.

As the parasite infects the epithelial cell in the midgut, we dissected midgut tissue for RNA-seq. Three bees per cup per day were dissected and pooled for RNA extraction using Trizol. The library was prepared using the NEBNext Ultra RNA Kit. In total, 90 RNA libraries (5 days * 6 treatments * 3 replicates) were sequenced on the Illumina NovaSeq 6000 platform.

The quality of RNA reads (150bp, paired-end) was first viewed using Fastqc (Andrews, 2010) and trimmed using the Seqtk package with default parameters (Li, 2022). The processed reads were aligned to N. ceranae (Version Ncer 3.0), A. cerana (version CC1.0), and A. mellifera (version Hav3.1) genome, respectively, with the Hisat2 package with default parameters (Kim et al., 2015; Diao et al., 2018; Wallberg et al., 2019; Huang et al., 2021). The read count per gene was retrieved using bedtools (Quinlan and Hall, 2010). The within-group dispersion was calculated from the three replicates to determine the significantly regulated genes with edgeR package and adjusted for FDR (false discovery rate) (Robinson et al., 2010). The genes with FDR < 0.05 were defined as significantly regulated ones. Gene Ontology (GO) terms were retrieved using EggNOG-mapper, and the enrichment analysis was performed using the TopGo package with an adjusted weighted ks test (Robinson et al., 2010; Alexa and Rahnenfuhrer, 2021). Bee survival was analyzed using the Kaplan-Meier estimate in the survival package, adjusted for FDR (R Core Team, 2013; Therneau, 2022). The variance of the spore load among the treatment groups was analyzed with the Wilcoxon rank test, and the p values were adjusted with FDR to reduce false positives. The impact of the spore source and day on the gene expression was analyzed using ANONA, where the day and parasite source were fixed factors, and the replicates were random factors. The number of up and down-regulated genes was analyzed using Pearson’s Chi-squared test.

The uninfected H_Amel shows the highest survival (99.4% survival), followed by H_Amel infected with P_Acer (H_Amel_P_Acer, 98.5% survival) and the H_Amel_P_Amel group (96.9% survival). The bees in the H_Acer_P_Acer group show the lowest survival (86.6% survival). The parasite source (P_Amel and P_Acer) shows a minor impact on the survival of H_Amel (H_Amel_P_Acer vs H_Amel_P_Amel, Kaplan-Meier test, P > 0.05) ( Figure 1A , Supplementary Table S2 ). Comparatively, the parasite P_Acer causes significantly higher mortality in H_Acer than P_Amel (Coxph test, P < 0.001). The uninfected H_Amel survives substantially better than the uninfected H_Acer (Coxph test, P < 0.01). Overall, the parasites cause a higher mortality in H_Acer than in H_Amel. Additionally, we normalize the mortality variance between the two honey bee species using uninfected_H_Amel and uninfected_H_Acer. Again, the parasites cause higher mortality in H_Acer than H_Amel (P < 0.001, Supplementary Figure S1 ).

The spores are not found in the uninfected bees. The spore load is not evenly distributed among the four infected honey bee groups (Kruskal-Wallis test, χ² = 43.1, df=3, P < 0.0001, Figure 1B ). The parasite produces more spores in H_Amel than H_Acer, when infected by either P_Amel or P_Acer (Wilcoxon rank sum test, df=1, P < 0.001) ( Supplementary Table S3 ). Thus, host species substantially impact the parasite proliferation (F=32.3, df=2, ANOVA, P < 0.0001).

To investigate the parasite gene expression profile, we quantify the parasite transcriptome in the two host species ( Figure 1C ). H_Acer shows a stronger tendency to suppress the parasite gene expression than H_Amel (F=77.1, df=1, ANOVA, P < 0.0001). A significantly lower number of up-regulated genes are observed in P_Amel than P_Acer when infecting the H_Amel (Pearson’s Chi-squared test, χ² = 20.8, df=2, P < 0.0001). The highest variance was at four dpi, and the number of up-regulated genes is threefold higher in P_Acer (36 genes) than in P_Amel (12 genes) in H_Amel ( Supplementary Table S4 ). Thus, the parasites express a higher number of genes and transcript levels in H_Amel than in H_Acer.

To infer how strong is the impact of the infection on the host, the honey bee transcriptomes responding to each of the two parasite sources are quantified as well. We find the expression of bee (H_Amel and H_Acer) genes is suppressed, where more down-regulated genes than up-regulated ones when infected by either P_Amel or P_Acer ( Table 1 , Pearson’s Chi-squared test, χ² = 51.0, df=4, P < 0.001). Comparatively, more genes were regulated by the infection in H_Acer than H_Amel. The bee hosts share a common subset of genes that respond to infection ( Figure 1D ). The immune response toward the infection is small in H_Amel, reflected by only two down-regulated immune genes in Toll pathway involved in pathogen recognition (LOC412536) and melanization (LOC406115). We found both up and down regulated immune genes in H_Acer. Particularly, the Toll pathway is generally suppressed, including the genes involved in parasite recognition PGRP (APICC_00292), signal transport SPZ (APICC_01852), and antimicrobial peptides defensin and lysozyme (APICC_08301, APICC_03572, APICC_08272). These suppressed genes have previously been confirmed using qPCR (Ke et al., 2022). A few genes are continuously down-regulated in H_Amel, reflecting genes inhibited by the infection ( Supplementary Figure S2 ). GO enrichment analysis indicates that apoptosis regulation (GO:0043066, FDR< 0.0001), cell cycle phase transition (GO:0044843, FDR < 0.0001), and a few other biological functions are altered in the infected honey bee genes ( Figure 2 ).