Transcriptome of Nosema ceranae and Upregulated Microsporidia Genes during Its Infection of Western Honey Bee (Apis mellifera)

Simple Summary In this study, the gene expression profile of a honey bee microsporidium, Nosema ceranae, was investigated at 5, 10, and 20 days post-infection. Based on transcriptome data, the common DEGs and stage-specific genes were identified and validated. Our data reveal that the gene expression of N. ceranae during infection is highly related to the mechanisms of gene transcription, protein synthesis, and structural proteins. This information could be a reference for the control of nosemosis at the genetic level in apiculture. Abstract Nosema ceranae is one of the fungal parasites of Apis mellifera. It causes physical and behavioral effects in honey bees. However, only a few studies have reported on gene expression profiling during A. mellifera infection. In this study, the transcriptome profile of mature spores at each time point of infection (5, 10, and 20 days post-infection, d.p.i.) were investigated. Based on the transcriptome and expression profile analysis, a total of 878, 952, and 981 differentially expressed genes (DEGs) (fold change ≥ 2 or ≤ −2) were identified in N. ceranae spores (NcSp) at 5 d.p.i., 10 d.p.i., and 20 d.p.i., respectively. Moreover, 70 upregulated genes and 340 downregulated genes among common DEGs (so-called common DEGs) and 166 stage-specific genes at each stage of infection were identified. The Gene Ontology (GO) analysis indicated that the DEGs and corresponding common DEGs are involved in the functions of cytosol (GO:0005829), cytoplasm (GO:0005737), and ATP binding (GO:0005524). Furthermore, the pathway analysis found that the DEGs and common DEGs are involved in metabolism, environmental information processing, and organismal systems. Four upregulated common DEGs with higher fold-change values, highly associated with spore proteins and transcription factors, were selected for validation. In addition, the stage-specific genes are highly involved in the mechanism of pre-mRNA splicing according to GO enrichment analysis; thus, three of them showed high expression at each d.p.i. and were also subjected to validation. The relative gene expression levels showed a similar tendency as the transcriptome predictions at different d.p.i., revealing that the gene expression of N. ceranae during infection may be related to the mechanism of gene transcription, protein synthesis, and structural proteins. Our data suggest that the gene expression profiling of N. ceranae at the transcriptomic level could be a reference for the monitoring of nosemosis at the genetic level.


Introduction
Microsporidia are considered primitive eukaryotic and obligate intracellular parasites; they have 70 S ribosomes of small subunit (SSU) and large subunit (LSU) rRNA, which is similar to eukaryotes, but they lack mitochondria, peroxisomes, and centrioles. Nonetheless,

-Honey bees
Worker honey bees were collected from A. mellifera colonies in apiaries at the National Chung Hsing University (NCHU, Taiching, Taiwan). Each honey bee colony had six to eight frames and the presence of a normal spawning queen. Brood cells of honeycombs from the six colonies, which were capped at more than 70-80%, were selected and incubated at 34 • C until honey bees emerged. The newly emerged adult honey bees were collected in a rearing cage and supplied with 50% syrup for the following experiments.
-Artificial infection of Nosema ceranae Mature N. ceranae spores were collected from the midgut tissues of honey bees from apiaries at the NCHU (National Chung Hsing University, Taichung, Taiwan) and purified by the protocol proposed by Huang et al., 2007 [34]. Briefly, the midgut tissues of honey bees were collected in 500 µL 1 × TE buffer (10 mM Tris, 1 mM EDTA, pH 7.5) (Bioman, New Taipei, Taiwan) and homogenized with a pestle. The homogenized tissues were then filtered with 100 mesh stainless steel wire cloth. The filtrate was centrifuged at 1000× g for 10 min, and the supernatant and tissue debris were removed. Then, the pellet was resuspended in 500 µL 1 × TE buffer. The wash step was repeated to remove any remaining tissue debris. The spore suspension was centrifuged at 16,000× g for 3 min at 4 • C, and the collected spores were counted with a haemocytometer (Paul Marienfeld GmbH & Co. KG, Baden-Württemberg, Germany). The concentration of the mature spore suspension was adjusted to 1 × 10 5 spore/2 µL for the infection of honey bees [35].
In the N. ceranae-infected group, the newly emerged honey bees were fed 2 µL of 50% sugar syrup with N. ceranae spores (1 × 10 5 spores per honey bee), while the honey bees in the control group were fed 2 µL of 50% sugar syrup. Each group contained fifty honey bees in one rearing cage, which were supplied with 50% syrup and incubated at 34 • C until sampling. Samples for both RNA-seq and the detection of N. ceranae infections were selected in triplicate, and qPCR validation for gene expression was repeated in four replicates. The survival rates of the infected and control groups were observed and recorded daily. A Kaplan-Meier survival curve analysis was performed using SPSS (Norman H. Nie, C. Hadlai (Tex) Hull and Dale H. Bent, IBM SPSS Statistics version 20, Palo Alto, CA, USA). -

