To obtain hybrid seeds, we performed controlled pollinations between L. formosana and L. styraciflua trees. Developing fruits (multiples of capsules) were collected from the L. styraciflua trees. Whole fruits were surface-disinfested using the following sequence: rinse with water for 30 min, soaked in 75% ethanol for 1min and sodium hypochlorite for 4 min on the ultra-clean table, and cleaned with sterile distilled water for 5–6 times. The fruits were then dissected aseptically with a scalpel to remove the immature seeds. Each seed was nicked with a scalpel and cultured, ten seeds per plate, on induction medium: A modified Blaydes medium (Merkle et al., 1998) with 1.0 mg/L 2,4-dichlorophenoxyacetic acid (2,4-D), 0.5 mg/L 6-benzylaminopurine (6-BA), and 1.0 g/L casein hydrolysate (CH). The medium was Explants were transferred to fresh media of the same composition monthly until NEC (callus that were completely NEC or had both NEC and EC) started to form. Callus were transferred to fresh medium every 3 weeks. EC were transferred to embryo development medium –modified Blaydes basal medium without plant growth regulators (PGRs) – and incubated in the dark for 3–4 weeks, until embryos began to emerge from the PEMs. Individual somatic embryos were selected and cultured (16 per 90-mm plastic petri plate) on fresh medium of the same composition, also in the dark. When the embryos had enlarged to 2–3 mm in length, they were transferred to medium lacking PGRs and CH and incubated under cool white fluorescent light (140 μmol m^–2 s^–1) for germination. The cultures at different developmental stages, including NEC, friable EC, PEM after maturation cultured for 20 days (PEM1), PEM after maturation cultured for 45 days (PEM2), GE, HE, TE, CE, roots (R), stems (S), and leaves (L) were sampled, as shown schematically in Figures 1, 2. We collected the NEC, EC, PEM1, PEM2, GE, HE, TE, CE, and vegetative organs of plantlet samples (roots, stems, and leaves) from tissues in each of the three replicates and stored these at –80°C for RNA extraction.

A total of 1 μg of RNA per sample was used as the input material for the RNA sample preparations. Three experimental replicates per sample were used for RNA library construction. RNA-Seq libraries were generated using the NEBNext Ultra RNA Library Prep Kit for Illumina (New England Biolabs, United States) following the manufacturer’s recommendations, and index codes were added to attribute sequences to each sample. After sequencing, the raw data were filtered to remove adaptor contamination and low-quality reads. All clean reads were then mapped to the hybrid sweetgum genome using the HISAT2 program. The mapped reads of each gene were extracted using HISAT2 (Kim et al., 2015). We filtered the significantly differentially expressed genes (DEGs) using | log₂ (fold change)| ≥ 1.5 and P-value < 0.01 between seven pairwise comparisons: NEC_vs_EC, EC_vs_PEM1, PEM1_vs_PEM2, PEM2_vs_GE, GE_vs_HE, HE_vs_TE, TE_vs_CE, CE_vs_R, CE_vs_L, and CE_vs_S. Gene expression levels were normalized using the [fragments per kilobase of transcript per million fragments mapped (FPKM)] (Florea et al., 2013) method. All samples had more genes expressed at the level of 10 to 100 FPKM, and only 1196 to 1426 genes were identified at levels higher than 100 FPKM (Supplementary Table 1).

K-means clustering with Euclidean distance in MeV 4.9 software yielded 50 clusters based on inputs of log₂ relative FPKM values. In total, 13,574 genes residing in clusters with even more obvious tissue- or stage-specific expression trends were selected to make the final 35 K-means clusters. Only 18 of the 35 clusters are shown in Figure 5. TFs were extracted from the respective clusters. The log₂ relative FPKM value of each TF was used to produce heat maps in TBtools (Chen et al., 2018) (Figure 5B).

