Chao Tianle
Yang Liuxu
Lou Delong2
Wang Guizhi- 1Shandong Provincial Key Laboratory of Animal Biotechnology and Disease Control and Prevention, College of Animal Science and Veterinary Medicine, Shandong Agricultural University, Tai’an, China
- 2Comprehensive Testing and Inspection Center, Shandong Provincial Animal Husbandry and Veterinary Bureau, Jinan, China
- 3Animal Husbandry Development Center of Linqu County, Weifang, China
Fluvalinate is a widely used and relatively safe acaricide for honeybees, but it still has a negative impact on honeybee colonies. Such negative effects may be related to fluvalinate-induced brain nerve tissue damage, but the detailed molecular regulatory mechanism of this phenomenon is still poorly understood. In this study, we analyzed the miRNA expression profile changes in the brain tissue of Apis mellifera ligustica by miRNA sequencing after fluvalinate treatment. A total of 1,350 miRNAs were expressed in Apis mellifera ligustica brain tissue, of which only 180 were previously known miRNAs in honeybees. Among all known and novel miRNAs, 15 were differentially expressed between at least two of the four time periods before and after fluvalinate administration. Further analysis revealed five significantly enriched KEGG pathways of the differentially expressed miRNA (DEM) potential target genes, namely, “Hippo signaling pathway-fly,” “Phototransduction-fly,” “Apoptosis-fly,” “Wnt signaling pathway,” and “Dorso-ventral axis formation,” which indicates that differentially expressed miRNA function may be related to cell apoptosis and memory impairment in the fluvalinate-treated Apis mellifera ligustica brain. Ame-miR-3477-5p, ame-miR-375-3p, and miR-281-x were identified as key miRNAs. Overall, our research provides new insights into the roles of miRNAs in brain tissue during the process of fluvalinate-induced Apis mellifera ligustica poisoning.
Introduction
Among all of the commonly used anti-mite insecticides for honeybees, fluvalinate is known to have low toxicity (Johnson et al., 2006; Johnson et al., 2009). However, there is growing evidence that fluvalinate is not entirely benign for Apis mellifera. It has been reported that honeybee workers that received field-level fluvalinate treatment showed a higher mortality rate than those that did not receive treatment (Medici et al., 2012; Berry et al., 2013). Apis mellifera individuals exposed to fluvalinate have increased mortality, lower body weight, reduced reproductive viability, and decreased learning and memory ability (Rinderer et al., 1999; Burley et al., 2008; Berry et al., 2013; Frost et al., 2013). Furthermore, as the primary molecular target for pyrethroid insecticides, the Apis mellifera voltage-dependent sodium channel has been reported to have a higher affinity than a variant of the Varroa destructor voltage-dependent sodium channel for fluvalinate (Soderlund, 2010; Gosselin‐Badaroudine and Chahine, 2017). Although there is no direct evidence to date, growing evidence of memory loss, learning ability decrease, and toxicity for voltage-gated sodium channels indicate that the lethal and sublethal effects of fluvalinate on honeybees may be partly due to nervous system damage.
MicroRNA is a type of nonprotein-coding endogenous small RNA molecule with a sequence length of 20–24 nt. By inhibiting or degrading mRNA, it is involved in many biological regulatory processes at the posttranscriptional level in plants, animals, and microorganisms and it participates in many pathways in vivo, such as regulating cell proliferation, differentiation, apoptosis, stress response, and metabolism (Griffiths-Jones 2004; Calin and Croce, 2006). Some effects of miRNAs on the nervous system of A. mellifera ligustica have been reported.
This study used transcriptome sequencing technology to reveal changes in the miRNA expression profiles of Apis mellifera ligustica brain tissue under the influence of fluvalinate, identify miRNAs that were differentially expressed before and during fluvalinate treatment, and study the potential function of these miRNAs at different time points after fluvalinate treatment. The miRNAs involved in the fluvalinate response were revealed, and functional enrichment analysis of the miRNA target genes was used to perform a preliminary assessment of the pathways and mechanisms by which miRNAs are involved in regulating this process. Our study provides new insights into miRNA function in fluvalinate-induced honeybee brain tissue molecular regulation changes.
Materials and Methods
Sample Collection
We took six colonies of Apis mellifera ligustica with similar colony potential and randomly divided them into two sets (six colonies in each set). Each colony lives in a standard sixteen-frame hive with one naturally mated queen inside, and the colony size was similar between the two sets (approximately 16,000–18,000 workers in each colony). One was set as the control and the other as the test group. The test groups were exposed to fluvalinate (fluvalinate strip, Sichuan Pengshan Wangshi Animal Health Co., Ltd.) at a normal dosage for administration to honeybees (according to the instructions of the drug, a single strip was placed in the beehive every 7 days, which contained 40 mg of fluvalinate). For all of the colonies, we sampled foragers that returned from 08:00 to 12:00. Sample honeybees were taken 1 day before and 10, 20, and 30 days after fluvalinate treatment began, and 20 adult worker honeybees were collected from each bee colony for each sampling time. All collected honeybees were immediately used for dissection under a microscope, and 20 bee brains collected from the same colony were pooled for subsequent RNA extraction. The royal jelly glands and pigments were removed from the brain, and the brain tissue was immediately stored in an RNA preservation solution (Beijing Tiangen Biochemical Technology Company). The sampling site was Fuxin Honeybee Farm in Tai’an City, Shandong Province, China, which provided a uniform breeding environment for the experimental animals.
RNA Extraction, Library Construction, and Sequencing
We used TRIzol to extract total RNA from each sample, and the RNA molecules in a size range of 18–30 nt were then enriched by polyacrylamide gel electrophoresis (PAGE). The RNA quality was evaluated with an Agilent 2100 Bioanalyzer. The total RNA samples satisfied the required quality of RNA integrity number (RIN) > 7.0. Then, 1 μg RNA per sample was used for library construction according to the Illumina TruSeq sample RNA Sample Prep Kit (Cat#RS-200-0012). Sequencing was performed by Guangzhou Gene Denovo Biotechnology Co., Ltd., using an Illumina HiSeq TM 2500. In particular, libraries from untreated brain tissue samples were named CK1, CK2, and CK3; libraries from samples collected 10 days after fluvalinate treatment were named D10-1, D10-2, and D10-3; libraries from samples collected 20 days after fluvalinate treatment were named D20-1, D20-2, and D20-3; and libraries from samples collected 30 days after fluvalinate treatment were named D30-1, D30-2, and D30-3.
Sequencing Data Quality Control
The low-quality reads were filtered out of the data to obtain high-quality reads. Then, reads without 3′ linkers and 5′ linkers, reads without inserts or with inserts less than 18 nt in length, and polyA reads were removed. Next, small RNA tags with a read frequency lower than two were also removed. Finally, we obtained clean small RNA tag sequences.
Sequencing Data Alignment
The software Blastall 2.2.25 (blastn) was used to align all the clean small RNA tags with the small RNAs in GenBank and compare them with the small RNAs in the Rfam database (release 11.0) (Griffiths-Jones et al., 2003) to identify and remove rRNA, scRNA, snoRNA, snRNA, and tRNA. Bowtie software (version 1.1.2) was used to compare all clean reads with the A. mellifera ligustica reference genome (GCF-003254395.2), and all the reads mapped to exons, introns, and repeat sequences were removed. According to the genome alignment results, the small RNA tag sequence set that characterized each group was determined.
We also used Bowtie (version 1.1.2) alignment to identify the known miRNAs of A. mellifera ligustica. The same method was used to compare the miRNA sequences of all the species in the miRbase database (V22) (Griffiths-Jones et al., 2006) to identify homologous miRNAs. Then, the clean reads that had not yet been annotated were compared with the reference genome, and novel miRNAs were predicted based on the genome position and hairpin structure predicted by MIREAP V0.2 software (Wang et al., 2015).
Identification of Differentially Expressed MicroRNA
The miRNA expression abundance was evaluated by the transcripts per million reads (TPM). The edgeR package was used to perform a differential expression analysis of the miRNAs, and miRNAs with |log2(FC)| > 1 and FDR < 0.05 were identified as DEMs.
MicroRNA Target Gene Prediction
RNAhybrid (version 2.1.2) (Krüger and Rehmsmeier, 2006) + svm-light (version 6.01) (Joachims, 1999), Miranda (version 3.3a), and TargetScan (version 7.0) (Turner, 1985; Lewis et al., 2005) were applied for miRNA target gene prediction. The intersection of the software results was taken as the potential miRNA target gene set.
Differentially Expressed MicroRNA Expression Profile Clustering
The z-score of each DEM was calculated using a tpm value. Using the MultiExperiment Viewer (MeV) (Saeed et al., 2003), hierarchical clustering was performed with the z-score.
Quantitative Real-Time PCR (qRT-PCR) Validation of Differentially Expressed MicroRNA
To further verify the accuracy of our sequencing data, 10 DEMs were selected for the qRT-PCR test. A LightCycler 480 real-time system (Roche, United States) was applied for the experiment. The detailed information of primers for qRT-PCR testing was listed in Supplementary Table S1, and U5 snRNA and 18S rRNA (Cristino et al., 2014; Mello et al., 2014) were selected as the reference gene. qRT-PCR was performed according to the instructions of the Mir-X™ miRNA qRT-PCR TB Green™ Kit (Clontech). All reactions were performed using three biological replicates, and each biological replicate had three technical replicates. The qRT-PCR testing results were calculated using 2−ΔΔCT ± SEM, in which the geometric mean of the threshold cycles (Ct) values of the two reference genes was used for the calculation, and the significance of the difference in miRNA expression was tested using Tukey’s honestly significant difference (HSD).
Functional Enrichment Analysis
In order to generate a comprehensive description of the biological characteristics of the DEMs, potential target genes of DEMs from three different expression profile clusters were respectively enriched into the Gene Ontology (GO) database (Ashburner et al., 2000) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al., 2008) with FDR ≤ 0.05 using the R package clusterProfiler (Yu et al., 2012).
Network Analysis
We used Cytoscape 3.7.0 to construct the miRNA-KEGG pathway network diagram. Based on the results of KEGG enrichment, we used the significantly enriched KEGG pathways and 15 DEMs as nodes to make a network diagram. DEM node sizes were determined by the expression values.
Results
Assessment of Sequencing Quality
A total of 184,051,902 clean reads were obtained by sequencing. The sample with the most clean reads was D30-3, with 18,694,375 clean reads, and the lowest was CK-1, with 9,844,859. The number of clean tag sequences in each sample ranged from 18,109,705 to 9,626,768, and the filtering rate ranged from 94% to 98%. The number of reads in each sample library was within a reasonable range, and there were no abnormal samples (Supplementary Table S2).
Alignment Analysis
Comparisons were used to determine the miRNA contents present in the samples and their base distributions and identify known and novel miRNAs. The proportion of previously known miRNAs, that is, the honeybee miRNAs included in miRBase, was above 80%. The proportions of the miRNAs known from other species but not annotated in the honeybee miRBase data were all approximately 2%, while some newly discovered miRNAs accounted for only 0.1%. Overall, 1,350 miRNAs were detected, of which only 180 were previously known miRNAs in honeybees, while 924 were known by comparison with other species’ known miRNAs, and 246 were newly discovered novel miRNAs (Supplementary Table S3).
Analysis of Differential MicroRNA Expression
Using |log2(FC)| > 1 and FDR < 0.05 as the thresholds, 15 miRNAs were identified as differentially expressed among at least two different time points before or after fluvalinate exposure (Supplementary Table S4). A heatmap of the expression levels of these 15 DEMs in each group was drawn, and expression profile cluster analysis was performed (Figure 1). All 15 DEMs were divided into three clusters, and DEMs included in each cluster had similar expression patterns. Among them, 10 DEMs were selected for qRT-PCR testing to validate the accuracy of the sequencing data. As shown in Figure 2, the qRT-PCR results of 9 miRNAs were consistent with the RNA-seq results. The miR-750-y qRT-PCR results did not show a trend of differences among the different time periods.
FIGURE 1. Expression pattern cluster analysis of differential expressed miRNAs. Each column represents a sample, and each row represents a miRNA. The expression level of miRNA in different samples is indicated by different colors. The redder the color, the higher the expression level, and the bluer the color, the lower the expression level.
FIGURE 2. RNA-Seq and qRT-PCR relative expression level of differential expressed miRNAs. (A) MiRNA expression TPM value of RNA-Seq, (B) MiRNA relative expression level of qRT-PCR. Significant differences are indicated by different symbols with FDR < 0.05 in RNA-Seq results and p < 0.05 in qRT-PCR results.
Two DEMs were contained in Cluster one, and their expression levels showed a continuous upward trend with increasing treatment time. In Cluster two, six DEMs showed the highest expression in the D10 stage. In Cluster three, the expression of the seven DEMs gradually decreased after treatment with fluvalinate. A reasonable guess is that miRNAs with the same expression pattern have similar functions. Therefore, to further explore the regulatory roles of these DEMs, we performed GO and KEGG enrichment analyses on the target genes of the three clusters of DEMs.
Gene Ontology Enrichment Analysis
To clarify the main functions and regulatory mechanisms associated with the DEMs, we performed GO analysis on their potential target genes (Supplementary Figure S1).
For Cluster one target genes, 40 biological process terms and six molecular function terms showed significant enrichment, while no significant result was obtained from the cellular component terms. Among all Cluster one enriched GO terms, the top 12 most significant terms were all related to RNA splicing and processing.
For Cluster two, 61 cellular component terms, 54 molecular function terms, and 71 biological process terms were significantly enriched. It is worth noting that, among the three types of terms, the terms with the highest significant enrichment level were all related to ion transmembrane transport and protein modification.
For Cluster three, 10 cellular component terms, 43 molecular function terms, and 44 biological process terms achieved significant enrichment levels. The top significantly enriched biological process terms were “cellular protein modification process, protein modification process, and macromolecule modification.”
Furthermore, we also noticed that 27 GO terms were significantly enriched in all three clusters. Among them, we found three biological processes related to advanced brain functions: “behavior,” “learning or memory,” and “cognition.”
Kyoto Encyclopedia of Genes and Genomes Enrichment Analysis
We performed KEGG pathway enrichment analysis on the target genes of the DEMs. The pathways that obtained significant enrichment are shown in Supplementary Table S5.
The results showed that, in Cluster one, 200 target genes were enriched in 53 KEGG pathways, and only five pathways were significantly enriched. Among them, the pathway with the most significant enrichment level and the largest number of target genes was “Dorsoventral axis formation.”
For Cluster two, 89 KEGG pathways were enriched, and 15 pathways achieved a significant level. Among them, the pathway with the most significant enrichment and the largest number of target genes was the “Wnt signaling pathway.”
In Cluster three, 749 target genes were enriched in 105 KEGG pathways, 11 of which were significantly enriched. The pathway with the most significant enrichment was “Phototransduction-fly,” and the pathway containing the largest number of target genes was “Wnt signaling pathway.”
To compare the KEGG pathway enrichment results of the three clusters of DEMs, we selected and combined the top five enriched pathways in each cluster (Figure 3). The results showed large differences among the three clusters. The KEGG pathways “Dorso-ventral axis formation,” “Wnt signaling pathway,” and “Hedgehog signaling pathway-fly” were significantly enriched in all three DEM clusters, which may indicate that these pathways are regulated by DEMs at multiple periods after fluvalinate treatment. The pathways “Phototransduction-fly” and “Hippo signaling pathway-fly” were significantly enriched in Clusters two and three, which may indicate that these pathways are regulated by DEMs immediately after fluvalinate treatment. The pathways “spliceosome,” “circadian rhythm-fly,” and “phagosome” were specifically significantly enriched in Clusters one, two, and three, respectively, which may indicate that these pathways are related to the unique functions of DEMs with different expression patterns.
FIGURE 3. KEGG enrichment for miRNA target genes of each cluster. The horizontal axis of the bubble graph represents the number of target genes in each KEGG pathway, the vertical axis represents different channel names, the different colors of bubbles represent the FDR values, and the size of the bubbles represents the gene ratio value.
To further confirm the pathways that each miRNA may regulate, KEGG enrichment analysis was also performed for all 15 DEM target genes, and the top five enriched pathways were selected for comparison (Figure 4). In Cluster one, target genes were enriched in nine KEGG pathways, and only “Dorso-ventral axis formation” achieved the top five significance levels in both DEMs. In Cluster two, target genes were enriched in 20 pathways, while the pathways “Wnt signaling pathway” and “Phototransduction-fly” contained the most target genes and DEMs. In Cluster three, target genes were enriched in 23 pathways, among which target genes of four DEMs were enriched in “Hippo signaling pathway-fly,” and target genes of three DEMs were enriched in “Phototransduction-fly” and “Starch and sucrose metabolism.” Among all of the top enriched pathways, the “Hippo signaling pathway-fly” achieved the top five significant levels in six DEMs with the largest DEM number, followed by the “Dorso-ventral axis formation,” “Wnt signaling pathway,” and “Phototransduction-fly” pathways; all achieved the top five significant levels in five DEMs. According to the above enrichment results, “Hippo signaling pathway-fly,” “Dorso-ventral axis formation,” “Wnt signaling pathway,” and “Phototransduction-fly” were considered the key KEGG pathways in the miRNA coordinated posttranscriptional regulation after fluvalinate treatment.
FIGURE 4. KEGG enrichment for target genes of each miRNA.
Network Analysis
Based on the KEGG enrichment results, we imported the DEMs and all KEGG pathways with significant enrichment as nodes into Cytoscape 3.7.0 to draw a network map, and the size of the miRNA nodes was determined by their expression level. In total, 15 miRNAs and 39 significantly enriched KEGG pathways were included in the network (Figure 5). The most highly expressed miRNA was ame-miR-3477-5p, which exceeded the sum of the other 14 miRNAs, and its potential target genes were significantly associated with eight pathways.
FIGURE 5. Differential expressed miRNA and KEGG pathway network. Square nodes represent miRNAs, and arrow-shaped nodes represent pathways. The miRNA nodes' sizes were determined by the expression value.
Among the 39 significantly enriched KEGG pathways, “Hippo signaling pathway-fly” was associated with the target genes of seven miRNAs. Both the “Phototransduction-fly” and “Apoptosis-fly” pathways were associated with the target genes of six miRNAs, while the “Wnt signaling pathway” and “Dorso-ventral axis formation” were associated with the target genes of five miRNAs. Therefore, “Apoptosis-fly” was also considered a key KEGG pathway, similar to the other four pathways. DEMs with the top three expression levels that participate in the regulation of at least one key pathway were determined to be key miRNAs. Therefore, ame-miR-3477-5p, ame-miR-375-3p, and miR-281-x were identified as key miRNAs, and their target genes enriched in key KEGG pathways are listed as potential key target genes (Supplementary Table S6).
Discussion
This study focused on identifying Apis mellifera ligustica brain tissue DEMs before and after fluvalinate treatment. As a pyrethroid, fluvalinate has been used to control Varroa mites in honeybees since the 1980s (Vu et al., 2020). In recent years, beekeepers have noted an increase in certain phenomena: honeybees gradually die and homing decreases after fluvalinate use, colonies generally decline, and empty nests often appear. Many studies have reported that fluvalinate has negative impacts on multiple survival traits of Apis mellifera, including learning and memory ability (Rinderer et al., 1999; Burley et al., 2008; Berry et al., 2013; Frost et al., 2013). However, to our knowledge, there is still no research on the function of miRNAs in this process.
In this study, 15 DEMs were identified at different times after fluvalinate exposure in Apis mellifera ligustica brain tissues. Among them, the five DEMs with the highest expression were ame-miR-3477-5p, ame-miR-375-3p, miR-281-x, miR-1-y, and ame-miR-281-3p. Ame-miR-3477-5p showed the highest expression value, which exceeded the sum of the other 14 miRNAs.
Expression pattern analysis showed that ame-miR-3477-5p reached the peak expression level at 10 days after fluvalinate treatment and then returned to the preadministration level. This phenomenon may suggest that ame-miR-3477-5p plays a role in the early brain damage or detoxification process caused by fluvalinate. However, we also found very few research reports related to ame-miR-3477-5p, and it has only been detected in three species: Apis mellifera, Tribolium castaneum, and Nasonia vitripennis. The only research related to ame-miR-3477-5p′s function reported that after worker honeybees perform visual pattern learning in a maze, the expression of miR-3477 in the brain is significantly upregulated (Qin et al., 2014). In our research, the potential target genes of ame-miR-3477-5p were significantly enriched in eight KEGG pathways, while the top enriched pathway was phototransduction-fly. This result suggests that the function of ame-miR-3477-5p may be related to vision, learning, or memory.
DEMs ame-miR-375-3p and miR-281-x, which were also identified as key miRNAs, showed similar expression patterns to ame-miR-3477-5p. Among them, the target genes involved in the key KEGG pathways “Phototransduction-fly” and “Wnt signaling pathway” of ame-miR-375-3p were highly similar to those of ame-miR-3477-5p. Both ame-miR-3477-5p and ame-miR-375-3p regulate the phototransduction pathway and Wnt signaling pathway mainly by targeting CaMKII (both target 42 transcripts of CaMKII). CaMKII is a multifunctional serine/threonine protein kinase that is highly expressed in brain tissue; it is also the main component of PSD and is autophosphorylated to participate in learning, memory, and cognitive ability and may be the molecular basis of nervous system plasticity (Joiner and Griffith, 1997; Bejar et al., 2002; Lisman et al., 2002; Zhang et al., 2005). Previous studies have confirmed that CaMKII is a molecular memory switch in vertebrates and insects and plays important roles in memory formation in Drosophila and Apis mellifera (Mehren and Griffith, 2004; Malik et al., 2013; Matsumoto et al., 2014; Scholl et al., 2015). Furthermore, activation of CaMKII was reported to be involved in the modulation of the voltage dependence of steady-state inactivation, and activation of the CaM-kinase II pathway by an increased intracellular calcium concentration might cause specific ion channels to be sensitive to the pesticide (Lavialle-Defaix et al., 2010). Several miRNAs, including miR-1, miR-30b-5p, miR-145, miR-148, and miR-152, have been reported to target the CaMKII gene, but miR-375-3p and miR-3477-5p are not among them (Terentyev et al., 2009; Liu et al., 2010; Cha et al., 2013; He et al., 2013). Our results reveal that ame-miR-375-3p and ame-miR-3477-5p are newly discovered potential regulators of CaMKII and may play roles in memory or cognitive impairment in the Apis mellifera ligustica brain after fluvalinate treatment. It is worth noting that the monkey DEMs bantam-x and miR-750-y were also predicted to target CaMKII, which may indicate the existence of a complex posttranscriptional regulatory network controlling CaMKII.
The key miRNA miR-281-x, which is also known as ame-miR-281-5p, has been reported to be differentially expressed in the honeybee brain after treatment with dinotefuran (Huang et al., 2021). In this research, the potential target genes of miR-281-x were significantly enriched in the key pathways Dorso-ventral axis formation and Hippo signaling pathway-fly. Our results showed that miR-281-x targets 27 transcripts of the LOC409780 gene, currently referred to as CELF2 (CUGBP Elav-like family member 2) in Apis mellifera. As a member of the CELF family and the dorsoventral axis formation pathway, CELF2 is widely expressed in various systems of animals and is implicated in the regulation of RNA biogenesis as a splicing factor (Dasgupta and Ladd, 2012; Ajith et al., 2016). MiR-281-x also targets 13 transcripts of the LOC408616 gene, currently known as TEAD1/TEF-1 (transcriptional enhancer factor TEF-1) in Apis mellifera. As a Hippo signaling effector, TEAD1 participates in various biological processes and plays important roles in animal tissue growth and maintenance (Dey et al., 2020; Heng et al., 2021).
The current results also show that the potential target genes of DEMs miR-11987-x, novel-m0194-5p, and novel-m0184-5p are enriched in the dorsoventral axis formation pathway, and both can target the Broad-Complex (Br-c) gene. As a member of the “Dorso-ventral axis formation” pathway, Br-C is an early ecdysone-responsive gene encoding a family of zinc-finger transcription factors and is involved in insect development regulation (Paul et al., 2006; Jiang et al., 2015). Br-C was reported to be relatively highly expressed in the brain tissue of honeybee and oriental river prawn, which may indicate that Br-C could participate in insect brain function regulation (Paul et al., 2006; Jiang et al., 2015). Interestingly, our results showed that DEMs miR-11987-x, novel-m0194-5p, and novel-m0184-5p have different target sites in the Br-C gene. miR-11987-x can target two transcripts, while novel-m0194-5p and novel-m0184-5p target seven different transcripts of Br-C.
In summary, the results of this study reveal that the miRNA expression profiles in the brains of A. mellifera ligustica changed before and after treatment with fluvalinate. The results of target gene functional enrichment analyses indicate that the DEMs may be related to cell apoptosis and memory damage in the brain. This study further confirms that abnormal activities observed among honeybees after fluvalinate administration may be caused by neurological disorders, including excessive arousal and memory loss caused by fluvalinate-induced brain damage. However, the specific roles of these miRNAs in brain function remain to be investigated, and their regulatory mechanisms also need to be further studied. This study provides new information on the molecular mechanisms and regulatory pathways of the fluvalinate response in A. mellifera ligustica brain tissue.
In summary, the results of this study reveal that the miRNA expression profiles in the brains of A. mellifera ligustica changed before and after treatment with fluvalinate. The results of target gene functional enrichment analyses indicate that DEMs may be related to cell apoptosis and memory damage in the brain. This study further confirms that abnormal activities in honeybees after fluvalinate administration may be caused by neurological disorders, including excessive arousal and memory loss caused by fluvalinate-induced brain damage. However, specific roles of these miRNAs in brain function remain to be investigated, and their regulatory mechanisms also need to be further studied. This study provides new information on the molecular mechanisms and regulatory pathways of fluvalinate response in A. mellifera ligustica brain tissue.
Data Availability Statement
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 below: https://www.ncbi.nlm.nih.gov/geo/, GSE162423.
Author Contributions
Conceptualization and methodology: WG and CT; experiments and data acquisition: WG, CT, YL, LD, XH, FY, HY, and SX; data analysis: YL and CT; writing: CT and YL. All authors have read and agreed to the published version of the manuscript.
Funding
This work was financially supported by the Funds of Shandong Province Agricultural Variety Project (no. 2016LZGC039) and Shandong “Double Tops” Program (no. SYL2017YSTD12). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
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.
Acknowledgments
The author would like to thank all the members of the Animal Protection Laboratory, School of Animal Science and Technology, Shandong Agricultural University, who collected the whole sample of our experiment. We are grateful to the livestock department and friends of apiculture in the sampling area for their help in sample collection on site/for providing samples. We would also thank Guangzhou Gene Denovo Biotechnology Co., Ltd. for their help.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.855987/full#supplementary-material
References
Ajith, S., Gazzara, M. R., Cole, B. S., Shankarling, G., Martinez, N. M., Mallory, M. J., et al. (2016). Position-dependent Activity of CELF2 in the Regulation of Splicing and Implications for Signal-Responsive Regulation in T Cells. RNA Biol. 13, 569–581. doi:10.1080/15476286.2016.1176663
Ashburner, M., Ball, C. A., Blake, J. A., Botstein, D., Butler, H., Cherry, J. M., et al. (2000). Gene Ontology: Tool for the Unification of Biology. Nat. Genet. 25, 25–29. doi:10.1038/75556
Bejar, R., Yasuda, R., Krugers, H., Hood, K., and Mayford, M. (2002). Transgenic Calmodulin-dependent Protein Kinase II Activation: Dose-dependent Effects on Synaptic Plasticity, Learning, and Memory. J. Neurosci. 22, 5719–5726. doi:10.1523/jneurosci.22-13-05719.2002
Berry, J. A., Hood, W. M., Pietravalle, S., and Delaplane, K. S. (2013). Field-Level Sublethal Effects of Approved Bee Hive Chemicals on Honey Bees (Apis mellifera L). PLoS One 8, e76536. doi:10.1371/journal.pone.0076536
Burley, L. M., Fell, R. D., and Saacke, R. G. (2008). Survival of Honey Bee (Hymenoptera: Apidae) Spermatozoa Incubated at Room Temperature from Drones Exposed to Miticides. J. Econ. Entomol. 101, 1081–1087. doi:10.1603/0022-0493(2008)10110.1093/jee/101.4.1081
Calin, G. A., and Croce, C. M. (2006). MicroRNA Signatures in Human Cancers. Nat. Rev. Cancer 6, 857–866. doi:10.1038/nrc1997
Cha, M.-J., Jang, J.-K., Ham, O., Song, B.-W., Lee, S.-Y., Lee, C. Y., et al. (2013). MicroRNA-145 Suppresses ROS-Induced Ca2+ Overload of Cardiomyocytes by Targeting CaMKIIδ. Biochem. Biophysical Res. Commun. 435, 720–726. doi:10.1016/j.bbrc.2013.05.050
Cristino, A. S., Barchuk, A. R., Freitas, F. C. P., Narayanan, R. K., Biergans, S. D., Zhao, Z., et al. (2014). Neuroligin-associated microRNA-932 Targets Actin and Regulates Memory in the Honeybee. Nat. Commun. 5, 5529. doi:10.1038/ncomms6529
Dasgupta, T., and Ladd, A. N. (2012). The Importance of CELF Control: Molecular and Biological Roles of the CUG-BP, Elav-like Family of RNA-Binding Proteins. WIREs RNA 3, 104–121. doi:10.1002/wrna.107
Dey, A., Varelas, X., and Guan, K.-L. (2020). Targeting the Hippo Pathway in Cancer, Fibrosis, Wound Healing and Regenerative Medicine. Nat. Rev. Drug Discov. 19, 480–494. doi:10.1038/s41573-020-0070-z
Frost, E. H., Shutler, D., and Hillier, N. K. (2013). Effects of Fluvalinate on Honey Bee Learning, Memory, Responsiveness to Sucrose, and Survival. J. Exp. Biol. 216, 2931–2938. doi:10.1242/jeb.086538
Gosselin‐Badaroudine, P., and Chahine, M. (2017). Biophysical Characterization of the Varroa Destructor Na V 1 Sodium Channel and its Affinity for τ‐fluvalinate Insecticide. FASEB j. 31, 3066–3071. doi:10.1096/fj.201601338r
Griffiths-Jones, S., Bateman, A., Marshall, M., Khanna, A., and Eddy, S. R. (2003). Rfam: An RNA Family Database. Nucleic Acids Res. 31, 439–441. doi:10.1093/nar/gkg006
Griffiths-Jones, S., Grocock, R. J., van Dongen, S., Bateman, A., and Enright, A. J. (2006). miRBase: microRNA Sequences, Targets and Gene Nomenclature. Nucleic Acids Res. 34, D140–D144. doi:10.1093/nar/gkj112
Griffiths-Jones, S. (2004). The microRNA Registry. Nucleic Acids Res. 32, 109D–111D. doi:10.1093/nar/gkh023
He, J., Jiang, S., Li, F.-l., Zhao, X.-j., Chu, E.-f., Sun, M.-n., et al. (2013). MicroRNA-30b-5p Is Involved in the Regulation of Cardiac Hypertrophy by Targeting CaMKIIδ. J. Investig. Med. 61, 604–612. doi:10.2310/JIM.0b013e3182819ac6
Heng, B. C., Zhang, X., Aubel, D., Bai, Y., Li, X., Wei, Y., et al. (2021). An Overview of Signaling Pathways Regulating YAP/TAZ Activity. Cell. Mol. Life Sci. 78, 497–512. doi:10.1007/s00018-020-03579-8
Huang, M., Dong, J., Guo, H., and Wang, D. (2021). Effects of Dinotefuran on Brain miRNA Expression Profiles in Young Adult Honey Bees (Hymenopptera: Apidae). J. Insect Sci. 21, 3. doi:10.1093/jisesa/ieaa131
Jiang, S. F., Zhang, Y. P., Sun, S. M., Gong, Y. S., Xiong, Y. W., Qiao, H., et al. (2015). Molecular Cloning, Characterization, and Expression Analysis of a Broad-Complex Homolog during Development in the oriental River Prawn Macrobrachium Nipponense. Genet. Mol. Res. 14, 5141–5152. doi:10.4238/2015.May.18.4
Joachims, T. (1999). Svmlight: Support Vector Machine. Available at: http//svmlight.joachims.org/https://akela.mendelu.cz/∼zizka/Machine_Learning/SVM/SVM-light/SVM-Light Support Vector Ma.pdf.
Johnson, R. M., Pollock, H. S., and Berenbaum, M. R. (2009). Synergistic Interactions between In-Hive Miticides in Apis mellifera. ec 102, 474–479. doi:10.1603/029.102.0202.479
Johnson, R. M., Wen, Z., Schuler, M. A., and Berenbaum, M. R. (2006). Mediation of Pyrethroid Insecticide Toxicity to Honey Bees (Hymenoptera: Apidae) by Cytochrome P450 Monooxygenases. J. Econ. Entomol. 99, 1046–1050. doi:10.1093/jee/99.4.1046
Joiner, M.-l. A., and Griffith, L. C. (1997). CaM Kinase II and Visual Input Modulate Memory Formation in the Neuronal Circuit Controlling Courtship Conditioning. J. Neurosci. 17, 9384–9391. doi:10.1523/jneurosci.17-23-09384.1997
Kanehisa, M., Araki, M., Goto, S., Hattori, M., Hirakawa, M., Itoh, M., et al. (2008). KEGG for Linking Genomes to Life and the Environment. Nucleic Acids Res. 36, D480–D484. doi:10.1093/nar/gkm882
Kruger, J., and Rehmsmeier, M. (2006). RNAhybrid: MicroRNA Target Prediction Easy, Fast and Flexible. Nucleic Acids Res. 34, W451–W454. doi:10.1093/nar/gkl243
Lavialle-Defaix, C., Moignot, B., Legros, C., and Lapied, B. (2010). How Does Calcium-dependent Intracellular Regulation of Voltage-dependent Sodium Current Increase the Sensitivity to the Oxadiazine Insecticide Indoxacarb Metabolite Decarbomethoxylated JW062 (DCJW) in Insect Pacemaker Neurons? J. Pharmacol. Exp. Ther. 333, 264–272. doi:10.1124/jpet.109.163519
Lewis, B. P., Burge, C. B., and Bartel, D. P. (2005). Conserved Seed Pairing, Often Flanked by Adenosines, Indicates that Thousands of Human Genes Are microRNA Targets. Cell 120, 15–20. doi:10.1016/j.cell.2004.12.035
Lisman, J., Schulman, H., and Cline, H. (2002). The Molecular Basis of CaMKII Function in Synaptic and Behavioural Memory. Nat. Rev. Neurosci. 3, 175–190. doi:10.1038/nrn753
Liu, X., Zhan, Z., Xu, L., Ma, F., Li, D., Guo, Z., et al. (2010). MicroRNA-148/152 Impair Innate Response and Antigen Presentation of TLR-Triggered Dendritic Cells by Targeting CaMKIIα. J.I. 185, 7244–7251. doi:10.4049/jimmunol.1001573
Malik, B. R., Gillespie, J. M., and Hodge, J. J. L. (2013). CASK and CaMKII Function in the Mushroom Body α′/β′ Neurons during Drosophila Memory Formation. Front. Neural Circuits 7, 52. doi:10.3389/fncir.2013.00052
Matsumoto, Y., Sandoz, J.-C., Devaud, J.-M., Lormant, F., Mizunami, M., and Giurfa, M. (2014). Cyclic Nucleotide-Gated Channels, Calmodulin, Adenylyl Cyclase, and Calcium/calmodulin-dependent Protein Kinase II Are Required for Late, but Not Early, Long-Term Memory Formation in the Honeybee. Learn. Mem. 21, 272–286. doi:10.1101/lm.032037.113
Medici, S. K., Castro, A., Sarlo, E. G., Marioli, J. M., and Eguaras, M. J. (2012). The Concentration Effect of Selected Acaricides Present in Beeswax Foundation on the Survival of Apis mellifera Colonies. J. Apicultural Res. 51, 164–168. doi:10.3896/IBRA.1.51.2.03
Mehren, J. E., and Griffith, L. C. (2004). Calcium-independent Calcium/calmodulin-dependent Protein Kinase II in the Adult Drosophila CNS Enhances the Training of Pheromonal Cues. J. Neurosci. 24, 10584–10593. doi:10.1523/JNEUROSCI.3560-04.2004
Mello, T. R. P., Aleixo, A. C., Pinheiro, D. G., Nunes, F. M. F., Bitondi, M. r. M. G., Hartfelder, K., et al. (2014). Developmental Regulation of Ecdysone Receptor (EcR) and EcR-Controlled Gene Expression during Pharate-Adult Development of Honeybees (Apis mellifera). Front. Genet. 5, 445. doi:10.3389/fgene.2014.00445
Paul, R. K., Takeuchi, H., and Kubo, T. (2006). Expression of Two Ecdysteroid-Regulated Genes,Broad-ComplexandE75, in the Brain and Ovary of the Honeybee (Apis mellifera L.). Zoolog. Sci. 23, 1085–1092. doi:10.2108/zsj.23.1085
Qin, Q.-H., Wang, Z.-L., Tian, L.-Q., Gan, H.-Y., Zhang, S.-W., and Zeng, Z.-J. (2014). The Integrative Analysis of microRNA and mRNA Expression inApis Melliferafollowing Maze-Based Visual Pattern Learning. Insect Sci. 21, 619–636. doi:10.1111/1744-7917.12065
Rinderer, T. E., De Guzman, L. I., Lancaster, V. A., Delatte, G. T., and Stelzer, J. A. (1999). Varroa in the Mating Yard: I. The Effects of Varroa Jacobsoni and Apistan on Drone Honey Bees. Am. Bee J. 139, 134–139.
Scholl, C., Kübert, N., Muenz, T. S., and Rössler, W. (2015). CaMKII Knockdown Affects Both Early and Late Phases of Olfactory Long-Term Memory in the Honeybee. J. Exp. Biol. 218, 3788–3796. doi:10.1242/jeb.124859
Saeed, A. I., Sharov, V., White, J., Li, J., Liang, W., Bhagabati, N., et al. (2003). TM4: A Free, Open-Source System for Microarray Data Management and Analysis. BioTechniques 34 (2), 374–378. doi:10.2144/03342mt01
Soderlund, D. M. (2010). State-dependent Modification of Voltage-Gated Sodium Channels by Pyrethroids. Pestic. Biochem. Physiol. 97, 78–86. doi:10.1016/j.pestbp.2009.06.010
Terentyev, D., Belevych, A. E., Terentyeva, R., Martin, M. M., Malana, G. E., Kuhn, D. E., et al. (2009). miR-1 Overexpression Enhances Ca 2+ Release and Promotes Cardiac Arrhythmogenesis by Targeting PP2A Regulatory Subunit B56α and Causing CaMKII-dependent Hyperphosphorylation of RyR2. Circ. Res. 104, 514–521. doi:10.1161/CIRCRESAHA.108.181651
Turner, D. A. (1985). Miranda: A Non-strict Functional Language with Polymorphic Types. Lecture Notes Comp. Sci. (including subseries Lecture Notes Artif. Intelligence Lecture Notes Bioinformatics) 201, 1–16. doi:10.1007/3-540-15975-4_26
Vu, P. D., Rault, L. C., Jenson, L. J., Bloomquist, J. R., and Anderson, T. D. (2020). Voltage-gated Chloride Channel Blocker DIDS as an Acaricide for Varroa Mites. Pestic. Biochem. Physiol. 167, 104603. doi:10.1016/j.pestbp.2020.104603
Wang, R., Xu, L., Zhu, X., Zhai, L., Wang, Y., Yu, R., et al. (2015). Transcriptome-Wide Characterization of Novel and Heat-Stress-Responsive microRNAs in Radish (Raphanus Sativus L.) Using Next-Generation Sequencing. Plant Mol. Biol. Rep. 33, 867–880. doi:10.1007/s11105-014-0786-1
Yu, G., Wang, L.-G., Han, Y., and He, Q.-Y. (2012). ClusterProfiler: An R Package for Comparing Biological Themes Among Gene Clusters. OMICS: A J. Integr. Biol. 16, 284–287. doi:10.1089/omi.2011.0118
Keywords: Apis mellifera ligustica, fluvalinate, brain, transcriptome, MicroRNA
Citation: Tianle C, Liuxu Y, Delong L, Yunhan F, Yu H, Xueqing S, Haitao X and Guizhi W (2022) Fluvalinate-Induced Changes in MicroRNA Expression Profile of Apis mellifera ligustica Brain Tissue. Front. Genet. 13:855987. doi: 10.3389/fgene.2022.855987
Received: 16 January 2022; Accepted: 17 March 2022;
Published: 12 April 2022.
Edited by:
Gary S. Stein, University of Vermont, United StatesReviewed by:
Pei Fan, Henan University of Technology, ChinaMichele Mortarino, University of Milan, Italy
Copyright © 2022 Tianle, Liuxu, Delong, Yunhan, Yu, Xueqing, Haitao and Guizhi. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Chao Tianle, chaotianle@sdau.edu.cn; Wang Guizhi, wgz@sdau.edu.cn