Abstract

Background

   The extract from Metasequoia glyptostroboides Hu et Cheng, a rare and
   endangered species native to China, exhibits numerous biological and
   pharmacological activities. The species is recalcitrant to rooting
   during micropropagation, a challenge that has yet to be resolved. In
   this study, transcriptomic and hormonal analyses were conducted to
   appreciate the molecular mechanism of adventitious root (AR) formation
   in optimized rooting conditions.

Results

   The use of 2/5-strength Woody Plant Medium (WPM) significantly promoted
   AR formation of M. glyptostroboides shoots while the content of
   endogenous auxin, cytokinins and gibberellins (GAs) varied at different
   stages of AR formation. Transcriptomic analysis showed the significant
   up- or down-regulation of differentially expressed genes (DEGs)
   associated with plant hormone signal transduction and the
   phenylpropanoid biosynthesis pathway in response to 2/5-strength WPM.
   DEGs related to the biosynthesis of indole-3-acetic acid, cytokinins
   and GAs were identified. Transcript factors involved in 13 families
   were also revealed. A weighted gene co-expression network analysis
   indicated a strong correlation between hormones and genes involved in
   plant hormone signal transduction and the phenylpropanoid biosynthetic
   pathway.

Conclusions

   These results indicate that the AR-promoting potential of 2/5-strength
   WPM in M. glyptostroboides was due to complex interactions between
   hormones and the expression of genes related to plant hormone signal
   transduction and the phenylpropanoid biosynthetic pathway.

Supplementary Information

   The online version contains supplementary material available at
   10.1186/s12864-024-10989-6.

   Keywords: Metasequoia glyptostroboides Hu et Cheng, Transcriptome,
   Endogenous hormones, Adventitious root

Background

   Roots are an essential element of plant growth as they provide water
   and nutrients while providing structural support. Adventitious roots
   (ARs), which form from non-root tissue post-embryonically, are crucial
   for the establishment of a rooting system for plants [[40]1, [41]2],
   and are the main developmental path by which plantlets form roots
   during vegetative propagation [[42]3]. ARs can form during vegetative
   propagation with or without stimulation by plant growth regulators
   (PGRs), typically auxins [[43]4, [44]5]. ARs also enhance plants’
   adaptation to adverse environmental conditions and abiotic stresses
   [[45]6–[46]8]. In vitro propagation is an important technology for
   plant conservation because it allows for the mass propagation and
   regeneration of elite germplasm. However, it is not uncommon to observe
   recalcitrance in rooting, especially of woody plants, as has been
   observed for Berberis aristata DC. [[47]9], magnolia [[48]10],
   Cariniana legalis [[49]11], and Pinus koraiensis Sieb. et Zucc.
   [[50]12]. Therefore, the ability to induce ARs is a key objective for
   improving the recalcitrance of in vitro rooting in difficult-to-root
   species.

   Several factors can optimize AR induction, including PGRs type and
   concentration [[51]13], light [[52]14], temperature [[53]15], and other
   factors. Nutrient deficiency can also serve as a cue for the initiation
   and elongation of ARs. For example, low phosphorus content increased AR
   formation in Solanum lycopersicum [[54]16], Populus ussuriensis
   [[55]17], and Malus domestica [[56]18]. The deficiency of Mn or Fe
   during the AR induction phase promoted more and longer roots in
   Eucalyptus globulus cuttings [[57]19]. While previous studies on
   nutrient deficiency-induced AR development centered on a single
   nutrient, little research has been conducted on the impact of multiple
   nutrient deficiency. In vitro propagation often involves a complex
   mixture of macro- and micro-nutrients, in the form of a basal medium,
   whose strength can impact the induction of ARs. Therefore, the
   adjustment of basal medium strength is a simple but effective way to
   adjust the molar concentrations of several nutrients simultaneously.
   For example, the use of different strengths of a popular basal medium,
   Murashige and Skoog (MS), showed an AR formation gradient from
   Echinodorus ‘Indian Red’ shoots: 1/2 > 3/4 > 1/4 > full-strength
   [[58]20], or 1/2 > 3/4 > full > 1/4-strength for Ledebouria revoluta
   bulblets [[59]21]. Half-strength Woody Plant Medium (WPM) [[60]22] free
   of PGRs induced significantly more ARs from Salix tetrasperma shoots
   than full-strength WPM [[61]23].

   The accumulation and distribution of endogenous hormones in plants tend
   to change in response to exogenous hormones or other factors, resulting
   in AR formation. In Malus hupehensis (Pamp.) Rehd., stem cuttings
   exogenously treated with indole-3-acetic acid (IAA) significantly
   accelerated the synthesis of endogenous IAA, gibberellic acid 3
   (GA[3]), and zeatin riboside (ZR), sped up the metabolism of
   carbohydrates and soluble proteins, leading to more rapid AR formation
   [[62]24]. After treating Vaccinium corymbosum L. cuttings with
   indole-3-butyric acid (IBA), endogenous IAA increased significantly at
   the root primordium formation stage whereas GA[3] content did not
   change significantly, while the contents of abscisic acid (ABA) and ZR
   were significantly lower than the control group, suggesting that higher
   IAA and lower ABA and zeatin contents might facilitate AR formation and
   development [[63]25]. In stem micro-cuttings of Malus prunifolia var.
   ‘Ringo’, an apple rootstock, darkness promoted AR formation by
   increasing or suppressing the endogenous content of IAA, ZR, ABA,
   GA[3], brassinolide and jasmonic acid (JA) at different stages
   [[64]26]. In summary, the ability to regulate endogenous hormone
   content is an essential strategy in the control of AR formation.

   Alongside the fluctuation in endogenous hormones, various key genes,
   proteins and pathways are involved in AR development [[65]27–[66]30].
   Transcriptome sequencing revealed that most candidate genes involved in
   AR formation encoded polar auxin transport carriers, auxin-regulated
   proteins, and cell wall remodeling enzymes, leading to an enrichment of
   the biosynthesis of secondary metabolites, plant hormone signal
   transduction, starch and sucrose metabolism, and other metabolic
   processes [[67]31–[68]34]. In the absence of PGRs and other growth
   regulators, genes related to the synthesis, transport, metabolism and
   recognition of plant hormones were involved in AR formation from
   Populus euramericana micro-cuttings [[69]35]. Red light promoted AR
   formation in Camellia gymnogyna Chang cuttings due to a complex
   interaction of hormones and regulation of the expression of genes
   participating in the response to hormones, cellular glucan metabolism,
   auxin, and other metabolic processes [[70]36]. Despite these advances
   in understanding basic and more complex molecular mechanisms underlying
   AR induction and formation, in vitro rooting remains recalcitrant in
   some woody species.

   Metasequoia glyptostroboides Hu et Cheng (Cupressaceae), a known living
   fossil endemic to southeast China, is a large deciduous conifer. Living
   individuals of M. glyptostroboides were first discovered in the early
   1940s in western Sichuan, China, and subsequent population studies
   revealed that the tree only grows in Sichuan, Hunan, and Hubei
   provinces [[71]37, [72]38]. However, given the importance and
   socio-economic value of its wood, M. glyptostroboides was successfully
   cultivated in Asia, Africa, Europe, and North America [[73]39]. M.
   glyptostroboides has remained genetically unchanged for millions of
   years and survived extreme ecological and climatic changes since the
   Cretaceous period [[74]40]. Extracts from M. glyptostroboides displayed
   numerous biological and pharmacological activities, including
   antioxidant, antidermatophytic [[75]41], antifungal [[76]42],
   anticancer [[77]43], and anti-amyloidogenic [[78]44]. Due to the
   limited distribution and declining number of native individual trees
   found in natural stands, M. glyptostroboides was listed as a Class I
   rare and endangered species in the List of Chinese Rare and Endangered
   Plants about 30 years ago [[79]45]. It is also classified as an
   endangered species on the IUCN Red List of Threatened Species [[80]46].
   M. glyptostroboides seeds are scarce due to the occurrence of male and
   female cones on different branches of adult trees, and a long juvenile
   phase before female buds appear [[81]47]. Therefore, the propagation of
   M. glyptostroboides is the main limiting factor for further research
   development of this living fossil tree.

   In a previous study, we developed a shoot-based micropropagation
   protocol for M. glyptostroboides. When treated with either IAA, IBA or
   1-naphthaleneacetic acid (NAA), M. glyptostroboides shoots showed
   limited rooting and root number on full-strength WPM after culture for
   2 months [[82]48]. Thus, the objectives of this study were to optimize
   AR induction of M. glyptostroboides shoots and attempt to reveal the
   mechanism underlying AR formation using transcriptomic and hormone
   analyses.

