Integrated Transcriptomics and Metabolomics Analysis Promotes the Understanding of Adventitious Root Formation in Eucommia ulmoides Oliver

As a primary approach to nutrient propagation for many woody plants, cutting roots is essential for the breeding and production of Eucommia ulmoides Oliver. In this study, hormone level, transcriptomics, and metabolomics analyses were performed on two E. ulmoides varieties with different adventitious root (AR) formation abilities. The higher JA level on the 0th day and the lower JA level on the 18th day promoted superior AR development. Several hub genes executed crucial roles in the crosstalk regulation of JA and other hormones, including F-box protein (EU012075), SAUR-like protein (EU0125382), LOB protein (EU0124232), AP2/ERF transcription factor (EU0128499), and CYP450 protein (EU0127354). Differentially expressed genes (DEGs) and metabolites of AR formation were enriched in phenylpropanoid biosynthesis, flavonoid biosynthesis, and isoflavonoid biosynthesis pathways. The up-regulated expression of PAL, CCR, CAD, DFR, and HIDH genes on the 18th day could contribute to AR formation. The 130 cis-acting lncRNAs had potential regulatory functions on hub genes in the module that significantly correlated with JA and DEGs in three metabolism pathways. These revealed key molecules, and vital pathways provided more comprehensive insight for the AR formation mechanism of E. ulmoides and other plants.


Introduction
Eucommia ulmoides Oliver is a perennial tree with important medicinal, industrial, and ecological values [1][2][3].E. ulmoides is an herbal medicine and edible oil resource and is also a natural high-quality rubber resource with great developmental potentials [4][5][6].Cutting propagation has the advantages of stable genetic characteristics, a short seedlingraising period, simple operation, and low cost [7].As a dioecious plant, the propagation of E. ulmoides mainly depends on cutting to speed up the growth process and maintain excellent characters in the genetic breeding.The formation of adventitious roots (AR) is closely related to the survival of cuttings [8,9].The improvement of AR formation ability is essential and will greatly promote the breeding and production of E. ulmoides.
Plant hormones play an indispensable role in AR formation [10,11].Jasmonic (JA) acts with auxin and other phytohormones through sophisticated crosstalk during AR development [12,13].JA might enhance auxin synthesis in the induction stage to promote AR formation [14].The reduction of AR numbers in JA-deficient cuttings indicated that JA played an active role in the AR formation of wild-type petunia [15].JA could improve the biomass, total phenolic content, and antioxidant activity of AR of Acmella radicans cultured in ahake flasks [16].Notably, some studies have denoted that low concentrations of methyl jasmonate (MeJA) and high concentrations of MeJA promoted and inhibited AR, respectively, which had opposite effects on AR formation [17,18].In addition, other plant hormones can also regulate AR development through a series of crosstalk with JA and/or auxin, including ethylene, abscisic acid (ABA), salicylic acid (SA), and so on.JA could also mutually interfere with ethylene signaling to facilitate auxin-induced AR in dark-grown Arabidopsis thaliana seedlings [19].In Boehmeria nivea L., JA could interact with ethylene to regulate the occurrence and number of AR [20].ABA regulated indole-3-acetic acid (IAA) oxidase activity and the AR development of mung bean seedlings under cadmium stress [21].Exogenous ABA inhibited the formation of AR at a dose of 0.25 µm or higher in Arabidopsis [22].Low-concentration SA promoted AR development and altered the architecture of the root apical meristem [23].The appropriate concentration of exogenous SA could improve IAA level and promote the AR formation of cucumber [24].
AR formation has been deeply studied in the model plants, including A. thaliana [44], Oryza [45], and Populus [46].In E. ulmoides, the mechanism of the systematic regulation of AR formation by various hormones and key molecules remains elusive.The release of the high-quality genome of E. ulmoides provided an opportunity to deeply analyze the E. ulmoides AR formation at the molecular level [47,48].In this study, physiological indexes, whole transcriptomics, and metabolomics during the AR development of E. ulmoides were comprehensively analyzed.The dynamic levels of IAA, ABA, JA, and SA hormones were examined.Weighted co-expression network analysis (WGCNA) was performed, combined with hormone traits and expressed genes.Differentially expressed genes (DEGs), differentially expressed metabolites (DEMs), and pathway enrichment were investigated by transcriptomic data and metabolomic data.The cis-acting long noncoding RNAs (lncR-NAs) with potential regulatory functions were identified.This study provides abundant information for further research on the molecular mechanism during AR formation and has an important significance for the breeding and application of new varieties of E. ulmoides.1A).This is consistent with our previous study on the average number of roots and average length of roots of H6 and H8 [8].The hormone dynamic changes of H6 and H8 displayed similar rising or falling patterns in three stages.However, compared to H8, H6 exhibited significant differences in hormone levels in its roots at each stage (Figure 1B-E).In the S1 stage, H6 exhibited significantly higher levels of SA and JA compared to H8, and H6 had a significantly lower level of ABA than H8.In the S2 stage, H6 showed significantly higher levels of IAA and ABA compared to H8, and H6 had significantly lower levels of SA and JA than H8.In the S3 stage, H6 exhibited a significantly higher level of ABA compared to H8, and H6 had significantly lower levels of IAA, SA, and JA than H8.For each variety, the level of JA decreased, while the levels of IAA and SA increased from the S1 stage to the S2 stage.
by transcriptomic data and metabolomic data.The cis-acting long noncoding RNAs (lncRNAs) with potential regulatory functions were identified.This study provides abundant information for further research on the molecular mechanism during AR formation and has an important significance for the breeding and application of new varieties of E. ulmoides.

