From spring to autumn 2021, seeds of Youbi No.5 (R. communis line) were sown and planted under standard field conditions in Yangluo experimental field at OCRI, CAAS. The plant materials came from National Infrastructure for Crop Germplasm Resources of China. Seed samples were extracted from the fruits after removal of the pericarp for experimentation. Totals of 15 unmatured seed from 7 days after flowering (DAF) (S1), 14 DAF (S2), 21 DAF (S3), 28 DAF (S4) and 35 DAF (S5) (covering the whole development stages of seeds) were collected with three biological replicates, which were promptly frozen and stored at -80°C for transcriptome sequencing and metabolite profiling.

The protocols of metabolome profiling were performed according to Wang et al. (2019b) and Cao et al. (2019) and implemented by Wuhan MetWare Biotechnology Co., Ltd. (www.MetWare.cn). Briefly, three biological samples are freeze-dried and crushed with a zirconia bead for 1.5 min at 30 Hz. Then, 100mg of lyophilized powder was dissolved with 1.2 mL 70% methanol solution. Following centrifugation, the extracts were filtrated and then metabolome analysis was performed using UPLC/MS/MS system (SHIMADZU Nexera X2, Japan; Applied Biosystems 4500 Q-TRAP, USA) as described by Guo et al. (2019). The identification of metabolites was achieved by integrating multiple dimensions of analytical data, including the exact mass, MS/MS fragmentation patterns, isotopic distributions of MS/MS fragment ions, and retention time (RT). Specifically, a proprietary developed intelligent MS/MS spectrum-matching algorithm was employed to align the MS/MS spectra and RTs of metabolites detected in the project samples with those recorded in the in-house reference database (GB-PLANT). During the matching process, mass tolerance thresholds were set at 20 ppm for both MS and MS/MS data, and the RT tolerance was defined as 0.2 min. Multi reaction monitoring triple quadrupole mass spectrometry and peak area were used to quantify metabolites.

Principal component analysis (PCA) and correlation coefficients among the sample were determined using R software (v.4.5.2). Normalized metabolite data were used to compare all samples. Differentially accumulated metabolites (DAMs) between groups were determined by variable importance of projection (VIP) ≥ 1 and fold change > 1.5. VIP values were extracted from the OPLS-DA results generated using the R package (Metabo Analyst R) (v.1.0.1). The Metware database (MWDB) and other public databases were used to annotate metabolites. Metabolic pathways were subsequently annotated based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) database.

Following the standard protocol of Oxford Nanopore Technologies (ONT), total RNA was first extracted, assessed for purity and concentration, and only high-quality samples were used for transcriptome sequencing. Library construction was then performed, which involved primer annealing and reverse transcription to generate full-length cDNA with the addition of a switch oligo, followed by synthesis of the complementary strand. The resulting DNA underwent damage repair and end-repair was subsequently purified using magnetic beads. An ONT sequencing adapter was ligated using T4 DNA ligase (NEB). Finally, for sequencing, the prepared cDNA libraries were loaded onto FLO-MIN109 flow cells before running on a PromethION platform at Biomarker Technology Company (Beijing, China). The low-quality reads (length less than 500bp, quality score less than 6) and ribosomal RNA sequences were filtered out from raw sequences. Then, full-length non-chimeric (FLNC) transcripts based on the presence of primers at both ends of the sequence were obtained. FLNC transcripts were polished to obtain consensus isoforms before using minimap2 (Li, 2018) to remove redundancy based on the comparison results with the reference genome of R. commons (Chan et al., 2010). The final transcripts were ready for subsequent analysis.

The alternative splicing events in each sample were obtained through Astalavista software (Foissac and Sammeth, 2007), and the occurrence of the above five alternative splicing events in the transcripts was statistically analyzed according to Barbazuk et al. (2008) and Li et al. (2021a).