Results

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. [83]1). The rooting percentage (Fig. [84]1b), average root number
   (Fig. [85]1c) and root length (Fig. [86]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. [87]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.

Fig. 1.

   [88]Fig. 1
   [89]Open in a new tab

   Adventitious root induction of Metasequoia glyptostroboides shoots
   during in vitro rooting after 30 d (a, phenotype of adventitious
   rooting in different strength of WPM medium; b, rooting percentage; c,
   average root number; d, average root length) treatment. Red
   bars = 1 cm. Treatments (X-axis) indicate the strength of WPM medium.
   Bars indicate means ± SD. Different letters indicate statistically
   significant differences between the groups based on Duncan's multiple
   range test (p < 0.05) for the designated treatments. Forty shoots were
   prepared for each treatment, and three replicates were performed for
   each treatment. Graph was generated by OriginPro 2021

Fig. 2.

   Fig. 2
   [90]Open in a new tab

   Adventitious root development of Metasequoia glyptostroboides shoots
   stem base (0.5–1.0 cm) after 8, 10, 12, 16 and 20 d in 2/5 (MG) and
   full-strength (CK) WPM treatment during in vitro rooting. Red bars in
   photos representing phenotypic changes = 1 cm. Red bars in photos of
   histological analyses = 0.1 mm. Red arrows indicate meristematic cells
   (8 d in MG), the root primordium (10 d in MG; 12 d in CK), and
   adventitious roots (12 d in MG; 16 d in CK)

Analysis of endogenous hormone content

   Nine endogenous hormones, IAA, isopentenyladenine (IP),
   isopentenyladenosine (IPA), trans-zeatin riboside (TZR), trans-zeatin
   (TZT), GA[1], GA[3], GA[4] and GA[7], at 8, 10, 12, 16, and 20 d in MG
   and CK were detected in the stem bases of M. glyptostroboides
   (Fig. [91]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 16 d.
   TZT, GA[1] and GA[3] accumulated rapidly at 16 and 20 d in MG and CK
   (GA[1] was not detected at 8–12 d). TZT content in MG was significantly
   higher than in CK at 20 d, whereas GA[1] and GA[3] contents were
   significantly lower in MG at 16 and 20 d. GA[4] 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 GA[4] at 16 and 20 d, corresponding to
   the timing of AR primordia formation and breaking through the
   epidermis. Compared with CK, GA[7] 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.

Fig. 3.

   [92]Fig. 3
   [93]Open in a new tab

   Determination of endogenous IAA, IP, IPA, TZT, TZR, GA[1], GA[3], GA[4]
   and GA[7] content at 8, 10, 12, 16, and 20 d in 2/5 (MG) and
   full-strength (CK) WPM treatment in stem bases of Metasequoia
   glyptostroboides. Bars indicate means ± SD. Asterisk indicates
   significant differences following one-way ANOVA and according to the
   Mann Whitney non-parametric test (p < 0.05) in comparisons of MG with
   CK at designed time points. Graph was generated by OriginPro 2021. The
   green color indicated MG group, and the orange color indicate CK group

   In brief, 2/5-strength WPM induced AR formation from M.
   glyptostroboides stems, displaying a lower endogenous content of IAA,
   TZT, GA[1], GA[3] and GA[7] 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 GA[4] content at 12 d. At the
   early stage of AR elongation (16 d), stems showed a lower content of
   GA[4] 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 20 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 [94]1). After filtering, the
   number of clean reads per library ranged from 48,871,938 to 69,828,774
   (Table [95]1).

Table 1.

   Quality of the Metasequoia glyptostroboides transcriptome
   Sample Raw reads Raw bases Clean reads Clean bases Error rate (%) Q20
   (%) Q30 (%) GC content (%) Total mapped Multiple mapped Uniquely mapped
   CK8_1 60,082,318 9,072,430,018 59,774,080 8,769,466,610 0.0248 98.23
   94.34 45.27 55,333,828(92.57%) 4,554,573(7.62%) 50,779,255(84.95%)
   CK8_2 58,119,546 8,776,051,446 57,774,384 8,450,597,743 0.0257 97.86
   93.43 44.78 53,178,728(92.05%) 2,866,686(4.96%) 50,312,042(87.08%)
   CK8_3 57,664,460 8,707,333,460 57,361,744 8,446,704,536 0.0248 98.24
   94.29 45.03 53,061,066(92.5%) 3,151,893(5.49%) 49,909,173(87.01%)
   CK10_1 57,155,996 8,630,555,396 56,851,942 8,363,762,156 0.0245 98.3
   94.57 45.01 52,573,479(92.47%) 3,064,712(5.39%) 49,508,767(87.08%)
   CK10_2 56,157,598 8,479,797,298 55,862,822 8,225,946,595 0.0249 98.16
   94.17 44.94 51,601,975(92.37%) 2,942,352(5.27%) 48,659,623(87.11%)
   CK10_3 54,395,208 8,213,676,408 54,110,662 7,959,165,645 0.0252 98.06
   93.92 44.8 49,835,785(92.1%) 2,600,890(4.81%) 47,234,895(87.29%)
   CK12_1 65,305,512 9,861,132,312 64,905,262 9,527,758,039 0.0249 98.15
   94.17 45 59,959,450(92.38%) 3,632,263(5.6%) 56,327,187(86.78%)
   CK12_2 53,578,760 8,090,392,760 53,313,598 7,821,609,247 0.0248 98.24
   94.33 44.97 49,315,718(92.5%) 3,049,087(5.72%) 46,266,631(86.78%)
   CK12_3 65,361,464 9,869,581,064 64,952,778 9,537,325,409 0.0253 98
   93.78 45.13 59,959,417(92.31%) 3,812,965(5.87%) 56,146,452(86.44%)
   CK16_1 49,398,180 7,459,125,180 49,099,194 7,197,946,881 0.0258 97.85
   93.29 44.87 45,391,609(92.45%) 2,435,551(4.96%) 42,956,058(87.49%)
   CK16_2 70,131,470 10,589,851,970 69,828,774 10,290,102,514 0.0245 98.34
   94.62 44.95 64,514,772(92.39%) 4,078,200(5.84%) 60,436,572(86.55%)
   CK16_3 55,947,022 8,448,000,322 55,610,250 8,183,099,410 0.0246 98.26
   94.48 44.98 51,456,751(92.53%) 3,269,546(5.88%) 48,187,205(86.65%)
   CK20_1 58,068,412 8,768,330,212 57,800,252 8,474,341,705 0.0248 98.23
   94.35 44.84 53,401,925(92.39%) 2,947,942(5.1%) 50,453,983(87.29%)
   CK20_2 58,827,058 8,882,885,758 58,530,182 8,605,083,535 0.0249 98.16
   94.19 44.95 54,079,930(92.4%) 3,275,795(5.6%) 50,804,135(86.8%)
   CK20_3 63,537,498 9,594,162,198 63,192,082 9,281,162,149 0.0248 98.22
   94.33 44.7 58,509,070(92.59%) 2,973,012(4.7%) 55,536,058(87.88%)
   MG8_1 49,120,012 7,417,121,812 48,871,938 7,187,278,001 0.0246 98.28
   94.51 45.05 45,149,064(92.38%) 2,678,993(5.48%) 42,470,071(86.9%)
   MG8_2 56,827,144 8,580,898,744 55,167,482 7,052,330,992 0.0267 97.32
   93.17 44.3 49,156,124(89.1%) 3,242,379(5.88%) 45,913,745(83.23%)
   MG8_3 58,709,874 8,865,190,974 58,429,364 8,606,828,896 0.0249 98.15
   94.17 45.11 53,848,855(92.16%) 3,663,939(6.27%) 50,184,916(85.89%)
   MG10_1 52,329,230 7,901,713,730 52,093,458 7,656,876,244 0.0246 98.3
   94.56 44.77 48,098,250(92.33%) 2,698,678(5.18%) 45,399,572(87.15%)
   MG10_2 54,553,552 8,237,586,352 54,320,766 7,999,361,115 0.0245 98.31
   94.56 44.9 50,234,532(92.48%) 2,981,073(5.49%) 47,253,459(86.99%)
   MG10_3 55,617,404 8,398,228,004 55,327,556 8,130,139,634 0.0247 98.23
   94.39 44.79 51,112,215(92.38%) 2,707,780(4.89%) 48,404,435(87.49%)
   MG12_1 62,748,826 9,475,072,726 62,389,216 9,075,067,389 0.0249 98.19
   94.29 44.97 57,730,219(92.53%) 3,621,721(5.81%) 54,108,498(86.73%)
   MG12_2 53,627,786 8,097,795,686 53,360,844 7,765,546,534 0.0244 98.39
   94.76 44.91 49,498,029(92.76%) 2,834,986(5.31%) 46,663,043(87.45%)
   MG12_3 54,641,228 8,250,825,428 54,339,608 7,978,760,330 0.0252 98.07
   93.97 44.92 50,285,195(92.54%) 2,873,572(5.29%) 47,411,623(87.25%)
   MG16_1 55,040,448 8,311,107,648 54,724,378 8,037,510,329 0.0247 98.27
   94.43 44.67 50,671,783(92.59%) 2,631,513(4.81%) 48,040,270(87.79%)
   MG16_2 57,240,268 8,643,280,468 56,916,708 8,367,806,053 0.0246 98.27
   94.5 44.82 52,756,590(92.69%) 2,816,102(4.95%) 49,940,488(87.74%)
   MG16_3 49,161,246 7,423,348,146 48,896,832 7,193,153,938 0.025 98.13
   94.11 44.66 45,220,062(92.48%) 2,333,811(4.77%) 42,886,251(87.71%)
   MG20_1 55,207,152 8,336,279,952 54,882,790 8,086,450,953 0.025 98.12
   94.08 44.83 50,320,178(91.69%) 2,927,434(5.33%) 47,392,744(86.35%)
   MG20_2 54,752,716 8,267,660,116 54,443,578 7,998,713,968 0.0248 98.2
   94.26 44.67 50,192,321(92.19%) 2,861,494(5.26%) 47,330,827(86.94%)
   MG20_3 57,440,612 8,673,532,412 57,096,700 8,388,980,293 0.0249 98.18
   94.26 44.88 52,781,360(92.44%) 3,447,756(6.04%) 49,333,604(86.4%)
   [96]Open in a new tab

   Note 1: Raw reads = Total number of reads in the raw sequencing data.
   Raw bases = Total raw sequencing data. Clean reads = Total number of
   reads of sequencing data after quality control. Clean bases = Total
   sequencing data after quality control. Error rate (%) = Average error
   rate of sequencing bases corresponding to quality control data. Q20 (%)
   and Q30 (%) refer to the percentage of bases with sequencing
   quality > 99% and 99.9%, respectively in the total base after quality
   control. GC content (%) = The percentage of the sum of G and C bases
   corresponding to the quality control data in the total base. Total
   mapped = Number of clean reads located on the genome. Multiple
   mapped = Number of clean reads with multiple alignment locations on the
   reference sequence. Uniquely mapped = Number of clean reads with unique
   alignment positions on the reference sequence

   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. [97]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. [98]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. [99]4b).