Morphological Characteristics and Hormone Content during AR Development
E. ulmoides Huazhong 6 (H6) had more AR numbers compared to E. ulmoides Huazhong 8 (H8) across three developmental stages, indicating obvious morphological differences (Figure 1A).This is consistent with our previous study on the average number of roots and average length of roots of H6 and H8 [8].The hormone dynamic changes of H6 and H8 displayed similar rising or falling patterns in three stages.However, compared to H8, H6 exhibited significant differences in hormone levels in its roots at each stage (Figure 1B-E).In the S1 stage, H6 exhibited significantly higher levels of SA and JA compared to H8, and H6 had a significantly lower level of ABA than H8.In the S2 stage, H6 showed significantly higher levels of IAA and ABA compared to H8, and H6 had significantly lower levels of SA and JA than H8.In the S3 stage, H6 exhibited a significantly higher level of ABA compared to H8, and H6 had significantly lower levels of IAA, SA, and JA than H8.For each variety, the level of JA decreased, while the levels of IAA and SA increased from the S1 stage to the S2 stage.

The DEGs and KEGG Pathways Enrichment during AR Development
The whole transcriptome strand-specific RNA sequencing of cuttings generated a total of 205.60 Gb clean bases (Q30 > 94.09%), with at least 15.67 Gb from each sample (Table S1).After mapping the clean reads to the E. ulmoides genome, 34,974 transcripts were finally obtained.In the principal component analysis (PCA), the first principal component (PC1) explained 26.3% of the variation, and the second principal component (PC2) explained 17.43% of the variation (Figure 2A).The samples were separated by cultivars along PC2.Correlation analysis explained the strong correlation of sample repetition (Figure 2B).
The number of DEGs in two E. ulmoides varieties at three AR development stages were 6571 in H8S1 vs. H6S1, 6251 in H8S2 vs. H6S2, and 3538 in H8S3 vs. H6S3, with the highest number of DEGs discovered in the S1 stage (Figure 2C).In each stage, the number of up- (E) SA level.S1, cutting on the 0th day.S2, cutting on the 18th day.S3, cutting on the 32nd day.Different letters denote significant differences at p < 0.05 level.

The DEGs and KEGG Pathways Enrichment during AR Development
The whole transcriptome strand-specific RNA sequencing of cuttings generated a total of 205.60 Gb clean bases (Q30 > 94.09%), with at least 15.67 Gb from each sample (Table S1).After mapping the clean reads to the E. ulmoides genome, 34,974 transcripts were finally obtained.In the principal component analysis (PCA), the first principal component (PC1) explained 26.3% of the variation, and the second principal component (PC2) explained 17.43% of the variation (Figure 2A).The samples were separated by cultivars along PC2.Correlation analysis explained the strong correlation of sample repetition (Figure 2B).
The number of DEGs in two E. ulmoides varieties at three AR development stages were 6571 in H8S1 vs. H6S1, 6251 in H8S2 vs. H6S2, and 3538 in H8S3 vs. H6S3, with the highest number of DEGs discovered in the S1 stage (Figure 2C).In each stage, the number of up-regulated DEGs was greater than the number of down-regulated DEGs.DEGs were used for KEGG enrichment analysis, and several pathways that may play key roles in AR development were identified.The top 20 enriched metabolic pathways in each stage were displayed (Figure 2D-F).In H8S1 vs. H6S1, metabolic pathways such as phenylpropanoid biosynthesis, starch and sucrose metabolism, glycolysis/gluconeogenesis, flavonoid biosynthesis, and the pentose phosphate pathway were enriched (Figure 2D).In H8S2 vs. H6S2, metabolic pathways such as ribosome, phenylpropanoid biosynthesis, glycolysis/gluconeogenesis, ABC transporters, and cysteine and methionine metabolisms were enriched (Figure 2E).In H8S3 vs. H6S3, metabolic pathways such as plant-pathogen interaction, phenylpropanoid biosynthesis, ribosome, MAPK signaling pathway-plant, and starch and sucrose metabolisms were enriched (Figure 2F).Furthermore, phenylpropanoid biosynthesis and flavonoid biosynthesis were enriched in the S1, S2, and S3 stages, and isoflavonoid biosynthesis was enriched in both the S1 and S2 stages.

Construction and Analysis of the WGCNA
All expressed genes were divided into four modules containing 685, 985, 270, and 60 genes in blue, turquoise, brown, and gray, respectively (Figure 3A).The genes in the blue module were positively correlated with JA and negatively correlated with SA, while the genes in the turquoise module were negatively correlated with JA and positively correlated with SA (Figure 3B).The significantly relevant genes were demonstrated by gene significance and module membership in MEblue-JA (Figure 3C) and MEturquoise-JA (Figure 3D).Twenty hub genes were identified in the blue and turquoise modules, respectively (Figure 3E,F).A total of 21 orthologous genes of these E. ulmoides hub genes were found in A. thaliana (Table S2).Compared with H8, 19 hub genes in the blue module were significantly up-regulated in the S1 stage in H6 (Figure 4A,B).The eight hub genes in the blue module were differentially expressed in three stages of AR development of H6, namely EU0109382, EU0106969, EU0100375, EU0120751, EU0109584, EU0116878, EU0125382, and EU0109541 (Figure 4A,B).The results of the functional annotation of proteins exposed that EU0120751 was an F-box protein and that EU0125382 was a SAUR-like auxin-responsive family protein.Compared with H8, 19 hub genes in the turquoise module were differentially expressed in H6, most of which were significantly down-regulated in the S1 stage and up-regulated in the S2 stage (Figure 4C,D).A total of 18 hub genes in the turquoise module were differentially expressed, containing EU0112695, EU0116412, EU0128499, EU0100696, EU0114710, EU0124232, EU0117900, EU0127605, EU0107185, EU0119332, EU0127354, EU0110081, EU0113433, EU0117354, EU0131180, EU0104849, EU0125430, and EU0102912 (Figure 4C,D).The results of the functional annotation of proteins indicated that EU0124232 belonged to the LOB protein and EU0128499 belonged to the B3 domaincontaining AP2/ERF family transcription factor.

