Dynamic Landscape of Alternative Splicing During Early Embryogenesis of the Rhesus Monkey

TL;DR

This study maps how alternative splicing changes during early development in rhesus monkeys, offering insights into human embryogenesis.

Contribution

The first systematic analysis of alternative splicing dynamics in early rhesus macaque embryogenesis.

Findings

Identified previously unannotated zygotic genome activation (ZGA) genes in rhesus macaque embryos.

Observed a significant increase in alternative splicing events from the 8-cell to morula stages.

Found 35 key splicing factors and potential regulatory relationships with specific splicing events.

Abstract

The rhesus macaque is one of the closest evolutionary relatives to humans, making the study of alternative splicing (AS) during its early embryonic development highly valuable for understanding human embryogenesis and related diseases. However, systematic studies in this context remain limited. Here, a comprehensive bioinformatic analysis of AS was performed using RNA-seq data spanning early rhesus macaque embryogenesis. We identified multiple previously unannotated zygotic genome activation (ZGA) genes, thereby refining the rhesus macaque ZGA gene repertoire. The landscape of AS and differential AS events (DASEs) across early stages was characterized, revealing dynamic and stage-specific regulation, with a marked increase in AS events from the 8-cell to morula stages. In addition, weighted gene co-expression network analysis identified 35 key splicing factors (SFs) involved in regulating early rhesus macaque embryonic development. Finally, we calculated the correlation between differentially expressed SFs and DASEs during the ZGA process, and identified potential regulatory relationships between several SFs (TRA2B, IGF2BP1, HNRNPAB, and MATR3) and specific DASEs. Collectively, this study provides the first systematic analysis of AS dynamics and regulation in early rhesus macaque embryogenesis, highlighting its critical role in development and offering a valuable reference for understanding AS in early human embryos.

Figures

Full text

1. Introduction

Early embryonic development in mammals is a complex process [1]. It begins with the zygote formed by the fusion of sperm and oocyte [2]. The zygote subsequently undergoes successive rounds of cell fate specification and extensive epigenetic reprogramming, including zygotic genome activation (ZGA) [3] ultimately developing into a mature blastocyst accompanied by dissolution of the surrounding zona pellucida [4]. Following blastocyst formation, implantation into the uterine endometrium is initiated [5]. As a pivotal event in early embryonic development, ZGA plays a fundamental role in maintaining embryonic totipotency and initiating developmental programs [6]. Notably, ZGA exhibits pronounced species specificity. In mice, major ZGA occurs at the 2-cell (2C) to 4-cell (4C) stage. In humans, it takes place between the 4C and 8C stages. In rhesus macaque, major ZGA is initiated from the 8C to the morula stage [7,8,9].

Alternative splicing (AS) is a post-transcriptional regulatory process through which introns are removed and exons are differentially joined to generate mature mRNA isoforms [10]. This process greatly enhances the diversity of the transcriptome and proteome [11]. Canonical AS events include skipped exon (SE), mutually exclusive exons (MXE), retained intron (RI), alternative 5′ splice site (A5SS), alternative 3′ splice site (A3SS), alternative first exon (AFE), and alternative last exon (ALE) [12] (Figure S1). AS is pervasive during mammalian embryonic development [13] and is implicated in approximately 15% of human genetic diseases [14,15,16,17]. Owing to strict ethical constraints, studies of human embryonic development are severely limited. Therefore, the use of model organisms such as mice and rhesus macaques represents a critical strategy for investigating early developmental mechanisms.