For co-expression network analysis, the WGCNA (Zhang and Horvath, 2005; Langfelder and Horvath, 2008) package was used. Based on the log₂ (1 + FPKM) values, a matrix of pairwise (spearman correlation coefficient) SCCs between all pairs of genes was generated and transformed into an adjacency matrix (a matrix of connection strengths) using the following formula: connection strength (adjacency value) = |(1 + correlation)/2|^β. Here, β represents the soft threshold for the correlation matrix, which gives greater weight to the strongest correlations while maintaining gene–gene connectivity. Based on the scale-free topology criterion, the β value was selected as 13 (Supplementary Figure 1A). The adjacent matrix was transformed into topological overlap (TO) matrix by TOM similarity algorithm, and the hierarchical clustering of genes was carried out based on TO similarity. The dynamic tree-cutting algorithm was used to cut the hierarchal clustering dendrogram, decomposed/merged branches and defined modules to achieve a stable number of clusters (Supplementary Figure 1B) (Langfelder and Horvath, 2008). For each module, a summary profile (module eigengene) was calculated via principal component analysis. Further, the modules with higher TO values (average TO for all genes in a given module) than those of modules comprised of randomly selected genes were retained. Finally, the gene interactions in key modules were visualized using Cytoscape software, and a gene co-expression network based on the hit number was constructed. Cytoscape software (Shannon et al., 2003) was used to visualize the gene interactions of key modules and constructed the gene co-expression network based on hit number. Considering the expression information, we marked these genes and confirmed the top 20 hub unigenes per group. The Cytoscape working data are available in Supplementary Table 2.

Gene Ontology (GO) enrichment analysis for DEG sets was performed using Cytoscape (BiNGO plug-in) (Maere et al., 2005). The P-value for enrichment was calculated for each represented GO term and corrected using the Benjamini–Hochberg error correction method. The GO terms exhibiting a corrected (after adjusting with the false discovery rate) P-value ≤ 0.05 were considered to be significantly enriched. Furthermore, pathway enrichment analysis of different sets of genes was performed using MapMan (v 3.6.0RC1) categories (significance value ≤ 0.05). MapMan bins of supercluster genes were assigned using the Mercator pipeline for automated sequence annotation (Thimm et al., 2004)¹.

Quantitative real-time polymerase chain reaction (qRT-PCR) was conducted on 33 RNA samples that were used in the preparation of sequencing libraries using a CFX Connect^® Real-Time PCR Detection System (CFX Connect, Bio-Rad, Munich, Germany). Each sample included three biological replicates. Nine genes that specifically highly expressed in non-embryogenic callus, embryogenic callus and somatic embryo development stage, respectively, were selected for qRT-PCR to verify the RNA-seq results, and then reverse transcribed of 1 μg total RNA into cDNA using a TransScript One-Step DNA Removal and cDNA Synthesis SuperMix Kit (Beijing Transgen Biotech Co. Ltd., Beijing, China). Each PCR reaction mixture contained 10 μl of Trans Start Top Green qPCR SuperMix (Beijing Transgen Biotech Co. Ltd., Beijing, China), 1 μl of template cDNA, 0.4 μl of forwarding primer, 0.4 μl of reverse primer (Supplementary Table 3), and 8.2 μl ddH₂O to a final volume of 20 μl with three technical replicates of each gene. The amplification condition for qRT-PCR was under the following program: denaturation at 94°C for 30 s, 40 cycles at 94°C for 5 s, 60°C for 15 s, 72°C for 15 s, with a melt cycle from 65 to 95°C. The apple EF1-α gene (Zheng et al., 2017) served as the internal standard. The relative gene expression in each sample was normalized to EF1 Ct value and calculated using the 2^{–ΔΔ Ct} method.

To study the SE of hybrid sweetgum, immature seeds of hybrid sweetgum (L. styraciflua × L. formosana) were tissue culture-induced to NEC. The NEC cells could be differentiated from the EC cells by visual discrimination. The NECs were dense and watery, white to brown in color, and had disorganized growth (Figure 1A). The ECs were friable, light yellow in color, and had organized globular structures that gave rise to SEs (Figure 1B). Once embryogenic cells have formed, they continue to proliferate and form PEMs. Browning takes place in PEMs when they are cultured on maturation medium for 20 days (Figure 1C), and the browning callus can regrow a new globular structure callus that continues to culture for 25 days (Figure 1D). Finally, globular, heart, torpedo, and cotyledon stage somatic embryos appear after 60 days (Figures 1E,F, 2).

To provide a clear understanding of hybrid sweetgum at the transcription level, we sequenced the 33 cDNA libraries constructed from 11 developmental stages: NEC, EC, PEM1, PEM2, GE, HE, TE, CE, roots (R), stems (S), and leaves (L) (Figure 1). In total, 998,501,629 clean reads (containing about 299.46 G of nucleotides) were obtained with quality checking. After alignment with the hybrid sweetgum reference genome, the clean data were mapped to the reference genome with a mapping ratio varying from 48,451,523 (76.36%) to 38,189,210 (81.49%). Among these, the unique mapping ratio varied from 47,112,326 (74.25%) to 37,112,232 (79.19%). In each sample, more than 91.08% had a base score of Q30 or higher. In our transcriptome, alternative splicing (AS) events were predicted by ASprofile. In total, 1,528,934 AS events, including alternative 5’ splicing, alternative 3’ splicing, exon skipping, intron retention, and the other 12 types of AS, were checked across the 33 samples. HE2 has the largest number of AS events (51,015), and the smallest was detected in L2 (37,151). A summary of the mapping statistics obtained for each sample is given in Supplementary Table 4.

