The mature seeds of fine single plants of O. henryi were collected in Mengguan town, Guiyang city in China (latitude: 26°14′23”N, longitude: 106°25′12″W) in Nov 2017. The authority responsible for the O. henryi resources is the Mengguan Forestry Station in Guizhou Province, China, who provided permission to collect the seeds of O. henryi. The formal identification of the plant material was undertaken by Prof. Mingtai An (Guizhou University). The matuer seeds were treated with concentrated H₂S0₄ for 30 min, washed with tap water for 30 min, disinfected in 75% alcohol for 1 min and then treated with 2% NaClO for 8 min, followed by 5 rinses in sterile distilled water. They were then soaked in sterile water for 24 h to make them swell and soften, and finally, mature embryos were obtained for somatic embryo induction. The somatic embryo induction methods referred to Wu et al., (2020). NEC ( Figure 1A ), EC ( Figure 1B ), GE ( Figure 1C ) and CE ( Figure 1D ) were collected. The culture medium on the surface of the calluses was washed with distilled water and placed into a 5 ml centrifuge tube, which was quickly frozen in liquid nitrogen and then transferred to a freezer at -80°C for storage. After the samples were collected, endogenous hormone determination and transcriptome sequencing were performed.

The test samples were sent to Metware (http://www.metware.cn/). The contents of endogenous hormones, including AUX (indole-3-acetic acid, IAA), CKs (N6-isopentenyladenine, IP), trans-zeatin (tZ), cis-zeatin (cZ), dihydrozeatin (DZ), GAs (gibberellin, GA1, GA3, GA4, GA7, GA9, GA15, GA19, GA20, GA24, GA53) and abscisic acid (ABA), were detected by liquid chromatography-tandem mass spectrometry (LC−MS/MS).

Total RNA was extracted from NEC, EC, GE and CE using TRIzol reagent (Metvivre, Wuhan, China). RNA purity was examined using a NanoPhotometer spectrophotometer (Implen, CA, USA), RNA concentration was measured with a high-accuracy Qubit^® 2.0 fluorometer (Life Technologies, CA, USA), and RNA integrity was assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA).

Poly A-tailed mRNA was enriched by Oligo(dT) magnetic beads, and then fragmentation buffer was added to break the RNA into short fragments. Using the short fragment RNA as a template, the first-strand cDNA was synthesized with six nucleobase random hexamers, buffer, dNTPs and DNA polymerase I to generate double-stranded cDNA, and then AMPure XP beads were used to purify the double-stranded cDNA. The purified double-stranded cDNA was then subjected to end repair, an A-tail was added, and a sequencing adapter was connected. Then, AMPure XP beads were used for fragment size selection, and finally, PCR enrichment was performed to obtain the final cDNA library.

After the library was constructed, Qubit 2.0 was used for preliminary quantification, and Agilent 2100 was used to detect the insert size of the library. The next experiment was performed only when the insert size met the expectation. Then, the effective concentration of the library was accurately quantified by the Q-PCR method (the effective concentration of the library was>2 nM) to complete the library inspection. Finally, the cDNA libraries were sequenced using the Novaseq 6000 platform (Illumina, San Diego, CA, USA).

After obtaining sufficiently high-quality clean reads, Trinity software was used to splice clean reads to obtain unigenes for subsequent analysis. BLAST software was used to align the unigene sequences with the KEGG, NR, Swiss-Prot, GO, COG/KOG, and Trembl databases (E ≤ 1e-5), and HMMER software was used to compare them with the Pfam database to obtain annotation information for unigenes after predicting the amino acid sequences of unigenes.

FPKM (Fragments Per Kilobase of transcript per Million mapped reads) was used as an indicator to measure transcript or gene expression levels, and the threshold standard of DEGs was |log2Fold Change|≥1 and FDR < 0.05.

Six genes related to the regulation of endogenous hormones were selected for qRT−PCR verification, and the total RNA of the test samples was extracted by reverse transcription using a Thermo kit to obtain cDNA. A QIAGEN kit was used for qPCR detection. The specific reaction procedure was as follows: preheating at 95°C for 2 min, followed by 40 cycles at 95°C for 5 s and 60°C for 30 s. The relative expression level was calculated using the 2^-ΔΔCt method. HLM-Actin was used as the internal reference gene, and the primers are shown in Table 1 . The relative expression levels were determined with 3 biological replicates, and each biological replicate was determined with 3 technical replicates.