Fig. 4.

   [100]Fig. 4
   [101]Open in a new tab

   DEGs identified in different comparisons during adventitious root
   development of Metasequoia glyptostroboides. a Number of up- and
   down-regulated DEGs at 8, 10, 12, 16 and 20 d compared 2/5-strength WPM
   treatment (MG) with the control (CK). The X-axis represents the groups
   that were compared during rooting, while the Y-axis represents the
   number of DEGs. b Venn diagram of DEGs (absolute numbers and relative
   percentage values in parentheses below) between different comparisons

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. [102]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).

Fig. 5.

   Fig. 5
   [103]Open in a new tab

   Expression analysis of DEGs enriched in the “plant hormone signal
   transduction” pathway. Different colors indicate DEG expression levels
   based on log10 (FPKM) values. Blue represents down-regulated expression
   while red represents up-regulated expression

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. [104]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.

Fig. 6.

   [105]Fig. 6
   [106]Open in a new tab

   Expression analysis of DEGs enriched in the “phenylpropanoid
   biosynthesis” pathway. Different colors indicate DEG expression levels
   based on log10 (FPKM) values. Blue represents down-regulated expression
   while red represents up-regulated expression

   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. [107]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. [108]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. [109]7b, Table S3).

Fig. 7.

   [110]Fig. 7
   [111]Open in a new tab

   Expression analysis of DEGs enriched in “tryptophan metabolism” (a),
   “diterpenoid biosynthesis” (b) and “zeatin biosynthesis” (c) pathways.
   Different colors indicate DEG expression levels based on log10 (FPKM)
   values. Blue represents down-regulated expression while red represents
   up-regulated expression

   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. [112]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 GA[1] was high (R = 0.737 and
   R = 0.771, respectively) (Fig. [113]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 GA[1]
   (R = -0.841, p = 0.00) (Fig. [114]8a). MEred had a highly positive
   correlation with GA[4] (Fig. [115]8a). MEbrown had positive
   correlations with IP, TZT, TZR, GA[1] and GA[3] (Fig. [116]8a). MEblack
   had positive correlations with IP, IPA, TZT, TZR, GA[1], GA[3] and
   GA[7](Fig. [117]8a). According to the related phytohormone components
   in other modules and correlations, MEblack, MEblue and MEturquoise were
   selected for network model analysis.

Fig. 8.

   [118]Fig. 8
   [119]Open in a new tab

   WGCNA-based co-expression analysis of genes in AR formation of
   Metasequoia glyptostroboides stems. a Correlation analysis was carried
   out between the co-expression modules of genes and endogenous hormone
   content. The numbers in the column on the left indicate the number of
   genes in each module, and each set of data on the right indicates the
   correlation coefficient and significance p-value of the module with the
   phenotype (in parentheses). b-d Heatmaps of top 500 connectivity genes
   corresponding to the MEblack, MEblue and MEturquoise modules,
   respectively. The different colors indicate DEG expression levels based
   on log10 (FPKM) values. Blue represents down-regulated expression while
   red represents up-regulated expression

   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. [120]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 [121]2). A total of 1080 genes with high
   expression at the AR elongation stage (CK20, MG16, MG20) (Fig. [122]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 [123]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. [124]8d) were identified in the MEturquoise module. Three
   candidate hub genes (PER, GLC, KS1) were involved in “phenylpropanoid
   biosynthesis” and “diterpenoid biosynthesis” pathways (Table [125]2).