The aim of this study was to reveal the potential key genetic factors related in the SE of hybrid sweetgum. Among these seven comparisons (Figure 3A), 3737, 5726, 7984, 6371, 136, 403, and 72 DEGs were identified, respectively. In comparison with NEC, EC had 1866 upregulated and 1871 downregulated genes. In comparison with EC, PEM1 had 3009 upregulated and 2717 downregulated genes. In comparison with PEM1, PEM2 had 3367 upregulated and 4307 downregulated genes. In comparison with PEM1, GE had 3433 upregulated and 2938 downregulated genes. In comparison with GE, HE had 49 upregulated and 87 downregulated genes. In comparison with HE, TE had 192 upregulated and 211 downregulated genes. In comparison with TE, CE had 36 upregulated and 36 downregulated genes. The Upset plot was also created to find the overlapped genes in the comparisons of eleven pairwise comparisons of hybrid sweetgum SE (Figure 3B). As shown in Figure 3B, 853 genes were exclusive to PEM1_vs_PEM2. 618 genes were exclusive to PEM2_vs_GE. Moreover, 598 genes were found in NEC_vs_EC. DEG analysis showed that the hybrid sweetgum transcriptome changed significant dynamically during SE, especially during the transition period from PEM1 to PEM2.

To evaluate the potential functions of the DEGs, we used GO term assignment to classify the functions of the DEGs in pairwise comparisons under three main GO categories: biological process, cellular component, and molecular function (Supplementary Figure 2). In all pairwise comparisons, the term with the largest proportion in biological process was metabolic process, followed by cellular process, single-organism process, biological regulation, respond to stimulus and localization, the term with the largest proportion in cellular component were cell and cell part, followed by membrane and organelle, the term with the largest proportion in molecular function was catalytic activity, followed by binding, transporter activity, structural molecular activity and nucleic acid binding transcription factor activity. In all pairwise comparisons among the functional groups, “metabolic process,” “cell,” “catalytic activity” terms were dominant. To investigate the biological pathways of the DEGs, we used the Kyoto Encyclopedia of Genes and Genomes (KEGG) database to classify the functions of the DEGs with emphasis on biological pathways (Supplementary Figure 3). According to the KEGG annotations, 701 DEGs (NEC_vs_EC) were assigned to 120 pathways, 1308 DEGs (EC_vs_PEM1) were assigned to 125 pathways, 1403 DEGs (PEM1_vs_PEM2) were assigned to 123 pathways, 1389 DEGs (PEM2_vs_GE) were annotated to 123 pathways, 18 DEGs (GE_vs_HE) were annotated to 21 pathways, 100 DEGs (HE_vs_TE) were assigned to 63 pathways, and 17 DEGs (TE_vs_CE) were annotated to 18 pathways. The pathway enrichment analysis showed that the annotated changes in NEC_vs_EC were mainly related to plant hormone signal transduction and fatty acid elongation. Among them, the genes annotated to fatty acid metabolism were up-regulated. The annotated changes in EC_vs_PEM1 were primarily involved in ribosomes and oxidative phosphorylation. The annotated changes in PEM1_vs_PEM2 were primarily involved in glutathione metabolism and phenylpropanoid biosynthesis. The annotated changes in PEM2_vs_GE were primarily involved in ribosomes and fatty acid elongation. Among them, the genes annotated to ribosome biogenesis in eukaryotes were up-regulated. The annotated changes in GE_vs_HE were mostly involved in brassinosteroid biosynthesis and starch and sucrose metabolism. Among them, the genes annotated to brassinosteroid biosynthesis pathway were significantly up-regulated, while the genes with alpha-linolenic acid metabolism were significantly down-regulated. The annotated changes in HE_vs_TE were primarily involved in phenylpropanoid biosynthesis and biosynthesis of amino acids. Among them, the genes annotated to phenylalanine metabolism were down-regulated. The annotated changes in TE_vs_CE were primarily involved in linoleic acid metabolism and alpha-linoleic acid metabolism. Among them, the genes annotated to alpha-linolenic acid metabolism and linolenic acid metabolism were down-regulated.

Based on the KEGG and GO annotations, plant hormone signal transduction, zeatin biosynthesis, and IAA biosynthesis were the representative pathways in our study. A large number of genes involved in auxin (95 DEGs) and cytokinin (40 DEGs) biosynthesis and signal transduction pathways were differentially expressed when comparing EC with NEC and different SE stages.

