For this study, we utilized seeds from the rubber tree variety GT-1 (Gondang Tapen 1), which is known for its excellent yield traits. Seeds were obtained from the experimental farm of the Chinese Academy of Tropical Agricultural Sciences (Danzhou, Hainan, China) and rubber seeds soaked in water until seed coat rupture stage, and then rubber seed of radicle protrusion and seedling emergence stages planted in sandy soils within a controlled chamber at the Rubber Research Institute, Chinese Academy of Tropical Agricultural Sciences, Haikou, China. The sandy soil was kept consistently moist, and the ambient temperature was maintained at 20°C (Yan and Cao, 2009). We collected samples from five tissues—the endosperm, cotyledon, embryonic axis, radicle, and plumule—during four stages of seed germination: seed imbibition, seed coat rupture, radicle protrusion, and seedling emergence. Three biological replicates were performed for each tissue.

Total RNA was promptly extracted on the same day as the samples were collected and stored at –80°C for further analysis. The rubber seeds of radicle protrusion and seedling emergence stage planted in sandy soil, and the tissue samples such as endosperm, cotyledon, embryonic axis, radicle or plumule separated from the seeds that take out from sandy soil and be cleaned, and flash-frozen in liquid nitrogen immediately after collection. Subsequently, the samples were placed in a mortar within liquid nitrogen and then grinded to powder in liquid nitrogen, and total RNA was extracted using TRIzol reagent (TransGen Biotech, Beijing, China) following the manufacturer’s instructions (Tang et al., 2010).

RNA extracted from various tissue samples at different germination stages was subjected to transcriptome sequencing. For the total RNA, three biological replicates were pooled. Transcriptome sequencing was conducted by Majorbio Technology Corporation using the NovaSeq 6000 system following the Illumina Sequencing Protocol. The quality of the raw reads was checked using ‘FastP’ (Chen et al., 2018). High-quality reads from each sample were aligned to the rubber tree reference genome available at https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/001/654/055 (Kim et al., 2015) using ‘HISAT2’. The aligned reads were subsequently utilized to construct a reference annotation-based transcript assembly, enabling the identification of uniquely mapped reads for each sample. The results were subjected to ‘featureCounts v1.5.0’ (Liao et al., 2014) to obtain read counts for all samples. Differences in gene expression were assessed using DESeq2 1.2.5R (Love et al., 2014), with the p-value and false discovery rate (FDR) adjusted p-value (or q-value) set at 0.05. Gene Ontology (GO) enrichment analyses, including the categories of biological process, cellular component, and molecular function, as well as Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses, were performed on the differentially expressed genes (DEGs) across the samples using the clusterProfiler R package (Yu et al., 2012).

The expression levels of each transcript across the 44 samples were calculated and normalized to fragments per kilobase per million mapped reads (FPKM). DESeq2 was used to identify DEGs according to the following criteria: FDR < 0.05 and fold change > 2. Subsequent enrichment analyses included KEGG and GO enrichment analyses, as well as reference to the Plant Transcription Factor Database (PlantTFDB 4.0) http://planttfdb.cbi.pku.edu.cn/.

Fresh samples of endosperm, embryonic axis, radicle and plumule from seeds of seedling emergence stage and cotyledon from freshly harvested seeds were fixed in 80% (v/v) ethanol for 24 h and then sliced by a vibratome (Leica, Germany). The vibratome parameters were as follows: speed = 1.00mm/s, amplitude =1.0mm, and thickness = 100 μm. Sections were stained with Oil Red O (Sigma, USA) and mounted in 60% (v/v) glycerol (Lin et al., 2022). Pictures were taken with a light microscope SZX16 (Olympus, Japan).