The rhesus macaque diverged from the human ancestral lineage approximately 25 million years ago and shares ~93% genetic homology with humans [18]. This high degree of evolutionary conservation confers unique advantages for modeling human developmental processes, disease pathogenesis, and gene regulatory mechanisms, making the rhesus macaque an especially relevant system for studying AS during early embryonic development [19]. However, AS during preimplantation embryonic development in the rhesus macaque has rarely been studied. Therefore, this study leveraged RNA-seq data spanning eight consecutive developmental stages from oocytes to blastocysts to construct a comprehensive AS landscape and identified 4542 differentially AS events (DASEs) across 2094 genes. Subsequently, a total of 638 splicing factors (SFs) were curated in the rhesus macaque. Using weighted gene co-expression network analysis (WGCNA), 35 core SFs were identified. Correlation analyses further revealed the regulation of these SFs in splicing modulation during embryonic development. Overall, this study systematically analyzes the molecular regulatory mechanisms of cell fate transitions during early embryogenesis through the dynamic features of AS.

2. Results

2.1. Global Landscapes of Gene Expression and AS During Early Embryonic Development in Rhesus Macaque

To delineate the global dynamics of gene expression and AS during early embryonic development in rhesus macaque, we analyzed RNA-seq data generated by paired-end sequencing on the Illumina platform, spanning eight developmental stages from the oocytes (GV) to the blastocyst (BLA).

In total, 21,585 protein-coding genes were detected across all samples (Figure 1a), of which constitutively expressed genes accounted for approximately 61% (12,302 genes). To assess the stage-specific transcriptional activity during embryogenesis, protein-coding genes were assigned to individual developmental stages (Figure 1b). The number of expressed genes decreased from 16,817 at the GV stage to 14,791 at the metaphase II oocytes (MII) stage, followed by a progressive increase that peaked at the 8-cell (8C) stage with 18,274 genes. Subsequently, gene numbers declined sharply to 16,595 at the MOR stage, and increased again at the blastocyst stage, reaching 17,651 genes. These dynamic patterns indicate that early rhesus macaque embryogenesis is governed by tightly orchestrated temporal regulation of gene expression.

Using rMATS, we systematically quantified AS events across preimplantation developmental stages (Table 1). Approximately 40% (8488) of protein-coding genes underwent AS during early development. Among the five major AS types, genes harboring SE events exhibited the most pronounced stage-dynamic changes. Overall, the number of AS-associated genes remained relatively low from the 1-cell (1C) to 4C stages, followed by a marked increase beginning at the 8C stage and reaching a maximum at the blastocyst stage (Figure 1c). These results suggest that the activation of AS is closely coupled to developmental progression.

To ensure the robustness of AS detection, we further filtered AS events by retaining only those consistently identified in at least 60% of biological replicates at each stage, and all subsequent analyses were based on the filtered splicing events. The resulting distribution of filtered AS events mirrored the trends observed prior to filtering, indicating that the observed developmental dynamics of AS are highly reproducible and not driven by AS sparse (Figure 1c,d).

To further explore global developmental trajectories, principal component analysis (PCA) was performed based on gene expression levels (TPM) and AS levels (IncLevel) (Figure 1e). PCA of gene expression profiles revealed clear stage-specific transitions across early embryonic development. Notably, PCA based on IncLevel values produced trajectories that were highly concordant with those derived from gene expression data, underscoring a tight coupling between transcriptional and post-transcriptional regulation. Together, these results highlight AS as an integral regulatory layer that coordinates with gene expression programs to ensure orderly progression of early embryogenesis.

2.2. Analysis of DEGs and DASEs During Early Embryonic Development in Rhesus Macaque

To more clearly delineate the dynamic changes in gene expression and AS during early embryonic development in rhesus macaque, DESeq2 and rMATS were used to identify differentially expressed genes (DEGs) and DASEs between adjacent developmental stages, respectively. Across all stages, a total of 10,837 DEGs and 4542 DASEs were detected. These DASEs occurred in 2094 genes, hereafter referred to as genes with DASE. Notably, 89% of genes with DASE exhibited differential expression (Figure 2a), indicating a strong association between differential gene expression and AS during early embryogenesis. This convergence suggests that gene regulation and splicing programs are coordinated during embryonic genome activation and cell fate determination.

