Integrative co-expression and deep learning reveal stage-specific regulators of Apis cerana drone embryogenesis

Integrative co-expression and deep learning reveal stage-specific regulators of Apis cerana drone embryogenesis | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Article Integrative co-expression and deep learning reveal stage-specific regulators of Apis cerana drone embryogenesis Xiang Ding, Dan Yue, Jiamin Huan, Runlang Su, Yehui Tan This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-7939377/v1 This work is licensed under a CC BY 4.0 License Status: Posted Version 1 posted You are reading this latest preprint version Abstract Stage-resolved regulators in Apis cerana drone embryos remain unclear. Here, RNA-seq sampled at 12/24/36/48/60 h was analyzed by WGCNA to build co-expression modules, refined with k-nearest neighbors on five-point trajectories, and benchmarked across 19 ML/DL classifiers (including ResNet,SVM,ExtraTrees,KNN and RandomForest). To balance accuracy and interpretability, we combined model benchmarking, consensus voting, and feature attribution (e.g., gradient/SHAP) and report parameter and sampling sensitivities. Expression co-regulated into 6 clusters aligned with developmental stages—germline/epigenetic control (12 h), ncRNA/rRNA-centric regulation (24 h), neurogenesis/morphogenesis (36–48 h), and organ specialization (60 h). ResNet achieved ~ 97% accuracy for gene-pattern classification; consensus reduced > 5,800 candidates to ~ 2,000 high-confidence stage-specific genes and hub PPI subnetworks. The WGCNA→kNN(unsupervised cluster analysis) →ML/DL pipeline thus resolves stage-specific regulators while addressing the “black-box” trade-off, and its performance generalizes with attention to module granularity (β, minModuleSize), k, model capacity/regularization, and temporal sampling density. Biological sciences/Computational biology and bioinformatics Biological sciences/Developmental biology Apis cerana drones embryonic deep learning ResNet Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 Figure 7 Figure 8 Figure 9 1 Introduction The Asian honeybee ( Apis cerana ) is an important social insect that not only provides us with honey and beeswax, but also plays a crucial role as a pollinator in ecosystems [ 1 ] . Embryonic development is the starting phase of honeybee lifecycle, during which the main organs of adult bees progressively develop (REF). So studying its developmental mechanisms is essential for understanding honeybee growth and reproduction [ 2 ] .Studies on the embryo of honeybee workers ( Apis mellifera ) have revealed that it undergoes 10 successive developmental phases, lasting 72 h [ 3 , 4 ] . And significant morphological changes occur and tissues and organs progressively develop from stage 5 to stage 9 [ 5 ] .The strongly expressed proteins and pathways related to transcriptional − translational machinery and morphogenesis at 24 h drone embryo relative to the worker, illustrating the earlier occurrence of morphogenesis in the drone than worker [ 6 ] .Proteomic analysis showed that in the 24-hour stage, proteins and pathways related to transcription-translation mechanisms and morphogenesis show significant expression, indicating that molecular preparation for morphogenesis likely begins before 24 hours [ 6 ] . Moreover, phosphorylated proteins associated with morphogenesis represent the highest proportion [ 7 ] . The 48-hour stage is the most active period of protein expression during embryonic development, primarily involving proteins related to carbohydrate metabolism and energy production [ 8 ] The haploid drones, which develop from unfertilized eggs (REF), are differ from the fertilized diploid the fertilized eggs giving rise to diploid workers are an ideal model for studying the morphological and biochemical characteristics of recessive mutations [ 9 ] .In previous transcriptomic studies of worker bee embryos, 24, 48, and 72 hours were typically chosen as sampling points, roughly corresponding to the early, middle, and late stages of the embryo's transition into a larva [ 10 ] . Based on these studies, we added three more time points: 12h, 36h, and 60h, aiming for a more detailed transcriptomic analysis of drone embryos. The 12h and 60h time points were selected from the longest developmental periods in worker bee development. The 12h stage is crucial for studying early gene expression in drone embryos, during which significant changes in gene and microRNA expression occur [ 11 ] , while 36h represents a key middle point during major morphological changes in worker embryos. To date, there has been little transcriptomic research focusing on these three time points. In this study, we conducted an in-depth analysis of gene expression in the drone embryos of Apis cerana at five key developmental stages (12h, 24h, 36h, 48h, 60h). Through RNA-seq analysis, we observed significant differences in gene expression across these stages, with the most notable changes occurring between 12h and 24h. Additionally, we employed temporal transcriptome analysis and weighted gene co-expression network analysis (WGCNA) to identify stage-specific genes at these five time points, providing deeper insights into their molecular mechanisms. We also discovered potential key genes and identified several important but functionally unknown genes in Apis cerana , which may be closely related to drone embryo development.The development of high-throughput sequencing technologies has made it possible to study the embryonic development of non-model insects at the genomic level, providing opportunities to identify new candidate genes [ 12 ] . While substantial genomic and transcriptomic data are available for the embryonic development of various insect species, including D. melanogaster [ 13 ] , G. bimaculatus [ 14 ] , Athetis lepigone [ 15 ] , T. castaneum [ 16 ] , Oncopeltus fasciatus [ 17 ] , Nilaparvata lugens [ 18 ] , and the silkworm (Bombyx mori) [ 19 ] , transcriptomic studies on the embryonic development of male honeybees remain relatively sparse. It was not until 2020 that the first high-quality genome assembly of Apis cerana was completed [ 20 ] , which laid the foundation for precise gene annotation. Existing research has shown that, within 24 hours, the expression of proteins and pathways related to transcription-translation mechanisms and morphogenesis significantly increases, suggesting that the molecular preparations for morphogenesis begin before the 24-hour mark. The highest proportion of phosphoproteins associated with morphogenesis was observed in the 24-hour phosphoproteome [ 21 ] . At 48 hours, embryonic protein expression peaks, primarily involving proteins related to carbohydrate metabolism and energy production. At 72 hours, proteins related to antioxidation and cytoskeletal functions are significantly upregulated [ 22 ] . Although these studies reveal certain characteristics of protein expression during male honeybee embryogenesis, systematic studies on the dynamic changes in gene expression remain insufficient. To address this gap, we performed a time-series transcriptomic analysis of A. cerana drone embryos at five key stages (12, 24, 36, 48, and 60 h) to delineate stage-specific regulatory frameworks at the transcriptional level. Methodologically, we first constructed co-expression modules using weighted gene co-expression network analysis (WGCNA), then refined within-module expression patterns with k-nearest neighbors (KNN) based on the five-timepoint trajectories. We next combined “module features + KNN subcluster labels/consistency” as inputs to benchmark 19 machine-/deep-learning models. This integrated strategy follows the recently advocated paradigm of “network-first → model benchmarking → interpretable validation”: WGCNA compresses high-dimensional expression into modules to denoise, after which models such as LASSO enable more precise screening of phenotype- or disease-relevant target genes, substantially reducing the search space [ 23 , 24 ] .In head-to-head comparisons, the residual network (ResNet) achieved the highest accuracy for gene-pattern classification (96.78%), with concordant training/validation curves and no signs of overfitting. To further enhance confidence, we applied a consensus-intersection of the top three performers (e.g., ResNet, ExtraTrees, KNN classifier) and integrated this with KNN subcluster consistency, thereby systematically narrowing the initial large candidate pool to > 2,000 high-confidence, stage-specific key genes. We then constructed protein–protein interaction (PPI) networks in Cytoscape, highlighting subnetworks with ≥ 20 nodes, and performed GO annotation, which resolved functional groups related to visual organs, wing discs, gonads, and metabolism. To experimentally support the computational predictions, we randomly selected 1–2 highly expressed representative genes at each stage for qPCR validation; results were consistent with both model predictions and transcriptomic trends. In sum, we establish and implement a closed-loop time-series framework—WGCNA → KNN(unsupervised cluster analysis) → multi-model ML/DL → consensus intersection → PPI/GO → qPCR—that maps stage-resolved gene networks underlying A.cerana drone embryogenesis, from early RNA/epigenetic regulation through mid-stage organ morphogenesis to late-stage tissue specialization and metamorphic preparation. This framework couples accuracy with interpretability and provides an actionable shortlist of testable regulators and pathways, offering a practical template for subsequent functional genetics and breeding applications. 2 Materials and Methods 2.2 Sample collection Embryo samples of Apis cerana were obtained from three queenright colonies in Cangyuan, Yunan Province, China. To collect right-aged embryonic samples, mated egg-laying queens were limited in two frames with empty drone comb, newly laid eggs were monitored at six-hour intervals and were marked with different colored paint pens, the time of the first observed egg deposition was recorded as 0 h. After marking, the comb frame was returned to the hive until eggs from the 12 h, 24 h, 36 h, 48 h, and 60 h stages were obtained. At each stage, 20–30 eggs were collected, lysed in lysis buffer (PRL), and immediately stored at -80°C. In total, samples related to five developmental stages with three biological replicates each were collected for RNA-seq. 2.3 RNA extraction and cDNA library construction Total RNA was extracted from the embryo samples according to the protocol of GeneBetter’s R516 FlashPure Total RNA micro Kit, and then were preserved in lysis buffer and stored at -80°C. The purity and concentration of RNA were measured using NanoDrop 2000, and quantity of the total RNA samples were determined by using Agilent 2100/4200 system. cDNA libraries were prepared according to standard protocols and sequenced as 150-bp paired-end reads on the Illumina HiSeq platform by Anhui Gaohe Biotechnology Co., Ltd. 2.4 Mapping and quality control of RNA-seq data Before mapping, raw data were filtered by the following steps: (1) reads with adapters were removed; (2) reads in which unknown bases represented more than 10% of the total sequence were removed; (3) low quality reads, in which the percentage of low-quality bases (quality no more than 10) exceeded 50%, were removed. The remaining reads were defined as “clean reads” and used for subsequent bioinformatics analyses [ 25 ] . Quality control on the clean data was performed by the production of a base composition chart and quality distribution chart using FastQC. We used Hisat2(v2.2.1) [ 26 ] to map clean reads to the reference genome (GCA_029169275.1). The SAM files were converted to BAM files and sorted using samtools (v1.6) [ 27 ] . These sorted BAM files were then quantified by passing them to Feature Counts [ 28 ] to count the number of reads mapped to each gene. Finally, differential expression analysis was conducted using the DESeq2 algorithm implemented in Trinity (v2.15.1 ) software [ 29 ] . 2.5 PCA analysis To evaluate the consistency of sample collection and investigate the transcriptomic relationships between 12h, 24h, 36h, 48h, and 60h embryos, we performed hierarchical clustering analysis and principal component analysis (PCA) on the full gene expression data of these samples. Both methods effectively captured the similarities and differences between samples. Before conducting these analyses, gene expression levels were standardized using the regularized log transformation method in the limma package. 2.6 GO and KEGG analysis The longest cds protein sequence from the reference Apis cerana genome was annotated for GO/KEGG using emapper.py (2.1.10) with the eggnog database. The annotated file was converted into org.db format, and GO and KEGG enrichment analysis was performed using the clusterProfiler package (v 4.4) in R. 2.7 Temporal clustering analysis of genes Using the ClusterGVis package in R, time-series transcriptome analysis was performed on significantly expressed genes ( |log₂FoldChange| >1 and padj < 0.05) from the expression matrix. Through clustering genes with similar expression patterns, the gene modules of different clusters were analyzed to determine their expression profiles across the various embryonic stages. 2.8 WGCNA, KNN refinement, and downstream deep-learning–based prioritization Raw counts were normalized, lowly expressed genes were filtered, and outliers were removed based on inter-sample correlations. A weighted co-expression network was constructed using WGCNA (R package WGCNA v1.xx): the soft-thresholding power (β) was chosen with pickSoftThreshold to approximate scale-free topology (target R² ≥ 0.8), a signed adjacency matrix was computed and transformed to the topological overlap matrix (TOM), modules were detected by dynamic tree cutting, and module eigengenes were correlated with the five developmental stages (12、24、36、48and60 h). After module detection, we first refined gene-level temporal patterns with KNN clustering within each module, using five-stage trajectories as features (K selected by silhouette/elbow criteria). Each gene received a KNN subcluster label and a neighborhood-consistency score. We then formed a composite feature set (module-level normalized profiles + KNN labels/scores) and fed it into a multi-model machine-/deep-learning pipeline for stage classification and gene prioritization. In total, 19 models were benchmarked (11 deep-learning models: RNN, LSTM, BiLSTM, GRU, CNN, ResNet, Transformer, AttentionRNN, TCN, WaveNet, InceptionTime; and 8 traditional models: e.g., Random Forest, Gradient Boosting, AdaBoost). Data were split 60/20/20 into training/validation/test sets and evaluated with 5-fold cross-validation. The best-performing ResNet used a two-block residual architecture (fully connected 128–64), batch normalization, dropout = 0.3, and was trained with Adam (lr = 0.001) and cross-entropy loss for 100 epochs. For high-confidence calls, we applied a consensus intersection of the top three models (e.g., ResNet, ExtraTrees, KNN classifier) and integrated this with KNN-subcluster consistency, yielding stage-specific key-gene sets across all five time points. Subsequent PPI/GO/qPCR analyses and visualizations are described in dedicated sections. 3 Result 3.1 Overall RNA-Seq results and quality control Illumina RNA sequencing generated an average of 66.52 million (M) raw reads per sample. After filtering, an average of 66.20 million clean reads per sample were retained for subsequent transcriptome analysis. The clean reads had high quality, with average Q20 and Q30 scores of 98.84% and 96.31%, respectively. Additionally, the average GC content was 39.49% (Supplemental Table S1 ). HISAT2 software was used for alignment against the reference genome. Sliding window density analysis showed that more than 80% of transcripts were uniquely aligned to the reference genome. Across 15 samples, a total of 61.82 million clean reads were uniquely mapped to the reference genome, with an average mapping rate of 89.86%. Only 10.14% of reads (7.97 million) failed to align with the genome (Supplemental Table S2 ). 3.2 Hierarchical clustering and principal component analysis of samples The results of hierarchical clustering and PCA (Fig. 1 a, b) demonstrated that the samples collected at the same time point had good homogeneity. The genes related to 12h were significantly different from the other four stages (Fig. 1 b), which exhibited the highest expression levels (Fig. 1 a). It indicated that gene expression profiles at 12h were markedly distinct from the other four stages. 3.3 Analysis of differentially expressed genes (DEGs) To study the differential expression of genes during the development of drone embryos, we conducted a statistical analysis of differentially expressed genes (DEGs) across five developmental stages of drone embryos, comparing gene expression differences between 12h and 24h, 24h and 36h, 36h and 48h, and 48h and 60h. The results (Fig. 2 a) showed the following: between 12h and 24h, 2538 genes were significantly upregulated, and 3282 genes were significantly downregulated; between 24h and 36h, 1004 genes were significantly upregulated, and 1723 genes were significantly downregulated; between 36h and 48h, 888 genes were significantly upregulated, and 574 genes were significantly downregulated; and between 48h and 60h, 328 genes were significantly upregulated, and 1268 genes were significantly downregulated. The most differentially expressed genes (both upregulated and downregulated) were observed between the 12h and 24h developmental stages (Fig. 2 b), while the degree of gene expression differences between 36h, 48h, and 60h gradually decreased. The heatmap of differentially expressed genes (Fig. 2 c) indicates clear stage-specific gene expression across the five developmental stages. These stage-specific genes are likely associated with physiological, biochemical processes, morphological changes, and behavioral activities, potentially playing an important role in the development of drone embryos. We selected the top 20 genes with the most significant expression differences at each stage, filtering out 21 non-coding genes ( Fig. 2 b ) . Through EggNOG database annotation and protein function prediction, we found that these genes are primarily involved in functions such as DNA-binding transcription factors, regulation of gene expression, G-protein-coupled receptors, cytoskeleton binding, methylation, negative regulation of programmed cell death, protein transport, phosphorylation, and molecular interactions during embryonic development. Additionally, we identified five genes ( LOC100578606, LOC100578389, LOC100577882, LOC551929, LOC725838 ) whose functions in Apis cerana are still unknown. These genes may play a critical role in drone embryo development, but further experiments are needed to verify their specific functions. 3.4 Gene expression in drone embryos is co-regulated across 10 clusters, which correspond to different developmental stages. To further explore the dynamic changes in gene expression during drone embryo development, we applied the k-means clustering algorithm, dividing 9492 genes from the reference genome into 10 clusters based on the similarity of their expression patterns. Detailed analysis of these 10 clusters revealed significant gene enrichment at specific developmental stages (Fig. 3 ). The C1 and C2 modules were highly expressed significantly enriched at the 12-hour time point of drone embryo development, with related GO terms mainly involving sperm cell differentiation, maintenance of germline stem cell populations, male gametogenesis, and functions associated with protein modifications such as phosphorylation, methylation, epigenetic regulation of gene expression, and mRNA regulation or degradation. On the other hand, genes in the C5 and C8 modules showed higher expression after 24 hours of embryo development, primarily related to processes such as cell growth and differentiation, transcription regulation, and rRNA and ncRNA metabolism. 3.5 Co-Expression Network Construction and Hub Module Identification Weighted gene co-expression network analysis (WGCNA) were conducted to explore the association of genes at five specific developmental stages using the all aligned genes (9492). The samples were clustered (Fig. 4 a), missing values were processed, and outliers were excluded. To ensure that the network exhibited scale-free characteristics, a soft-threshold power of 24 was selected (R² for the scale-free network = 0.80) (Fig. 4 b,c). A co-expression matrix was constructed using the one-step method, and dynamic tree cutting identified nine gene modules (Fig. 4 a,d), which demonstrated their correlation with drone embryo development. The results showed that the turquoise module (containing 2643 genes) had the highest correlation with the 12-hour embryo development stage (cor = 0.69; P = 0.004) (Fig. 4 d); the red and blue modules had the highest correlation with the 24-hour stage; the black module was most correlated with the 36-hour stage (cor = 0.69; P = 0.004); the brown and grey modules were associated with the 48-hour stage (cor = 0.69; P = 0.004); and the yellow, pink, and green modules were highly correlated with the 60-hour stage (cor = 0.69; P = 0.004) (Fig. 4 d). The differences in the correlation of these nine modules across the developmental stages of drone embryos indicate that gene expression at each stage exhibits clear specificity. 3.6 Integrative WGCNA and Deep Learning Identify Stage-Specific Genes in Apis cerana Embryogenesis In this study, we first conducted a systematic analysis of gene co-expression modules underlying Apis cerana embryonic development using WGCNA. Although only five developmental time points were sampled, heterogeneity in temporal expression dynamics enabled further refinement: by combining K-means clustering with module–trait association analysis, we partitioned the data into six representative expression modules (Fig. 5 ). Building on these results, we applied machine-learning models to prioritize key genes, constructing and evaluating 19 models in total—11 deep learning models (RNN, LSTM, BiLSTM, GRU, CNN, ResNet, Transformer, AttentionRNN, TCN, WaveNet, InceptionTime) and 8 traditional models (e.g., Random Forest, Gradient Boosting, AdaBoost). The ResNet model achieved the best performance for gene-expression pattern classification, reaching an accuracy of 96.78% (Fig. 6 –8,Supplemental Fig. 1). Notably, training and validation learning curves were highly concordant, indicating no overfitting and demonstrating strong generalization and robustness. Among all candidates, we further selected the top three performers and intersected their predictions, yielding a robust set of > 2,000 genes with putative stage-specific roles, which were subsequently categorized for downstream functional analyses (Fig. 6 –8,Supplemental Fig. 1). (a) The number of correctly predicted genes by ResNet, ExtraTrees, and KNN. ResNet achieved the highest number of correctly predicted genes (2251), followed by ExtraTrees (2186) and KNN (2148).(b) Gene overlap matrix showing the number of shared correctly predicted genes between models. Diagonal values represent the total correctly predicted genes for each model, while off-diagonal values represent the overlaps between two models (e.g., ResNet–ExtraTrees: 2156; ResNet–KNN: 2118; ExtraTrees–KNN: 2098). The color scale indicates the number of shared genes, with darker colors representing larger overlaps. 3.7 PPI network construction and hub genes identification Building on the WGCNA-derived co-expression modules, we applied KNN clustering to further refine temporal expression subgroups within modules and, based on these refined clusters, selected candidate genes for downstream analyses. We then constructed protein–protein interaction (PPI) networks in Cytoscape and visualized subnetworks with ≥ 20 nodes (Fig. 8). Using network topology metrics, stage-specific hub genes were prioritized; these hubs aggregated into functional modules that collectively drive key physiological and biochemical processes during honeybee embryogenesis (Fig. 8). In addition, 1–2 highly expressed representative genes per stage were randomly selected for qPCR validation, and the results were concordant with both model predictions and transcriptomic trends, supporting the reliability and reproducibility of the identified hubs and modules. Finally, Gene Ontology (GO) annotation of these genes delineated their stage-specific functional roles (Fig. 9 ; Supplemental Tables 3 and 5). Fugure 8 . PPI networks across five embryonic stages built from ResNet-prioritized genes 4 Discussion At the 12h early developmental stage, temporal transcriptomic dynamics were closely associated with cytoskeletal remodeling, precise chromosomal positioning, and the formation of male germ cells (Fig. 8a). In parallel, epigenetic and phosphorylation-associated layers likely contribute to stage-specific control, consistent with systems-level evidence implicating TOR/epigenetic axes and phospho-regulatory drivers in developmental regulation [ 30 , 31 ] . Concordant results from WGCNA, deep-learning–guided analyses, and temporal clustering (Supplemental Tables 3a, 3b; Table 4a) further support the importance of these early programs. Collectively, these observations indicate that multiple regulatory layers act in concert to stabilize and specify developmental trajectories at 12 h. In 24h embryos (Fig. 3 ,Supplemental Table 4b), GO enrichment of Cluster 5 and Cluster 8 highlighted potentially critical roles for ncRNAs and rRNAs. Foundational studies have established that ncRNAs—including microRNAs and lncRNAs—govern chromatin states, transcriptional outputs, RNA processing, and translation, thereby shaping cell division, differentiation, and developmental timing [ 32 – 35 ] . Moreover, gene-set enrichment comparing 12 h with 24 h module genes showed higher 24 h expression in sets related to compound-eye photoreceptor differentiation, ncRNA processing, and rRNA processing (Fig. 8a,b;Supplemental Table 5a,b), reinforcing stage-specific deployment of RNA-centric regulation. At 36h, hub-gene GO enrichment indicated DNA binding, neuronal projection development, cell morphogenesis, and differentiation (Fig. 8c; Supplemental Tables 4c and 5c), consistent with accelerated organogenesis and functional specialization. Nutrient- and growth-sensing pathways (e.g., TOR) intersect with transcriptional/epigenetic circuits, providing a plausible framework for coupling metabolic state to morphogenesis at this stage [ 30 ] . GSEA comparing 24 hand 36 h further showed that genes related to male gamete generation and RNA splicing were significantly upregulated at 24 h (Supplemental Fig. 2). Given that alternative splicing (AS) broadly expands transcriptome/proteome diversity and profoundly impacts developmental phenotypes—with recent studies underscoring rapid, organism-wide AS remodeling [ 34 , 36 – 38 ] 。We infer that heightened AS activity at 24 h promotes expression control of sex-differentiation and developmental genes, with important consequences for drone physiology and biochemistry. At 48h, hub-gene GO terms centered on anterior–posterior patterning in wing discs, adult walking behavior, gonad development, eye morphogenesis, sensory organ morphogenesis, and positive regulation of the Wnt pathway (Fig. 8d; Supplemental Tables 4d and 5d). This constellation is consistent with coordinated deployment of core developmental signaling during organogenesis in holometabolous insects, supporting polarity establishment, proliferation, and differentiation [ 39 , 40 ] . At 60h, temporal clustering (Cluster 9) and hub-gene enrichment jointly indicated robust organ morphogenesis and skeletal muscle development (Supplemental Table 4e; Table 5e). Homology analysis mapped LOC410326 to Drosophila cbt ; given the known roles of cbt in cell-cycle control and tissue growth, the bee homolog likely contributes to stem-cell maintenance, rapid tissue proliferation, and tissue specializationorgan development, thereby laying the groundwork for subsequent metamorphosis. Thus, late embryogenesis is characterized by intensified, directionally patterned growth consistent with imminent metamorphic transitions [ 41 , 42 ] . We further dissected the deep-learning outputs with a composite scoring system to identify genes with pronounced temporal dynamics. For dynamic genes, Dynamic Score = R² × 20 + |slope| × 15 + log₂(fold_change + 1) × 5; for stable genes, Stability Score = |slope| × 10 + CV × 10 + Δmax × 2 (lower scores indicate greater stability). Genes were stratified into four tiers: Tier 1 (strongly dynamic)—Dynamic Score > 30 and R² >0.8 with monotonic trends (|slope| >1.0) and substantial changes (fold_change > 4); Tier 2 (transient response)—Peak/Valley patterns with fold_change > 2; Tier 3 (complex regulation)—Complex patterns or R² < 0.5; Tier 4 (homeostatic)—Stability Score < 10 with CV < 0.2 and |slope| < 0.1 (Supplemental Figure. 3). To validate stratification, we randomly selected the most pronounced Tier-1 genes and the lowest-scoring Tier-4 genes for box-plot visualization and qPCR (Supplemental Fig. 4a,b, 5a,b; Supplemental Tables 6, 7). Tier-1 genes exhibited consistent linear trends across time (mean R² = 0.85 ± 0.08), whereas Tier-4 genes remained stably expressed (mean CV = 0.15 ± 0.05). This framework strengthens reliability assessment via R², enhances biological interpretability via pattern label (Increasing/Decreasing/Peak/Valley/Stable/Complex), and optimizes identification of both dynamic and stable genes via dual scoring, thereby improving functional categorization and providing robust annotation for downstream model validation. 5 Conclusion Altogether, our results provide new insights into the molecular underpinnings of A. cerana drone embryo genesis. The discovery of previously uncharacterized but network-central genes points to additional developmental regulators. More broadly, we delineate gene networks mediating the transition from early germ-line/organ primordia to mature organ systems. By integrating WGCNA, deep learning (ResNet), and temporal transcriptomics, we defined stage-specific key genes, modules, and hubs, highlighting histone-linked/epigenetic regulation, ncRNAs, and sex-differentiation programs in drone embryogenesis and offering a theoretical foundation for bee developmental biology, breeding, and conservation. Declarations 6 Funding This research is supported by the Yiwu Industrial and Commercial College Program under Grant Nos. XJKJ2502YB. This work was supported by Jinhua Key Laboratory of Robot Intelligent Welding Technology and by Jinhua Public Welfare Technology Application Research Project in 2024. 7 Availability of data and materials The short-read RNA-Seq data involved in this study were submitted to the NCBI BioProject database (https://www.ncbi.nlm.nih.gov/bioproject/ accessed on 22 January 2025), with accession numbers PRJNA1186594. 8 Authors ’ contributions The genome assembly and analyses were conducted by Xiang Ding and Yehui Tan. The manuscript was written by Runlang Su and Dan Yue. 9 Ethics approval and consent to participate Not Applicable. 10 Consent for publication Not Applicable. 11 Competing interests The authors declare that they have no competing interests. References Pan, L-X., Hu, W-W., Cheng, F-P., Hu, X-F. & Wang, Z-L. Transcriptome analysis reveals differentially expressed genes between the ovary and testis of the honey bee Apis mellifera. Apidologie 53 , 15 (2022). Fang, Y. et al. In-depth Proteomics Characterization of Embryogenesis of the Honey Bee Worker. Mol. Cell. Proteom. 13 , 2306–2320 (2014). DuPraw, E. J. The honeybee embryo, (1967). Hu, X., Ke, L. & Wang, Z. Zeng Z Dynamic transcriptome landscape of Asian domestic honeybee ( Apis cerana ) embryonic development revealed by high-quality RNA sequencing. BMC Dev. Biol. 18 , 11 (2018). Cridge, A. G. et al. Dearden PK The honeybee as a model insect for developmental genetics. genesis 55 , e23019 (2017). Fang, Y. et al. Analysis Unravels Mechanism Underling the Embryogenesis of the Honeybee Drone and Its Divergence with the Worker ( Apis mellifera lingustica ). J. Proteome Res. 14 , 4059–4071 (2015). Ma, B., Ma, C., Li, J. & Fang, Y. Revealing phosphorylation regulatory networks during embryogenesis of honey bee worker and drone ( Apis mellifera ). Front. Cell. Dev. Biol. 10 , 1006964 (2022). Tan, K. et al. & Gudrun Geographic distribution of the eastern honeybee, Apis cerana (Hymenoptera: Apidae), across ecological zones in China: morphological and molecular analyses. Syst. Biodivers. : 473–482, (2006). Beye, M. The dice of fate: the csd gene and how its allelic composition regulates sexual development in the honey bee. Apis mellifera BioEssays . 26 , 1131–1139 (2004). Hu, X. F. & Ke, L. Zeng ZJ Histomorphological study on embryogenesis of the honeybee Apis cerana . J. Asia. Pac. Entomol. 22 , 860–867 (2019). C, M. W. D., Z, H. & Z, L. & L W Sex-Specific Development in Haplodiploid Honeybee Is Controlled by the Female-Embryo-Specific Activation of Thousands of Intronic LncRNAs. Front. (2021). Cell Dev. Biol,(. Li, J., Fang, Y. & Zhang, L. Begna D Honeybee ( Apis mellifera ligustica ) drone embryo proteomes. J. Insect. Physiol. 57 , 372–384 (2011). Oppenheim, S. J., Baker, R. H. & Simon, S. DeSalle R We can't all be supermodels: the value of comparative transcriptomics to the study of non-model insects. Insect Mol. Biol. 24 , 139–154 (2015). Pires, C. V., Freitas, F. C., Cristino, A. S. & Dearden, P. K. Simoes ZL Transcriptome Analysis of Honeybee ( Apis mellifera ) Haploid and Diploid Embryos Reveals Early Zygotic Transcription during Cleavage. PLoS One . 11 , e0146447 (2016). Simon, S. et al. DeSalle R Comparative transcriptomics reveal developmental turning points during embryogenesis of a hemimetabolous insect, the damselfly Ischnura elegans. Sci. Rep. 7 , 13547 (2017). Sandler, J. E. & Stathopoulos, A. Quantitative Single-Embryo Profile of Drosophila Genome Activation and the Dorsal–Ventral Patterning Network. Genetics 202 (4), 1575–1584 (2016). Bando, T. et al. Mito T Analysis of RNA-Seq data reveals involvement of JAK/STAT signalling during leg regeneration in the cricket Gryllus bimaculatus . Development 140 , 959–964 (2013). Li, L-T., Zhu, Y-B., Ma, J-F. & Li, Z-Y. Dong Z-P An Analysis of the Athetis lepigone Transcriptome from Four Developmental Stages. PLoS One . 8 , e73911 (2013). Pridöhl, F. et al. Transcriptome sequencing reveals maelstrom as a novel target gene of the terminal system in the red flour beetle Tribolium castaneum. Development 144 (7), 1339–1349 (2017). Ginzburg, N. & Cohen, M. & Chipman AD Factors involved in early polarization of the anterior-posterior axis in the milkweed bug Oncopeltus fasciatus . genesis 55, (2017). Fan, X-B. et al. An Overview of Embryogenesis: External Morphology and Transcriptome Profiling in the Hemipteran Insect Nilaparvata lugens. Front. Physiol. 11 , (2020). Wang, Z. L. et al. A Chromosome-Scale Assembly of the Asian Honeybee Apis cerana Genome. Front. Genet. 11 , 279 (2020). Wang, Y. et al. Integrative machine learning models predict prostate cancer diagnosis and biochemical recurrence risk: Advancing precision oncology. npj Digit. Med. 8 , 524 (2025). Zhu, Z. et al. 2025 Exploration of differential expression and biological significance of amino acid metabolism genes in osteoarthritis. Front. Immunol. Volume 16 , (2025). de Sena Brandine, G. & Smith, A. D. Falco: high-speed FastQC emulation for quality control of sequencing data. F1000Res Nov . 7 , 8:1874 (2019). Kim, D., Langmead, B. & Salzberg, S. HISAT: a fast spliced aligner with low memory requirements. Nat Methods (2015). Li, H. et al. Genome Project Data Processing S: The Sequence Alignment/Map format and SAMtools. Bioinformatics 25 (16), 2078–2079 (2009). Liao, Y., Smyth, G. K. & Shi, W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 30 (7), 923–930 (2014). Haas, B. J., Papanicolaou, A. & Yassour, M. A De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. Nat. Protoc. 8 , 1494–1512 (2013). Sun, J. et al. TOR signaling, epigenetics, and tumor development. Front. Genet. 13 , (2022). Zhang, T. et al. Drivers of protein phosphorylation (multi-omics). Genome Biol. 24 (1), 52 (2023). Bartel, D. P. MicroRNAs:genomics,biogenesis, mechanism, and function. Cell 116 (2), 281–297 (2004). Malone, C. D. Hannon GJ Small RNAs as Guardians of the Genome. Cell 136 (4), 656–668 (2009). Fu, X-D. & Non-coding RNA: a new frontier. Natl. Sci. Rev. 1 , 190–204 (2014). Qu, Z. & Adelson, D. L. Conservation and roles of ncRNA. Front. Genet. 3 . (2012). Calixto, C. P. G. et al. Rapid, dynamic AS in Arabidopsis cold response. Plant. Cell. 30 , 1424–1444 (2018). Zhang, H. et al. Population-level AS and agronomic traits of rice. Plant Cel.(2024) l. Zhong, Y. et al. Plant.Pan-transcriptomic AS and cold tolerance in rice. (2024). Nusse, R. & Clevers, H. Wnt/β-Catenin Signaling, Disease, and Emerging Therapeutic Modalities. Cell 169 (6), 985–999 (2017). Tripathi, B. K. & Irvine, K. D. The wing imaginal disc. Genetics. Apr 4;220(4):iyac020. PMID: 35243513; PMCID: PMC8982031. (2022). Zhang, P. & Edgar, B. A. The Krüppel-like factor Cabut has cell cycle regulatory properties similar to E2F1.118,( (2025). Blanco, E. et al. Gene expression following induction of regeneration in Drosophila wing imaginal discs. Expression profile of regenerating wing discs(2010). Additional Declarations No competing interests reported. Supplementary Files SupplementalFigure1.pdf SupplementalFigure2.pdf SupplementalFigure3.pdf SupplementalFigure4.pdf SupplementalFigure5.pdf SupplementalTable1.xlsx SupplementalTable2.xlsx SupplementalTable3.xlsx SupplementalTable4a.csv SupplementalTable4b.csv SupplementalTable4c.csv SupplementalTable4d.csv SupplementalTable4e.csv SupplementalTable5a.csv SupplementalTable5b.csv SupplementalTable5c.csv SupplementalTable5e.csv SupplementalTable6.xlsx SupplementalTable7.xlsx Cite Share Download PDF Status: Posted Version 1 posted You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-7939377","acceptedTermsAndConditions":true,"allowDirectSubmit":true,"archivedVersions":[],"articleType":"Article","associatedPublications":[],"authors":[{"id":550059202,"identity":"ef6dbacd-cfc8-433a-bdd0-c597c0df0da8","order_by":0,"name":"Xiang Ding","email":"","orcid":"","institution":"University of Chinese Academy of Sciences","correspondingAuthor":false,"prefix":"","firstName":"Xiang","middleName":"","lastName":"Ding","suffix":""},{"id":550059203,"identity":"a5e0da03-fd2f-473b-9d6a-18cc2fa1cc9d","order_by":1,"name":"Dan Yue","email":"","orcid":"","institution":"College of Animal Science and Technology, Yunnan Agricultural University","correspondingAuthor":false,"prefix":"","firstName":"Dan","middleName":"","lastName":"Yue","suffix":""},{"id":550059204,"identity":"2830945b-7d5e-4b72-891a-9f6f581370d1","order_by":2,"name":"Jiamin Huan","email":"","orcid":"","institution":"University of Chinese Academy of Sciences","correspondingAuthor":false,"prefix":"","firstName":"Jiamin","middleName":"","lastName":"Huan","suffix":""},{"id":550059205,"identity":"53792d23-001a-41bc-9c0a-5450e6134bcc","order_by":3,"name":"Runlang Su","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA1klEQVRIiWNgGAWjYNCCAiBmbwASBhbEajEAYp4DIIYEKVokEkAsIrTIz8h9Jl1gYJcnH/n86oYfBRIM/O3dCfjNv5FuJj3DILnY8HZO2c0eoMMkzpzdgF+LRBqbNI8Bc+LG2TlpN3iAWgwkcvFrkZ8B1lKfuHHmmbSbf4jRwnADrOVw4nwJ9mO3ibLF4MwzZmseg+OJG3hy2G7LGEjwEPSLfHsa422eiurE+e3Hn91888dGjr+9l4DD4NYd4AHFDgMPccrB1jWwPyBe9SgYBaNgFIwoAACu2UGGvyXSUQAAAABJRU5ErkJggg==","orcid":"","institution":"College of Animal Science and Technology, Yunnan Agricultural University","correspondingAuthor":true,"prefix":"","firstName":"Runlang","middleName":"","lastName":"Su","suffix":""},{"id":550059206,"identity":"eba1b59a-239f-44fa-af52-c3af2aee16e9","order_by":4,"name":"Yehui Tan","email":"","orcid":"","institution":"University of Chinese Academy of Sciences","correspondingAuthor":false,"prefix":"","firstName":"Yehui","middleName":"","lastName":"Tan","suffix":""}],"badges":[],"createdAt":"2025-10-24 14:59:27","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-7939377/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-7939377/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":96711000,"identity":"a37f8473-0555-45d6-8171-f71135f1f7f5","added_by":"auto","created_at":"2025-11-25 10:11:29","extension":"docx","order_by":2,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":6816797,"visible":true,"origin":"","legend":"","description":"","filename":"IntegrativecoexpressionanddeeplearningrevealstagespecificregulatorsofApisceranadroneembryogenesis.docx","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/d9dac1f157dd8cf57d02cec0.docx"},{"id":96710380,"identity":"0a69e86c-5d73-4fc5-9326-d8dbdc8af814","added_by":"auto","created_at":"2025-11-25 10:10:33","extension":"json","order_by":10,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":6236,"visible":true,"origin":"","legend":"","description":"","filename":"9e709631659943ec87068723e1aeaff5.json","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/4c00f0a432b74ee2399ffa16.json"},{"id":96689296,"identity":"ba494605-8338-4b15-8053-0c2b1cb729cb","added_by":"auto","created_at":"2025-11-25 06:09:04","extension":"pdf","order_by":11,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":87372,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementalFigure1.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/d626e5a8da7a81a0e29999b3.pdf"},{"id":96689340,"identity":"20327bbd-895b-4bbb-a537-216c468cdbdf","added_by":"auto","created_at":"2025-11-25 06:09:06","extension":"pdf","order_by":12,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":10503377,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementalFigure2.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/7bc427e5a7122ef2037ece8a.pdf"},{"id":96689319,"identity":"530440c1-e562-4edc-bdbc-9af3cff86dbf","added_by":"auto","created_at":"2025-11-25 06:09:05","extension":"pdf","order_by":13,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":461687,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementalFigure3.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/d653de07c66530b7c5485ddb.pdf"},{"id":96710573,"identity":"cc83a25d-0f84-4ce1-b4bd-e90b01a3417e","added_by":"auto","created_at":"2025-11-25 10:10:56","extension":"pdf","order_by":14,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":1872614,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementalFigure4.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/33ec4a96d7446deb13574f35.pdf"},{"id":96710386,"identity":"a57acf42-0bdd-482b-a1bc-033d2f113ca0","added_by":"auto","created_at":"2025-11-25 10:10:34","extension":"pdf","order_by":15,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":485338,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementalFigure5.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/7fd7ec55a3e9259f7bf2c698.pdf"},{"id":96689309,"identity":"646ae367-35d2-40e2-8b0b-c9086ec20773","added_by":"auto","created_at":"2025-11-25 06:09:05","extension":"xlsx","order_by":16,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":10738,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementalTable1.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/8aaa69807f19ab409f4dd235.xlsx"},{"id":96689332,"identity":"4641fb1e-483b-4ab1-91d8-990ac58e8267","added_by":"auto","created_at":"2025-11-25 06:09:06","extension":"xlsx","order_by":17,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":10309,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementalTable2.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/6a46aa22dc18126916e10ad2.xlsx"},{"id":96689300,"identity":"116a8dcf-3ad6-46bd-84d0-61a4e658ee82","added_by":"auto","created_at":"2025-11-25 06:09:04","extension":"xlsx","order_by":18,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":15774,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementalTable3.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/014d209e71832e2447008610.xlsx"},{"id":96689286,"identity":"74792ef4-a905-4e6e-a1db-7121d0f24644","added_by":"auto","created_at":"2025-11-25 06:09:04","extension":"csv","order_by":19,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":87362,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementalTable4a.csv","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/40eaf0fe1bf15f936a319da7.csv"},{"id":96689310,"identity":"9f903e0f-d5e4-440b-bf35-69cbc1300888","added_by":"auto","created_at":"2025-11-25 06:09:05","extension":"csv","order_by":20,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":17173,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementalTable4b.csv","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/15e851f2be4604754441a98a.csv"},{"id":96710581,"identity":"0e83501f-d7a7-483d-8702-bfc01557c70a","added_by":"auto","created_at":"2025-11-25 10:10:57","extension":"csv","order_by":21,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":12162,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementalTable4c.csv","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/f8d383c482cec64f8f61ebc5.csv"},{"id":96710802,"identity":"80855d71-ecb7-419d-9134-4f347a3f953b","added_by":"auto","created_at":"2025-11-25 10:11:12","extension":"csv","order_by":22,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":26785,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementalTable4d.csv","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/629b4301a8c39eb1fa1361f0.csv"},{"id":96710799,"identity":"f7ea317a-79a8-4bf8-aa6f-a7361a94e47e","added_by":"auto","created_at":"2025-11-25 10:11:11","extension":"csv","order_by":23,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":13281,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementalTable4e.csv","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/cea124fadb4cf76ae9f9948b.csv"},{"id":96689311,"identity":"392f8d67-7391-4f00-a894-ee91abc03cb6","added_by":"auto","created_at":"2025-11-25 06:09:05","extension":"csv","order_by":24,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":26669,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementalTable5a.csv","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/0494f23df6757597b932d8b2.csv"},{"id":96689301,"identity":"3ceb7f13-faf4-448a-abc4-d85a69553be9","added_by":"auto","created_at":"2025-11-25 06:09:05","extension":"csv","order_by":25,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":27157,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementalTable5b.csv","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/55a88c72d270a26b35a261b6.csv"},{"id":96689333,"identity":"8c0ddd52-0909-46ba-b4de-f31d170ca03b","added_by":"auto","created_at":"2025-11-25 06:09:06","extension":"csv","order_by":26,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":9282,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementalTable5c.csv","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/18ff5fbfb785a4f6817d0307.csv"},{"id":96710790,"identity":"f080d137-e74e-4b99-bdf8-e4036d9b93eb","added_by":"auto","created_at":"2025-11-25 10:11:11","extension":"csv","order_by":27,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":47103,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementalTable5e.csv","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/8abbbffa807d3b8aacaa4187.csv"},{"id":96710930,"identity":"80d4a09d-d0f1-400e-9394-20e00dc931d1","added_by":"auto","created_at":"2025-11-25 10:11:26","extension":"xlsx","order_by":28,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":16493,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementalTable6.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/82a8806754d1e0052d14bd46.xlsx"},{"id":96710943,"identity":"16c2c45b-2c48-47cf-85c0-5126a2c46e7e","added_by":"auto","created_at":"2025-11-25 10:11:26","extension":"xlsx","order_by":29,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":13072,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementalTable7.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/deb74ceb97de047f9c4c80f2.xlsx"},{"id":96689326,"identity":"f5b2f498-c392-4df8-9cfe-a69a8cad465e","added_by":"auto","created_at":"2025-11-25 06:09:05","extension":"xml","order_by":30,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":101971,"visible":true,"origin":"","legend":"","description":"","filename":"9e709631659943ec87068723e1aeaff51enriched.xml","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/9bd7ba83e06674bc9f72ca4f.xml"},{"id":96689330,"identity":"24c46e4c-c764-4324-a302-99b3c032cc47","added_by":"auto","created_at":"2025-11-25 06:09:06","extension":"pdf","order_by":31,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":38510,"visible":true,"origin":"","legend":"","description":"","filename":"Fig1.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/d9d36a60c42b53124f69b4ae.pdf"},{"id":96689307,"identity":"f6e58626-724e-4939-b044-3b136b06cb81","added_by":"auto","created_at":"2025-11-25 06:09:05","extension":"pdf","order_by":32,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":1231448,"visible":true,"origin":"","legend":"","description":"","filename":"Fig2.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/ed93d6d76eac47e84361fb41.pdf"},{"id":96711616,"identity":"0f419e9d-2df1-43ee-8b4c-7ebb8f607bcd","added_by":"auto","created_at":"2025-11-25 10:12:20","extension":"pdf","order_by":33,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":1492095,"visible":true,"origin":"","legend":"","description":"","filename":"Fig3.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/85b045c72fdecd942340385d.pdf"},{"id":96710928,"identity":"7ccf2012-1861-4cc4-9f4e-c34f8656069c","added_by":"auto","created_at":"2025-11-25 10:11:25","extension":"pdf","order_by":34,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":61819,"visible":true,"origin":"","legend":"","description":"","filename":"Fig4.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/eb47a150a2a77c4aa422b591.pdf"},{"id":96710455,"identity":"d0e34bf4-ee47-438b-bc1a-1563dc6109ae","added_by":"auto","created_at":"2025-11-25 10:10:41","extension":"pdf","order_by":35,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":33572,"visible":true,"origin":"","legend":"","description":"","filename":"Fig5.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/256c217bd0fc2888d92f9ed2.pdf"},{"id":96710788,"identity":"617d55c4-19d4-4b00-ad38-8b4c8bf83cde","added_by":"auto","created_at":"2025-11-25 10:11:11","extension":"pdf","order_by":36,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":471127,"visible":true,"origin":"","legend":"","description":"","filename":"Fig6.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/508d1fbd1eddb1d6ed8eae21.pdf"},{"id":96710569,"identity":"757b4e15-7be0-4db1-afa3-be5ff904dfeb","added_by":"auto","created_at":"2025-11-25 10:10:56","extension":"pdf","order_by":37,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":765026,"visible":true,"origin":"","legend":"","description":"","filename":"Fig7.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/ae428d7c9f3184f7d7b3f454.pdf"},{"id":96689303,"identity":"5a8f4cab-243e-49d9-80d9-921c9a7f9cb2","added_by":"auto","created_at":"2025-11-25 06:09:05","extension":"pdf","order_by":38,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":416073,"visible":true,"origin":"","legend":"","description":"","filename":"Fig8.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/9da58caf73b6b2ec53d009d3.pdf"},{"id":96711227,"identity":"7c280c94-3f35-405c-8851-2fe18391ac8f","added_by":"auto","created_at":"2025-11-25 10:11:47","extension":"pdf","order_by":39,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":10558,"visible":true,"origin":"","legend":"","description":"","filename":"Fig9.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/ce9f81d03e78679e0ed60776.pdf"},{"id":96710464,"identity":"9d494c1b-9e31-46c2-9be3-07217749d4a1","added_by":"auto","created_at":"2025-11-25 10:10:41","extension":"jpeg","order_by":40,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":1144477,"visible":true,"origin":"","legend":"","description":"","filename":"floatimage1.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/764cf2e566b6b308db1d98eb.jpeg"},{"id":96689322,"identity":"834e8ce3-5334-4bca-91d8-9662a79d4343","added_by":"auto","created_at":"2025-11-25 06:09:05","extension":"jpeg","order_by":41,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":1478794,"visible":true,"origin":"","legend":"","description":"","filename":"floatimage2.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/44315c3c8378eecc831e22cb.jpeg"},{"id":96689346,"identity":"c6d7af48-151c-4604-b69c-b8a2f55be0fd","added_by":"auto","created_at":"2025-11-25 06:09:06","extension":"jpeg","order_by":42,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":4334936,"visible":true,"origin":"","legend":"","description":"","filename":"floatimage3.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/022776a10c6199323b381949.jpeg"},{"id":96689329,"identity":"a09bd9be-2ebe-45f1-810b-80c9837be41f","added_by":"auto","created_at":"2025-11-25 06:09:05","extension":"jpeg","order_by":43,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":1386935,"visible":true,"origin":"","legend":"","description":"","filename":"floatimage4.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/dce848a953fb7e5583e8f670.jpeg"},{"id":96689347,"identity":"46c9fdee-ec20-4ddd-90c5-ab7dcd0bb6b1","added_by":"auto","created_at":"2025-11-25 06:09:06","extension":"jpeg","order_by":44,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":2486454,"visible":true,"origin":"","legend":"","description":"","filename":"floatimage5.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/5d84485103c31f3eb3012854.jpeg"},{"id":96711466,"identity":"8b3281a5-bd71-4b3c-b66d-a0b86a508c6c","added_by":"auto","created_at":"2025-11-25 10:12:02","extension":"jpeg","order_by":45,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":2357068,"visible":true,"origin":"","legend":"","description":"","filename":"floatimage6.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/a40aae686ce06a09cd787333.jpeg"},{"id":96689315,"identity":"7e4cd9d0-ee3e-418e-9ac6-37335cb2d94c","added_by":"auto","created_at":"2025-11-25 06:09:05","extension":"jpeg","order_by":46,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":3625602,"visible":true,"origin":"","legend":"","description":"","filename":"floatimage7.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/77fff32cdd3f5d2bf6b4b980.jpeg"},{"id":96689336,"identity":"16269914-b0da-48ad-8e62-3afaba278dac","added_by":"auto","created_at":"2025-11-25 06:09:06","extension":"jpeg","order_by":47,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":5004411,"visible":true,"origin":"","legend":"","description":"","filename":"floatimage8.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/3fb3c75df057a689a8e8f36b.jpeg"},{"id":96711612,"identity":"86271899-7b85-4d3d-baa5-bbd6ada93b15","added_by":"auto","created_at":"2025-11-25 10:12:19","extension":"jpeg","order_by":48,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":3540246,"visible":true,"origin":"","legend":"","description":"","filename":"floatimage9.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/c300fc580173b59583c04b08.jpeg"},{"id":96710805,"identity":"72ab80dd-9b00-4eb0-a4a5-8aa01f6c6dbd","added_by":"auto","created_at":"2025-11-25 10:11:12","extension":"png","order_by":49,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":77333,"visible":true,"origin":"","legend":"","description":"","filename":"Onlinefloatimage1.png","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/e8eb31be2b084c1b43177754.png"},{"id":96710841,"identity":"dc699ac8-bdf9-490e-b53f-f7ea0b3f9b90","added_by":"auto","created_at":"2025-11-25 10:11:16","extension":"png","order_by":50,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":40394,"visible":true,"origin":"","legend":"","description":"","filename":"Onlinefloatimage2.png","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/c1981c7d6a7f319878826118.png"},{"id":96710991,"identity":"f3705608-f4ff-41c7-9b4f-d832908c914a","added_by":"auto","created_at":"2025-11-25 10:11:29","extension":"png","order_by":51,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":41589,"visible":true,"origin":"","legend":"","description":"","filename":"Onlinefloatimage3.png","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/76c791ebde007323ab3afc35.png"},{"id":96689348,"identity":"7a114449-2609-46f7-8c05-5e72d06e023e","added_by":"auto","created_at":"2025-11-25 06:09:06","extension":"png","order_by":52,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":161579,"visible":true,"origin":"","legend":"","description":"","filename":"Onlinefloatimage4.png","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/67ab53798c471de2321b109d.png"},{"id":96689317,"identity":"51f00a45-c520-49f0-b9b7-f5f357aa6000","added_by":"auto","created_at":"2025-11-25 06:09:05","extension":"png","order_by":53,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":28088,"visible":true,"origin":"","legend":"","description":"","filename":"Onlinefloatimage5.png","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/6eb2dd6725f0a11493d4caf8.png"},{"id":96689335,"identity":"3dfcf1a6-67ee-4aaa-aaea-4125596eded2","added_by":"auto","created_at":"2025-11-25 06:09:06","extension":"png","order_by":54,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":16329,"visible":true,"origin":"","legend":"","description":"","filename":"Onlinefloatimage6.png","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/9ae85512a7ef1d3dad0e6691.png"},{"id":96689344,"identity":"211c519b-deb2-4b0d-bd7f-9d6378fc04d0","added_by":"auto","created_at":"2025-11-25 06:09:06","extension":"png","order_by":55,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":53011,"visible":true,"origin":"","legend":"","description":"","filename":"Onlinefloatimage7.png","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/86c4e5b8977fd017b59454b7.png"},{"id":96710606,"identity":"ca15dfb2-a62d-49bb-9d55-1168c9365732","added_by":"auto","created_at":"2025-11-25 10:10:58","extension":"png","order_by":56,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":1060671,"visible":true,"origin":"","legend":"","description":"","filename":"Onlinefloatimage8.png","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/76c866d90f6947a1b2059065.png"},{"id":96689349,"identity":"941a5100-2e13-499e-ab7c-8462302b9ea2","added_by":"auto","created_at":"2025-11-25 06:09:06","extension":"png","order_by":57,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":26496,"visible":true,"origin":"","legend":"","description":"","filename":"Onlinefloatimage9.png","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/c4446f631b43655b603f8a51.png"},{"id":96711046,"identity":"fce4f585-0264-4918-9624-aaaa57f0edd9","added_by":"auto","created_at":"2025-11-25 10:11:33","extension":"xml","order_by":58,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":102401,"visible":true,"origin":"","legend":"","description":"","filename":"9e709631659943ec87068723e1aeaff51structuring.xml","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/59fc8e7d707ff756dda4ddee.xml"},{"id":96689351,"identity":"e6ac4f7e-9664-4491-bc59-73c4b021988a","added_by":"auto","created_at":"2025-11-25 06:09:06","extension":"html","order_by":59,"title":"","display":"","copyAsset":false,"role":"acdc-reference","size":116342,"visible":true,"origin":"","legend":"","description":"","filename":"earlyproof.html","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/cc3573006aef3f3cc85b3783.html"},{"id":96710688,"identity":"0455b827-5d5a-4fd7-a0d1-2f336c3d70fd","added_by":"auto","created_at":"2025-11-25 10:11:04","extension":"png","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":27363,"visible":true,"origin":"","legend":"\u003cp\u003eSample clustering for drone embryos of \u003cem\u003eA.cerana\u003c/em\u003e. (a) Heatmap and hierarchical clustering of the embryo samples using the whole transcriptome data. (b) Principal component analysis (PCA) of the embryo samples using the whole transcriptome data. (D:Drone)\u003c/p\u003e","description":"","filename":"Fig1.png","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/ae1ce3f680e0b931f98ee4b0.png"},{"id":96711087,"identity":"292a5204-2611-447a-9c79-598f757447f3","added_by":"auto","created_at":"2025-11-25 10:11:38","extension":"png","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":53561,"visible":true,"origin":"","legend":"\u003cp\u003eDifferential gene expression analysis. (a) Statistical data of differentially expressed genes. (b) Volcano plots of the five stages of drone embryonic development. (c) Heatmap of differentially expressed genes across the five stages of drone embryonic development (base_mean: absolute gene expression level; length: gene length).\u003c/p\u003e","description":"","filename":"Fig2.png","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/0bbd722bfc7b6cf0bb668d71.png"},{"id":96689273,"identity":"fd695e45-7a23-4eb9-aa9a-1763407f8fb0","added_by":"auto","created_at":"2025-11-25 06:09:03","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":143263,"visible":true,"origin":"","legend":"\u003cp\u003eillustrates the dynamic changes in gene expression during the embryonic development of drones. Using k-means clustering analysis, the gene expression patterns of drones were divided into 10 clusters. The heatmap on the left side of the Figure shows the expression levels of genes at each developmental stage, while the line plots on the right depict the trends in gene expression across five key time points during embryonic development (12h, 24h, 36h, 48h, and 60h). The clusters are labeled as C1 to C10, with X representing the drone embryo.\u003c/p\u003e","description":"","filename":"Fig3.png","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/05021ab683424aed8be4f748.png"},{"id":96710823,"identity":"2c328dd2-b847-4aa7-8f26-18f51de04996","added_by":"auto","created_at":"2025-11-25 10:11:13","extension":"png","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":34677,"visible":true,"origin":"","legend":"\u003cp\u003eConstruction of WGCNA modules. \u003cstrong\u003e(a)\u003c/strong\u003e Shows the distribution of gene numbers corresponding to the nine color modules.\u003cstrong\u003e(b, c)\u003c/strong\u003e Present the analysis of the scale-free fit index and mean connectivity for various soft-thresholding powers. The red line indicates where the correlation coefficient reaches 0.8, corresponding to a soft-thresholding power of 24. \u003cstrong\u003e(d)\u003c/strong\u003e Displays the heatmap analysis of the correlation between the five developmental stages and gene expression. The thickness of the lines reflects the strength of the correlation, with brown lines representing highly significant P-values and green lines representing significant P-values; the stars indicate the correlation between each color module, with larger stars corresponding to higher absolute correlation values.\u003c/p\u003e","description":"","filename":"Fig4.png","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/e96e7ab588411943307584ca.png"},{"id":96689277,"identity":"bd886584-5cd7-401e-b56e-32b2f6461054","added_by":"auto","created_at":"2025-11-25 06:09:04","extension":"png","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":213555,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cem\u003eGene expression trends of six clusters across five embryonic developmental stages in \u003c/em\u003e\u003cem\u003e\u003cstrong\u003eA. cerana\u003c/strong\u003e\u003c/em\u003e\u003cstrong\u003e.\u003c/strong\u003e(Each panel represents one gene cluster obtained from K-means clustering of WGCNA-derived modules, with the number of genes indicated in parentheses. The gray lines denote expression trajectories of individual genes, while the colored solid line represents the cluster centroid (average expression pattern). The red dashed line indicates the fitted linear trend of the cluster, with the slope shown in parentheses. Positive slopes reflect overall upregulation, negative slopes reflect downregulation, and slopes near zero indicate stable expression across developmental stages (12 h, 24 h, 36 h, 48 h, and 60 h)).\u003c/p\u003e","description":"","filename":"Fig5.png","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/261b55130fea97d3b67b3846.png"},{"id":96711604,"identity":"84416459-0030-4a75-8dde-e88ae7299dbb","added_by":"auto","created_at":"2025-11-25 10:12:18","extension":"png","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":143307,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eComparison of correctly predicted genes and overlaps among models.\u003cbr\u003e\n \u003c/strong\u003e(a) The number of correctly predicted genes by ResNet, ExtraTrees, and KNN. ResNet achieved the highest number of correctly predicted genes (2251), followed by ExtraTrees (2186) and KNN (2148).(b) Gene overlap matrix showing the number of shared correctly predicted genes between models. Diagonal values represent the total correctly predicted genes for each model, while off-diagonal values represent the overlaps between two models (e.g., ResNet–ExtraTrees: 2156; ResNet–KNN: 2118; ExtraTrees–KNN: 2098). The color scale indicates the number of shared genes, with darker colors representing larger overlaps.\u003c/p\u003e","description":"","filename":"Fig6.png","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/f788acea0802c0cee2429011.png"},{"id":96711036,"identity":"542d32ec-7097-49e2-9eee-0b684402ea4b","added_by":"auto","created_at":"2025-11-25 10:11:32","extension":"png","order_by":7,"title":"Figure 7","display":"","copyAsset":false,"role":"figure","size":673574,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eComprehensive performance comparison of 19 machine learning models for gene expression pattern classification.(a)\u003c/strong\u003e Accuracy comparison across 19 models, including 11 deep learning models (red) and 8 traditional machine learning models (blue). ResNet achieved the highest accuracy (0.9678), followed by ExtraTrees, KNN, SVM, and RandomForest, while WaveNet and InceptionTime showed the lowest performance.\u003cstrong\u003e(b)\u003c/strong\u003e F1 score comparison of the same models. ResNet achieved the best F1 score (0.958), confirming its superior classification ability, whereas WaveNet and InceptionTime again performed the worst.\u003cstrong\u003e(c)\u003c/strong\u003e Precision–recall trade-off plot for all models. The best-performing models (ResNet, ExtraTrees, KNN, SVM, RandomForest) cluster near the top-right corner, indicating a good balance between precision and recall.\u003cstrong\u003e(d)\u003c/strong\u003e Radar chart of the top five models (ResNet, ExtraTrees, KNN, SVM, RandomForest) showing performance across four evaluation metrics (accuracy, precision, recall, and F1 score). ResNet demonstrated consistently superior performance across all metrics.\u003c/p\u003e","description":"","filename":"Fig7.png","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/7c88b559e5f03e44b9d16383.png"},{"id":96711014,"identity":"acf93560-0b04-407d-a389-537fd6e57275","added_by":"auto","created_at":"2025-11-25 10:11:30","extension":"png","order_by":8,"title":"Figure 8","display":"","copyAsset":false,"role":"figure","size":235417,"visible":true,"origin":"","legend":"\u003cp\u003ePPI networks across five embryonic stages built from ResNet-prioritized genes\u003c/p\u003e","description":"","filename":"Fig8.png","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/a1071e1491ab3e707035dd70.png"},{"id":96689289,"identity":"5c72f5ee-49ab-4f35-9523-d6b74b64783c","added_by":"auto","created_at":"2025-11-25 06:09:04","extension":"png","order_by":9,"title":"Figure 9","display":"","copyAsset":false,"role":"figure","size":78810,"visible":true,"origin":"","legend":"\u003cp\u003eQuantitative polymerase chain reaction (qPCR) results demonstrate the expression changes of genes in drone embryos, as identified through WGCNA and Deep learning analysis.(t-test, p \u0026lt; 0.05). (* p \u0026lt; 0.05, ** p \u0026lt; 0.01, *** p \u0026lt; 0.001)\u003c/p\u003e","description":"","filename":"Fig9.png","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/2b3e0f94134c3b8b7049ab6a.png"},{"id":98622290,"identity":"b1e0fb54-994c-4cad-bce7-9ff4eb3e3cd5","added_by":"auto","created_at":"2025-12-19 16:51:00","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":2759223,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/d581096b-b775-44fb-ac82-47e92fff1b50.pdf"},{"id":96689282,"identity":"f6f7fff2-1501-4367-83d3-4bf10489d896","added_by":"auto","created_at":"2025-11-25 06:09:04","extension":"pdf","order_by":1,"title":"","display":"","copyAsset":false,"role":"supplement","size":87372,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementalFigure1.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/aeeb4801bdc13204e5660824.pdf"},{"id":96710731,"identity":"81e5d925-5037-480c-bfa6-81a37574e7b1","added_by":"auto","created_at":"2025-11-25 10:11:08","extension":"pdf","order_by":2,"title":"","display":"","copyAsset":false,"role":"supplement","size":10503377,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementalFigure2.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/7b8cca1bc1bf7706c52490b0.pdf"},{"id":96689292,"identity":"acca6e6f-f50f-4030-9d8f-3ff164f876c8","added_by":"auto","created_at":"2025-11-25 06:09:04","extension":"pdf","order_by":3,"title":"","display":"","copyAsset":false,"role":"supplement","size":461687,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementalFigure3.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/e10c297dda932e325c3fc9ee.pdf"},{"id":96689283,"identity":"88b605b6-24aa-4c76-adc9-a9ef4717ed0d","added_by":"auto","created_at":"2025-11-25 06:09:04","extension":"pdf","order_by":4,"title":"","display":"","copyAsset":false,"role":"supplement","size":1872614,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementalFigure4.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/570898b785be8192c3c4af6e.pdf"},{"id":96689295,"identity":"c3c38b3e-51b6-4133-9379-64b68b9631ad","added_by":"auto","created_at":"2025-11-25 06:09:04","extension":"pdf","order_by":5,"title":"","display":"","copyAsset":false,"role":"supplement","size":485338,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementalFigure5.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/4b96b42c7e633398e47e340c.pdf"},{"id":96710747,"identity":"9cc776c8-6a56-4678-ab6a-e0ce7143ec31","added_by":"auto","created_at":"2025-11-25 10:11:09","extension":"xlsx","order_by":6,"title":"","display":"","copyAsset":false,"role":"supplement","size":10738,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementalTable1.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/257ab7c574d3b8729f4c6a89.xlsx"},{"id":96689274,"identity":"fcc44a3a-a859-45cb-8583-c3359fa98b72","added_by":"auto","created_at":"2025-11-25 06:09:04","extension":"xlsx","order_by":7,"title":"","display":"","copyAsset":false,"role":"supplement","size":10309,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementalTable2.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/ced59b9aed85ad04ff0c2f4b.xlsx"},{"id":96689275,"identity":"9e35d636-49b9-4258-bc61-b552633b7b0b","added_by":"auto","created_at":"2025-11-25 06:09:04","extension":"xlsx","order_by":8,"title":"","display":"","copyAsset":false,"role":"supplement","size":15774,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementalTable3.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/bd1d67cfe372fad0610b631c.xlsx"},{"id":96689299,"identity":"71e08668-d9e6-4630-968d-2d5bca88fcc2","added_by":"auto","created_at":"2025-11-25 06:09:04","extension":"csv","order_by":9,"title":"","display":"","copyAsset":false,"role":"supplement","size":87362,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementalTable4a.csv","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/ca1a00d91047ebeaa5f43614.csv"},{"id":96711054,"identity":"a3f127cd-d101-4887-a858-a6a937d702ff","added_by":"auto","created_at":"2025-11-25 10:11:34","extension":"csv","order_by":10,"title":"","display":"","copyAsset":false,"role":"supplement","size":17173,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementalTable4b.csv","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/a893c0b24287c2b55d659b26.csv"},{"id":96689279,"identity":"96f8f597-916c-4298-954b-4384b3bd7105","added_by":"auto","created_at":"2025-11-25 06:09:04","extension":"csv","order_by":11,"title":"","display":"","copyAsset":false,"role":"supplement","size":12162,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementalTable4c.csv","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/6dbcd06605bd2b59e33c6679.csv"},{"id":96710980,"identity":"6a8667c5-e020-4b12-b350-1b60c458c145","added_by":"auto","created_at":"2025-11-25 10:11:28","extension":"csv","order_by":12,"title":"","display":"","copyAsset":false,"role":"supplement","size":26785,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementalTable4d.csv","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/c8b5d708127ec71fefdb1362.csv"},{"id":96689325,"identity":"62fabaa9-89d8-453f-a7e8-75fb41beb9d5","added_by":"auto","created_at":"2025-11-25 06:09:05","extension":"csv","order_by":13,"title":"","display":"","copyAsset":false,"role":"supplement","size":13281,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementalTable4e.csv","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/8d1227a07fb93d3648542453.csv"},{"id":96689291,"identity":"b12e5684-a80f-4727-bb64-1caffc95370c","added_by":"auto","created_at":"2025-11-25 06:09:04","extension":"csv","order_by":14,"title":"","display":"","copyAsset":false,"role":"supplement","size":26669,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementalTable5a.csv","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/4f7afaae9e8866b359f57248.csv"},{"id":96689304,"identity":"afcf05fe-3847-4ac7-b6b1-6fe358f5bf90","added_by":"auto","created_at":"2025-11-25 06:09:05","extension":"csv","order_by":15,"title":"","display":"","copyAsset":false,"role":"supplement","size":27157,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementalTable5b.csv","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/f07cb2450203fbb3288a9ae4.csv"},{"id":96711034,"identity":"0dc8425f-b246-419e-9697-16ad897216cf","added_by":"auto","created_at":"2025-11-25 10:11:32","extension":"csv","order_by":16,"title":"","display":"","copyAsset":false,"role":"supplement","size":9282,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementalTable5c.csv","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/efb0827c9cb5df4788b408d4.csv"},{"id":96689276,"identity":"ce5763bf-c7b7-4a21-85c3-8a5009a3d158","added_by":"auto","created_at":"2025-11-25 06:09:04","extension":"csv","order_by":17,"title":"","display":"","copyAsset":false,"role":"supplement","size":47103,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementalTable5e.csv","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/18c92212ea82f8770af350e8.csv"},{"id":96689321,"identity":"d583bf5c-0476-4365-8a87-b5cb025b9769","added_by":"auto","created_at":"2025-11-25 06:09:05","extension":"xlsx","order_by":18,"title":"","display":"","copyAsset":false,"role":"supplement","size":16493,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementalTable6.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/a7fb2e1e3c702b840a4c5755.xlsx"},{"id":96689341,"identity":"3b583082-c5f3-414a-a3cd-5dbab691736c","added_by":"auto","created_at":"2025-11-25 06:09:06","extension":"xlsx","order_by":19,"title":"","display":"","copyAsset":false,"role":"supplement","size":13072,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementalTable7.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-7939377/v1/2301d920dfd650eee20a7ddd.xlsx"}],"financialInterests":"No competing interests reported.","formattedTitle":"Integrative co-expression and deep learning reveal stage-specific regulators of Apis cerana drone embryogenesis","fulltext":[{"header":"1 Introduction","content":"\u003cp\u003eThe Asian honeybee (\u003cem\u003eApis cerana\u003c/em\u003e) is an important social insect that not only provides us with honey and beeswax, but also plays a crucial role as a pollinator in ecosystems \u003csup\u003e[\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e]\u003c/sup\u003e. Embryonic development is the starting phase of honeybee lifecycle, during which the main organs of adult bees progressively develop (REF). So studying its developmental mechanisms is essential for understanding honeybee growth and reproduction \u003csup\u003e[\u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e]\u003c/sup\u003e.Studies on the embryo of honeybee workers (\u003cem\u003eApis mellifera\u003c/em\u003e) have revealed that it undergoes 10 successive developmental phases, lasting 72 h \u003csup\u003e[\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e, \u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e]\u003c/sup\u003e. And significant morphological changes occur and tissues and organs progressively develop from stage 5 to stage 9 \u003csup\u003e[\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e]\u003c/sup\u003e.The strongly expressed proteins and pathways related to transcriptional\u0026thinsp;\u0026minus;\u0026thinsp;translational machinery and morphogenesis at 24 h drone embryo relative to the worker, illustrating the earlier occurrence of morphogenesis in the drone than worker \u003csup\u003e[\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e]\u003c/sup\u003e.Proteomic analysis showed that in the 24-hour stage, proteins and pathways related to transcription-translation mechanisms and morphogenesis show significant expression, indicating that molecular preparation for morphogenesis likely begins before 24 hours \u003csup\u003e[\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e]\u003c/sup\u003e. Moreover, phosphorylated proteins associated with morphogenesis represent the highest proportion \u003csup\u003e[\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e]\u003c/sup\u003e. The 48-hour stage is the most active period of protein expression during embryonic development, primarily involving proteins related to carbohydrate metabolism and energy production \u003csup\u003e[\u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e8\u003c/span\u003e]\u003c/sup\u003eThe haploid drones, which develop from unfertilized eggs (REF), are differ from the fertilized diploid the fertilized eggs giving rise to diploid workers are an ideal model for studying the morphological and biochemical characteristics of recessive mutations \u003csup\u003e[\u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e]\u003c/sup\u003e.In previous transcriptomic studies of worker bee embryos, 24, 48, and 72 hours were typically chosen as sampling points, roughly corresponding to the early, middle, and late stages of the embryo's transition into a larva \u003csup\u003e[\u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e10\u003c/span\u003e]\u003c/sup\u003e. Based on these studies, we added three more time points: 12h, 36h, and 60h, aiming for a more detailed transcriptomic analysis of drone embryos. The 12h and 60h time points were selected from the longest developmental periods in worker bee development. The 12h stage is crucial for studying early gene expression in drone embryos, during which significant changes in gene and microRNA expression occur \u003csup\u003e[\u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e11\u003c/span\u003e]\u003c/sup\u003e, while 36h represents a key middle point during major morphological changes in worker embryos. To date, there has been little transcriptomic research focusing on these three time points. In this study, we conducted an in-depth analysis of gene expression in the drone embryos of \u003cem\u003eApis cerana\u003c/em\u003e at five key developmental stages (12h, 24h, 36h, 48h, 60h). Through RNA-seq analysis, we observed significant differences in gene expression across these stages, with the most notable changes occurring between 12h and 24h. Additionally, we employed temporal transcriptome analysis and weighted gene co-expression network analysis (WGCNA) to identify stage-specific genes at these five time points, providing deeper insights into their molecular mechanisms. We also discovered potential key genes and identified several important but functionally unknown genes in \u003cem\u003eApis cerana\u003c/em\u003e, which may be closely related to drone embryo development.The development of high-throughput sequencing technologies has made it possible to study the embryonic development of non-model insects at the genomic level, providing opportunities to identify new candidate genes \u003csup\u003e[\u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e12\u003c/span\u003e]\u003c/sup\u003e. While substantial genomic and transcriptomic data are available for the embryonic development of various insect species, including \u003cem\u003eD. melanogaster\u003c/em\u003e \u003csup\u003e[\u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e13\u003c/span\u003e]\u003c/sup\u003e, \u003cem\u003eG. bimaculatus\u003c/em\u003e\u003csup\u003e[\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e]\u003c/sup\u003e, \u003cem\u003eAthetis lepigone\u003c/em\u003e \u003csup\u003e[\u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e]\u003c/sup\u003e, \u003cem\u003eT. castaneum\u003c/em\u003e \u003csup\u003e[\u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e]\u003c/sup\u003e, \u003cem\u003eOncopeltus fasciatus\u003c/em\u003e \u003csup\u003e[\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e]\u003c/sup\u003e, \u003cem\u003eNilaparvata lugens\u003c/em\u003e \u003csup\u003e[\u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e18\u003c/span\u003e]\u003c/sup\u003e, and the silkworm (Bombyx mori) \u003csup\u003e[\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e]\u003c/sup\u003e, transcriptomic studies on the embryonic development of male honeybees remain relatively sparse. It was not until 2020 that the first high-quality genome assembly of \u003cem\u003eApis cerana\u003c/em\u003e was completed\u003csup\u003e[\u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e]\u003c/sup\u003e, which laid the foundation for precise gene annotation.\u003c/p\u003e\u003cp\u003eExisting research has shown that, within 24 hours, the expression of proteins and pathways related to transcription-translation mechanisms and morphogenesis significantly increases, suggesting that the molecular preparations for morphogenesis begin before the 24-hour mark. The highest proportion of phosphoproteins associated with morphogenesis was observed in the 24-hour phosphoproteome \u003csup\u003e[\u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e]\u003c/sup\u003e. At 48 hours, embryonic protein expression peaks, primarily involving proteins related to carbohydrate metabolism and energy production. At 72 hours, proteins related to antioxidation and cytoskeletal functions are significantly upregulated \u003csup\u003e[\u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e22\u003c/span\u003e]\u003c/sup\u003e. Although these studies reveal certain characteristics of protein expression during male honeybee embryogenesis, systematic studies on the dynamic changes in gene expression remain insufficient. To address this gap, we performed a time-series transcriptomic analysis of \u003cem\u003eA. cerana\u003c/em\u003e drone embryos at five key stages (12, 24, 36, 48, and 60 h) to delineate stage-specific regulatory frameworks at the transcriptional level. Methodologically, we first constructed co-expression modules using weighted gene co-expression network analysis (WGCNA), then refined within-module expression patterns with k-nearest neighbors (KNN) based on the five-timepoint trajectories. We next combined \u0026ldquo;module features\u0026thinsp;+\u0026thinsp;KNN subcluster labels/consistency\u0026rdquo; as inputs to benchmark 19 machine-/deep-learning models. This integrated strategy follows the recently advocated paradigm of \u0026ldquo;network-first \u0026rarr; model benchmarking \u0026rarr; interpretable validation\u0026rdquo;: WGCNA compresses high-dimensional expression into modules to denoise, after which models such as LASSO enable more precise screening of phenotype- or disease-relevant target genes, substantially reducing the search space\u003csup\u003e[\u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e23\u003c/span\u003e, \u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e24\u003c/span\u003e]\u003c/sup\u003e.In head-to-head comparisons, the residual network (ResNet) achieved the highest accuracy for gene-pattern classification (96.78%), with concordant training/validation curves and no signs of overfitting. To further enhance confidence, we applied a consensus-intersection of the top three performers (e.g., ResNet, ExtraTrees, KNN classifier) and integrated this with KNN subcluster consistency, thereby systematically narrowing the initial large candidate pool to \u0026gt;\u0026thinsp;2,000 high-confidence, stage-specific key genes. We then constructed protein\u0026ndash;protein interaction (PPI) networks in Cytoscape, highlighting subnetworks with \u0026ge;\u0026thinsp;20 nodes, and performed GO annotation, which resolved functional groups related to visual organs, wing discs, gonads, and metabolism. To experimentally support the computational predictions, we randomly selected 1\u0026ndash;2 highly expressed representative genes at each stage for qPCR validation; results were consistent with both model predictions and transcriptomic trends.\u003c/p\u003e\u003cp\u003eIn sum, we establish and implement a closed-loop time-series framework\u0026mdash;WGCNA \u0026rarr; KNN(unsupervised cluster analysis) \u0026rarr; multi-model ML/DL \u0026rarr; consensus intersection \u0026rarr; PPI/GO \u0026rarr; qPCR\u0026mdash;that maps stage-resolved gene networks underlying \u003cem\u003eA.cerana\u003c/em\u003e drone embryogenesis, from early RNA/epigenetic regulation through mid-stage organ morphogenesis to late-stage tissue specialization and metamorphic preparation. This framework couples accuracy with interpretability and provides an actionable shortlist of testable regulators and pathways, offering a practical template for subsequent functional genetics and breeding applications.\u003c/p\u003e"},{"header":"2 Materials and Methods","content":"\u003cdiv id=\"Sec3\" class=\"Section2\"\u003e\u003ch2\u003e2.2 Sample collection\u003c/h2\u003e\u003cp\u003eEmbryo samples of \u003cem\u003eApis cerana\u003c/em\u003e were obtained from three queenright colonies in Cangyuan, Yunan Province, China. To collect right-aged embryonic samples, mated egg-laying queens were limited in two frames with empty drone comb, newly laid eggs were monitored at six-hour intervals and were marked with different colored paint pens, the time of the first observed egg deposition was recorded as 0 h. After marking, the comb frame was returned to the hive until eggs from the 12 h, 24 h, 36 h, 48 h, and 60 h stages were obtained. At each stage, 20\u0026ndash;30 eggs were collected, lysed in lysis buffer (PRL), and immediately stored at -80\u0026deg;C. In total, samples related to five developmental stages with three biological replicates each were collected for RNA-seq.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec4\" class=\"Section2\"\u003e\u003ch2\u003e2.3 \u003cem\u003eRNA\u003c/em\u003e extraction and cDNA library construction\u003c/h2\u003e\u003cp\u003eTotal RNA was extracted from the embryo samples according to the protocol of GeneBetter\u0026rsquo;s R516 FlashPure Total RNA micro Kit, and then were preserved in lysis buffer and stored at -80\u0026deg;C. The purity and concentration of RNA were measured using NanoDrop 2000, and quantity of the total RNA samples were determined by using Agilent 2100/4200 system. cDNA libraries were prepared according to standard protocols and sequenced as 150-bp paired-end reads on the Illumina HiSeq platform by Anhui Gaohe Biotechnology Co., Ltd.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec5\" class=\"Section2\"\u003e\u003ch2\u003e2.4 Mapping and quality control of RNA-seq data\u003c/h2\u003e\u003cp\u003eBefore mapping, raw data were filtered by the following steps: (1) reads with adapters were removed; (2) reads in which unknown bases represented more than 10% of the total sequence were removed; (3) low quality reads, in which the percentage of low-quality bases (quality no more than 10) exceeded 50%, were removed. The remaining reads were defined as \u0026ldquo;clean reads\u0026rdquo; and used for subsequent bioinformatics analyses \u003csup\u003e[\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e]\u003c/sup\u003e. Quality control on the clean data was performed by the production of a base composition chart and quality distribution chart using FastQC. We used Hisat2(v2.2.1) \u003csup\u003e[\u003cspan citationid=\"CR26\" class=\"CitationRef\"\u003e26\u003c/span\u003e]\u003c/sup\u003e to map clean reads to the reference genome (GCA_029169275.1). The SAM files were converted to BAM files and sorted using samtools (v1.6) \u003csup\u003e[\u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e27\u003c/span\u003e]\u003c/sup\u003e. These sorted BAM files were then quantified by passing them to Feature Counts \u003csup\u003e[\u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e28\u003c/span\u003e]\u003c/sup\u003e to count the number of reads mapped to each gene. Finally, differential expression analysis was conducted using the DESeq2 algorithm implemented in Trinity (v2.15.1\u003cb\u003e)\u003c/b\u003e software \u003csup\u003e[\u003cspan citationid=\"CR29\" class=\"CitationRef\"\u003e29\u003c/span\u003e]\u003c/sup\u003e.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec6\" class=\"Section2\"\u003e\u003ch2\u003e2.5 PCA analysis\u003c/h2\u003e\u003cp\u003eTo evaluate the consistency of sample collection and investigate the transcriptomic relationships between 12h, 24h, 36h, 48h, and 60h embryos, we performed hierarchical clustering analysis and principal component analysis (PCA) on the full gene expression data of these samples. Both methods effectively captured the similarities and differences between samples. Before conducting these analyses, gene expression levels were standardized using the regularized log transformation method in the limma package.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec7\" class=\"Section2\"\u003e\u003ch2\u003e2.6 GO and KEGG analysis\u003c/h2\u003e\u003cp\u003eThe longest cds protein sequence from the reference \u003cem\u003eApis cerana\u003c/em\u003e genome was annotated for GO/KEGG using emapper.py (2.1.10) with the eggnog database. The annotated file was converted into org.db format, and GO and KEGG enrichment analysis was performed using the clusterProfiler package (v 4.4) in R.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec8\" class=\"Section2\"\u003e\u003ch2\u003e2.7 Temporal clustering analysis of genes\u003c/h2\u003e\u003cp\u003eUsing the ClusterGVis package in R, time-series transcriptome analysis was performed on significantly expressed genes ( |log₂FoldChange| \u0026gt;1 and padj\u0026thinsp;\u0026lt;\u0026thinsp;0.05) from the expression matrix. Through clustering genes with similar expression patterns, the gene modules of different clusters were analyzed to determine their expression profiles across the various embryonic stages.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec9\" class=\"Section2\"\u003e\u003ch2\u003e2.8 WGCNA, KNN refinement, and downstream deep-learning\u0026ndash;based prioritization\u003c/h2\u003e\u003cp\u003eRaw counts were normalized, lowly expressed genes were filtered, and outliers were removed based on inter-sample correlations. A weighted co-expression network was constructed using WGCNA (R package WGCNA v1.xx): the soft-thresholding power (β) was chosen with pickSoftThreshold to approximate scale-free topology (target R\u0026sup2; \u0026ge; 0.8), a signed adjacency matrix was computed and transformed to the topological overlap matrix (TOM), modules were detected by dynamic tree cutting, and module eigengenes were correlated with the five developmental stages (12、24、36、48and60 h).\u003c/p\u003e\u003cp\u003eAfter module detection, we first refined gene-level temporal patterns with KNN clustering within each module, using five-stage trajectories as features (K selected by silhouette/elbow criteria). Each gene received a KNN subcluster label and a neighborhood-consistency score. We then formed a composite feature set (module-level normalized profiles\u0026thinsp;+\u0026thinsp;KNN labels/scores) and fed it into a multi-model machine-/deep-learning pipeline for stage classification and gene prioritization. In total, 19 models were benchmarked (11 deep-learning models: RNN, LSTM, BiLSTM, GRU, CNN, ResNet, Transformer, AttentionRNN, TCN, WaveNet, InceptionTime; and 8 traditional models: e.g., Random Forest, Gradient Boosting, AdaBoost). Data were split 60/20/20 into training/validation/test sets and evaluated with 5-fold cross-validation. The best-performing ResNet used a two-block residual architecture (fully connected 128\u0026ndash;64), batch normalization, dropout\u0026thinsp;=\u0026thinsp;0.3, and was trained with Adam (lr\u0026thinsp;=\u0026thinsp;0.001) and cross-entropy loss for 100 epochs.\u003c/p\u003e\u003cp\u003eFor high-confidence calls, we applied a consensus intersection of the top three models (e.g., ResNet, ExtraTrees, KNN classifier) and integrated this with KNN-subcluster consistency, yielding stage-specific key-gene sets across all five time points. Subsequent PPI/GO/qPCR analyses and visualizations are described in dedicated sections.\u003c/p\u003e\u003c/div\u003e"},{"header":"3 Result","content":"\u003cdiv id=\"Sec11\" class=\"Section2\"\u003e\u003ch2\u003e3.1 Overall RNA-Seq results and quality control\u003c/h2\u003e\u003cp\u003eIllumina RNA sequencing generated an average of 66.52\u0026nbsp;million (M) raw reads per sample. After filtering, an average of 66.20\u0026nbsp;million clean reads per sample were retained for subsequent transcriptome analysis. The clean reads had high quality, with average Q20 and Q30 scores of 98.84% and 96.31%, respectively. Additionally, the average GC content was 39.49% (Supplemental Table \u003cspan refid=\"MOESM1\" class=\"InternalRef\"\u003eS1\u003c/span\u003e). HISAT2 software was used for alignment against the reference genome. Sliding window density analysis showed that more than 80% of transcripts were uniquely aligned to the reference genome. Across 15 samples, a total of 61.82\u0026nbsp;million clean reads were uniquely mapped to the reference genome, with an average mapping rate of 89.86%. Only 10.14% of reads (7.97\u0026nbsp;million) failed to align with the genome (Supplemental Table \u003cspan refid=\"MOESM2\" class=\"InternalRef\"\u003eS2\u003c/span\u003e).\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec12\" class=\"Section2\"\u003e\u003ch2\u003e3.2 Hierarchical clustering and principal component analysis of samples\u003c/h2\u003e\u003cp\u003eThe results of hierarchical clustering and PCA (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003ea, b) demonstrated that the samples collected at the same time point had good homogeneity. The genes related to 12h were significantly different from the other four stages (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eb), which exhibited the highest expression levels (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003ea). It indicated that gene expression profiles at 12h were markedly distinct from the other four stages.\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec13\" class=\"Section2\"\u003e\u003ch2\u003e3.3 Analysis of differentially expressed genes (DEGs)\u003c/h2\u003e\u003cp\u003eTo study the differential expression of genes during the development of drone embryos, we conducted a statistical analysis of differentially expressed genes (DEGs) across five developmental stages of drone embryos, comparing gene expression differences between 12h and 24h, 24h and 36h, 36h and 48h, and 48h and 60h. The results (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ea) showed the following: between 12h and 24h, 2538 genes were significantly upregulated, and 3282 genes were significantly downregulated; between 24h and 36h, 1004 genes were significantly upregulated, and 1723 genes were significantly downregulated; between 36h and 48h, 888 genes were significantly upregulated, and 574 genes were significantly downregulated; and between 48h and 60h, 328 genes were significantly upregulated, and 1268 genes were significantly downregulated. The most differentially expressed genes (both upregulated and downregulated) were observed between the 12h and 24h developmental stages (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eb), while the degree of gene expression differences between 36h, 48h, and 60h gradually decreased. The heatmap of differentially expressed genes (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ec) indicates clear stage-specific gene expression across the five developmental stages. These stage-specific genes are likely associated with physiological, biochemical processes, morphological changes, and behavioral activities, potentially playing an important role in the development of drone embryos.\u003c/p\u003e\u003cp\u003eWe selected the top 20 genes with the most significant expression differences at each stage, filtering out 21 non-coding genes \u003cb\u003e(\u003c/b\u003eFig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eb\u003cb\u003e)\u003c/b\u003e. Through EggNOG database annotation and protein function prediction, we found that these genes are primarily involved in functions such as DNA-binding transcription factors, regulation of gene expression, G-protein-coupled receptors, cytoskeleton binding, methylation, negative regulation of programmed cell death, protein transport, phosphorylation, and molecular interactions during embryonic development. Additionally, we identified five genes (\u003cem\u003eLOC100578606, LOC100578389, LOC100577882, LOC551929, LOC725838\u003c/em\u003e) whose functions in \u003cem\u003eApis cerana\u003c/em\u003e are still unknown. These genes may play a critical role in drone embryo development, but further experiments are needed to verify their specific functions.\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003e\u003cb\u003e3.4 Gene expression in drone embryos is co-regulated across 10 clusters, which correspond to different developmental stages.\u003c/b\u003e\u003c/p\u003e\u003cp\u003eTo further explore the dynamic changes in gene expression during drone embryo development, we applied the k-means clustering algorithm, dividing 9492 genes from the reference genome into 10 clusters based on the similarity of their expression patterns. Detailed analysis of these 10 clusters revealed significant gene enrichment at specific developmental stages (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003e). The C1 and C2 modules were highly expressed significantly enriched at the 12-hour time point of drone embryo development, with related GO terms mainly involving sperm cell differentiation, maintenance of germline stem cell populations, male gametogenesis, and functions associated with protein modifications such as phosphorylation, methylation, epigenetic regulation of gene expression, and mRNA regulation or degradation. On the other hand, genes in the C5 and C8 modules showed higher expression after 24 hours of embryo development, primarily related to processes such as cell growth and differentiation, transcription regulation, and rRNA and ncRNA metabolism.\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec14\" class=\"Section2\"\u003e\u003ch2\u003e3.5 Co-Expression Network Construction and Hub Module Identification\u003c/h2\u003e\u003cp\u003eWeighted gene co-expression network analysis (WGCNA) were conducted to explore the association of genes at five specific developmental stages using the all aligned genes (9492). The samples were clustered (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003ea), missing values were processed, and outliers were excluded. To ensure that the network exhibited scale-free characteristics, a soft-threshold power of 24 was selected (R\u0026sup2; for the scale-free network\u0026thinsp;=\u0026thinsp;0.80) (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eb,c). A co-expression matrix was constructed using the one-step method, and dynamic tree cutting identified nine gene modules (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003ea,d), which demonstrated their correlation with drone embryo development.\u003c/p\u003e\u003cp\u003eThe results showed that the turquoise module (containing 2643 genes) had the highest correlation with the 12-hour embryo development stage (cor\u0026thinsp;=\u0026thinsp;0.69; P\u0026thinsp;=\u0026thinsp;0.004) (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003ed); the red and blue modules had the highest correlation with the 24-hour stage; the black module was most correlated with the 36-hour stage (cor\u0026thinsp;=\u0026thinsp;0.69; P\u0026thinsp;=\u0026thinsp;0.004); the brown and grey modules were associated with the 48-hour stage (cor\u0026thinsp;=\u0026thinsp;0.69; P\u0026thinsp;=\u0026thinsp;0.004); and the yellow, pink, and green modules were highly correlated with the 60-hour stage (cor\u0026thinsp;=\u0026thinsp;0.69; P\u0026thinsp;=\u0026thinsp;0.004) (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003ed). The differences in the correlation of these nine modules across the developmental stages of drone embryos indicate that gene expression at each stage exhibits clear specificity.\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec15\" class=\"Section2\"\u003e\u003ch2\u003e3.6 Integrative WGCNA and Deep Learning Identify Stage-Specific Genes in \u003cem\u003eApis cerana\u003c/em\u003e Embryogenesis\u003c/h2\u003e\u003cp\u003eIn this study, we first conducted a systematic analysis of gene co-expression modules underlying \u003cem\u003eApis cerana\u003c/em\u003e embryonic development using WGCNA. Although only five developmental time points were sampled, heterogeneity in temporal expression dynamics enabled further refinement: by combining K-means clustering with module\u0026ndash;trait association analysis, we partitioned the data into six representative expression modules (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003e). Building on these results, we applied machine-learning models to prioritize key genes, constructing and evaluating 19 models in total\u0026mdash;11 deep learning models (RNN, LSTM, BiLSTM, GRU, CNN, ResNet, Transformer, AttentionRNN, TCN, WaveNet, InceptionTime) and 8 traditional models (e.g., Random Forest, Gradient Boosting, AdaBoost). The ResNet model achieved the best performance for gene-expression pattern classification, reaching an accuracy of 96.78% (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003e\u0026ndash;8,Supplemental Fig.\u0026nbsp;1). Notably, training and validation learning curves were highly concordant, indicating no overfitting and demonstrating strong generalization and robustness. Among all candidates, we further selected the top three performers and intersected their predictions, yielding a robust set of \u0026gt;\u0026thinsp;2,000 genes with putative stage-specific roles, which were subsequently categorized for downstream functional analyses (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003e\u0026ndash;8,Supplemental Fig.\u0026nbsp;1).\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003e(a) The number of correctly predicted genes by ResNet, ExtraTrees, and KNN. ResNet achieved the highest number of correctly predicted genes (2251), followed by ExtraTrees (2186) and KNN (2148).(b) Gene overlap matrix showing the number of shared correctly predicted genes between models. Diagonal values represent the total correctly predicted genes for each model, while off-diagonal values represent the overlaps between two models (e.g., ResNet\u0026ndash;ExtraTrees: 2156; ResNet\u0026ndash;KNN: 2118; ExtraTrees\u0026ndash;KNN: 2098). The color scale indicates the number of shared genes, with darker colors representing larger overlaps.\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec16\" class=\"Section2\"\u003e\u003ch2\u003e3.7 PPI network construction and hub genes identification\u003c/h2\u003e\u003cp\u003eBuilding on the WGCNA-derived co-expression modules, we applied KNN clustering to further refine temporal expression subgroups within modules and, based on these refined clusters, selected candidate genes for downstream analyses. We then constructed protein\u0026ndash;protein interaction (PPI) networks in Cytoscape and visualized subnetworks with \u0026ge;\u0026thinsp;20 nodes (Fig.\u0026nbsp;8). Using network topology metrics, stage-specific hub genes were prioritized; these hubs aggregated into functional modules that collectively drive key physiological and biochemical processes during honeybee embryogenesis (Fig.\u0026nbsp;8). In addition, 1\u0026ndash;2 highly expressed representative genes per stage were randomly selected for qPCR validation, and the results were concordant with both model predictions and transcriptomic trends, supporting the reliability and reproducibility of the identified hubs and modules. Finally, Gene Ontology (GO) annotation of these genes delineated their stage-specific functional roles (Fig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e9\u003c/span\u003e; Supplemental Tables\u0026nbsp;3 and 5).\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003e\u003cb\u003eFugure 8\u003c/b\u003e. PPI networks across five embryonic stages built from ResNet-prioritized genes\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003c/div\u003e"},{"header":"4 Discussion","content":"\u003cp\u003eAt the 12h early developmental stage, temporal transcriptomic dynamics were closely associated with cytoskeletal remodeling, precise chromosomal positioning, and the formation of male germ cells (Fig.\u0026nbsp;8a). In parallel, epigenetic and phosphorylation-associated layers likely contribute to stage-specific control, consistent with systems-level evidence implicating TOR/epigenetic axes and phospho-regulatory drivers in developmental regulation \u003csup\u003e[\u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e30\u003c/span\u003e, \u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e31\u003c/span\u003e]\u003c/sup\u003e. Concordant results from WGCNA, deep-learning\u0026ndash;guided analyses, and temporal clustering (Supplemental Tables\u0026nbsp;3a, 3b; Table\u0026nbsp;4a) further support the importance of these early programs. Collectively, these observations indicate that multiple regulatory layers act in concert to stabilize and specify developmental trajectories at 12 h.\u003c/p\u003e\u003cp\u003eIn 24h embryos (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003e,Supplemental Table\u0026nbsp;4b), GO enrichment of Cluster 5 and Cluster 8 highlighted potentially critical roles for ncRNAs and rRNAs. Foundational studies have established that ncRNAs\u0026mdash;including microRNAs and lncRNAs\u0026mdash;govern chromatin states, transcriptional outputs, RNA processing, and translation, thereby shaping cell division, differentiation, and developmental timing \u003csup\u003e[\u003cspan additionalcitationids=\"CR33 CR34\" citationid=\"CR32\" class=\"CitationRef\"\u003e32\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR35\" class=\"CitationRef\"\u003e35\u003c/span\u003e]\u003c/sup\u003e. Moreover, gene-set enrichment comparing 12 h with 24 h module genes showed higher 24 h expression in sets related to compound-eye photoreceptor differentiation, ncRNA processing, and rRNA processing (Fig.\u0026nbsp;8a,b;Supplemental Table\u0026nbsp;5a,b), reinforcing stage-specific deployment of RNA-centric regulation.\u003c/p\u003e\u003cp\u003eAt 36h, hub-gene GO enrichment indicated DNA binding, neuronal projection development, cell morphogenesis, and differentiation (Fig.\u0026nbsp;8c; Supplemental Tables\u0026nbsp;4c and 5c), consistent with accelerated organogenesis and functional specialization. Nutrient- and growth-sensing pathways (e.g., TOR) intersect with transcriptional/epigenetic circuits, providing a plausible framework for coupling metabolic state to morphogenesis at this stage \u003csup\u003e[\u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e30\u003c/span\u003e]\u003c/sup\u003e. GSEA comparing 24 hand 36 h further showed that genes related to male gamete generation and RNA splicing were significantly upregulated at 24 h (Supplemental Fig.\u0026nbsp;2). Given that alternative splicing (AS) broadly expands transcriptome/proteome diversity and profoundly impacts developmental phenotypes\u0026mdash;with recent studies underscoring rapid, organism-wide AS remodeling \u003csup\u003e[\u003cspan citationid=\"CR34\" class=\"CitationRef\"\u003e34\u003c/span\u003e, \u003cspan additionalcitationids=\"CR37\" citationid=\"CR36\" class=\"CitationRef\"\u003e36\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR38\" class=\"CitationRef\"\u003e38\u003c/span\u003e]\u003c/sup\u003e 。We infer that heightened AS activity at 24 h promotes expression control of sex-differentiation and developmental genes, with important consequences for drone physiology and biochemistry.\u003c/p\u003e\u003cp\u003eAt 48h, hub-gene GO terms centered on anterior\u0026ndash;posterior patterning in wing discs, adult walking behavior, gonad development, eye morphogenesis, sensory organ morphogenesis, and positive regulation of the \u003cem\u003eWnt\u003c/em\u003e pathway (Fig.\u0026nbsp;8d; Supplemental Tables\u0026nbsp;4d and 5d). This constellation is consistent with coordinated deployment of core developmental signaling during organogenesis in holometabolous insects, supporting polarity establishment, proliferation, and differentiation\u003csup\u003e[\u003cspan citationid=\"CR39\" class=\"CitationRef\"\u003e39\u003c/span\u003e, \u003cspan citationid=\"CR40\" class=\"CitationRef\"\u003e40\u003c/span\u003e]\u003c/sup\u003e.\u003c/p\u003e\u003cp\u003eAt 60h, temporal clustering (Cluster 9) and hub-gene enrichment jointly indicated robust organ morphogenesis and skeletal muscle development (Supplemental Table\u0026nbsp;4e; Table\u0026nbsp;5e). Homology analysis mapped \u003cem\u003eLOC410326\u003c/em\u003e to Drosophila \u003cem\u003ecbt\u003c/em\u003e; given the known roles of \u003cem\u003ecbt\u003c/em\u003e in cell-cycle control and tissue growth, the bee homolog likely contributes to stem-cell maintenance, rapid tissue proliferation, and tissue specializationorgan development, thereby laying the groundwork for subsequent metamorphosis. Thus, late embryogenesis is characterized by intensified, directionally patterned growth consistent with imminent metamorphic transitions\u003csup\u003e[\u003cspan citationid=\"CR41\" class=\"CitationRef\"\u003e41\u003c/span\u003e, \u003cspan citationid=\"CR42\" class=\"CitationRef\"\u003e42\u003c/span\u003e]\u003c/sup\u003e.\u003c/p\u003e\u003cp\u003eWe further dissected the deep-learning outputs with a composite scoring system to identify genes with pronounced temporal dynamics. For dynamic genes, Dynamic Score\u0026thinsp;=\u0026thinsp;R\u0026sup2; \u0026times; 20 + |slope| \u0026times; 15\u0026thinsp;+\u0026thinsp;log₂(fold_change\u0026thinsp;+\u0026thinsp;1) \u0026times; 5; for stable genes, Stability Score = |slope| \u0026times; 10\u0026thinsp;+\u0026thinsp;CV \u0026times; 10\u0026thinsp;+\u0026thinsp;Δmax \u0026times; 2 (lower scores indicate greater stability). Genes were stratified into four tiers: Tier 1 (strongly dynamic)\u0026mdash;Dynamic Score\u0026thinsp;\u0026gt;\u0026thinsp;30 and R\u0026sup2; \u0026gt;0.8 with monotonic trends (|slope| \u0026gt;1.0) and substantial changes (fold_change\u0026thinsp;\u0026gt;\u0026thinsp;4); Tier 2 (transient response)\u0026mdash;Peak/Valley patterns with fold_change\u0026thinsp;\u0026gt;\u0026thinsp;2; Tier 3 (complex regulation)\u0026mdash;Complex patterns or R\u0026sup2; \u0026lt; 0.5; Tier 4 (homeostatic)\u0026mdash;Stability Score\u0026thinsp;\u0026lt;\u0026thinsp;10 with CV\u0026thinsp;\u0026lt;\u0026thinsp;0.2 and |slope| \u0026lt; 0.1 (Supplemental Figure. 3). To validate stratification, we randomly selected the most pronounced Tier-1 genes and the lowest-scoring Tier-4 genes for box-plot visualization and qPCR (Supplemental Fig.\u0026nbsp;4a,b, 5a,b; Supplemental Tables\u0026nbsp;6, 7). Tier-1 genes exhibited consistent linear trends across time (mean R\u0026sup2; = 0.85\u0026thinsp;\u0026plusmn;\u0026thinsp;0.08), whereas Tier-4 genes remained stably expressed (mean CV\u0026thinsp;=\u0026thinsp;0.15\u0026thinsp;\u0026plusmn;\u0026thinsp;0.05). This framework strengthens reliability assessment via R\u0026sup2;, enhances biological interpretability via pattern label (Increasing/Decreasing/Peak/Valley/Stable/Complex), and optimizes identification of both dynamic and stable genes via dual scoring, thereby improving functional categorization and providing robust annotation for downstream model validation.\u003c/p\u003e"},{"header":"5 Conclusion","content":"\u003cp\u003eAltogether, our results provide new insights into the molecular underpinnings of \u003cem\u003eA. cerana\u003c/em\u003e drone embryo genesis. The discovery of previously uncharacterized but network-central genes points to additional developmental regulators. More broadly, we delineate gene networks mediating the transition from early germ-line/organ primordia to mature organ systems. By integrating WGCNA, deep learning (ResNet), and temporal transcriptomics, we defined stage-specific key genes, modules, and hubs, highlighting histone-linked/epigenetic regulation, ncRNAs, and sex-differentiation programs in drone embryogenesis and offering a theoretical foundation for bee developmental biology, breeding, and conservation.\u003c/p\u003e"},{"header":"Declarations","content":"\u003cp\u003e\u003cstrong\u003e6 Funding\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThis research is supported by the Yiwu Industrial and Commercial College Program under Grant Nos. XJKJ2502YB. This work was supported by Jinhua Key Laboratory of Robot Intelligent Welding Technology and by Jinhua Public Welfare Technology Application Research Project in 2024.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e7 Availability of data and materials\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe short-read RNA-Seq data involved in this study were submitted to the NCBI BioProject database (https://www.ncbi.nlm.nih.gov/bioproject/\u0026nbsp;accessed on 22 January 2025), with accession numbers PRJNA1186594.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e8 Authors\u003c/strong\u003e\u003cstrong\u003e\u0026rsquo;\u003c/strong\u003e\u003cstrong\u003econtributions\u0026nbsp;\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe genome assembly and analyses were conducted by Xiang Ding and Yehui Tan. The manuscript was written by Runlang Su and Dan Yue.\u0026nbsp;\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e9 Ethics approval and consent to participate\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eNot Applicable.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e10 Consent for publication\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eNot Applicable.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e11 Competing interests\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe authors declare that they have no competing interests.\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\u003cli\u003e\u003cspan\u003ePan, L-X., Hu, W-W., Cheng, F-P., Hu, X-F. \u0026amp; Wang, Z-L. Transcriptome analysis reveals differentially expressed genes between the ovary and testis of the honey bee Apis mellifera. \u003cem\u003eApidologie\u003c/em\u003e \u003cb\u003e53\u003c/b\u003e, 15 (2022).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eFang, Y. et al. In-depth Proteomics Characterization of Embryogenesis of the Honey Bee Worker. \u003cem\u003eMol. Cell. Proteom.\u003c/em\u003e \u003cb\u003e13\u003c/b\u003e, 2306\u0026ndash;2320 (2014).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eDuPraw, E. J. The honeybee embryo, (1967).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eHu, X., Ke, L. \u0026amp; Wang, Z. Zeng Z Dynamic transcriptome landscape of Asian domestic honeybee (\u003cem\u003eApis cerana\u003c/em\u003e) embryonic development revealed by high-quality RNA sequencing. \u003cem\u003eBMC Dev. Biol.\u003c/em\u003e \u003cb\u003e18\u003c/b\u003e, 11 (2018).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eCridge, A. G. et al. Dearden PK The honeybee as a model insect for developmental genetics. \u003cem\u003egenesis\u003c/em\u003e \u003cb\u003e55\u003c/b\u003e, e23019 (2017).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eFang, Y. et al. Analysis Unravels Mechanism Underling the Embryogenesis of the Honeybee Drone and Its Divergence with the Worker (\u003cem\u003eApis mellifera lingustica\u003c/em\u003e). \u003cem\u003eJ. Proteome Res.\u003c/em\u003e \u003cb\u003e14\u003c/b\u003e, 4059\u0026ndash;4071 (2015).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eMa, B., Ma, C., Li, J. \u0026amp; Fang, Y. Revealing phosphorylation regulatory networks during embryogenesis of honey bee worker and drone (\u003cem\u003eApis mellifera\u003c/em\u003e). \u003cem\u003eFront. Cell. Dev. Biol.\u003c/em\u003e \u003cb\u003e10\u003c/b\u003e, 1006964 (2022).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eTan, K. et al. \u0026amp; Gudrun Geographic distribution of the eastern honeybee, \u003cem\u003eApis cerana\u003c/em\u003e (Hymenoptera: Apidae), across ecological zones in China: morphological and molecular analyses. \u003cem\u003eSyst. Biodivers.\u003c/em\u003e : 473\u0026ndash;482, (2006).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eBeye, M. The dice of fate: the csd gene and how its allelic composition regulates sexual development in the honey bee. \u003cem\u003eApis mellifera BioEssays\u003c/em\u003e. \u003cb\u003e26\u003c/b\u003e, 1131\u0026ndash;1139 (2004).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eHu, X. F. \u0026amp; Ke, L. Zeng ZJ Histomorphological study on embryogenesis of the honeybee \u003cem\u003eApis cerana\u003c/em\u003e. \u003cem\u003eJ. Asia. Pac. Entomol.\u003c/em\u003e \u003cb\u003e22\u003c/b\u003e, 860\u0026ndash;867 (2019).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eC, M. W. D., Z, H. \u0026amp; Z, L. \u0026amp; L W Sex-Specific Development in Haplodiploid Honeybee Is Controlled by the Female-Embryo-Specific Activation of Thousands of Intronic LncRNAs. Front. (2021). Cell Dev. Biol,(.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eLi, J., Fang, Y. \u0026amp; Zhang, L. Begna D Honeybee (\u003cem\u003eApis mellifera ligustica\u003c/em\u003e) drone embryo proteomes. \u003cem\u003eJ. Insect. Physiol.\u003c/em\u003e \u003cb\u003e57\u003c/b\u003e, 372\u0026ndash;384 (2011).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eOppenheim, S. J., Baker, R. H. \u0026amp; Simon, S. DeSalle R We can't all be supermodels: the value of comparative transcriptomics to the study of non-model insects. \u003cem\u003eInsect Mol. Biol.\u003c/em\u003e \u003cb\u003e24\u003c/b\u003e, 139\u0026ndash;154 (2015).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003ePires, C. V., Freitas, F. C., Cristino, A. S. \u0026amp; Dearden, P. K. Simoes ZL Transcriptome Analysis of Honeybee (\u003cem\u003eApis mellifera\u003c/em\u003e) Haploid and Diploid Embryos Reveals Early Zygotic Transcription during Cleavage. \u003cem\u003ePLoS One\u003c/em\u003e. \u003cb\u003e11\u003c/b\u003e, e0146447 (2016).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eSimon, S. et al. DeSalle R Comparative transcriptomics reveal developmental turning points during embryogenesis of a hemimetabolous insect, the damselfly Ischnura elegans. \u003cem\u003eSci. Rep.\u003c/em\u003e \u003cb\u003e7\u003c/b\u003e, 13547 (2017).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eSandler, J. E. \u0026amp; Stathopoulos, A. Quantitative Single-Embryo Profile of Drosophila Genome Activation and the Dorsal\u0026ndash;Ventral Patterning Network. \u003cem\u003eGenetics\u003c/em\u003e \u003cb\u003e202\u003c/b\u003e (4), 1575\u0026ndash;1584 (2016).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eBando, T. et al. Mito T Analysis of RNA-Seq data reveals involvement of JAK/STAT signalling during leg regeneration in the cricket \u003cem\u003eGryllus bimaculatus\u003c/em\u003e. \u003cem\u003eDevelopment\u003c/em\u003e \u003cb\u003e140\u003c/b\u003e, 959\u0026ndash;964 (2013).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eLi, L-T., Zhu, Y-B., Ma, J-F. \u0026amp; Li, Z-Y. Dong Z-P An Analysis of the \u003cem\u003eAthetis lepigone\u003c/em\u003e Transcriptome from Four Developmental Stages. \u003cem\u003ePLoS One\u003c/em\u003e. \u003cb\u003e8\u003c/b\u003e, e73911 (2013).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003ePrid\u0026ouml;hl, F. et al. Transcriptome sequencing reveals maelstrom as a novel target gene of the terminal system in the red flour beetle Tribolium castaneum. \u003cem\u003eDevelopment\u003c/em\u003e \u003cb\u003e144\u003c/b\u003e (7), 1339\u0026ndash;1349 (2017).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eGinzburg, N. \u0026amp; Cohen, M. \u0026amp; Chipman AD Factors involved in early polarization of the anterior-posterior axis in the milkweed bug \u003cem\u003eOncopeltus fasciatus\u003c/em\u003e. genesis 55, (2017).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eFan, X-B. et al. An Overview of Embryogenesis: External Morphology and Transcriptome Profiling in the Hemipteran Insect Nilaparvata lugens. \u003cem\u003eFront. Physiol.\u003c/em\u003e \u003cb\u003e11\u003c/b\u003e, (2020).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eWang, Z. L. et al. A Chromosome-Scale Assembly of the Asian Honeybee \u003cem\u003eApis cerana\u003c/em\u003e Genome. \u003cem\u003eFront. Genet.\u003c/em\u003e \u003cb\u003e11\u003c/b\u003e, 279 (2020).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eWang, Y. et al. Integrative machine learning models predict prostate cancer diagnosis and biochemical recurrence risk: Advancing precision oncology. \u003cem\u003enpj Digit. Med.\u003c/em\u003e \u003cb\u003e8\u003c/b\u003e, 524 (2025).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eZhu, Z. et al. 2025 Exploration of differential expression and biological significance of amino acid metabolism genes in osteoarthritis. \u003cem\u003eFront. Immunol. Volume\u003c/em\u003e \u003cb\u003e16\u003c/b\u003e, (2025).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003ede Sena Brandine, G. \u0026amp; Smith, A. D. Falco: high-speed FastQC emulation for quality control of sequencing data. \u003cem\u003eF1000Res Nov\u003c/em\u003e. \u003cb\u003e7\u003c/b\u003e, 8:1874 (2019).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eKim, D., Langmead, B. \u0026amp; Salzberg, S. HISAT: a fast spliced aligner with low memory requirements. \u003cem\u003eNat Methods\u003c/em\u003e (2015).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eLi, H. et al. Genome Project Data Processing S: The Sequence Alignment/Map format and SAMtools. \u003cem\u003eBioinformatics\u003c/em\u003e \u003cb\u003e25\u003c/b\u003e (16), 2078\u0026ndash;2079 (2009).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eLiao, Y., Smyth, G. K. \u0026amp; Shi, W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. \u003cem\u003eBioinformatics\u003c/em\u003e \u003cb\u003e30\u003c/b\u003e (7), 923\u0026ndash;930 (2014).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eHaas, B. J., Papanicolaou, A. \u0026amp; Yassour, M. A De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. \u003cem\u003eNat. Protoc.\u003c/em\u003e \u003cb\u003e8\u003c/b\u003e, 1494\u0026ndash;1512 (2013).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eSun, J. et al. TOR signaling, epigenetics, and tumor development. \u003cem\u003eFront. Genet.\u003c/em\u003e \u003cb\u003e13\u003c/b\u003e, (2022).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eZhang, T. et al. Drivers of protein phosphorylation (multi-omics). \u003cem\u003eGenome Biol.\u003c/em\u003e \u003cb\u003e24\u003c/b\u003e (1), 52 (2023).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eBartel, D. P. MicroRNAs:genomics,biogenesis, mechanism, and function. \u003cem\u003eCell\u003c/em\u003e \u003cb\u003e116\u003c/b\u003e (2), 281\u0026ndash;297 (2004).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eMalone, C. D. Hannon GJ Small RNAs as Guardians of the Genome. \u003cem\u003eCell\u003c/em\u003e \u003cb\u003e136\u003c/b\u003e (4), 656\u0026ndash;668 (2009).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eFu, X-D. \u0026amp; Non-coding RNA: a new frontier. \u003cem\u003eNatl. Sci. Rev.\u003c/em\u003e \u003cb\u003e1\u003c/b\u003e, 190\u0026ndash;204 (2014).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eQu, Z. \u0026amp; Adelson, D. L. Conservation and roles of ncRNA. \u003cem\u003eFront. Genet.\u003c/em\u003e \u003cb\u003e3\u003c/b\u003e. (2012).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eCalixto, C. P. G. et al. Rapid, dynamic AS in Arabidopsis cold response. \u003cem\u003ePlant. Cell.\u003c/em\u003e \u003cb\u003e30\u003c/b\u003e, 1424\u0026ndash;1444 (2018).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eZhang, H. et al. Population-level AS and agronomic traits of rice. Plant Cel.(2024) l.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eZhong, Y. et al. Plant.Pan-transcriptomic AS and cold tolerance in rice. (2024).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eNusse, R. \u0026amp; Clevers, H. Wnt/β-Catenin Signaling, Disease, and Emerging Therapeutic Modalities. \u003cem\u003eCell\u003c/em\u003e \u003cb\u003e169\u003c/b\u003e (6), 985\u0026ndash;999 (2017).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eTripathi, B. K. \u0026amp; Irvine, K. D. The wing imaginal disc. Genetics. Apr 4;220(4):iyac020. PMID: 35243513; PMCID: PMC8982031. (2022).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eZhang, P. \u0026amp; Edgar, B. A. The Kr\u0026uuml;ppel-like factor Cabut has cell cycle regulatory properties similar to E2F1.118,( (2025).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eBlanco, E. et al. Gene expression following induction of regeneration in Drosophila wing imaginal discs. Expression profile of regenerating wing discs(2010).\u003c/span\u003e\u003c/li\u003e\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":false,"hideJournal":true,"highlight":"","institution":"","isAcceptedByJournal":false,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true},"keywords":"Apis cerana, drones, embryonic, deep learning, ResNet","lastPublishedDoi":"10.21203/rs.3.rs-7939377/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-7939377/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003eStage-resolved regulators in Apis cerana drone embryos remain unclear. Here, RNA-seq sampled at 12/24/36/48/60 h was analyzed by WGCNA to build co-expression modules, refined with k-nearest neighbors on five-point trajectories, and benchmarked across 19 ML/DL classifiers (including ResNet,SVM,ExtraTrees,KNN and RandomForest). To balance accuracy and interpretability, we combined model benchmarking, consensus voting, and feature attribution (e.g., gradient/SHAP) and report parameter and sampling sensitivities. Expression co-regulated into 6 clusters aligned with developmental stages\u0026mdash;germline/epigenetic control (12 h), ncRNA/rRNA-centric regulation (24 h), neurogenesis/morphogenesis (36\u0026ndash;48 h), and organ specialization (60 h). ResNet achieved\u0026thinsp;~\u0026thinsp;97% accuracy for gene-pattern classification; consensus reduced\u0026thinsp;\u0026gt;\u0026thinsp;5,800 candidates to ~\u0026thinsp;2,000 high-confidence stage-specific genes and hub PPI subnetworks. The WGCNA\u0026rarr;kNN(unsupervised cluster analysis) \u0026rarr;ML/DL pipeline thus resolves stage-specific regulators while addressing the \u0026ldquo;black-box\u0026rdquo; trade-off, and its performance generalizes with attention to module granularity (β, minModuleSize), k, model capacity/regularization, and temporal sampling density.\u003c/p\u003e","manuscriptTitle":"Integrative co-expression and deep learning reveal stage-specific regulators of Apis cerana drone embryogenesis","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2025-11-25 06:08:58","doi":"10.21203/rs.3.rs-7939377/v1","editorialEvents":[{"type":"communityComments","content":0}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"c8a80b71-d94a-42c6-a6d2-59726ce8da08","owner":[],"postedDate":"November 25th, 2025","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"posted","subjectAreas":[{"id":58523734,"name":"Biological sciences/Computational biology and bioinformatics"},{"id":58523735,"name":"Biological sciences/Developmental biology"}],"tags":[],"updatedAt":"2025-12-15T07:10:25+00:00","versionOfRecord":[],"versionCreatedAt":"2025-11-25 06:08:58","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-7939377","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-7939377","identity":"rs-7939377","version":["v1"]},"buildId":"8U1c8b4HqxoKbykW_rLl7","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}