For example, the expression levels of PIN1-like, PIN2, AUX1, ARFs (ARF5, ARF17, ARF18, ARF19), GH3 (GH3.6, GH3.1, GH3.17), and five SAUR genes—all involved in auxin signal transduction—were significantly upregulated from NEC to EC. Nevertheless, AUX1, GNOM, GNOM-like, DRM1, PIN, PIN2, AUX4-like, TIR1, IAA9, ARFs (ARF4, ARF6, ARF19-like), GH3 (GH3.9, GH3.1, GH3, GH3.17-like), and SAUR (SAUR50, SAUR72) were mainly expressed in the NEC stage and downregulated in EC. From the EC to PEM and GE stages, GNOM, AFB2, GH3.6, IAA33, ARFs (ARF1, ARF1-like, ARF6, ARF18, ARF19, ARF19-like), and two SAUR32-like genes showed noteworthy upregulated expression (Figure 4A). In indole-3-acetic acid (IAA) biosynthesis, NIT and TSB were upregulated in EC. CYP714A1, CYP90A1, CYP1A1, CYP724B1, and three NIT genes showed an NEC-specific expression pattern. One YUCCA3, three CYP1A1, and one PAI gene were upregulated in PEM, and CYP724B1, CYP714A1, CYP90A1, ASA, and NIT remained upregulated during the development stage of SEs (Figure 4B).

As shown in Figure 4C, Three CKS, one IPT, and one Cis-ZOG family genes involved in cis-zeatin O-glycosylation were highly expressed in NEC and significantly downregulated from NEC to EC and stages of somatic embryo development. Two TPIT1 and one CKS gene were notably downregulated from NEC to EC and PEM and interestingly decreased during somatic embryogenesis, and the expression trends of these genes were consistent in studies on somatic embryogenesis of Longan (Chen et al., 2020). One CYP735A was upregulated in somatic embryo development. For the cytokinin signaling pathway, four A_ARR, two AHP, three B_ARR, and CRE1 genes were mainly expressed in NEC but downregulated in EC. Two A_ARR, one B_ARR, and three AHP genes were upregulated in EC. One A_ARR, four B_ARR, and three AHP genes were upregulated in PEM, and one CRE1, one B_ARR, and three A_ARR genes showed upregulated expression during the development stage of SEs (Figure 4D). Analysis of differential gene expression patterns indicated that genes related to cytokinin signal transduction and synthesis were differentially expressed during somatic embryogenesis, which may be involved in the regulation of somatic embryogenesis in hybrid sweetgum. Multiple sequences of the same hormone-related genes were compared, and the highest alignment rate was 97.21%, while the lowest was 17.58%. Most of the gene sequences were different, indicating that these genes were different isoforms (Supplementary Figure 4).

Thirty-five K-means clusters (containing 19,957 genes) with the most obvious tissue- and stage-specific expression trends were selected. Among these, 18 clusters representing 6815 genes are shown in Figure 5A and Supplementary Table 5. The 18 clusters can be divided into four large clusters, which are NEC, early SE (EC, PEM1, PEM2, GE), somatic embryo development (GE, HE, TE, CE) and vegetative organ (L, R, S). Early SE-specific gene clusters (C23, C20, C26, and C19) and embryo development-specific gene clusters (C17, C19, C34, and C35) were identified (Figure 5A). The NEC and vegetative organ categories had highly unique transcriptomes; almost 60% of the genes with tissue-specific expression were found in the NEC (C2 and C14) and vegetative organ (C1, C4, C5, C7, C8, C12, C13, C16, and C32) groups. However, a unique transcriptome was also found during the somatic embryo development stage; almost 10% of the genes with stage-specific expression were in embryo development.

Clusters with similar expression trends were combined into 12 superclusters (Figure 5A), and specific genes or pathway enriched in specific superclusters were identified by MapMan Bins (Figure 5C and Supplementary Table 5). The top abundant MapMan bin for the NEC-specific supercluster (supercluster 1) was polyamine metabolism, protein homeostasis, and secondary metabolism. As shown in Figure 5C, the biological functions associated with early somatic embryogenesis, nutrient uptake, enzyme classification, multi-process regulation, protein modification, and RNA processing were enriched in EC, PEM, and GE. Genes that participated in major vesicle trafficking, redox homeostasis, and regulation of proteins, including biosynthesis, homeostasis, and translocation, were particularly enriched in the somatic embryo development stage.