Stage-wise comparison revealed a pronounced increase in both DEGs and DASEs during the transition from the 8C to the MOR stage (Figure 2b), coinciding with the major wave of ZGA period. This sharp increase indicates that large-scale transcriptional activation is accompanied by extensive remodeling of the splicing landscape, underscoring coordinated regulation at both transcriptional and post-transcriptional levels during this critical developmental window. Among the five major types of AS events, SE events predominated across all developmental stages, consistently accounting for at least 65% of total DASEs (Figure 2c). Although SE events were dominant, the relative proportions of distinct splicing types varied across developmental stages, indicating pronounced temporal regulation of AS patterns.

To explore the functional relevance of developmentally regulated splicing, we performed Gene Ontology (GO) enrichment analysis on genes undergoing differential SE events at each stage. Enrichment biological processes were closely associated with early embryonic development (Figure 2d), including translation, cell division, autophagy, methylation, and chromatin remodeling [16,20,21,22,23,24,25]. Notably, splicing-related biological processes were specifically and significantly enriched during the 8C-MOR transition, including “RNA splicing”, “mRNA processing”, and “mRNA splicing, via spliceosome,” whereas no significant enrichment of these terms was observed at other stages. This stage-restricted enrichment pattern is consistent with previous observations [26] and further supports the notion that large-scale activation of AS is a defining molecular feature of the ZGA stage.

ZGA represents a developmental transition in which control of embryonic development shifts from maternally deposited transcripts to autonomous transcription from the zygotic genome [27,28]. To delineate the contribution of AS to this process, we identified genes significantly upregulated during the 8C-MOR transition and defined a ZGA gene set comprising 3618 genes (Table S1). This set includes well-established ZGA regulators reported in humans and mice, such as DUXB and GATA6 [29,30], supporting the reliability of the dataset. Notably, 608 ZGA genes also underwent differential AS (Figure 2e), suggesting that splicing regulation is extensively integrated into the ZGA transcriptional program. Several ZGA-associated genes with DASE have been implicated in embryonic development and genome regulation. For instance, POU5F1 plays a critical role in dorsoventral patterning during zebrafish embryogenesis by regulating VOX and FGF8A expression [31]. In addition, NLRP7 has been shown to modulate AS of homologous recombination-related genes in human embryonic stem cells through interactions with RNA splicing and the DNA damage response factors, and its dysfunction compromises genomic integrity [32]. Together, these observations indicate that coordinated transcriptional activation and AS of ZGA genes constitute an essential regulatory axis underlying early embryonic development in rhesus macaque.

2.3. Identification of Key SFs During ZGA in Rhesus Macaque

To identify core regulators underlying the extensive transcriptional and splicing reprogramming during ZGA, we performed WGCNA using DEGs identified during the 8C to morula (8C-MOR) transition. Using a soft-thresholding power of 8, a gene co-expression network was constructed, yielding 15 distinct co-expression modules (Figure 3a). Each module was assigned a unique color and contained 39 to 2449 genes. The turquoise module was the largest, comprising 2449 genes, whereas the midnightblue module was the smallest, containing 39 genes. Genes that could not be confidently assigned to any module were grouped into the grey module. Hub-genes within each module were subsequently identified based on high module eigengene connectivity (|KME| ≥ 0.8) and strong intramodular connectivity (kWithin ≥ 0.5). Notably, hub-genes from the brown, green, pink, and black modules were almost exclusively composed of ZGA genes (Figure 3b). Therefore, these four modules are referred to as ZGA-associated modules. Moreover, the majority of genes within these four modules exhibited coordinated upregulation at the MOR stage (Figure 3c), indicating that these modules represent ZGA-associated transcriptional programs that are synchronously activated during the 8C-MOR transition.

To investigate the biological functions of these ZGA-associated modules, we performed GO enrichment analysis on the genes in the four ZGA-associated modules. All four modules showed significant enrichment for AS-related biological processes, including mRNA processing, RNA splicing, and mRNA splicing via spliceosome (Figure 3d), suggesting that splicing regulation is a prominent and coordinated feature of the ZGA transcriptional network.