Confirmation of Nosema ceranae infection
Before the experiment, multiplex PCR [36] was performed to confirm the absence of Nosema spp. infection in newly emerged adult honey bees (Supplementary Figure S1). To confirm the N. ceranae infection, three midgut tissues of N. ceranae-infected bees and noninfected honey bees were collected in 1 × TE buffer. The midgut tissues were homogenized with a pestle, and the tissue suspension was subjected to observation, spore counting, and DNA extraction. Mature spores were observed under 400× light microscopy (WHITED, Taiwan) and counted by a haemocytometer (Paul Marienfeld GmbH & Co. KG, Baden-Württemberg, Germany). A total of 100 µL of homogenized tissue suspension was subjected to DNA extraction by the EDNA HISPEX (Easy DNA High-Speed Extraction) tissue kit (Fisher Biotec, Stirling, Australia) according to the user manual. The DNA sample was then subjected to PCR detection with the microsporidia universal primer set 18f/1537r, and the Hb18f/Hb18r primer set was used as an internal control (Supplementary Table S1). For the PCR, each 20 µL PCR contained 10 µL 2× SuperRed PCR Master mix (BIOTOOLs, New Taipei, Taiwan), 1 µL 100 mM of each primer, and 1 µL template DNA. PCR amplification was performed by a MiniAmp™ Thermal Cycler (Thermo Scientific, Waltham, MA, USA) with the following program: the thermal cycle was preheated to 95 • C for 10 min, followed by 35 cycles of 95 • C for 15 s, 55 • C for 30 s, and 72 • C for 1 min 30 s, and then a final extension was performed of 72 • C for 10 min and storage at 20 • C. The PCR product was analyzed by electrophoresis on a 1% agarose gel in 1 × Tris-acetate-EDTA (TAE) buffer.

-Preparation of RNA samples for next-generation sequencing
Based on the data of the survival curve and mature spore counts, midgut tissues from three honey bees were collected at 0 (3 h), 5, 10, and 20 days post-infection (d.p.i.) and purified N. ceranae mature spores (1 × 10 8 ) were homogenized by a Tissue Lyser II (QIAGEN, Solingen, Germany) at f = 30/s for 1 min. This procedure was repeated at three time points. The total RNA was extracted using TRIzol reagent (Invitrogen Life Technologies, Waltham, MA, USA) following the user's manual. The extracted RNA sample was then purified using a Locking Gel (BIOTOOLS, New Taipei, Taiwan). The RNA purity and integrity were checked by a Nanodrop 1000 Spectrophotometer V3.5 (Thermo Scientific, Waltham, MA, USA) and Qubit 2.0 Fluorometer (Invitrogen Life Technologies, Waltham, MA, USA), respectively. Total RNA samples of three biological repeats were pooled and used for library construction. The size distribution of each RNA sample was determined by Fragment Analyzer Systems (Agilent, Santa Clara, CA, USA). The RNA sample was then used to construct a library for transcriptome sequencing. -

