Animals and experimental design
Three honeybee (Apis mellifera) colonies were used in this study, and each colony had a mated queen and about 30,000 bees. They were kept in the Honeybee Research Institute, Jiangxi Agricultural University, China.
Queens were caged onto a drone frame to lay eggs for 6 h. After hatching for 6 h, drone larvae were moved to worker cells and commercial plastic queen cells, and the rest of larvae remained in drone cells. We transplanted these drone larvae to drone cells after feeding with royal jelly in QCs for 48 h. The DC and QC drones finally were emerged from drone cells, and the WC drones were emerged from WCs.
We collected 3rd drone instars (around 54 h old larvae) of QCs, WCs and DCs, respectively, for RNA-SEq. Each sample contained 6 larvae, and each group had 3 biological replicates from three colonies. Newly emerged drones from these three groups were collected for morphological measurement and RNA-SEq. For RNA-Seq, each sample contained 2 newly emerged drones, and each group had three biological replicates. All samples for RNA-Seq were put into liquid nitrogen, frozen for 30 min and then stored at -80℃ before RNA extraction. For morphological measurement, each group had 15 newly emerged drones from three colonies. For weight measurement, each group had 90 replicates.
Morphological Measurement
Morphological data included birth weight, length and width of drone wings, width of thoraxes and horizontal area of head. Drone frames were placed into an incubator for 12 h before emerging. Newly emerged drones were captured and weighed using an analytical balance (0.01 mg, AUY120, Shimazu Co. Ltd., Japan). Subsequently, drones were anesthetized by CO2, and their wing length and width, thorax width and head horizontal area were measured by a microscopic imaging system (microscope: GL99TI, Guiguang Co. Ltd., China; Troup view software: x64, ToupTek Co. Ltd., China). The data of head horizontal area were obtained by Image J (1.52a, Wayne Rasband National Institutes of Health, USA) using images taken by the above microscopic imaging system (microscope: GL99TI, Guiguang Co. Ltd., China; Troup view software: x64, ToupTek Co. Ltd., China).
The drone reproductive organ, including seminiferous tubules length, horizontal area of seminal vesicles and mucous glands, were measured. Firstly, the drones were dissected, and the testes, seminal vesicles and mucous glands were taken out and washed with saline solution and placed on glass slides. After removing redundant tissues, the location of the tissues was adjusted on glass slides, and the images were taken under the microscopic imaging system. The tunica testis of the testes was removed, and seminiferous tubules were divided into many small aggregations according to the method developed by Tavares et al. [20]. The horizontal area of seminal vesical and mucous glands was measured by Image J, and the length of seminiferous tubules was measured by Troup view.
Rna Extraction And Sequencing
Total RNA of newly emerged drones and 3rd instars of those three groups was extracted in accordance with the standard protocol of the TRlzol Reagent (Life Technologies, CA, USA), respectively. RNA integrity and concentration were tested by an Agilent 2100 Bioanalyzer (Agilent Technologies, Inc., CA, USA).
RNA-Seq was done according to our previous study [36]. Briefly, mRNA of each sample was isolated from total RNA by a NEBNext Poly (A) mRNA Magnetic Isolation Module (NEB, E7490; New England Biolabs Inc., USA) and broken randomly using Fragmentation Buffer. Afterward, a cDNA paired-end library of each sample was constructed using a NEBNext Ultra RNA Library Prep Kit (NEB, E7530; New England Biolabs Inc.,USA) and the NEBNext Multiplex Oligos (NEB, E7500 ; New England Biolabs Inc.,USA). Purified double-stranded cDNA was isolated by AgencourtAMPure XP beads (Beckman Coulter, Inc.) and acquired by PCR. The effective concentration of the library (> 2 nM) was accurately quantified by qRT-PCR to ensure the quality of the library. All constructed cDNA libraries were sequenced by the HiSeq X Ten sequencing platform (Illumina Inc., USA).
Analysis Of Pearson`s Correlation
The reliability of biological replicates in each group was evaluated using Pearson’s correlation coefficient analysis. The rate of Pearson’s correlation coefficient over 0.8 was considered as a conventionally accepted threshold for valid replicates [37]. One 3rd drone instar sample with a very low rate of Pearson’s correlation coefficient (< 0.8) was removed before gene expression analysis (Fig. S1).
Co-expression Network Analysis
WGCNA was done using R package [38] (4.0.2) to determine the correlations of gene expression among the three drone groups. The read count of each gene from each sample was used, and the WGCNA method was done according to Langfelder and Horvath [38]. The gene clustering tree from this analysis is shown in Fig. 3.
Gene Expression Analysis
Like our previous study [36], low-quality reads were filtered and those with a sequencing error rate of less than 1% (Q20 > 98%) were retained. The clean reads were mapped to the newest version of honeybee genomics (Amel_HAv3.1). Gene expression levels were estimated using fragments per kilobase of exon per million fragments mapped (FPKM) values by the Cufflinks software [37].
DESeq2 [39] was used to analyze the differential expression among three drone groups using gene read counts rather than FPKM values. Fold Change ≥ 1 and p < 0.05 were used as the screening criteria to identify significant DEGs among three groups. We selected 61 and 50 interesting genes that are involved in caste differentiation and development regulation [9, 25–30] from comparisons at larval and adult stages respectively. The log10 fold change values of these genes were used for heatmap analysis in R package [40] (4.0.2), and results were shown in Figs. 4 and 5.
Enrichment Of Go And Kegg
Sequences of DEGs from all comparisons were against various protein and nucleotide sequence databases by BLASTX, including the National Center for Biotechnology Information (NCBI) non-redundant protein (Nr) database and Swiss-Prot database and non-redundant nucleotide sequence (Nt) database with a cut-off E-value of 10− 5. DEGs were mapped in terms of the GO database and a hypergeometric test [41] (p < 0.05 indicates the significance) in GO enrichment analysis to identify their significantly enriched GO terms.
Similarly, all DEGs from each comparison were mapped to the KEGG protein database (http://www.genome.jp/kegg/kegg1.html) using BLAST (E-value < 1e-5). The statistical enrichment of DEGs in KEGG pathways was analyzed by a hypergeometric test (Q-value < 0.05) using the KOBAS 2.0 software [42].
Qrt-pcr Analysis Of Ten Selected Genes
Total RNA of 3 d larvae and newly emerged drones from DCs, QCs and WCs were extracted and used in qRT-PCR validation of RNA-Seq data. The purity (260 nm/280 nm ratio between 1.8 and 2.0 for RNA) and the concentration of each RNA sample were measured according to our previous study [43]. RNA samples were standardized for reverse transcription. cDNA was synthesized using MLV reverse transcriptase (Takara, Japan), according to the manufacturer’s instructions. The β-action gene was used as an internal housekeeping gene. A total of 10 genes from RNA-Seq results were randomly selected as target genes for qRT-PCR analysis (Fig. S2 and S3). Primers for these genes were designed using Primer 5.0 software (Table S12). qRT-PCR cycling conditions were as follows: 94 °C for 2 min, 40 cycles, followed by 94 °C for 15 s, 60 °C for 30 s, and 72 °C for 30 s. The specificity of the PCR products was verified by melt curve analysis for each sample. For each gene, three biological replicates (with five technical replicates for each biological replicate) were performed. Control and target genes for each sample were run on the same plate to control for interpolate variation. The Ct value for each biological replicate was obtained by calculating the mean of three technical replicates. The relative expression level among DC, QC and WC larvae or newly emerged drones was calculated using the 2−ΔΔCt formula reported by Liu and Saint [44].
Data analysis
Data of morphological indexes and reproductive tissues among three groups (Figs. 1 and 2) were analyzed with one-way ANOVA (Statview 5.01 package, SAS Institute Inc., USA) followed by Fisher’s PLSD test. The relative expression levels of 10 genes in qRT-PCR experiment were calculated using 2−ΔΔCt format, and then compared with RNA-Seq results.