To further pinpoint key splicing regulators, we intersected hub-genes from the four ZGA-associated modules with the curated SF dataset. This analysis identified 35 SFs as hub genes (Table 2). These SFs included members of the heterogeneous nuclear ribonucleoprotein (hnRNP) family (e.g., HNRNPA3, HNRNPAB, and HNRNPD) and serine/arginine-rich (SR) protein family (e.g., SRSF3 and TRA2B) [33], as well as several SFs implicated in early embryonic development, such as SNRPD and SNRPB2 [34]. Collectively, these findings indicate that a defined set of SFs occupies central positions within ZGA-associated co-expression networks and may serve as core regulators of AS programs during early rhesus macaque embryogenesis.

2.4. Expression Dynamics and Regulatory Roles of SFs During the Major ZGA Period

To systematically interrogate the regulatory roles of SFs in AS during the ZGA, we focused on the 8C-MOR transition. A total of 229 differentially expressed SFs (DESFs) were identified at this stage. We assessed the association between SF expression levels (TPM) and splicing changes by correlating DESF expression with IncLevel values of SE type DASEs. Correlation analysis revealed extensive coupling between SF expression and SE splicing dynamics (Figure 4a). Specially, 137 DESFs exhibited strong positive correlations with differential SE events (r > 0.8), whereas 151 DESFs showed strong negative correlations (r < −0.8), indicating that SFs participate in both activating and repressing exon inclusion during ZGA.

In addition, we focused on the changes in the expression levels of DESFs and the changes in splicing levels during embryonic development (Figure 4b). The results showed that during the 8C-MOR stage, 115 out of the 137 DESFs showed an increase (Figure 4b, upper panel), and the splicing levels of the corresponding 69 differential SE events also increase (Figure 4b, low panel). Similarly, during the 8C-MOR stage, 22 out of 137 DESFs showed a decrease (Figure 4b, upper panel), and the splicing levels of the corresponding 31 SE events also decreased (Figure 4b, low panel). Genes harboring these differential SE events were significantly enriched for spliceosome-mediated mRNA splicing-related biological processes (Figure 4c), suggesting that coordinated upregulation of SFs reinforces global activation of the splicing machinery during ZGA. Next, we focused on 35 core DESFs identified as hub genes by WGCNA. Notably, 24 of these hub SFs displayed strong positive correlations with SE-type DASEs (r > 0.8, Figure 4a) and all 24 were classified as ZGA genes (Figure 4d). These observations indicate that a subset of ZGA-activated SFs occupies central regulatory positions in the splicing network and may exert stage-specific post-transcriptional control during early embryonic development.

SFs regulate AS by recognizing sequence-specific RNA-binding motifs (RBMs) on pre-mRNAs [34]. To explore potential direct regulatory interactions, we retrieved RBMs for the 24 highly correlated hub SFs from the CISBP-RNA database. Eight high-confidence RBMs were identified, including those for TRA2B and IGF2BP1 (Figure 5a,d, Table S2). TRA2B, a member of the SR protein family, has been implicated in early embryonic development [33]. During the ZGA stage, the splicing levels of development-related genes such as PCBP2 and ANLN [35,36] are significantly upregulated and show a positive correlation with the expression level of TRA2B (log_2_FC = 2.2) (Figure 5c). Consistently, FIMO analysis identified TRA2B binding sites located 76 bp upstream of exon 12 of PCBP2 and within exon 21 of ANLN (Figure 5c,b), supporting a direct association between TRA2B binding and exon inclusion in these genes. Similarly, IGF2BP1 exhibited pronounced upregulation during ZGA (log_2_FC = 4.2). Motif analysis identified a putative IGF2BP1 binding site within exon 28 of DOT1L (Figure 5e,f). DOT1L encodes an epigenetic regulator responsible for H3K79 methylation and plays a critical role in transcriptional activation and embryonic viability [37]. These findings suggest that IGF2BP1 may modulate the splicing of chromatin regulators to fine-tune transcriptional programs during ZGA.