Metabolome Profile of H8 and H6 during AR Development
Based on the LC-QTOF platform, qualitative and quantitative analyses of the metabolome were conducted on 24 samples, resulting in the detection of 13,088 peaks, with the annotation of 2396 metabolites.PCA results indicated that the metabolites in these samples exhibited obvious differences, consistent with the AR phenotype and hormone levels (Figure 5A).The PC1 can explain 35.69% of the features of the original data set.According to PC1, the separation of the same E. ulmoides variety at different times (S1 stage and other two stages) was uncovered, indicating that the metabolites in the S2 and S3 stages were quite different from those in the S1 stage.Conspicuously, the treatments of different E. ulmoides varieties were separated on the PC2, and PC2 explained 21.8% of the characteristics of the original data set.Furthermore, OPLS-DA results indicated Q 2 > 0.9, suggesting a highly reliable model (Figure 5B).
The differences in accumulation patterns of metabolites in different varieties and stages were manifested (Figure 5C).All metabolites were divided into six clusters.The content of metabolites in cluster 1 was the highest in H8S2 and H8S3, and the content in other groups was low.The content of metabolites in cluster 2 was the highest in H6S2 and H6S3 but lower in other groups.The content of metabolites in cluster 3 was the highest in H8S1 and H6S1 groups but lower in other groups.Different biological repeats were also clustered into a cluster, which indicated the good homogeneity of biological repeats and the high reliability of metabolic data.Above all, cluster results and PCA results indicated significant differences of metabolites in different E. ulmoides varieties and AR development stages.

Metabolome Profile of H8 and H6 during AR Development
Based on the LC-QTOF platform, qualitative and quantitative analyses of the metabolome were conducted on 24 samples, resulting in the detection of 13,088 peaks, with the annotation of 2396 metabolites.PCA results indicated that the metabolites in these samples exhibited obvious differences, consistent with the AR phenotype and hormone levels (Figure 5A).The PC1 can explain 35.69% of the features of the original data set.According to PC1, the separation of the same E. ulmoides variety at different times (S1 stage and other two stages) was uncovered, indicating that the metabolites in the S2 and S3 stages were quite different from those in the S1 stage.Conspicuously, the treatments of different E. ulmoides varieties were separated on the PC2, and PC2 explained 21.8% of the characteristics of the original data set.Furthermore, OPLS-DA results indicated Q 2 > 0.9, suggesting a highly reliable model (Figure 5B).The differences in accumulation patterns of metabolites in different varieties and stages were manifested (Figure 5C).All metabolites were divided into six clusters.The content of metabolites in cluster 1 was the highest in H8S2 and H8S3, and the content in other groups was low.The content of metabolites in cluster 2 was the highest in H6S2 and H6S3 but lower in other groups.The content of metabolites in cluster 3 was the highest in H8S1 and H6S1 groups but lower in other groups.Different biological repeats were also clustered into a cluster, which indicated the good homogeneity of biological repeats and the high reliability of metabolic data.Above all, cluster results and PCA results indicated significant differences of metabolites in different E. ulmoides varieties and AR development stages.
A large number of DEMs were revealed in different comparison groups.In the comparison group of H6S1 vs. H8S1, 1361 DEMs were identified, of which 686 were up-regulated and 675 were down-regulated (Figure 6A).These DEMs were mainly enriched in phenylpropanoid biosynthesis, porphyrin and chlorophyll metabolism, histidine A large number of DEMs were revealed in different comparison groups.In the comparison group of H6S1 vs. H8S1, 1361 DEMs were identified, of which 686 were up-regulated and 675 were down-regulated (Figure 6A).These DEMs were mainly enriched in phenylpropanoid biosynthesis, porphyrin and chlorophyll metabolism, histidine metabolism, nicotinate and nicotinamide metabolism, and alpha-linolenic acid metabolism pathways (Figure 6B).In the comparison group of H6S2 vs. H8S2, 1395 DEMs were exposed, of which 536 were up-regulated and 859 were down-regulated (Figure 6C).These DEMs were mainly enriched in ABC transporters, phenylpropanoid biosynthesis, indole alkaloid biosynthesis, purine metabolism, and nicotinate and nicotinamide metabolism (Figure 6D).In the comparison group of H6S3 vs. H8S3, 1402 DEMs were unveiled, of which 730 were up-regulated and 672 were down-regulated (Figure 6E).These DEMs were mainly enriched in phenylpropanoid biosynthesis, carbon metabolism, flavonoid biosynthesis, ubiquinone and other terpenoid-quinone biosynthesis, and isoflavonoid biosynthesis (Figure 6F).Phenylpropanoid biosynthesis was enriched in every stage of AR development.