Among the 6815 genes that make up the aforementioned 12 superclusters (Figure 5A), 142 are TFs belonging to 60 families (Figure 5B and Supplementary Table 5). These TFs exhibit different stage- or tissue-specific expression patterns, and the dynamic and stage- or tissue-specific expression patterns of these TFs may reflect the pivotal roles that they play. For instance, the EC-specific (supercluster 2) cluster has four RLK/Pelle family kinases, one C2C2 gene, and one ethylene response factor gene with potential roles in embryogenic competence. Interestingly, the embryo development-specific (superclusters 5–8) cluster has one GRF gene, indicating that this gene plays a key role in the morphogenesis of SEs.

As an alternative analysis tool, weighted gene coexpression network analysis (WGCNA), was adopted. WGCNA is a systems biology approach that aims to understand entire gene networks rather than single genes. In this study, we constructed a co-expression network based on the pairable correlation of gene co-expression trends in all sample tissues. Modules were defined as clusters of genes that were highly interrelated; The correlation coefficient between genes in the same cluster was high. This analysis resulted in 22 different modules, color-coded in the dendrogram (Figure 6A and Supplementary Figure 1B), where each branch constitutes a module and each branch is a gene (Figure 6A and Supplementary Table 2). The first principal component of a given module was the module eigengene, which represented the gene expression profile of the module. Due to the tissue-specific expression profile of the eigengene, 22 characteristic genes of 22 different modules were associated with different tissue types. In comparison with other co-expression modules, for instance, the module with the highest correlation with NEC was the lightcyan module (p = 8 × 10^–22, r = 0.99), whereas genes in the green module had a higher correlation with EC (p = 3 × 10^–2, r = 0.42) than NEC, except for PEM2. Genes in the blue module were comparatively highly correlated with early SE (EC, PEM1, PEM2). Interestingly, three co-expression modules, the pink, salmon, and purple modules, showed relatively high correlation (r ≥ 0.15) with somatic embryo development stages (GE, HE, TE, CE) (Figure 6B). Many of the modules were correlated with more than one vegetative organ of the plantlet; however, a few were correlated with only a specific vegetative organ of the plantlet. For example, the yellow module was specifically correlated (p = 7 × 10^–14, r = 0.95) with leaf tissue, whereas the midnight-blue and red modules were specifically correlated with root (p = 2 × 10^–20, r = 0.98) and stem (p = 5 × 10^–13, r = 0.94) tissues, respectively (Figure 6B).

Figures 7A,C,E,G show the eigengene expressions for the lightcyan (NEC), thistle (EC) modules, blue (EC, PEM1, PEM2) and pink (GE, HE, TE, CE) modules. The heatmap shows the expression profiles of all co-expressed genes from the blue, pink, lightcyan, and thistle modules. WGCNA can also be used to drawn gene networks, where each node represents a gene and the edges are the connecting lines between nodes, represent co-expression-related genes. The node gene with the highest connectivity, as the hub gene, may play a key role in different modules. The lightcyan, thistle, blue and pink module networks are shown in Figures 7B,D,F,H. Some of the more interconnected node genes are those of TFs. Only the top 20 genes with high connectivity were selected, and 10 TFs representing different families were obtained, including the WRKY, RLK-Pelle, MYB, and STE_STE11 TFs. Another nine hub genes were also identified as key proteins or enzymes that play a key role in early SE (blue module) (Figure 7F and Supplementary Table 2). In the pink module network, 77 of the 965 genes encode TFs. Selecting only the top 23 genes showing high connectivity weight resulted in 11 TFs representing distinct families, including the AP2/ERF-ERF, C2C2, C3H, bHLH, GRF, and bZIP TFs. Another 12 hub genes, namely cyclin-dependent protein kinase inhibitor (SMR4), seed biotin-containing protein (SBP65), and late embryogenesis abundant protein 6 (LEA6), are involved in regulation of the cell cycle, seed germination capacity, and embryo development ending in seed dormancy (Figure 7H and Supplementary Table 2). In the lightcyan module network (Figure 7B), 37 of the 416 genes encode TFs. Selecting only the top 19 genes showing high connectivity weight indicated 11 TFs representing distinct families, including the AUX/IAA (EVM0009101; EVM0009101), C2C2-Dof, GNAT, bHLH, AP2/ERF-ERF, and NAC TFs (Supplementary Table 2). Another eight hub genes, namely protein NRT1/PTR FAMILY 4.3 (NPF4.3), cytokinin dehydrogenase 1 (CKX1), and brassinosteroid-related acyltransferase 1 (BAT1), are involved in regulation of transmembrane transporter activity, cytokinin metabolic processes, and response to brassinosteroids. Interestingly, after we performed the co-expression analysis of the data again, we observed that one module displayed opposite expression patterns between the two types of calluses (NEC and EC) having different embryogenic potentials. For example, the thistle module showed low transcript abundance in NEC but high abundance in EC. In this module, 38 of the 212 genes encode TFs. Selecting only the top 15 genes showing high connectivity weight revealed 10 TFs representing distinct families, including the HB-WOX, AP2/ERF-AP2, MYB, B3-ARF, AP2/ERF-ERF, FAR1, and C2H2 TFs. Another five hub genes, namely, peptidyl-prolyl cis-trans isomerase (FKBP43), CASP-like protein 1E1, beta-fructofuranosidase, insoluble isoenzyme 1 (INV1), GDSL esterase/lipase (At5g45910), and WAT1-related protein (At1g25270) were involved in regulation of peptidyl-prolyl cis-trans isomerase activity, protein histidine kinase binding, carbohydrate storage, lipid catabolic processes, and transmembrane transporter activity (Figure 7D). Interestingly, the FKBP43 gene is a fundamental protein that regulates cell division, adhesion and elongation throughout the plant development and embryogenesis and required for the spatial organization of apical meristems. These results suggest that the FKBP43 gene may play a key role in the proliferation and maintenance of embryogenic ability of embryogenic callus.