In addition, the ZGA gene KDM5B not only underwent SE events but was also significantly enriched in methylation-related biological processes [38] (Figure 4c). KDM5B expression showed a strong positive correlation with the splicing factor SRSF2, and a putative SRSF2 binding site was detected within exon 23 of KDM5B (Figure S2a–c). Given the established role of KDM5B in H3K4 demethylation and ZGA initiation, these observations suggest a potential mechanism by which splicing regulation intersects with epigenetic remodeling during early embryogenesis.

In parallel, 151 DESFs exhibited strong negative correlations with DASEs. Among these, 120 were upregulated and 31 were downregulated, corresponding to 58 upregulated and 143 downregulated SE events, respectively (Figure S3a,b). GO enrichment analysis indicated that the genes associated with these differential SE events were significantly enriched in spliceosome-related biological processes, such as mRNA splicing (Figure S3c).

Among the 35 hub DESFs, 25 displayed significant negative correlations with DASEs. RNA-binding motifs were identified for seven of these SFs (Figure S3d, Table S3). For example, HNRNPAB, a core member of the hnRNP family, specifically recognizes pyrimidine-rich sequences within pre-mRNAs or mRNAs via its RNA-binding domains [39], exhibited a binding site immediately downstream of exon 14 of MUTYH. MUTYH plays a crucial role in repairing oxidative DNA damage during early embryogenesis, when reactive oxygen species levels are elevated [40]. In addition, MATR3 binding motifs were detected downstream of exon 4 of FOXR1, a transcription factor essential for neural development [41] (Figure S3d–f). Given MATR3’s known affinity for AU-rich and structured RNA elements [42,43], these findings provide a molecular basis for its involvement in exon repression during ZGA.

Together, these integrative analyses demonstrate that SFs play central and multifaceted regulatory roles during early rhesus macaque embryogenesis. SFs function not only as ZGA-activated genes themselves but also as key modulators of AS in developmentally and epigenetically important genes. By coordinating both positive and negative regulation of exon inclusion, SFs contribute to precise post-transcriptional control during ZGA, thereby shaping the molecular trajectory of early embryonic development.

3. Discussion

In this study, we performed a systematic and comprehensive analysis of AS dynamics during early embryonic development in rhesus macaque, leveraging time-resolved RNA-seq data spanning all preimplantation stages. To our knowledge, this work provides the first global landscape of AS across primate preimplantation development, thereby establishing a foundational resource for dissecting post-transcriptional regulatory mechanisms during early embryogenesis.

Our analyses revealed pronounced temporal remodeling of AS throughout development, which plays a pivotal role in the key events of the maternal-to-zygotic transition. Similar studies have also been conducted in other species, particularly in humans and mice [16]. Principal component analysis based on IncLevel recapitulated the developmental trajectory inferred from gene expression profiles, underscoring the notion that AS constitutes a regulatory layer that is coordinated with, and comparably influential to, transcriptional control. These findings highlight AS as an integral component of the gene regulatory architecture governing early embryonic development.

ZGA represents a pivotal developmental milestone characterized by a shift from maternal to zygotic transcriptional control. Our results demonstrate that a substantial fraction of ZGA genes undergo differential AS, suggesting that AS is tightly coupled to transcriptional activation during this stage. By selectively processing pre-mRNAs into distinct isoforms, AS may enable fine-tuning of protein function, dosage, and interaction networks, thereby contributing to the rapid reorganization of developmental programs required for successful embryonic progression. These observations are consistent with accumulating evidence that AS plays critical roles in cell fate specification, developmental plasticity, and lineage commitment [17].

