The Comparison of Juvenile Hormone and Transcriptional Changes between Three Different Juvenile Hormone Analogs Insecticides on Honey Bee Worker Larval’s Development

Juvenile Hormone and Transcriptional Changes between Three Different Juvenile Hormone Analogs Insecticides on Honey Bee Abstract: Juvenile hormones (JHs) play a crucial role in the development of honey bee ( Apis mellifera ) worker larvae. Juvenile hormone analogs (JHAs), insecticides widely used in pest control, have been reported to affect the health and survival of honey bee worker larvae. However, the molecular mechanisms of JHAs in the honey bee remain unclear. In this study, we treated honey bee worker larvae with pyriproxyfen, fenoxycarb, and methoprene, three different JHAs. We monitored the changes in the transcription of genes encoding major JH response enzymes (CYP15A1, CYP6AS5, JHAMT, and CHT1) using RT-qPCR and analyzed the transcriptome changes in worker larvae under JHA stress using RNA-seq. We found that the enrichment pathways differed among the treatment groups, but the classiﬁcation of each pathway was generally the same, and fenoxycarb affected more genes and more pathways than did the other two JHAs. Notably, treatment with different JHAs in the honey bee changed the JH titers in the insect to various extents. These results represent the ﬁrst assessment of the effects of three different JHAs on honey bee larvae and provide a new perspective and molecular basis for the research of JH regulation and JHA


Introduction
The honey bee is widely distributed all over the world [1]. It plays a crucial role in the pollination of crops [2,3], and a healthy and strong colony of honey bees well enhances the pollination area of local crops. Thus, the honey bee greatly contributes to the social economy and ecology. The honey bee has developed into a new model organism, which, because of its extremely complex and delicate social structure and variety of behaviors, allows for the biology and behaviors of insects to be studied [4][5][6][7]. Moreover, exploring the effects of insecticides on the honey bee has always been an important component in honey bee science.
Juvenile hormones (JHs) are a group of sesquiterpenoids [8] that are indispensable endocrine hormones in insects. They are synthesized and secreted by the corpora allata [9,10] and play a crucial role in the regulation of vitellogenesis, reproduction, and ovary development of insects [11][12][13]. They also function in insect migration, energy metabolism, and caste differentiation [14][15][16]. JHAs are compounds similar to JH in structure and mechanism [17][18][19][20]. JHAs can effectively mimic JHs in insects to exercise its function or interfere with JH synthesis. Consequently, changes occur in JH titers in the insect hemolymph, thereby interfering with the normal physiological and biochemical processes of insects and achieving the purpose of controlling and eliminating pests because of their similarities [21]. JHA insecticides include pyriproxyfen, fenoxycarb, and methoprene. Previous research has shown that pyriproxyfen can reduce the growth rate of Hippodamia convergens and inhibit their population. Fenoxycarb reduces the longevity of pear psylla (Cacopsylla pyricola). Methoprene is an important mediator of female Drosophila melanogaster, which, when treated with methoprene, exhibits sexual inhibition to male D. melanogaster after mating [22][23][24]. Other studies have suggested that JHA insecticides can affect the social activity of the honey bee and change its growth and development processes [25,26]. Genes such as cytochrome CYP6AS5 (CYP6AS5), chitinase (CHT1), cytochrome CYP15A1 (CYP15A1), and juvenile hormone acid methyltransferase (JHAMT) reportedly induce responses when insects are under JHA stress [27][28][29]. Currently, with the wider use of JHA insecticides in pest control, it is becoming increasingly important to test the toxicity of pyriproxyfen, fenoxycarb, and methoprene in honey bees as well as their regulation and molecular mechanisms.
In the present study, changes in the JH titers inside the honeybees' hemolymphs were measured after treatment with three different JHAs. We also present an analysis of the treatment transcriptome, for which high-throughput RNA-seq technology was utilized to identify and quantify the expression levels of the genes transcribed in honey bees' larvae, which had been exposed to three different JHAs. The genes related to metabolism processes and physiological changes were then identified and analyzed. Moreover, the differentially expressed genes (DEGs) possibly involved in the physiological metabolism and detoxification responses were identified and verified using the RT-qPCR method. The resulting transcriptome library may help to elucidate the regulation and molecular mechanisms of JHA on honey bees and screen the most effective JHA pesticide among these three JHAs in honey bees.