RNA-Seq library construction and sequencing
The mRNA from each sample was further purified for sequencing library construction using an Ultra RNA Library Prep Kit (NEB, Essex, MA, USA) according to the manufacturer's protocol. Briefly, poly-A-containing mRNA was purified by oligo d (T) magnetic beads. During the final step of mRNA purification, RNA was fragmented and primed with NEBNext random primers for cDNA synthesis. After mRNA was reverse-transcribed, end-repair and adaptor ligation (NEBNext Adaptor) were performed. Finally, DNA fragments were selectively enriched by PCR and then qualified and quantified by Fragment Analyser Systems (Agilent, Santa Clara, CA, USA) and qPCR, respectively, for the DNA library quality control analysis. The libraries were then pooled and sequenced by the Illumina MiSeq sequencer with a paired-end (PE) technology of 301 bp in length to generate high-throughput transcriptome sequencing data. Sequencing adaptors were trimmed using Trimmomatic [37] from raw PE reads, and then a read quality check was performed by PRINseq [38] from trimmed PE reads. To confirm the pathogen load in the sequencing data of the newly emerged adult honey bees, the trimmed PE reads were mapped to the genomes of Nosema apis and 6 other viruses.
-Differentially expressed gene (DEG) analysis The raw PE reads were trimmed by Trimmomatic [37], and then the clean reads were obtained. The quality reads were mapped to the N. ceranae genome index (ASM98816v1) using HISAT2 with default parameters, and the corresponding gene annotation information (GCF_000988165.1_ASM98816v1_genomic.gff) from the National Biotechnology Information Center (NCBI) (https://www.ncbi.nlm.nih.gov/ (accessed on 19 January 2021)) was used to obtain the N. ceranae-related PE reads for further analysis. Differentially expressed genes (DEGs) among experimental groups were analyzed via Cuffdiff, and the gene expression heatmap was generated by the R package [39] (Supplementary Figure S2). For the DEG analysis, genes with a zero FPKM value from each group of the analysis were filtered out. To facilitate logarithmic transformation, a value of 1 was added to the FPKM value of the remaining genes. Raw signals (FPKM t 1) were then transformed using log2-based transformation. A statistical test of fold changes was performed, and a significant difference in the expression of these two groups was defined as a log2 ratio ≥ 1 or ≤ −1 (fold change ≥ 2 or ≤ −2). The DEGs from each comparison (NcSp/5 d.p.i., NcSp/10 d.p.i., and NcSp/20 d.p.i.) were further compared by Venny 2.1.0 (Juan Carlos Oliveros, Madrid, Spain) [40] to obtain the common DEGs of N. ceranae in different infection stages. The DEGs and common DEGs were further used in the GO enrichment analysis and submitted to the KEGG Automatic Annotation Server (KAAS, http://www.genome.jp/tools/kaas/ (accessed on 23 April 2021)) for gene orthologue assignment and pathway mapping [41]. The selected N. ceranae common DEGs were further analyzed by EMB BLAST (European Molecular Biology Network-Swiss node) for a functional homology search. The results with an E-value < 10 −4 were recorded for further discussion of functional proteins. -

qRT-PCR validation
For validation, three midgut tissues of infected honey bees were collected into 300 µL of TE buffer at 0.5 (12 h), 1, 3, 5, 10, 15, and 20 d.p.i., with four replicates, and were homogenized for RNA extraction. Four replicates of 1 × 10 8 purified mature N. ceranae spores were also subjected to RNA extraction. The total RNA for qRT-PCR validation was extracted from 200 µL of homogenized midgut tissue and the N. ceranae mature spore sample using the GENEzol™ TriRNA Pure Kit (Geneaid, New Taipei, Taiwan) as described in the user manual. The concentration and quantity of the total RNA were measured using a Nanodrop™ spectrophotometer (Thermo Scientific, Waltham, USA). Each RNA sample was adjusted to 1 µg/mL for cDNA synthesis by BioDiamond RT-PCR mix (Origin Pure BioSci & Tech.Co., Yilan, Taiwan) as described in the user manual. Four highly upregulated common DEGs, AAJ76_1600052943, AAJ76_2000149161, AAJ76_2000141845, and AAJ76_500091146, and three highly expressed stage-specific genes, AAJ76_5000026485, AAJ76_2400020190, and AAJ76_2300035467, were selected for RT-qPCR validation. β-Tubulin from N. ceranae served as the internal control. Real-time PCR was performed by iQ™ SYBR ® Green Supermix (Bio-Rad, Hercules, CA, USA) in a CFX Connect Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA). qRT-PCR was performed as follows: 95 • C for 3 min, followed by 40 cycles of 95 • C for 10 s, 60 • C for 30 s, and 72 • C for 45 s. The relative gene expression levels for common DEGs were calculated by the 2-Ct method, and stage-specific genes were calculated by the 2-Ct method [42]. -

Statistical analysis
Statistical analysis and plotting were performed with R 4.0.3 [43] and the heatmap, ggplot2, and ggpubr [39,44,45] packages. The statistical analysis of the Kaplan-Meier survival curves was performed using the R packages survminer [46] and survival [47]. The statistical analysis of the spore count and qRT-PCR data was performed using one-way ANOVA with SPSS ® [48].