SFs are central executors of AS regulation, acting through sequence- and structure-specific recognition of pre-mRNAs to modulate spliceosome assembly and exon selection. In the absence of a SF database for rhesus macaques, we constructed a macaque SF database by homology-based conversion of human and mouse SF data. Building upon the AS landscape defined here, we systematically investigated the involvement of SFs in ZGA-associated splicing regulation. Using WGCNA, we identified 35 core SFs occupying hub positions within ZGA-associated modules. Correlation analyses further revealed strong associations between SF expression levels and changes in skipped exon inclusion, accompanied by enrichment of SF binding motifs in the vicinity of differentially spliced exons. Together, these results support a model in which SF-mediated AS regulation is embedded within the core regulatory network of early embryonic development, acting in concert with transcriptional programs to shape developmental trajectories.

This study has several limitations. The conclusions are derived from integrative bioinformatic analyses of high-throughput sequencing data. Although our predictions and models are theoretically well supported at the level of individual AS events, we have not yet elucidated the interaction mechanisms between SFs and RNA, nor have we fully clarified the functional effects mediated by AS isoforms of core genes involved in ZGA. The above analyses are all theoretical inferences based on bioinformatics, and the corresponding experimental validation cannot be carried out in our laboratory at present. In future work, we will collaborate with other laboratories to validate and extend the regulatory relationships identified in this study through the integration of perturbation assays, single-cell long-read sequencing and functional validation approaches.

4. Materials and Methods

4.1. Dataset

4.1.1. RNA-Seq Data Processing

In this study, we collected transcriptomic sequencing data from 26 rhesus macaque samples spanning eight stages of early embryonic development, including GV, MII, 1C, 2C, 4C, 8C, MOR, and BLA, which were obtained from the public dataset GSE86938 [9] (Table S4). Each developmental stage was represented by at least two biological replicates. On average, approximately 10.23 million paired-end reads (125 bp) were generated per sample. Quality control was conducted using FastQC (v0.12.1), adapters and low-quality bases at the 5′ ends were trimmed with fastp (v0.24) [44]. After filtering, each sample retained an average of ~10.04 million clean reads with a mean read length of 110 bp, which were used for all downstream analyses.

4.1.2. Collection of SFs

As no curated SFs database is currently available for the rhesus macaque, SF annotations were inferred through cross-species orthology mapping. A total of 597 human and 598 mouse SFs were collected from published literature. Orthology mapping from human and mouse to rhesus macaque was independently performed using the getLDS function implemented in the R package biomaRt (v2.62.1) [45]. This approach yielded 638 putative rhesus macaque SFs (Table S5). After excluding SFs that were not expressed across all early embryonic developmental stages, a total of 550 SFs were retained for downstream analyses (Table S6).

4.2. Quantification of Gene Expression

Based on the rhesus macaque reference genome annotation (Mmul_10) obtained from Ensembl [46], gene expression quantification was performed with Salmon (v1.10.3) [47] and tximport (v1.34.0) [48]. During quantification, the library type parameter was set to -l A, and the more stringent mapping validation mode (--validateMappings) was enabled to improve quantification accuracy. Tximport generated a gene expression matrix comprising 22,519 genes across 26 samples. After excluding non-coding genes, a final expression matrix containing 21,585 protein-coding genes was obtained for subsequent differential expression analysis.

4.3. Identification of Differentially Expressed Genes

Based on the protein-coding gene expression matrix, DESeq2 (v1.46.0) [49] was used to analyze DEGs between consecutive developmental stages. DEGs were defined with |log_2_ fold change| ≥ 1 and p < 0.05.

4.4. Identification of AS Events

rMATS (v4.3.0) [50] is a specialized tool for detecting and quantifying alternative splicing events from RNA-seq data. IncLevel value which quantifies the inclusion level of splicing events is defined as follows:

$$[eqn]$$