Table 2.

   Candidate hub genes in MEblack, MEblue and MEturquoise modules
   Module Gene id Gene name Connectivity KEGG pathway
   MEblue Mgl0016480 PER39 215.16 Phenylpropanoid biosynthesis
   MEblue Mgl0017520 PER4 200.46 Phenylpropanoid biosynthesis
   MEblue Mgl0017890 PER1 131.69 Phenylpropanoid biosynthesis
   MEblue Mgl0125510 PER3 120.20 Phenylpropanoid biosynthesis
   MEblue Mgl0141010 PER39 137.28 Phenylpropanoid biosynthesis
   MEblue Mgl0178090 PER4 220.23 Phenylpropanoid biosynthesis
   MEblue Mgl0178150 PER4 192.39 Phenylpropanoid biosynthesis
   MEblue Mgl0240770 PR1 193.44 Plant hormone signal transduction;MAPK
   signaling pathway—plant;Plant-pathogen interaction
   MEblue Mgl0287740 XTH1 92.85 Plant hormone signal transduction
   MEblue Mgl0289590 XTH5 126.04 Plant hormone signal transduction
   MEblue Mgl0097340 XTH34 117.56 Plant hormone signal transduction
   MEblue Mgl0116960 ASMT 116.90 Tryptophan metabolism
   MEblue Mgl0049820 CISZOG1 163.26 Zeatin biosynthesis pathway
   MEturquoise Mgl0017750 PER1 306.48 Phenylpropanoid biosynthesis
   MEturquoise Mgl0037470 PER25 301.31 Phenylpropanoid biosynthesis
   MEturquoise Mgl0271670 GLC 345.56 Phenylpropanoid biosynthesis
   MEturquoise Mgl0229440 KS1 259.03 Diterpenoid biosynthesis
   MEblack Mgl0299640 Receptor-like protein EIX2, EIX2 12.13
   Plant-pathogen interaction
   MEblack Mgl0072980 L-gulonolactone oxidase 5; GGLO5 11.13 Ascorbate and
   aldarate metabolism;;Biosynthesis of cofactors
   MEblack Mgl0074800 Alpha-amylase, AMYA 10.67 Starch and sucrose
   metabolism
   MEblack Mgl0060640 Glutathione S-transferase U7, GSTU7 8.45 Glutathione
   metabolism
   MEblack Mgl0231410 Endochitinase 4, CHI4 6.97 Amino sugar and
   nucleotide sugar metabolism
   MEblack Mgl0145540 Xanthotoxin 5-hydroxylase CYP82C4, C82C4 5.22
   Isoflavonoid biosynthesis
   [126]Open in a new tab

   Note 2: Connectivity is a concept that is similar to the degree of
   networking. The connectivity of each gene is the sum of the edge
   attributes of the connected genes, and a greater the connectivity value
   indicates the greater importance of the node in the network

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. Among those genes, four candidate genes (Mgl0287740,
   Mgl0289590, Mgl0141010, Mgl0049820) were also identified as hub genes
   by WGCNA. The expression trends of the DEGs from qRT-PCR and RNA-seq
   analysis were largely consistent (Fig. [127]9). In addition,
   correlation analysis showed that the RNA-seq results were highly
   consistent (R^2 = 0.8083) with the RT-qPCR results (Fig. [128]9),
   demonstrating that the transcriptome data was reliable to reflect the
   response of 2/5-strength WPM-induced AR formation of M.
   glyptostroboides plantlets.

