Transcriptome Analysis and Genes Function Verification of Root Development of Paeonia suffruticosa under Sandy Loam Cultivation

1  Introduction

Tree peony, which is a deciduous subshrub in the family Paeoniaceae, produces large and brightly colored flowers. In China, it is a popular traditional flower that has high ornamental and medicinal values [1,2]. Traditionally, peony is mainly propagated by cuttings and grafting, both of which can take root effectively, and the best performance is in sandy loam soil [3,4]. The advantages of plant tissue culture include rapid reproduction, high multiplication coefficient, short cycle, and retention of the excellent traits of the mother plant. Accordingly, it has recently been widely used for the rapid reproduction of tree peony [5,6]. However, there are several problems associated with tree peony test-tube seedlings (e.g., poor induction of roots, low rooting rate, low quality of the roots, and low survival rate after transplantation), which seriously affect the development of rapid propagation methods for the commercial cultivation of peony [6,7]. Therefore, to establish an efficient tissue culture system for tree peony, the problems related to poor rooting, especially root formation, must be addressed [4,7]. The explants used in tissue culture were equivalent to cuttings and scions. Therefore, explored the rooting mechanism of seedlings, and combined with the evidence of cuttings and grafting, could provide theoretical support for the rooting of tree peony in vitro. Several studies have been conducted on tree peony root development, but they mainly focused on the effects of the growth medium, hormone levels, and endogenous phenolics [5,7,8]. Unfortunately, the underlying molecular mechanism is still unclear.

Root development is a complex process regulated by nutrients, oxidase activity, environmental conditions, and various endogenous hormones [9–12]. Of these factors, endogenous hormones, especially auxin, are particularly important for plant rooting [13,14]. Auxin has been confirmed as the main hormone that promotes root formation in Vigna radiata, Arabidopsis, Mangifera and other plants [15–18]. More specifically, auxin directly regulates plant root growth and development, but it also indirectly affects the expression of root development-related genes, thereby regulating root formation [19,20].

High-throughput transcriptome sequencing (RNA-seq) can analyze gene expression in specific plant parts and during particular stages [21,22]. It has been commonly used to identify differentially expressed genes (DEGs) and key genes that regulate plant growth and development [4,23]. Many new rooting-related genes were identified in Saccharum sp., Chrysanthemum, and Zea mays by transcriptome sequencing, and their relationships with auxin and other related hormones have been clarified in follow-up studies [21,24,25]. Thus, RNA-seq may be useful for identifying rooting-related genes to elucidate the mechanism mediating root formation in tree peony. The transcriptome sequencing of tree peony to date has primarily focused on tissue culture browning, genetic breeding, and the floral fragrance mechanism [18,26,27]. There are no reports describing research on the mechanism underlying root formation and the relationship with endogenous hormones at the transcription level.

In this study, a transcriptome sequencing analysis was performed to screen and identify DEGs and the associated metabolic pathways at different stages of tree peony root development under sandy loam cultivation. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses revealed many of the DEGs were involved in auxin metabolism. Additionally, the rooting-related genes were screened, after which their expression levels in seedlings, cuttings and grafted seedlings were analyzed by quantitative real-time polymerase chain reaction (qRT-PCR) and their functions were verified. The results of these analyses may be useful for characterizing the main molecular mechanism controlling tree peony root differentiation. Furthermore, the data presented in this manuscript may provide the theoretical basis for future investigations of the tree peony rooting mechanism and the establishment of viable peony tissue culture and rapid propagation systems.

2  Materials and Methods

2.1 Plant Materials and Treatments

Tree peony (Paeonia suffruticosa ‘Fengdanbai’) was used as the study material. All plants were grown in a nursery at the Henan Agricultural University in Zhengzhou, China. Seedlings were approximately 10 years old, cultivated in sandy loam soil, and their growth condition was good. The external morphology and the anatomical structure of the roots were examined to determine the three critical rooting period (Fig. 1): ungerminated root primordium (R1), induced root primordium (R2), and lateral root formation (R3). The roots which were positive and negative 0.5 cm with the occurrence point of lateral roots were intercepted as the materials for RNA-Seq and qRT-PCR. The screening method of the critical period of rooting and the sampling method of cuttings and grafted seedlings referred to our previous research [28]. The current year semi-lignified twigs of ‘Fengdanbai’ which was 10 years old used above were collected as the cuttings and scions for qRT-PCR materials. They were planted in sandy loam soil, and cultivated in a glass greenhouse at a temperature of 15°C ± 2°C and an air humidity of 50%. The phloem of the rootstock interface and roots during the rooting critical stage were collected as qRT-PCR materials [28]. For each sample, the root materials were examined in triplicate, and placed in a sterile 50 mL centrifuge tube, immediately frozen in liquid nitrogen, and stored at −80°C until analyzed.