Transcriptome and Metabolome of Key Pathways duing AR Development
The key metabolic pathways in the AR development of E. ulmoides were constructed, including phenylpropanoid biosynthesis, flavonoid biosynthesis, and isoflavonoid biosynthesis.Functional genes and key metabolites involved in three pathways were identified.The difference profiles of gene expression and metabolite content between two E. ulmoides cultivars were exhibited in H8S1 vs. H6S1, H8S2 vs. H6S2, and H8S3 vs. H6S3 (Figure 7).
In the phenylpropanoid biosynthesis pathway, compared with H8, the metabolite content of H6 showed significant differences in three stages of AR development.The metabolite content of p-coumaryl alcohol and 4-hydroxycinnamyl alcohol 4-D-glucoside of H6 was significantly up-regulated in all three stages.The metabolite content of the coumarine, mlethyilchavicol, and syringin of H6 was significantly up-regulated at the S2 stage and the S3 stage.In addition, the content of some metabolites of H6 was also up-regulated or down-regulated in different stages, including phenylalanine, cinnamic acid, and p-coumaralde hyde (Figure 7A).The DEG profiles uncovered that six phenylalanine ammonialyase (PAL), six 4-coumarate: CoA ligase (4CL), two cytochrome P450 73A (CYP73A), 21 cinnamoyl-CoA reductase (CCR), 30 cinnamyl-alcohol dehydrogenase (CAD), and one coniferyl-alcohol glucosyltransferase (UGT72E) genes were up-regulated or down-regulated in three stages of AR development.In the S2 and S3 stages, more than half of the CCR genes were significantly up-regulated.The two CCR genes (EU0105592 and EU0106052) were significantly up-regulated in three stages.Most CAD genes were significantly up-regulated in the S2 stage and down-regulated in the S1 and S3 stages.As the up-regulated DEGs, EU0113086 (CYP73A), EU0122718 (CAD), and EU0124962 (CAD) had the largest change folds, which were 6.3, 5.8, and 7.0 respectively (Figure 7B).
In the flavonoid biosynthesis pathway, compared with H8, the metabolite content of H6 also had significant differences.The metabolite content of sakurane tin, luteoforol, liquiritigenin, and 7,4-dihydroxyflavone of H6 was significantly up-regulated in all three stages.In addition, the metabolite content of H6, such as naringin, neoheperidin, and delohiridin, also suggested significant up-regulation (Figure 7A).A total of 20 bifunctional dihydroflavonol 4-reductase (DFR) genes were differentially expressed during AR development.The most DFR genes were up-regulated at S2 stage, and the most DFR genes were down-regulated at S1 and S3 stages.EU0100163 (DFR) demonstrated the largest change fold in the S2 stage (Figure 7B).
In the isoflavonoid biosynthesis pathway, significant differences in metabolite content and gene expression between H6 and H8 were also observed.The content of four metabolites of H6 were significantly up-regulated in all three stages, including malonylgenistin, isoformononetin, malonyldaidzin, and (-) Medic arpin-3-O-glucoside-6 ′ -O-malonate (Figure 7A).The seven 2-hydroxyisoflavanone dehydratase (HIDH) genes were differentially expressed during AR development.The vast majority of HIDH genes conveyed the same expression pattern in the S1 and S2 stages, i.e., down-regulated expression in the S1 stage and up-regulated expression in the S2 stage.EU0132738 (HIDH) and EU0120511 (HIDH) showed the largest change fold in the S2 stage.Additionally, cytochrome P450 81E (CYP8IE) genes were also differentially expressed in different stages of AR development (Figure 7B).
The 130 lncRNAs with potential regulatory capability were screened, and cis-targeting hub genes in the modules significantly were correlated with JA and DEGs in phenylpropanoid biosynthesis, flavonoid biosynthesis, and isoflavonoid biosynthesis pathways (Figure 8).Altogether, 17 lncRNAs exhibited differential expressions, and 113 lncRNAs possessed potential regulatory effects during the AR development of E. ulmoides (Figure S2).Eminently, MSTRG.50256.1 was the cis-acting lncRNA of the cytochrome P450 (CYP450) superfamily gene (EU0127354) in the module, which significantly correlated with JA, and it performed simultaneously as the cis-acting lncRNA of the CCR gene (EU0127349) in the phenylpropanoid biosynthesis pathway (Figure 8).

Discussion
AR formation is the crucial process of cutting propagation and the primary bottleneck restricting the survival of cuttings for horticultural and forestry plants [49,50].Therefore, understanding the regulation of AR formation is extremely necessary in cultivation and breeding.Complex interactions among genotypes, hormone crosstalk, metabolic states, epigenetic regulation, and other factors are considered to shape the AR characteristics of woody plants [51].In this study, the key molecules and important pathways of E. ulmoide AR formation were comprehensively analyzed by integrating physiological data and whole transcriptomes and metabolomes.What needs to be emphasized is that these results only revealed preliminary associations or patterns between the key molecules and important pathways and AR development.These potential causal relationships provide

Discussion
AR formation is the crucial process of cutting propagation and the primary bottleneck restricting the survival of cuttings for horticultural and forestry plants [49,50].Therefore, understanding the regulation of AR formation is extremely necessary in cultivation and breeding.Complex interactions among genotypes, hormone crosstalk, metabolic states, epigenetic regulation, and other factors are considered to shape the AR characteristics of woody plants [51].In this study, the key molecules and important pathways of E. ulmoide AR formation were comprehensively analyzed by integrating physiological data and whole transcriptomes and metabolomes.What needs to be emphasized is that these results only revealed preliminary associations or patterns between the key molecules and important pathways and AR development.These potential causal relationships provide rich genetic and metabolite resources for future research, and the determination of causal relationships needs to be explored from multiple perspectives in subsequent work.This study significantly promotes the comprehension of the regulation of AR formation in E. ulmoides and other plants.