Based on Fold Change≥2 and FDR<0.01, DESeq2 (1.6.3) was performed for differential expression analysis. The p-value obtained from the null hypothesis test was adjusted using the Benjamini and Hochberg’s approach. Gene Ontology (GO, http://geneontology.org/), Non-redundant (Nr, http://ftp.ncbi.nih.gov/blast/db/FASTA/nr.gz), and Kyoto Encyclopedia of Genes and Genomes (KEGG, https://www.kegg.jp/) databases were used for gene annotation.

Twelve overlapping DEGs (28166.t000046, 29200.t000004, 29200.t000006, 29601.t000014, 29620.t000004, 29646.t000070, 29841.t000070, 29881.t000005, 29929.t000288, 30055.t000018, 30131.t000233, 30147.t000716, 30174.t000160, 30183.t000035, 60629.t00001, 60637.t00001), 6 DEGs belong to Cluster I related to castor oil accumulation (28035.t000007, 27985.t000040, 29673.t000033, 29726.t000092, 29682.t000014 and 27810.t000017) and 5 candidate DEGs from correlation analysis of transcriptome and metabolome (29794.t000071, 29917.t000061, 30147.t000162, 30008.t000036 and 28035.t000007) were verified by qRT-PCR analysis. With a reverse transcriptase (Vazyme, Nanjing, China), total RNA was reverse transcribed to cDNA. QRT-PCR reaction was performed in ABI PRISM7500 sequencing system (Applied Biosystems, USA) with ChamQ Universal SYBR qPCR Master Mix (Vazyme, Nanjing, China). QRT-PCR reaction procedure was performed as follows: 95 °C for 5 minutes, then 40 cycles at 95 °C for 15 seconds, 60 °C for 15 seconds, and finally 72 °C for 30 seconds. All reactions were executed using one biological sample with four technical replicates. Aata analysis adopted the comparative Ct method according to Schmittgen and Livak (2008). RcActin gene was used as the internal standard (Cao et al., 2019). The qRT-PCR primers were designed using Primer 6.0 software (Premier, Canada) (Supplementary Table S1).

Correlation analysis between the screened transcription factors (TFs) and lipid-related genes was conducted based on Pearson correlation coefficient using the built-in cor function in R software (v.4.5.2). The statistical significance of the correlation coefficients was determined by calculating the corresponding p-values with the corPvalueStudent function. Only the gene pairs meeting the stringent thresholds (coefficient of determination (R²) > 0.8 and p-value < 0.05) were retained. After removing duplicate gene pairs from the filtered results, the interaction network of the remaining gene pairs was visualized using Cytoscape software (v.3.10.4) (Shannon et al., 2003).

Correlation analysis between differential genes and differential metabolites were calculated based on Pearson’s correlation coefficients method. Before calculating the correlation, the z-value transformation method was used for data preprocessing, and then filtering was performed based on the correlation coefficient (r) and the p-value of the correlation. The filtering threshold was |r| > 0.80 and p-value < 0.05. The relationships between the metabolome and transcriptome data were illustrated using R software (v.4.5.2).

To gain deeper insights into the metabolic changes during seed development in R. communis, we conducted metabolite profiling across these five stages (S1-S5). Three biological replicates of immature seeds from each stage (S1-S5) were subjected to metabolomic analysis. Principal component analysis (PCA) revealed a clear separation among the five distinct seed developmental stages (S1-S5) (Supplementary Figure S1). The results demonstrated high correlation coefficients among the sample replicates, indicating the reliability of the metabolomic data (Supplementary Figure S2). Following quality filtering, totals of 790 known metabolites were identified in S1 and S5. Among these, lipids (17.1%), Phenolic acid (15.9%), organic acids (12.8%), flavonoids (10.0%), and amino acids and derivatives (9.7%) accounted for a substantial proportion (Supplementary Figure S3). Detailed information on all identified metabolites, including ionization models, molecular weights, KEGG pathways, and the names, categories and quantities of compounds across the five seed developmental stages, was provided in Supplementary Table S2.

Metabolomic profiling across five seed developmental stages (S1-S5) revealed dynamic changes in metabolite accumulation. A total of 10 pairwise comparisons were performed (Figures 1A-J; Supplementary Figure S4). The number of differentially accumulated metabolites (DAMs) varied among comparisons. Notably, stage S5 emerged as a critical phase for metabolic process, showing the highest number of common DAMs shared across comparisons involving this stage (Supplementary Figure S5). Three DAMs (isohyperoside, 1-O-(3,4-Dihydroxy-5-methoxy-benzoyl)-glucoside, and perillyl alcohol) were common to all ten comparison groups, suggesting their central roles throughout seed development.

All DAMs were categorized into 19 expression trend clusters via K-means clustering (Supplementary Figure S6A-S). KEGG enrichment analysis of DAMs indicated that metabolic shifts during seed development were prominently associated with amino acid biosynthesis, phenylpropanoid metabolism, and oxidative carboxylic acid metabolism across multiple comparisons (Supplementary Figure S7; Supplementary Tables S3-S6). Focusing on lipids relevant to castor oil composition, we identified several DAMs corresponding to major fatty acids, including hydroxy ricinoleic acid (Lmbn007891), ricinoleic acid methyl ester (Lmbn009747), ricinoleic acid (Zmyn004714), stearic acid (mws1489), palmitic acid (mws1488) and linoleic acid (mws1491). Notably, three ricinoleic acid-related metabolites (Lmbn007891, Lmbn009747, and Zmyn004714), based on MS/MS spectral matching, clustered together (subcluster 16) and exhibited a progressive increase in abundance from stage S1 to S5 (Figure 2; Supplementary Figure S6P), aligning with the active phase of oil accumulation during late seed maturation.

To identify candidate factors associated with oil accumulation in castor seeds, a total of 15 RNA-seq libraries were constructed from immature seeds collected at five developmental stages (S1-S5). High-quality sequencing data were obtained with mean read lengths spanning 845 to 1,555 bp across all samples. Clean reads showing high alignment rates (96.84-98.62%) were successfully aligned to the R. communis reference genome (Supplementary Table S7). A total of 32,194 genes were annotated for subsequent analysis. Differential expression analysis identified 9,940 genes that were differentially expressed in at least one pairwise comparison between stages (Supplementary Table S8). The greatest number of differentially expressed genes (DEGs) was observed in comparisons involving the earliest stage (S1), indicating a major shift in gene expression during initial development. In contrast, fewer DEGs were detected between the later stages (e.g., S3-S5), suggesting increased transcriptional stability as maturation progresses.

To characterize the expression dynamics of DEGs throughout seed development, K-means clustering was performed on the 9,940 DEGs based on log₂(FPKM + 1) values, using a similarity regulation model combined with hierarchical clustering. The DEGs were partitioned into five distinct clusters according to their expression trends. Cluster I (1,127 DEGs) exhibited a continuous upward trend from S1 to S5 (Supplementary Figure S8A). Cluster II (2,640 DEGs) displayed a bell-shaped profile, increasing from S1 to S3 and subsequently declining through S5 (Supplementary Figure S8B). Cluster III (1,519 DEGs) showed an initial decrease from S1 to S2, followed by stable expression in later stages (Supplementary Figure S8C). Cluster IV (2,865 DEGs) demonstrated a consistent down-ward trend across all stages (Supplementary Figure S8D). Finally, Cluster V (1,789 DEGs) exhibited a fluctuating pattern, characterized by a decline from S1 to S2, stability from S2 to S3, an increase from S3 to S4, and a final decrease from S4 to S5 (Supplementary Figure S8E).

Based on Oxford Nanopore Technologies (ONT) sequencing data, a total of 3,556 alternative splicing (AS) events were identified during seed development from S1 to S5. Intron retention (IR; 1,235 events) and alternative 3’splice sites (A3SS; 1,054 events) were the most common types, followed by exon skipping (ES; 625 events), alternative 5’splice sites (A5SS; 622 events), and mutually exclusive exons (MEE; 20 events). The number of AS events was highest at stage S1 and lowest at S5, indicating dynamic AS regulation during seed development The distribution of AS types varied across developmental stages (Figures 3A-E). IR, A3SS, and ES consistently represented most events in S1-S3. In S4, A3SS and A5SS proportions increased slightly, while IR and ES decreased. In S5, the profile shifted markedly: A3SS became the most abundant type (34.17%), followed by ES (29.17%), while IR declined to 23.33%. MEE also increased to 2.5% in S5 but remained below 1% in earlier stages. These shifts suggest that specific splicing patterns are associated with different phases of seed development, and phenotypic differences between S5 and the other four stages may be partially influenced by shifts in alternative splicing profiles.

We further investigated differentially expressed alternative splicing (DAS) events by overlapping DEGs and different AS events across stages S1-S5. Most DAS events occurred in comparisons involving S1 (Figures 3F-I), with fewer events detected between later stages (Figures 3J-O). This indicates that splicing-mediated regulatory changes are most active during the transition from S1 to S2, whereas later development (S3-S5) is characterized by relative splicing stability. Functional annotation of DAS genes highlighted several genes associated with lipid metabolism and hormone response. Specifically, SCAD (30032.t000002) was identified in the S1 vs S2 comparison; LACS1 (30076.t000180) in S1 vs S3; SBH2 (29813.t000027) and LIP1P (29804.t000062) in S1 vs S4; and ARF2 (29647.t000040) in S3 vs S5 (Supplementary Table S9). These findings suggest that alternative splicing plays an important role during seed development by generating distinct isoforms in a stage-specific manner.

To identify candidate genes potentially involved in castor oil accumulation, we intersected DEGs across multiple seed developmental stages to pinpoint common DEGs shared among ten comparative groups. This analysis revealed only 16 overlapping DEGs exhibiting distinct expression profiles (Table 1). Among these, the expression of five genes, 29200.t000006, 29620.t000004, 30174.t000160, 60637.t00001 and 60629.t00001, increased gradually during seed development (Figure 4A), suggesting that they might play a significant role in seed development. Functional annotation indicated that these genes encode proteins including glutelin type-A 3 precursor (LEGB, 29200.t000006), non-specific lipid-transfer protein A-like (nsLTP, 30174.t000160; 29620.t000004), ricin-like (Ricin, 60629.t00001) and agglutinin (RCA, 60637.t00001). GO enrichment analysis further highlighted their predominant involvement in metabolic processes (Figure 4B). The expression patterns of eight randomly selected common DEGs were validated by qRT-PCR, confirming consistency with the RNA-seq data (Supplementary Figure S9).

Key enzymes involved in castor oil biosynthesis, such as ACBP, ACCase, ACP, ACS, DGAT, EAR, FATA, FATB, FAH12, GPAT, KAS, KAR, LPAAT, LPCAT, PAP, PDCT, PDAT, PLA2, PLC2, PLD, SAD, as well as oil body-related proteins (OLE and PXG), were surveyed among the DEGs. A total of 30 genes potentially encoding enzymes associated with ricinoleic acid biosynthesis and six oil body-related genes were identified (Table 2). Based on their expression trends across five developmental stages, these genes were categorized into 5 clusters (Table 2; Supplementary Table S6). Thirteen genes, including OLEs, PDATs, PXGs, ACP1, ACCase, FATA1, FAH12 and DGAT2, were assigned to Cluster I, showing a continuously rising expression pattern, consistent with the accumulation pattern of three ricinoleic acid-related metabolites (hydroxy ricinoleic acid, ricinoleic acid methyl ester and ricinoleic acid). Ten genes, such as ACPs, FATB, PLA2, LPAAT2 and PLDs, fell into Cluster II, exhibiting a bell-shaped expression profile. The remaining genes were distributed among Cluster III (three genes: PDCT1, LPAAT4, and PLD), Cluster IV (six genes: PDAT1, ACP4, DGAT1, DGAT3, PLC2, and PLD), and Cluster V (four genes: SAD, PLD, PLC2, and FAFB), all displaying declining expression trends. Expression heatmap analysis of the 30 key enzyme and six oil body-related genes (Figure 5) revealed that one gene encoding ACP2 (29929.t000054) maintained relatively high expression throughout five seed development stages, implying a sustained role in oil synthesis. Additionally, seven genes including ACP1 (29709.t000042), DGAT2 (29682.t000014), ACP1 (29726.t000092), FATA1 (30217.t000013), LPAAT2 (27810.t000017), PDCT (29841.t000129) and SAD (29929.t000018), exhibited moderate expression across all stages, suggesting their involvement in maintaining oil biosynthesis. In contrast, eight genes, including PDAT1 (29706.t000035), LPAAT4 (30170.t000402), PLC2 (28833.t000008), FATB (30147.t000739), PLD (28725.t000003), ACCase (ONT.13196), PLD (29848.t000187), and ACP4 (230147.t000696), showed consistently low expression, indicating limited contribution to oil synthesis during the observed stages.

To verify transcriptome reliability, five randomly selected DEGs from Cluster I were subjected to qRT-PCR. Their expression trends aligned well with the RNA-seq data (Supplementary Figure S10), supporting the robustness of our transcriptomic findings.

Transcription factors (TFs), encoded by an extensive gene family in plants, serve as core regulators of gene expression. In this study, a total of 1,701 genes belonging to 69 TF families were identified through alignment of annotated transcripts against the AGRIS database (Supplementary Table S11). The top 20 TF families, each containing more than 30 genes, included C2H2, MYB-related, bHLH, AP2/ERF-ERF, NAC, MYB, WRKY, C3H, bZIP, GRAS, B3, GARP-G2-like, LOB, MADS-MIKC, Trihelix, TCP, MADS-M-type, HB-HD-ZIP, FAR1 and B3-ARF (Figure 6A). Enrichment analysis of DEG expression patterns demonstrated that a substantial number of lipid metabolism-related DEGs were predominantly clustered in clusters I and II (Supplementary Table S9). To further identify key TFs involved in lipid synthesis, we performed co-expression analyses between lipid-related genes and TFs within clusters I and II (Figures 6B, C). The results showed that 35 TFs, categorized into 17 distinct types, were identified in cluster I. Among these, the AP2/ERF family was the most abundant with 6 members, followed by the LOB family (5 members), MYB family (4 members), and bZIP family (3 members) (Supplementary Table S12). Notably, the AP2/ERF (e.g., WRI1, DREB) and LOB family members exhibited co-expression with multiple lipid-related genes in cluster I, such as FAH12, PDAT2, ACP1, FATA1, identifying them as candidate genes for functional validation in the oil accumulation process of castor seeds (Figure 6C; Supplementary Table S12). A total of 24 TFs, categorized into 20 distinct types, were identified in cluster II (Supplementary Table S12). Among these, only bHLH and Tify families contained three and two members, respectively (Supplementary Table S12). PLATZ, GARP-G2-like, MYB-related, GeBP and GRAS exhibited co-expression with multiple lipid-related genes in cluster II, such as FATB, LPAAT2, PLD (Figure 6C; Supplementary Table S12), indicating their potential important roles in lipid metabolism. Collectively, the identification of specific TF families, such as AP2/ERF and LOB in cluster I and PLATZ, GARP-G2-like, and others in cluster II, as co-expressed with lipid metabolism genes highlights them as strong candidate genes implicated in oil accumulation.

Pearson correlation analysis revealed 266 very strong (|r| > 0.9) and 338 strong (0.8 < |r| ≤ 0.9) correlation events between the top 100 DEGs and the top 50 DAMs (p-value < 0.05; Supplementary Figure S11). Co-enrichment analysis identified 9 KEGG pathways shared by transcripts and metabolites during seed development, including: biosynthesis of unsaturated fatty acids (ko01040), isoquinoline alkaloid biosynthesis (ko00950), tropane, piperidine and pyridine alkaloid biosynthesis (ko00960), phenylpropanoid biosynthesis (ko00940), ubiquinone and other terpenoid-quinone biosynthesis (ko00130), tyrosine metabolism (ko00350), phenylalanine, tyrosine and tryptophan biosynthesis (ko00400), phenylalanine metabolism (ko00360), and biosynthesis of amino acids (ko01230).

Furthermore, to examine the associations between 9,940 DEGs and 3 key DAMs, Pearson correlation analysis identified 1,982, 1,913, and 395 DEGs strongly correlated (0.8 < |r| ≤ 0.9, p-value < 0.05) with hydroxy ricinoleic acid (Lmbn007891), ricinoleic acid methyl ester (Lmbn009747) and ricinoleic acid (Zmyn004714), respectively (Supplementary Table S13). Notably, only six of these correlated DEGs were annotated as key enzymes involved in lipid biosynthesis or oil body-related proteins. Specifically, two DEGs (FAH12 and PLD) were key enzymes for ricinoleic acid biosynthesis and four DEGs (OLE1, OLE1, OLE and PXG) encoded as oil body-related proteins (Table 3). Among these, PLD (29841.t000111) was classified into Cluster IV, and its expression level shows a decreasing trend. Notably, five DEGs (OLE1: 29794.t000071, OLE: 29917.t000061, OLE1: 30147.t000162, PXG: 30008.t000036 and FAH12: 28035.t000007) were assigned to Cluster I, exhibiting an upward expression trend consistent with the accumulation patterns of the three key DAMs (Figure 5). The expression patterns of the five candidate genes were validated by qRT-PCR, showing strong agreement with the RNA-seq data (Supplementary Figure S10). FAH12 is a key gene influencing ricinoleic acid synthesis during seed development (Chen et al., 2007). Our findings provide further evidence confirming this earlier report. Oleosins (OLEs) and caleosins (PXGs) play a crucial role in encapsulating and stabilizing lipid droplets during triacylglycerol assembly (Badami and Kudari, 1970; Hanano et al., 2006). Collectively, these integrated transcriptomic and metabolomic findings highlight the essential role of further substantiate the pivotal function of FAH12 in ricinoleic acid synthesis and oil body protein-related genes (OLEs and PXGs) in oil stabilization and storage during castor oil accumulation process in Ricinus communis.

Metabolite profiling across different developmental stages identified several differentially accumulated metabolites (DAMs) associated with key fatty acids: including hydroxy ricinoleic acid, ricinoleic acid methyl ester, ricinoleic acid, stearic acid, palmitic acid and linoleic acid. Three ricinoleic acid-related DAMs (hydroxy ricinoleic acid, ricinoleic acid methyl ester and ricinoleic acid) exhibited a progressive increase from stages S1 to S5, indicating significant changes in ricinoleic acid content during seed development. Expression analysis of DAMs corresponding to stearic, palmitic, and linoleic acids revealed a declining trend for mws1489, whereas mws1488 and mws1491 showed upward trends (Supplementary Table S3). These metabolomic findings are consistent with previous reports. Zhang et al. (2016b) observed that oleic, linoleic, and ricinoleic acid levels increased steadily during castor seed development, aligning with the upward trends of mws1488 and mws1491 in our study. Furthermore, Chen et al. (2007) and Wu et al. (2021) reported that ricinoleic acid was undetectable in early developmental stages but accumulated substantially during mid-to-late stages, peaking at maturity. This pattern corroborates the expression profiles of the three ricinoleic acid-related DAMs identified in our metabolome analysis (Figure 2).

Transcriptomic analysis of developing castor seeds enabled the identification of key genes involved in castor oil accumulation. Wu et al. (2021) identified 10 key candidate genes (including OLEs, PXG, LACS7, PLA2, ACCase, DGAT, and PDAT) that exhibited increasing expression patterns during seed development. Similarly, Tian et al. (2019) reported 12 ricinoleic acid-related genes in Hiptage benghalensis, further supporting the functional importance of OLE, PXG, ACCase, DGAT, and PDAT in ricinoleic acid accumulation. In this study, transcriptomic analysis identified 30 genes encoding key enzymes in the ricinoleic acid biosynthesis pathway and six oil body-related proteins involved in castor oil accumulation and storage (Table 2; Supplementary Table S10). Among these, 13 genes, including OLEs, PDATs, PXGs, ACP1, ACCase, FATA1, FAH12, and DGAT2, were grouped into Cluster I, characterized by an upward expression trend during seed development (Table 2). Notably, approximately 85% of these genes, such as OLEs, PDATs, PXGs, FAH12, and DGAT2, function in the later stages of castor oil synthesis and accumulation. Their rising expression paralleled the increasing ricinoleic acid content. This co-occurrence suggests these genes are implicated in two key processes: ricinoleic acid biosynthesis (as enzymes) and castor oil accumulation/storage (as oil body-associated proteins). Additionally, the Cluster II genes with a “bell-shaped” expression trend also played a crucial role in the accumulation of ricinoleic acid in castor oil. Ten genes, including ACPs, FATB, PLA2, LPAAT2 and PLDs, were classified into this cluster (Supplementary Table S10).