All data were assessed for significant differences using Tukey’s test with SPSS 21.0 software. All data are presented as the mean ± standard deviation (SD) of three replicates. The graphs were created with Origin 9.0 and Photoshop 2019 software.

There were significant differences in endogenous hormone contents at different SE stages of O. henryi (P<0.05, Figure 2 ). The content of IAA in NEC was the lowest, while the contents of CKs, GAs and ABA in EC were the lowest. The content of IAA in EC was higher than that in NEC. The contents of GAs and ABA in NEC were 5.5 and 12 times those in EC, respectively. The results showed that relatively high IAA and low CKs, GAs and ABA contents promoted EC formation. With EC development, the contents of IAA, CKs, GAs and ABA gradually increased and were the highest in CE, which were 3.7, 10.4, 10.7 and 12.7 times those in EC, respectively. These results indicated that these four endogenous hormones contents were increased with somatic embryo development in O. henryi, and each tissue stage had its specific hormone content.

In this study, a total of 12 samples were sequenced, and 105.67 GB of clean data was obtained. The clean data, Q30, Q20, GC, clean reads and transcripts of each sample are shown in Table 2 . The sequence and assembly quality levels of the 12 samples were high. The overall expression trends of the qRT−PCR and RNA-Seq expression levels of the 6 genes were consistent ( Figure 3 ), indicating that the transcriptome data were reliable and could be further analysed.

The functions of unigenes were annotated by BLAST software and were predicted by comparative analysis with seven databases (KEGG, NR, Swiss-Prot, Trembl, KOG, GO and Pfam). Among the 235394 unigenes distributed to each of the databases, 47.93%, 59.47%, 39.58%, 59.05%, 35.75%, 49.27% and 40.44% of the genes were annotated ( Additional file 1 ). The highest species match was found in Lupinus angustifolius (18.74%), followed by Cajanus cajan (17.21%) and Glycine max (4.23%) ( Additional file 2 ).

A total of 115981 genes were annotated in the GO database, among which there were a large number of genes in cell (76165), cell component (76047), organelle (59197), cell membrane (36686) and organelle component (27642) for the cellular component; binding activity (72917) and catalytic activity (61800) for molecular function; and cellular process (73879) and metabolic process (64693) for biological process ( Additional file 3 ).

A total of 84165 genes were annotated in the KOG database, and a total of 25 classification functions were obtained. The top three KOG functions were “General function prediction only” (28881), “Posttranslational modification, protein turnover, chaperones” (7702) and “Signal transduction mechanisms” (7492), and the lowest was “Cell motility” (37) ( Additional file 4 ).

In total, 112816 unigenes were mapped into 141 KEGG database pathways. The top eleven KEGG pathways were “Metabolic pathways” (20565), “Biosynthesis of secondary metabolites” (9661), “Carbon metabolism” (2930), “Protein processing in endoplasmic reticulum” (2296), “RNA transport” (2283), “Plant−pathogen interaction” (2232), “Biosynthesis of amino acids” (2207), “Amino sugar and nucleotide sugar metabolism” (2169), “Endocytosis” (2151), “Plant hormone signal transduction” (2018) and “Starch and sucrose metabolism” (1961). In summary, the annotated gene information could provide data support for exploring its biological functions and revealing the molecular mechanism involved in SE at different stages for O. henryi ( Additional file 5 ).

The numbers of DEGs in NEC vs. EC, EC vs. GE, and GE vs. CE were 11589, 8999, and 27982, respectively, totalling 38100, and the number of common DEGs was 876 ( Figure 4A ). In NEC vs. EC, 7716 DEGs were downregulated and 3873 were upregulated, of which the number of upregulated DEGs was 2 times that of downregulated DEGs. In EC vs. GE, 4358 DEGs were downregulated and 4641 were upregulated. GE vs. CE had the largest number of DEGs, with 14060 upregulated and 13922 downregulated ( Figure 4B ).