Gene Ontology (GO) enrichment analysis was performed for each module using Cytoscape (BiNGO plug-in) (Maere et al., 2005) highlighted important biological processes represented by a range of co-expressed genes. The blue module formed a cluster of 2632 DEGs enriched in functions related to phosphorus metabolic processes (phosphate metabolic process, phosphorylation, and protein amino acid phosphorylation) and amine metabolic processes (aminoglycan catabolic process, chitin metabolic process, and aminoglycan catabolic process) (Supplementary Figure 5C). The pink module formed a cluster of 965 DEGs enriched in functions related to RNA metabolic processes (RNA processing, rRNA metabolic process, rRNA processing, and ncRNA processing) and responses to stress (response to cold, response to temperature stimulus, and response to abiotic stimuli) (Supplementary Figure 5D). The lightcyan module formed a cluster of 416 DEGs enriched in functions related to the catabolic process, carbohydrate catabolic process, polysaccharide catabolic process, pectin catabolic process, and cotyledon vascular tissue pattern formation (Supplementary Figure 5A). The thistle module obtained from the other WGCNA analysis formed a cluster of 211 DEGs enriched in functions related to regulation of metabolic processes (transcription, macromolecule metabolic process, RNA metabolic process, and nucleic acid metabolic process) and developmental processes (reproductive developmental process, anatomical structure development, and multicellular organismal development) (Supplementary Figure 5B).

To verify the reliability of transcriptome data, twenty genes related to SE were selected.

The GLP1-13 (EVM0020840), ARF5 (EVM0026553), ABI3-1 (EVM0020968), LTP3 (EVM0021586), WOX9 (EVM0013624), LEA6 (EVM0000925), SMR4 (EVM0007351), ECP40 (EVM0003818), and MYB36 (EVM0015487) genes with high specific expression in various stages of hybrid sweetgum transcriptome were selected as marker genes to measure expression by qRT-PCR (Supplementary Table 3). The qRT-PCR analysis showed that expression levels during different SE phases changed between the genes (Figure 8). WOX9 and MYB36, as selected molecular marker genes, were almost not expressed in NEC, while they mainly expressed during early SE, with the highest expression at EC stage and then decreased in the mature somatic embryos of SE. This showed that those molecular markers played a key role in the maintenance and acquisition of embryogenic ability. Meanwhile, ECP40, ABI3-1, and LEA6 were highly expressed or up-regulated during somatic embryo development processes, and minimally or undetected in NEC. Which suggested that those marker genes may promote the hybrid sweetgum somatic embryo development. LTP3, SMR4, and ARF5 were highly expressed or up-regulated during the whole SE processes, and low expression in NEC and vegetative organ (R, L, S), indicated that those molecular markers played a key role in SE. In addition, the transcription level of GLP1-13 was highly and specific expressed in NEC, and this gene may inhibit the transition from non-embryonic to embryonic cells. qRT-PCR verification of SE-related genes also demonstrated a high correlation between RNA-seq and qRT-PCR data.