Fig. 9.

   [129]Fig. 9
   [130]Open in a new tab

   Analysis of the fold changes of eight candidate genes expressed in
   Metasequoia glyptostroboides stems during adventitious root induction
   and formation, determined by RNA-seq (green) and RT-qPCR (orange). Bars
   indicate means ± SD. The X-axis represents the time point after each
   treatment while the Y-axis represents the log2 (fold change). In
   correlations analysis of the expression levels of the DEGs measured via
   RNA-seq and RT-qPCR, the X-axis represents log2 (fold change) in
   RT-qPCR and the Y-axis represents log2 (fold change) in RNA-seq. Graph
   was generated by OriginPro 2021

Discussion

   Rooting recalcitrance has been a major limitation in the vegetative
   propagation of elite plant germplasm [[131]49–[132]51]. In general,
   PGRs can effectively induce or promote AR formation, combined with or
   without other regulators, and are the dominant factor affecting AR
   induction in the clonal propagation of many plant species
   [[133]52–[134]54]. However, the inducible ability of PGRs is
   insufficient in some difficult-to-root species. The effect of nutrient
   deficiency on AR formation has been reported in many plant species,
   especially crops [[135]55, [136]56]. Nutrient deficiency-induced stress
   stimulated AR growth, by increasing the available root surface area and
   contributing to nutrient uptake [[137]3]. In the present study,
   low-strength WPM significantly promoted AR formation from M.
   glyptostroboides stems. Compared with full-strength WPM, after
   elongation on WPM medium, plantlets likely experienced nutrient
   deficiency-induced stress during low strength WPM-induced AR formation.
   Thus, 2/5-strength WPM-promoted AR formation in M. glyptostroboides
   plantlets might serve as a response to nutrient deficiency. The low
   cost and efficient rooting promoted by the use of low ionic strength of
   the basal culture medium is a candidate strategy for the induction and
   formation of ARs in the micropropagation systems of difficult-to root
   species.

   Phytohormones, as physiological signals that regulate plant growth and
   development, are very important for AR development. There have been
   many reports dedicated to the accumulation patterns of endogenous IAA,
   cytokinins and GAs during AR formation. It is widely accepted that IAA
   is the central player that induces AR in competent cells, and an
   increase in IAA content is always observed over the course of rooting
   [[138]57–[139]59], as a result of the expression of multiple genes
   involved in the maintenance of auxin homeostasis [[140]60–[141]62].
   Cytokinins, including IP, IPA, TZR and TZT, are thought to have an
   early positive regulative function at the stage of dedifferentiation
   during AR formation, whereas high concentrations of cytokinins act
   antagonistically to auxin during the induction phase [[142]58, [143]63,
   [144]64]. GAs may have an AR phase-dependent effect, being inhibitory
   to root induction and stimulatory to AR formation in some species
   [[145]52, [146]63]. For M. glyptostroboides plantlets, 2/5-strength
   WPM-promoted rooting displayed a lower level of endogenous IAA and
   cytokinins than the control. Similar findings were observed during
   adventitious rooting of Camellia sinensis cuttings [[147]65] and Picea
   abies seedlings [[148]66]. In addition, a lower level of GA[1], GA[3]
   and GA[7] over the entire rooting process of M. glyptostroboides
   plantlets was revealed in 2/5-strength WPM-promoted rooting, confirming
   the hypothesis that high levels of GA negatively regulate AR formation
   [[149]67], and suggesting species- or treatment-dependent phytohormone
   accumulation patterns during AR formation. In this study, during AR
   formation in M. glyptostroboides plantlets, genes in tryptophan
   metabolism, diterpenoid biosynthesis and zeatin biosynthesis pathways
   showed differential expression patterns, alongside a fluctuation in
   IAA, IP, IPA, TZR, TZT and GA content. These results suggest a module
   in which low-strength WPM induced ARs in M. glyptostroboides plantlets
   by regulating the level of phytohormones and the expression of
   phytohormone-related genes.

   Generally, AR formation is inseparable from the regulation of plant
   hormone signal transduction [[150]68–[151]70]. Multiple plant hormone
   signal transduction-related genes are involved in adventitious rooting.
   GH3, AUX/IAA and SAUR are early auxin response proteins that are active
   during AR formation. GH3 plays a crucial role in conjugating IAA to
   amino acids, and is critical in maintaining auxin homeostasis
   [[152]71–[153]73]. AUX/IAA always participates in the auxin signaling
   pathway by interacting with auxin response factor protein or other
   genes [[154]74–[155]76]. SAURs mediate the regulation of several
   aspects of plant growth and development, and have positive or negative
   effects on primary, lateral and adventitious root development
   [[156]77–[157]80]. In M. glyptostroboides, GH3, IAA, and SAUR DEGs were
   significantly enriched in plant hormone signal transduction, implying a
   link to AR induction of M. glyptostroboides plantlets. Furthermore,
   numerous XTH DEGs that were significantly enriched in plant hormone
   signal transduction were also differentially expressed. XTH, belonging
   to glycoside hydrolase family 16, is a key enzyme in plant cell wall
   remodeling and a cell-wall-loosening protein involved in cell expansion
   or cell-wall weakening [[158]81]. In marigold explants, methane-rich
   water and 2,4-epibrassinolide significantly promoted the activity of
   two major cell wall-relaxing factors, XTH and peroxidase, which in turn
   promoted AR formation [[159]82]. However, an understanding of the
   direct relationship between XTH expression and AR formation is still
   limited. XTH might be an excellent candidate for future
   biotechnological efforts to improve rooting in some plant species.

   Phenylpropanoid biosynthesis has been identified as a remarkable
   pathway of AR formation in several species [[160]70, [161]83, [162]84].
   Multiple genes, such as, 4-Coumaroyl Co-A ligase (4CL), C4H, PER and
   phenylalanine ammonia lyase (PAL), involved in phenylpropanoid
   biosynthesis have been reported vital roles in plant growth,
   development, adaptation against various stresses, etc.
   [[163]85–[164]87]. The expressions of those genes contributed the
   differential accumulation of multiple phenolic compounds, such as
   chlorogenic acid, caffeic acid and ferulic acid, that promoted cell
   division and provided energy for the rooting process [[165]88–[166]90].
   Meanwhile, the negative role of some members of the above genes in root
   growth was also reported. In transgenic Nicotiana benthamiana, the
   overexpression of Ocimum kilimandscharicum 4CL11 resulted in a rootless
   or reduced root growth phenotype, which are likely involved in
   inhibition of auxin transport and signaling [[167]91]. In this study,
   phenylpropanoid biosynthesis was also analyzed as a significantly
   enriched pathway during the stage of AR emergence from the epidermis
   and elongation, with a number of PER DEGs having been identified. PER
   encodes peroxidase (E1.11.1.7), which plays an important role in a
   number of physiological processes in plants, are heme-containing
   enzymes with catalytic action on diverse organic compounds, including
   IAA, and their activity has been used as a biochemical marker of
   rooting [[168]92]. In basal tissue of Saccharum spp. interspecific
   hybrids microshoot, the altered expression of PER genes involved in
   phenylpropanoid biosynthesis during AR formation might cause
   alterations in the production of lignin and many functional compounds
   involved in various processes associated with cell division and
   organogenesis [[169]70]. In cucumber hypocotyls, AAA ATPase
   domain-containing protein CsARN6.1, promoted waterlogging-triggered AR
   formation by activating the expression of CsPrx5, which is a
   waterlogging-responsive class-III peroxidase that enhanced hydrogen
   peroxide production and increased adventitious rooting [[170]93].
   Therefore, PER could be considered as a crucial DEG in AR formation.
   However, the information of promotive or negative role of genes
   involved in phenylpropanoid biosynthesis pathway on AR formation is
   still limited, which need further researches to illustrate.

   WGCNA is a useful approach for identifying sample-specific modules and
   candidate hub genes in plant growth and development [[171]94–[172]96].
   The hub genes in adventitious rooting have been revealed by WGCNA
   analysis for some plant species. In Camellia sinensis (L.) O. Kuntze
   cuttings, hub genes involved in NAA-mediated AR formation were
   identified by WGCNA, such as Helix-loop-helix DNA-binding Domain
   transcription factor, AP2 domain (AP2)-containing protein, SAUR, etc.
   [[173]65]. In C. gymnogyna Chang, WGCNA analysis identified three
   highly correlated modules and hub genes mainly participated in
   “response to hormone”, “cellular glucan metabolic progress” and
   “response to auxin” during AR development under red light [[174]36]. A
   total of nine co-expression modules were constructed in this study,
   including three modules that were significantly related to phytohormone
   content in AR development of M. glyptostroboides plantlets. The
   candidate hub genes were involved in plant hormone signal transduction
   and phenylpropanoid biosynthesis pathway. These findings suggest a
   complex interaction between hormones and the expression of genes
   related to the above two pathways in AR formation of M.
   glyptostroboides plantlets.