N. Ceranae Infection
The Kaplan-Meier survival curve results showed that survival fitness was significantly reduced (p < 0.01) in the N. ceranae-infected A. mellifera ( Figure 1A). Mature spores were observed from 5 to 20 d.p.i. in the homogenized midgut tissues of the N. ceranae-infected honey bees under microscopy; while only a few mature spores were observed at 5 d.p.i., the spore counts significantly increased from 10 to 20 d.p.i ( Figure 1B,C). Molecular detection was performed to detect the infection process of N. ceranae on A. mellifera. Although the earliest time that mature spores were observed was at 5 d.p.i., the results showed that N. ceranae infection was first detected at 3 d.p.i. by PCR ( Figure 1D).

Sequencing Data Summary
The sequencing data are summarized in Supplementary Table S2 The mappability of each sample to the N. ceranae genome was 0%, 4%, 3%, 3%, and 12% for 0 (3 h), 5, 10, and 20 d.p.i. and mature N. ceranae spores, respectively. In addition, the newly emerged adult honey bees were confirmed to be free of N. apis and common honey bee viral infections (Supplementary Table S3). These mapped reads were subjected to further analysis of the expression of different genes.

DEG Analysis
For the DEG analysis, a total of 878, 952, and 981 genes were identified to have significant fold changes (fold change ≥ 2 or ≤−2) in NcSp/5 d.p.i., NcSp/10 d.p.i., and NcSp/20 d.p.i., respectively. Among these DEGs, 429, 117, and 113 genes were upregulated and 449, 835, and 868 genes were downregulated in NcSp/5 d.p.i., NcSp/10 d.p.i., and NcSp/20 d.p.i., respectively ( Figure 2A; Tables 1, 2, and S4-S9). The common DEGs were selected from all three groups of common DEGs by set intersection, and stage-specific genes were defined as the genes that were only expressed at 5, 10, and 20 d.p.i., respectively (Tables 3 and S10). There was a total of 70 upregulated and 340 downregulated common DEGs ( Figure 2B; Tables 1, 2, and S4-S9), and the fold-change distribution at each infection stage is shown in Figure 2B. The results also showed that a greater upregulation of N. ceranae common DEGs and stage-specific genes was found in the initial infectious stage, whereas more genes of N. ceranae were suppressed in the late stage of infection ( Figure 2). For the stage-specific genes, a total of 25, 71, and 70 genes were identified at each time post-infection (Supplementary Table S10), and the top 10 expression values (FPKM) of the stage-specific genes are summarized in Table 3.

KEGG Analysis
The results of the KEGG analysis show the pathway annotation of the DEGs and common DEGs at different stages and indicate the mechanisms of N. ceranae infection ( Figure 5 and Supplementary Tables S17-S22). In the comparison of NcSp at 5 d.p.i., there were 211 pathways that included upregulated genes and 200 pathways that included downregulated genes ( Figure 5A). Among these pathways, highly abundant N. ceranae DEGs were involved in the metabolism (KO:09100), environmental information processing (KO:09130), and organismal systems (KO:09150) categories, including the carbohydrate metabolism (KO:09101), signal transduction (KO:09132), and immune system (KO:09151) pathways ( Figure 5A and Supplementary Tables S17 and S18). In the NcSp at 10 d.p.i. group, 112 pathways included upregulated genes and 1224 pathways included downregulated genes. Among these pathways, three categories, including the metabolism pathway (KO:09100), human diseases (KO:09160), and environmental information processing (KO:09130) were the most mapped, and they contained the pathways of carbohydrate metabolism (KO:09101) and signal transduction (KO:09132) ( Figure 5B; Supplementary Tables S19 and S20). For the NcSp/20 d.p.i. group, there were 129 pathways that included upregulated genes and 267 pathways that included downregulated genes, and the pathway categories that highly regulated genes belonged to, including metabolism pathways (KO:09100), environmental information processing (KO:09130), human diseases (KO:09160), and organismal systems (KO:09150), were identified. These pathways also included carbohydrate metabolism (KO:09101) and signal transduction (KO:09132) ( Figure 5C; Supplementary Tables S21 and S22). The KEGG analysis results of the common DEGs were similar to those of the DEGs ( Figure 5D; Supplementary Tables S23 and S24). A total of 80 pathways included upregulated genes and 199 pathways included downregulated genes during N. ceranae infection ( Figure 5D). A high number of downregulated common DEGs were involved in the pathways of signal transduction (KO:09132); protein families: genetic information processing (KO:09182) and carbohydrate metabolism (KO:09101) within the environmental information process (KO:09130); brite hierarchies (KO:09180); and metabolism (KO:09100) categories.