Figure 1: Morphological and anatomical observation of tree peony root (A/a, Root primordium ungerminated period/R1; B/b, Root primordium induction period/R2; C/c, Lateral root formation period/R3, Rp. Lateral root primordium; Lr. Lateral root)

2.2 RNA Extraction, Construction of the cDNA Library, and Transcriptome Sequencing

Total RNA was isolated from fresh root tissue of R1, R2, and R3 according to a modified CTAB procedure [29]. The RNA concentration and quality were determined using the 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) and the NanoDrop 1000 UV/Visible spectrophotometer (Thermo Scientific, Waltham, MA, USA), respectively. The cDNA library was constructed from the high-quality RNA samples and then sequenced by Shenzhen Hengchuang Biotechnology Co., Ltd. (China). The preliminary quantitative analysis was performed using the Qubit® 2.0 Fluorometer (Life Technologies, Grand Island, NY, USA). The insert fragment size of the library was determined using the 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). The double-stranded cDNA of the above-mentioned samples was sequenced using the GA IIx system (Illumina, CA, USA).

2.3 De Novo Transcriptome Assembly and Annotation

The raw transcriptome sequencing data were filtered to obtain the original clean data. A de novo sequence assembly was used because of a lack of an available fully sequenced tree peony reference genome. This sequence was compared with the sequences of other model plants available in the NCBI database using the BLASTn algorithm. The quality of the library was assessed (e.g., insert length check and matched fragment randomness check). The fragments per kilobase of transcript per million fragments mapped value [30] was used as an index for measuring gene expression levels. The EBseq software [31] was used to analyze the data for the different sample groups. The following criteria were used to identify DEGs (i.e., differentially expressed genes): expression level fold-change ≥2 and false discovery rate <0.01. On the basis of the COG (i.e., cluster of orthologous groups of proteins), GO (i.e., gene ontology), KEGG (i.e., kyoto encyclopedia of genes and genomes), Swiss-Prot, NR (i.e., non-redundant amino acids), and NT (i.e., non-redundant nucleotides) databases, the functions, types, and enriched metabolic pathways of the DEGs were determined.

2.4 Transcriptome Data Verification and Differentially Expressed Gene Analysis