Chemicals and Honey Bee Larvae Rearing
Pyriproxyfen, fenoxycarb, methoprene standard, JH III, HPLC-grade methanol, acetonitrile, and isooctane were purchased from ANPEL Laboratory Technologies (Shanghai) Inc. Cell culture plates, D-glucose, D-fructose, yeast extract, and other reagents (analytical reagents) were obtained from Sangon Bioengineering (Shanghai Co., Ltd., Shanghai China). Royal jelly was supplied by a bee keeper in Nanning.
Our experiments were conducted during the winter of 2020 at the GuangXi University Insect Application Laboratory (22 • 50 N, 108 • 17 E), Nanning Guangxi, China. All honey bee larvae samples were collected from A. mellifera ligustica colonies. Three treatment groups and one control group were established, and three different JHAs were used separately in different treatment groups.
In order to obtain precisely aged larvae, a queen was captured and locked in an empty comb for 4 h for egg laying. Newly emerged larvae (<12 h old) were transferred to a clean 48-well cell culture plates, which were placed in sterile queen cell cups containing 25 µL of diet A, as shown in Table 1. All plates were placed in an incubator (34.5 • C, 95% Relative Humidity), and the day of larval removal was defined as day 1; no diet was served on day 2. A total of 25 µL of diet B was provided on day 3, and 30 and 40 µL of diet C were provided on days 4 and 5, respectively. All our rearing methods were conducted following the standard beekeeping methods for breeding [30]. The concentration of each JHA for treatment was previously described by Tasei and Milchreit [31,32]. In this report, the maximum fenoxycarb concentration found in fieldrealistic doses ranged between 7.5 and 217 ng/µL. Our former research (not published) shows that most honey bee larvae will not survive when they are exposed to 100 ng/µL of these three JHAs. Here, we chose 50 ng/µL for each JHA (dissolved in 0.1% v/v DMSO) [32]. Generally, different treatments and times may lead to discrepant conclusions, even when experimental subjects are exposed to identical concentrations in toxicology tests. Thus, a necessary pre-experiment was conducted to demonstrate that there was no apparent difference in mortality between control and 50 ng/µL of each JHA-treated group. According to the sequence of pyriproxyfen, fenoxycarb, methoprene, their p-values = 0.230, 0.422, and 0.768, respectively, but all of their adjusted p > 0.05 (independent sample t-test) (Table S1). Hence, it was considered suitable to apply a 50 ng/µL concentration of these three JHAs to honey bees in our experiments. We used 1 µL of different JHAs (dissolved in 0.1% v/v DMSO) and added them directly to diet C at a concentration of 50 ng/µL on day 4 (50 ng/bee); the control group was provided with a total of 1 µL DMSO (0.1% v/v) without any JHAs in diet C, and only healthy six-day-old (48 h after day 4) larvae (raised using our method and no food left in the plate) were used for future tests. Three biological replicates were prepared for each group. All six-day-old bee larva samples were used for the extraction of JH, RT-qPCR experiments, and transcriptome analysis.

