Comparative Identification of MicroRNAs in Apis cerana cerana Workers’ Midguts in Response to Nosema ceranae Invasion

Here, the expression profiles and differentially expressed miRNAs (DEmiRNAs) in the midguts of Apis cerana cerana workers at 7 d and 10 d post-inoculation (dpi) with N. ceranae were investigated via small RNA sequencing and bioinformatics. Five hundred and twenty nine (529) known miRNAs and 25 novel miRNAs were identified in this study, and the expression of 16 predicted miRNAs was confirmed by Stem-loop RT-PCR. A total of 14 DEmiRNAs were detected in the midgut at 7 dpi, including eight up-regulated and six down-regulated miRNAs, while 12 DEmiRNAs were observed in the midgut at 10 dpi, including nine up-regulated and three down-regulated ones. Additionally, five DEmiRNAs were shared, while nine and seven DEmiRNAs were specifically expressed in midguts at 7 dpi and 10 dpi. Gene ontology analysis suggested some DEmiRNAs and corresponding target mRNAs were involved in various functions including immune system processes and response to stimulus. KEGG pathway analysis shed light on the potential functions of some DEmiRNAs in regulating target mRNAs engaged in material and energy metabolisms, cellular immunity and the humoral immune system. Further investigation demonstrated a complex regulation network between DEmiRNAs and their target mRNAs, with miR-598-y, miR-252-y, miR-92-x and miR-3654-y at the center. Our results can facilitate future exploration of the regulatory roles of miRNAs in host responses to N. ceranae, and provide potential candidates for further investigation of the molecular mechanisms underlying eastern honeybee-microsporidian interactions.


Introduction
Honeybees play vital roles not only in the pollination of crops and wild flora, but also in the support of critical ecosystem balance [1,2]. In addition, honeybees serve as key models for studying development, social behavior, and disease transmission [3,4]. Western honeybees (Apis mellifera) and eastern honeybees (Apis cerana), the two best-known honeybee species, have been used for honey production and crop pollination. Compared with A. mellifera, A. cerana has some prominent advantages, such as long flight duration, adaptation to extreme weather conditions, frequent grooming and hygienic behaviors, and cooperative group-level defenses [5]. Apis cerana cerana, a subspecies of A. cerana, is an important endemic species that is crucial to ecosystem balance and environmental improvement in China [6].
Nosemosis is a serious disease in adult honeybees caused by several Nosema species. Nosema infection occurs through the ingestion of spores in contaminated food or water, followed by the

Small RNA Isolation, cDNA Library Construction and Deep Sequencing
Firstly, total RNA of each midgut sample in N. ceranae-treated and control groups were extracted using TRIzol Reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer's protocols. Secondly, DNA contaminants were removed with RNase-freeDNase I (TaKaRa, Beijing, China). The purified RNA quantity and quality were checked using Nanodrop 2000 spectrophotometer (Thermo Fisher, Waltham, MA, USA), and the integrity of the RNA samples were evaluated using Agilent 2100 bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) and only values of 28S/18S ≥ 0.7 and RIN ≥ 7.0 were considered qualified for the subsequent small RNA library construction. Thirdly, RNA molecules in the size range of 18-30 nt were enriched by agarose gel electrophoresis (AGE) and then ligated with 3 and 5 RNA adaptors, and fragments with adaptors on both ends were enriched by PCR after reverse transcription. Fourthly, the subsequent cDNAs were purified and enriched by 3.5% AGE to isolate the expected size (140-160 bp) fractions and eliminate unincorporated primers, primer dimer products, and dimerized adaptors. Ultimately, the 12 cDNA libraries were sequenced on Illumina sequencing platform (HiSeq TM 4000) using the single-end technology by GENE DENOVO Biotechnology Co. (Guangzhou, China). The libraries were as follows: AcCK1-1, AcCK1-2 and AcCK1-3 as replicate libraries for normal midguts of workers at 7 dpi with sucrose solution; AcT1-1, AcT1-2 and AcT1-3 as replicate libraries for midguts of workers at 7 dpi with sucrose solution containing N. ceranae spores; AcCK2-1, AcCK2-2 and AcCK2-3 as replicate libraries for normal midguts of workers at 10 dpi with sucrose solution; AcT2-1, AcT2-2 and AcT2-3 as replicate libraries for midguts of workers at 10 dpi with sucrose solution containing N. ceranae spores.