12 differentially expressed genes were randomly selected from the transcriptome sequencing data, and detected the expression levels by qRT-PCR technology in three different developmental stages of tree peony. By analyzing the transcriptome sequencing data, two significant differentially expressed rooting-related genes were selected and cloned. BioXM2.7 was used to translate the open reading frame (ORF). The physicochemical properties and conserved domains were estimated using analytical tools at the NCBI (https://www.ncbi.nlm.nih.gov/). The online software SMART (http://smart.embl-heidelberg.de) was used for domain analysis. DNAMAN8 was used for multiple sequence alignment and MEGA7.0 was used to construct phylogenetic trees with the neighbor-joining method. Primer 5.0 was used for designing the qRT-PCR primers, which were synthesized by Guangzhou Saizhe Biotechnology Co., Ltd. (China). A peony β-tubulin gene (GenBank: EF608942.1) was used as the internal reference gene. The related gene IDs and primer sequences are listed in Table 1.

2.5 Construction of the Arabidopsis Overexpression Vector

Primers containing restriction sites were designed (Table 1) for the amplification of the full-length ORF of the target gene. The PCR product was purified (Tiangen DP204, Beijing, China) and inserted into the pGEM-T Easy vector. The recombinant plasmid was then inserted into Escherichia coli cells. A single clone was sequenced. After confirming the sequence was accurate, the plasmid was extracted and then digested by restriction enzymes (Nco I, Spe I, Takara Biomedical Technology Co., Ltd., Beijing, China). The target fragment was inserted into a pCAMBIA1302 vector.

2.6 Functional Verification of Transformed Arabidopsis Plants

The constructed pCAMBIA1302 overexpression vector was inserted into wild-type Arabidopsis plants according to an Agrobacterium tumefaciens-mediated inflorescence dip method [32]. Transgenic seedlings were selected on 1/2MS medium (containing 100 mg⋅L−1 Cef and 25 mg⋅L−1 Hyg). The number of PsIAA27 and PsARF19 overexpression positive transgenic Arabidopsis obtained were 15 and 18, respectively. They were then transplanted into matrix (perlite: ver-miculite = 1:1) and analyzed. The phenotypic examination and identification were initiated starting with the T1 generation plants. Continuing to cultivate transgenic Arabidopsis to harvest T2 homozygous transgenic Arabidopsis, 3 plants with the same growth status were selected for qRT-PCR and index determination. Details regarding the primers are listed in Table 1.

2.7 Statistical Analysis

All data in this research underwent a two-way analysis of variance (ANOVA), and significant differences between the means were tested using Duncan’s post-hoc test (p ≤ 0.05). All statistical analyses were conducted using the Microsoft Office 2010 and SPSS 19.0 (SPSS, Chicago, USA) program for Windows.

3  Results

3.1 Morphological and Anatomical Observation of Adventitious Root Formation Process

Morphological anatomic observation of the growth process of root primordium germination and lateral roots was shown in Fig. 1. The root surface was smooth in the Fig. 1A, and there was no sign of root primordium sprouting (Fig. 1a). The root surface in Fig. 1B was smooth and without convex, but it had obvious root primordium germination phenomenon (Fig. 1b). A large number of bulges had been found on the root surface in Fig. 1C and indefinite root had been formed (Fig. 1c). The above three periods were named as R1 (ungerminated root primordium), R2 (induced root primordium) and R3 (lateral root formation), respectively, and these roots were used as test materials for transcriptome sequencing analysis.

3.2 Illumina Sequencing and de Novo Transcriptome Assembly

A total of 31.63 Gb raw data (i.e., raw sequence data) were obtained during the transcriptome sequencing analysis in the three periods of R1, R2 and R3 (Table 2). Additionally, 207,856 unigenes were assembled after eliminating the adapter sequences, ambiguous reads, and low-quality reads, with a minimum length (i.e., minimum sequence length) of 200 bp, a maximum length (i.e., maximum sequence length) of 16,605 bp, and a mean length (i.e., mean sequence length) of 911.98 bp. The N50 (i.e., contig N50 and scaffold N50) length, N90 (i.e., contig N90 and scaffold N90) length, and GC (i.e., guanine and cytosine) content of the assembled unigenes were 1,522 bp, 358 bp, and 41.05%, respectively.

3.3 Functional Annotation of the Tree Peony Rooting Transcriptome

A total of 120,188 unigenes were annotated using one nucleotide database and five functional protein databases (Table 3). Of these unigenes, 111,369 (92.66%) were annotated on the basis of the NR database, 96,232 (80.07%) were annotated according to the NT database, 72,879 (60.64%) were annotated using the Swiss-Prot database, 45,216 (37.62%) were annotated on the basis of the COG database, 71,752 (59.70%) were annotated according to the KEGG database, and 62,769 (52.23%) were annotated using the GO database.

3.4 Analyses of DEGs

In the R1 vs. R2 comparison, 41,162 DEGs were identified, of which 10,642 were up-regulated and 30,520 were down-regulated. In the R1 vs. R3 comparison, 53,436 DEGs were detected, of which 14,636 were up-regulated and 38,800 were down-regulated. The R2 vs. R3 comparison revealed 34,538 DEGs, of which 14,496 were up-regulated and 20,042 were down-regulated. The comparisons indicated that as the root primordium developed, an increasing number of unigenes had a down-regulated expression level. Thus, the number of down-regulated genes exceeded the number of up-regulated genes (Fig. 2).

Figure 2: Differentially expressed genes in three evelopmental stages of tree peony rooting (R1 vs. R2: R1/R2 stage, R1 vs. R3: R1/R3 stage, and R2 vs. R3: R2/R3 stage. The numbers of DEGs are marked in the column chart)

The differential gene expression in each metabolic pathway varied among the three stages of tree peony root development (Fig. 3). The comparison of differential gene expression in different stages indicated that most of the DEGs in the pathways were down-regulated. Of the examined metabolic pathways, carbohydrate metabolism had the most DEGs. There were 314 DEGs between R1 and R2, of which 218 were down-regulated (i.e., 69%). There were 371 DEGs between R2 and R3, most of which were down-regulated. There were 336 DEGs between R1 and R3, of which 102 were up-regulated and 234 were down-regulated. Moreover, the expression level of the down-regulated genes was double that of the up-regulated. The number of DEGs in hormone was second more than others. There were 226 DEGs between R2 and R3, of which 129 were down-regulated (i.e., 57.1%) and 97 were up-regulated (i.e., 42.9%). There were 321 DEGs between R1 and R2, of which 256 were down-regulated (i.e., 79.8%) and 65 were up-regulated (i.e., 20.2%). There were 227 DEGs between R1 and R3, of which 56 were up-regulated (i.e., 24.7%) and 171 were down-regulated (i.e., 75.3%).

Figure 3: Analysis of the DEGs of main metabolic pathways in the R1/R2 stage, R1/R3 stage, and R2/R3 stage (The numbers of DEGs are marked in the column chart)

Regarding the R1 vs. R2 comparison, the most and least common pathways associated with the DEGs were the auxin and gibberellin pathways, respectively (Fig. 4). An examination of the expression levels of the genes related to the hormone metabolism pathway revealed there were more down-regulated genes than up-regulated genes. In the R2 vs. R3 comparison, there was roughly the same number of up-regulated and down-regulated genes involved in the auxin pathway. For the R1 vs. R3 comparison, there were substantially more down-regulated genes than up-regulated genes. The down-regulated DEGs were mainly related to the zeatin, gibberellin, and abscisic acid metabolic pathways. Of the examined hormone pathways, auxin pathway had the most DEGs. There were 127 DEGs between R1 and R2, of which 9 were up-regulated (i.e., 7.1%) and 118 were down-regulated (i.e., 92.9%). There were 87 DEGs between R2 and R3, of which 42 were up-regulated (i.e., 48.3%) and 45 were down-regulated (i.e., 51.7%). There were 39 DEGs between R1 and R3, of which 9 were up-regulated (i.e., 23.1%) and 30 were down-regulated (i.e., 76.9%).

Figure 4: Analysis of the DEGs of four kinds of hormonal pathways in the R1/R2 stage, R1/R3 stage, and R2/R3 stage (The numbers of DEGs are marked in the column chart)

3.5 GO Term Enrichment Analysis of DEGs

An analysis of the annotated peony differential unigenes using the GO database divided the sequences into 55 subcategories in the three main functional categories (i.e., biological process, cellular component, and molecular function) (Fig. 5). The most common annotated GO terms in the biological process category were regulation of transcription, transcription, carbohydrate metabolic process, metabolic process, and defense response. The predominant cellular component GO terms were integral component of membrane, nucleus, cytoplasm, plasma membrane, and intracellular. The main molecular function GO terms were ATP binding, DNA binding, RNA binding, metal ion binding, and sequence-specific DNA binding transcription factor activity.

Figure 5: Differential gene GO classification in the R1/R2 stage, R1/R3 stage, and R2/R3 stage

3.6 KEGG Pathway Enrichment Analysis of DEGs

The active metabolic pathways associated with the tree peony rootings in the three examined differentiation stages were clarified by a KEGG pathway enrichment analysis. The 10,687, 8,232, and 6,907 DEGs detected for the R1, R2, and R3 periods were involved in 133, 133, and 132 metabolic pathways, respectively (Fig. 6). The top three metabolic pathways were plant hormone signal transduction, starch and sucrose metabolism, and RNA transport. Other enriched metabolic pathways among the DEGs were mRNA surveillance pathway, endocytosis, and phenylpropanoid biosynthesis. The plant hormone signal transduction played an important regulatory role in the three root development stages of tree peony.

Figure 6: KEGG pathway enrichment scatter map of differentially expressed transcripts in the R1/R2 stage, R1/R3 stage, and R2/R3 stage

3.7 Rooting-Related Transcription Factor Analysis

A total of 58 transcription factor families were annotated following the sequence alignment and screening of known transcription factors (Fig. 7). The top five transcription factor families were MYB (393 genes), bHLH (271 genes), bZIP (212 genes), C3H (212 genes), and ERF (190 genes). Additionally, there were many transcription factors belonging to auxin transcription factor ARF (94 genes), and transcription factors AP2 (97 genes), GATA (56 genes), LBD (62 genes) and NAC (172 genes) involved in auxin regulation, the transcription factors ARF (94 genes), bHLH (271 genes), bZIP (212 genes), NAC (172 genes), WOX (17 genes), WRKY (166 genes), etc., were closely related to root genesis.

Figure 7: Transcription factor families (The numbers of transcription factor are marked at the top of the column chart)

3.8 qRT-PCR Verification of Transcriptome Data

In order to verify the reliability of the transcriptome data, 12 differentially expressed genes were randomly selected and detected the expression levels by qRT-PCR in three different developmental stages, and then compared with the RNA-seq results. The specific results are shown in the Fig. 8, there were some differences in the expression levels of the 12 genes between qRT-PCR and RNA-seq, but the overall expression trends were analogous in the three periods. Therefore, the transcriptome data was reliable.

Figure 8: qRT-PCR verification of transcriptome data (Random selection of 12 differential genes, Each column represents the mean of three independent measurements. Error bars represent the SD of the mean values. Diferent letters represent signifcant diference at p ≤ 0.05 using Duncan’s post-hoc test)

3.9 Heatmap Analysis of DEGs in Auxin Pathway

The cluster annotation and heatmap analysis of DEGs in the auxin pathway (Fig. 9) showed that Unigene13430_All was highly expressed in the R2 period, extremely low in the R3 period, and centered in the R1 period. CL10096.Contig1_ All also expressed very low in R3, but high in R1, and relatively low in R2. Unigene 61813_All and CL14666.Contig4_All were highly expressed only in R1 period, CL12740.Contig1_All was highly expressed in R1 and R2. The DEGs were mainly annotated in IAA27, ARF19, IAA9, ARF8 and other families. Among them, Unigene13430_All and CL10096.Contig1_All showed the most significant differences in expression levels at R1, R2 and R3 periods.

Figure 9: Heatmap analysis of the IAAs and ARFs gene expression levels in the R1, R2, and R3 stages (A: IAAs; B: ARFs. The colour scale indicates a relative fold-change, with red indicates high expression and green low expression)

3.10 Cloning and Analysis of DEGs in Auxin Pathway

Since the most DEGs were clustered into auxin pathway, mainly IAA27 and ARF19, the related genes Unigene13430_All and CL10096.Contig1_All were cloned by combining transcriptome database. Their sequences were amplified by PCR. The full-length Unigene13430_All coding sequence consists of 843 bp, encoding 280 amino acids. The full-length CL10096.Contig1_All coding sequence comprises 1,470 bp, encoding 489 amino acids. A smart BLAST protein sequence analysis revealed that Unigene13430_All encode a complete AUX-IAA domain, which was similar with the IAA27 of Macadamia ntegrifolia, Telopea speciosissima, Vitis riparia, Juglans regia, Carya illinoinensis, Ziziphus jujuba, Prunnus avium, Vigna unguiculata and Cajanus cajan (Fig. 10A). CL10096.Contig1_All also encode a complete AUX-IAA domain, which was similar with the ARF19 of Vitis inifera, Prosopis alba, Carica papaya, Pistacia vera, Jatropha curcas, Quercus lobata, Hevea brasiliensis, Camellia sinensis and Ricinus communis (Fig. 10B).

Figure 10: Protein conservative domain analysis (A: Unigene13430_All/PsIAA27, B: CL10096.Contig1_All/PsARF19)

Phylogenetic tree clustering analysis of Unigene13430_All and CL10096.Contig1_All with all IAA and ARF sequences in Arabidopsis thaliana showed that Unigene13430_All was clustered with IAA27 in Arabidopsis thaliana, and CL10096.Contig1_All was clustered with ARF19 in Arabidopsis thaliana. Therefore, Unigene13430_All and CL10096.Contig1_All were designated as PsIAA27 and PsARF19 in tree peony, respectively (Fig. 11).

Figure 11: Evolutionary tree analysis (The black triangle marker gene is the gene used in this research, Unigene13430_All: PsIAA27, CL10096.Contig1_All: PsARF19)

3.11 RT-PCR of PsIAA27 and PsARF19 in Different Tree Peony Roots

The expression of PsIAA27 and PsARF19 were detected in the rooting critical period of seedlings, cuttings and grafted scion by RT-PCR (Fig. 12). As shown in Fig. 13, their expression levels were differences at different rooting stages of seedlings, cuttage seedlings and graft seedlings. Among them, the expression level of PsIAA27 showed a increased first and then decreased trend in seedlings and it was highest in R2 stage. However, the expression level of PsIAA27 was gradual upward in cuttings and grafted scion, and both achieved the maximum value in R3 stage. The PsIAA27 expression abundance in cuttings was higher than that of grafted scion in R3 stage, but lower than that of seedlings in R2 stage. The expression level of PsARF19 showed a upward trend during the three periods, and all got the highest value in R3 stage. The PsARF19 expression abundance in seedlings was the highest in R3 stage, but lower than that of cuttings in R2 stage, and no significant difference.

Figure 12: Rooting critical period of tree peony seedlings, cuttings and grafted scion for qRT-PCR (A: Seedlings, B: Cuttage seedlings, C: Graft seedlings. The black scale bar is 1 cm)

Figure 13: RT-PCR of PsIAA27 and PsARF19 during roooting of tree peony seedlings, cuttings and grafted scion. (A: PsIAA27, B: PsARF19. Diferent letters represent signifcant diference at p ≤ 0.05 using Duncan’s post-hoc test)

3.12 Functional Analysis of PsIAA27 and PsARF19 in Transgenic Arabidopsis

The overexpression vectors pD1301S-PsIAA27 and pD1301S-PsARF19 were constructed and inserted into wild-type Arabidopsis plants. The PsIAA27 and PsARF19 expression levels were significantly higher in the T1 transgenic Arabidopsis plants than in the wild-type control plants (Fig. 14A), indicating that the overexpression vectors were successfully incorporated into the T1 Arabidopsis plants. Continuing to cultivate transgenic Arabidopsis to harvest T2 homozygous transgenic Arabidopsis. Further analyses of the T2 Arabidopsis generation revealed that the plants transformed with the PsIAA27 overexpression vector had significantly longer roots and had less lateral roots than the wild-type control plants (Fig. 14B). Furthermore, the lateral root number and lehgth of the T2 Arabidopsis plants transformed with the PsARF19 overexpression vector was significantly greater than that of the wild-type control plants (Fig. 14B). These observations suggested that PsIAA27 promoted the elongation of roots, but inhibited the development of lateral roots, whereas PsARF19 significantly promoted roots elongation and the development of lateral roots (Fig. 14C).

Figure 14: Positive Arabidopsis RT-PCR verification and phenotypic traits. (A: Positive identification, B: Morphological variation, C: Phenotypic traits. Diferent letters represent signifcant diference at p ≤ 0.05 using Duncan’s post-hoc test)

4  Discussion

The development of tree peony industrialization urgently needs the establishment of a rapid propagation system by tissue culture. However, in the process of tissue culture of tree peony, there are problems such as difficult to induce explants to form roots, and a large proportion of the roots that do form are poor quality, with a disconnect between the root and the stem. This results in low survival rates of plantlets at the trans-planting stage [5,7,33,34]. Therefore, the transcriptome analysis of the critical period of rooting of seedlings, and combine with the relevant verification of cuttings and grafting, can provide a certain reference for the rooting of tree peony in vitro. In this study, a transcriptome sequence database of tree peony seedlings under sandy loam cultivation at different rooting periods was established. A total of 31.63 Gb raw data were obtained for the sequenced transcriptome (Table 2). Additionally, 207,856 unigenes were assembled, with an average length of 911.98 bp (Table 3). Moreover, 45,216 unigenes were annotated according to at least one database (i.e., NR, NT, Swiss-Prot, KEGG, and GO; Table 3).

The identified DEGs at the different tree peony rooting stages were mainly associated with the auxin pathway, but many genes had down-regulated expression levels in R2 (i.e., induced root primordium; Fig. 4). This may explain why there were relatively few tree peony roots. During the transition from R2 to R3 (i.e., from induced root primordium to lateral root formation), there was no major difference in the number of up-regulated and down-regulated genes in the auxin pathway (Fig. 4). This indicates that after the root primordium forms, auxin likely plays a regulatory role and does not function as an inducer. From R1 to R2 (i.e., ungerminated root primordium to induced root primordium), there were considerably fewer DEGs in the gibberellin and zeatin pathways than in the auxin pathway (Fig. 4). During this period, the DEGs in the abscisic acid pathway were mostly down-regulated, implying that abscisic acid and auxin coordinately influence root development. From R2 to R3, there was an equal number of up-regulated genes and down-regulated genes in the auxin pathway (Fig. 4). Additionally, there were more DEGs in the auxin pathway than in the other three pathways, suggesting that the auxin content increases during the entire rooting process, possibly reflecting a major regulatory role. The DEGs in auxin pathway with the most significant changes in expression were mainly annotated in IAA27 and ARF19 (Fig. 9). Both Unigene13430_All and CL10096.Contig1_All were cloned; on the basis of a bioinformatics analysis, these genes were named PsIAA27 and PsARF19 (Figs. 10 and 11), respectively.

Bassa et al. [35] reported that the expression pattern of SI-IAA27 differed from that of traditional auxin genes. The expression of most auxin genes is regulated by auxin. The expression of SI-IAA27 alters auxin sensitivity. The down-regulated expression or silencing of SI-IAA27 leads to increased auxin sensitivity, thereby improving root development. In contrast, the overexpression of this gene significantly changes the endogenous auxin gradient in roots and inhibits lateral root formation, ultimately resulting in relatively few roots [36]. In this study, the PsIAA27 expression level was higher in R2 than in R1 and R3 (Figs. 8 and 13), which may have affected the concentration of endogenous auxin in the roots and inhibited the formation of lateral roots. This might also help to explain why there were few tree peony lateral roots. We speculated that PsIAA27 influences auxin sensitivity and polar transport. High expression of PsIAA27 could promote main root elongation, but had a certain inhibitory effect on lateral root development.

The protein encoded by the ARF gene can bind to auxin and activate the production of auxin-related factors, which then regulate the expression of specific genes [37,38]. ARF19 is important for root formation, it has shown that ARF19 can coordinate with ARF7 to regulate the expression of CHI and IAA14, thereby regulating lateral root development, and arf7 arf19 can cause severely impaired lateral root formation and abnormal gravitropism in both hypocotyl and root [39–41]. A previous study by Okushima et al. [42] confirmed that ARF19 contributes to the early stage of root formation, possibly by modulating the induction of the root primordium. Other investigations proved that ARF19 regulates the formation of Arabidopsis adventitious roots, indicating that ARF19 has a critical role related to root development [39,43,44]. However, other ARFs can stimulate root development in the absence of ARF19 [41]. In this study, PsARF19 was more highly expressed in R3 than in R2 and R1 (Fig. 13), and gradually increased with the rooting process, which may have influenced the formation of tree peony rooting. This gene is expected to be expressed during the tree peony root induction period.

5  Conclusions

In summary, our RNA-seq analysis involving one nucleotide database and five functional protein databases resulted in the annotation of 120,188 unigenes. We mainly focused on identifying genes that regulate root formation under sandy loam cultivation. In this study, auxin and the related genes were revealed to be important for the induction and formation of roots in tree peony. The expression of PsIAA27 and PsARF19 significantly affected roots development of tree peony, and had similar effects in seedlings, cuttings and grafted seedlings, but the specific mechanisms were different. The data presented herein indicated that PsARF19 might be a regulatory gene related to the formation of tree peony roots, whereas PsIAA27 might be involved in regulating auxin sensitivity during root development. However, the specific mechanisms underlying the functions of these genes will need to be thoroughly investigated in future studies. The findings of this study provide a theoretical basis for future research on the mechanism mediating peony root formation. They may also be relevant for optimizing tree peony tissue culture, cuttings and grafting. The generated transcriptome data represent an excellent resource for exploring the genes involved in tree peony root development.

Authorship: Study conception and design: Z.W. and S.H.; data collection: Y.S., X.M., L.S.; analysis and interpretation of results: Y.S., D.H., X.L., Y.S.; draft manuscript preparation: Y.S., W.S. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement: Data supporting the reported results will be available and provided upon request.

Funding Statement: This study was financially supported by the National Key Research and Development Program of China (Grant No. 2020YFD1000503), the Natural Science Foundation of China (Grant Nos. 31870698, 32001353), the Key Scientific and Technological Projects of Henan Province (Grant No. 202102110083), and the Science and Technology Program of Shanghai (Grant No. 21DZ1202000).

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.

Appendix

Supplementary Table 1: Sequence information