Analysis of miRNAs in the Heads of Different Castes of the Bumblebee Bombus lantschouensis (Hymenoptera: Apidae)

Abstract

Bumblebees are important insect pollinators for many wildflowers and crops. MicroRNAs (miRNAs) are endogenous non-coding small RNAs that regulate different biological functions in insects. In this study, the miRNAs in the heads of the three castes of the bumblebee Bombus lantschouensis were identified and characterized by small RNA deep sequencing. The significant differences in the expression of miRNAs and their target genes were analyzed. The results showed that the length of the small RNA reads from males, queens, and workers was distributed between 18 and 30 nt, with a peak at 22 nt. A total of 364 known and 89 novel miRNAs were identified from the heads of the three castes. The eight miRNAs with the highest expressed levels in males, queens, and workers were identical, although the order of these miRNAs based on expression differed. The male vs. queen, male vs. worker, and worker vs. queen comparisons identified nine, fourteen, and four miRNAs with significant differences in expression, respectively. The different castes were clustered based on the differentially expressed miRNAs (DE miRNAs), and the expression levels of the DE miRNAs obtained by RT-qPCR were consistent with the read counts obtained through Solexa sequencing. The putative target genes of these DE miRNAs were enriched in 29 Gene Ontology (GO) terms, and catalytic activity was the most enriched GO term, as demonstrated by its association with 2837 target genes in the male vs. queen comparison, 3535 target genes in the male vs. worker comparison, and 2185 target genes in the worker vs. queen comparison. This study highlights the characteristics of the miRNAs in the three B. lantschouensis castes and will aid further studies on the functions of miRNAs in bumblebees.

Keywords: Bombus lantschouensis; Solexa sequencing; microRNA; pollinators.

Figure 1

Distribution of the lengths of the small RNA reads in the three B. lantschouensis castes. The black vertical lines represent error bars. The size distribution of the 18–30 nt clean reads obtained from all three groups was assessed. Most sequences in the three libraries were 21–23 nt in length, and the most abundant size was 22 nt.

Figure 2

Distribution of known microRNAs (miRNAs) in the three B. lantschouensis castes. The Venn diagram displays the distribution of 364 unique known miRNAs among the male, queen, and worker libraries.

Figure 3

Heat map diagram showing the differentially expressed miRNAs (DE miRNAs) in the three B. lantschouensis castes. Nineteen differentially expressed miRNAs with |fold change| ≥1 and p < 0.05 were screened using the DESeq R package. The rows represent the different miRNAs, and the columns represent males (M1, M2, M3), workers (W1, W2, W3), and queens (Q1, Q2, Q3). The expression data for each miRNA were calculated from three biological replicates.

Figure 4

Gene Ontology (GO) enrichment analysis of potential target genes of the DE miRNAs of B. lantschouensis obtained in the male vs. queen, male vs. worker, and worker vs. queen comparisons. The x-axis shows the GO category, and the y-axis indicates the number of genes.

Figure 5

Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways of potential target genes of DE miRNAs of B. lantschouensis obtained in the male vs. queen, male vs. worker, and worker vs. queen comparisons. The x-axis shows the pathway terms, and the y-axis indicates the number of genes.

Figure 6

RT-qPCR validation of the DE miRNAs in the three B. lantschouensis castes identified by Solexa small RNA sequencing. RT-qPCR was performed using nine B. lantschouensis heads, which included three biological replicates of the three castes. The log2 (fold-change) relative expression values obtained from the male vs. queen, male vs. worker, and worker vs. queen comparisons are shown. The black vertical lines represent error bars. In the RT-qPCR assays, the expression of the miRNAs was normalized to the U6 snRNA levels. The data were statistically analyzed using the comparative quantity (2^−ΔΔCT) method.