The DEGs of NEC vs. EC, EC vs. GE and GE vs. CE were annotated into the KEGG database to obtain 137, 135 and 141 metabolic pathways, respectively ( Additional file 6 ). Among them, metabolic pathways such as “plant hormone signal transduction”, “metabolic pathways”, “biosynthesis of secondary metabolites”, “circadian rhythm - plant”, “flavonoid biosynthesis”, “isoflavonoid biosynthesis”, “isoquinoline alkaloid biosynthesis” and “phenylpropanoid biosynthesis” were significantly enriched ( Figure 5 ).

The above analysis showed that there were significant differences in the contents of endogenous hormones at different SE stages. To understand how the changes in endogenous hormone contents were regulated by the gene expression of the plant hormone biosynthesis and signal transduction pathways, we focused on five metabolic pathways involved in the biosynthesis and signal transduction of auxin, cytokinin, gibberellin and abscisic acid. The results showed that (1) tryptophan metabolism (ko00380) was annotated with 42 DEGs; (2) zeatin biosynthesis (KO00908) was annotated to 27 DEGs; (3) diterpenoid biosynthesis (KO00904) was annotated to 24 DEGs; (4) carotenoid biosynthesis (KO00906) was annotated to 48 DEGs; and (5) in plant hormone signal transduction (KO04075), there were 192, 90, 97 and 69 differentially expressed signal transduction genes annotated with auxin, cytokinin, GAs and abscisic acid, respectively ( Additional file 7 ). The results showed that the differences in endogenous hormone contents at different SE stages were regulated by related genes encoding plant hormones.

I-III represent KEGG enrichment of differential genes at NEC vs. EC, EC vs. GE and GE vs. CE, respectively.

To further explore the related genes that affect the changes in endogenous hormones during SE for O. henryi, DEGs related to plant hormone biosynthesis and signal transduction were screened from the transcriptome data, as shown in Figure 6 . In tryptophan metabolism (ko00380), there was a DEG, YUCCA, that was not expressed in NEC, was expressed at low levels in EC and GE, and was expressed at high levels in CE. In diterpene biosynthesis (ko00904), one GA3ox and two GA20ox genes were highly expressed in NEC and low in EC, and these expression levels gradually increased with EC development. In carotenoid biosynthesis (ko00906), the expression patterns of CYP97A3 (cytochrome P450 enzymes), ZEP (zeaxanthin epoxidase), ABA2 and AAO3 (Arabidopsis aldehyde oxidase 3) were the same as those of GA3ox and two GA20ox genes during SE for O. henryi. In plant hormone signal transduction (ko04075), 1 small auxin-up RNA (SAUR), 2 type-B ARA-bidopsis response regulators (B-ARR), 2 GA-insensitive dwarf 1 (GID1), 3 DELLA, 1 PYR1-like (PYL) and 1 ABRE-binding factor (ABF) were screened. Their expression patterns were consistent with the trend of endogenous hormone contents in SE at different stages for O. henryi. It was speculated that these genes were involved in the biosynthesis of corresponding endogenous hormones and affected the SE process of O. henryi.

A large number of differential TFs were found during SE of O. henryi, with 713 in NEC vs. EC, 750 in EC vs. GE, and 1959 in GE vs. CE ( Figure 7 ). Among them, there were 66, 62 and 188 differential TFs involved in plant hormone regulation, respectively, 55 AUX/IAAs and 58 ARFs in auxin; 8 CER1s, 30 B-ARRs and 5 A-ARRs in cytokinin; 35 DELLAs and 32 TFs in gibberellin; and 15 ABFs in abscisic acid, for a total number of 238 differential TFs during SE for O. henryi ( Additional file 8 ). It is worth noting that most AUX/IAA TFs were downregulated during EC formation and GE differentiation into CE and highly expressed in NEC, while other TFs were upregulated and downregulated at different SE stages. These results suggest that AUX/IAA played important roles in EC formation and CE differentiation of O. henryi ( Figure 8 ), and the regulation of its development by other TFs is a complex process.