Differences in Plant Hormone Content between Two E. ulmoides Cultivars
The IAA level of H6 (easy-to-root cultivar) was significantly higher than that of H8 (difficult-to-root cultivar), which indicated that IAA is an important endogenous hormone that regulates AR formation.The level of IAA increased during the critical period of AR formation, significantly promoting the occurrence of AR.Auxin is a key hormone affecting the formation of AR in plants, and it can interact with other hormones to jointly affect the formation of AR [13].The level of IAA can be used as one of the important indexes to judge the ability of cutting the roots of improved varieties of E. ulmoides.
Notably, the JA level of two E. ulmoides cultivars decreased significantly from the 0th day (S1 stage) to the 18th day (S2 stage).The JA level of petunia, carnation, and pea cuttings reached the highest level within a period of time after cutting and then decreased [52][53][54].JA is a hormone involved in abiotic stress response [55][56][57][58].The higher JA level in the 0th day and the lower JA level in the 18th day of H6 variety may be important factors for superior AR formation ability.The dual effects of JA on AR development may be related to level and stage.This view is supported by some studies [10].For example, the lower wound-induced JA level of transgenic plants reduced the numbers of AR in petunia [15].At the same time, the SA level also changed reversely in the two varieties from the 0th day to the 18th day with a different order than the JA level, which implies that there may be distinct regulatory mechanisms and interfering signal transduction between JA and SA during the early stages of AR formation [59][60][61].

Several Hub Genes Regulate AR Development
JA can interact with various hormones to regulate AR formation [10,62].In this study, several hub genes that can directly or indirectly regulate hormone levels were found in the gene module and were significantly correlated with JA.The gene modules with positive correlation and negative correlation of JA content were manifested by WGCNA.Twenty hub genes were screened from the two modules, and most of them were differentially expressed in different stages of AR development in E. ulmoides.EU0120751 and EU0125382 were hub genes in modules with significant positive correlation with the JA level.EU0120751 was an F-box protein.In A. thaliana, TIR1 and auxin-signaling F-box 2 protein could modulate JA homeostasis and control AR initiation by specific sensing complexes with IAA6, IAA9, and/or IAA17 [63].EU0125382 was a SAUR-like auxin-responsive family protein.The WOX11-SAUR-auxin signaling regulatory module was essential for AR development in poplar [64].EU0120751 and EU0125382 may directly or indirectly affect the accumulation of JA and control the formation of AR of E. ulmoides.EU0124232 and EU0128499 were hub genes in modules with a significantly negative correlation with JA level.EU0124232 belonged to the LOB protein.In apple, the MdLBDs promoter contained many cis-acting regulatory elements related to the response of growth regulators, and the over-expression of MdLBD16a in tomato increased the number of AR [65].Our study revealed that AT5G06080 was the orthologous gene of EU0124232.It has been demonstrated that AtLBD33 (AT5G06080) played an important regulatory role in A. thaliana root development [66,67].EU0128499 belonged to the B3 domain-containing AP2/ERF family transcription factor.AP2/ERF protein played a vital role in signal pathways mediated by different plant hormones [68].OsAP2/ERF-40 could promote AR development in rice [69].EU0120751, EU0125382, EU0124232, and EU0128499 may modulate the AR development of E. ulmoides by regulating the plant hormone network composed of JA, IAA, ethylene, and other hormones, which needs to be further studied by molecular biology methods, such as CRISPR, in the future [70].Previous studies initially defined the AR development phase of E. ulmoides [8], and the functional mechanism of each key gene at every developmental stage can be scrutinized in prospective research endeavors.

Metabolic Pathway Enrichment of DEGs and DEMs during AR Development
The integrated analysis of transcriptomics and metabolomics revealed that DEGs and DEMs were co-enriched in phenylpropanoid biosynthesis, flavonoid biosynthesis, and isoflavonoid biosynthesis, suggesting that these three metabolic pathways were essential to form AR for E. ulmoides.The phenylpropanoid biosynthesis pathway was also enriched during AR formation in other plants [42,71,72].Many genes were differentially expressed, including six PAL, six 4CL, two CYP73A, 21 CCR, 30 CAD, one UGT72E, 20 DFR, seven HIDH, and seven CYP8IE.The PAL genes have been demonstrated to affect the formation of AR in walnut, wheat, and Echinacea purpurea [73][74][75].Conspicuously, the genes of E. ulmoides H6 were generally significantly up-regulated.Therefore, the difference in gene expression level on the 18th day may be an important reason for the difference in AR development between the two E. ulmoides varieties.Some genes with the same expression pattern were observed, such as CAD genes (EU0103834, EU0103838, and EU0122718).At the same time, genes with opposite expression patterns also were found, such as CCR genes (EU0121857 and EU0123001).This suggests that these genes may have manifold effects in AR formation.