where IJC (Inclusion Junction Counts) is the number of RNA-seq reads mapped to splice junctions with the occurrence of target alternative splicing events, and SJC (Skipped Junction Counts) is the number of reads mapped to junctions with the absence of target events. It identifies five major types of AS events: SE, A5SS, A3SS, MXE, and RI. The workflow consists of the following steps: (i) indexing the rhesus macaque reference genome FASTA file using HISAT2-build (v2.2.1) [51]. (ii) aligning RNA-seq reads to the reference genome using HISAT2. (iii) converting SAM files to sorted BAM files using SAMtools (v1.21) [52]. (iv) applying rMATS to the resulting BAM files to identify AS events and DASEs between adjacent developmental stages. DASEs were defined as events with an |IncLevelDifference| > 0.05 and a false discovery rate (FDR) < 0.05.

4.5. GO Enrichment Analysis

GO functional enrichment analysis was performed using clusterProfiler (v4.14.6), with org.Mmu.eg.db (http://bioconductor.org/packages/org.Mmu.eg.db/, accessed on 29 November 2024) (v3.20.0) serving as the annotation database. GO terms with an adjusted p-value (p.adjust) < 0.05 were considered significantly enriched [53].

4.6. WGCNA

4.6.1. Construction of Co-Expression Modules

WGCNA [54] analysis was performed on the gene expression data of DEGs from the 8C-MOR period. First, missing value checks and outlier sample identification were conducted (Figure S4a). Then, standardized Z-scores were used for cluster analysis. Finally, WGCNA was used for module classification. The soft threshold power was determined to be 8 based on the scale-free topology fit index (R^2^) and mean connectivity (Figure S4b). The network type was set to “unsigned” (considering only the absolute value of correlation), and the correlation type was “pearson”. The minimum module size (number of DEGs) was set to 30, and the module merging height threshold was set to 0.25.

4.6.2. Identification of Hub Genes

Within each WGCNA module, intramodular connectivity was evaluated using the intramodular connectivity function. Two metrics were calculated: module eigengene-based connectivity (kME) and intramodular connectivity (kWithin). kME reflects the correlation between a gene and the corresponding module eigengene, whereas kWithin measures the extent of co-expression between a gene and other genes within the same module. Genes with |kME| ≥ 0.8 and kWithin ≥ 0.5 were defined as hub genes.

4.7. Correlation Analysis Between DESFs and DASEs

To explore the relationship between DESFs and DASEs from the 8C to MOR stage, a systematic correlation analysis was conducted. First, the TPM values of all DESFs and the IncLevels of all DASEs during this developmental period were extracted. Pearson correlation coefficients were then calculated between each DESF and all DASEs to assess the strength of linear associations. Based on the correlation coefficient matrix, a clustered heatmap was generated to visualize the global association patterns between DESFs and DASEs. To further investigate the potential regulatory roles of SFs in DASEs, a stringent threshold of absolute correlation coefficient (|r| > 0.8) was applied to identify strongly correlated DESF-DASE pairs.

4.8. Binding Motif Analysis for SFs

To further investigate the regulation of AS by SFs, we first extracted the sequence of the alternative exon and its flanking regions (500 bp upstream and downstream) for DASEs (SE type) that were highly correlated with SFs. Secondly, we retrieved the binding motifs for the relevant SFs from the CISBP-RNA database [55]. Thirdly, the extracted sequences and SF motifs were imported into the Find Individual Motif Occurrences (FIMO) tool (v 5.3.0) [56] within the MEME Suite (v5.3.0) to identify potential binding sites for SFs in the DASE sequences, using a p-value < 0.01 as the screening threshold. Finally, the Integrative Genomics Viewer (IGV) (v2.19.4) [57] was used to visualize the SF binding locations.

5. Conclusions

In summary, this work provides a comprehensive view of AS dynamics and regulation during early rhesus macaque embryogenesis, revealing AS as a critical and coordinated layer of gene regulation during ZGA. By elucidating the interplay between SFs, AS, and transcriptional activation, our findings offer important insights into the molecular logic of early primate development and provide a valuable framework for understanding AS regulation during early human embryogenesis.