Extraction and Determination of JH Titers in Honey Bees
JH III is known as the only hormone found in the worker honey bee [33]. Herein, the method used to extract JH III was based on the method described by Furuta [34]. Six-day-old honey bee samples weighing 1 g (about 11 6-day-old larvae; the actual weight of each repeat was recorded for later JH III level calculations) were anesthetized at a low temperature, placed in a grinding pestle, and cooled by liquid nitrogen. After adding liquid nitrogen, the honey bees were continuously ground into powder by automatic grinding. The bee powder was transferred into a 10 mL centrifuge tube, and then 0.2 mL of acetonitrile and 1 mL of methanol were added after grinding was completed. The sample was mixed in a shaker for 3 min to which 2 mL of isooctane was added, and it was shaken again for another 3 min for complete mixing. Centrifugation was performed at 4000 rpm for 5 min, after which the supernatant was carefully poured into a new centrifuge tube. The isooctane extraction and supernatant absorption steps were repeated twice. Finally, a supernatant of about 6 mL was dried with nitrogen-blowing apparatus. After drying, 1 mL of acetonitrile (HPLC-grade) was added to the centrifuge tube, and the solution was passed through a green organic membrane (purification equipment). Finally, 0.2 mL of solution was poured into a small sample bottle for future determination.
Agilent 1260 series high-performance liquid chromatography-tandem mass spectrometry (Agilent Technologies, Inc., Colorado Springs, CO, USA) was used for chromatographic analysis. The chromatographic column was CHIRALPAK AY-3R (150 mm × 2.1 mm, 3 µm; Daicel Chiral Technologies, Shanghai, China). The column temperature was 30 • C, the mobile phase was 0.1% formic acid and acetonitrile (v:v = 55:45), the flow rate was 0.3 mL/min, the injection volume was 3 µL, the mobile phase was 0.1% formic acid water-acetonitrile solution, and the detection wavelength was 214 nm. The mass spectrometry conditions were as follows: electrospray ion source and positive ion scanning mode (ESI+); capillary voltage, 4.0 kv; air temperature, 330 • C; and atomizing gas pressure, 103.425 kPa. Multiple response monitoring mode was used with isocratic elution. Qualitative analysis was carried out based on retention time and ion information comparison. The parent ion and ion with the highest response value were quantitatively analyzed for the R-JH III and S-JH III retention time, ion monitoring, pyrolysis voltage, and collision energy.

RNA Extraction, Library Preparation, and Sequencing
The total RNA of honey bee larvae was extracted with a Trizol (MRC Cincinnati, OH, USA) reagent kit according to the manufacturer's instructions. The biological repli-cates of pyriproxyfen, fenoxycarb, and methoprene treatment groups were pyriproxyfen-1, pyriproxyfen-2, and pyriproxyfen-3; fenoxycarb-1, fenoxycarb-2, and fenoxycarb-3; and methoprene-1, methoprene-2, and methoprene-3, respectively. The biological replicates of the control group were control-1 (CK-1), control-2 (CK-2), and control-3 (CK-3). For all groups, three specimens were collected for each sample, three replicates were performed, and the quantity and quality of RNA were assessed using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific). The used reverse transcription kit was Hifair ® III 1st Stand cDNA Synthesis Supermix for qPCR (gDNA digester plus) (Yeasen Biotechnology (Shanghai) Co., Ltd., Shanghai, China). Quantitative real-time PCR was carried out using a QuantStudio 6 Flex instrument (Thermo Fisher Scientific) at 95 • C for 30 s for pre-denaturation and 40 cycles of 95 • C for 5 s and 60 • C for 30 s for PCR. The total reaction volume was 10 mL. Primer sequences for genes encoding JHA stress enzymes and internal controls were designed using Primer-BLAST (NCBI). RNA samples were prepared for a subsequent study.
RNA samples were collected for cDNA library construction. Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) and Qubit fluorometer (Invitrogen, Carlsbad, CA, USA) were used to determine the RNA quantity and quality, respectively. Three replicates each of the control and treatment groups were analyzed. Total RNA extraction steps were the same. Library preparation and whole-transcriptome sequencing were performed on the Illumina sequencing platform developed by Genedenovo Biotechnology Co., Ltd. (Guangzhou, China). RNA integrity was checked using standard agarose gel electrophoresis with TAE buffer. First-strand cDNA was synthesized in the M-MuLV reverse transcriptase system with fragmentized mRNA as a template and random oligonucleotides as primers. The RNA strand was then degraded with RNaseH. Second-strand cDNA was synthesized from dNTPs in the DNA polymerase I system. Random hexamer primers were used for library preparation and amplification. The AMPure XP System (Beckman Coulter, Brea, CA, USA), was used to purify double-stranded cDNA, which was then eluted with EB buffer for end repair and the addition of poly (A). Illumina HiSeqTM2000 (Illumina, San Diego, CA, USA) was used to sequence the cDNA library and generate about 200 bp paired-end reads.