LncRNAs Regulate Plant Hormone Signal Transduction during AR Development
LncRNAs can regulate mRNAs by co-localization [76][77][78], and they exhibit a critical regulatory role for AR formation in forest trees and other plants [79][80][81].For example, the lncWOX11a suppressed AR development by regulating the expression of PeWOX11a in poplar [82].The lncWOX5 negatively regulated WOX5, and lncWOX11 positively regulated WOX11 during the AR formation of poplar [46].The overexpression of lncWOX5 negatively regulated the development of AR in Populus [83].A total of 130 cis-acting lncRNAs were identified associated with the hub genes of modules significantly correlated with JA and DEGs in the phenylpropanoid biosynthesis, flavonoid biosynthesis, and isoflavonoid biosynthesis pathway, with 17 of them exhibiting significant differential expression during root cuttings in favorable conditions of other growth environment factors in E. ulmoides.
MSTRG.50256.1, as an exemplar of cis-acting lncRNA, potentially regulates the CYP450 superfamily gene (EU0127354) in modules significantly correlated with JA and the CCR family gene (EU0127349) in the phenylpropanoid biosynthesis pathway.CYP450 is a huge enzyme superfamily with heme as a cofactor, which is known as a universal biocatalyst and performs an important feature in the metabolic process associated with plant growth and development, such as fatty acids, phytosterols, plant hormones, and phenylpropanoids [84][85][86][87].In the phenylpropanoid biosynthesis pathway, CYP450 monooxygenase 98 catalyzed a rate-limiting step [88].In soybean, CYP728H1 was proven to be co-expressed with the genes related to the phenylpropanoid metabolism [89].In Linum usitatissimum L., seven gene families, namely CYP73, 74,75,76,77,84, and 709, encoded enzymes associated with the phenylpropanoid metabolism [90].AR can be formed in many situations, including abiotic stresses and exogenous hormones, in addition to cuttings [91,92].When cuttings are accompanied by floods and exogenous hormones, the functions of these lncRNAs may be stimulated.MSTRG.50256.1, as the cis-acting lncRNA of the CYP450 gene in the module significantly correlated with JA and the CCR gene in the phenylpropanoid biosynthesis pathway, emphasizes a global regulatory function that needs to be revealed urgently in the future during the AR development of E. ulmoides.

Plant Materials and Sampling
The semi-lignified cuttings were taken from three-year-old E. ulmoides H6 and E. ulmoides H8, cultivated from the Eucommia Engineering Research Center of State Forestry and Grassland experimental base (34 • 500 ′ N, 112 • 330 ′ E, altitude: 108 m) located in the Mengzhou, Henan province, China (Figure 9).The semi-lignified spring shoots, measuring 10-15 cm in length, were chosen as cuttings, and 2-3 mature leaves on the top tender shoots were retained.The cuttings were broken manually from the base without scissors.These cuttings were dipped in the special rooting agent (patent number: CN202011236585.9)for 10 s, and they were inserted into the cutting bed that was disinfected with a 500-fold dilution of 40% carbendazim wettable powder before cutting.The cutting bed was 8 m long and 1.5 m wide, laid with 0.15 m of coarse river sand at the bottom, and covered with 0.15 m of matrix (turfy soil:perlite = 3:1) at the top.The ambient temperature was controlled at 20~29 • C, and the relative humidity was controlled above 90%.E. ulmoides cuttings were studied after 0, 18, and 32 days, which were recorded as S1, S2, and S3 respectively.These sampling times were considered by referring to our previous study on the primary definition of the AR development phase of E. ulmoides [8].Additionally, for convenience, S1-S3 stages of H6 and H8 cuttings were named H6S1-H6S3 and H8S1-H8S3.Stem segments with lengths of 1-3 cm at the base of cuttings were collected, and two biological replicates were made for whole transcriptome strand-specific RNA sequencing and four biological replicates for LC-MS/MS non-targeted determination, respectively.All samples were immediately frozen in liquid nitrogen and stored at −80

Determination of the Content of Plant Endogenous Hormones
High performance liquid chromatography-tandem mass spectrometry (HPLC-MS/MS) was used to determine the content of endogenous hormones, including IAA, ABA, JA, and SA.The determination conditions were as follows: BEH C18 column (2.1 × 100 mm, 1.8 µm), 0.1% formic acid aqueous solution as mobile phase A, acetonitrile as mobile phase B, linear gradient elution, flow rate 0.25 mL/min, injection volume 5 µL, and column temperature 40 °C.The mass spectrometry parameters were as follows.Electrospray ionization source (ESI) was in positive and negative ion ionization mode.The ion source temperature was 500 °C.The ion source voltage was 5500 V-4500 V.The curtain gas was 30 psi.The atomization gas and auxiliary gas were both 50 psi.Multiple Reaction Monitoring (MRM) was adopted for scanning.Comparisons were analyzed using the oneway analysis of variance (ANOVA) test.

RNA Extraction, Library Construction and Sequencing
Total RNA was extracted using TRIzol (Invitrogen, Waltham, MA, USA) following the manufacturer's protocol.RNA quality was determined by 1% agarose gel electrophoresis.The concentration and purity were measured using Nanodrop 2000 (Thermo Fisher

Determination of the Content of Plant Endogenous Hormones
High performance liquid chromatography-tandem mass spectrometry (HPLC-MS/MS) was used to determine the content of endogenous hormones, including IAA, ABA, JA, and SA.The determination conditions were as follows: BEH C18 column (2.1 × 100 mm, 1.8 µm), 0.1% formic acid aqueous solution as mobile phase A, acetonitrile as mobile phase B, linear gradient elution, flow rate 0.25 mL/min, injection volume 5 µL, and column temperature 40 • C. The mass spectrometry parameters were as follows.Electrospray ionization source (ESI) was in positive and negative ion ionization mode.The ion source temperature was 500 • C. The ion source voltage was 5500-4500 V.The curtain gas was 30 psi.The atomization gas and auxiliary gas were both 50 psi.Multiple Reaction Monitoring (MRM) was adopted for scanning.Comparisons were analyzed using the one-way analysis of variance (ANOVA) test.