As we all known that auxin and cytokinin are important factors in plant cell division, differentiation, and SE induction. In this study, auxin and cytokinin signal transduction and biosynthesis related genes were significantly expressed during somatic embryogenesis, among which there were many auxin signal transduction genes, such as SUAR (small auxin up RNA) family genes had the largest number. However, the levels of endogenous auxin and cytokinin are influenced by their exogenous application (Vondráková et al., 2011; Xu et al., 2013). Auxin is considered to be the central regulator of SE, which may be due to the establishment of an auxin gradient during SE induction (Yang et al., 2012). To date, the effect of exogenous auxin during SE has been reported in several studies (Yang et al., 2012; Chu et al., 2017). In accordance with these results, we found a mass of auxin-responsive genes differentially expressed between tissues with contrasting embryogenic potential and at particular stages in hybrid sweetgum SE.

In samples of EC, contents of many auxin related genes, such as ARF17, ARF18 and AUX1, were significantly higher than those of NEC. The expression levels of ARF9 and ARF2 were significant during somatic embryo development, and increased gradually with somatic embryo development, and reached the highest in cotyledon embryo maturation. It is well-known that members of the ARF and AUX/IAA transcription regulator/signaling families act in concert to regulate the expression of responsive genes (Vanneste and Friml, 2009; Kieffer et al., 2010). The level of endogenous IAA has been correlated with PEM formation and high-frequency SE competency (Zhou et al., 2017). The balance of free IAA and IAA conjugates during SE induction of Coffea canephora was critical for embryogenic potential (Ayil-Gutierrez et al., 2013). The conjugation of auxin is synthesized by the GH3 family (Chen et al., 2010), and PIN1 plays a fundamental role in maintaining the embryonic auxin gradients (Feher, 2019). To summarize, we found that some of these family genes (PIN1-like, PIN2, AUX1, ARFs, GH3, and SAUR) had lower expression in NEC. Most of them were dramatically upregulated in EC and downregulated during SE, indicating that the auxin-responsive genes played an important role in hybrid sweetgum SE. Tryptophan (Trp)-dependent IAA biosynthesis is an important pathway in higher plants, and exogenous application of doses of Trp and IAA produced similar enhancements during rice SE (Siriwardana and Nabors, 1983). In our study, the expression levels of TSB genes, key genes in Trp synthesis, were drastically upregulated in EC and downregulated during SE. No Trp synthesis gene showed NEC specificity with high FPKM, suggesting that the level of Trp was higher during SE than during NEC.

Aside from auxin being a main inducer of SE, the use of exogenously supplied cytokinin to induce SE has been well established in many species (Sagare et al., 2000; Cabrera-Ponce et al., 2018). Large numbers of transcripts involved in zeatin biosynthesis and signal transduction are differentially expressed during cotton SE (Xu et al., 2013). Furthermore, endogenous cytokinin levels have been shown to be higher in SE than in NEC (Lai and Chen, 2002). In our study, from NEC to EC, a total of 31 DEGs were implicated in cytokinin signal transduction, including four CRE1 (two upregulated and one downregulated), seven B-ARR (two upregulated and five downregulated), and 11 A-ARR (two upregulated and four downregulated) genes. In the zeatin synthesis pathway, TRIT1 was upregulated from NEC to EC, and most of the CisZOG, IPT, and CKS were downregulated in EC and remained low during SE. However, IAA and zeatin biosynthesis- and signal transduction-related genes regulate SE in a complex and comprehensive way. Additional study of these genes in hybrid sweetgum SE is required.

Many genes and TFs have been reported as key factors in SE induction. Some TFs, including LEC1, LEC2, FUS3, ABI3, WOX9, WOX2, and BBM, were identified as molecular marker genes for early microspore embryogenesis of Brassica napus (Malik et al., 2007). In our study, 11 molecular markers were TFs (AUX22, ABI3, ARF3, ARF5, AIL1, AIL5, AGL15, WOX11, WOX9, IAA29, BBM1, and MYB36), and some of these gene functions on embryogenesis have been well characterized in various plants (Figure 9 and Supplementary Table 6). LEC2, ABI3, and FUS3 are B3 domain-containing transcription factors; ABI3 is highly expressed in embryogenic masses of Coffea arabica, and it has been speculated that its activity is related to embryogenic potential (Freitas et al., 2019). The AP2/ERF family has two very important members, namely BBM and PLT2, which have been reported to play important roles in microspore, somatic, and zygotic embryogenesis (Tsuwamoto et al., 2007). Overexpression of BBM can trigger spontaneous formation of SEs in Arabidopsis thaliana and B. napus. Therefore, BBM1 can serve as a marker gene for embryogenic callus cells in B. napus (Boutilier et al., 2002). In this assay, BBM and ABI3 were highly expressed during the SE processes. However, BBM1 was more highly expressed during the embryonic processes, and it can be used as a remarkable marker for maintenance of totipotent potential in hybrid sweetgum. To date, The MADS-box gene of AGL15 is expressed preferentially during embryonic development, and the protein accumulation reaches the highest level; it is expressed at lower levels at the later stage of germination (Fernandez et al., 2000; Yang and Zhang, 2010). Notably, AGL15 was shown to be beneficial to somatic embryo development and promoted the production of secondary embryonic callus from cultured zygotic embryos (Harding et al., 2003). The homologous domain transcription factor WUS, which can regulate the formation and maintenance of stem cells, was an early marker gene for SAM initiation in embryos (Weigel and Jürgens, 2002). STIMPY/WOX9 played an key role in increasing the proliferation rate and preventing premature differentiation in emerging seedlings (Wu et al., 2005). WOX2 and WOX9 were highly expressed at the early stage of SE in Picea abies, and may play a role together in conifer embryo model (Palovaara et al., 2010). WOX11 and WOX12 genes promoted the division of procambium cells and played a key role in callus and adventitious root formation (Liu et al., 2014). In this assay, WOX9 was specifically expressed in EC and downregulated during SE. WOX11 was highly expressed throughout SE and almost not expressed in NEC or vegetative organs. Hence, these genes might play an important role in hybrid sweetgum SE.