Transcriptomic Data Analysis
Raw data (raw reads) in FASTA format were processed through in-house perl scripts (NCBI number: PRJNA734139). In this step, clean data (clean reads) were obtained by removing reads containing an adapter, reads with all A bases, reads containing ploy-N, and low-quality reads from raw data. Meanwhile, the contents of GC, Q20, and Q30 of clean data were calculated. All downstream analyses were based on high-quality clean data and HISAT2 was used to count the read numbers mapped to each gene. The FPKM of each gene was then calculated based on the gene length and read count mapped to this gene.

Functional Analyses of DEGs
The relative expression levels of genes among the control, pyriproxyfen-treated, fenoxycarb-treated, and methoprene-treated groups were analyzed using the DEGseq method [35], which provided statistical routines to determine differential expression in digital gene expression data by using a model based on the negative binomial distribution. Corrected p-value < 0.05 and log2 (fold change) ≥1 were set as the threshold for significantly differential expression. GO terms with corrected p-values less than 0.05 were considered as significantly enriched by DEGs. KOBAS software was used to determine the statistical enrichment of DEGs in Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways.

Validation of Four DEGs
Four candidate differentially regulated JHA-sensitive genes with various biological functions were analyzed by RT-qPCR by using three biological replicates to verify the transcriptome results. We selected four DEGs related to JHA stress, CHT1, CYP6AS5, CYP15A1, and JHAMT genes, which were expressed in three different treatments (Table S2). About 10 µL of PCR reaction volume contained 5 µL of TB Green TM Pre mix Ex taqTMII, 3 µL of DEPC water, 1 µL of diluted cDNA template, and 0.5 µL of each primer. The PCR reaction was performed on the Step One Plus TM Real-Time System (Bio Rad, USA) under the following conditions: 95 • C for 30 s for the holding stage; 40 cycles of 95 • C for 5 s, 60 • C for 34 s for the PCR stage. A final melting curve analysis was performed. The relative expression levels of the selected transcripts normalized to the housekeeping gene (actin) were calculated using the 2 −∆∆CT method [36]. ∆∆Ct = (Ct, Target − Ct, Actin) treat − (Ct, Target − Ct, Actin) control.

Statistics
Student t-tests were used to determine the significance of DEGs expression differences between control and treatment groups (RNA-seq and RT-qPCR). Two-tailed probabilities were adopted in the tests. Multistrata parametric tests were used to determine the significance of JH III titer changes between control and treatment groups. All analyses were performed with SPSS Statistics 21.0 software (SPSS Inc., Chicago, IL, USA).

Determination of JH Titers in Six-Day-Old Worker Larvae
The JH content in honey bee larvae after different JHA treatments was determined by LC-MS. Three kinds of JHA changed the JH titers in honey bee larvae. Among them, the JH titers in honey bee larvae treated with fenoxycarb were the lowest. The JH titers treated with methoprene were slightly higher than those treated with fenoxycarb. The JH titers treated with pyriproxyfen were the highest (Figure 1).
Agronomy 2021, 11, x FOR PEER REVIEW 6 Figure 1. Titers of JH Ⅲ in 6-days-old larvae after different JHAs exposure and in the control group. Data presented in the table are mean ± standard error. X-axis represents the control group and different-treated groups. Y-axis shows the titers level of JH Ⅲ. With significant differences between different treatment groups and control groups are indicated using asterisks: * p < 0.05, ** p < 0.01 (post-hoc parametric test) (Table S3).

Quality Analysis of Transcriptome Sequencing Results
After filtering low-quality data, clean reads were obtained. The effective read RNA-seq ranged from 40,097,124 to 58,162,002, with the lowest unique matching rat 85.10% and the highest 88.98%. Q30 values were also above 90.14% (Table S4). This find showed that the quality of the RNA-seq data was high. The final sequencing results w Data presented in the table are mean ± standard error. X-axis represents the control group and different-treated groups. Y-axis shows the titers level of JH III. With significant differences between different treatment groups and control groups are indicated using asterisks: * p < 0.05, ** p < 0.01 (post-hoc parametric test) (Table S3).

Quality Analysis of Transcriptome Sequencing Results
After filtering low-quality data, clean reads were obtained. The effective reads of RNA-seq ranged from 40,097,124 to 58,162,002, with the lowest unique matching rate of 85.10% and the highest 88.98%. Q30 values were also above 90.14% (Table S4). This finding showed that the quality of the RNA-seq data was high. The final sequencing results were reliable.

