The analyses in this study are based on publicly available sequence data from the NCBI GenBank database under accession number SRP131292 (https://www.ncbi.nlm.nih.gov/sra/?term=SRP131292) (Chen et al., 2016). The datasets used and analyzed were obtained from this publication and can be requested from the corresponding author of the original study upon reasonable request.

Before alignment, the raw data were quality-checked and cleaned using FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Specifically, we used the default parameters of FastQC to assess quality scores (Q values), GC content, and the proportion of N bases in the reads. The filtered reads were then aligned to the reference genome (GCA_029169275.1) using Hisat2 (v2.2.1) with the following parameters: -max-intronlen 10,000 (maximum intron length set to 10,000 bp), --dta (optimization for downstream transcript assembly), and--phred33 (Phred+33 quality score encoding) (Wang et al., 2021a). The resulting SAM files were converted to BAM format, sorted, and indexed with Samtools (v1.6) (Ma et al., 2022). Quantification of aligned reads was performed using FeatureCounts (Patro et al., 2017). Finally, differential expression analysis was conducted using Trinity (v2.15.1) with the DESeq2 algorithm (Rogers et al., 2021).

Using emapper. py (v2.1.10) with the eggNOG database (Roth et al., 2023), GO and KEGG annotations were performed on the protein sequence of the longest csd sequence in the reference genome of A. cerana. The annotated file was then converted into org. db format to facilitate GO and KEGG enrichment analysis using the clusterProfiler (v4.0) package in R (Setton and Sharma, 2018).

Time-series transcriptome analysis was conducted on significant genes in the expression matrix using the ClusterGVis package in R (https://github.com/junjunlabClusterGVis). By clustering similar genes, gene modules in different clusters were analyzed to understand gene expression patterns across various embryonic stages. The ClusterGVis package integrates fuzzy C-means (Mfuzz) and K-means clustering algorithms, allowing for efficient clustering analysis of gene expression matrices and providing insights into gene expression trends at each time point.

For alternative splicing analysis, filtered sequencing reads were aligned to the reference genome (GCA_029169275.1) using the Salmon software, which provides a fast and accurate tool for quantifying transcript expression. Salmon aligns the reads pseudo-selectively, allowing for quantification without the need to generate large alignment files. This approach is efficient in handling complex transcriptomes, especially when exploring alternative splicing events across multiple samples (Slabaugh et al., 2019). Next, alternative splicing events were identified using the SUPPA (Simulation of Uncertainty and Prediction of Alternative splicing) tool (Song et al., 2020; Tao et al., 2021). SUPPA leverages transcript-level quantification data from Salmon and annotation files (e.g., GTF format) to generate alternative splicing events, such as exon skipping, alternative 5′ and 3′ splice sites, and intron retention. By applying SUPPA’s splicing event generation functionality (suppa.py), the annotated files were processed to identify splicing variations across conditions, providing insights into how alternative splicing patterns may vary between developmental stages or treatment groups (Song et al., 2020; Tao et al., 2021).

The Illumina RNA sequencing generated an average of 86.89 million (M) raw reads per sample. After filtering, an average of 86.59 million clean reads per sample was retained for subsequent transcriptome analysis. The average quality scores for clean reads showed Q20 and Q30 ratios of 100%, indicating high quality of the obtained clean reads. Additionally, the average GC content was 38.24% (Supplementary Table S1). We used HISAT2 software to perform alignment analysis with the reference genome. The sliding window density analysis showed that over 80% of the transcripts were uniquely aligned to the reference genome, with an average alignment rate of 83.33% across nine samples. Only 15.66% of the reads (7.97 million) failed to align to the genome (Supplementary Table S2).

To assess the consistency of sample collection and explore the transcriptomic relationships among embryos at 24, 48, and 72 h, we conducted hierarchical clustering and principal component analysis (PCA) on the global gene expression data of these samples. These two analytical methods summarize the similarities and differences among the samples. Prior to these analyses, we normalized the gene expression levels using the regularized log transformation method implemented in the limma package. The results showed that the 48-hour and 72-hour samples were more similar to each other, while both time points exhibited significant differences compared to the 24-hour samples, consistent with previous studies (Figure 1a). In the PCA plot (Figure 1b), embryo samples collected at the same time point clustered together, further validating the results of the hierarchical clustering analysis. Notably, the 24-hour samples were distinctly separated from the 48-hour and 72-hour samples along the first principal component (PC1), indicating that PC1 primarily reflects the gene expression differences between the 24-hour stage and subsequent developmental stages. Meanwhile, the 48-hour and 72-hour samples were differentiated along the second principal component (PC2), suggesting additional transcriptomic differences between these two time points. Moreover, we found that the Apd−2 and HEX70b gene had high loadings on both PC1 and PC2 (Supplementary Figure S1), indicating its significant role in distinguishing samples from different developmental stages. Both analyses demonstrate that gene expression at 24 h of embryonic development was significantly different from the subsequent two stages, suggesting that the most substantial changes and regulation in gene expression occur between 24 and 48 h.

The volcano plots (Figure 2b) visually illustrate the overall trends in gene expression levels across stages. Each point represents a gene, with the horizontal axis showing the logarithmic fold change in gene expression (log2 Fold Change) and the vertical axis showing the negative logarithm of the significance level (−log10 p-value). Significantly upregulated genes were distributed in the upper right of the plot, while significantly downregulated genes appear in the upper left. It is evident that the number of differentially expressed genes is highest during the 24-hour to 48-hour stage; the volcano plot displayed a more dispersed pattern, indicating the most dramatic changes in gene expression occur at this stage. This was consistent with the previous hierarchical clustering and principal component analysis results, further confirming that the period from 24 to 48 h was a critical turning point in embryonic development, during which significant reprogramming of gene expression occurs.

In contrast, from 48 h to 72 h, the number of significantly differentially expressed genes in the volcano plot decreases markedly, and the plot becomes more concentrated. This suggested that the degree of gene expression changes was reduced during this stage, and the embryonic development process enters a relatively stable period.

The genes showing significant differential expression in these volcano plots may be closely related to important physiological and biochemical processes, morphological formation, and functional differentiation at each embryonic stage. Upregulated genes may be involved in cell proliferation, differentiation, and the acquisition of specific functions, while downregulated genes may be associated with the regulation of the cell cycle and adjustments in metabolic activities. Further functional annotation and pathway enrichment analyses of these differentially expressed genes will help reveal their specific roles and regulatory mechanisms in embryonic development, providing important insights into the molecular basis of embryogenesis.

To further explore the dynamic changes in gene expression during embryonic development, we employed a k-means clustering algorithm to partition 9,492 genes into 10 clusters based on the similarity of their expression patterns. Detailed analysis of these clusters revealed significant enrichment of genes at specific developmental stages (Figure 3).

At the 24-hour embryonic stage, genes in clusters C3, C5, and C6 exhibited significant enrichment. The associated Gene Ontology (GO) terms mainly include fatty acid derivative metabolic processes, biosynthetic processes of organic acids and carboxylic acids, key complexes related to transcription and translation, mRNA metabolism and splicing, and the metabolic processing of rRNA and ncRNA (Supplementary Tables S3, S4, S5). These results indicated that at the 24-hour stage, gene expression primarily involves fundamental metabolic processes and the establishment of gene expression machinery, emphasizing the construction of basic cellular components and functions required for early development.

Overall, these genes play crucial roles in multiple key processes and functions of cell biology and developmental biology throughout various stages of embryonic development. As development progresses into the 48-hour stage, genes in clusters C7 and C8 showed high expression levels. The related GO terms were concentrated on the positive regulation of biological and cellular processes, post-embryonic development, larval or pupal development, appendage morphogenesis, and RNA synthesis, modification, and processing (Supplementary Tables S6, S7). This suggested that at the 48-hour stage, genes were involved in developmental programs that initiated morphological changes and cell differentiation, promoting the formation of specific structures while continuing RNA processing activities. Further, at the 72-hour stage, genes in clusters C9 and C10 become predominant. The associated GO terms mainly pertain to protein synthesis processes, the structure and function of ribosomes, appendage development and morphogenesis (including imaginal disc-derived appendage development, appendage morphogenesis, and leg disc development), and specific molecular binding and enzymatic activities (such as signal receptor binding, heme binding, tetrapyrrole binding, and nucleotide monophosphate kinase activity) (Supplementary Tables S8, S9).

These resulted reflect that at the 72-hour stage, gene expression focuses on protein synthesis and the role of ribosomes as the core machinery for protein production, promoting further development of appendages and demonstrating diverse functions in molecular binding and catalysis, indicating an advanced stage of tissue formation and functional specialization.

Alternative splicing analysis validated the expression levels of different gene transcripts at three time points during embryonic development. The results showed that the numbers of alternatively spliced genes and splicing events at these three time points significantly overlapped with the distribution of differentially expressed genes. During the critical developmental stage from 24 to 48 h, the counts of all seven types of splicing events reached their highest levels, with Alternative First Exon (AF) events being the most abundant, totaling 255 (Figures 4a,b).

The high frequency of AF events may be related to the complexity and diversity of gene expression regulation in the early stages of embryonic development. Between 24 and 48 h a pivotal developmental period embryonic cells undergo rapid proliferation and differentiation, necessitating precise and flexible mechanisms of gene expression regulation. AF events, by selecting different initial exons, can produce protein isoforms with varying N-terminal structures, thereby influencing protein localization, function, and interactions with other molecules. This mechanism allowed the same gene to perform diverse functions in different cell types or developmental stages.

The results of GO enrichment analysis (Figure 5) indicate that significantly differentially spliced transcripts in embryos are primarily involved in biological processes such as morphogenesis, growth, cytoskeleton organization, and the formation of specific organs and appendage structures. Specifically, these functions encompass biological processes such as cytoskeleton organization, morphogenesis, neuromuscular junction formation, axon development, cell polarity, protein polymerization, cell junction organization, organ development, and cell structure formation (Supplementary Table S10). These processes include core functions such as cell differentiation, signal transduction, structural assembly, and cell-cell junctions, revealing the role of complex regulatory networks in maintaining and modulating the morphology and function of cells and tissues during development.

These findings suggest that alternative splicing events may play an essential role in regulating the formation and maintenance of specific cellular structures at critical stages of embryonic development. For instance, the dynamic regulation of cytoskeleton and cell junctions is crucial for maintaining cell morphology, differentiation, and migration, while regulation of protein polymerization and axon development may be pivotal in the formation of the nervous system. Additionally, these differentially spliced transcripts may have significant roles in the morphogenesis and physiological functions of specific organs and appendages (e.g., muscle and nervous system), highlighting the highly refined gene expression regulatory mechanisms during embryonic development. These insights provide important clues for further exploration of how alternative splicing mediates functional diversity of genes and precise regulation in specific developmental processes.

To investigate the key genes involved in alternative splicing events during embryonic development, we selected three genes with the highest connectivity in GO enrichment analysis: Dll, CaMKII, and LOC408563. Our results showed that the Dll gene exhibited exon-skipping events at both the 24-h and 48-h time points, with each time point featuring unique exon skipping events. Additionally, at 72 h, an additional exon skipping event was observed compared to the 48-h time point (Figure 6).

In the CaMKII gene, compared to 48 h, the 24-hour time point had two additional exon skipping events and one intron retention event. At 72 h, CaMKII showed one unique intron retention event but had three fewer exon skipping events and one fewer intron retention event. The LOC408563 gene exhibited one more exon skipping event and one more intron retention event at 24 h compared to 48 h. Conversely, at 48 h compared to 24 h, it had one unique intron retention event. Compared to 72 h, the 48-hour time point had two more exon skipping and one more intron retention event, while at 72 h, a unique intron retention event was present (Figure 6).

The changes in these splicing events at different developmental stages may suggest distinct roles for these genes in regulating embryonic development. Time-specific splicing patterns likely relate to the precise regulatory requirements of gene expression, with such dynamic splicing mechanisms enabling the production of different protein isoforms to meet the varying functional needs of cells during development. Thus, these alternative splicing events could represent an important aspect of the regulatory complexity in embryonic development.

To further validate the reliability of our public database-based analysis results, we performed independent verification of the identified differentially spliced genes using quantitative real-time PCR (qRT-PCR). Specifically, we selected three representative genes, Dll (XM_006560304, XM_006560302.3), CaMKII (XM_00650514.3), and LOC408563 (XM_006563138.3), and validated the transcripts corresponding to the relatively unique splicing events at the 24-hour, 48-hour, and 72-hour time points. The results demonstrated that the alternative splicing patterns of these three genes were highly consistent with those predicted by our bioinformatics analysis, confirming the accuracy and reliability of our analytical pipeline (Figure 7; Supplementary Table S11).