Conclusion

   In this research, the use of 2/5-strength WPM promoted AR development
   from M. glyptostroboides stems. The contents of endogenous auxins,
   cytokinins and GAs changed dramatically at different stages of AR
   formation. The AR-promoting potential of 2/5-strength WPM for M.
   glyptostroboides is due to complex interactions between hormones and
   the expression of genes related to plant hormone signal transduction
   and phenylpropanoid biosynthesis pathway. A hypothetical model of
   adventitious root formation of M. glyptostroboides plantlets based on
   our results is shown in Fig. [175]10. The knowledge gained from this
   study will help researchers better understand the molecular traits of
   low strength basic culture medium-based regulation of adventitious
   rooting of M. glyptostroboides plantlets.

Fig. 10.

   [176]Fig. 10
   [177]Open in a new tab

   A possible network of AR-promoting potential of 2/5-strength WPM for
   Metasequoia glyptostroboides. Red arrows indicated up-regulated DEGs,
   and blue arrows indicated down-regulated DEGs. Green arrows indicated
   an interaction between the corresponding pathway. Orange arrows
   indicated the decrease in hormone accumulation. The complex
   interactions between hormones biosynthesis and the expression of genes
   related to plant hormone signal transduction and phenylpropanoid
   biosynthesis pathway contributed to the AR formation of M.
   glyptostroboides induced by 2/5-strength WPM under in vitro culture