Validation
Four highly upregulated common DEGs (AAJ76_1600052943, AAJ76_2000149161, AAJ76_2000141845, and AAJ76_500091146) and three highly expressed stage-specific genes (AAJ76_5000026485, AAJ76_2400020190, and AAJ76_2300035467) at each time postinfection were selected to validate the RNA-seq results during infection by N. ceranae. The results showed that all of the upregulated common DEGs increased at 3 d.p.i. and reached a high peak at 15 to 20 d.p.i., while no or a low abundance of gene expression was detected at 0.5 to 1 d.p.i. (Figure 6). The validation result matched the detection of N. ceranae infection, which was detectable at 3 d.p.i. (Figure 1D), suggesting that N. ceranae begins the infection process and initiates gene expression at 3 d.p.i. In addition, the upregulated common DEGs, reaching a high peak at 15-20 d.p.i., also showed coordination with the increasing spore count of N. ceranae ( Figure 1C). For the stage-specific gene validations, the stage-specific genes were only expressed at specific times post-infection, while the gene expression of AAJ76_5000026485, AAJ76_2400020190, and AAJ76_2300035467 was detected at all times post-infection, including in mature spores ( Figure 7). However, the results of three stage-specific genes showed higher expression levels at 5, 10, and 20 d.p.i., revealing a similar expression trend as the transcriptomic prediction ( Figure 7).

