Comparative Transcriptome Analysis Reveals Stem Secondary Growth of Grafted Rosa rugosa ‘Rosea’ Scion and R. multiflora ‘Innermis’ Rootstock

Grafted plant is a chimeric organism formed by the connection of scion and rootstock through stems, so stem growth and development become one of the important factors to affect grafted plant state. However, information regarding the molecular responses of stems secondary growth after grafting is limited. A grafted Rosa plant, with R. rugosa ‘Rosea’ as the scion (Rr_scion) grafted onto R. multiflora ‘Innermis’ as the stock (Rm_stock), has been shown to significantly improve stem thickness. To elucidate the molecular mechanisms of stem secondary growth in grafted plant, a genome-wide transcription analysis was performed using an RNA sequence (RNA-seq) method between the scion and rootstock. Comparing ungrafted R. rugosa ‘Rosea’ (Rr) and R. multiflora ‘Innermis’ (Rm) plants, there were much more differentially expressed genes (DEGs) identified in Rr_scion (6887) than Rm_stock (229). Functional annotations revealed that DEGs in Rr_scion are involved in two Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways: the phenylpropanoid biosynthesis metabolism and plant hormone signal transduction, whereas DEGs in Rm_stock were associated with starch and sucrose metabolism pathway. Moreover, different kinds of signal transduction-related DEGs, e.g., receptor-like serine/threonine protein kinases (RLKs), transcription factor (TF), and transporters, were identified and could affect the stem secondary growth of both the scion and rootstock. This work provided new information regarding the underlying molecular mechanism between scion and rootstock after grafting.


Introduction
Grafting is an ancient plant asexual propagation technique. It is practiced in forest trees, fruit trees, vegetable crops, and in many ornamentals. There are many advantages of using grafted plants, such as yield increase, stress tolerance, and successive cropping, which have been well-studied for several decades [1,2]. The rootstock, i.e., the entire root system of a grafting plant, is often selected to enhance nutrient uptake and alter various physiological processes in the scion, such as biomass accumulation [3], fruit quality [4], and response to abiotic stresses (e.g., water deficit and salinity) [5,6]. Little published information regarding the effect of the scion on rootstock development after grafting is available, but the scion genotypes that confer rootstock vigor and root patterns have been identified [3]. Grafting techniques are increasing for many plants around the world. scion), Rm (R. multiflora 'Innermis' ungrafted), and Rr (R. rugosa 'Rosea' ungrafted). Stems were separately harvested from the seedlings, Stem samples of Rm_stock and Rm were collected at the same position as well as Rr_scion and Rr: Rm_stock (5 cm distance from the branch point, about 5-8 cm distance from the graft union), Rm (5 cm distance from the branch point), Rr_scion (between the tenth and twelfth leaf, about 5-8 cm above the graft union), Rr (between the tenth and twelfth leaf). Rr_scion, Rr, Rm_stock, and Rm as samples for transcriptome analysis. These were quickly frozen in liquid nitrogen and then stored at −80 • C until further use. All experiments included three to five grafted plants under each treatment, and all were repeated three times. The stem diameter (the sixth internode) was measured in Rr_scion and Rr; for Rm_stock and Rm, the stem (about 5 cm below the grafting point) and root neck were measured. In addition, the lignin and cellulose content of stem dry matter in Rm_stock, Rm, Rr, and Rr_scion were determined by visible spectrophotometry (Solarbio, Beijing, China). Between Rr_scion and Rr, Rm_stock and. Rm, statistical analysis to find significant differential content of both lignin and cellulose was determined using a two-tailed Student's t-test in Microsoft Office Excel 2017 (p-values < 0.05 α-level).

Histological Analysis
The stems of Rr_scion, Rr, Rm_stock, and Rm were cut at the same measurement position. Tissue sections were prepared using a slicer (VT1200, Leica, Wetzlar, Germany) and then stained with 0.05% toluidine blue for 1 min to visualize secondary xylem tissues. The images were recorded using a Zeiss Axioplan light microscope (Zeiss, Jena, Germany), and images were captured using an axiocam digital camera (Zeiss) and AXIOVISION v.4.5 software (Carl Zeiss Vision GmbH, Hallbergmoos, Germany). The xylem width was measured from the border of the cambium flanking the xylem to the outermost cells of the pith [18].