Methods

Adventitious root induction

   This study employed the micropropagation system and culture conditions
   for M. glyptostroboides that were previously established [[178]48]. In
   vitro plantlets were maintained and propagated on WPM supplemented with
   5.0 μM 6-benzyladenine (BA) (Solarbio, Beijing, China). The composition
   of WPM was described in previous study [[179]97]. Single shoots (1 cm
   in height) cut from multiple in vitro shoots were inoculated on
   PGR-free WPM for 30 d. ARs were not induced during this period. For the
   rooting treatment, after trimming 2 mm from the base of single shoots
   (each containing four leaves and two nodes) cultured on PGR-free WPM
   for 30 d, trimmed shoots were inoculated onto 0, 1/5, 2/5, 4/5 or
   full-strength WPM supplemented with 0.05 μM NAA (Macklin, Shanghai,
   China), 10 g/L sucrose and 0.5% (w/v) agar. The 0-strength WPM
   treatment contained only 10 g/L sucrose and 0.5% (w/v) agar. The
   full-strength WPM treatment was used as the control (CK). Ten shoots
   were placed in each jar (four jars per treatment). Three replicates
   were performed for each treatment (n = 12 jars; total = 120 shoots).
   Rooting percentage, average root number and root length were calculated
   for each treatment after 30 d. Significant differences in each
   parameter between treatments were analyzed by Duncan's multiple range
   test (p ≤ 0.05) after mean separation using one-way analysis of
   variance.

Histological analysis

   At least 15 bases of M. glyptostroboides shoots (0.5–1.0 cm) from the
   control and optimized rooting treatments were collected at 8, 10, 12,
   16, and 20 d and fixed in formalin/acetic acid/alcohol (50% ethanol:
   37% formalin: glacial acetic acid, 18:1:1, v/v/v) at 25 ± 2 °C
   for > 24 h, then treated with protocol described in a previous study
   [[180]52] with minor modifications. After dehydrating shoot bases in a
   70–100% alcohol series, fixed material was infiltrated with, and
   embedded in molten paraffin wax (Sinopharm, Beijing, China).
   Sections (4–6 μm thick) for microscopy were made with a Leica RM2255
   (Wetzlar, Germany) rotary microtome and stained in 0.5% fast green
   (Sinopharm). Micrographs were obtained with a Nikon Eclipse E200
   microscope (Nikon, Tokyo, Japan) and HQimage C630 digital camera
   (Hengqiao, Hangzhou, China).