Screening of DEGs
Based on the results of differential analyses, we screened the false discovery rate (FDR) < 0.05 and log2FC > 1 as a DEGs. A total of 34, 532, and 222 DEGs (Table S5) were identified among the three different JHA treatment and control groups. Among them, 21 were upregulated and 13 were downregulated in the pyriproxyfen-treated group compared with the control group; 222 were upregulated and 310 were downregulated in the fenoxycarb-treated group compared with the control group; and 155 were upregulated and 67 were downregulated in the methoprene-treated group compared with the control group ( Figure 2A). Moreover, the results for the DEGs as visualized via a Venn diagram ( Figure 2B) revealed that 34, 532, and 222 genes were specifically expressed in the different JHA treatment groups, and they had two DEGs in common.

GO Functional Annotation and KEGG Pathway
GO annotations were classified into the category of biological processes, molecular functions, and cellular components. When the p-value was <0.05, it was considered as an enrichment item. Among the DEGs in the pyriproxyfen-treated group, 91 genes with GO annotation were enriched to 20 GO terms; 1567 gens with GO annotation were enriched to 38 GO terms in the fenoxycarb-treated group; and 796 genes with GO annotation were enriched to 33 GO terms in the methoprene-treated group. The GO annotations for the DEGs of different JHA treatment groups are shown in Figure S1. The top 20 significantly enriched GO terms of each JHA treatment group are shown in Figure S2.
The biochemical pathways of the DEGs were investigated using the KEGG database. KEGG analysis suggested that 49 DEGs were divided into 31 pathways in the pyriproxyfen-treated group. In the fenoxycarb-treated group, a total of 601 DEGs were assigned to 227 pathways. In the methoprene-treated group, 409 DEGs were involved in 171 pathways. The top 20 significantly enriched KEGG pathways are shown in Figure S3.

Validation of DEGs by Using RT-qPCR
To verify the reliability of the RNA-Seq data, four DEGs involved in the JH response were selected and verified using the RT-qPCR method. The results showed a similar trend of expression changes for RT-qPCR and RNA-Seq data, indicating that the RNA-Seq results were reliable (Figure 3), and the FPKM mean value of RNA-seq and RT-qPCR data had been shown in Tables 2 and 3.

GO Functional Annotation and KEGG Pathway
GO annotations were classified into the category of biological processes, molecular functions, and cellular components. When the p-value was <0.05, it was considered as an enrichment item. Among the DEGs in the pyriproxyfen-treated group, 91 genes with GO annotation were enriched to 20 GO terms; 1567 gens with GO annotation were enriched to 38 GO terms in the fenoxycarb-treated group; and 796 genes with GO annotation were enriched to 33 GO terms in the methoprene-treated group. The GO annotations for the DEGs of different JHA treatment groups are shown in Figure S1. The top 20 significantly enriched GO terms of each JHA treatment group are shown in Figure S2.
The biochemical pathways of the DEGs were investigated using the KEGG database. KEGG analysis suggested that 49 DEGs were divided into 31 pathways in the pyriproxyfentreated group. In the fenoxycarb-treated group, a total of 601 DEGs were assigned to 227 pathways. In the methoprene-treated group, 409 DEGs were involved in 171 pathways. The top 20 significantly enriched KEGG pathways are shown in Figure S3.

Validation of DEGs by Using RT-qPCR
To verify the reliability of the RNA-Seq data, four DEGs involved in the JH response were selected and verified using the RT-qPCR method. The results showed a similar trend of expression changes for RT-qPCR and RNA-Seq data, indicating that the RNA-Seq results were reliable (Figure 3), and the FPKM mean value of RNA-seq and RT-qPCR data had been shown in Tables 2 and 3. Figure 3. Confirmation of RNA-Seq results by RT-qPCR. Data presented in the figure are mean ± standard error. X-axis represents different treatment groups. Y-axis shows the relative level of gene expression. With significantly expression differences between different treatment groups and control groups are indicated using asterisks: ** p < 0.01 (independentsample t-test).

