Skip to content
Frontiers
Search

Original Research

Frontiers in Insect Science ·

Volume 6 - 2026 | https://doi.org/10.3389/finsc.2026.1828903

Regulation of Nosema ceranae gene expression by Bidens pilosa phytogenic treatment in infected honey bees

  • Yi-Hsuan Li1
  • Kai-Di Chen2
  • Ju-Chun Chang2
  • Zih-Ting Chang3,4
  • Chong-Yu Ko4
  • Shang-Tse Ho5
  • Li-Hung Chen6,7
  • Ming-Cheng Wu1,2
  • Wen-Chin Yang8
  • Yu-Liang Yang8
  • Yue-Wen Chen4
  • Yu-Shin Nai1,2*corresponding author
  1. 1Doctoral Program in Microbial Genomics, National Chung Hsing University and Academia Sinica, Taichung, Taiwan
  2. 2Department of Entomology, National Chung Hsing University, Taichung, Taiwan
  3. 3Department of Entomology, National Taiwan University, Taipei, Taiwan
  4. 4Department of Biotechnology and Animal Science, National Ilan University, Yilan, Taiwan
  5. 5Department of Wood-Based Materials and Design, National Chiayi University, Chiayi, Taiwan
  6. 6Department of Plant Pathology, National Chung Hsing University, Taichung, Taiwan
  7. 7Advanced Plant Biotechnology Center, National Chung Hsing University, Taichung, Taiwan
  8. 8Agricultural Biotechnology Research Center, Academia Sinica, Taipei, Taiwan
Article metrics0 viewsView details
PDFCiteMetrics

Abstract

Nosema ceranae is an obligate microsporidian parasite that causes nosemosis in honey bees (Apis mellifera). In controlled laboratory cage experiments, newly emerged worker bees were experimentally infected with N. ceranae spores and treated with Bidens pilosa phytogenic extract (BP). Each treatment group consisted of 50 bees per replicate, and experiments were repeated independently. To explore the molecular basis of BP-associated effects, we conducted transcriptome profiling of N. ceranae from infected honey bee midguts at multiple time point. RNA samples from infected bee midguts with BP treatment were collected at 5, 10, and 20 days post-infection (dpi) for transcriptomic analysis. BP treatment significantly improved survival probability and reduced pathogen load compared to infected controls. During the early infection stage (5 dpi), BP treatment was associated with extensive downregulation of parasite genes, including components of the V-type ATPase pathway. Gene Ontology and KEGG analysis suggested suppression of metabolic and ion transport processes. To further evaluate the potential role of V-type ATPase, RNAi-mediated knockdown resulted in reduced gene expression and showed a trend toward decreased pathogen load and modest improvement in host survival. Although the RNAi results do not provide definitive evidence of causality, they support a potential involvement of the V-type ATPase in parasite proliferation. Overall, BP altered the transcriptomic profile of N. ceranae in a stage-dependent manner and may influence parasite development by affecting key metabolic pathways.

Keywords

1 Introduction

The health of honey bees has been a significant concern for both beekeepers and scientists. Multiple factors negatively impact honey bee health, including environmental pollution, food scarcity, poor nutrition, pesticide exposure, and various pathogens (1–4). Among these threats, diseases caused by bacteria, viruses, and fungi are major contributors to colony decline. Notably, nosemosis is highly prevalent in honey bee populations (5–7). This disease is caused by two species of microsporidia, Nosema apis and Nosema ceranae, which primarily infect adult honey bees (8, 9). Infection with N. ceranae weakens honey bees and, in severe cases, can lead to substantial declines in adult bee populations, reduced honey yields, and ultimately, colony collapse. The infection impairs hormone production, lipid synthesis, and midgut tissue integrity (10). It also induces nutritional and energetic stress, suppresses the immune system, and interferes with host cell apoptosis (11). Infected bees exhibit impaired olfactory learning and memory, premature foraging behavior, reduced homing ability, and decreased flight capacity (10, 12). N. ceranae infection is associated with increased energetic demand and may contribute to malnutrition and elevated mortality (12, 13).

Given the widespread occurrence of this disease and the limited efficacy of existing treatments, there is a continued need to identify safe and effective control strategies. Researchers have devoted substantial efforts to developing strategies to control N. ceranae and mitigate its impact on honey bee populations. While the antibiotic fumagillin was previously considered an effective treatment (14, 15), recent studies have demonstrated its limited efficacy against N. ceranae infections (15). Regardless of dosage or method of administration, fumagillin treatment in heavily infected colonies has shown no significant improvement in winter colony size or survival rates (14). Consequently, ongoing research is focused on identifying new, safe, and effective alternatives for the prevention and treatment of N. ceranae infections in honey bees. Based on the above-mentioned issue, an alternative approach has been investigated for the management of N. ceranae infections. In addition, numerous studies have evaluated the effectiveness of organic extracts and natural supplements, such as commercial products and plant-based supplements (16–18). These include gastric pentadecapeptide BPC 157, BEEWELL AminoPlus, and other commercial products (18, 19). Our previous study also showed that dimethyl sulfoxide (DMSO) has a positive therapeutic effect on honey bees infected with N. ceranae (20).

Among these, extracts from a variety of plants have been extensively studied and show potential in controlling N. ceranae. Extracts from Annona squamosa, Ocimum basilicum, Psidium guajava, and Syzygium jambos have exhibited strong anti-microsporidian activity, inhibiting the development of N. ceranae spores, with effect comparable to fumagillin (21). Furthermore, a herbal extract mixture comprising Rumex acetosella, Achillea millefolium, Plantago lanceolata, Salvia officinalis, Thymus vulgaris, Rosmarinus officinalis, and Laurus nobilis has demonstrated efficacy in comparison to fumagillin (22).

One such plant, B. pilosa, is of particular interest due to its known antiparasitic properties and its role as a food source for bees. B. pilosa is an erect, perennial herb belonging to the Asteraceae family. Native to South America, it has spread widely across tropical regions worldwide. In many areas, BP serves as a major pollen source for A. mellifera, contributing significantly to bee bread production (23, 24). In Taiwan, B. pilosa is commonly used in the production of capsules, decoctions, and tinctures derived from its dried powder, which are marketed as dietary supplements or health food products (25–27). Crude root extracts of B. pilosa have demonstrated in vitro antiparasitic activity against the protozoan parasite Plasmodium falciparum (28, 29). Moreover, phytogenic extracts from B. pilosa have shown effectiveness against Eimeria parasites, the causative agents of coccidiosis in poultry, while also helping to reduce the risk of drug resistance (27, 30).

As mentioned above, BP has shown inhibitory activity against protistan pathogens; however, limited data are available regarding its effects on nosemosis. Therefore, this study aims to evaluate the effects of BP on N. ceranae infection in A. mellifera. Based on the transcriptomic analysis, particular attention was given to genes involved in the V-type ATPase pathway. Furthermore, RNA interference (RNAi) was employed to investigate the potential role of this pathway during infection and to further assess BP-associated effects.

2 Materials and methods

2.1 Honey bees

Honey bees (A. mellifera) used in all experiments were collected from colonies maintained at the apiaries of National Ilan University (NIU, Taiwan) and National Chung Hsing University (NCHU, Taiwan). Each selected colony consisted of six to eight frames and was confirmed to have a healthy, laying queen to ensure normal brood development. Colonies were routinely monitored for general health status, and those with low or negligible N. ceranae infection levels were selected. Standard beekeeping management practices, including regular syrup feeding and routine Varroa mite control, were maintained throughout the study. Brood frames with more than 70-80% of the cells capped were collected and incubated at 34 °C until adult emergence. Newly emerged worker bees were collected and transferred into rearing cages (90 mm in diameter × 55 mm in height). Bees were maintained under laboratory conditions and provided ad libitum access to 50% (w/v) syrup as their food resource.

2.2 Extraction of Bidens pilosa

BP was extracted with methanol at a 10:1 (v/v, methanol:BP) ratio. Whole parts of B. pilosa were collected, washed, air-dried, and ground into powder prior to extraction. The extraction was performed at room temperature, followed by filtration and solvent evaporation under reduced pressure to obtain the crude extract. The extract was stored at −20 °C until use. After extraction, the BP extract was freeze-dried as described previously (30). The dried BP extract was weighed and then dissolved in DMSO to a final concentration of 100 ppm for further use (Sigma Chemical Co., St. Louis, MO).

2.3 N. ceranae infection and BP treatment