Discussion
The pattern of N. ceranae infection in this study was similar to a previous report that N. ceranae was detectable in an early infection stage at 3 d.p.i., while mature spores were only observed at 5 d.p.i. [18,49]. Negative impacts on the survival fitness of N. ceranae-infected honey bees have been reported, and they explain the changes to immune suppression, nutritional deficiencies, and hormone regulation [18,[50][51][52]. These data indicate that infected N. ceranae might begin to affect and have consistent negative effects on honey bees at earlier infectious stages. Altogether, these data reveal that N. ceranae infections can be detected as early as 3 d.p.i. and the infection level dramatically increases from 10 to 20 d.p.i.
The gene expression profiling of N. ceranae was described at the transcriptomic level in this study. The results indicated that a greater upregulation of N. ceranae's common and stage-specific genes was found in the initial infectious stage, whereas more genes of N. ceranae were suppressed in the late stage of infection. This tendency was similar to the study of N. ceranae gene expression in Apis cerana, suggesting that the N. ceranae genes were expressed at an early infection stage, followed by the gene expression profile of N. ceranae switching from the merogony to sporogony stages [53]. Microsporidia lack mitochondria, and the loss of most metabolic pathways might involve an alteration of gene expression during the course of infection; therefore, host-generated ATP would be imported to support their metabolism [54,55]. Several highly expressed stage-specific genes (AAJ76_5000026485, AAJ76_2400020190, AAJ76_2300035467, and AAJ76_600074061) were identified as pre-mRNA splicing factors, indicating that post-transcriptional regulation was continually activated during N. ceranae infection.
For the GO annotation of common DEGs, cytosol (GO:0005829) and cytoplasm (GO:0005737) were also found in the upregulated common DEG group. The cytoplasm mainly comprises microsporidia, containing many ribosomes in a regular array [56]. During the replication of sporonts, division by binary fission produces two sporoblasts and infects the hosts through karyokinesis and cytokinesis [2]. Therefore, the upregulation of genes belonging to the cytoplasm with corresponding GO terms might contribute to the early propagation of N. ceranae. At the initiation of microsporidial infection, sporoplasm was injected into the cytoplasm of the host cell by the germination of the polar tube [2]. Because of the lack of genes for many metabolic pathways in the microsporidia, the microsporidia acquire nutrients from the host. Microsporidia in contact with the host cytoplasm could provide more opportunities to export cytotoxic compounds across the plasma membrane as an alternative to compensate for development [57]. In addition, upregulated genes belonging to the ATP-binding (GO:0005524) pathway were found in the early stage of N. ceranae infection (NcSp/5 d.p.i.). In the ATP-binding (GO:0005524) pathway, the transmembrane protein superfamily of microsporidia, which includes ATP-binding cassette (ABC) transporters, plays an important role in substrate transportation and pathogen development, indicating that energy transport is abundant in this stage [57].
It was noted that metal ion binding (GO:0046872), GTP binding (GO:0005525), and the ubiquitin-dependent protein catabolic process (GO:0006511) showed high gene counts as observed in the downregulated common DEGs. Metal ions are vital elements for proper folding and the ribozyme catalysis of nucleic acids [58][59][60]. GTP-binding proteins regulate cellular processes for multiple protein biosynthesis and intracellular membrane traffic [61]. The role of the ubiquitin-proteasome pathway has been reported to be the selective degradation of soluble cellular proteins under most conditions [62]. The downregulation of the genes belonging to the GO term indicate the requirement of nucleotide and protein bases from host cells in N. ceranae infections. In contrast, many stage-specific genes of N. ceranae were involved in pre-mRNA-associated pathways, including "RNA binding" (GO:0003723) and "mRNA splicing, via spliceosome" (GO:0000398), suggesting the impacts on the transcriptional level for N. ceranae infection.
Based on the result of KEGG pathway analysis on DEGs, the categories of metabolism and environmental information processing included a high number of downregulated genes. Moreover, in the pathways of metabolism, the "carbohydrate metabolism" included a high number of downregulated genes at each infection stage. It has been reported that Nosema acquires carbohydrates from the epithelial cells of the honey bee gut lining, leading to a series of metabolic changes in the host [3,20]; therefore, this point was also proven by the downregulation of genes belonging to "carbohydrate metabolism".
Instead of mitochondria, microsporidia contain a genome-less organelle called a mitosome, which encodes only 20 mitosomal proteins in microsporidia, indicating the loss of function of the mitosome organ compared to that of normal mitochondria in the yeast Saccharomyces cerevisiae [32,63]; therefore, the mature spores could still maintain basic metabolism by their own system. In the beginning of infection, genes belonging to the "folding, sorting, and degradation (K09123)" and "replication and repair (K09124)" pathways in the genetic information processing (K09120) categories, the cellular processes pathway (KO:09140) of "cell growth and death" (K09143), and the environmental information processing pathway (KO:09130) of "signal transduction" (K09132) were highly upregulated in the NcSp/5 d.p.i. group. Meanwhile, the genes involved in these pathways were downregulated at a late stage of N. ceranae infection (NcSp/10 d.p.i. and NcSp/20 d.p.i. groups), suggesting that the merogony stage needs to initiate more functions of genetic and information processes than mature spores to facilitate the propagation of N. ceranae and maintain the activities of infection. Lastly, from our data, the downregulated genes involved in the pathway of carbohydrate metabolism suggest that N. ceranae takes energy from honey bees during infection [8].
The validation of four highly upregulated common DEGs and three highly expressed stage-specific genes at each time post-infection showed similar tendencies to those of the DEG predictions. The AAJ76_1600052943 and AAJ76_2000141845 showed a functional homology to the spore wall proteins SWP25 and SWP30, respectively. In Nosema bombycis, the two spore wall proteins have been characterized as an endosporal protein (SWP30) and an exosporal protein (SWP32) [64,65]. These two spore wall proteins bind to the chitin component of N. bombycis, implying that these proteins are the structural support proteins in the spore wall. SWP25 was first found in N. bombycis; the localization of SWP25 is the endospore, and this protein possesses a signal peptide and a heparin-binding motif. In addition, SWP25 is critical for host-to-host transmission and persistence in the environment [66]. Furthermore, SWP25 has been reported to be related to spore wall integrity and germination, and recent research has shown that the induction of SWP25 expression after binding with haematopoietin may change the integrity of pathogens and affect the clearance of immune cells or the transmission of pathogens in haemolymph [65,67]. Therefore, the upregulation of AAJ76_1600052943 and AAJ76_2000141845 might also play roles similar to those mentioned above to facilitate spore germination and mature spore development. The highly expressed SPW was consistent with our previous study on N. ceranae gene expression in A. mellifera [33] and the research on the N. ceranae gene expression in A. cerana, suggesting a significant role for SPW during N. ceranae infection [53].
AAJ76_2000149161 was mapped to a DNA-dependent protein kinase catalytic subunit (DNA-PKcs) by SWISS-MODEL [68]. DNA-PKcs is a key member of the phosphatidylinositol-3 kinase-like (PIKK) family of protein kinases and plays an essential role in the repair of DNA double-strand breaks by nonhomologous end joining (NHEJ) [69]. The upregulation of this gene might contribute to the relative function of N. ceranae in using host DNA fragments for its replication. AAJ76_500091146 is described as a forkhead hnf3 transcription factor based on protein functional prediction. HNF-3/FKH genes are a large family of transcriptional activators, and all forkhead proteins contain a highly conserved DNA binding domain and play numerous regulatory roles in eukaryotes. At least four genes in yeast are involved in a wide range of biological functions, such as apoptosis, cell cycle, immunity, development, growth, stress resistance, metabolism, and reproduction [70,71]. Based on the above functions, the upregulated AAJ76_500091146 might be involved in multiple regulatory roles in the development and growth of N. ceranae.
Interestingly, the functional prediction of all three stage-specific genes indicated posttranscriptional modification. AAJ76_5000026485 is a functional homologue of the pre-mRNA-splicing factor cwf14, which is a type of spliceosome that catalyzes the excision of introns from pre-mRNA [72]. A study from yeast showed that cwf14 functions in the regulation of cell growth and mRNA splicing and is also a key factor regulating cell wall biogenesis [73]. The function of AAJ76_2400020190 was found to match that of the putative H/ACA ribonucleoprotein complex subunit 3, which is known as a small nuclear RNP (snRNP). snRNP predominantly guides the site-directed pseudouridylation of target RNAs, processes ribosomal RNA, and stabilizes telomerase RNA [74]. Moreover, the function of H/ACA ribonucleoproteins is essential for ribosome biogenesis, pre-mRNA splicing, and telomere maintenance [74]. AAJ76_2300035467 is functionally similar to U6 snRNAassociated Sm-like protein LSm6. The Sm proteins are RNA-protein complex proteins, and LSm6 is an Sm-like protein (LSM) in S. cerevisiae. A study of S. cerevisiae suggested that LSm6 is associated with U6 snRNA and is required for the maintenance of normal U6 snRNA levels and for pre-mRNA splicing [75]. In summary, the primary transcription products (pre-mRNAs) undergo a variety of post-transcriptional modifications to create functional mRNAs [76]. These three stage-specific genes were identified at different infection stages, and all associations with pre-mRNA splicing suggest that they are highly involved in the activation of post-transcriptional modifications of N. ceranae at each infectious stage.

Conclusions
In this study, the transcriptomic analysis results indicated alterations in the N. ceranae gene expression status during A. mellifera infection. The identification of the common and stage-specific genes expressed during the infection process revealed the interaction between N. ceranae and A. mellifera. Downregulated genes during the infection, especially at a late period of infection, were found. Regarding RT-qPCR validation, the common DEGs showed higher expression levels from 5 d.p.i. to 20 d.p.i., indicating the importance of common DEGs. In our study, we attempted to identify the DEGs resulting from gene expression profiling based on transcriptomic analysis and proposed the common and stage-specific genes of N. ceranae during its infection of A. mellifera. Our data also describe the relatively comprehensive gene expression of the N. ceranae infection stage and analyze the details during infection, thus providing basic information on the N. ceranae mechanism during infection and furthering the application, including the early detection of nosemosis.
Supplementary Materials: The following supporting information can be downloaded at: https://www. mdpi.com/article/10.3390/insects13080716/s1, Figures S1 and S2; Tables S1-S24. Figure S1. The detection of Nosema spp. of newly emerged adult honey bees by multiplex PCR before experiment. The information of primer sets were listed in Supplementary Tables S1; Figure S2. The pipeline of transcriptome analysis. The commend lines and analysis purposes are noted pairwise the item. DEG = differential expressed gene; FPKM = fragment per kilobase per million; GO = gene ontology; KEGG = Kyoto Encyclopedia of Genes and Genomes; Table S1. Primer sets list; Table S2. RNA sequencing summary; Table S3. Transcriptomic data for honey bee infected with Nosema ceranae were screened by mapping to the genomes of Nosema apis and 6 honey bee infecting viruses;