Immature seeds of P. mongolica were obtained every 10 days after pollination between June and September 2019, from adult trees in the Baiyinaobao Nature Reserve of Inner Mongolia, China. Cones were periodically collected every week for all desired developmental stages of ZEs. The megagametophytes containing different developmental stages of ZEs were divided for RNA sequencing based on stereomicroscope observations (von Arnold et al., 2002), as follows: SDYS1 collected after 1 week of fertilization, which likely included early embryos in the proembryogenic stage, SDYS2 including the cleavage polyembryony in the early embryogenic stage, SDYS3 including the dominant embryos, SDYS4 including the pre-cotyledonary embryos, SDYS5 including the cotyledonary embryo, and SDYS6 including the mature embryo (Figure 1).

For transcriptome sequencing, total RNA was obtained using the Plant RNA kit (Sigma-Aldrich, St. Louis, MO, United States) following the manufacturer’s instructions. RNA degradation and contamination were detected in 1% agarose gels. RNA quantity and integrity were assessed using the NanoPhotometer spectrophotometer (NanoPhotometer purchased from Implen, Munich, Germany) and RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, United States). Total RNA (1 µg) was used to generate sequencing libraries using the NEBNext^® Ultra™ RNA Library Prep Kit for Illumina (NEB, Ipswich, MA, United States) sequencing.

Differential expression analysis was performed using the DESeq2 R package based on all obtained reads per kilo base per million mapped reads (RPKM) (1.16.1) (Love et al., 2014). Pairwise differential expression analysis was conducted among six group samples in sequence. Each set of samples contains three replicates. Transcripts that have log2FC ≥ 1 or log2FC ≤ 1 and p-value < 0.05 were used to determine DEGs. Further analysis will be based on the fold-change value of EEGs using Expander 7 software (Ulitsky et al., 2010) with the K-means algorithm (Shamir et al., 2005).