Spores of N. ceranae were obtained and purified from the midgut tissues of forager honey bees collected at the apiaries, following protocols described in previous studies (31, 32). The purified spores were quantified using a hemocytometer (Paul Marienfeld GmbH & Co. KG, Baden-Württemberg, Germany). To assess the effects of BP on N. ceranae infection, newly emerged adult honey bees were randomly divided into four groups: Blank (honey bees fed with 50% syrup), Nc (honey bees infected with N. ceranae and fed with 50% syrup), DMSO+Nc (honey bees infected with N. ceranae and fed with 100 ppm DMSO in 50% syrup), and BP+Nc (honey bees infected with N. ceranae and fed with 100 ppm BP (dissolved in DMSO) in 50% syrup). For infection, each honey bee was individually fed with 2 μL of N. ceranae spore suspension containing 1×105 spores/μL. A total of 50 honey bees were included per group. After inoculation, honey bees were maintained in rearing cages at 34 °C and monitored daily. This experiment was repeated four times for genomic copy number quantification, three times for RNA sequencing (RNA-seq), and four times for quantitative real-time PCR (qPCR) validation. For the RNA-seq, midgut tissues from three honey bees per replicate in each group were collected at 0 (3 hours post-infection, hpi), 5, 10, and 20 days post-infection (dpi) for RNA extraction and library construction.

2.4 Quantification of N. ceranae genome copy number

Confirmation of N. ceranae infection followed the protocol described in our previous study (32). At 0, 1, 3, 5, 10, 15, and 20 dpi, three midguts from each group were collected into 1×TE buffer, homogenized with a pestle, and examined under a light microscope (WHITED, Taiwan). Total DNA was extracted using the EDNA HISPEX tissue kit (Fisher Biotec, Stirling, Australia), then diluted 100-fold. Absolute quantitation of N. ceranae genome copy number was conducted via qPCR using primer set: SSU-F/SSU-R (Supplementary Table 1) on a CFX Connect Real-Time PCR Detection System (Bio-Rad, California, USA). Each 20 μL reaction contained: 10 μL 2× iQ™ SYBR® Green Supermix, 1 μL of each forward (10 mM) and reverse (10 mM) primer, and 1 μL of 100-fold diluted DNA template. The qPCR reactions were performed using a two-step amplification protocol by the following procedure: initial denaturation at 95 °C for 3 min, 40 cycles of denaturation at 95 °C for 10 s and annealing at 59 °C for 30 s, a final extension at 72 °C for 10 s, and a melting curve at 65 °C to 95 °C with an increase of 0.5 °C in 5 s. All samples were performed in quadruplicate. A standard curve was generated from DNA extracted from a known spore suspension (1×107 spores/100μL), serially diluted 10-fold down to 1×100 spores/mL. Threshold cycle (Ct) values were plotted against the log10 of the initial spore concentration, yielding a linear regression equation:: y = 3.0383x+42.084 (R2 = 0.997; Supplementary Figure 1). This approach is consistent with established absolute qPCR quantification techniques using standard curves to derive copy number from Ct values. To correct for the sample preparation, calculated copy numbers were multiplied by the factors: 2 (DNA extraction dilution correction) ×100 (template dilution) ×100 (volume factor, 1 μL of diluted DNA). This experiment was repeated four times.

2.5 RNA sequencing and raw data processing

Midgut tissues from honey bees were homogenized using a Tissue Lyser II (QIAGEN, Solingen, Germany). Total RNA was extracted following the manufacturer’s instructions for TRIzol reagent (Life Technologies, USA) in combination with Locking Gel (BIOTOOLS, Taiwan). RNA purity and concentration were assessed using a Nanodrop 1000 Spectrophotometer V3.5 (Thermo Scientific) and a Qubit 2.0 Fluorometer (Life Technologies, USA), respectively. RNA from three biological replicates was pooled for library construction to obtain sufficient RNA quantity and represent an averaged transcriptomic profile (33, 34). While this approach limits statistical inference, therefore, the data was interpreted as exploratory. The mRNA enrichment was performed using poly-A mRNA with oligo-d (T) magnetic beads included in an Ultra RNA Library Prep Kit (NEB, Essex, MA, USA). The mRNA was fragmented, and first- and second-strand cDNA synthesis was conducted using NEBNext random primers followed by end-repair, A-tailing, and adaptor ligation (NEBNext Adaptor). Libraries were enriched via PCR, and quality and quantity were confirmed using an Agilent Fragment Analyzer and qPCR. Sequencing was performed on an Illumina MiSeq platform (301 bp paired-end reads) by Genomics BioSci. & Tech. Co., Ltd. (Taipei, Taiwan). Raw reads were processed by adapter removal using Trimmomatic (35), and quality trimming and filtering were performed with PRINseq (36).

2.6 Differentially expressed gene (DEG), GO, and KEGG analysis