RNA Extraction, Library Construction and Sequencing
Total RNA was extracted using TRIzol (Invitrogen, Waltham, MA, USA) following the manufacturer's protocol.RNA quality was determined by 1% agarose gel electrophoresis.The concentration and purity were measured using Nanodrop 2000 (Thermo Fisher Scientific, Bothell, WA, USA).RNA integrity was checked by Agient 2100 (Agilent Technologies, Santa Clara, CA, USA).rRNA was removed using the Epicentere Ribo-zero™ rRNA Removal Kit (Epicentere, Middletown, DE, USA).Sequencing libraries were prepared by using the NEBNextR Ultra™ Directional RNA Library Prep Kit for Illumina ® (NEB, San Diego, CA, USA) following the manufacturer's instructions.Finally, 12 libraries were constructed and sequencing was performed on an Illumina NovaSeq 6000 platform, and 150-bp paired-end reads were obtained.

Transcriptome Assembly, lncRNA Identification, and Differential Expression Analysis
Clean reads were obtained by removing reads containing adapter, poly-N, and lowquality reads from raw data.At the same time, the Q20, Q30, GC-content, and sequence duplication level of the clean data were calculated.All downstream analyses were based on clean data with high quality.The clean reads were mapped to the E. ulmoides genome [48] using HISAT2 v2.0.4 [93].Then, the transcriptome was assembled by StringTie v1.3.1 [94] based on the reads mapped to the reference genome.The gffcompare program was used to annotate the assembled transcripts.
The unknown transcripts were used to screen for putative lncRNAs.Transcripts with a length of ≥200 bp, number of exons ≥ 2, and expression levels of fragments per kilobase of exons per million fragments mapped (FPKM) ≥ 0.1 were selected as lncRNA candidates and further screened using four computational approaches (CPC/CNCI/Pfam/CPAT) that could distinguish protein-coding genes from non-coding genes.Using Perl script, the neighboring genes within 100 kb upstream and downstream of lncRNAs were used as cis-target genes of lncRNAs.The network diagram of the relationship between cis-acting lncRNAs and genes was shown by ChiPlot (https://www.chiplot.online/(accessed on 20 October 2023)).

Metabolite Profiling and Data Analysis
Lyophilized powder (50 mg) was weighed and extracted using 1 mL of extract solution (methanol: acetonitrile: water = 2:2:1) and was vortexed for 30 s.It was homogenized in ball mill for 10 min at 45 Hz, then put in an ultrasonic 10 min ice-water bath, and incubated at −20 • C for 1 h.The extracts were centrifuged at 12,000 rpm for 15 min at 4 • C, and supernatant (500 µL) was transferred into EP tubes, then dried in a vacuum concentrator, and 160 µL of extraction solvent (acetonitrile: water= 1:1) was added for reconstitution.It was vortexed for 30 s and sonicated for 10 min (ice-water bath).Finally, the supernatant was obtained by centrifugation at 4 • C at 12,000 rpm for 15 min and transferred to a 2 mL glass vial through a 0.22 µm membrane.The filtrate was analyzed using a UPLC-MS system (UPLC: Acquity I-Class PLUS, Waters, Milford, MA, USA; MS: Xevo G2-XS QTOF, Waters, Milford, MA, USA).The ESI source conditions were set as follows: capillary voltage: 2500 V (positive ion mode) or −2000 V (negative ion mode); cone voltage: 30 V; ion source temperature: 150 • C; desolvent gas temperature 500 • C; backflush gas flow rate: 50 L/h; desolventizing gas flow rate: 800 L/h.
Metabolites were characterized based on secondary spectral information using the METLIN and MVDBV database.They were quantified by triple quadruple mass spectrometry in the multiple reaction monitoring mode.PCA analysis, OPLS-DA analysis, and expression profiling of all metabolites of 24 metabolome samples were performed using the Metware Clouds (https://cloud.metware.cn(accessed on 10 July 2023)).The differently expressed metabolites (DEMs) were obtained from 3 combinations of H6S1 vs. H8S1, H6S2 vs. H8S2, and H6S3 vs. H8S3 groups.The criteria were the threshold of variable importance in projection ≥ 1, |log 2 fold change| ≥ 1, and p-value ≤ 0.05.The KEGG enrichment analyses of DEMs were carried out using clusterProfiler v3.14.3 R package.

Conclusions
In this study, the pivotal factors and main pathways related to the AR formation of E. ulmoides were identified by comparing the physiological indexes, transcriptomes, and metabolomes of two E. ulmoides varieties with different AR formation abilities.Significant differences of IAA, ABA, JA, and SA levels between two varieties were found.The hub genes significantly related to JA were demonstrated, such as F-box, SAUR-like, LOB, and AP2/ERF.DEGs and DEMs were enriched in multiple pathways, including phenylpropanoid biosynthesis, flavonoid biosynthesis, and isoflavonoid biosynthesis.Most PAL, CCR, CAD, DFR, and HIDH genes were significantly up-regulated on the 18th day in easy-to-root cultivar.The 130 cis-targeting lncRNAs with potential regulatory roles in AR formation were revealed.The role of these key molecules and pathways in AR development can be further researched through CRISPR, single-cell, and spatial omics technologies in future studies.In conclusion, these results provide an important theoretical foundation for understanding the molecular mechanism of AR formation of E. ulmoides.

Figure 1 .
Figure 1.The morphological characteristics and hormone levels during the AR formation of two E. ulmoides varieties.(A) Morphological characteristics.(B) IAA level.(C) ABA level.(D) JA level.(E) SA level.S1, cutting on the 0th day.S2, cutting on the 18th day.S3, cutting on the 32nd day.Different letters denote significant differences at p < 0.05 level.

Figure 1 .
Figure 1.The morphological characteristics and hormone levels during the AR formation of two E. ulmoides varieties.(A) Morphological characteristics.(B) IAA level.(C) ABA level.(D) JA level.(E)SA level.S1, cutting on the 0th day.S2, cutting on the 18th day.S3, cutting on the 32nd day.Different letters denote significant differences at p < 0.05 level.

Figure 2 .
Figure 2. DEGs and KEGG enrichment analysis during AR development.(A) PCA of samples in three developmental stages.(B) Correlation analysis of samples.The scale bar represents the size of correlation.(C) DEGs of H8 and H6 in three developmental stages.(D) KEGG enrichment of DEGs in H8S1 vs. H6S1.(E) KEGG enrichment of DEGs in H8S2 vs. H6S2.(F) KEGG enrichment of DEGs in H8S3 vs. H6S3.The q-value is p-value after multiple hypothesis testing corrections based on false discovery rate measures.

Figure 2 .
Figure 2. DEGs and KEGG enrichment analysis during AR development.(A) PCA of samples in three developmental stages.(B) Correlation analysis of samples.The scale bar represents the size of correlation.(C) DEGs of H8 and H6 in three developmental stages.(D) KEGG enrichment of DEGs in H8S1 vs. H6S1.(E) KEGG enrichment of DEGs in H8S2 vs. H6S2.(F) KEGG enrichment of DEGs in H8S3 vs. H6S3.The q-value is p-value after multiple hypothesis testing corrections based on false discovery rate measures.

Figure 3 .
Figure 3. WGCNA of expressed genes with FPKM > 1 during E. ulmoides H6 AR development.(A) Hierarchical clustering tree depicting four co-expressed gene modules.(B) Heatmap of the correlation between modules and hormones.(C) Gene significance and module membership in MEblue-JA.(D) Gene significance and module membership in MEturquoise-JA.(E) Correlation network of 20 hub genes in the blue module.(F) Correlation network of 20 hub genes in the turquoise module.

Figure 3 .
Figure 3. WGCNA of expressed genes with FPKM > 1 during E. ulmoides H6 AR development.(A) Hierarchical clustering tree depicting four co-expressed gene modules.(B) Heatmap of the correlation between modules and hormones.(C) Gene significance and module membership in MEblue-JA.(D) Gene significance and module membership in MEturquoise-JA.(E) Correlation network of 20 hub genes in the blue module.(F) Correlation network of 20 hub genes in the turquoise module.

Figure 4 .
Figure 4.The differential expression of hub genes in the blue and turquoise modules significantly correlated with JA. (A) Expression levels of hub genes in the blue module.(B) Differential changes of hub genes in the blue module.(C) Expression levels of hub genes in the turquoise module.(D) Differential changes of hub genes in the turquoise module.The red circles, blue circles, and gray circles represent significant up-regulation, significant down-regulation, and no significant change for hub genes under corresponding comparison groups, respectively.The color depth of red circles and blue circles represents the size of log2 fold change.

Figure 4 .
Figure 4.The differential expression of hub genes in the blue and turquoise modules significantly correlated with JA. (A) Expression levels of hub genes in the blue module.(B) Differential changes of hub genes in the blue module.(C) Expression levels of hub genes in the turquoise module.(D) Differential changes of hub genes in the turquoise module.The red circles, blue circles, and gray circles represent significant up-regulation, significant down-regulation, and no significant change for hub genes under corresponding comparison groups, respectively.The color depth of red circles and blue circles represents the size of log 2 fold change.

Figure 7 .
Figure 7.The differential expression profile of transcriptomics and metabolics in phenylpropanoid biosynthesis, flavonoid biosynthesis, and isoflavonoid biosynthesis in three stages during E. ulmoides AR development.(A) Metabolic pathway map and differential expression profile of metabolite contents.The purple boxes, green boxes, and gray boxes represent significant up-regulation, significant down-regulation, and no significant change for metabolite contents in H6S1 vs. H8S1, H6S2 vs. H8S2, and H6S3 vs. H8S3, respectively.(B) The differential expression profile of genes in the pathways.The red boxes, blue boxes, and gray boxes represent significant up-regulation, significant down-regulation, and no significant change for gene expression in H6S1 vs. H8S1, H6S2 vs. H8S2, and H6S3 vs. H8S3, respectively.The color depth of boxes represents the size of log2 fold change.

Figure 7 .
Figure 7.The differential expression profile of transcriptomics and metabolics in phenylpropanoid biosynthesis, flavonoid biosynthesis, and isoflavonoid biosynthesis in three stages during E. ulmoides AR development.(A) Metabolic pathway map and differential expression profile of metabolite contents.The purple boxes, green boxes, and gray boxes represent significant up-regulation, significant down-regulation, and no significant change for metabolite contents in H6S1 vs. H8S1, H6S2 vs. H8S2, and H6S3 vs. H8S3, respectively.(B) The differential expression profile of genes in the pathways.The red boxes, blue boxes, and gray boxes represent significant up-regulation, significant down-regulation, and no significant change for gene expression in H6S1 vs. H8S1, H6S2 vs. H8S2, and H6S3 vs. H8S3, respectively.The color depth of boxes represents the size of log 2 fold change.

Figure 8 .
Figure 8.The 130 lncRNA cis-target hub genes in modules that significantly correlated with JA and DEGs in phenylpropanoid biosynthesis, flavonoid biosynthesis, and isoflavonoid biosynthesis pathway.Each type of factor is marked with various colors.The lncRNA MSTRG.50256.1 is marked with a pentagram.

Figure 8 .
Figure 8.The 130 lncRNA cis-target hub genes in modules that significantly correlated with JA and DEGs in phenylpropanoid biosynthesis, flavonoid biosynthesis, and isoflavonoid biosynthesis pathway.Each type of factor is marked with various colors.The MSTRG.50256.1 is marked with a pentagram.