Total RNA was reverse-transcribed into cDNA using the PrimeScript 1st Strand cDNA Synthesis Kit (Takara, Dalian, China). RT−qPCR amplification was carried out using SYBR^® Premix Ex Taq™ II (Perfect Real Time) (Takara, Dalian, China) and the CFX96 Touch™ Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA) following the manufacturer’s protocol. The mitosis protein YLS8 (Long et al., 2016) served as an endogenous control for the normalization of RT−qPCR Ct values, and relative gene expression levels were determined using the ΔΔCt method (Livak et al., 2001). Gene-specific primers were designed utilizing the NCBI Primer-BLAST tool (http://www.ncbi.nlm.nih.gov/tools/primer-blast/). The primer sequences for the three candidate genes are listed in Supplementary Table 1 .

The results are presented as the mean ± standard error (SE), and three technical replicates used in qPCR, and the statistical analyses were performed using SPSS v17.

A total of 44 RNA-Seq libraries were constructed for five tissues of the rubber tree ( Supplementary Figure 1A ), including the endosperm, cotyledon, embryonic axis, radicle, and plumule, at four different germination stages ( Supplementary Figure 1B ). We obtained 306.32 Gb of clean data from the 44 libraries representing five tissues: 12 endosperm samples, 11 cotyledon samples, 12 embryonic axis samples, 6 radicle samples, and 3 plumule samples at 4 stages. The clean data from each sample exceeded 5.95 Gb, with more than 92.53% of the reads being Q30 bases. More than 92.22% of the clean reads were successfully mapped to the rubber tree reference genome (Hevea brasiliensis).

Differentially expressed genes (DEGs) were identified in the endosperm, cotyledon, and embryonic axis samples at various stages of seed germination. The number of DEGs ranged from 255 (C1_vs_D1) to 16,673 (A3_vs_B3) ( Supplementary Table 2 ). Across all four stages, 41 genes were expressed in the cotyledonous tissue ( Supplementary Figure 2B ), 24 genes were expressed in the embryonic axis ( Supplementary Figure 2C ), and no genes were enriched in the endosperm ( Supplementary Figure 2A ). Similarly, compared with the other groups, the A1 vs. B1, A1 vs. C1, and A1 vs. D1 groups shared the largest number of DEGs (3451). This trend was also observed in the A2 vs. B2, A2 vs. C2, and A2 vs. D2 groups, which included 7875 DEGs, as well as in the A3 vs. B3, A3 vs. C3, and A3 vs. D3 groups, which included 8686 DEGs. These results indicated that most of the genes were expressed during the early stages of seed germination compared to the later three stages.

GO and KEGG pathway enrichment analyses of the DEGs were also conducted. A total of 34 highly significant terms were identified, including 11 for molecular function, 9 for cellular component, and 14 for biological process ( Figure 1A ).

In the biological process category, pathways related to auxin signaling, cellular response to auxin stimulus, and hormone-mediated signaling were significantly enriched. This finding indicates the critical involvement of hormones, particularly auxin, throughout all stages of seed germination. Furthermore, KEGG pathway analysis demonstrated enrichment in 22 pathways, including terpenoid backbone biosynthesis, ubiquinone and other terpenoid-quinone biosynthesis, peroxisome, and MAPK signaling ( Figure 1B ). These results indicate that active terpenoid metabolism occurs during seed germination and potentially serves as the foundation for rubber synthesis.

Phytohormones, including auxin, abscisic acid (ABA), gibberellin (GA), and cytokinin (CTK), are critical regulators of seed germination and crop development (Liu et al., 2007; Miransari and Smith, 2014; Oracz and Karpinski, 2016; Xue et al., 2021; Yang et al., 2022). They also have the potential to influence laticifer development. Therefore, we obtained phytohormone-related genes http://hevea.catas.cn/home/index, specifically auxin, ABA, GA, and CTK, from in different tissues and at four stages of seed germination to analyze the relationship between laticifer differentiation and plant hormones.

Natural rubber, an essential hydrocarbon polymer used in industry, primarily consists of cis-1,4-polyisoprene (Cherian et al., 2019; Yamashita and Takahashi, 2020; Calvo-Polanco et al., 2021). Two rubber biosynthesis pathways exist (Long et al., 2021). The 2-C-methyl-D-erythritol-4-phosphate (MEP) pathway supplies isopentenyl pyrophosphate (IPP) monomers for rubber synthesis, while the MVA pathway is the primary biosynthetic route for rubber and steroids in the cytosol (Chow et al., 2012). With respect to the transcriptome dataset, we identified genes associated with these two pathways in various tissues and at various seed developmental stages. A previous study indicated that latex vessels were already formed in the embryonic axis and cotyledon (Wu et al., 1990), and the expression of rubber synthesis-related genes in the cotyledon and embryonic axes was verified in this study ( Figure 3A ).

In our study, genes from these pathways exhibited differential expression, with high expression levels in the embryonic axis and cotyledon ( Figure 3 ). MVA, the primary rubber synthesis pathway, exhibited increased gene activation throughout seed development, especially in the late stages of rubber biosynthesis involving REF/SRPP genes, where numerous genes, such as SRPP and CPT, were highly expressed in the embryonic axis and cotyledon ( Figure 3A ). scaffold0992_72818 is an acetyl-CoA acetyltransferase (AACT) that catalyzes acetyl-CoA to produce acetoacetyl-CoA in cotyledons at later three stages, and another AACT scaffold0001_3383134, is expressed in only all tissues at the seed imbibition stage.

The number and expression level of genes related to latex synthesis were greater in cotyledons than in the embryonic axis. Interestingly, several CPTs (scaffold0385_30655 and scaffold0387_920391), CBPs (scaffold1160_94208 and scaffold1368_123062) and REF/SRPP (scaffold0824_400587, scaffold0916_24536 and scaffold2538_3915) were expressed in the radicle or plumule at later two stages.

In the MEP pathway, the rate-limiting enzyme 1-deoxy-D-xylulose-5-phosphate reductoisomerase (DXR) catalyzes IPP production for rubber synthesis, and DMAPP catalyzes steroid synthesis. Two DXR homologs exhibited high expression during the seed imbibition stage but relatively low expression in the other three developmental stages ( Figure 3B ). These results highlight the differential expression of genes involved in rubber biosynthesis during seed germination and indicate that laticifer species are likely to be present in cotyledon and embryonic axis tissues.

We found that laticifers may exist only in the cotyledon and cotyledon tissues of rubber seeds; therefore, to classify and analyze genes in cotyledon and embryonic axis tissues with similar expression trends, we utilized the short time-series expression miner (STEM) approach in this study. This analysis identified eight genes (profiles 0, 1, 2, 7, 12, 17, 18, and 19) among the significantly differentially expressed genes (DEGs) with a significance threshold of P < 0.05. Among these cotyledon profiles, Profiles 0, 2, 17, and 19 exhibited distinct and noticeable gene expression trends, warranting further investigation. In Profiles 0 and 2, gene expression showed a downward trend during seed germination, indicating downregulation during seed germination. Conversely, the expression of profiles 17 and 19 demonstrated an overall upward trend, indicating upregulation during seed germination. Profile 17 revealed the activation of pathways related to nutrient supply and recycling, including the TCA cycle, endocytosis, and nitrogen metabolism in the cotyledon. Moreover, Profile 19 indicated the activation of energy metabolism pathways, such as oxidative phosphorylation, photosynthesis, and carbon fixation pathways, in the cotyledons ( Figure 4 ). Similarly, on the embryonic axis, the expression of the genes in profiles 0 and 2 tended to decrease, similar to what was observed in the cotyledons. In contrast, profiles 17 and 19 in the embryonic axis were associated with the enrichment of pathways related to terpenoid biosynthesis, plant hormone signal transduction, phenylpropanoid biosynthesis, flavonoid biosynthesis, and the phagosome. These findings suggest the activation of secondary metabolism pathways, including the MVA and MEP pathways, in the embryonic axis during seed germination ( Figure 5 ).

Transcription factors (TFs) play significant roles in regulating tissue development, natural rubber biosynthesis, and laticifer differentiation (Iglesias-Fernandez et al., 2013; Lee et al., 2015; Wang et al., 2017; Yu et al., 2017; Kumpeangkeaw et al., 2019). To analyze the possible relationships between transcription factors and laticifer differentiation in cotyledons and cotyledons, according to the Plant Transcription Factor Database (PlantTFDB), a total of 2,276 TFs were predicted for the rubber tree. Among these genes, 713 were upregulated, and 932 were downregulated at all stages of cotyledon development ( Figure 6A ). The top 10 transcription factor families in the cotyledons included the MYB, AP2/ERF, bHLH, NAC, WRKY, C2C2, bZIP, B3 superfamily, HB, and GRAS ( Figure 6B ). In all stages of embryonic axis development, 589 TFs were upregulated, and 652 were downregulated ( Figure 6C ). The top 10 transcription factor families observed in the embryonic axis were similar to those in the cotyledon but with varying numbers ( Figure 6D ). Furthermore, pathways containing both up- and downregulated TFs exhibited enrichment in processes related to plant hormone signal transduction, MAPK signaling, plant−pathogen interactions, and circadian rhythm ( Figure 7 ). According to the laticifer-specific expression data (unpublished), scaffold0675_129662 (auxin response factor), scaffold0778_314155 (ABSCISIC ACID-INSENSITIVE), scaffold0941_312542 (DELLA protein, GA pathway), scaffold1421_86909 and scaffold0017_1455021 (AP2) and scaffold0332_497007, scaffold0148_1097530, and scaffold1297_134611 (MYB) were upregulated in the cotyledons. In addition, scaffold0675_129662, scaffold1230_90996 (auxin response factor), scaffold0778_314155 (ABSCISIC ACID-INSENSITIVE), scaffold0941_312542 (DELLA protein, GA pathway), scaffold1421_86909, scaffold0017_1455021, scaffold0037_175787, scaffold0013_532383 (AP2), scaffold0332_497007, scaffold0148_1097530, and scaffold0073_653695 (MYB) were upregulated in the cotyledons. These TFs, which are specifically expressed in laticifer, may regulate laticifer differentiation and rubber synthesis.

There are few reports about laticifer cells during seed germination. We identified laticifer cells of the endosperm, cotyledon and embryonic axis before germination and the endosperm, cotyledon, embryonic axis, radicle, and plumule after germination by staining with Oil Red O. The result showed that laticifer cells were observed in the cotyledon and embryonic axis but not in endosperm before seed germination ( Figures 8A–C ). Laticifer cells were observed in all tissues except endosperm after seed germination ( Figures 8D–H ). The results reveled that the original laticifer cells are likely to develop with seed formation.

We validated the expression patterns of six DEGs identified through RNA-Seq analysis in two tissues (endosperm and cotyledon) and four seed developmental stages (A, B, C, and D). These six genes, namely, HDR (scaffold0309_1225356), GGPPS (scaffold0434_538628), GPPS (scaffold0221_308815), REF/SRPP (scaffold1222_100110) are involved in the MEP pathway, and SS (scaffold1368_123062), HMGCR (scaffold0419_722144) are involved in the MVA pathway. The qPCR analysis results closely matched the transcriptome data ( Figure 9 ), confirming the reliability of the DEG identification results.

Study reported that hormone-related genes promoted the differentiation of secondary laticifer cells (Loh et al., 2019). IAA was shown to promote laticifer differentiation (Suri and Ramawat, 1995). In our study, we observed enrichment of phytohormone pathways during seed germination, with a particular emphasis on the auxin signaling pathway, which exhibited significant changes. Only seven genes, including six auxin-responsive genes and one auxin transporter gene (scaffold0520_237508), were consistently expressed across all tissues during the four seed germination stages. The auxin transporter gene is likely involved in the development of laticifers and other tissues. Previous research has shown that ARF5 and WOX13, which are associated with the auxin signaling pathway, are highly expressed in newly developed stem tissues of high latex–yield rubber tree clones, suggesting that increased expression of ARF5 and WOX13 may promote laticifer differentiation (Sae-Lim et al., 2019). Most of these genes were expressed in the embryonic axis, where laticifer cells exist. The expression of the scaffold0675_129662 and scaffold1230_90996 TFs (involved in the auxin pathway), which are specifically expressed in laticifer, increased in both the cotyledon and embryonic axes, suggesting that auxin likely plays a role in laticifer differentiation or rubber synthesis.

As a stress-related hormone, the ABA signaling pathway mediates laticifer differentiation in rubber trees (Tian et al., 2015). In the present study, ABA-related genes exhibited specific expression patterns in different tissues and at different stages. For example, scaffold0050_2638274 (the PYL6 ABA receptor) was exclusively expressed in the endosperm, while scaffold0047_2077855 (an ABA receptor), scaffold0245_1434801 (associated with abscisic acid biosynthesis), scaffold4755_4240 (ABSCISIC ACID-INSENSITIVE 5), and scaffold0919_21529 (the abscisic acid G protein-coupled receptor) were expressed during the seed coat rupture stage. G protein-coupled receptor (GPCR) is a crucial developmental regulator that acts as a regulator of G protein signaling and participates in various developmental processes, including seed germination (Kansup et al., 2013). We observed increased expression of ABA-related genes in the endosperm during the seedling emergence stage, suggesting their vital role in nutrient supply during late seed germination (the seedling emergence stage). Additionally, several ABA-related genes were expressed in the embryonic axis during the seed imbibition stage, and scaffold0778_314155, a TF involved in the ABA pathway that is specifically expressed in laticifer, was upregulated in the cotyledon and embryonic axes, suggesting that ABA likely plays a role in laticifer differentiation or rubber synthesis.

GA-related genes exhibited low expression during the seed imbibition stage, with two genes, scaffold0288_1337849 (gibberellin 2-beta-dioxygenase 1) and scaffold5194_1884 (gibberellin 20-oxidase), exclusively expressed in the endosperm. GA is known to promote starch catabolism and reduce the accumulation of stored proteins (Hu et al., 2021). These two genes likely promote gibberellin synthesis and accumulation. Scaffold0629_69288 (gibberellin 2-beta-dioxygenase) was expressed only in the plumule, suggesting its role in plumule development. GA-related genes were predominantly expressed during the seed coat rupture stage in three tissues: endosperm, cotyledon, and embryonic axis. scaffold0941_312542, a TF involved in the GA pathway that is specifically expressed in the laticifer, was upregulated in the cotyledon and embryonic axes, indicating that GA may regulate laticifer differentiation or rubber synthesis.

The differentiation of laticifers depends on cytokinin levels (Suri and Ramawat, 1995). However, the expression pattern of CTK differed from that of the aforementioned hormones. CTK-related genes were not expressed during the seed imbibition stage, indicating that CTK likely does not function in the early stages of rubber tree seed germination. While CTK is abundant in the endosperm of rice and maize during seed germination (Rijavec and Dermastia, 2010), CTK-related genes are scarcely activated in the endosperm of rubber plants at all seed germination stages. This finding suggested that CTK levels are likely low. CTK-related genes were highly expressed in the radicle and cotyledon ( Figure 5C ), suggesting that CTK may promote the development of the radicle and cotyledon during seed germination. Furthermore, in vitro laticifer differentiation in Calotropis procera is dependent on cytokinin (Suri and Ramawat, 1995), and high expression of CTK-related genes likely promotes laticifer differentiation in cotyledons.

Studies have shown that brassinosteroids (BRs) pathway is related to plant morphogenesis (Wang et al., 2012). In this study, we analyzed the gene expression changes of BRs pathway during seed germination, but only three genes were identified in DEGs. And just two genes, scaffold0123_1858764 and scaffold0441_941076, differentially expressed in cotyledon at early period of seed germination. It is suggested that BRs may not be the main hormone in the differentiation of early laticifer cells ( Supplementary Figure 4 ).

The response to different hormones during seed germination may link hormones to seed germination rate and laticifer development, and provide farmers with convenient breeding methods. Meanwhile, the analysis of hormone related genes expression may provide effective clues for us to promote the differentiation of laticifer cells in seeds by certain hormones in the future.

According to the structure, laticifers in rubber tree is divided into articulated and non-articulated laticifers (Sando et al., 2009). Non-articulated laticifers derive from single cell, and articulated laticifers formed by the fusion of multiple cells (Castelblanque et al., 2016).

Exogenous PSK promoted the differentiation of laticifers in the epicormic shoots of rubber tree(Chao et al., 2023). While there are limited reports on the relationship between seed germination and rubber synthesis, our study revealed differential expression of genes related to rubber biosynthesis in specific tissues and at specific seed germination stages. Previous studies suggested that laticifer initiation occurs in tissues such as cotyledons and on the embryonic axis (Bruni, A et al., 1978). Similarly, laticifer initiation in the early stages of laticifer cell differentiation was found in the cotyledon and embryonic axes of rubber plants (Castelblanque et al., 2016), and our laticifer cells staining confirms that laticifer cells were already present in the cotyledon and embryonic axes of seeds before germination ( Figures 8B, C ). We hypothesized that the initiation of differentiation of laticifer cells began with seed germination. Laticifer initiation involves invasive growth between adjacent cells and subsequent branching and elongation until mature laticifer formation occurs (Tan et al., 2019; Johnson et al., 2021). In our study, during seed germination, SRPP, CBP, CPT, and genes related to rubber synthesis were preferentially expressed in the cotyledon and embryonic axes and subsequently expressed in the radicle, while little was expressed in the radicle, indicating that laticifer differentiation occurred first in the cotyledon and embryonic axes, followed by the radicle and plumule. Moreover, more CPT and REF/SRPP genes were expressed in the cotyledons than in the other cotyledons at all stages. Moreover, REF/SRPP was primarily expressed in cotyledons during the first three stages and then significantly elevated in expression in the fourth stage, indicating that cotyledons may contain more laticifers than embryos at the early stage. However, CPTs (scaffold0385_30655 and scaffold0387_920391), CBPs (scaffold1160_94208 and scaffold1368_123062), and REF/SRPP (scaffold0824_400587, scaffold0916_24536, and scaffold2538_3915) were expressed in the radicle or plumule at the radicle protrusion or seedling emergence stages, indicating that these genes may play a role in the differentiation of different tissues in laticifers. Previous studies did not detect latex synthesis during seed germination, but in the present study, we showed that CPT and REF/SRPP, which are specifically expressed in laticifers, were expressed at the transcriptional level in the cotyledon and embryonic axes during seed germination. These results verified that laticifers are likely present in the seed. In the MEP pathway, HDR can catalyze the transformation of HMBPP to IPP, which is the substrate of rubber synthesis, but HDR has a very low expression level in the latter three stages, which is contrary to the expression levels of a large number of rubber synthesis-related genes in the MVA pathway. These findings indicate that the IPP produced by the MEP pathway is not amenable to rubber synthesis. Study reported that the laticifer density can be used as a screening marker for high-yield rubber tree breeding and may be shorten the rubber tree breeding period(Tan et al., 2019). Therefore, the number of original laticifer cells may determine its production after maturity. In this study, we found that IAA, GA and CTK related genes expressed highly in laticifer cells containing tissues. It is possible that these hormones are involved in the initial development of laticifer cells. Hormone-treated seeds may be a new breeding method to enhance the initial development of laticifer cells and thus enhance their latex-producing potential in mature rubber tree.

A variety of transcription factors have been identified in the cotyledon and embryonic axes. As important regulators, TFs play crucial roles in cell differentiation and embryogenesis (Long et al., 2006). In this study, TFs associated with the cotyledon and embryonic axes were enriched, with the three genes being MYB, AP2/ERF, and bHLH ( Figure 6 ). Previous studies have shown that the AP2/ERF family is involved in cell differentiation and dedifferentiation (Iwase et al., 2011; Marsch-Martinez et al., 2006), while MYB and bHLH are involved in cell differentiation (Gómez et al., 2009; Katayama et al., 2015; Tominaga et al., 2014). Screening for the specific expression of these transcription factors in laticifer cells may be related to laticifer differentiation or rubber synthesis. We screened several laticifer-specifically expressed TFs, including those upregulated both in cotyledons and in the embryonic axis. As MYB-type transcription factors, scaffold0332_497007 and scaffold0148_1097530 are both upregulated in the laticifer cotyledon and embryonic axes. Similarly, scaffold0675_129662 (auin response factor), scaffold0778_314155 (ABSCISIC ACID-INSENSITIVE) and scaffold0941_312542 (GA pathway) were upregulated in the laticifer cotyledon and embryonic axes. The above genes, which are specifically expressed in laticifer, are likely involved in laticifer differentiation or rubber biosynthesis.