Construction of the cDNA Library and Solexa Sequencing for Transcriptomic Analysis
For all samples, total RNA was extracted using TRIzol (Invitrogen, Carlsbad, CA, USA) as per the manufacturer's protocol. The quality and concentration of each RNA sample was assayed using 1% agarose gel electrophoresis, a spectrophotometer (K5500, Kaeo, China), and a bioanalyzer (2100 RNA Nano 6000 assay kit, Agilent Technologies, Santa Clara, CA, USA). Poly(A) mRNA was isolated using magnetic oligo(dT) beads, and the isolates were then divided into short fragments (Ambion, Austin, TX, USA). The first-strand cDNA was synthesized using random hexamer primers, and the second-strand cDNA was synthesized using dNTPs, buffer, RNaseH (Invitrogen), and DNA polymerase I (New England Biolabs, Ipswich, MA, USA). Short fragments were purified, subjected to end repair and the addition of poly(A), and were then ligated to sequencing adapters. The fragments of interest were purified with agarose gel electrophoresis and enriched using polymerase chain reaction (PCR) amplification. Finally, all cDNA libraries were sequenced using an Illumina Hiseq4000 platform with a PE150 sequencing strategy.

Sequencing, Assembly, and Functional Annotation of cDNA
Raw data (Raw reads) was firstly processed by with Perl Scripts. In this step, clean data (reads) of fastq format were obtained by removing reads with any of the following: (1) adapter contamination > 5 bp; (2) both more than a 15% base calling accuracy, a Phred quality score Q ≤ 15; (3) reads with number of N base accounting for more than 5%; As for paired-end sequencing data, both reads will be filtered out if any read paired-end reads are adaptor-polluted. Transcriptome assembly was accomplished using Trinity method (20140717) [19]. The obtained clean data after filtering will be carried out on statistical analysis, Q30, GC-content and sequence duplication level of the clean data were calculated . We predicted open reading frames (ORFs) of the assembled transcripts using TransDecoder (v20140717) based on the following criteria: a minimum length of 100 amino score and greater than 0 is reported; if a shorter ORF is fully encapsulated by a longer ORF, the longer one is reported; Trinotate (20140717) was used for performing the functional annotation of ORFs on the base of database. We predicted the unigene homology using the BLAST database (2.2.28), protein signal peptide and transmembrane domain prediction using SignalP (4.1) and TmHMM database (2.0), and then compared these outputs to current annotation databases, UniProt protein database, eggNOG database (4.5.1). The results were controlled by an e-value threshold of 10 −5 . Gene ontology (GO) annotations of contigs were determined using Blast2GO (http://www.blast2go.com/) according to molecular function, biological process, and cellular component ontologies (http://www.geneontology.org/) [20]. Here, we used an e-value threshold of 10 −5 . Annotated contigs were used to query the KEGG to define the KEGG orthologs (KOs). The KEGG mapping tool was used to plot these KOs into the complete metabolic atlas.

Identification of DEGs and Functional Analysis
The read counts for each gene in each sample was counted by high-throughput sequencing data (HTseq-count) (v.0.6.0). We applied the DESeq2 (V.1.4.5) method to analyze the expression of contigs between the test and control treatments of our experiment, and a model based on the negative binomial distribution was performed for normalization [21]. The p-value was assigned using Benjamini and Hochberg's approach for controlling the false discovery rate [22], The DEGs were selected if their |log 2 Ratio| ≥ 1, and q < 0.05.
DEGs were functionally annotated using the UniProt, Pfam, GO, and KEGG databases. A GO enrichment analysis was performed using Blast2GO, and we detected which of the DEGs were significantly enriched in GO terms using q < 0.05 as the threshold of significance. The cellular metabolisms, biochemical pathways, and biological potential of DEGs in the KEGG pathway were analyzed, and the enrichment pathways of DEGs were also assessed using a significance threshold of q < 0.05.

A Quantitative Reverse-Transcription PCR (qRT-PCR) for the Validation and Analysis of Expression Patterns
Total RNA was extracted from the stems of the different samples (Rm_stock, Rm, Rr_scion, and Rr), as described above. Reverse transcription was performed using Superscript II reverse transcriptase (Invitrogen) according to the manufacturer's instructions. Reverse transcription was performed using oligo (dT) primers and 1 µg of total RNA. The primer sequences are listed in Supplementary  Table S1. To detect transcript abundance, qRT-PCR was performed using an ABI Prism 7500 system (Applied Biosystems, Foster City, CA, USA) and a SYBR Premix Ex Taq kit (TaKaRa, Tokyo, Japan). The qRT-PCR was performed in 20-µL volumes containing 2 µL first-strand cDNA, 200 nM of each primer, and 10 µL of the 2× SYBR PCR mixture, with the following cycling parameters: 95 • C for 30 s, 40 cycles of denaturation at 95 • C for 3 s, and annealing and extension at 60 • C for 30 s. Three replicates were conducted in parallel, and the results were normalized differentially to the expression level of constitutive actin (Rm) or GAPDH (Rr). A relative quantitative method (∆∆Ct) was used to evaluate quantitative variation [23].

Effects of Rosa Grafting on the Stem Growth of Scion and Stock
To reveal the effect of grafting on the growth of scion and rootstock and on the stem thickness of Rr_scion (R. rugosa 'Rosea' as scion) and Rm_stock (R. multiflora 'Innermis' as stock), a phenotype analysis of a cross-section and of the lignin and cellulose content was conducted 1 year after grafting. The stem thickness of Rr_scion (the sixth internode) increased by 47.32% compared with Rr, and the stem (about 5 cm below the grafting point) and stem neck of Rr_scion increased by 291.98% and 206.12%, respectively, compared with Rm ( Figure 1). There were no differences in the xylem morphological characteristics in Rr_scion vs. Rr and Rm_stock vs. Rm, but the xylem width in the stem increased by 127.27% and 49.50% in Rr_scion and Rm_stock, respectively ( Figure 1). The lignin and cellulose content increased in Rr_scion and Rm_stock; the lignin content was significantly increased, by 42

Illumina Sequencing and De Novo Assembly of the Grafted Rosa Transcriptome
To investigate the genes associated with grafting response, four cDNA libraries were constructed from total RNA extracted from stems of Rr, Rr_scion, Rm and Rm_stock. The libraries were then sequenced by the Illumina Hiseq4000 platform. Summary of sequence assembly after illumina squencing was shown in Table 1. A tatal of 184,747 transcripts from Rr and 154,572 transcripts from Rm were obtained from clean data. In addition, 136,293 unigenes from Rr and 184,747 unigens from Rm were obtained with an average length of 816 bp and 859bp ( Table 2; Table  S2). The RNA-seq data can be found in the NCBI Sequence Read Archive (SRA) under ID number PRJNA406996, and the accession of the assembly cDNA data in Rm and Rr are SRR11091586 and SRR11091585. The seedlings of R. rugosa 'Rosea' (Rr), R. multiflora 'Innermis' (Rm), and a grafted plant (Rr grafted on Rm). (B) An analysis of stems between R. rugosa 'Rosea' and R. rugosa 'Rosea' as the scion, R. multiflora 'Innermis' and R. multiflora 'Innermis' grafted, and the root neck between R. multiflora 'Innermis'(Rm) and R. multiflora 'Innermis' grafted (Rm_stock). Rr (R. rugosa 'Rosea'); Rr_scion (R. rugosa 'Rosea' grafted); Rm-s (the stem of R. multiflora 'Innermis'); Rm-stock-s (the stem of R. multiflora 'Innermis' grafted); Rm-g (the root-neck of R. multiflora 'Innermis'); Rm-stock-g (the root-neck of R. multiflora

Illumina Sequencing and De Novo Assembly of the Grafted Rosa Transcriptome
To investigate the genes associated with grafting response, four cDNA libraries were constructed from total RNA extracted from stems of Rr, Rr_scion, Rm and Rm_stock. The libraries were then sequenced by the Illumina Hiseq4000 platform. Summary of sequence assembly after illumina squencing was shown in Table 1. A tatal of 184,747 transcripts from Rr and 154,572 transcripts from Rm were obtained from clean data. In addition, 136,293 unigenes from Rr and 184,747 unigens from Rm were obtained with an average length of 816 bp and 859bp ( Table 2; Table S2). The RNA-seq data can be found in the NCBI Sequence Read Archive (SRA) under ID number PRJNA406996, and the accession of the assembly cDNA data in Rm and Rr are SRR11091586 and SRR11091585. Compared with ungrafted samples, there were 6877 and 229 DEGs in Rr_scion and Rm_stock, respectively (q < 0.05). In the Rm_stock vs. Rm comparison, 89 DEGs showed increased transcript abundance, and 140 DEGs exhibited decreased transcript abundance (q < 0.05). In the Rr_scion vs. Rr comparison, 4600 DEGs were upregulated and 2277 DEGs were downregulated (q < 0.05). There were 4511 upregulated DEGs and 2137 decreased transcripts in Rr_scion, more than in Rm_stock, when compared with their control samples. (Figure 2A,B, Table S3). The heatmap of DEGs in Rr_scion vs. Rr and Rm_stock vs. Rm was listed in Figure S1. We also found that 67.25% of DEGS (154 of 229) were successfully annotated, and 22.75% (75 of 229) were unannotated in Rm_stock vs. Rm, whereas 71.62% of DEGs (4925 of 6877) were annotated and 28.38% (1945 of 6877) were unannotated in Rr_scion vs. Rr (Table S4).

GO Enrichment Analysis of DEGs
After conducting a GO analysis using blast2GO, the contigs of Rr and Rm were classified into 67 terms, and "cellular process," "cell part," and "binding" were all dominant among the each of the three main GO classifications "biological process", "cellular component" and "molecular function", ( Figure S2). To better understand the function of genes that affected growth and development after grafting, we further analyzed the GO terms of the DEGs enriched between Rr_scion vs. Rr and Rm_stock vs. Rm (Table S5). In the biological process category, these biological processes were mainly associated with "response to stimulus or stress," "developmental process," and

GO Enrichment Analysis of DEGs
After conducting a GO analysis using blast2GO, the contigs of Rr and Rm were classified into 67 terms, and "cellular process," "cell part," and "binding" were all dominant among the each of the three main GO classifications "biological process", "cellular component" and "molecular function", ( Figure S2). To better understand the function of genes that affected growth and development after grafting, we further analyzed the GO terms of the DEGs enriched between Rr_scion vs. Rr and Rm_stock vs. Rm (Table S5). In the biological process category, these biological processes were mainly associated with "response to stimulus or stress," "developmental process," and "multicellular organismal process" between Rr_scion and Rr (p < 0.05) ( Figure 3; Table S5). We also found that the GO terms "oxidoreductase activity," "oxygen oxidoreductase activity," "cellulose synthase," and "protein serine/threonine kinase activity" in the molecular function category and "plasma membrane" in the cellular component category were enriched in Rr_scion vs. Rr. Among Rm_stock and Rm (control) samples, most of the DEGs enriched in biological process terms were involved in "metabolic process" (GO:0008152, 70 DEGs) and "cellular process" (GO:0009987, 67 DEGs) ( Figure 3; Table S5). These cellular components were associated with "intrinsic component of membrane" (GO:0031224, 43 DEGs), and the molecular functions "cytokinin dehydrogenase activity" (GO:0019139, 2DEGs) and "transmembrane transporter activity" (GO:0022857,19DEGs) were observed in Rm_stock vs. Rm. The enriched GO analysis also showed that the scion and rootstock responded differently to grafting (Table S5).
Genes 2020, 11, x FOR PEER REVIEW 7 of 18 "multicellular organismal process" between Rr_scion and Rr (p < 0.05) ( Figure 3; Table S5). We also found that the GO terms "oxidoreductase activity," "oxygen oxidoreductase activity," "cellulose synthase," and "protein serine/threonine kinase activity" in the molecular function category and "plasma membrane" in the cellular component category were enriched in Rr_scion vs. Rr. Among Rm_stock and Rm (control) samples, most of the DEGs enriched in biological process terms were involved in "metabolic process" (GO:0008152, 70 DEGs) and "cellular process" (GO:0009987, 67 DEGs) ( Figure 3; Table S5). These cellular components were associated with "intrinsic component of membrane" (GO:0031224, 43 DEGs), and the molecular functions "cytokinin dehydrogenase activity" (GO:0019139, 2DEGs) and "transmembrane transporter activity" (GO:0022857,19DEGs) were observed in Rm_stock vs. Rm. The enriched GO analysis also showed that the scion and rootstock responded differently to grafting (Table S5).

Protein Kinases
In our study, 103 DEGs were predicted to encode protein kinases according to the functional annotation, of which 96 and 7 contigs encoding protein kinases were differentially expressed in Rr_scion vs. Rr and Rm_stock vs. Rm samples, respectively ( Figure 5, Table S7). In Rr_scion vs. Rr, we found that 89.58% of DEGs (86 of 96) encoding protein kinase were upregulated, and 10.42% (10 of 96) were downregulated. We further identified 26 DEGs encoding leucine-rich repeat receptor-like serine/threonine protein kinase (LRR-RLK). Another 32 DEGs encoding cysteine-rich receptor-like protein kinase (CRK), wall-associated receptor kinase (WAK), and receptor-like protein kinase were all differently upregulated; of these, 15 DEGs (including c160972_g4, c177139_g1, c112384_g1, and c112384_g1) were upregulated more than 20-fold ( Figure 5, Table S7). For Rm_stock vs. Rm, three categories of DEGs encoding LRR-RLK, receptor protein kinase, and receptor-like cytosolic serine/threonine-protein kinase, were identified, and the four DEGs encoding LRR-RLK were downregulated differently ( Figure 5, Table S7). The contrasting expression patterns of LRR-RLK indicated that they were differentially modulated in Rr_scion vs. Rr and Rm_stock vs. Rm.
vs. Rm, three categories of DEGs encoding LRR-RLK, receptor protein kinase, and receptor-like cytosolic serine/threonine-protein kinase, were identified, and the four DEGs encoding LRR-RLK were downregulated differently ( Figure 5, Table S7). The contrasting expression patterns of LRR-RLK indicated that they were differentially modulated in Rr_scion vs. Rr and Rm_stock vs. Rm.   vs. Rm, three categories of DEGs encoding LRR-RLK, receptor protein kinase, and receptor-like cytosolic serine/threonine-protein kinase, were identified, and the four DEGs encoding LRR-RLK were downregulated differently ( Figure 5, Table S7). The contrasting expression patterns of LRR-RLK indicated that they were differentially modulated in Rr_scion vs. Rr and Rm_stock vs. Rm.

Transcription Factors (TFs)
A number of TFs have been reported to play important roles in the secondary growth of stems. In our study, 96 differentially expressed transcripts were identified belonging to 15 TF categories in Rr_scion vs. Rr, and 16 transcripts belonging to 12 TF categories were identified in Rm_stock vs. Rm ( Figure 6, Table S8). In Rr_scion vs. Rr, the largest number of DEGs was found for WRKY TFs (24 DEGs), followed by zinc finger protein (19 DEGs) and MYB (14 DEGs). We also identified 76 DEGs encoding TFs as being upregulated in Rr_scion vs. Rr, whereas 24 were found to be downregulated. In addition, 94% of DEGs (15 of 16) encoding MYB TFs, DOF TFs, and others were upregulated in Rm_stock vs. Rm ( Figure 6, Table S8).

Transcription Factors (TFs)
A number of TFs have been reported to play important roles in the secondary growth of stems. In our study, 96 differentially expressed transcripts were identified belonging to 15 TF categories in Rr_scion vs. Rr, and 16 transcripts belonging to 12 TF categories were identified in Rm_stock vs. Rm ( Figure 6, Table S8). In Rr_scion vs. Rr, the largest number of DEGs was found for WRKY TFs (24 DEGs), followed by zinc finger protein (19 DEGs) and MYB (14 DEGs). We also identified 76 DEGs encoding TFs as being upregulated in Rr_scion vs. Rr, whereas 24 were found to be downregulated. In addition, 94% of DEGs (15 of 16) encoding MYB TFs, DOF TFs, and others were upregulated in Rm_stock vs. Rm ( Figure 6, Table S8).

Transporter Genes
Many DEGs encoding different types of transport factors were differentially expressed. In Rm_stock, 19 DEGs were related to the transport of vacuolar amino acids, ABC, sugar, boron, inorganic phosphate, potassium, metals, folate-biopterin, xyloglucan endotransglucosylases/hydrolases, and lipids. The expression of a vacuolar amino acid transporter gene (contig c130686_g1) showed a 15-fold increase in Rm_stock relative to Rm. In Rr_scion, 69 DEGs encoded 18 different categories of transport proteins (Table S9). A total of 16 ABC transporters were significantly differentially expressed in Rr_scion vs. Rr.

qRT-PCR Validation of Differentially Expressed Transcripts from RNA-seq
To confirm the accuracy and reproducibility of our Illumina RNA-seq results, a small number of unigenes, including TFs and hormone-associated genes, were randomly chosen for qRT-PCR amplification. The primer sequences for these 26 unigenes are listed in Table S1. The relative expression levels of these unigenes were analyzed, yielding general agreement in their transcript abundance, as determined by RNA-seq (Figure 7), thus corroborating our RNA-seq transcriptomic data. For example, the contigs c120633_g1 and c141086_g1 in Rm_stock and Rm and the contigs c112414_g1 and c52312_g1 in Rr_scion and Rr showed similar relative differences following the

Transporter Genes
Many DEGs encoding different types of transport factors were differentially expressed. In Rm_stock, 19 DEGs were related to the transport of vacuolar amino acids, ABC, sugar, boron, inorganic phosphate, potassium, metals, folate-biopterin, xyloglucan endotransglucosylases/hydrolases, and lipids. The expression of a vacuolar amino acid transporter gene (contig c130686_g1) showed a 15-fold increase in Rm_stock relative to Rm. In Rr_scion, 69 DEGs encoded 18 different categories of transport proteins (Table S9). A total of 16 ABC transporters were significantly differentially expressed in Rr_scion vs. Rr.

qRT-PCR Validation of Differentially Expressed Transcripts from RNA-seq
To confirm the accuracy and reproducibility of our Illumina RNA-seq results, a small number of unigenes, including TFs and hormone-associated genes, were randomly chosen for qRT-PCR amplification. The primer sequences for these 26 unigenes are listed in Table S1. The relative expression levels of these unigenes were analyzed, yielding general agreement in their transcript abundance, as determined by RNA-seq (Figure 7), thus corroborating our RNA-seq transcriptomic data. For example, the contigs c120633_g1 and c141086_g1 in Rm_stock and Rm and the contigs c112414_g1 and c52312_g1 in Rr_scion and Rr showed similar relative differences following the qRT-PCR and RNA-seq. Overall, we found similarity in the expression trends between our qRT-PCR and RNA-seq results (r 2 = 060), as shown in Figure 7.
Genes 2020, 11, x FOR PEER REVIEW 10 of 18 qRT-PCR and RNA-seq. Overall, we found similarity in the expression trends between our qRT-PCR and RNA-seq results (r 2 = 060), as shown in Figure 7.

Stem Secondary Growth Response for Grafted R. Multiflora 'Innermis'/R. Rugosa 'Rosea' Plants
Rosa is one of the most popular ornamental plants in the world, and the mechanisms of growth, development, and resistance in Rosa are currently the foci of several research projects. The genotypes of the rootstock or scion are often selected for their ability to alter the growth and development of a grafted plant. In our study, the secondary growth (stem thickness) of the scion R. rugosa 'Rosea' significantly increased compared to an ungrafted control. This result was similar to that of previous studies on grafted apples, grapevines, and watermelons compared with a self-grafted control [11,24,25]. In addition, the stem thickness of the rootstock was significantly increased in the Rm_stock vs. Rm, suggesting that the scion biomass also conferred vigor to the rootstock, similar to the scion genotype's effect on the shoot and root in grated grape [3].

Stem Secondary Growth Response for Grafted R. Multiflora 'Innermis'/R. Rugosa 'Rosea' Plants
Rosa is one of the most popular ornamental plants in the world, and the mechanisms of growth, development, and resistance in Rosa are currently the foci of several research projects. The genotypes of the rootstock or scion are often selected for their ability to alter the growth and development of a grafted plant. In our study, the secondary growth (stem thickness) of the scion R. rugosa 'Rosea' significantly increased compared to an ungrafted control. This result was similar to that of previous studies on grafted apples, grapevines, and watermelons compared with a self-grafted control [11,24,25]. In addition, the stem thickness of the rootstock was significantly increased in the Rm_stock vs. Rm, suggesting that the scion biomass also conferred vigor to the rootstock, similar to the scion genotype's effect on the shoot and root in grated grape [3].

DEGs in Key Pathway about Lignin and Cellulose Biosynthesis in the Scion and Rootstock
Lignins are complex phenolic polymers of plant cell walls, and reductions in lignin during stem secondary development can cause growth defects and affect water and nutrient transport in plants [26]. Key insights into the molecular aspects of gene regulation in grafted plants can be attained by analyzing DEGs [24]. The phenylpropanoid biosynthesis pathway (map00940), which is related to lignin biosynthesis, was significantly enriched, and DEGs encoding enzymes implicated in lignin synthesis were identified in Rr-scion vs. Rr (Table S10). Peroxidase (PRX), AtPrx2, AtPrx25, and AtPrx71 are involved in lignin biosynthesis, and a significant decrease in the total lignin content has been reported in ATPRX2-and ATPRX25-deficient mutants [27]. These DEGs were confirmed in Rr_scion vs. Rr and Rm_stock vs. Rm (Supplementary Table S10) and were upregulated by 1.34-5.39-fold in both of them. Interestingly, DEGs (60%, 6 DEGs) encoding PRX were downregulated in Rr_scion vs. Rr, which may have regulated the stem lignin structure of the scion [28]. Cinnamyl alcohol dehydrogenase (CAD) catalyzes the last step in the monolignol biosynthesis pathway, and the expression of At4g34230 encoding CAD in Arabidopsis has been reported to increase eight-fold, with very high signal intensity from immature to mature stems [29]. In our study, two DEGs (c158813_g2 and c133940_g1) encoding CAD were upregulated in Rr_scion vs. Rr, however, downregulated in Rm_stock vs. Rm. In addition, we also found that many DEGs encoding 4-coumarate-CoA ligase (CCoAOMT), 4-coumarate-CoA ligase (4-CL) were upregulated in Rr_scion vs. Rr (Table S10). The lignin content increase of Rr_scion also implies that was maybe affected by these genes after grafting (Figure 1) [30].
Cellulose is a major important component of plant cell wall as well as lignin, which are responsible for both oriented cell elongation during plant growth and the strength required to maintain an upright growth habit [31]. Cellulose is a simple polymer of unbranched β-1,4-linked glucan chains [32]. In our study, five DEGs encoding α, α-trehalase, 1,4-α-glucan branching enzyme, and β-amylase were enriched in the "Starch and sucrose mechanism" pathway. It implies the transformation of starch to soluble sugars at the transcriptional level to involved in cellulose synthesis in Rm_stock vs. Rm. In addition, we found DEG (c147623_g4) encoding xyloglucan endotransglucosylase/hydrolase protein (XTH) were upregulated 2.39-fold in Rm_stock vs. Rm, which has been recognized as a cell wall-modifying enzyme. Addition to purified XTH enzyme to Arabidopsis will act predominantly to strengthen or "tighten" cells, or loosen cell walls [33,34]. These DEGs were maybe regulating the stem thickness in Rm_stock vs. Rm, by changing the content of energy production including sugar, cellulose and starch, which need further to be identified [34].

DEGs in Response to Phytohormone Signal Transduction in the Scion and Rootstock
Hormones have been widely applied to regulate secondary vascular growth. Auxin, cytokinins, brassinosteroids, gibberellins, abscisic acid, and ethylene have all been found to control cambial growth and differentiation through a complex regulatory network [35]. In our study, a number of DEGs related to the six hormones mentioned above were identified, and their expression levels were also analyzed in Rr_scion vs. Rr and Rm_stock vs. Rm (Table 3. We found many more DEGs related to hormone signal transduction (e.g., auxin-responsive GH3, EIN3-binding F-box protein EBF, and ethylene receptor ETR) than involved in hormone synthesis pathways (e.g., gibberellin 20 oxidase, cytokinin, and dehydrogenase). This implies that enhanced hormone signal transport and crosstalk may be the main outcome of grafting, rather than increased hormone content. In addition, auxin appeared to play a primary role in vascular patterning in Rr_scion vs. Rr, and the transcriptomes of 18 genes (including AUX/IAA, auxin binding protein, and auxin-induced protein) were significantly changed. They were found to be highly and preferentially expressed in xylem and during the differentiation of tracheary elements [36], and the loss and gain of function resulted in reduced and discontinuous vascular formation, respectively [37,38]. For Rm_stock vs. Rm, the DEGs were more related to cytokinin and indole-3-acetic acid synthesis in the zeatin biosynthesis pathway (map00908). Cytokinin has been recognized as key regulator of cambial activity, procambium maintenance, and development of the vasculature [39]. Degradation of cytokinin is catalyzed by cytokinin oxidase/dehydrogenase (CKX) enzymes. It was identified that ckx3 ckx5 and ckx7 regulated the activity of the reproductive meristems of Arabidopsis [40,41], It was found the inflorescence stem of ckx3 ckx5 mutants is thicker than the wild type, and CKX7-overexpressing in transgenic Arabidopsis reduce the meristem size of roots and induce the cessation of root growth [41]. In our study, two DEGs (c143974_g1 and c145655_g1) were upregulated or downregulated in the stem of Rm_stock vs. Rm, which implied to regulate the development of stem, and the specific regulating should be identified in the future study. In addition, some findings have indicated an important role for long-distance basipetal transport of cytokinin through the phloem in controlling vascular patterning in roots via inhibitory interaction with auxin [42,43], and whether the stem vascular development was affected by cytokinin transportation should be explored. The different expressions of hormone-related DEGs suggested distinct responses in Rr_scion and Rm compared to the control.

Protein Kinases Are Related to Stem Vascular Development in the Scion and Rootstock
Protein kinases (PKs) play an important role in cellular singular transduction processes during plant growth and development [44,45]. In Arabidopsis, 223 LRR-RLK gene expression patterns were identified that had specific characteristics at various developmental stages [46]. We compared the DEGs in Rr_scion vs. Rr and Rm_stock vs. Rm with 14 LRR-RLK subfamilies in Arabidopsis (Figure 7). For Rr_scion vs. Rr, the DEGs (C157420_g1, C152554_g1) were separately assembled into LRRX and LRRIX subfamilies in Arabidopsis, and were upregulated by 3.92 and 12.58-fold ( Figure S3A), resulting in increased stem growth and vascular development and increasing the length of petioles and hypocotyl growth due to brassinosteriod activity when they were overexpressed in transgenic Arabidopsis [47,48]. Genes from the LRR-RLKXI subfamily were found to reduce plant growth and affect vascular development when loss of either RLK PHLOEM INTERCALATED WITH XYLEM (pxy) or RLK XYLEM INTERMIXED WITH PHLOEM1 (xip1) occurred. Interestingly, the homologous gene (c144124_g2) was significantly downregulated, by 4.76-fold, in Rm_stock vs. Rm, and was not consistent with the stem thickening growth of the rootstock ( Figure S3B). The other DEGs encoding LRR-RLKs (except c9455_g1) were also all downregulated in Rm_stock vs. Rm, suggesting that these genes were negatively related to plant growth and development. In addition, many other protein kinases, such as WAKs, MAPK, LecRKs, GsSRK, and CRKs, were identified as DEGs. Their function in stem secondary growth requires further study.

Transcription Factors Involved in Stem Secondary Growth in the Scion and Rootstock
Furthermore, there is increasing evidence that various TF families, such as MYB, NAC, and WRKY, are involved in the regulation of stem secondary growth [26,49]. Our data provide evidence that some candidate TFs were involved in stem secondary growth after grafting. The NAC and MYB TF genes are the key players regulating the complex transcriptional network leading to wall-thickening cell differentiation [50,51]. We further analyzed NAC and MYB DEGs in our study according to the phylogenetic tree domain in Arabidopsis. Two DEGs (c157108_g1 and c130114_g1) were homologous with AtVND1, AtVND7, and AtXND1 in Rr_scion vs. Rr ( Figure S3C); these play important roles in vessel formation [52,53]. Interestingly, the DEG c109220_g1 encoding MYB TF in Rr_scion vs. Rr and c119611_g1 in Rm_stock vs. Rm were both found to be homologous with AtMYB75 ( Figure  S3E,F), which negatively regulated the secondary cell wall (SCW) [54]. It has been suggested that the NAC-MYB-based transcriptional network can regulate the stem secondary growth of Rr_scion and Rm_stock after grafting [55]. In addition, TFs such as the AP2-EREBP, bHLH, C2H2, C2C2-GATA, and GRAS gene families have been found to regulate secondary cell wall metabolic genes through the protein-DNA network in Arabidopsis thaliana [56]. In our study, 9 bHLH family, 10 C2H2 family, and 7 HD-ZIP family TFs were identified in Rr_scion vs. Rr, and the C2C2-GATA and GRAS gene families were identified in Rm_stock vs. Rm, providing a number of candidate regulators of stem growth after grafting.

Transporter Proteins Are Important for Stem Secondary Growth in the Scion and Rootstock
Substantial DEGs encoding transporters such as ABC transporter, sugar, boron, inorganic phosphate, and potassium have been found in Rr_scion vs. Rr and Rm_stock vs. Rm. In plants, they are an important membrane transport protein family, relevant to the transportation of hormones, lipids, metal ions, and exogenous substances [57]. Some ABC transporters, especially AtABCB (AtABCB1 and AtABCB19) and AtABCG (AtABCG14, AtABCG25, and AtABCG40) have been reported to serve a function in stem lignification by auxin, tZ-type cytokinins, and ABA transport and distribution [58,59]. Most DEGs (16 DEGs in Rr_scion; 3 DEGs in Rm_stock) encoding ABC transporter were identified in Rr_scion vs. Rr and Rm_stock vs. Rm. The DEGs (c124372_g1, c150319_g3, and c154902_g1) were homologous with AtABCG according to the phylogeny tree ( Figure S3H) in Rr_scon vs. Rr; they seemed to function as hormone transporters, affecting the secondary growth of stems after grafting, and should, therefore, be further studied. In addition, sugar transporters (16 DEGs) were more strongly differentially expressed in Rr_scion vs. Rr and Rm_stock vs. Rm. Many sugar transporters localized in phloem companion cells and the associated parenchyma in maturing stems were found to affect vascular development [60,61]. Transporter DEGs may suggest that the growth and development of Rosa grafted plants were altered by the transportation of complex organic materials in Rr_scion and Rm_stock. This has been verified in grafted watermelon plants, where many DEGs are responsible for water transport [24]. Five DEGs were enriched in the KEGG pathway "Starch and sucrose metabolism" (map00500) in Rm_stock vs. Rm, which implies that changes in sugar content may affect the stem growth of rootstock.
Recently, regulatory networks activated in response to different developmental stages and various abiotic and biotic stimuli have been identified [62,63]. Stress responses are likely integrated into the gene regulatory network that determines xylem cell specification and differentiation [56]. Hydrogen peroxide (H 2 O 2 ), one of the stress-induced reactive oxygen species (ROS), also plays an important regulatory role in lignin biosynthesis [64]. Ascorbate peroxidase (APX) and CuZn-superoxide dismutase (CuZn-SOD) genes have been found to positively regulate secondary wall biosynthesis and promote growth in Arabidopsis [65]. In addition, the stem lignin metabolism is known to be related to plant disease resistance, insect resistance, and tolerance of drought, salt, heat, cold, heavy metals, and other stresses [66]. In our study, 845 DEGs significantly enriched in 25 GO terms were identified in the biological process "Response to stimulus or stress" in Rr_scion vs. Rr (p < 0.05), and 17 DEGs were also involved in Rm_stock vs. Rm. How these stress DEGs are regulated by grafting to improve stem secondary growth requires further investigation. In addition, many related abiotic stimuli DEGs have been identified in other grafted plants such as watermelon and grapevine [11,24].

Conclusions
Grafting is particularly important for plant development. This study used transcriptome data for the first time to analyze the stem secondary growth of the scion and rootstock after grafting. The scion R. rugosa 'Rosea' and rootstock R. multiflora 'Innermis' stem morphology were observed, and the lignin and cellulose content were measured. These were found to be significantly changed compared to an ungrafted control. A total of 136,293 unigenes in Rr and 108,651 unigenes in Rm were obtained, and 6,877 and 229 DEGs were detected, respectively. The data revealed important pathways for stem secondary growth of the scion and rootstock after grafting, such as "Plant hormone signal transduction" and "Phenylpropanoid biosynthesis," consistent with morphological changes in the scion, and "starch and sucrose metabolism" in the rootstock. Substantial signal transduction genes such as PK, TF, and transporters were found to regulate the secondary cell wall of the stem and may be important determinants of the underlying stem secondary growth. These results will facilitate future analysis of the roles of these genes in stem secondary growth after grafting and will also be useful for enabling rootstock breeding with thicker stems.