Discussion
JHs are important hormones in honey bees, as they regulate the growth and development of honey bee larvae. JHs can increase food consumption and reduce fat amounts in worker and queen bumble bees (Bombus terrestris) [37]. Former researchers discovered that there are JH titer differences between honey bee workers and queens [38,39], and they also proved that JHs are indispensable hormones in the developmental determination of queen-like characters [40,41]. In our study, the application of JHAs at environmental concentrations affected the JH titers of honey bees. Bee larvae with such hormonal changes probably fail to molt or become abnormal workers or queens. In addition, we do not know how long the negative effects of this elevated JH level will last and whether it will continue to affect bees after metamorphosis; thus, more work studying the effects of JH III abnormalities induced by JHAs in the larva stage on the life cycle of adults is necessary.
The synthesis and metabolism of JHs in insects are well understood, but the mechanism regulating their levels remains unclear. In our experiment, when exposure to different JHAs downregulated JH III levels, the RT-qPCR result of all related metabolic enzymes showed the same expression patterns as those of the transcriptome results, and the RT-qPCR results of some genes showed significant changes when compared with their transcriptome results in different treatments. Our results indicate that the expressions of CYP15A1 and JHAMT in different treatments were upregulated. As important enzymes related to JH metabolism [41][42][43], CYP15A1 and JHAMT should increase JH titers in insects when they are upregulated, but our results suggest that the JH titers show downregulation when CYP15A1 and JHAMT are upregulated. This may be related to the expression of other metabolic enzymes of JHs, such as juvenile hormone esterase and juvenile hormone epoxide hydrolase (Table S6). In addition, hormone regulation in an insect is complex and involves a large number of factors. Ion levels, insulin signaling [44], and allatostatins [9,45] are also considered influencing factors for JH synthesis.
Cytochrome P450 enzymes are well-known for their metabolism of pesticides [46,47], not only CYP15A1 as discussed above but also the CYP6AS5. In our experiment, we demonstrated that the CYP6AS subfamily is unique to hymenopterans. Wang et al. [48] showed that carbendazim decreases JH titer levels in honey bee larvae, accompanied by the downregulation of CYP6AS5, which shares the same trend with our results. A previous study conducted by our group demonstrated that CYP6AS5 is upregulated when honey bees are exposed to thiamethoxam. In addition, JH levels show an increasing trend [49]. Thus, a positive correlation between the level of JHs in insect larvae and CYP6AS5 transcription in honey bees was observed. We speculate that CYP6AS5 is consistent with CYP15A1 and CYP4C7, and that it plays an important role in regulating JH titers [50,51].
There are large differences in the number of DEGs in different treatment groups. Figures S2 and S3 indicate that the top 20 significantly enriched GO terms and KEGG pathways were different among the different treatments. This result shows that the effects of the molecular mechanisms of pyriproxyfen, fenoxycarb, and methoprene on honey bees are different, which may explain the results of JH titers. Figure S1 shows that fenoxycarb affected more genes and more pathways than did the other two JHAs, and it also showed a stronger effect on the control JH titers.
In conclusion, the results of this study show that the environmental concentration of these three JHAs significantly reduced the level of JH III in honey bee worker larvae, which may harm honey bee caste differentiation. The data of our RNA-seq results reveal that these JHAs affected the catalytic activity, biosynthetic process, and metabolism pathways. We speculate that some of these alternations are related to the regulation of JH levels. Our results provide a new perspective and molecular basis for the study of the effect of JHA mechanisms on honey bees.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/10 .3390/agronomy11122497/s1. Figure S1: Correlation coefficient of samples, a histogram representation of the GO terms in different treatment groups. Figure S2: Top 20 significantly enriched GO terms of different treatments. Figure S3: Top 20 significantly enriched KEGG pathways of different treatments. Table S1: Survival rate. Survival rates of control and 50 ng/ul of each treatment groups (t-Tests). Table S2: Gene-specific primers. Specific primers of functional gene used for RT-qPCR.

Conflicts of Interest:
The authors declare that there is no competing interest and that the article is submitted without any commercial or economic interest that could be generated as a potential conflict of interest.