Determination of endogenous hormones

   To analyze the content of endogenous hormones, the same material
   employed for histological analysis was used, and stored at -80℃. Three
   biological replicates (20 cut stem bases in each) were harvested (0.2 g
   fresh weight (FW) for each base) to assess endogenous IAA, GA[1],
   GA[3], GA[4], GA[7], TZR, TZT, IP, and IPA content by liquid
   chromatography-tandem mass spectrometry (LC–MS/MS). The quantification
   endogenous hormones contents were performed according to previously
   described methods [[181]98] with minor modifications by Wuhan ProNets
   Biotechnology Co, Ltd. (Wuhan, China). All standards for this analysis
   were purchased from Sigma (St. Louis, MI, USA), except GA[1] and GA[7],
   which were purchased from Toronto Research Chemicals (Toronto, Canada).
   The internal standard mixture was consisted of D-IAA and D-SA, with the
   concentration of 1 μg/mL. The samples were ground to a powder using
   liquid nitrogen, then 10 volumes of acetonitrile solution and 2 μL of
   internal standard solution were added to the powder. The mixture was
   incubated at 4 °C overnight, then centrifuged at 12,000 g for 5 min
   with high-speed refrigerated centrifuge (Ke Cheng H3-18KR, Hunan,
   China). Five volumes of acetonitrile solution were added to the
   precipitate, which was extracted twice, then supernatants were pooled.
   After purified with CNW C18 chromatographic packing (CNW, Duesseldorf,
   Germany), the supernatant was collected and concentrated until dry with
   a CV600 centrifugal concentrator (Jiaimu, Beijing, China). For
   analysis, the dried supernatant was suspended in 150 μL methanol and
   filtered through a 0.22 μm organic phase filter membrane (Jingteng,
   Tianjin, China). Hormones were measured by high-performance liquid
   chromatography (HPLC) (Agilent 1290; Santa Clara, CA, USA) with WATERS
   ACQUITY UPLC HSS T3 (2.1 × 100 mm, 1.8 µm) (Waters, MA, USA)
   chromatographic column and AB Qtrap6500 mass spectrometer (AB Sciex,
   Foster City, CA, USA). The analytical mobile phase of HPLC consisted of
   0.1% formic acid (mobile A) and acetonitrile (mobile B). The gradient
   began at 95% A and 5% B for 1.0 min, moving to 5% A and 95% B at
   8.0 min, then returning to 95% A and 5% B at 8.1 min, and remaining
   there to 10.0 min. The flow rate was 0.3 mL/min and the column
   temperature was set at 30 °C. The MS parameters are indicated in Table
   S5. Following mean separation using one-way ANOVA, significant
   differences in the content of endogenous hormones between treatments
   were assessed with the Mann Whitney non-parametric test (p ≤ 0.05). The
   chromatogram of IAA, GA1, GA3, GA4, GA7, TZR, TZT, IP, and IPA
   standards were shown in Figure S4. The linear regression of the
   standards curves was detailed in Table S6.

Library preparation, sequencing and read mapping

   The samples used to analyze the content of endogenous hormones were
   also employed for RNA-seq analysis. Total RNA was extracted from each
   sample using TRIzol® reagent according the manufacturer’s instructions
   (Invitrogen, Carlsbad, CA, USA). The quality of total RNA was
   determined by a 5300 Bioanalyser (Agilent) and quantified using
   NanoDrop 2000 (Thermo Fisher Scientific, Wilmington, NC, USA). RNA
   purification, reverse transcription, library construction and
   sequencing were performed by Shanghai Majorbio Bio-pharm Biotechnology
   Co., Ltd. (Shanghai, China) according to the manufacturer’s
   instructions (Illumina, San Diego, CA, USA). The paired-end RNA-seq
   sequencing library was sequenced with the Illumina NovaSeq 6000. Raw
   paired end reads were trimmed and quality control was achieved by fastp
   [[182]99] using default parameters. Clean reads were separately aligned
   to the reference M. glyptostroboides genome [[183]100] in orientation
   mode using HISAT2 software [[184]101]. The mapped reads of each sample
   were assembled by String Tie [[185]102] using a reference-based
   approach.

Differential expression analysis and functional enrichment

   The expression levels of genes were calculated as reads per kilobase of
   transcript sequence per million base pairs sequenced (FPKM). Gene
   abundance was quantified by RSEM
   ([186]http://deweylab.biostat.wisc.edu/rsem/) [[187]103]. DESeq2
   [[188]104] was used to perform differential expression analysis. DEGs
   with |log[2]FC|≥ 1 and a false discovery rate (FDR) ≤ 0.05 were
   considered as significantly differently expressed genes. In addition,
   functional-enrichment analysis, including the use of Gene Ontology (GO)
   and Kyoto Encyclopedia of Genes and Genomes (KEGG), was performed to
   identify which DEGs were significantly enriched, as assessed by GO
   terms and KEGG pathways, compared with the whole-transcriptome
   background (Bonferroni-correction; p ≤ 0.05). GO terms and KEGG pathway
   enrichment analysis were carried out by Goatools
   ([189]https://github.com/tanghaibao/Goatools) and KOBAS
   ([190]http://kobas.cbi.pku.edu.cn/home.do) [[191]105]. Identification
   of TF families was performed by using BLASTX similarity from the Plant
   Transcription Factor Database (PlantTFDB v5.0;
   [192]http://planttfdb.gao-lab.org/) with the default parameters. Gene
   co-expression networks were constructed with the WGCNA, v 1.63 package
   [[193]106]. The WGCNA parameters settings were: average FPKM > 1;
   variable coefficient > 0.1; soft threshold = 10; minimum module
   size = 30; merge cut height = 0.1; connection value (|KME|) > 0.5.

RT-qPCR analysis

   For RT-qPCR analysis, eight candidate DEGs in the enriched KEGG pathway
   were selected to validate the RNA-seq results. Gene-specific primers
   are listed in Table S7. The qTower^3 (Analytik Jena, Jena, Germany) and
   PerfectStart Green qPCR Supermix (TransGen Biotech, Beijing, China)
   were employed for RT-qPCR. Three biological replicates and three
   technical replicates were performed for each candidate gene. M.
   glyptostroboides actin was used as the internal control and the
   differential expression of candidate DEGs was analyzed by the 2^−ΔΔCt
   method [[194]107].

Supplementary Information

   [195]Supplementary Material 1.^ (2.6MB, zip)

Acknowledgements