Adventitious root formation
During in vitro rooting, the absence of WPM or the use of low-strength WPM significantly promoted adventitious rooting of M. glyptostroboides shoots, when compared with the control group (full-strength WPM) (Fig. 1). The rooting percentage (Fig. 1b), average root number (Fig. 1c) and root length (Fig. 1d) increased significantly in low-strength WPM, and 2/5-strength WPM most effectively induced ARs, resulting in the highest rooting percentage (84.17%), root number (5.67), and root length (3.30 cm) after 30 d of treatment. Histological analysis (Fig. 2) showed that AR initiation was observed after 8 d in 2/5-strength WPM, evidenced by clusters of meristematic cells in the outer phloem. AR primordia were evident at 10 d, and ARs emerged from the epidermis at 12 d and elongated after 16 d. In the control, no AR primordia were observed between 8 and 10 d, only after 12 d. AR primordia emerged from the epidermis at 16 d and elongated at 20 d. Thus, the base of shoots (0.5-1.0 cm) collected at 8, 10, 12, 16, and 20 d from 2/5-strength WPM (marked as MG) and control (CK) were employed as the samples for endogenous hormone content and RNA-seq analysis.
Analysis of endogenous hormone content
Nine endogenous hormones, IAA, isopentenyladenine (IP), isopentenyladenosine (IPA), trans-zeatin riboside (TZR), trans-zeatin (TZT), GA1, GA3, GA4 and GA7, at 8, 10, 12, 16, and 20 d in MG and CK were detected in the stem bases of M. glyptostroboides (Fig. 3). Endogenous IAA mainly accumulated at 12 d, and levels at 8, 10, 12, 16, and 20 d in the control were significantly higher than in all MG treatments. There was no significant difference in IP content between MG and CK, except at 10 d. IPA content in MG was significantly lower than in CK, except at 16 d. TZR in stems accumulated from 8 to 20 d in MG and CK, and in MG was significantly higher than in CK at 20 d. TZT, GA1 and GA3 accumulated rapidly at 16 and 20 d in MG and CK (GA1 was not detected at 8–12 d). TZT content in MG was significantly higher than in CK at 20 d, whereas GA1 and GA3 contents were significantly lower in MG at 16 and 20 d. GA4 content increased sharply at 12 d, corresponding to the timing of ARs breaking through the epidermis in MG, then decreased rapidly at 16 and 20 d, corresponding to the timing of AR elongation. Unlike MG, CK displayed a significantly higher level of GA4 at 16 and 20 d, corresponding to the timing of AR primordia formation and breaking through the epidermis. Compared with CK, GA7 content in MG showed no significant difference at 8 and 10 d, but was significantly lower at 12, 16 and 20 d, corresponding to the stage of ARs breaking through the epidermis and extension.
In brief, 2/5-strength WPM induced AR formation from M. glyptostroboides stems, displaying a lower endogenous content of IAA, TZT, GA1, GA3 and GA7 at different stages, except TZT at the late stage of AR extension (20 d). At the stage of AR primordia formation (10 d) and breaking through the epidermis (12 d), stems showed lower contents of IP, IPA and TZR, and higher GA4 content at 12 d. At the early stage of AR elongation (16 d), stems showed a lower content of GA4 and higher contents of IPA and TZR. At the late stage of AR extension (20 d), a higher content of TZT was observed in stems.
Analysis of RNA sequences and identification of differentially expressed genes (DEGs)
To identify genes involved in AR formation from M. glyptostroboides stem bases, 30 cDNA libraries were prepared from three repeat mRNA samples collected from the 2/5-strength WPM treatment (MG) and the control (CK) at 8, 10, 12, 16 and 12 d, marked as MG8, MG10, MG12, MG16, MG20 and CK8, CK10, CK12, CK16, CK20, respectively. The sequenced raw data was submitted to the SRA at the NCBI database with the following accession number: PRJNA999065. The total number of raw reads obtained for each library ranged from 49,120,012 to 70,131,470 with Q20 > 97.32% and Q30 > 93.17% (Table 1). After filtering, the number of clean reads per library ranged from 48,871,938 to 69,828,774 (Table 1).
When the stems of M. glyptostroboides placed in 2/5-strength WPM and control medium were compared, 253 (MG8 vs CK8), 239 (MG10 vs CK10), 417 (MG12 vs CK12), 540 (MG16 vs CK16), and 364 (MG20 vs CK20) DEGs were identified on days 8, 10, 12, 16 and 20, respectively (Fig. 4a). The highest number of up-regulated genes was observed in the MG16 vs CK16 library (330) and fewest in the MG10 vs CK10 library (99) (Fig. 4a). The Venn diagram indicates that three genes maintained significant differential expression in the 2/5-strength WPM treatment from between 8 and 20 d (Fig. 4b).
GO and KEGG enrichment analysis of DEGs
The significant enrichment GO terms of DEGs showed a few differences in the five compared libraries (Fig. S1). In MG8 vs CK8, the significantly enriched terms were “glutamate catabolic process”, “dicarboxylic acid catabolic process”, “xyloglucan: xyloglucosyl transferase activity”, “xyloglucan metabolic process”, etc. In MG10 vs CK10, the most significantly enriched terms were “sexual reproduction”, “reproduction”, “mannose binding”, etc. In MG12 vs CK12, “dipeptide transmembrane transporter activity”, “tripeptide transmembrane transporter activity” and “cell wall”, etc. were indicated as being the most significantly enriched terms. In MG16 vs CK16, “lactoperoxidase activity”, “hydrogen peroxide catabolic process”, “hydrogen peroxide metabolic process”, etc. were the main significantly enriched terms. In MG20 vs CK20, the most significantly enriched terms were “sexual reproduction”, “reproduction”, “xyloglucan metabolic process”, etc. The KEGG pathways significantly enriched in each stage displayed a significant difference between each treatment (Fig. S2). In MG8 vs CK8, the significantly enriched DEG pathways were “taurine and hypotaurine metabolism”, “butanoate metabolism”, “alanine, aspartate and glutamate metabolism”, “plant hormone signal transduction”, etc. In MG10 vs CK10, DEGs in “arginine and proline metabolism” and “vitamin B6 metabolism” pathways were significantly enriched. In MG12 vs CK12 and MG16 vs CK16, “phenylpropanoid biosynthesis” was the most significantly enriched pathway. In MG20 vs CK20, “plant hormone signal transduction” and “pentose and glucuronate interconversions” pathways were the most significantly enriched.
DEGs enriched in the plant hormone signal transduction pathway
KEGG pathway enrichment analysis showed that 34 DEGs were enriched in the “plant hormone signal transduction” pathway, including genes encoding IAA-amido synthetase Gretchen Hagen3 (GH3), auxin/IAA (AUX/IAA), Protein TIFY 10b (TI10B), pathogenesis-related protein (PR), auxin-responsive proteins (SMALL AUXIN UP RNA, SAUR), and xyloglucan endotransglucosylase protein (XTH) (Fig. 5). Three GH3 genes were identified, one of which (Mgl0097120) was exclusively up-regulated in MG16 vs CK16 and two (Mgl0098450 and Mgl0097090) that were differentially expressed in MG20 vs CK20. Only one IAA27 gene (Mgl0010460) was identified; it was significantly up-regulated in MG12 vs CK12 and MG16 vs CK16. Four DEGs encoding jasmonate ZIM domain-containing protein (JAZ) were significantly up-regulated in MG8 vs CK8 and MG16 vs CK16, respectively, including three TI10B DEGs (Mgl0152960, Mgl0211930, Mgl0211910) and one hypothetical protein (Mgl0211950). Four SAUR and three PR1 DEGs showed differential expression in MG8 vs CK8 (Mgl0056630, SAUR12), MG12 vs CK12 (Mgl0056630, SAUR12; Mgl0266300, Mgl0240770, PR1), MG16 vs CK16 (Mgl0056790, SAUR; Mgl0222810, SAUR71; Mgl0266300; Mgl0240770; Mgl0240810, PR1) and MG20 vs CK20 (Mgl0038910, SAUR32; Mgl0240770, PR1). Nineteen XTH DEGs were identified, but only four XTH1 (Mgl0287740, Mgl0287830, Mgl0287700, Mgl0287730) and one XTH34 (Mgl0097340) were significantly up-regulated in MG8 vs CK8, two XTH1 (Mgl0287710 and Mgl0287780) in MG20 vs CK20, one XTH1 (Mgl0282260) and one XTH5 (Mgl0289590) in MG16 vs CK16, nine XTH1 (Mgl0287680, Mgl0287690, Mgl0287720, Mgl0287750, Mgl0287760, Mgl0287870, Mgl0287860, Mgl0287850, and Mgl0287880) in MG8 vs CK8 and MG20 vs CK20, and one XTH1 (Mgl0287840) in MG8 vs CK8, MG16 vs CK16 and MG20 vs CK20 (Table S1).
DEGs enriched in the phenylpropanoid biosynthesis pathway
A total of 52 peroxidase (PER) DEGs were enriched in the “phenylpropanoid biosynthesis” pathway, most of which were differentially expressed in MG12 vs CK12 and MG16 vs CK16, corresponding to the stage of AR emergence from the epidermis and elongation (Fig. 6, Table S2). Ten out of 52 PER DEGs were significantly down-regulated in MG8 vs CK8 (Mgl0266600, PER12), MG10 vs CK10 (Mgl0017650, PER4), MG12 vs CK12 (Mgl0266600; Mgl0004180, PER2; Mgl0037470, PER25; Mgl0017750, PER1; Mgl0017650, PER4), MG16 vs CK16 (Mgl0266600; Mgl0004180; Mgl0017210, PER12; Mgl0272220, PER5; Mgl0017750; Mgl0017760, PER70), and MG20 vs CK20 (Mgl0004180; Mgl0173900, PER51; Mgl0202410, PER66). PER19 (Mgl0077730) was initially significantly up-regulated in MG12 vs CK12, then decreased significantly in MG20 vs CK20, while two PER5 (Mgl0272180 and Mgl0272230) were down-regulated in MG16 vs CK16 then up-regulated in MG20 vs CK20. Thirty-nine out of 52 PER DEGs were up-regulated at different stages, 10 of which were exclusively significantly up-regulated (Mgl0272300, PER5; Mgl0017540, PER4; Mgl0016520, PER27; Mgl0184200, PER1; Mgl0125520, PER3; Mgl0167670, PER56; Mgl0294200, PER5; Mgl0016070, PER3; Mgl0167660, PER27; Mgl0219990, PER3) in MG12 vs CK12, nine (Mgl0017510, PER4; Mgl0141010, PER39; Mgl0008240, PER47; Mgl0114430, PER2; Mgl0016530, PER56; Mgl0141000, PER3; Mgl0056040, PER51; Mgl0121310, PER55; Mgl0125510, PER3) in MG16 vs CK16, and five (Mgl0263220, PER5; Mgl0017630, PER52; Mgl0017200, PER12; Mgl0272210, PER5; Mgl0017790, PER4) in MG20 vs CK20. Three PER DEGs (Mgl0017780, PER1; Mgl0016080, PER3; Mgl0272310, PER5) were significantly up-regulated at 8 d, and Mgl0272310 (PER5) and Mgl0016080 (PER3) were significantly up-regulated at 10 and 12 d, respectively. Seven PER DEGs (Mgl0068470, PER5; Mgl0017520, PER4; Mgl0305250, PER1; Mgl0134470, PER27; Mgl0037420, PER28; Mgl0178090, PER4; Mgl0016510, PER39) were significantly up-regulated at 12 and 16 d in MG, relative to CK. Mgl0016480 (PER39) and Mgl0178150 (PER4) were significantly up-regulated in MG8 vs CK8, MG12 vs CK12 and MG16 vs CK16. Similarly, Mgl0017890 (PER1) was significantly up-regulated in MG16 vs CK16 and MG20 vs CK20. Mg10016500 (PER56) was significantly up-regulated at 10, 12 and 16 d in MG, compared with CK. Mgl0282530 (PER72) was significantly up-regulated from 8 to 16 d in MG vs CK.
In addition, one DEG for an aldehyde dehydrogenase family 2 member B4 gene (ALDH2B4), one caffeic acid 3-O-methyltransferase gene (COMT1), one trans-cinnamate 4-monooxygenase C4H2 gene (C4H2), one acyltransferase GLAUCE (GLC) gene, and two cinnamoyl-CoA reductase 1 (CCR) genes (Fig. 6, Table S2) were enriched in the “phenylpropanoid biosynthesis” pathway. ALDH2B4 (Mgl0101560) was significantly up-regulated in MG16 vs CK16. C4H2 (Mgl0012640) was significantly up-regulated in MG8 vs CK8 and MG20 vs CK20, and down-regulated in MG12 vs CK12. GLC (Mgl0271670) was significantly down-regulated in MG16 vs CK16. CCR (Mgl0131050 and Mgl0131240) was significantly up-regulated in MG8 vs CK8, and COMT1 (Mgl0247750) in MG12 vs CK12.
DEGs enriched in tryptophan metabolism, diterpenoid biosynthesis and zeatin biosynthesis pathways
Four DEGs were enriched in the “tryptophan metabolism” pathway, which is involved in IAA biosynthesis, including genes encoding proteins for indole-3-pyruvate monooxygenase (YUC2) (Mgl0088070), acetylserotonin O-methyltransferase (ASMT) (Mgl0116960 and Mgl0209550), and COMT1 (Mgl0247750) (Fig. 7a, Table S3). Mgl0088070 was significantly down-regulated in MG12 vs CK12. Mgl0116960 (ASMT) was significantly up-regulated in MG16 vs CK16 whereas Mgl0209550 (ASMT) was significantly down-regulated in MG10 vs CK10. One DEG (Mgl0229440, ent-kaur-16-ene synthase, KS1) involved in gibberellin biosynthesis was enriched in the “diterpenoid biosynthesis” pathway, and was significantly up-regulated in MG16 vs CK16 (Fig. 7b, Table S3).
One DEG related to zeatin O-xylosyltransferase (ZOX1) (Mgl0264320), two cis-zeatin O-glucosyltransferase 1 (CISZOG1) DEGs (Mgl0049820 and Mgl0049830) and three scopoletin glucosyltransferase (TOGT1) DEGs (Mgl0307290, Mgl0063030 and Mgl0062950) were enriched in the “zeatin biosynthesis” pathway; these regulated the biosynthesis of TZT, TZR, IP and IPA (Fig. 7c and Table S3). ZOX1 (Mgl0264320) was significantly down-regulated in MG8 vs CK8 and MG16 vs CK16. Two CISZOG1 were significantly up-regulated in MG12 vs CK12, and Mgl0049830 (CISZOG1) was also significantly up-regulated in MG10 vs CK10. Three TOGT1 DEGs showed considerable differential expression patterns in different comparison groups: Mgl0307290 (TOGT1) was up-regulated in MG16 vs CK16; Mgl0063030 (TOGT1) was up-regulated in MG12 vs CK12 and down-regulated in MG20 vs CK20; Mgl0062950 (TOGT1) was up-regulated in MG16 vs CK16 and MG20 vs CK20.
Transcript factors (TFs) analysis
A total of 31 DEGs were identified as TF, and were classified into 13 families with reference to the PlantTFDB, including 9 MYB (v-MYB avian myeloblastosis viral oncogene homolog), 4 NAC (NAM/ATAF/CUC protein), 3 AP2 (APETALA2), 3 bHLH (basic helix-loop-helix domain protein), 2 ERF (Ethylene responses factor), 2 MYB-related, 2 WRKY, 1 GRAS (Glycopeptide resistance-associated protein), 1 HB-other (Homeobox-other), 1 HSF (Heat stress transcription factor), 1 LBD (Lateral organ boundaries domain protein), 1 YABBY and 1 bZIP (basic leucine zipper domain protein) (Fig. S3 and Table S4) TF. Most of those DEGs were not annotated in a certain pathway, except the bZIP. Among those DEGs, eight out of nine MYB were significantly up-regulated in MG8 vs CK8, MG12 vs CK12, MG16 vs CK16 and MG20 vs CK20 group, respectively, indicating a potential positive role in meristematic cells and AR formation, and AR elongation of M. glyptostroboides.
Co-expression network analysis by weighted gene co-expression network analysis (WGCNA)
To explore the relationship between phytohormone content and gene expression in M. glyptostroboides AR formation, a co-expression network was constructed by WGCNA. The modules in WGCNA are defined as clusters of highly interconnected genes, and genes in the same cluster have a high correlation coefficient. Five modules had a highly positive correlation with phytohormone content. Among them, the correlation coefficient between MEblue with TZT and GA1 was high (0.737 and 0.771, respectively) (Fig. 8a). MEturquoise had a highly positive correlation (R = 0.565, p = 0.00) with IAA, but had highly negative correlations with TZT (R = -0.856, p = 0.0) and GA1 (R = -0.841, p = 0.00) (Fig. 8a). MEred had a highly positive correlation with GA4 (Fig. 8a). MEbrown had positive correlations with IP, TZT, TZR, GA1 and GA3 (Fig. 8a). MEblack had positive correlations with IP, IPA, TZT, TZR, GA1, GA3 and GA7 (Fig. 8a). According to the related phytohormone components in other modules and correlations, MEblack, MEblue and MEturquoise were selected for network model analysis.
The genes with top 500 connectivity were selected to perform further analysis. The hub genes are those that show high connectivity in a network. In the MEblack module, 118 genes were identified with high expression in the control group (Fig. 8b). Genes encoding Receptor-like protein EIX2, L-gulonolactone oxidase 5, α-amylase, glutathione S-transferase U7, endochitinase 4 and xanthotoxin 5-hydroxylase CYP82C4 were identified as candidate hub genes for this module. Those genes were involved in “plant-pathogen interaction”, “ascorbate and aldarate metabolism”, “biosynthesis of cofactors”, “starch and sucrose metabolism”, glutathione metabolism”, “amino sugar and nucleotide sugar metabolism” and “isoflavonoid biosynthesis” pathways, respectively (Table 2). A total of 1080 genes with high expression at the AR elongation stage (CK20, MG16, MG20) (Fig. 8c) were categorized into the MEblue module, and several genes involved in “phenylpropanoid biosynthesis”, “plant hormone signal transduction”, “MAPK signaling pathway - plant”, “plant-pathogen interaction”, “tryptophan metabolism” and “zeatin biosynthesis” pathways, including PER, PR1, XTH, ASMT and CISZOG1, were identified as candidate hub genes for this module (Table 2). A total of 1467 genes with high expression at the AR initiation phase (MG8) or the stages in which no AR primordia with potential meristematic cells were observed (CK8, CK10) (Fig. 8d) were identified in the MEturquoise module. Three candidate hub genes (PER, GLC, KS1) were involved in “phenylpropanoid biosynthesis” and “diterpenoid biosynthesis” pathways (Table 2).
Reverse transcription quantitative polymerase chain reaction (RT-qPCR) analysis of gene expression
To further validate the RNA-seq results, eight candidate DEGs (four DEGs (Mgl0010460, IAA27; two XTH1, Mgl0287720 and Mgl0287740; Mgl0289590, XTH5) enriched in the “plant hormone signal transduction” pathway; two DEGs (Mgl0141010, PER39; Mgl0167670, PER56) enriched in the “phenylpropanoid biosynthesis” pathway; TOGT1 (Mgl0063030) and CISZOG1 (Mgl0049820) enriched in the “zeatin biosynthesis” pathway) were selected for RT-qPCR analysis of M. glyptostroboides stems that were collected from 2/5-strength WPM and control medium at 8, 10, 12, 16 and 20 d. The expression trends of the DEGs from qRT-PCR and RNA-seq analysis were largely consistent (Fig. 9). In addition, correlation analysis showed that the RNA-seq results were highly consistent (R2 = 0.8083) with the RT-qPCR results (Fig. 9), demonstrating that the transcriptome data was reliable to reflect the response of 2/5-strength WPM-induced AR formation of M. glyptostroboides plantlets.