The agriGO (version 2.0) (http://systemsbiology.cau.edu.cn/agriGOv2/) was used in Gene Ontology (GO) enrichment analysis of obtained DEGs, which was corrected by a p value less than 0.05. The cluster Profiler R package was used to test the enrichment of DEGs in KEGG pathways. Venn diagrams were constructed using online software (http://bioinformatics.psb.ugent.be/webtools/Venn/). Hierarchical clustering heat maps were built, and cluster analyses were conducted using MultiExperiment Viewer (http://www2.heatmapper.ca/expression/).

The R package WGCNA V1.41-1 was used for WGCNA analysis (Langfelder and Horvath, 2008). After filtering, the automatic network construction function blockwiseModules was used to construct co-expression modules based on gene expression values. Cluster correlation was used to construct cluster families within the network. The network image was created using Cytoscape software.

Total RNA was obtained from each sample using an RNA extraction kit (Takara, Japan), and cDNA was obtained using the PrimeScript™ RT Reagent Kit (Takara, Japan). PCRs were performed using AceQ qPCR SYBR Green Master Mix with the following thermal cycling conditions: 95°C for 10  min, 40 cycles of 95°C for 10 s, and 60°C for 60  s (de Vega-Bartol et al., 2013b). Relative gene expression was quantified using the 2−ΔΔCt method (Cao et al., 2016).

To evaluate the molecular mechanisms governing the embryogenesis of P. mongolica, we isolated ZEs from June to September after pollination (Figure 1). Gene expression profiles were developed using RNA-seq as described in the methods section. After trimming adapter sequences and removing low-quality and multi-matched reads, more than 136.7 million unique clean reads and 20.5- to 30.4-G clean bases remained. A total of 76–80% clean reads were mapped to the reference genome of the Norway spruce. All samples were highly correlated and therefore clustered together (Figure 2A). Moreover, a principal component analysis (PCA) result demonstrated that the three replicates of each developmental stages can be clearly distinguished (Figure 2B).

We further analyzed DEGs through pairwise comparisons between consecutive stages of development. A two-fold difference (log2FC ≥ 1 or ≤ −1) and false discovery rate of 0.05 or less were used as thresholds. Each pairwise comparison identified a large proportion of DEGs, including 12,473 in proembryogenic versus cleavage polyembryony stages, 8,868 in cleavage polyembryony versus dominant embryo stages, 5,389 in dominant embryos versus pre-cotyledonary embryos, 9,077 in pre-cotyledonary embryos versus cotyledonary embryos, 15,878 in cotyledonary embryos versus mature embryos, and 16,987 in mature embryos versus proembryo embryos (Figure 2C). Next, we selected eight transcripts as candidate genes for further investigation of their expression patterns using RT-qPCR to evaluate their transcriptome data. The expression of these genes was consistent with the transcriptome data (Figure 3).

Based on the significantly expressed genes, a heat map was generated after background subtraction and normalization to represent the global dynamics of gene expression during zygotic embryogenesis in P. mongolica (Figure 4A). Furthermore, Venn diagram analysis showed that 5,039 transcripts were differentially expressed from stages 1 to 2, 1,606 from stages 2 to 3, 639 DEGs from stages 3 to 4, 1,646 from stages 4 to 5, and 4,425 from stages 5 to 6. Additionally, 1,279 genes appeared in five comparisons across all six stages of embryo development (Figure 4B). In all, the greatest difference in the number of DEGs was in the early embryogenic stage and later to mature embryonic stage. In contrary, little difference was observed between stages 3 and 4. Based on these results, embryonic developmental stages can be divided into the following: stage 1 is the proembryogenic phase; stage 2 is the early embryogenic phase; stages 3–4 are the mid-embryogenic phase; stage 5 is the latter embryogenic phase, and stage 6 is the mature embryo in P. mongolica.

Furthermore, we classified these DEGs based on KEGG pathways, which were mainly enriched in protein processing of the endoplasmic reticulum, fatty acid metabolism, starch and sucrose metabolism, phenylpropanoid biosynthesis, plant hormone signal transduction, and MAPK signaling pathway (Figure 4C). Subsequently, GO enrichment analysis was performed to further characterize their biological roles. Most biological process (BP) terms were enriched in cellulose metabolic and biosynthetic processes, polysaccharide metabolic processes, and glucan metabolic processes. In the molecular function (MF) category, DNA-binding transcription factor activity and cellulose synthase activity were significantly enriched. The significantly enriched cellular component (CC) terms were cell wall and apoplast (Figure 4D). In brief, many DEGs were involved in the regulation of embryonic development, and these pathways have great significance for further revealing the molecular events involved in zygotic embryogenesis in P. mongolica.

Among the DEGs identified, 1,279 genes appeared in five comparisons across all six stages of embryo development. Based on functional annotations and accumulated data on model plant systems, statistical results showed that the expressions of all six LEA genes, seven seed storage proteins, three seed maturation protein-encoding genes, and one seed dormancy gene were different in every developmental stage (Figure 5A). These associated genes were enriched in post-embryonic development, seed development, and fruit development (Supplementary Figure S1). In addition, many epigenetic regulation components were found during embryogenesis. Several putative histone subunit homologs and methyltransferases had different expression profiles across developmental stages. Additionally, genes associated with glucose transport and stress protein were also differentially expressed (Figures 5B–D). Five sugar efflux transporters and three sugar transporter genes involved in glucose metabolism and the transport pathway were expressed during ZE maturation. Moreover, stress response genes and HSP family genes were obviously differentially expressed in all six developmental stages during ZE maturation (Supplementary Figure S1). Importantly, several tyrosine phosphatases and PP2C were downregulated in later embryogenic stages, suggesting that the ABA signal transduction cascade plays an important role in later embryogenesis.

Because multiple phytohormone-related genes are involved in the acquisition of embryogenic competence, we further analyzed the transcriptome data to gain insights into the function of plant hormones in zygotic embryogenesis (Figure 6A). In the early embryogenic stage, auxin responsive genes, auxin transport genes AUX/IAA, and gibberellin response genes were upregulated in early embryogenesis. Conversely, ethylene-responsive protein kinase Le-CTR1 and ABA-related protein PP2C phosphatase were downregulated in this stage. In addition, several PP2C was still upregulated in early to middle embryogenesis. EIN3, which could initiate downstream transcriptional cascades for ethylene responses, was upregulated in the pre-cotyledon embryo stage. Auxin transporter including AUX/IAA and PIN genes and auxin responsive genes were upregulated in the cotyledon embryo stage. Subsequent seed maturation, PP2C phosphatase, and AUX/IAA family genes were significantly increased. Taken together, these results suggested that auxin, GA, and ethylene may participate in early embryonic initiation and cotyledon embryonic formation. Auxin and ABA synergistically regulate late embryonic development and seed maturation in P. mongolica.

Furthermore, transcripts related to epigenetic regulation were also found in different developmental phases (Figure 6B). SAM-dependent carboxyl methyltransferase was found during early embryogenesis, suggesting that DNA methylation has an important role in embryonic development of P. mongolica. In addition, chromatin modification-related genes, such as Jumonji C-domain histone demethylase-encoding gene, showed maximum expression during early embryogenesis. The expression of three transcripts related to histone deacetylase, which is usually associated with transposon silencing, were also increased in early embryogenesis. Therefore, histone methylation and acetylation modification play a relevant role in the early embryogenic phase of P. mongolica.

Most changes in transcript levels occurred in the early embryonic pattering determination and mature embryos formation. DEGs between stage 6 and stage 1 were identified to further investigate transcriptional regulation in ZE development. In the megagametophytes including proembryos and auxin polar transport (AUX/IAA and PIN3), response TFs (MYB, AP2, bZIP, and B3) were highly abundant. More upregulated transcripts were highly enriched in response to stimulus, such as salt, drought, cold, and wound stress. Moreover, histone H2A, lysine-specific demethylase JMJ706 and JMJ25, and histone RNA stem-loop-binding two genes were up-expressed in early proembryos. On the contrary, histone acetylation-related genes (HAC1 and histone deacetylase 6) and chromatin formation or remodeling genes (MBD9, RRC1, and SNF2-related transcripts) were highly expressed at the mature phases. ABA response, dehydration, photosynthesis initiation, and storage/energy-related genes were associated with mature stages of embryo development. In addition, many DEGs were annotated as containing “terminal inverse repeats (TIR),” “leucine rich repeats,” or “NB-ARC domains,” which are molecular switches implicated in immune signal transduction mechanisms (Figure 6C). These genes were highly abundant at stage 1 but declined successively during later developmental stages. These results are consistent with the previous report of transcriptome in maritime pine (Gonçalves et al., 2005; Rodrigues et al., 2018).

Because all DEGs represent different functions, a WGCNA assay was used to define clusters and obtain gene modules with specific patterns of expression during the six-embryo development stages (Figure 7). Based on expression pattern and tissue specificity, 20 gene modules were identified. Six stage-specific modules were significantly correlated with one specific stage of embryonic development (r > 0.8 and p value < 0.001). Next, we analyzed the functions of these genes in each model and found that the translation initiation factor, CCAAT-box binding factor, and mRNA splicing factor were significantly associated with early embryogenic stage, suggesting that this period may be the transcriptional activation stage. Heat shock protein was significantly correlated with dominant embryonic developmental stage. Auxin-responsive protein was mainly enriched in the pre-cotyledon stage. Several late-embryogenesis abundant-like proteins were found during late embryogenesis and embryonic maturation stages. Therefore, LEA are likely key genes for embryo maturation in P. mongolica. In addition, the top 50 genes associated with MEdarkred were of unknown function. This may be due to incomplete genome annotation in P. mongolica.

To further identify embryonic markers in P. mongolica, the MGFR tool (El Amrani et al., 2015) was used to analyze the DEGs of all samples. Fragments per kilobase of transcript per million mapped reads > 30 in at least one embryonic stage were used as the standard for DEG evaluation. We identified 184 transcripts in pro-embryos, 20 in early embryos, 7 in dominant embryos, 18 in pre-cotyledon embryo, 29 in cotyledon embryos, and 164 in mature seed embryos. These transcripts seem to be candidate genes for embryo-specific markers (Figure 8A). Considering that TFs are major determinants of embryogenesis, they are best for use as embryonic markers. We investigated transcripts likely TFs, whose levels were at least fourfold higher than other transcripts in a set of phase-enriched markers. There were MYB, GATA, TCP, and AP2 TFs associated with pro-embryo initiation. MYB was a key TF family in early embryo formation; HLH and AP2 TFs were highly expressed in the embryonic maturation stage (Figure 8B).

In addition, we also identified the known transcription factor marker gene WOX during embryonic morphogenesis (Figure 8C). Pm1910, Pm116716, and Pm10344604, likely putative orthologs of WOX8, were found to be downregulated in early embryogenesis. Pm10344604, a putative ortholog of WOX11, was upregulated in the dominant embryonic stage. Pm658686 was likely WOX12 and showed high expression in this stage. Pm81351 (a putative ortholog of WOX8) was strongly induced in the pre-cotyledonary stage. Pm74653 (WOX5) was upregulated in the cotyledonary stage.

TFs play a significant role in cell fate determination during embryogenesis. Hence, we further identified TFs involved in embryonic development in P. mongolica based on model plant systems (Figure 9). In all TFs, 144 were specifically upregulated and 90 were downregulated from pro-embryo to early embryo formation; 175 TFs were specifically upregulated and 71 were downregulated from early embryo to dominant embryo generation; 94 TFs (62 upregulated and 32 downregulated) in dominant embryo stages compared with pre-cotyledon embryonic stages; and 129 TFs were upregulated and 66 downregulated in the cotyledon embryonic period. One hundred fifty-four were upregulated and 101 downregulated from the cotyledon embryonic stage to the mature seeds (Figures 9A,C). In addition, we further analyzed the types of transcription factors in all DEGs (Figure 9B). They are mainly enriched in MYB, AP2, WRKY, HLH, bZIP, and NAC TF families. Furthermore, we analyzed the TFs in differentially expressed transcripts shared by the six groups (Figure 9D). We found that MYB TFs were upregulated in early embryogenesis. A putative WRKY transcript was strongly upregulated in middle embryogenesis. The AP2 and MYB transcripts were crucial in the transition from early embryogenesis to the pre-cotyledonary stage. Additionally, GATA, TCP, BES1/BZR1, and WRKY TFs were upregulated in mature embryos.

Auxin plays a considerable role in cellular patterning during embryogenesis. Its biosynthesis and polar transport-related localization are crucial in plant embryogenesis (Pérez-Pastrana et al., 2019). A previous study found that auxin-induced genes, such as ARF16, significantly increase in early cotyledonary embryonic stages in P. pinaster (de Vega-Bartol et al., 2013a). ARF16 represses WOX5 transcription, which is required to maintain pluripotent columella stem cells (Radoeva et al., 2019). In our study, we found that the auxin transport carrier genes AUX and PIN3 were upregulated in early embryonic stages. Moreover, its expression was also increased from the early cotyledonary to later embryonic stages. Thus, these results suggest that auxin biosynthesis, auxin transport, and auxin response are more prominent in early embryogenesis and affect cotyledonary embryonic pattern formation. This result is consistent with a previous report that auxin influx (AUX1, LAX1, and LAX2) is needed for cellular patterning in early embryogenesis in Arabidopsis (Robert et al., 2015). In addition, we also found that WOX8 and several auxin response TFs were specifically expressed in early embryogenesis. These findings provide a framework for processes that the transportation and accumulation of auxin likely activate the expression of these TFs during early embryogenesis in P. mongolica.

In addition, abscisic acid (ABA) plays a critical role in maintaining embryogenesis by regulating tissue differentiation and determining the physiological characteristics of late embryos (Pérez et al., 2015). Moreover, protein accumulation is generally under the control of ABA in conifer (Högberg et al., 2001). In this study, we statistically analyzed the expression of PP2C, an ABA signal transducer and key factor protein phosphatase in different ZE developmental stages. Six transcripts annotated to PP2C phosphatase were downregulated in the cotyledonary embryonic stage. Moreover, genes associated with ABA response and dehydration are more highly expressed at the embryonic mature phases. These results suggest that ABA play a critical role in the embryonic maturation processes of P. mongolica.

Epigenetic reprogramming, such as histone posttranslational modifications, DNA methylation, and chromatin-remodeling maintenance, has co-regulated functions during early embryogenesis in pine (de Vega-Bartol et al., 2013b; Trontin et al., 2016). HAD2, an acetylation regulation gene, is particularly increased in the early cotyledonary stage in P. pinaster (Rodrigues et al., 2018). We found that histone H2A histone subunits and lysine-specific demethylases show a higher abundance in early embryogenesis. Thus, histone posttranslational modifications might drive the expression of the transcriptome of early embryos in P. mongolica. Furthermore, histone acetylation-related genes and chromatin formation or remodeling genes were found in early embryogenesis. Similar results have been reported during maritime pine and Norway spruce embryogenesis (Uddenberg et al., 2011; Rodrigues et al., 2018). In addition, our results showed that several SAM-dependent carboxyl methyltransferase and other methyltransferases show different expression profiles. In general, DNA methylation is important for embryogenesis in P. mongolica.

A large number of TFs have been reported to regulate cell fate differentiation events, which determine the first zygotic division during embryogenesis. For example, the MYB TF family is expressed in early seed development (Avilés-Viñas et al., 2021). Moreover, MYB interacts with bHLH to regulate seed maturation. Network analysis showed that these TFs also associated with FUS3 and WOX12 in seed maturation (Tian et al., 2020). In here, we focused on TFs associated with plant development, phytohormones, and stress responses. The plant development-related MYB, AP2, and bZIP TF families were upregulated during early ZE development. The expressions of LEA, HSP, and HLH genes were increased during late ZE development. Additionally, homologues of the PLATZ and SRF family of plant-specific TFs, whose roles are poorly known, were also observed.

WOX genes (WOX2, WOX8, WOX1, and WOX3) play an important role during the ZE development. WOX2 is expressed in the apical daughter cell, and WOX8 is specifically expressed in the basal daughter cells of the zygote (Sakakibara et al., 2014; Armenta-Medina et al., 2021; Avilés-Viñas et al., 2021). The YDA-MAP signaling cascade phosphorylates and activates WRKY2 to promote WOX8 transcription during zygote development (Ueda et al., 2011; Horstman et al., 2017). We also found that WOX8 was specifically expressed in the early embryogenic stage, WOX11/12 were specifically expressed in the dominant embryonic stage, and WOX5 was upregulated in the cotyledonary stage. Additionally, a homolog of the WRKY family was also upregulated during this period. This may be consistent with regulation patterns in Arabidopsis that auxin response leads to activation of WRKY2 and expression WOX during embryogenesis.

Several genes are sufficient for embryogenesis when ectopically expressed. For example, when LEC1 and LEC2 are expressed in seedlings, SE is induced (Osorio-Montalvo et al., 2020). SERK1 encodes a leucine-rich repeat (LRR) transmembrane RLK, which induces SE when ectopically expressed. In addition, several TFs, such as AGAMOUS-Like15 (AGL15) and AGAMOUS-Like18 (AGL18), may also induce SE development (Zheng et al., 2016; Wójcik et al., 2020). The APETALA2 (AP2) domain TFs BABY BOOM/PLETHORA4 (BBM/PLT4), PLETHORA2 (PLT2), and WOUND INDUCED DEDIFFERENTIATION (WIND) promote SE in seedlings (DuBuc et al., 2020). Our transcriptome profiling data represent a valuable foundation for identifying genes, which are specifically expressed in particular stages of embryogenesis. We established a set of stringently defined temporal markers for embryogenesis. Transcription factors, including MYB, AP2, GATA, and TCP, could be markers during morphogenesis in early embryonic stages. AP2 and bHLH transcription factors are markers in the later embryonic stages to seed maturation. Heat shock proteins are specifically expressed in the domain embryo stage and serval auxin response genes in the pre-cotyledon embryo stage; late-embryogenesis abundant proteins are abundant during later embryogenesis. Therefore, these genes may be sufficient for embryogenesis and SE formation when ectopically expressed in P. mongolica.

In summary, our results establish a set of stringently defined temporal markers for identifying the key players involved in embryogenesis of P. mongolica. Our datasets and approach are expected to facilitate the discovery of molecular mechanisms underlying patterns of embryonic development. Moreover, functional understanding of these shared and distinct expression patterns of signaling, transcriptional, and epigenetic factors will help to address how embryonic development shapes the divergence of seed development in P. mongolica. Furthermore, our data lay a foundation for the protection and utilization of P. mongolica resources.