Because of N. ceranae is an obligate intracellular parasite, it is unculturable, therefore, high-quality reads were aligned to the N. ceranae reference genome (ASM98816v1) using HISAT2 with default parameters. Only reads mapped to the N. ceranae reference genome were retained for downstream analyses, allowing selective analysis of parasite-associated transcripts while excluding the majority of host-derived sequences. Gene annotation information (GCF_000988165.1_ASM98816v1_genomic.gff) was obtained from the National Biotechnology Information Center (NCBI) (https://www.ncbi.nlm.nih.gov/ accessed on 19 January 2021) for identifying N. ceranae-related PE reads. Differential gene expression analysis across experimental groups was performed using Cuffdiff. Genes with zero Fragments Per Kilobase per Million (FPKM) values across all samples were excluded; for the remaining genes, a pseudocount of 1 was added for log transformation (log2(FPKM + 1)). Genes exhibiting a log2 ratio ≥ 1 or <-1 (fold change ≥ 2 or <- 2) were considered putative differentially expressed between compared groups (DMSO+Nc and BP+Nc) (D/B) for 5, 10, and 20 dpi. Significantly altered genes (DEGs) were subjected to Gene Ontology (GO) enrichment analysis using the hypergeometric test, with a significance defined as P < 0.05. For gene orthologue assignment and pathway mapping, they were submitted to the KEGG Automatic Annotation Server (KAAS, http://www.genome.jp/tools/kaas/ accessed on 23 April 2021) (37).

2.7 Quantitative PCR (qPCR) validation

To validate the expression patterns of selected DEGs, midgut tissues from three honey bees per group were collected at time points 0, 0.5, 1, 3, 5, 10, 15, and 20 dpi, with four biological replicates per time point. Each sample (~300 μL) was homogenized in Trizol reagent, and total RNA was extracted using the GENEzolTM TriRNA Pure Kit (Geneaid, New Taipei, Taiwan) following the manufacturer’s protocol. RNA purity and concentration were assessed using a Nanodrop 1000 Spectrophotometer V3.5 (EZdrop 1000, New Taipei, Taiwan). cDNA was synthesized from 1 ng/µL of RNA using the BioDiamond RT-PCR mix (Origin Pure BioSci & Tech. Co., Yilan, Taiwan). Quantitative RT-PCR (qRT-PCR) was then performed using gene-specific primer sets (Supplementary Table 1). N. ceranae β-tubulin was used as an internal control to validate the expression levels of selected DEGs associated with V-type ATPase, including ATPeV1A (AAJ76_4700027584), ATPeV1B (AAJ76_700039738), ATPeV1D (AAJ76_1600056373), ATPeV1E (AAJ76_1700024327), ATPeV0B (AAJ76_1900038297), ATPeV0C (AAJ76_300040219), and ATPeV0D (AAJ76_630003653) (Supplementary Table 1). The qRT-PCR reactions and thermal cycling were conducted as described above; the annealing temperature of each primer set is summarized in Supplementary Table 1. The relative gene expression level was determined by the 2-ΔCt method (log2-transformed).

2.8 Gene knockdown test

To evaluate the role of ATPeV0B in the N. ceranae infection, a 212 bp dsRNA (dsVATP), covering nucleotides 253–464 of the gene (Supplementary Figure 2), was synthesized. This was achieved via Golden Gate cloning into the pSIG1 vector with the primer sets listed in Supplementary Table 1 following the protocol (38), and expressed in E. coli HT115. The dsVATP was purified with 8 M RNase‐free LiCl (38), air‐dried, and resuspended in 750 µL RNase‐free water for experimental use (39). Four treatment groups (n= 50 honey bees each) were established, including Blank (honey bees fed with 50% syrup), Nc (honey bees infected with N. ceranae and fed with 50% syrup), dsGFP+Nc (honey bees infected with N. ceranae and fed with GFP dsRNA in 50% syrup), and dsVATP+Nc (honey bees infected with N. ceranae and fed with dsVATP in 50% syrup). Infections were performed by feeding newly emerged honey bees 2 μL of a 50% syrup containing 5×103 spores/μL (N. ceranae). A low spore concentration was used in the RNAi experiment to avoid excessive infection intensity and allow clearer observation of gene knockdown effects (40, 41). Infections were performed by feeding newly emerged honey bees 2 μL of a 50% syrup containing 5×103 spores/μL (N. ceranae). Honey bees were kept at 34 °C with 95% RH in incubators. Starting at 5 dpi, bees in the dsRNA groups received 2 μg dsRNA in syrup every two days. Daily survival was monitored for Kaplan–Meier analysis up to 15 dpi (42). Midgut tissues (three honey bees per group) were collected at 5, 10 (5 days post-dsRNA treatment), and 15 dpi (10 days post-dsRNA treatment) for DNA and RNA extraction as described above. The spore counts, N. ceranae genome copies, and expression of ATPeV1A, ATPeV1B, ATPeV1D, ATPeV1E, ATPeV0B, ATPeV0C, and ATPeV0D were quantified by qPCR, using the method described above. The relative gene expression level of the dsGFP and dsVATP was determined by the 2-ΔCt method (43, 44). For statistical analysis, values were log2-transformed prior to performing Student’s t-test.

2.9 Statistical analysis

All statistical analyses and visualization were performed using R version 4.0.3 (45) using various packages, including heatmap, ggplot2, and ggpubr (46–48). Kaplan–Meier survival curves were analyzed using the survminer and survival packages in R, with statistical significance assessed by log-rank (Mantel–Cox) tests and Bonferroni correction for multiple comparisons (42, 49). Spore count and qPCR data, including genome copy and gene expression, were expressed as mean ± standard deviation (SD) from four separate experiments. The differences between groups were analyzed using one-way ANOVA followed by Tukey’s multiple comparisons test. Differences were deemed statistically significant at p<0.05.

3 Results

3.1 N. ceranae infection

Kaplan–Meier survival analysis revealed that the Nc, DMSO+Nc, and BP+Nc exhibited significantly reduced survival compared to the Blank. Among the infected groups, the BP+Nc showed higher survival probability compared to other infected groups, while no significant difference was observed between the BP+Nc and DMSO+Nc (Figure 1A). Spores appeared earlier in the Nc (3 dpi) than in the DMSO+Nc and BP+Nc (5 dpi). By 10 dpi, spore counts were significantly highest in the Nc (5.4×106 spores per bee), while the BP+Nc had the lowest count (1.9×106 spores per bee), which was significantly lower than that of the Nc (Figure 1B). At later time points, spore counts remained higher in the Nc than in the DMSO+Nc and BP+Nc, although these differences were not statistically significant. Genome copy numbers were detected as early as 0 dpi in all infection groups. Genome copy numbers increased markedly after 10 dpi in all infection groups and then declined thereafter, with relatively lower levels observed in the BP+Nc compared with the Nc (Figure 1C).

Observation and molecular identification of *N. ceranae* infection. **(A)** Kaplan-Meier survival curves for Blank, Nc, DMSO+Nc and BP+Nc. Pairwise log-rank (Mantel-Cox) tests with Bonferroni correction for multiple comparisons were used to compare each treatment to Nc (*P < 0.05; ***P < 0.001). **(B)** spore counts and **(C)** genome copies. dpi = days post-infection. For **(B, C)**, statistical significance among treatment groups was assessed separately for each time point using one-way ANOVA followed by Tukey’s multiple comparisons test. Data are shown as mean ± SD (n = 4).
Figure 1. Observation and molecular identification of N. ceranae infection. (A) Kaplan-Meier survival curves for Blank, Nc, DMSO+Nc and BP+Nc. Pairwise log-rank (Mantel-Cox) tests with Bonferroni correction for multiple comparisons were used to compare each treatment to Nc (*P < 0.05; *P < 0.001). (B) spore counts and (C) genome copies. dpi = days post-infection. For (B, C)**, statistical significance among treatment groups was assessed separately for each time point using one-way ANOVA followed by Tukey’s multiple comparisons test. Data are shown as mean ± SD (n = 4).

3.2 Sequencing data summary

The relatively low mapping rate reflects the dominance of host-derived RNA in midgut samples. RNA sequencing generated between 1.47 million (M) to 6.85 M across the time points. After quality filtering, the number of clean reads ranged from 1.35 M to 2.09 M. The proportion of reads mapped to the N. ceranae genome ranged from 0 to 9.6% in the DMSO+Nc and BP+Nc, while the spore sample showed a higher mapping rate of up to 12% (Supplementary Table 2). This low mapping rate likely reflects the predominance of host-derived RNA in whole midgut samples.

3.3 DEG analysis

In the D/B comparison, the majority of DEGs were downregulated at 5 dpi, whereas more genes were upregulated at 20 dpi, indicating a time-dependent transcriptional response to BP treatment (Figure 2A; Table 1). A total of 1006, 431, and 1013 DEGs were identified at 5, 10, and 20 dpi, respectively (Figure 2B; Supplementary Table 3 and S4).

Differential expression analysis between BP+Nc and DMSO+Nc across infection time points. **(A)** Scatter plots of log<sub>2</sub> FPKM values for BP+Nc versus DMSO+Nc at 5, 10, and 20 days post-infection (dpi). Each dot represents a gene. Genes with |fold change| ≥ 2 are colored: red = significantly upregulated in B/D, blue = significantly downregulated in B/D; non-significant genes are shown in black. **(B)** Histogram showing the number of DEGs in defined fold-change bins at each time point. Fold change of -16 to -8, -8 to -4, -4 to -2 for downregulated; 2 to 4, 4 to 8, 8 to 16, 16 to 32 for upregulated. Numbers above bars indicate DEG counts per bin.
Figure 2. Differential expression analysis between BP+Nc and DMSO+Nc across infection time points. (A) Scatter plots of log2 FPKM values for BP+Nc versus DMSO+Nc at 5, 10, and 20 days post-infection (dpi). Each dot represents a gene. Genes with |fold change| ≥ 2 are colored: red = significantly upregulated in B/D, blue = significantly downregulated in B/D; non-significant genes are shown in black. (B) Histogram showing the number of DEGs in defined fold-change bins at each time point. Fold change of -16 to -8, -8 to -4, -4 to -2 for downregulated; 2 to 4, 4 to 8, 8 to 16, 16 to 32 for upregulated. Numbers above bars indicate DEG counts per bin.

Table 1. List of the top 10 DEGs by fold change for up- and down-regulated genes at each dpi.

RegulationdpiGene IDDMSO+Nc FPKMBP+Nc FPKMFold changeNCBI protein (gene) name
Upregulated5AAJ76_20001515760.2937120.7964.543928079dna-directed rna subunit e
AAJ76_600092331232.924467.9763.922844096polycomb enhancer epc
AAJ76_1000131874144.58292.4133.903315287trna pseudouridine synthase a
AAJ76_200015201123.726250.5593.787484541wd40 domain-containing poly + rna export protein
AAJ76_400003711950.4602102.7573.378423061histone acetyltransferase
AAJ76_439000349.2658102.8953.270066557endonuclease-reverse transcriptase, partial
AAJ76_450001529989.9253192.8343.03162225hypothetical protein
AAJ76_1400060351100.089216.4953.008928865hypothetical protein
AAJ76_3500010795552.6371200.292.98541191hypothetical protein
AAJ76_800026843198.724443.0992.977600082dna replication complex gins protein sld5
10AAJ76_200014945836.0558432.63511.99906428dna repair protein rad4
AAJ76_400015050031.7441356.56611.23252214iron hydrogenase
AAJ76_17000612054515.31360767.98969261616S ribosomal RNA
AAJ76_90008231144.0685351.9157.985650874mrna capping enzyme
AAJ76_100094159901.737005.927.769436678aaa atpase domain-containing protein
AAJ76_200004276118.2023131.2457.210353541hypothetical protein
AAJ76_178000505128.8713206.97.16630775hypothetical protein
AAJ76_390003313013.209693.77727.09916822blood-stage membrane protein ag-1
AAJ76_100010033011.154678.52057.039290378endoribonuclease dicer
AAJ76_130005962627.746172.0866.20219979hypothetical protein
20AAJ76_1000013104203.3383317.1416.31344985hypothetical protein
AAJ76_150004364941.1108586.58514.2683835720s proteasome subunit
AAJ76_80000127281264.652849.9714.00276669ring zinc finger domain-containing protein
AAJ76_730007221153.752397.08213.30097409pyruvate dehydrogenase e1 component beta subunit
AAJ76_60004577417.737511.39312.60068274cullin like protein
AAJ76_200012567977.9985969.30612.42720542skt5-like protein
AAJ76_70008988689.3968182.2849.790550625rrp4-like rna-binding protein
AAJ76_2300055285271.1812644.719.752553449dna-directed rna polymerase i
AAJ76_3300035292219.4381977.179.01014298ras-related cell division protein 42-like protein
AAJ76_28000418922.8019191.078.379591993subtilisin-like endopeptidase
Downregulated5AAJ76_160006850216700.71134.28-14.7236147ccaat binding transcription factor subunit a
AAJ76_50001167882187.76160.32-13.646225hypothetical protein
AAJ76_10001667271975.56159.947-12.3513771u3 small nucleolar ribonucleoprotein
AAJ76_1500062863548.68745.0673-12.1748323hypothetical protein
AAJ76_20001045325865.4486.924-12.0458141d-tyrosyl-trna deacylase
AAJ76_22000193231212.12102.504-11.8250961hypothetical protein
AAJ76_2000840471003.1397.6616-10.2714646chromatin assembly factor 1 p60 subunit
AAJ76_18000161402000.5198.067-10.1001163hypothetical protein
AAJ76_38000259332213.85219.989-10.0634285g2 mitotic specific cyclin 1
AAJ76_100058174465.53546.6012-9.98975987hypothetical protein
10AAJ76_1600024122493.18749.1088-10.0427328pp-loop atpase
AAJ76_830001211364.22947.3456-7.69296636hypothetical protein
AAJ76_14000637303071.49406.382-7.55813349ubiquitin
AAJ76_3000392391875.57265.234-7.07137123hypothetical protein
AAJ76_4200011647701.31999.3674-7.05785602hypothetical protein
AAJ76_1000062515372.98654.0805-6.89686152hypothetical protein
AAJ76_90000126431454.27212.903-6.83068375hypothetical protein
AAJ76_14000566762517.17456.658-5.51213033transmembrane adaptor erv26-like protein
AAJ76_56000183348856.261753.21-5.05145852protein transport protein yos1
AAJ76_1300025816116.53223.184-5.02638089hypothetical protein
20AAJ76_1260005835296.07650.0031-5.92116066chitooligosaccharide deacetylase-like protein
AAJ76_260003138143.09724.7229-5.78806049hypothetical protein
AAJ76_300094826280.73248.8358-5.74847788tyrosyl-trna synthetase
AAJ76_860009304381.2876.3399-4.99449911pol
AAJ76_93000341187.460818.0246-4.8523042hypothetical protein
AAJ76_90002423223.96715.13648-4.66604429dynein heavy chain
AAJ76_1400045858761.568169.459-4.49412421transcription initiation factor tfiid subunit
AAJ76_1600062669345.30686.2473-4.00366149exosome rnase ph-like protein
AAJ76_100016514323.28280.9325-3.99445866hypothetical protein
AAJ76_1900051261199.20650.4844-3.94588798hypothetical protein

3.4 GO annotation and KEGG pathway analysis

GO annotation analysis revealed time-dependent functional changes in DEGs, indicating dynamic effects of BP treatment on N. ceranae gene function (Table 2; Supplementary Table 5, S6). At 5 dpi, downregulated DEGs were primarily associated with mitochondrion and hydrolase activity, whereas at 20 dpi, upregulated DEGs were enriched in cytosol and cytoplasm (Figures 3A, B). These patterns suggest a shift in functional activity during infection progression. The top GO terms for up- and downregulated DEGs at each time point are summarized in Figure 3A and Figures 3B, respectively, spanning molecular function, cellular component, and biological process categories. KEGG pathway analysis further revealed dynamic changes in metabolic modules and ion-transport processes under BP treatment (Figures 3C–E; Supplementary Table 7–S9). At 5 dpi, genes involved in the V-type ATPase pathway were predominantly downregulated. At 10 dpi, both up- and downregulated pathways were observed, including metabolic modules such as glycolysis and V-type ATPase. By 20 dpi, transcriptional activity increased markedly, with widespread upregulation of metabolic pathways, including V-type ATPase and glycolysis. Overall, genes associated with the V-type ATPase pathway showed downregulation at early stages and upregulation at later stages, suggesting a time-dependent response under BP treatment.

Table 2. List of GO terms associated with up- and down-regulated genes, showing the top 10 terms by DEG count and at least two genes per term.

RegulationdpiCategoryGO termDMG countP-value
Upregulated5Molecular functionzinc ion binding40.001577812
endonuclease activity20.0000244
nuclease activity20.000194442
ubiquitin-protein transferase activity20.000420177
ubiquitin protein ligase activity20.001262618
10Cellular componentmembrane230.001634744
endoplasmic reticulum140.002486499
cellular_component110.00138483
Molecular functionzinc ion binding160.002562502
mismatched DNA binding30.003683504
acetate-CoA ligase activity20
Biological processbiological_process100.001590311
phospholipid biosynthetic process40.0000799
signaling30.001724481
sporulation resulting in formation of a cellular spore30.001724481
20Cellular componentCytoplasm1610.012354641
Cytosol1440.008864022
Ribosome320.002205132
Golgi apparatus210.032790389
Molecular functionnucleotide binding800.049951026
ATP hydrolysis activity540.014344335
structural constituent of ribosome410.000158804
GTP binding340.00000136
GTPase activity260.0000157
Biological processTranslation460.005211188
Downregulated5Cellular componentMitochondrion80.001810113
extracellular region30.00016989
cell projection20.001429974
Molecular functionhydrolase activity90.000602789
nuclease activity30.000421653
double-stranded RNA binding20.000838472
endonuclease activity20.000838472
Biological processproteasomal protein catabolic process30.000615712
regulation of protein localization to chromatin20
positive regulation of canonical Wnt signaling pathway20.0000453
10Cellular componentU4/U6 x U5 tri-snRNP complex30.004526475
COPII-coated ER to Golgi transport vesicle20.000302933
CIA complex20.000302933
Lsm1-7-Pat1 complex20.000302933
sno(s)RNA-containing ribonucleoprotein complex20.001152103
lipid droplet20.001152103
U6 snRNP20.002738914
Molecular functionoxidoreductase activity, acting on the CH-OH group of donors, NAD or NADP as acceptor20.001152103
Biological processmRNA cis splicing, via spliceosome20.001152103
peptidyl-diphthamide biosynthetic process from peptidyl-histidine20.002738914
20Cellular componentTRAMP complex20.00001
Molecular functionhydrolase activity, acting on ester bonds20.00001
poly(A) RNA polymerase activity20.0000394
Biological processtranscription by RNA polymerase II20.004608625
Functional annotation of differentially expressed genes (DEGs) using gene ontology (GO) and KEGG. **(A)** GO annotation of downregulated DEGs; **(B)** GO annotation of upregulated DEGs; **(C)** KEGG pathway annotation of DEGs at 5, 10, and 20 dpi, respectively. GO terms shown are limited to those with the top 10 highest DEG counts and a minimum of 2 associated genes. dpi = days post-infection.
Figure 3. Functional annotation of differentially expressed genes (DEGs) using gene ontology (GO) and KEGG. (A) GO annotation of downregulated DEGs; (B) GO annotation of upregulated DEGs; (C) KEGG pathway annotation of DEGs at 5, 10, and 20 dpi, respectively. GO terms shown are limited to those with the top 10 highest DEG counts and a minimum of 2 associated genes. dpi = days post-infection.

3.5 Validation

To validate the RNA-seq results, the expression of seven V-type ATPase pathway-related genes was examined using qRT-PCR (Figure 4A). Among these genes (Figures 4B–H), ATPeV1A expression was upregulated at 3 dpi and 10 dpi, followed by downregulation by 20 dpi (Figure 4B). For ATPeV1B, expression was downregulated at most time points, with the exceptions of 5 and 20 dpi. At 0 dpi, it was significantly downregulated (Figure 4C). Similarly, ATPeV1D was consistently downregulated at all time points except for 5 dpi and showed significant downregulation at 0 dpi, just like ATPeV1B (Figure 4D). The expression of ATPeV1E was upregulated during the early stages in the BP+Nc, but it began to show significant downregulation at 5 dpi, a trend that aligns with our transcriptome data (Figure 4E). ATPeV0B showed a similar pattern to ATPeV1D, being downregulated at nearly every time point except 5 dpi (Figure 4F). ATPeV0C was significantly upregulated early on. Consistent with ATPeV1E, its expression was downregulated at 5 dpi and matched the transcriptome trends. Following this, ATPeV0C showed significant downregulation at the late stage of the experiment (Figure 4G). ATPeV0D exhibited a distinct pattern from the other V-type ATPase genes, as it was downregulated at 0, 3, 10, 15, and 20 dpi, with significant downregulation observed at 15 and 20 dpi (Figure 4H).

Structural overview of the V-type ATPase and qPCR validation of selected subunit expression. **(A)** Schematic representation of the V-type ATPase complex. The V1 domain (cytosolic side) comprises subunits **(A, B, E, G)**, etc., and the V0 domain (membrane-embedded) comprises subunits a, c, d, etc. Subunits validated in this study are highlighted in blue, and the subunit used for dsRNA treatment is highlighted in red. The structural model was modified from Chen et al., 2022. Relative expression levels of **(B)** *ATPeV1A* (AAJ76_4700027584), **(C)** *ATPeV1B* (AAJ76_700039738), **(D)** *ATPeV1D* (AAJ76_1600056373), **(E)** *ATPeV1E* (AAJ76_1700024327), **(F)** *ATPeV0B* (AAJ76_1900038297), **(G)** *ATPeV0C* (AAJ76_300040219) and **(H)** *ATPeV0D* (AAJ76_630003653) measured by qRT-PCR at the indicated days post-infection (dpi). The results are expressed as log2 relative gene expression (2<sup>-ΔCt</sup>). Relative gene expression levels are presented as mean ± SD (n = 4 biological replicates). At each time point, differences between the DMSO+Nc and BP+Nc were analyzed using an unpaired two-tailed t-test. *P < 0.05. Triangles indicate genes showing lower expression in the BP+Nc compared with the DMSO+Nc, while dots indicate genes showing the same directional trend as observed in the transcriptome data.
Figure 4. Structural overview of the V-type ATPase and qPCR validation of selected subunit expression. (A) Schematic representation of the V-type ATPase complex. The V1 domain (cytosolic side) comprises subunits (A, B, E, G), etc., and the V0 domain (membrane-embedded) comprises subunits a, c, d, etc. Subunits validated in this study are highlighted in blue, and the subunit used for dsRNA treatment is highlighted in red. The structural model was modified from Chen et al., 2022. Relative expression levels of (B) ATPeV1A (AAJ76_4700027584), (C) ATPeV1B (AAJ76_700039738), (D) ATPeV1D (AAJ76_1600056373), (E) ATPeV1E (AAJ76_1700024327), (F) ATPeV0B (AAJ76_1900038297), (G) ATPeV0C (AAJ76_300040219) and (H) ATPeV0D (AAJ76_630003653) measured by qRT-PCR at the indicated days post-infection (dpi). The results are expressed as log2 relative gene expression (2-ΔCt). Relative gene expression levels are presented as mean ± SD (n = 4 biological replicates). At each time point, differences between the DMSO+Nc and BP+Nc were analyzed using an unpaired two-tailed t-test. *P < 0.05. Triangles indicate genes showing lower expression in the BP+Nc compared with the DMSO+Nc, while dots indicate genes showing the same directional trend as observed in the transcriptome data.

ATPeV0B displayed the most consistent and pronounced downregulation, suggesting a potential association with the response of N. ceranae to BP treatment. Because ATPeV0B forms part of the central membrane ring of the V-ATPase complex and exists as a single-copy gene, it was selected for RNAi analysis to further investigate its functional role during infection.

3.6 RNAi feeding bioassays

In the feeding bioassays under a lower infection dose, dsVATP treatment successfully altered the expression patterns of several V-type ATPase subunits (Figures 5A–G). However, the knockdown effect was not consistently observed at all time points, particularly at 10 and 15 dpi, indicating that the RNAi effect may be transient or limited. As expected, the expression of ATPeV0B showed a consistent decrease across time points. Survival analysis showed that the Blank had the highest, whereas all infected groups showed reduced survival (Figure 5H). The dsVATP+Nc showed a trend toward higher survival compared with the Nc and dsGFP+Nc, although these differences were not statistically significant. Spore counts and genome copy number increased over time in all infection groups. Even though the dsVATP+Nc showed lower pathogen levels at later time points, these differences were not statistically significant (Figures 5I, J). These results suggest that dsVATP treatment may influence the infection dynamics of N. ceranae, although the observed effects were modest.

Relative gene expression of V-type ATPase subunit gene and validation in dsVATP-treated samples. **(A–G)** Relative expression of **(A)** *ATPeV1A* (AAJ76_4700027584), **(B)** *ATPeV1B* (AAJ76_700039738), **(C)** *ATPeV1D* (AAJ76_1600056373), **(D)** *ATPeV1E* (AAJ76_1700024327), **(E)** *ATPeV0B* (AAJ76_1900038297), **(F)** *ATPeV0C* (AAJ76_300040219) and **(G)** *ATPeV0D* (AAJ76_630003653) were measured by qRT-PCR at the indicated days post infection (dpi). The results are expressed as log<sub>2</sub>(2<sup>–ΔCt</sup>), and statistical analyses were performed using Student’s *t*-test. Triangles mark genes significantly downregulated by BP compared with the DMSO + Nc. **(H)** Kaplan–Meier survival curves of *N. ceranae* infection in dsGFP+Nc and dsVATP+Nc. Pairwise log-rank (Mantel–Cox) tests with Bonferroni adjustment were used to compare each treatment with Nc. **(I)** spore counts and **(J)** genome copies were analyzed across time points. Statistical significance among time points was assessed by one-way ANOVA followed by Tukey’s multiple comparisons with Bonferroni adjustment. Data are shown as mean ± SD (n = 4).
Figure 5. Relative gene expression of V-type ATPase subunit gene and validation in dsVATP-treated samples. (A–G) Relative expression of (A) ATPeV1A (AAJ76_4700027584), (B) ATPeV1B (AAJ76_700039738), (C) ATPeV1D (AAJ76_1600056373), (D) ATPeV1E (AAJ76_1700024327), (E) ATPeV0B (AAJ76_1900038297), (F) ATPeV0C (AAJ76_300040219) and (G) ATPeV0D (AAJ76_630003653) were measured by qRT-PCR at the indicated days post infection (dpi). The results are expressed as log2(2–ΔCt), and statistical analyses were performed using Student’s t-test. Triangles mark genes significantly downregulated by BP compared with the DMSO + Nc. (H) Kaplan–Meier survival curves of N. ceranae infection in dsGFP+Nc and dsVATP+Nc. Pairwise log-rank (Mantel–Cox) tests with Bonferroni adjustment were used to compare each treatment with Nc. (I) spore counts and (J) genome copies were analyzed across time points. Statistical significance among time points was assessed by one-way ANOVA followed by Tukey’s multiple comparisons with Bonferroni adjustment. Data are shown as mean ± SD (n = 4).

4 Discussion

This study showed that treating honey bees with BP provided a noticeable benefit against nosemosis caused by N. ceranae. Kaplan-Meier survival analysis showed that BP treatment improved the survival rate of honey bees during infection. Our analysis revealed a discrepancy between spore counts and genome copy numbers at certain time points. While genome copy numbers reflect the total amount of pathogen DNA (including immature or non-viable stages), spore counts measure only the mature, infectious particles. This suggests that BP treatment might have a differential effect on spores and immature or non-viable stages. BP appears to slow down the growth of N. ceranae, potentially by influencing gene expression during the early stage of infection. Specifically, BP treatment was associated with changes in genes related to basic metabolic functions, including those in the V-type ATPase pathway. Such changes may affect energy utilization and internal pH regulation, which are important for spore development and replication. B. pilosa is recognized as a medicinal plant and a beneficial forage source for honey bees. It contains bioactive compounds such as phenolics and flavonoids (25, 26). Several phenolic compounds have been reported to inhibit N. ceranae viability, while flavonoids have also been associated with anti-microsporidian activity in honey bees (50). Additionally, B. pilosa has shown anti-coccidial properties in chickens (27). Consistent with these findings, our results indicate that BP not only improves bee survival but also reduces pathogen load. It should also be noted that DMSO, used as the solvent for BP, exhibited a measurable effect on infection dynamics and may partially confound the interpretation of BP-specific effects. Although DMSO exerted a transient effect, BP treatment was associated with a more sustained change in infection dynamics. However, transcriptomic comparisons were performed between the BP+Nc and DMSO+Nc, which helps to partially isolate BP-associated effects, allowing the observed expression differences to be primarily attributed to BP treatment, although the potential influence of DMSO should still be considered.

The time-dependent changes in N. ceranae gene expression observed in response to BP treatment suggest a dynamic interaction between the host, pathogen, and compounds. Early in the infection, many genes were downregulated, which may reflect the suppression of key biological processes required for N. ceranae infection (12, 32). However, as the infection progressed, this suppression appeared to lessen, and more genes were upregulated. This temporal shift in gene expression suggests that N. ceranae may initiate a compensatory or stress-induced response in the later stages. This may indicate that N. ceranae activates adaptive mechanisms to mitigate the effects of both BP treatment and RNAi exposure. In fungi such as yeast, stress responses often involve the activation of cellular pathways that help maintain energy homeostasis and cellular integrity under environmental stress (51, 52). Although the stress response mechanisms of microsporidia remain poorly understood, a comparable adaptive response may be triggered in N. ceranae to maintain essential cellular functions under treatment-induced stress (53, 54).

This time-dependent pattern implies that BP exerts a stage-specific effect that may inhibit key functions necessary for early infection, while later allowing distinct transcriptional responses as N. ceranae attempts to adapt. Similar time-dependent or stage-specific responses have also been reported in other studies (12, 22). These stage-specific dynamics may highlight potential intervention windows, such as targeting the early downregulated pathways to hinder initial establishment, or interfering with later upregulated responses to limit pathogen recovery. Collectively, these findings suggest that BP does not exert a uniform inhibitory effect throughout N. ceranae infection, but may modulate its gene expression in distinct and dynamic ways depending on the stage of infection.

GO annotation suggests that BP treatment impacts N. ceranae functions in a stage-specific manner. At the early stage of infection (D/B 5 dpi), we observed downregulation of terms such as mitochondrion (GO:0005739) and hydrolase activity (GO:0016787). Although microsporidia lack classical mitochondria, they retain mitosome-derived proteins (55), and hydrolases are involved in nutrient processing and other cellular functions. This pattern suggests that BP may influence energy-related or metabolic processes during early infection, although its direct effect on mitosome components remains unclear. At the late stage (D/B 20 dpi), upregulated GO terms included cytoplasm (GO:0005737) and cytosol (GO:0005829). Typically, cellular formation and replication processes are activated early in microsporidia development. However, their delayed upregulation here may indicate that BP treatment alters the timing of these processes. The cytoplasm in microsporidia contains numerous ribosomes essential for protein synthesis (56), and during sporogony, binary fission leads to sporoblast formation, followed by karyokinesis and cytokinesis (57). The delayed expression of cytoplasm-related genes may therefore be associated with changes in replication dynamics, consistent with the observed delay in infection progression. Altogether, the GO results suggest that BP may influence early-stage metabolic functions and be associated with delay cell division during infection, potentially contributing to slower N. ceranae development. The observed downregulation of the V-type ATPase pathway provides a potential link to how BP affects N. ceranae physiology. As V-type ATPases are crucial for maintaining intracellular homeostasis and proton transport, changes in their expression may influence conditions required for spore formation and replication. Based on the KEGG analysis, BP treatment was associated with time-dependent changes in metabolic activity. In the early stage of infection, modules involved in energy generation and biosynthesis, such as V-type ATPase (M00160), nucleotide synthesis, glycolysis (M00001), were downregulated. This may reflect reduced ATP-dependent ion transport and metabolic activity during early infection. At 20 dpi, several of these pathways, including the V-type ATPase pathway, were upregulated again, suggesting a possible compensatory response by N. ceranae under BP-induced stress. V-type ATPase functions in proton transport across membranes via ATP hydrolysis (58), and may be involved in ion transport and intracellular homeostasis in microsporidia (59). The early downregulation of V-type ATPase-related genes may affect proton transport and ion balance, while the later upregulation may represent an adaptive response to restore cellular equilibrium. Overall, the KEGG results indicate that BP treatment is associated with an initial reduction in metabolic and ion transport activity, followed by a later transcriptional rebound, suggesting a dynamic response of N. ceranae to BP treatment across infection stages.

To further evaluate the effect of BP treatment on N. ceranae infection and its potential impact on the V-type ATPase pathway, we analyzed the expression levels of seven key V-type ATPase subunits identified as DEGs. This analysis aimed to assess the potential functional relevance of the V-type ATPase pathway in the context of infection and treatment. The timing of this expression peaks roughly coincided with increases in spore counts and genome copy numbers, suggesting that V-type ATPase activity may be associated with N. ceranae development. Although a causal relationship remains to be confirmed, the observed expression patterns are consistent with this possibility, particularly given the established role of V-type ATPases in energy metabolism, ion homeostasis, and intracellular trafficking in other eukaryotic systems (60, 61). V-type ATPases are multi-subunit pumps involved in intracellular acidification and various cellular processes. Some subunits contribute to structural stability, while others form the proton-conducting channel. For example, the V0 proteolipid ring has been shown to function as a membrane-spanning pore (62). Based on its central role in proton transport and its consistent expression pattern, ATPeV0B was selected for RNAi-mediated knockdown to explore its potential function during infection. Our qPCR results showed that ATPeV0B expression was consistently lower under BP treatment. Although ATPeV1D exhibited a similar tread, ATPeV0B was prioritized due to its role in the V0 domain, which is directly involved in proton transport, whereas ATPeV1D primarily contributes to the catalytic V1 domain (60). The infection dose used in the RNAi experiment (5×10³ spores per bee) was intentionally lower than that used in the BP efficacy experiment (1×105 spores per bee) to avoid excessive mortality and allow longitudinal observation of dsRNA treatment effects. While RNAi-mediated knockdown of ATPeV0B was associated with a partial reduction in N. ceranae load and a modest improvement in host survival, these differences were not statistically significant. Therefore, although the results suggest a possible involvement of ATPeV0B in infection dynamics, they do not provide definitive evidence of its functional role. Instead, these findings support the idea that the V-type ATPase pathway may be relevant to parasite development and warrants further investigation.

While the dsRNA treatment specifically targeted ATPeV0B, several V-type ATPase-related genes showed expression changes similar to those observed under BP treatment. This suggests that BP may affect multiple pathways, including the V-type ATPase pathway, in N. ceranae. The upregulation of the ATPeV0C subunit, a major component of the c-ring that functions together with ATPeV0B (62), may represent a compensatory response to changes in ATPeV0B expression. This observation suggests that N. ceranae may exhibit a degree of transcriptional flexibility under treatment conditions.

Although RNAi treatment was designed to target ATPeV0B, the observed effects on gene expression and pathogen load were modest and not statistically significant. Therefore, while changes in related subunits were observed, these patterns should be interpreted with caution, and do not provide clear evidence of a compensatory mechanism. Instead, they may reflect broader transcriptional responses to treatment. Despite these limitations, the altered expression of multiple V-type ATPase-related genes suggests that perturbation of this pathway may influence overall cellular function. Although dsVATP treatment was associated with a reduction in pathogen load, the differences were not statistically significant, which may be due to factors such as insufficient dsRNA dosage or transient knockdown effects. These findings suggest that targeting a single gene or pathway may not be sufficient to achieve a strong therapeutic effect against a complex pathogen. In contrast, BP treatment, which likely affects multiple pathways simultaneously, may exert broader effects on infection dynamics.

Generative AI statement

The author(s) declared that generative AI was used in the creation of this manuscript. Generative AI was used to assist with grammar checking and minor language editing.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/finsc.2026.1828903/full#supplementary-material

The standard curve of genome copies. Supplementary material at the source

The position of dsVATP on ATPeV0B. Supplementary material at the source

Primer list. Supplementary material at the source

RNA sequencing summary. Supplementary material at the source

Up-regulated DEG list of D/B comparison group. Supplementary material at the source

Down-regulated DEG list of D/B comparison group. Supplementary material at the source

Up-regulated GO list of D/B comparison group. Supplementary material at the source

Down-regulated GO list of D/B comparison group. Supplementary material at the source

Up-regulated KEGG list of D/B comparison group. Supplementary material at the source

Down-regulated DEG list of D/B 20 dpi comparison group. Supplementary material at the source

KEGG module number of the D/B comparison group. Supplementary material at the source

Footnotes

Edited by: Terd Disayathanoowat, Chiang Mai University, Thailand

Reviewed by: Eui-Joon Kil, Andong National University, Republic of Korea

Wenfeng Li, Guangdong Provincial Academy of Chinese Medical Sciences, China

Statements

Conflict of interest

The author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Funding

The author(s) declared that financial support was received for this work and/or its publication. This research was supported by 113-2313-B-005-024 -MY3 from the National Science and Technology Council (NSTC).

Data availability

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics approval

The manuscript presents research on animals that do not require ethical approval for their study.

Author contributions

YL: Writing – original draft, Conceptualization. KC: Conceptualization, Investigation, Writing – original draft. JC: Investigation, Software, Writing – review & editing. ZC: Investigation, Writing – review & editing. CK: Methodology, Writing – review & editing. SH: Methodology, Writing – review & editing. LC: Methodology, Writing – review & editing. MW: Methodology, Writing – review & editing. WY: Methodology, Writing – review & editing. YY: Methodology, Writing – review & editing. YC: Methodology, Writing – review & editing. YN: Writing – review & editing, Methodology.

Acknowledgements

We thanks to Fang-Min Chang for the experiment support. An additional edit was made and confirmed by the corresponding author.

References

  1. 1 Evans JD, Schwarz RS. Bees brought to their knees: microbes affecting honey bee health. Trends Microbiol. (2011) 19: 614–20. doi: doi:10.1016/j.tim.2011.09.003
  2. 2 Tauber JP, Collins WR, Schwarz RS, Chen Y, Grubbs K, Huang Q. Natural product medicines for honey bees: perspective and protocols. Insects. (2019) 10: 356. doi: doi:10.3390/insects10100356
  3. 3 Goulson D, Nicholls E, Botías C, Rotheray EL. Bee declines driven by combined stress from parasites, pesticides, and lack of flowers. Science. (2015) 347: 1255957. doi: doi:10.1126/science.1255957
  4. 4 Potts SG, Biesmeijer JC, Kremen C, Neumann P, Schweiger O, Kunin WE. Global pollinator declines: trends, impacts and drivers. Trends Ecol Evol. (2010) 25: 345–53. doi: doi:10.1016/j.tree.2010.01.007
  5. 5 Li Y-H, Chen Y-H, Chang F-M, Wu M-C, Nai Y-S. Monitoring the season–prevalence relationship of Vairimorpha ceranae in honey bees (Apis mellifera) over one year and the primary assessment of probiotic treatment in Taichung, Taiwan. Insects. (2024) 15: 204. doi: doi:10.3390/insects15030204
  6. 6 Gisder S, Schüler V, Horchler LL, Groth D, Genersch E. Long-term temporal trends of Nosema spp. infection prevalence in Northeast Germany: continuous spread of Nosema ceranae, an emerging pathogen of honey bees (Apis mellifera), but no general replacement of Nosema apis. Front Cell Infect Microbiol. (2017) 7: 301. doi: doi:10.3389/fcimb.2017.00301
  7. 7 Jabal-Uriel C, Barrios L, Bonjour-Dalmon A, Caspi-Yona S, Chejanovsly N, Erez T. Epidemiology of the microsporidium Nosema ceranae in four Mediterranean countries. Insects. (2022) 13: 844. doi: doi:10.3390/insects13090844
  8. 8 Chen Y, Evans JD, Zhou L, Boncristiani H, Kimura K, Xiao T. Asymmetrical coexistence of Nosema ceranae and Nosema apis in honey bees. J Invertebr Pathol. (2009) 101: 204–9. doi: doi:10.1016/j.jip.2009.05.012
  9. 9 Martín‐Hernández R, Bartolomé C, Chejanovsky N, Le Conte Y, Dalmon A, Dussaubat C. Nosema ceranae in Apis mellifera: a 12 years postdetection perspective. Environ Microbiol. (2018) 20: 1302–29. doi: doi:10.1111/1462-2920.14103
  10. 10 Goblirsch M. Nosema ceranae disease of the honey bee (Apis mellifera). Apidologie. (2018) 49: 131–50. doi: doi:10.1007/s13592-017-0535-1
  11. 11 Higes M, Juarranz Á, Dias‐Almeida J, Lucena S, Botías C, Meana A. Apoptosis in the pathogenesis of Nosema ceranae (Microsporidia: Nosematidae) in honey bees (Apis mellifera). Environ Microbiol Rep. (2013) 5: 530–6. doi: doi:10.1111/1758-2229.12059
  12. 12 Holt HL, Aronstein KA, Grozinger CM. Chronic parasitization by Nosema microsporidia causes global expression changes in core nutritional, metabolic and behavioral pathways in honey bee workers (Apis mellifera). BMC Genomics. (2013) 14: 799. doi: doi:10.1186/1471-2164-14-799
  13. 13 Mayack C, Naug D. Individual energetic state can prevail over social regulation of foraging in honeybees. Behav Ecol Sociobiol. (2013) 67: 929–36. doi: doi:10.1007/s00265-013-1517-6
  14. 14 Iorizzo M, Letizia F, Ganassi S, Testa B, Petrarca S, Albanese G. Recent advances in the biocontrol of Nosemosis in honey bees (Apis mellifera L.). J Fungi. (2022) 8: 424. doi: doi:10.3390/jof8050424
  15. 15 Huang W-F, Solter LF, Yau PM, Imai BS. Nosema ceranae escapes fumagillin control in honey bees. PloS Pathog. (2013) 9: e1003185. doi: doi:10.1371/journal.ppat.1003185
  16. 16 Glavinic U, Jovanovic NM, Dominikovic N, Lakic N, Ćosić M, Stevanovic J. Potential of wormwood and oak bark-based supplement in health improvement of Nosema ceranae -infected honey bees. Animals. (2024) 14: 1195. doi: doi:10.3390/ani14081195
  17. 17 Glavinic U, Blagojevic J, Ristanic M, Stevanovic J, Lakic N, Mirilovic M. Use of thymol in Nosema ceranae control and health improvement of infected honey bees. Insects. (2022) 13: 574. doi: doi:10.3390/insects13070574
  18. 18 Tlak Gajger I, Ribarić J, Smodiš Škerl M, Vlainić J, Sikirić P. Stable gastric pentadecapeptide BPC 157 in honeybee (Apis mellifera) therapy, to control Nosema ceranae invasions in apiary conditions. J Vet Pharmacol Ther. (2018) 41: 614–21. doi: doi:10.1096/fasebj.31.1_supplement.1068.13
  19. 19 Glavinic U, Stankovic B, Draskovic V, Stevanovic J, Petrovic T, Lakic N. Dietary amino acid and vitamin complex protects honey bee from immunosuppression caused by Nosema ceranae. PloS One. (2017) 12: e0187726. doi: doi:10.1371/journal.pone.0187726
  20. 20 Ho ST, Nai YS, Chang ZT, Chang JC, Hsu WC, Ko CY. Dimethyl sulfoxide, an alternative for control of Nosema ceranae infection in honey bees (Apis mellifera). Arch Insect Biochem Physiol. (2024) 116: e22099. doi: doi:10.1002/arch.22099
  21. 21 Chaimanee V, Kasem A, Nuanjohn T, Boonmee T, Siangsuepchart A, Malaithong W. Natural extracts as potential control agents for Nosema ceranae infection in honeybees, Apis mellifera. J Invertebr Pathol. (2021) 186: 107688. doi: doi:10.1016/j.jip.2021.107688
  22. 22 Özkırım A, Küçüközmen B. Application of herbal essential oil extract mixture for honey bees (Apis mellifera L.) against Nosema ceranae and Nosema apis. J Apicultural Sci. (2021) 65 (1): 163–75. doi: doi:10.2478/JAS-2021-0010
  23. 23 Cannizzaro C, Keller A, Wilson RS, Elliott B, Newis R, Ovah R. Forest landscapes increase diversity of honeybee diets in the tropics. For Ecol Manage. (2022) 504: 119869. doi: doi:10.1016/j.foreco.2021.119869
  24. 24 Hsu P-S, Wu T-H, Huang M-Y, Wang D-Y, Wu M-C. Nutritive value of 11 bee pollen samples from major floral sources in Taiwan. Foods. (2021) 10: 2229. doi: doi:10.3390/foods10092229
  25. 25 Xuan TD, Khanh TD. Chemistry and pharmacology of Bidens pilosa: an overview. J Pharm Invest. (2016) 46: 91 – 132. doi: doi:10.1007/s40005-016-0231-6
  26. 26 Bartolome AP, Villaseñor IM, Yang W-C. Bidens pilosa L.(Asteraceae): botanical properties, traditional uses, phytochemistry, and pharmacology. Evidence-Based Complementary Altern Med. (2013) 2013. doi: doi:10.1155/2013/340215
  27. 27 Yang W-C, Yang C-Y, Liang Y-C, Yang C-W, Li W-Q, Chung C-Y. Anti-coccidial properties and mechanisms of an edible herb, Bidens pilosa, and its active compounds for coccidiosis. Sci Rep. (2019) 9: 1 – 11. doi: doi:10.1038/s41598-019-39194-2
  28. 28 Oliveira F, Andrade-Neto V, Krettli A, Brandão M. New evidences of antimalarial activity of Bidens pilosa roots extract correlated with polyacetylene and flavonoids. J Ethnopharmacol. (2004) 93: 39 – 42. doi: doi:10.1016/j.jep.2004.03.026
  29. 29 Brandão M, Krettli A, Soares L, Nery C, Marinuzzi H. Antimalarial activity of extracts and fractions from Bidens pilosa and other Bidens species (Asteraceae) correlated with the presence of acetylene and flavonoid compounds. J Ethnopharmacol. (1997) 57: 131–8. doi: doi:10.1016/S0378-8741(97)00060-3
  30. 30 Yang W, Tien Y, Chung C, Chen Y, Chiou W, Hsu S. Effect of Bidens pilosa on infection and drug resistance of Eimeria in chickens. Res Vet Sci. (2015) 98: 74 – 81. doi: doi:10.56397/jimr/2022.07.04
  31. 31 Huang W-F, Jiang J-H, Chen Y-W, Wang C-H. A Nosema ceranae isolate from the honeybee Apis mellifera. Apidologie. (2007) 38: 30–7. doi: doi:10.1051/apido:2006054
  32. 32 Li Y-H, Chang Z-T, Yen M-R, Huang Y-F, Chen T-H, Chang J-C. Transcriptome of Nosema ceranae and upregulated microsporidia genes during its infection of western honey bee (Apis mellifera). Insects. (2022) 13: 716. doi: doi:10.1007/978-3-030-80836-5_6
  33. 33 Ko B, Van Raamsdonk JM. RNA sequencing of pooled samples effectively identifies differentially expressed genes. Biology. (2023) 12: 812. doi: doi:10.3390/biology12060812
  34. 34 Takele Assefa A, Vandesompele J, Thas O. On the utility of RNA sample pooling to optimize cost and statistical power in RNA sequencing experiments. BMC Genomics. (2020) 21: 312. doi: doi:10.1186/s12864-020-6721-y
  35. 35 Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. (2014) 30: 2114–20. doi: doi:10.1093/bioinformatics/btu170
  36. 36 Schmieder R, Edwards R. Quality control and preprocessing of metagenomic datasets. Bioinformatics. (2011) 27: 863–4. doi: doi:10.1093/bioinformatics/btr026
  37. 37 Moriya Y, Itoh M, Okuda S, Yoshizawa AC, Kanehisa M. KAAS: an automatic genome annotation and pathway reconstruction server. Nucleic Acids Res. (2007) 35: W182–5. doi: doi:10.1093/nar/gkm321
  38. 38 Nino-Sanchez J, Chen LH, De Souza JT, Mosquera S, Stergiopoulos I. Targeted delivery of gene silencing in fungi using genetically engineered bacteria. J Fungi (Basel). (2021) 7 (2): 125. doi: doi:10.3390/jof7020125
  39. 39 Wu C-F, Chang L-P, Lee C, Stergiopoulos I, Chen L-H. pSIG plasmids, MoClo-compatible vectors for efficient production of chimeric double-stranded RNAs in Escherichia coli HT115 (DE3) strain. Plant Methods. (2025) 21: 96. doi: doi:10.1186/s13007-025-01413-5
  40. 40 He N, Zhang Y, Duan XL, Li JH, Huang W-F, Evans JD. RNA interference-mediated knockdown of genes encoding spore wall proteins confers protection against Nosema ceranae infection in the European honey bee, Apis mellifera. Microorganisms. (2021) 9: 505. doi: doi:10.3390/microorganisms9030505
  41. 41 Kim IH, Kim DJ, Gwak WS, Woo SD. Increased survival of the honey bee Apis mellifera infected with the microsporidian Nosema ceranae by effective gene silencing. Arch Insect Biochem Physiol. (2020) 105: e21734. doi: doi:10.1002/arch.21734
  42. 42 Kassambara A, Kosinski M, Biecek P, Fabian S. Package ‘survminer’. Drawing survival curves using ‘ggplot2’ (R package version 03 1). (2017).
  43. 43 Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2-ΔΔCT method. Methods. (2001) 25: 402–8. doi: doi:10.1006/meth.2001.1262
  44. 44 Schmittgen TD, Livak KJ. Analyzing real-time PCR data by the comparative CT method. Nat Protoc. (2008) 3: 1101–8. doi: doi:10.1038/nprot.2008.73
  45. 45 Team RC. R: A language and environment for statistical computing (R Version 4.0. 3, R Foundation for Statistical Computing, Vienna, Austria, 2020). (2021). doi: doi:10.32614/r.manuals
  46. 46 Kolde R, Kolde MR. Package ‘pheatmap’. R package. (2015) 1: 790.
  47. 47 Wickham H, Navarro D, Pedersen T. ggplot2: elegant graphics for data analysis. New York, NY: Springer-Verlag (2016).
  48. 48 Kassambara A, Kassambara MA. Package ‘ggpubr’. (2020).
  49. 49 Moore DF. Applied survival analysis using R. Cham, Switzerland: Springer (2016).
  50. 50 Burnham AJ. Scientific advances in controlling Nosema ceranae (Microsporidia) infections in honey bees (Apis mellifera). Front Vet Sci. (2019) 6: 79. doi: doi:10.3389/fvets.2019.00079
  51. 51 Gasch AP, Spellman PT, Kao CM, Carmel-Harel O, Eisen MB, Storz G. Genomic expression programs in the response of yeast cells to environmental changes. Mol Biol Cell. (2000) 11: 4241–57. doi: doi:10.1091/mbc.11.12.4241
  52. 52 Święciło A. Cross-stress resistance in Saccharomyces cerevisiae yeast—new insight into an old phenomenon. Cell Stress Chaperones. (2016) 21: 187 – 200. doi: doi:10.1007/s12192-016-0667-7
  53. 53 Han B, Takvorian PM, Weiss LM. Invasion of host cells by microsporidia. Front Microbiol. (2020) 11: 172. doi: doi:10.3389/fmicb.2020.00172
  54. 54 Tamim El Jarkass H, Reinke AW. The ins and outs of host‐microsporidia interactions during invasion, proliferation and exit. Cell Microbiol. (2020) 22: e13247. doi: doi:10.1111/cmi.13247
  55. 55 Burri L, Williams BA, Bursac D, Lithgow T, Keeling PJ. Microsporidian mitosomes retain elements of the general mitochondrial targeting system. Proc Natl Acad Sci. (2006) 103: 15916–20. doi: doi:10.1073/pnas.0604109103
  56. 56 Vávra J, Larsson JR. Structure of microsporidia. In: Microsporidia: Pathogens of Opportunity. Ames, IA: Wiley-Blackwell (2014). p. 1 – 70.
  57. 57 Bigliardi E, Sacchi L. Cell biology and invasion of the microsporidia. Microbes Infection. (2001) 3: 373–9. doi: doi:10.1016/s1286-4579(01)01393-4
  58. 58 Beyenbach KW, Wieczorek H. The V-type H+ ATPase: molecular structure and function, physiological roles and regulation. J Exp Biol. (2006) 209: 577–89. doi: doi:10.1242/jeb.02014
  59. 59 Chetia H, Kabiraj D, Sharma S, Bora U. Comparative insights to the transportome of Nosema: a genus of parasitic microsporidians. bioRxiv. (2017), 110809. doi: doi:10.1101/110809
  60. 60 Vasanthakumar T, Rubinstein JL. Structure and roles of V-type ATPases. Trends Biochem Sci. (2020) 45: 295 – 307. doi: doi:10.1016/j.tibs.2019.12.007
  61. 61 Banerjee S, Kane PM. Regulation of V-ATPase activity and organelle pH by phosphatidylinositol phosphate lipids. Front Cell Dev Biol. (2020) 8: 510. doi: doi:10.3389/fcell.2020.00510
  62. 62 Couoh-Cardel S, Hsueh Y-C, Wilkens S, Movileanu L. Yeast V-ATPase proteolipid ring acts as a large-conductance transmembrane protein pore. Sci Rep. (2016) 6: 24774. doi: doi:10.1038/srep24774

Citation

Li, Y., et al. (2026). Regulation of Nosema ceranae gene expression by Bidens pilosa phytogenic treatment in infected honey bees. Frontiers in Insect Science, Volume 6. https://doi.org/
Received
Accepted
Published
Volume
6 - 2026

Copyright

Copyright © 2026 Li, Chen, Chang, Chang, Ko, Ho, Chen, Wu, Yang, Yang, Chen and Nai.

License: https://creativecommons.org/licenses/by/4.0/

*Correspondence: Yu-Shin Nai

Article metrics

Views
0
Downloads
0
Citations
0
Shares
0

Related articles

Original ResearchSeasonal gut microbiota dynamics and colony performance under stationary and short-distance migratory beekeeping in Apis mellifera ligusticaMassimo Iorizzo, Gianluca Albanese, Cosimo Tedino, Mariantonietta Succi, Licia Maria Di Donato, Antonio de Nigris, Dalila Di Criscio, Sonia Ganassi, Luigi Ambrosone, Antonio De CristofaroOriginal ResearchBiocontrol potential of hive-associated Bacillus species against American Foulbrood disease in honeybees (Apis mellifera L.)Anuja Shrestha, Yva Eline, Kenzie Sundstrom, James T. Van LeuvenBrief Research ReportCumulative adult production and aerobic culture-dependent recovery of gut-associated bacteria in Galleria mellonella maintained on polymer-containing dietsStephany Eileen Keferstein-Borges, Gabriela Alcalá-Gómez, Claudia Ivette Cisneros-Reyes, Jaime Alcalá-Gómez