As an effective data mining method, WGCNA can screen genes related to specific traits and conduct modular classification on large samples to obtain co-expression modules with high biological significance (Zhao et al., 2010). Weighted gene co-expression network analysis has been used widely for plants in recent years (Zhao et al., 2020; Liu et al., 2021; Panahi and Hejazi, 2021). Therefore, we constructed a hybrid sweetgum gene co-expression network using a WGCNA approach and identified co-expression modules using transcriptome data. In total, according to the gene expression pattern, 22 modules were identified, several of which showed functional specificity at different stages (Figure 6). Then, to investigate the functional biological roles of these DEGs, BiNGO (Maere et al., 2005) software was used (Supplementary Table 2).

Our analysis identified many genes that were differentially expressed in the thistle module (EC). Of those with annotations, we found that the expression of many TFs, such as the HB-WOX, AP2/ERF-AP2, MYB, B3-ARF, AP2/ERF-ERF, FAR1, and C2H2 TFs, increased in EC. TFs played an important role in regulating plant differentiation and development. For instance, MYB TFs may regulate cell proliferation and elongation (Wen et al., 2020). Several R2R3-MYB TFs are barely expressed in NEC, significantly upregulated from NEC to EC, and maintained at high levels during early SE in longan (Chen et al., 2020). WUSCHEL-related homeobox (WOX) genes encoding plant-specific homeobox TFs play key roles in plant development and growth (Li et al., 2019), can regulate stem cell maintenance and formation, and are an early molecular marker for shoot apical meristem (SAM) initiation in embryos (Weigel and Jürgens, 2002). AP2/ERF TFs function in the response to phytohormones and/or biotic or abiotic stresses. BBM and ESR2, which are known to be associated with cellular pluripotency in Arabidopsis (Matsuo et al., 2011; Kulinska-Lukaszek et al., 2012), are also involved in pluripotency acquisition in hybrid sweetgum callus. Therefore, these genes may play a key role in the acquisition and maintenance of embryogenic ability during somatic embryogenesis.

In the pink module (somatic embryo development), 23 hub genes that showed high connectivity weight included 11 TFs representing distinct families, such as the AP2/ERF-ERF, C2C2, C3H, bHLH, GRF, and bZIP TFs (Supplementary Table 2). Another 12 hub genes included SMR4, SBP65, LEA6, LIP1, and others. Based on their BiNGO annotations, some genes associated with SEs development are involved in stress signaling pathways and RNA processes, revealing the relationship between SE development and stress signaling. The rRNA gene are an important part of protein synthesis and ribosome production, and the change of rRNA can cause the change of nucleolus volume (Koukalova et al., 2005). Gaudino and Pikaard (1997) reported that exogenous application of cytokinin stimulated polymerase I transcription of rRNA genes in Arabidopsis (Gaudino and Pikaard, 1997). LEA proteins were first identified in cotton seeds more than 30 years ago (Danyluk et al., 1998) and were later found to be abundantly expressed in late zygotic embryogenesis in many plant species, including cotton, barley, rice, oilseed rape, and wheat (Dure et al., 1989). LEA gene was highly expressed in carrot callus and somatic embryo stage (Kiyosue et al., 1992). We speculated these genes may regulate the synthesis and accumulation of cellular proteins and small RNA through some stress signals, thus affecting somatic embryo development and promoting somatic embryo maturation.