Quality Control and Sequencing Data Analysis
The raw data were pre-processed to exclude low-quality reads (length < 20 nt and ambiguous N), 5 adapter, 3 adapter and poly(A) sequences to gain clean reads, which were aligned against NCBI GeneBank and Rfam databases to remove ncRNA including rRNA, scRNA, snoRNA, snRNA and tRNA. Then, the obtained sequences were compared with exons and introns in the A. cerana genome (assembly ACSNU-2.0) to classify mRNA degradation products and the repeat associate miRNA sequences. All the downstream analyses were based on clean reads with high quality.
Currently, only the miRNA information of western honeybee was recorded in miRBase database, while the related information of eastern honeybees, including A. c. cerana, was not included. By utilizing Bowtie (v1.1.0) [38], the filtered sequences were analyzed with BLAST search against miRBase 21.0 by allowing at most two mismatches outside of the seed region [39], and small RNAs that matched exist miRNAs of other animal species in miRBase were identified as known miRNAs. The sequences that did not match known miRNAs were used to identify potentially novel miRNA candidates using RNAfold software [40]. Only sequences with typical Stem-loop hairpins and free energy lower than -20 kcal/mol were considered as potential novel miRNAs. The suffixes "-x" and "-y" respectively mean a certain miRNA deriving from the processing of the 5 and 3 arms of its precursor, while the suffix "-z" means a certain miRNA with unknown processing direction. Size distribution and saturation analysis were performed using sequencing libraries after all annotation steps. Hierarchical clustering of novel miRNA expression was conducted using heatmap tool in OmicShare (http://www.omicshare.com/tools/).

Analysis of DEmiRNAs
The miRNA expression levels in each sample were normalized to the total number of sequence tags per million (TPM) following formula: normalized expression = mapped read count/total reads×10 6 . Differential expression analysis of two samples was performed using the DEGseq R package [41] and p values were adjusted using q value. The criteria of p value (FDR) < 0.05 and |log 2 (Fold change)| > 1 were set as the threshold for statistically significant differential expression. A positive value indicated upregulation of a miRNA, while a negative value indicated down-regulation.

Prediction of the Target mRNAs of DEmiRNAs and Construction of miRNA-mRNA Regulation Networks
Target mRNAs for miRNAs were predicted with miRanda (v3.3a) [42], RNAhybrid (v2.1.2) +svm_light (v6.01) [43] and TargetFinder (Version: 7.0) [44] software. The input files were miRNA FASTA sequences files. Intersections of the results from the three above mentioned programs comprised the final predicted gene targets. Then, miRNA-mRNA regulation networks were visualized using Cytoscape v.3.2.1 software [45], following parameters: the free energy of the target site/miRNA duplex needs to be lower than −35 kcal/mol and p < 0.05.

GO and KEGG Pathway Analyses of DEmiRNA Target mRNAs
To determine their main biological functions, DEmiRNA target mRNAs were annotated by terms in the Gene Ontology (GO) database (http://www.geneontology.org/) using Blast2GO [46], following their numeric orders in the NCBI nr database. Subsequently, GO functional categorizations for all target mRNAs were determined with WEGO software [47], and GO enrichment analysis of functional significance terms in the GO database was performed using the hypergeometric test to identify significantly enriched GO terms for the target mRNAs compared to the genome background. KEGG pathway analyses of the predicted target mRNAs were performed using the KEGG pathway database (http://www.genome.jp/kegg/pathway.html) [48].

Stem-Loop RT-PCR Confirmation of miRNAs
Total RNAs from AcCK1 and AcCK2 samples were isolated using RNAiso plus kit (TaKaRa) and then treated with DNase I (TaKaRa) to remove remaining DNA. According to the previously described method [49], Stem-loop primers, specific forward primers and universal reverse primers (presented in Table S2) were designed using DNAMAN software on basis of the sequences of the randomly selected nine A. c. cerana miRNAs, including miR-1943-x, miR-3793-x, miR-252-y, miR-3963-x, miR-8516-x, miR-149-y, miR-7311-y, miR-9008-x, novel-m0014-3p. The primers were synthesized by Sangon Biotech Co., Ltd. (Shanghai, China). One microgram of total RNA was reverse transcribed to cDNA using RevertAid First Strand cDNA Synthesis Kit (TaKaRa) and Stem-loop primers (Table S2). The PCR amplification of randomly selected miRNAs was conducted on a T100 thermo cycler (BIO-RAD) using Premix (TaKaRa) under the following conditions: pre-denaturation step at 94 • C for 5 min; 30 amplification cycles of denaturation at 94 • C for 50 s, annealing at 55 • C for 30 s, and elongation at 72 • C for 1 min; followed by a final elongation step at 72 • C for 10 min. The PCR products were detected on 2% AGE with Genecolor (Gene-Bio, Shenzhen, China) staining.

RT-qPCR Validation of DEmiRNAs
Four DEmiRNAs (miR-6547-x, miR-7311-x, novel-m0008-3p and miR-2779-y) within AcCK1 vs AcT1, and four DEmiRNAs (miR-3726-x, miR-1788-y, miR-3319-y and miR-1672-x) within AcCK2 vs AcT2 were randomly selected for RT-qPCR validation. Total RNA of two control midgut samples (AcCK1, AcCK2) and two N. ceranae-infected midgut samples (AcT1, AcT2) were respectively isolated using RNAiso plus kit (TaKaRa), followed by inverse transcription with corresponding Stem-loop primers using the method mentioned above. The generated cDNA was used for RT-qPCR validation with specific forward primer and universal reverse primer. Stem-loop primers and specific forward primers (presented in Table S2) were designed and synthesized using the aforementioned method. RT-qPCR was carried out in an Applied Biosystems QuantStudio 3 (Thermo Fisher, Waltham, MA, USA) in a 20 µL reaction volume containing 1 µL cDNA, 10 µL SYBR Premix (Vazyme, Nanjing, China), 0.2 µL specific forward primer (20 µM), 0.2 µL reverse primer (20 µM), and 8.6 µl DEPC water. The reaction was performed at 95 • C for 5 min, followed by 45 cycles of 94 • C for 15 s, 60 • C for 15 s, and 72 • C for 15 s. The abundance of miRNAs was normalized relative to that of endogenous control snRNA U6. All reactions were performed in triplicate. The threshold cycle (Ct) was determined using the default threshold settings, and the data were analyzed using 2 -∆∆Ct program [50]. The experiment was performed three times using three independent biological samples.

Statistical Data Analysis
All statistical analyses were performed using SPSS software (IBM, Armonk, NY, USA) and GraphPad Prism 6.0 software (GraphPad Software Inc., San Diego, CA, USA). Data were presented as mean ± standard deviation (SD). Statistical analysis was calculated using independent-samples t-test and one-way ANOVA. Fisher's exact test was employed to filter the significant GO terms and KEGG pathways using R software 3.3.1. p < 0.05 was considered statistically significant.

Overview of the Small RNA Library from A. c. cerana Workers' Guts
To identify the miRNA profiles and DEmiRNAs during the responses of A. c. cerana workers' midguts to N. ceranae invasion, 12 small RNA libraries were constructed using Illumina sequencing. In total, 127,523,419 raw reads with an average of~106,269,512 reads per sample were produced, and after filtering out low-quality sequences, 5 and 3 adaptors and reads 18 < nt, 122,104,443 clean reads were obtained for further analysis (Table 1). In addition, the Pearson correlation between every sample in each group was above 0.9619 ( Figure S2), suggesting sufficient reproducibility and rationality of sampling. A BLAST run against the NCBI GenBank and RFam databases identified 3,611,375 (36.20%) to 5,054,110 (52.00%) unique small RNAs as rRNA, 101,602 (1.02%) to 221,407 (1.97%) as tRNAs, 5123 (0.05%) to 11,137 (0.11%) as snRNAs, and 176 (0.002%) to 584 (0.005%) as snoRNAs (Table S3).
Base composition, which is a fundamental feature of miRNA sequences, influences miRNA physiochemical and biochemical properties [51,52]. The sequence lengths of the clean tags were almost all distributed between 16-30 nt, with peaks at 22 nt accounting for 14.32% of all clean reads, followed by 20 nt and 21 nt, which accounted for 13.79% and 13.98% of clean reads, respectively ( Figure 1A). The nucleotide bias of miRNAs in the midgut of the A. c. cerana worker was further investigated, and the results showed that U was the most common nucleotide at the 5 end of miRNAs ( Figure 1B), which was consistent with other studies that also found U to be the most common base at the extreme 5 end of miRNAs [53][54][55]. Additionally, analysis of the nucleotide bias at each position showed that U and A were mainly located at the beginning and end of reads of A. c. cerana miRNAs ( Figure 1C), suggesting that AU base pairing may affect miRNA secondary structure or target recognition [56]. (1.97%) as tRNAs, 5123 (0.05%) to 11,137 (0.11%) as snRNAs, and 176 (0.002%) to 584 (0.005%) as snoRNAs (Table S3). Base composition, which is a fundamental feature of miRNA sequences, influences miRNA physiochemical and biochemical properties [51,52]. The sequence lengths of the clean tags were almost all distributed between 16-30 nt, with peaks at 22 nt accounting for 14.32% of all clean reads, followed by 20 nt and 21 nt, which accounted for 13.79% and 13.98% of clean reads, respectively ( Figure 1A). The nucleotide bias of miRNAs in the midgut of the A. c. cerana worker was further investigated, and the results showed that U was the most common nucleotide at the 5′ end of miRNAs ( Figure 1B), which was consistent with other studies that also found U to be the most common base at the extreme 5′ end of miRNAs [53][54][55]. Additionally, analysis of the nucleotide bias at each position showed that U and A were mainly located at the beginning and end of reads of A. c. cerana miRNAs ( Figure 1C), suggesting that AU base pairing may affect miRNA secondary structure or target recognition [56].   To explore known miRNAs and novel miRNAs in the four groups (AmCK1, AmT1, AmCK2, AmT2), the mapped sequences were compared to reference genomes and aligned with known mature miRNAs in the miRBase 21.0 database. In total, 340, 338, 306, and 302 known miRNAs were identified from the AmCK1, AmT1, AmCK2 and AmT2 samples, respectively. These known miRNAs had a broad range of expression levels in A. c. cerana, ranging from TPM 90513.14 to TPM 0.38. Among them, bantam-y, miR-184-y, miR-1-y, miR-276-y and miR-750-y were the most abundant known miRNAs in both AcCK1 and AcCK2 (Table S4). We also predicted 25 novel miRNAs not previously found in A. c. cerana, including 25, 22, 23 and 17 in AmCK1, AmT1, AmCK2 and AmT2 (Table S5), respectively. Among these, novel-m0019-3p, novel-m0001-3p, novel-m0012-5p, novel-m0004-3p and novel-m0012-3p were the highest-expressed novel miRNAs (Table S5). Moreover, expression clustering analysis showed that these novel miRNAs had various expression levels in different groups ( Figure 2). Interestingly, 242 known miRNAs and 23 novel miRNAs were shared by AmCK1 and AmCK2, implying their developmental stage-specific functions. Secondary structures of precursors of three known miRNAs (miR-7-x, miR-9895-y and miR-750-y) and three novel miRNAs (novel-m0001-3p, novel-m0004-3p and novel-m0019-5p) are displayed in Figure 3, and they all have classical Stem-loop structures.    clustering analysis showed that these novel miRNAs had various expression levels in different groups ( Figure 2). Interestingly, 242 known miRNAs and 23 novel miRNAs were shared by AmCK1 and AmCK2, implying their developmental stage-specific functions. Secondary structures of precursors of three known miRNAs (miR-7-x, miR-9895-y and miR-750-y) and three novel miRNAs (novel-m0001-3p, novel-m0004-3p and novel-m0019-5p) are displayed in Figure 3, and they all have classical Stem-loop structures. Furthermore, nine miRNAs were selected for Stem-loop RT-PCR validation, the results showed that eight of them yielded signal bands (approximately 60 bp) in both the AcCK1 and AcCK2 samples ( Figure 4, see also Figure S3), indicating that most predicted miRNAs were expressed in the A. c. cerana worker's midgut. However, more than one fragment was amplified from miR-1943-x and miR-252-y, suggestive of a non-specific amplification ( Figure 4, see also Figure S3). Furthermore, nine miRNAs were selected for Stem-loop RT-PCR validation, the results showed that eight of them yielded signal bands (approximately 60 bp) in both the AcCK1 and AcCK2 samples ( Figure 4, see also Figure S3), indicating that most predicted miRNAs were expressed in the A. c. cerana worker's midgut. However, more than one fragment was amplified from miR-1943-x and miR-252-y, suggestive of a non-specific amplification ( Figure 4, see also Figure S3).

Differentially Expressed A. c. cerana miRNAs Induced by N. ceranae Invasion
A total of 14 DEmiRNAs were detected in AcCK1 vs AcT1 (Table S6), including eight upregulated and six downregulated miRNAs, while 12 miRNAs with differential expression were identified in AcCK2 vs AcT2 (Table S7), including nine upregulated and three downregulated ones Negative control (sterile water was used as PCR template); snRNA U6 was used as positive control. The upper panel presents the AGE of PCR amplification products of miR-1943-x, miR-3793-x, miR-252-y, and novel-m0014-3p; the lower panel displays the AGE of PCR amplification products of miR-3963-x, miR-8516-x, miR-149-y, and miR-7311-y. The complete gel was presented in Figure S3.

Differentially Expressed A. c. cerana miRNAs Induced by N. ceranae Invasion
A total of 14 DEmiRNAs were detected in AcCK1 vs AcT1 (Table S6), including eight up-regulated and six downregulated miRNAs, while 12 miRNAs with differential expression were identified in AcCK2 vs AcT2 (Table S7), including nine upregulated and three downregulated ones ( Figure 5A). Venn diagram analysis was performed, and the result suggested that five DEmiRNAs, including miR-60-y and miR-8462-x, were shared in AcCK1 vs AcT1 and AcCK2 vs AcT2, while nine (e.g., miR-1-x and miR-980-y) and seven (e.g., novel-m0003-3p and miR-92-x) DEmiRNAs were unique for the two comparison groups, respectively ( Figure 5B, see also Tables S6 and S7). In addition, more than 85.71% (12) of the DEmiRNAs in AcCK1 vs AcT1 showed a >7-fold difference in gene expression (Table S6), whereas approximately 91.67% (11) of the DEmiRNAs in AcCK2 vs AcT2 displayed a >9-fold differential expression (Table S7). We speculate that these DEmiRNAs are likely to play key roles in the regulation of host responses to N. ceranae infection. Negative control (sterile water was used as PCR template); snRNA U6 was used as positive control. The upper panel presents the AGE of PCR amplification products of miR-1943-x, miR-3793-x, miR-252-y, and novel-m0014-3p; the lower panel displays the AGE of PCR amplification products of miR-3963-x, miR-8516-x, miR-149-y, and miR-7311-y. The complete gel was presented in Figure S3.

Differentially Expressed A. c. cerana miRNAs Induced by N. ceranae Invasion
A total of 14 DEmiRNAs were detected in AcCK1 vs AcT1 (Table S6), including eight upregulated and six downregulated miRNAs, while 12 miRNAs with differential expression were identified in AcCK2 vs AcT2 (Table S7), including nine upregulated and three downregulated ones ( Figure 5A). Venn diagram analysis was performed, and the result suggested that five DEmiRNAs, including miR-60-y and miR-8462-x, were shared in AcCK1 vs AcT1 and AcCK2 vs AcT2, while nine (eg., miR-1-x and miR-980-y) and seven (eg., novel-m0003-3p and miR-92-x) DEmiRNAs were unique for the two comparison groups, respectively ( Figure 5B, see also Tables S6 and S7). In addition, more than 85.71% (12) of the DEmiRNAs in AcCK1 vs AcT1 showed a >7-fold difference in gene expression (Table S6), whereas approximately 91.67% (11) of the DEmiRNAs in AcCK2 vs AcT2 displayed a >9fold differential expression (Table S7). We speculate that these DEmiRNAs are likely to play key roles in the regulation of host responses to N. ceranae infection.

Prediction and Enrichment Analysis of the Target mRNAs of A. c. cerana DEmiRNAs
DEmiRNAs in AcCK1 vs AcT1 and AcCK2 vs AcT2 were respectively used to search A. c. cerana 3'-UTR sequences for predicting potential target mRNAs via a combination of RNAhybrid (v2.1.2) +svm_light (v6.01), Miranda (v3.3a), and TargetScan (v7.0) softwares. In total, 2615 target mRNAs of DEmiRNAs were predicted from AcCK1 vs AcT1, while 1905 DEmiRNA target mRNAs in AcCK2 vs AcT2 were identified. GO classification was conducted, and the results indicated that the target mRNAs of DEmiRNAs in AcCK1 vs AcT1 and AcCK2 vs AcT2 were involved in various biological processes (16 and 15 terms), cellular components (11 and 11 terms), and molecular functions (7 and 7 terms) ( Figure 6A,B, see also Tables S8 and S9). Cellular process was the most enriched class in the biological processes; cell and cell part were highly represented in the cellular component categories, and binding and catalytic activity were the top enriched items in molecular functions (Tables S8 and  S9). Further investigation indicated that a large quantity of the shared DEmiRNA target mRNAs were engaged in binding, cellular process, single-organism process, metabolic process, biological regulation, and catalytic activity (Tables S8 and S9). However, target mRNAs with the function of extracellular matrix and immune system process were related only in AcCK1 vs AcT1, accounting for three DEmiRNAs including miR-1-x, miR-6313-y, and miR-252-y. Supramolecular fiber function was associated only with AcCK2 vs AcT2, accounting for novel-m0003-3p.
KEGG enrichment analysis for DEmiRNA-targeted genes was conducted, and the results demonstrated that the target mRNAs of shared DEmiRNAs in AcCK1 vs AcT1 were connected with 104 pathways; amongst them, the highest enriched pathways were endocytosis, the phosphatidylinositol signaling system, the Wnt signaling pathway, inositol phosphate metabolism, and neuroactive ligand-receptor interaction ( Figure 7A, see also Table S10). The target mRNAs of DEmiRNAs in AcCK2 vs AcT2 were involved in 92 pathways, and the most enriched pathway was the phosphatidylinositol signaling system, followed by the Wnt signaling pathway, phototransduction, endocytosis and the FoxO signaling pathway ( Figure 7B, see also Table S11). Interestingly, 62 and 52 metabolism-related pathways were respectively enriched by 234 and 195 target mRNAs in AcCK1 vs AcT1 and AcCK2 vs AcT2, and, in addition, there were five and six target mRNAs enriched in development in these two comparison groups (Table S10 and Table S11). These results demonstrated that DEmiRNAs could regulate the growth, development and metabolism of A. c. cerana workers' midguts. extracellular matrix and immune system process were related only in AcCK1 vs AcT1, accounting for three DEmiRNAs including miR-1-x, miR-6313-y, and miR-252-y. Supramolecular fiber function was associated only with AcCK2 vs AcT2, accounting for novel-m0003-3p.
KEGG enrichment analysis for DEmiRNA-targeted genes was conducted, and the results demonstrated that the target mRNAs of shared DEmiRNAs in AcCK1 vs AcT1 were connected with 104 pathways; amongst them, the highest enriched pathways were endocytosis, the phosphatidylinositol signaling system, the Wnt signaling pathway, inositol phosphate metabolism, and neuroactive ligand-receptor interaction ( Figure 7A, see also Table S10). The target mRNAs of DEmiRNAs in AcCK2 vs AcT2 were involved in 92 pathways, and the most enriched pathway was the phosphatidylinositol signaling system, followed by the Wnt signaling pathway, phototransduction, endocytosis and the FoxO signaling pathway ( Figure 7B, see also Table S11). Interestingly, 62 and 52 metabolism-related pathways were respectively enriched by 234 and 195 target mRNAs in AcCK1 vs AcT1 and AcCK2 vs AcT2, and, in addition, there were five and six target mRNAs enriched in development in these two comparison groups (Table S10 and Table S11). These results demonstrated that DEmiRNAs could regulate the growth, development and metabolism of A. c. cerana workers' midguts.

Regulatory Networks between A. c. cerana DEmiRNAs and Their Target mRNAs
Cytoscape was employed for visualization of the regulatory networks between DEmiRNA and their target mRNAs. A single mRNA can be targeted by various miRNAs, whereas a single miRNA is also capable of targeting different mRNAs [57]. In these regulatory networks, it can be clearly seen that DEmiRNAs and their target mRNAs in AcCK1 vs AcT1 formed complex networks. MiR-598-y had as many as 45 target mRNAs, and miR-252-y bound to three targets including XM_017049297.1, XM_017049298.1 and XM_017058024.1 ( Figure 8A). In contrast, miR-6313-y, miR-3726-x, miR-9204-x and miR-6717-x could bind to only one target ( Figure 8A). In addition, as shown in Figure 8B, for DEmiRNAs in AcCK2 vs AcT2, miR-92-x, miR-3654-y, novel-m0003-3p and miR-6313-y linked to six, four, two and one target mRNAs, respectively.

Verification of A. c. cerana DEmiRNAs via RT-qPCR
To verify the high-throughput sequencing data in our study, the expression levels of nine randomly selected DEmiRNAs in AcCK1 vs AcT1 and AcCK2 vs AcT2 were examined by RT-qPCR. The expression trends of eight DEmiRNAs in the two comparison groups were roughly the same as the result of the sRNA-seq (Figure 9), indicative of the accuracy and reliability of our sequencing data.

Regulatory Networks between A. c. cerana DEmiRNAs and Their Target mRNAs
Cytoscape was employed for visualization of the regulatory networks between DEmiRNA and their target mRNAs. A single mRNA can be targeted by various miRNAs, whereas a single miRNA is also capable of targeting different mRNAs [57]. In these regulatory networks, it can be clearly seen that DEmiRNAs and their target mRNAs in AcCK1 vs AcT1 formed complex networks. MiR-598-y had as many as 45 target mRNAs, and miR-252-y bound to three targets including XM_017049297.1, XM_017049298.1 and XM_017058024.1 ( Figure 8A). In contrast, miR-6313-y, miR-3726-x, miR-9204-x and miR-6717-x could bind to only one target ( Figure 8A). In addition, as shown in Figure 8B, for DEmiRNAs in AcCK2 vs AcT2, miR-92-x, miR-3654-y, novel-m0003-3p and miR-6313-y linked to six, four, two and one target mRNAs, respectively.

Verification of A. c. cerana DEmiRNAs via RT-qPCR
To verify the high-throughput sequencing data in our study, the expression levels of nine randomly selected DEmiRNAs in AcCK1 vs AcT1 and AcCK2 vs AcT2 were examined by RT-qPCR. The expression trends of eight DEmiRNAs in the two comparison groups were roughly the same as the result of the sRNA-seq (Figure 9), indicative of the accuracy and reliability of our sequencing data.   , while the right two charts showed the differential expression (log2(Fold change)) of the same selected DEmiRNA from sRNA-seq data (in dark gray). RT-qPCR data are represented as the mean ± SD.

Discussion
Previous studies have mainly been focused on western honeybee miRNAs [58][59][60][61][62]. However, relevant studies on the eastern honeybee were extremely limited. Shi et al. previously performed Illumina sequencing of miRNAs in the royal jelly from A. mellifera and A. cerana, and the results showed that there were differences between them; they further found that the transcriptomes of A. mellifera adults were affected by the two kinds of royal jelly [63]. Over the years, increasing evidence has demonstrated that miRNAs play vital roles in host-pathogen interactions via regulation of the gene expression of hosts or pathogens [64][65][66]. Although advances have been made regarding the actions of miRNAs in the interactions between the western honeybee and N. ceranae [67,68], the miRNA profiles and DEmiRNAs in eastern honeybee-N. ceranae interactions remains completely unknown. In a previous study, Huang et al. deep-sequenced A. mellifera miRNAs daily across the sixday reproductive cycle of N. ceranae and identified 17 DEmiRNAs that target more than 400 genes [28]. Here, in order to identify miRNAs involved in the eastern honeybee responding to N. ceranae infection, we constructed for the first time 12 small RNA libraries from untreated and N. ceranaetreated A. c. cerana workers' midgut samples. In total, 529 known and 25 novel miRNAs were RT-qPCR validation of miR-3319-y and miR-1672-x; snRNA U6 was used as a positive control in Stem-loop RT-PCR and reference gene in RT-qPCR. The experiment was performed three times using three independent biological samples. The left two charts showed the result of RT-qPCR of selected DEmiRNA (in light gray), while the right two charts showed the differential expression (log 2 (Fold change)) of the same selected DEmiRNA from sRNA-seq data (in dark gray). RT-qPCR data are represented as the mean ± SD. RT-qPCR validation of miR-3319-y and miR-1672-x; snRNA U6 was used as a positive control in Stem-loop RT-PCR and reference gene in RT-qPCR. The experiment was performed three times using three independent biological samples. The left two charts showed the result of RT-qPCR of selected DEmiRNA (in light gray), while the right two charts showed the differential expression (log 2 (Fold change)) of the same selected DEmiRNA from sRNA-seq data (in dark gray). RT-qPCR data are represented as the mean ± SD.

Discussion
Previous studies have mainly been focused on western honeybee miRNAs [58][59][60][61][62]. However, relevant studies on the eastern honeybee were extremely limited. Shi et al. previously performed Illumina sequencing of miRNAs in the royal jelly from A. mellifera and A. cerana, and the results showed that there were differences between them; they further found that the transcriptomes of A. mellifera adults were affected by the two kinds of royal jelly [63]. Over the years, increasing evidence has demonstrated that miRNAs play vital roles in host-pathogen interactions via regulation of the gene expression of hosts or pathogens [64][65][66]. Although advances have been made regarding the actions of miRNAs in the interactions between the western honeybee and N. ceranae [67,68], the miRNA profiles and DEmiRNAs in eastern honeybee-N. ceranae interactions remains completely unknown. In a previous study, Huang et al. deep-sequenced A. mellifera miRNAs daily across the six-day reproductive  , while the right two charts showed the differential expression (log2(Fold change)) of the same selected DEmiRNA from sRNA-seq data (in dark gray). RT-qPCR data are represented as the mean ± SD.

Discussion
Previous studies have mainly been focused on western honeybee miRNAs [58][59][60][61][62]. However, relevant studies on the eastern honeybee were extremely limited. Shi et al. previously performed Illumina sequencing of miRNAs in the royal jelly from A. mellifera and A. cerana, and the results showed that there were differences between them; they further found that the transcriptomes of A. mellifera adults were affected by the two kinds of royal jelly [63]. Over the years, increasing evidence has demonstrated that miRNAs play vital roles in host-pathogen interactions via regulation of the gene expression of hosts or pathogens [64][65][66]. Although advances have been made regarding the actions of miRNAs in the interactions between the western honeybee and N. ceranae [67,68], the miRNA profiles and DEmiRNAs in eastern honeybee-N. ceranae interactions remains completely unknown. In a previous study, Huang et al. deep-sequenced A. mellifera miRNAs daily across the sixday reproductive cycle of N. ceranae and identified 17 DEmiRNAs that target more than 400 genes [28]. Here, in order to identify miRNAs involved in the eastern honeybee responding to N. ceranae infection, we constructed for the first time 12 small RNA libraries from untreated and N. ceranaetreated A. c. cerana workers' midgut samples. In total, 529 known and 25 novel miRNAs were identified from sRNA-seq data using bioinformatics, which enriches the miRNA reservoir of the eastern honeybee and provides candidate miRNAs for functional studies in the future. In addition, eight predicted miRNAs were confirmed to be expressed in the midguts of A. c. cerana workers via The experiment was performed three times using three independent biological samples. The left two charts showed the result of RT-qPCR of selected DEmiRNA (in light gray), while the right two charts showed the differential expression (log 2 (Fold change)) of the same selected DEmiRNA from sRNA-seq data (in dark gray). RT-qPCR data are represented as the mean ± SD. The experiment was performed three times using three independent biological samples. The left two charts showed the result of RT-qPCR of selected DEmiRNA (in light gray), while the right two charts showed the differential expression (log 2 (Fold change)) of the same selected DEmiRNA from sRNA-seq data (in dark gray). RT-qPCR data are represented as the mean ± SD.

Discussion
Previous studies have mainly been focused on western honeybee miRNAs [58][59][60][61][62]. However, relevant studies on the eastern honeybee were extremely limited. Shi et al. previously performed Illumina sequencing of miRNAs in the royal jelly from A. mellifera and A. cerana, and the results showed that there were differences between them; they further found that the transcriptomes of A. mellifera adults were affected by the two kinds of royal jelly [63]. Over the years, increasing evidence has demonstrated that miRNAs play vital roles in host-pathogen interactions via regulation of the gene expression of hosts or pathogens [64][65][66]. Although advances have been made regarding the actions of miRNAs in the interactions between the western honeybee and N. ceranae [67,68], the miRNA profiles and DEmiRNAs in eastern honeybee-N. ceranae interactions remains completely unknown. In a previous study, Huang et al. deep-sequenced A. mellifera miRNAs daily across the six-day reproductive RT-qPCR validation of miR-3319-y and miR-1672-x; snRNA U6 was used as a positive control in Stem-loop RT-PCR and reference gene in RT-qPCR. The experiment was performed three times using three independent biological samples. The left two charts showed the result of RT-qPCR of selected DEmiRNA (in light gray), while the right two charts showed the differential expression (log 2 (Fold change)) of the same selected DEmiRNA from sRNA-seq data (in dark gray). RT-qPCR data are represented as the mean ± SD.

Discussion
Previous studies have mainly been focused on western honeybee miRNAs [58][59][60][61][62]. However, relevant studies on the eastern honeybee were extremely limited. Shi et al. previously performed Illumina sequencing of miRNAs in the royal jelly from A. mellifera and A. cerana, and the results showed that there were differences between them; they further found that the transcriptomes of A. mellifera adults were affected by the two kinds of royal jelly [63]. Over the years, increasing evidence has demonstrated that miRNAs play vital roles in host-pathogen interactions via regulation of the gene expression of hosts or pathogens [64][65][66]. Although advances have been made regarding the actions of miRNAs in the interactions between the western honeybee and N. ceranae [67,68], the miRNA profiles and DEmiRNAs in eastern honeybee-N. ceranae interactions remains completely unknown. In a previous study, Huang et al. deep-sequenced A. mellifera miRNAs daily across the six-day reproductive cycle of N. ceranae and identified 17 DEmiRNAs that target more than 400 genes [28]. Here, in order to identify miRNAs involved in the eastern honeybee responding to N. ceranae infection, we constructed for the first time 12 small RNA libraries from untreated and N. ceranae-treated A. c. cerana workers' midgut samples. In total, 529 known and 25 novel miRNAs were identified from sRNA-seq data using bioinformatics, which enriches the miRNA reservoir of the eastern honeybee and provides candidate miRNAs for functional studies in the future. In addition, eight predicted miRNAs were confirmed to be expressed in the midguts of A. c. cerana workers via Stem-loop RT-PCR ( Figure 4). Additionally, we also validated the expression of another eight DEmiRNAs (Figure 9A-C,G-I,M,N). Deep sequencing has become the mainstream method to explore novel miRNAs at present [69]. Notably, all of the 25 novel miRNAs predicted in this work were weakly expressed, in accordance with the results for some lower invertebrates such as Drosophila and Cyprinus carpio [57,70]. In our study, the most abundant miRNAs in AcCK1 and AcCK2 were bantam-y, miR-184-y, miR-1-y, miR-276-y, miR-750-y, miR-3477-x, miR-283-x, miR-12-z, miR-31-x and miR-9-z with TPM ≥ 31,625.97 each (Table S4). Surprisingly, these miRNAs were also highly expressed in AcT1 and AcT2 (Table S4). As a versatile player in many biological processes in Drosophila, bantam regulates fly development [71,72], cell proliferation and apoptosis in many cell types [73,74], and the production of the molting hormone ecdysone [75]. Through systematic manipulation of Dpp signaling in the Drosophila wing, Zhang et al. discovered that Dpp promoted proliferation in the lateral wing disc and repressed proliferation in the medial wing disc via omb, which controlled the regional proliferation rate by oppositely regulating transcription of bantam in the medial versus lateral wing disc [71]. Wu et al. previously investigated the role of the bantam miRNA in the regulation of neuroblast homeostasis in the Drosophila brain and found that bantam is a direct transcriptional target of the Notch signaling pathway; in addition, bantam feedback regulates the pathway by negatively regulating its target mRNA numb [72]. In the current work, bantam-y was the miRNA with the highest expression level in both AcCK1 and AcCK2 (Table S4), implying its key role in the regulation of the development of the A. c. cerana workers' midguts. MiR-2, an invertebrate-specific miRNA family that has been predicted in fruit flies, is highly expressed in the heads of Drosophila and Bombyx mori [76,77] and the neurons in Caenorhabditis [78]. It has also been shown to participate in the regulation of apoptosis of neuroblasts during the normal development of the nervous system in Drosophila [79]. Recent evidence shows that miR-2 regulates functions essential for the normal wing morphogenesis of B. mori via targeting of awd and fng [80]. In the embryonic development of Drosophila, miRs-2/13, including miR-2a, miR-2b, miR-13a, and miR-13b, play an important regulatory part in the embryonic development of Drosophila by targeting nine genes, including Sos and Myd88 [81]. In the present study, miR-2-y and miR-13-y were detected to be highly expressed in the normal midguts of A. c. cerana workers (Table S4), which suggests these two miRNAs may play a fundamental role in regulating the development and morphogenesis of the midgut and apoptosis of the midgut cells. However, the functions of other miRNAs that were abundant in A. c. cerana midgut samples, such as miR-750-y and miR-3477-x, remain largely unknown.
DEmiRNAs in the midguts of A. c. cerana workers challenged by N. ceranae are particularly worthy of attention since they are likely to be direct regulators of host-pathogen interactions. In this research, 14 and 12 DEmiRNAs were identified in AcCK1 vs AcT1 and AcCK2 vs AcT2, respectively. Among them, miR-60-y, miR-8462-x, miR-2965-y, miR-676-y, and miR-6313-y were shared by the two comparison groups, implying their primary roles in host responses to N. ceranae infections. In addition, nine DEmiRNAs, including miR-1-x, miR-980-y, miR-965-x, miR-598-y, miR-6717-x, miR-4635-y, miR-9204-x, miR-3726-x, and miR-252-y, were specific in AcCK1 vs AcT1, and all of them were known miRNAs ( Figure 5B, see also Tables S6 and S7). Additionally, seven DEmiRNAs, including two novel miRNAs (novel-m0003-3p and novel-m0019-5p), were specific in AcCK2 vs AcT2 ( Figure 5B, see also Tables S6 and S7). We believe these unique DEmiRNAs play specific parts during different stages of host responses under the stress of N. ceranae.
To gain further insight into the potential genes associated with N. ceranae infection, we conducted GO analysis for DEmiRNA target mRNAs to identify biologically important terms. The mostly enriched terms for the targets in AcCK1 vs AcT1 and AcCK2 vs AcT2 were binding, cellular process, biological regulation, single-organism process, and catalytic activity ( Figure 6, see also Tables S8 and S9), indicating the regulation of host metabolism and cellular activity via DEmiRNAs in response to N. ceranae. Huang and colleagues found that the target mRNAs of DEmiRNAs in the N. ceranae-infected midgut of the A. mellifera workers were mostly involved in binding, signaling, nucleus, transmembrane transport, and DNA binding, which differ from our findings in the present study [28]. This difference reflects the different responses of the two honeybee species with different resistances to the same fungal pathogen, N. ceranae. Considering the fact that A. cerana is more resistant to N. ceranae compared to A. mellifera, we inferred that DEmiRNAs may in part play a special role in the N. ceranae-resistance difference. In addition, we found 88 and 89 DEmiRNA-targeted mRNAs in AcCK1 vs AcT1 and AcCK2 vs AcT2 were respectively enriched in response to stimulus, and there was one target mRNA for DEmiRNAs in AcCK1 vs AcT1 enriched in immune system process, which indicates the regulatory roles of the corresponding DEmiRNAs in host immunity defenses against N. ceranae. Furthermore, KEGG analysis for the target mRNAs was performed to investigate the pathways influenced by N. ceranae. The result showed that out of the 104 pathways enriched by target mRNAs of DEmiRNAs in AcCK1 vs AcT1, 62 were related to material and energy metabolism, such as carbohydrate metabolism including the pentose phosphate pathway and citrate cycle, lipid metabolism including arachidonic acid metabolism and fatty acid biosynthesis, nucleotide metabolism including purine metabolism and pyrimidine metabolism, energy metabolism including sulfur metabolism and oxidative phosphorylation (Table S10). A total of 92 pathways were enriched by DEmiRNA target mRNAs in AcCK2 vs AcT2; among them, 52 pathways were associated with material and energy metabolism (Table S11). These results demonstrated that host DEmiRNAs regulate metabolism-related genes in response to N. ceranae invasion. Meanwhile, immunity-related pathways were further analyzed, and the results showed that endocytosis, phagosome, lysosome, ubiquitin mediated proteolysis, Jak-STAT signaling pathway, and MAPK signaling pathways were enriched by the target mRNAs in AcCK1 vs AcT1 (Table S10). Surprisingly, these immune pathways were also enriched by DEmiRNA target mRNAs in AcCK2 vs AcT2 (Table S11). These results demonstrated that DEmiRNAs and their target mRNAs were involved in cellular and humoral immune responses of the host to N. ceranae invasion. Together, these findings indicate that host material and energy metabolism, cellular activity, and immunity defenses were significantly impacted by N. ceranae infection, and DEmiRNAs played a comprehensive role in host-pathogen interactions.
MiR-1 is a muscle-specific miRNA that plays important roles in regulating heart development and muscle differentiation [82,83]. It has also been known to be a tumor suppressor gene that not only inhibits cancer cell proliferation, metastasis, and invasion, but also enhances apoptosis by regulating oncogenic targets in many cancer cells [84,85]. In insects, miR-1 was found to be involved in immune processes in Drosophila melanogaster [86] and Aedes aegypti [87]. Besides, Huang et al. suggested a significant down-regulation of ame-miR-1 in Apis mellifera workers at 6 dpi with N. ceranae [25]. In this work, miR-1-x was dramatically down-regulated in AcCK1 vs AcT1, suggesting that the host and fungal pathogen can interact with each other via miR-1-x by regulating host proliferation and apoptosis and immune response. Considering the conservation between miR-1-x and ame-miR-1, we speculated that miR-1 family participates in honeybee response to N. ceranae invasion as a general regulator. Wu et al. carried out Illumina sequencing of dengue virus-2 (DENV-2)-infected and uninfected Aedes albopictus and observed a series of DEmiRNAs including miR-252 [88]. By using RT-qPCR, Yan et al. detected a high expression level of miR-252 in DENV-2-infected C6/36 cells and mosquitoes compared with uninfected cells and mosquitos and further revealed that DENV-2 envelope protein expression can be affected by the overexpression of miR-252 with a mimic or down-regulation using an inhibitor [89]. Recently, Lim and colleagues discovered that miR-252-5p can control the cell cycle by directly repressing Abelson interacting protein (Abi) in Drosophila S2 cells by using cross-linking immunoprecipitation and deep sequencing of endogenous Argonaute 1 (Ago1) protein [90]. In our study, miR-252-y had a weak up-regulation in the midgut 7 dpi with N. ceranae, suggesting that miR-252-y may participate in N. ceranae defense in A. c. cerana. As a member of the miR-17-92 cluster, miR-92a inhibits apoptosis and promotes the proliferation of other cell types, including endothelial cells [91,92]. Kurze et al. first showed the importance of apoptosis for N. ceranae infections [93]. In addition, apoptosis is of great importance for pathogen defense in multicellular organisms including B. mori [94]. Previous studies have demonstrated that Spodoptera frugiperda larvae and B. mori larvae battle baculovirus infection by selective apoptosis of infected cells from the midgut epithelium [95,96]. N. ceranae is able to enhance its development during the infection process by preventing apoptosis in epithelial cells of infected A. mellifera [97]. A previous study showed that miR-194 expression was markedly decreased in A549 alveolar epithelial cells following infection with influenza A virus (IAV), and miR-194 could suppress fibroblast growth factor 2 (FGF2) expression at the mRNA and protein levels [98]. Xie and colleagues suggested that miR-194 increases cell apoptosis by inhibiting the NF-κB pathway in WI38 cells, and it may be used as a potential targeted therapy for the treatment of infantile pneumonia [99]. In this present study, miR-92-x had a sharp down-regulation, while miR-194-y was significantly up-regulated in the A. c. cerana worker's midgut at 10 dpi with N. ceranae. Therefore, we inferred that A. c. cerana could fight N. ceranae by promoting apoptosis of host midgut cells via regulation of the expression levels of miR-92-x and miR-194-y. Further studies are needed to confirm the mechanism underlying the roles of above-mentioned DEmiRNAs in host-pathogen interactions. Recently, Evans and Huang discovered both N. ceranae and A. mellifera gene expression were altered by knockdown of Dicer gene of N. ceranae, and further found parasite miRNAs may regulate host metabolism and immune response via cross-talk with A. mellifera mRNAs [68]. It's known that complex interactions exist between honeybee and microsporidian [28,67,68]. Thus, whether A. c. cerana miRNAs identified in this study can regulate N. ceranae mRNAs, and whether N. ceranae miRNAs could regulate A. c. cerana mRNAs, is a very interesting research direction.
A single miRNA can regulate several target mRNAs at the same time to inhibit their expression and vice versa [100]. MiR-598 has been found to play inhibitory role in osteosarcoma progression in vivo and in vitro by modulating osteoblastic differentiation in the microenvironment and targeting PDGFB and MET. Based on expressed sequence tags (ESTs) resources from the LNCaP cells, Saravanan et al. identified an hsa-miR-3654 with a higher expression level in LNCaP cells than the normal and androgen insensitive prostate cancer cell lines (PNT1A, PC-3) [101]. To our knowledge, there is no documentation of miR-598 and miR-3654 in insects until now. In the current work, miR-598 lies in the center of the regulation networks of DEmiRNAs in AcCK1 vs AcT1, linking to as many as 45 target mRNAs such as XM_017048388.1, XM_017049101.1, and XM_017051459.1 ( Figure 8A), while miR-3654-y was a key miRNA in the regulation networks of DEmiRNAs in AcCK2 vs AcT2, binding to four targets including XM_017057343.1, XM_017057344.1, XM_017057345.1, and XM_017057346.1 ( Figure 8B). These results indicated that miR-598 and miR-3654-y are likely to be key regulators of the interactions between A. c. cerana and N. ceranae, and these two miRNAs would be our candidates for future studies of host immune response and resistance to N. ceranae.
In conclusion, using high-throughput sequencing technology and bioinformatics, we systematically investigated the highly expressed miRNAs in the normal A. c. cerana midgut, and DEmiRNAs as well as their target mRNAs in the N. ceranae-infected midgut. This study provides the first global view of the miRNA profiles and DEmiRNAs in the midgut of the A. c. cerana worker invaded by N. ceranae. Findings in the present study demonstrated that N. ceranae invasion causes alterations in the expression of miRNAs that regulate metabolism, cellular activity, and immune response in the host. Taken together, our results not only provide novel insights into understanding A. c. cerana-N. ceranae interactions but also lay a foundation for deciphering the molecular mechanisms underlying the resistance of the eastern honeybee to N. ceranae.

Conclusions
In a nutshell, this is the first study on eastern honeybee miRNA response to N. ceranae infection from a genome-wide perspective, which demonstrates the expression pattern of miRNAs in A. c. cerana workers was altered by N. ceranae invasion. We identified 579 known miRNAs and 25 novel miRNAs in the midguts of A. c. cerana workers, and showed that 14 and 12 DEmiRNAs were N. ceranae-responsive in midguts at 7 dpi and 10 dpi, respectively. DEmiRNAs may play an important role in the host stress responses including immune responses via negatively regulate target mRNAs. Our results provide a rich genetic resource and potential candidates for further investigation of the regulatory roles of miRNAs in host responses to N. ceranae. Furthermore, the study may help elucidating the complex molecular mechanisms underlying the eastern honeybee-microsporidian interactions.