Integrated analysis of transcriptome, small RNA, and phytohormonal content changes between Artemisia annua Linn. and Nicotiana benthamiana Domin in heterogeneous grafting

As one of the most important techniques in horticulture, in China, grafting dates from 7000 BC^[1]. During the natural grafting process, people found that scion and rootstock from two different species could join together, the interface junction named the callus, which subsequently differentiates into vascular tissue to refactor material translocation^[2,3]. With the improvement in the success rate of plant grafting, the technique is now widely used in fruit, crop and vegetables, for propagating superior strains^[4], reducing disease and insect infection rates, and increasing crop yields^[5].

Population compatibility between scion and rootstock is the most important prerequisite for survival of grafted species, which includes the growing and regeneration of callus, and establishment of symbiotic relationships between scion and rootstock^[6]. In general, grafting survival rate of intraspecies is higher than interspecies, which means grafting incompatibility increases with taxonomic distance^[7]. Compared with incompatible graft combination, compatible combination could induce genes expression involved in auxin signal transduction, wound healing, and secondary metabolism^[8]. Besides, proteomics differences can also reflect the physiological status between compatible-rootstock grafting and incompatible-rootstock grafting. Proteins related to the Calvin cycle, carbohydrate metabolism, photosynthetic organs are relatively enriched in graft-compatible cucumber/pumpkin^[6], which would contribute to superior photosynthetic capacity, growth performance and biomass. During the grafting process, the accumulation of phenolic compounds, especially coumaric acids and flavonoid, may result in incompatible graft phenomenon. Actually, successful graft combination could manifest several markers like transcripts changes, metabolite differences, and redox status^[9], and the markers can be used for screening affinity plants in horticultural research.

Micrografting experiments^[10] have revealed the healing process of graft union, including callus formation and then generating new vascular connections. During the first 3 d after grafting, scion and rootstock may attach to each other through adhesion of parenchyma. As the cells increased in value, callus formed between rootstock and scion, and then producing cell differentiation and organ regeneration^[11]. Through formed plasmodesmata^[3,12], scion and rootstock may produce a clear cell-to-cell transfer action, including not just RNA, small RNAs, metabolites, and protein, which are the focus point we have perpetually focused.

Messenger RNA (mRNA) is a critical role in the transmission of genetic information and regulating physiological function. In the tobacco/tomato grafted model^[13], the mobile mRNAs in phloem can be easily identified owing to its large phylogenetic distance. Meanwhile, in the grafting process of cucumber and pumpkin^[14], four mRNAs related to chilling-induced in pumpkin could move to chilling-tolerant cucumber, and then regulate the fatty acid β-oxidative degradation metabolism in its internal, thus promoting the ability to resist cold. It is expected that through the transfer of mRNAs between scion and rootstock, a new allopolyploid species may be created, thus providing a method to asexual reproduction.

Plant microRNAs (miRNAs) are a class of small non-coding RNAs, 21 to 23 nucleotides in length, which belong to pivotal regulatory elements of gene expression^[15]. In heterografted plants, miRNA could produce unidirectional or bidirectional movement, especially from shoot to root, which included metabolic pathways and phytohormone signal transduction^[16,17]. The directional long distance movement of miRNAs could induce phenotypic changes of scion^[16], which illustrates that miRNAs play an important role in epigenetic inheritance. Through the movement of miRNAs, we could use the pattern for cultivating the modified crops of higher tolerance to drought and salt^[18]. Moreover, in grafted avocado^[19], a miR156-SPL4-miR172 model could act as a marker for plant affinity in response.

Whatever the model of homo- or hetero-grafting, it will both induce the damage response mechanisms in scion and rootstock. During the grafting process, endogenous hormones in plants always play its role in callus formation, and respond to biotic or abiotic stress^[20]. There are dozens of hormones and their derivatives in plants, which performs their respective functions in concert or antagonism. In grafted apples, the development of rootstock was not only affected by sugar metabolism, but the signal pathway of auxin and cytokinin also played an important role.

Artemisia annua, Artemisia of the Compositae family, it is well-known for producing artemisinin, a sesquiterpene lactone with an endoperoxide bridge^[21,22]. Artemisinin is a special compound for treating malaria used in Artemisinin-based Combination Therapy (ACT), specially synthesized in glandular trichomes of leaves and flowers^[23]. However, the content of artemisinin in dry weight of A. annua is 0.5%−1.2%, so it is still a scarce resource. At present, the artemisinin biosynthetic pathway has only been very basically understood. Through the action of the MVA pathway^[24] and the MEP pathway^[25], farnesyl diphosphate can be recognized as a key compound for producing artemisinin. Subsequently, ADS, CYP71AV1, DBR2, ALDH1 were reported as four critical genes for participating in artemisinin synthesis. In order to devise powerful methods to improve content of artemisinin, researchers mostly focused on various types of transcription factors^[26], phytohormone^[27]. Owing to the enormity and complexity of the A. annua genome, it is hard to obtain targeted transgenic plants or ideal mutants, which makes it difficult to research plant resistance or physiology.

Nicotiana benthamiana is a model plant of Solanaceae and commonly used for grafting analysis due to its detailed genomic information^[13,28]. Owing to the mature transient transformation system of Agrobacterium rhizogenes, transgenic N. benthamiana lines can be easily obtained for investigating mRNAs transmission in the case of graft. Through the interfamily grafting results of N. benthamiana with a variety of plants, β-1,4-glucanases produced by N. benthamiana was found to accomplish cell-cell adhesion and facilitate success rate of grafting^[29], thus N. benthamiana can be expected to be an excellent material for interspecific or interfamily grafting.

In this study, we developed a heterograft system using A. annua as scion and N. benthamiana as rootstock to investigate possible elements involved in graft compatibility. We used this model to carry out three parts: mRNA, miRNA, and phytohormone. The expression level and movement of mRNAs before and after grafting were compared, we can easily obtain the crucial genes that affect the normal growth and development of plants. The mobile mRNAs may play an important role in host plants, which may act as a basis for predicting its physiological function. Further investigation of miRNAs could supplement evidence of physiological activity in the Aa/Nb graft model. The regulatory network of miRNAs may help us better understand the physiological changes in heterograft plants. Meanwhile, changes in gene expression associated with phytohormones along with content analysis would provide a more accurate and detailed information for compatibility grafting. We summarized the changing trend of Aa/Nb heterografted plants at the molecular and metabolic level, not just find genes that make plants obtain better physiological status, but provide research materials in heterograft plants.

Seedlings of Artemisia annua cultivar (Aa)^[30] named as 'Huhao 1' with high artemisinin content and Nicotiana benthamiana (Nb) were sown on nutrient soil under a 16/8 h light/dark photoperiod at 25 °C of plant incubator. The nutrient soils are formulated in a ratio of 4:1:1 with peat, perlite and vermiculite and are slightly acidic, where the pH value is 6.5. Considering the long growth cycle of A. annua, the seeds of N. benthamiana were cultivated immediately when A. annua sprouted. The apical wedge grafting technique was performed (Fig. 1) as described by Buoso & Loschi^[31]. In addition, the leaves of A. annua scion and N. benthamiana rootstock were all cut off to avoid interference of old leaves during sampling. The grafted plants were still placed in plant incubator for growing together for 30 d. Through 30 d after grafting (DAG) growth, 16 successful Aa/Nb grafted plants were obtained.

Figure 1. 

A heterograft model of A. annua scion and N. benthamiana rootstock. (a) The upper half of A.annua and bottom half of N. benthamiana were cut using a sharp knife. Through the apical wedge grafting technique, the two part were combined and co-grown for a month. (b) The density of gladular trichome in leaves and (c) artemisinin content were calculated as shown in the histogram, no less than three biological replicates were performed for each data.

A. annua scion and N. benthamiana rootstock were respectively sampled from growing vigorously Aa/Nb grafted plants for RNA-seq^[32] and miRNA-seq^[33] (Supplemental Fig. S1). The samples were collected and frozen immediately in liquid nitrogen and stored at −80 °C to avoid RNA degradation. Each of the samples was derived from one single plant with at least three biological replicates.

Meanwhile, in order to elaborate the changes of phytohormone content in each tissue, the leaves and stem of A. annua scion, the leaves, stem and root of N. benthamiana rootstock were sampled for determining the changes of phytohormone content. The samples of phytohormone testing were firstly washed in distilled water for preventing disturbance of constituents, and then similarly deposited at −80 °C to prevent degradation. Besides, non-grafted A. annua and N. benthamiana plants with the same growth cycles were sampled as the control group. Each of the samples was derived from one single plant with at least three biological replicates.

Calculation of glandular trichome density and artemisinin content determination

After 30 DAG, we observed the whole plant phenotype, including glandular trichome of newborn Aa leaves, Nb leaves and Nb roots by fluorescence microscope^[34]. The density of trichomes of Aa leaves is calculated by ImageJ. The same part of the leaf was taken for each count, and there were at least three biological replicates. Due to the related artemisinin content, we used Aa leaves and stems to measure the content. Leaves and stems were deposited into the vacuum freeze drier for draining water. We used 100 mg of dry materials to extract the compound with methanol as a solvent, with ultrasound twice at 40 w. After centrifugation, the supernatant was taken and sampled through 0.22 um filter membrane. Samples were analyzed by HPLC-MS (1200-G6410A) with Aglient system in positive ion models^[35]. Besides, the mobile phase includes acetonitrile and water solution containing 0.1% formic acid. The content of artemisinin is based on 1, 10, 50, 100, 1,000 μg/ml. MassHunter software was used for identifying the sample correctness and calculating relative content.

RNA extraction, sequencing of mRNA and data processing

Total RNA of each sample was extracted from frozen materials using TRIzol reagent (Invitrogen, USA) and the quality of isolated RNA was guaranteed by Bioanalyzer 2100 (Agilent, USA) and NanoDrop ND-1000 (NanoDrop, USA). The integrity of RNA was confirmed through Bioanalyzer 2100 (Agilent, CA, USA) with RIN number > 7.0 and electrophoresis with denaturing agarose gel. And then mRNA contained PolyA were specifically captured using Dynabeads Oligo (dT)25-61005 (Thermo Fisher, CA, USA) by two rounds of purification. Through Magnesium RNA Fragmentation Module (NEB, cat.e6150, USA) under 94 °C 5−7 min, the poly(A) RNA was fragmented. After reverse transcribing by SuperScript^TM II Reverse Transcriptase, we got the final cDNA library with 300 ± 50 bp fragments. Finally, based on IIIumina Novaseq^TM 6000 (LC-Bio Technology, China), the library was sequenced through 2× 150 bp paired-end sequencing. The preparation of library and sequencing experiments of miRNA are identical with mRNA protocol as mentioned above. The reads contained adaptor contamination, low quality bases and undetermined bases were removed by https://github.com/OpenGene/fastp. Using HISAT2 software^[36] (https://daehwankimlab.github.io/hisat2/), the reads were mapped to the reference genome like A. annua and N. benthamiana. The mapped reads of each sample were assembled using StringTie (https://ccb.jhu.edu/software/stringtie) with default parameters. All mapped transcriptomes were merged to construct a comprehensive transcriptome by gffcompare (https://github.com/gpertea/gffcompare/). Through StringTie, the expression levels of all transcripts was estimated by calculating FPKM (FPKM = [total_exon_fragments / mapped_reads (millions) × exon_length (kB)]) and R package edgeR^[37]. Differentially expressed mRNAs were chosen with fold change > 2 or < 0.5, along with the parametric F-test comparing nested linear models (p value < 0.05) by R package edgeR (https://bioconductor.org/packages/release/bioc/html/edgeR.html). And then GO enrichment and KEGG enrichment were analyzed based on differentially expressed mRNAs by DAVID (https://david.ncifcrf.gov/).

In this project, we use fold change (FC ≥ 2 or FC ≤ 0.5, i.e. absolute value of log2FC ≥ 1) as change threshold, p < 0.05 as the standard for screening differential genes. Results of differential expression gene analysis, GO enrichment analysis and KEGG pathway enrichment analysis of differential expression gene were obtained in the set comparison group. To better explore the differential genes and potential mobile mRNAs, measured transcriptome data were compared with two databases: A. annua genome and N. benthamiana genome. Mobile mRNAs were confirmed seriatim in the way that genes were never present in the normal plant but were detected in grafted plants. For every sequence a comparison was made through Local BLAST and NCBI database.

Data analysis of miRNA and regulatory network of target genes

The validate library of miRNA was constructed using TruSeq Small RNA Sample Prep Kits (IIIumina, USA) and then sequenced through 1× 50 bp single-end method based on Illumina Hiseq2000/2500. The raw reads of miRNA were subjected into an in-house program named ACGT101-miR (LC, USA) to obtain clean reads after removing adapter dimers, junk reads and repeats. Subsequently, length screening was performed to preserve the unique sequence with base length between 18 and 25 nucleotides in miRbase 22.0 (www.mirbase.org). The mapping procedure was performed on pre-miRNA against A. annua and N. benthamiana genome data, respectively. Sequences were compared with mRNA, Repbase, and RFAM databases to filter to obtain valid data, and the precursor and genome were compared for miRNA identification. To identify the results of putative miRNAs in A. annua and N. benthamiana, all the obtained miRNAs were used to predict the secondary structures using RNAfold software (http://rna.tbi.univie.ac.at/cgi-bin/RNAWebSuite/RNAfold.cgi).

After normalizing deep-sequencing counts, differentially expressed miRNAs were analyzed using Fisher's exact test. The significance threshold was set to be 0.05 in the test.

To predict the targeted genes by differentially expressed miRNAs, PsRobot^[38] (http://omicslab.genetics.ac.cn/psRobot/) was used to make sure the binding sites of miRNA. The software predicts targets using a plant-based target penalty strategy (the default threshold is Score ≤ 2.5). Meanwhile, GO terms^{[39, 40]} and KEGG pathway^[41] of miRNA were also calculated and annotated. According to screened the most different expression miRNAs and its target gene, a correlation network diagram was manufactured through OmicStudio tools at www.omicstudio.cn/tool.

Transcriptome analysis of hormone transduction genes, phytohormones extraction and quantitative analysis

In order to learn more about the hormonal changes in the plant during grafting, we firstly collected seven hormones for statistical analysis, including auxin, cytokinin, gibberellin, ABA, ethylene, JA and SA. Genes associated with phytohormone signal transduction were extracted owing to the pathways involved. The changes in the amount of expression during each transduction process were determined by the genes in that module.

The extraction protocol of phytohormones was followed using the simultaneous targeted profiling of Simura et al.^[42] Five plant tissues were contained in the quantitative experiment, including Aa leaf, Aa stem, Nb leaf, Nb stem, and Nb root with non-grafted samples as control. After pretreating of solid phase extraction, samples were analyzed by HPLC-MS with Aglient system in positive and negative ion models, respectively. Besides, the mobile phase includes acetonitrile and water solution containing 0.05% formic acid. Each tissue sample contains three biological repeats. Abscisic acid (ABA, neg263), jasmonic acid-isoleucine (JA-L-ILE, neg322.1), trans-Zeatin-riboside (tZR, pos352.1), and N6-isopentenyladenine riboside (IPR, pos336.1) were chosen for representative indicators with warfarin (307.1, 308.9) as the internal standard. MassHunter software was used for identifying the type of phytohormone and calculating relative content.

We firstly constructed a grafting model involved with A. annua as scion and N. benthamiana as rootstock to investigate heterograft changes at molecular and metabolic levels (Fig 1, Supplemental Fig. S1).

After observing the morphology of Aa and Nb in both heterograft and non-graft plants, we found no significant difference in phenotype, which may be the cause of short growth time. And then, we observed the trichomes by OLYMPUS microscope (Fig. 1b). Through the contrast the density of leaves, we found after grafting, the newborn leaves showed a marked state of stunting, and decreased the number of trichomes, which was about half with before. Through quantitative measurement of artemisinin, we found that the content showed a significant drop in Aa leaves as expected. However, the content in Aa stem was slightly increased (Fig. 1c).

Gene differential expression and mRNA migration

In order to distinguish the differentially expressed genes during graft process, we analyzed transcriptome database of scion and rootstock part with non-grafted samples, respectively. Meanwhile, as the development of graft union, genes in the plant infiltrated into each other, thus resulting in changing of gene regulation and physiological activity. Therefore, we used the measured transcriptome database of scion and rootstock to compare with non-graft Nb and Aa genome to search potential mobile mRNAs.

In the graft process, a total of 7,794 DEGs (Differentially Expressed Genes) of A. annua scion were detected. Among those DEGs, 4,754 genes were up-regulated and 3,040 genes were down-regulated (Fig. 2). Through Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of DEGs (Supplemental Fig. S2), mostly GO-enriched genes were distributed in 'DNA-templated', 'protein phosphorylation', 'defense response' and 'signal transduction'. Besides, KEGG pathways were significantly enriched in plant-pathogen interaction, plant hormone signal transduction and metabolism process.

Figure 2. 

DEGs and potential mobile mRNAs identified from the Aa/Nb heterograft model. The upper scion part belonged to A. annua, and the bottom rootstock part was N. benthamiana as separated with horizontal line. The left half of the information was the number of up-regulated and down-regulated of different expression genes related to comparing with the same plant. The right half of the photo exhibited the movement process of mRNAs and identification of non-homologous genes.

To identify root-to-shoot mobile mRNAs from N. benthamiana grown in A. annua, measured transcriptome were compared with tabacco genome again, 50 genes of N. benthamiana were identified as mobile mRNA. GO analysis (Supplemental Fig. S2) of these mobile mRNAs reveal that the most overrepresented proportion is molecular function (42.6%), which contained binding of ATP, GTP and calcium ion, 'structural constituent of ribosome', 'calcium ion binding', activity involved in ATPase, catalytic, and hydrolase accounted as the major part. Within the biological process (36.1%) category, 'cell wall macromolecule catabolic process', 'translational elongation', and 'oxidation-reduction process' occupied the main position. Moreover, intracellular, ribosome and membrane were overrepresented in the cellular component (21.3%) category. Meanwhile, the most significant pathways (Table 1), in which movement genes are primarily located, were ko03010 (ribosome), ko04141 (protein processing in endoplasmic reticulum), ko03040 (spliceosome) and ko01200 (carbon metabolism). It is worth nothing that there was an ABC transport gene named Niben101Scf01719g08010 obtained in mobile mRNA. Owing to the important role in the transport and accumulation of secondary metabolites in plants, we could transfer that the mobile Niben101Scf01719g08010 may also be involved in the diterpenes metabolic process in A. annua.

In addition, a total of 8,214 DEGs of N. benthamiana rootstock were identified after grafting, in which 4,980 genes were up-regulated, whereas 3,234 DEGs were down-regulated. GO and KEGG analysis of the DEGs (Supplemental Fig. S3) revealed that the gene type and pathway involved in rootstock were both significantly different with scion, indicating that within the metabolic flux of the two exists a wide discrepancy. The top3 GO enrichment terms were oxidation-reduction process, oxidoreductase activity and regulation of transcription, DNA-templated. Furthermore, substance metabolism took up the majority in KEGG pathways, which includes metabolism of terpenoids and polyketides, metabolism of cofactors and vitamins, and lipid metabolism. During the formation process of Aa/Nb, we found 20 mobile genes of A. annua in the rootstock. Through GO analysis (Supplemental Fig. S3), 39.5% of mobile genes located in molecular function, such as 'protein binding', 'fatty-acyl-CoA reductase (alcohol-forming) activity' and 'transporter activity'. In the category of biological process, above 35.8% genes have the function, in which 'oxidation-reduction process' and 'response to cadmium ion' were enriched. The most remarkable terms of cellular component (24.7%) were cytosol, nucleus and cytoplasm. The mobile Aa transcripts may be involved in a variety of biosynthesis and substance metabolism (Table 2) in N. benthamiana, photosynthesis (ko00195), cutin, suberine and wax biosynthesis (ko00073) and oxidative phosphorylation (ko00190) for instance.

Identification and regulatory network research of miRNA

miRNA sequencing, screening and identification

We constructed four miRNA libraries named Aa, Aa scion, Nb, Nb rootstock from grafted union. As shown in Table 3, the raw reads obtained from four database were 14,833,433 14,976,061, 14,057,479, 13,596,666, respectively. After removing 3' adaptor and junk reads, following screening and retaining reads with base length from 18 to 25 bp and RNA database comparison, a total of 10,043,761, 9,486,686, 8,143,546, 5,473,716 valid reads were obtained.

The overall distribution of differentially expressed miRNA was analyzed using a volcano figure (Fig. 3). Most miRNAs were not differentially expressed during grafting. Additionally, in grafted union, the total number of differentially expressed miRNAs in Aa scion is more than Nb rootstock. Besides, the number of up-regulated miRNAs were both higher down-regulated ones, whether for scion or rootstock. The top three most abundant miRNAs in scion and rootstock are miR159, miR396, miR166 and miR159, miR6149, miR166, respectively.

Figure 3. 

Volcano photo of differentially expressed miRNAs in scion and rootstock. Overall distribution of differentially expressed miRNA in (a) Aa scion and (b) Nb rootstock. The abscissa represents the differential expression multiple changes of miRNA in different samples. The ordinate represents the statistical significance of miRNA expression changes. Red dots represent significantly up-regulated differentially expressed genes, blue dots represent significantly down-regulated differentially expressed genes, and gray dots represent non-significantly differentially expressed genes.

Regulatory network research of transcription factors in differentially expression miRNAs

To determine types of differently expressed miRNAs, we chose p values < 0.01 and higher expression level miRNAs to analysis its potential function. Through psRobot^[38] software, we can predict the target genes of miRNAs with significant differences.

As expected, most of the difference miRNAs were located in plant-pathogen interaction (ko04626) and monoterpenoid biosynthesis (ko00902) of Aa scion, aminoacyl-tRNA biosynthesis (ko00970), plant hormone signal transduction (ko04075) and pentose and glucuronate interconversions (ko00040) of Nb rootstock, respectively. Thus, we can infer that most miRNAs we obtained may be involved in the recovery of grafted union physiological function and transduction of information.

Considering the crucial physiological function of Transcription Factors (TF) in plant development^[43], we focused the miRNA-TF network as a key point to analyze the regulation mode. According to the significant expression difference, 142 miRNAs predicted to be transcription factors of Aa scion were obtained firstly by screening. Among the miRNAs, seven types of TF were classified as shown in Fig. 4a, which contain R2R3-MYB, bHLH, GRAS, GAMYB, SBP-box, MASD-box, IIS.

Figure 4. 

Regulatory network research of transcription factors in differentially expression miRNAs and heatmap analysis of target gene among Aa scion and Nb rootstock. (a) Network plot of different expression miRNAs in A. annua scion. (b) Heatmap of different expression miRNAs targets in A. annua scion. (c) Network plot of different expression miRNAs in N. benthamiana rootstock. (d) Heatmap of different expression miRNAs targets in N. benthamiana rootstock. The legend at the bottom of the left figure marks the miRNA name, target gene name and gene type in different colors and shapes, respectively. Heatmap of target genes was used FPKM as expression quantity. Different colors indicate different gene expression levels, from blue to white to red, indicating low to high expression levels, with red indicating high expression genes and blue indicating low expression genes. The group of each gene corresponds to the color bar on the left.

In the Aa scion, the most difference expression miRNAs belong to hormone transduction and plant resistance. We combined the expression quantity of target genes to analyze the regulatory trends. Fifteen genes related to miRNAs were obtained for constructing the heatmap and mostly of genes demonstrating the up-regulated after grafting, which is the reason why miRNA promoted its expression.

Additionally, the expression level of MYC2 gene, named AA518540 and AA477190, was always very high throughout, indicating the jasmonate signaling factor could play an important role in the graft process or development of grafted plants. Based on this, we can primarily speculate that the 15 target genes are essential for the normal growth and development of A. annua, which is worth further research.

By contrasting the miRNA database of Nb rootstock with non-grafted Nb, we preliminarily obtained 52 different expression miRNAs related to transcription factor, which actually belonged to only one type miRNA, named miR169 (Fig. 4). Through the heatmap of target genes, we can clearly find that the expression level of each gene was increased, which may be related to the cleavage effect of miR169. During grafting, miR169 and its target gene NF-Y in N. benthamiana expressed the most outstanding role, which told us the prominent impact of miR169 in the grafting model.

Hormone signal transduction during graft process

Phytohormones, found in auxin, cytokinin, gibberellin, and abscisic acid, play a key role in graft union and act as signal molecules related to graft development. Considering the different effects related to graft union of every plant hormone, we chose seven kinds of hormones to investigate its changing mechanism as shown in Fig. 5.

Figure 5. 

Changes of genes involved in phytohormone signal transduction. (a) Phytohormone transduction pathway of A. annua scion. (b) Phytohormone transduction pathway of N. benthamiana rootstock. Changes in DEGs were mapped in boxes, green and red boxes represent down-regulated and up-regulated, respectively. Orange boxes represents genes in the pathway indicting bidirectional regulation. Purple boxes illustrate no DEG.

Auxin transduction

In the grafted model we built, DEGs of auxin transduction pathway basically showed a consistent change trend of expression including A. annua scion and N. benthamiana rootstock except TIR1 and GH3. In the heatmap, we can infer that genes of A. annua in GH3 were subject to positive regulation, however, some of GH3 genes in rootstock display down-regulated, this may be related to negative feedback regulation. In the other two primary-response genes, AUX/IAA and SAUR both show bidirectional regulation. AUX1 is an important auxin influx carrier, which mainly regulates root hair development and root gravitropism. In the process of grafting, scion and rootstock reveal the down-regulation of AUX1 expression. In conclusion, ARF family consists of up-regulated and down-regulated genes, which may be concerned with their role in the metabolism pathway.

Cytokinin and abscisic acid transduction

CRE1 is a cytokinin receptor, it can be activated by CKs to initiate the phosphorylation signal. The down-regulation of CRE1 of scion and rootstock may result in vascular bundle cells differentiate into xylem cells. AHP is a histidine-containing phosphate transfer protein, which may involve in drought stress or cold signal regulation through redundantly negative manner. In scion part, AHP is up-regulated, whereas down-regulated in rootstock. During the graft process, B-ARRs reveal identical expression model in scion and rootstock, in which members may exhibit different transcription function, whereas A-ARR is down-regulated in scion and up-regulated in rootstock. On the one hand, it may be affected by the regulatory effect of B-ARR on it. On the other hand, it may be closely related to the growth condition and development of grafted plants.

In the Aa/Nb grafted model, genes involved in ABA transduction reveal a clear upward trend, especially PP2C and ABF in scion, PYR/PYL and ABF in rootstock. PP2C, actually a kind of serine-threonine phosphatases type 2C protein, act as a negative regulator of ABA transduction.

Gibberellin and ethylene transduction

In the scion part, the expression level of GID1 and GID2 was significantly up-regulated. However, the expression trend of these two pathways was not obvious in rootstock. Besides, GID2 could repress DELLA, which may influence the expression of downstream gibberellin signal during grafting process.

The ETR gene family as an ethylene receptor, can play a role in seed germination and be induced by physiological process such as plant senescence through Ca²⁺ and ABA signal^[44,45]. We have found surprising discovered that the controlling model of ETR is reverse in graft union, representing the scion and rootstock reveal enormous differences in ethylene signal transduction. Moreover, ethylene insensitive 3 (EIN3) is a crucial transcription factor in ethylene signal transduction and biosynthesis. According to our measured data, EIN3 reveals identical expression models of up-regulation, which may be consistent with its positive effect.

JA and SA transduction

In scion tissue, MYC2 is up-regulated, which indicates that the intensity of JA signal may improve during grafting. Additionally, JAZ is up-regulated in both scion and rootstock. The JAZ gene family is a series of Jasmonate-zim domain protein, which can act as transcription repressors and JA co-receptors. The expression model of JAZ may be closely related to its various effects in physiological function.

SA signal transduction did not show a strong trend of change. In the scion part, the three pathways in SA contain up-regulated genes and down-regulated genes, whereas in rootstock, NPR1, the critical role in SA transduction, demonstrate down-regulated trends leading to overexpression of TGA. Besides, NPR1 could perform antagonistic effect of JA and SA. Down-regulated of NPR1 may result in enhancement of JA signal. The up-regulated of TGA can lead to a positive response to SA signals, and then respond to plant pathology as soon as possible.

Quantitative analysis of endogenous hormones

According to our focus phytohormones and testing instrument, we finally obtained three types of hormones to analysis its regulation process in scion and rootstock, including ABA, JA and CK in Fig. 6. In order to accurately measure changes in each plant part, five tissues were obtained for quantitative determination, including two scion parts and three rootstock parts.

Figure 6. 

The content of measured ABA, JA and CK in A. annua scion and N. benthamiana rootstock. The abscissa represents the sampling site, and the ordinate represents the phytohormone concentration. The error bar represents the standard error for three independent experiments.

N6-isopentenyladenine riboside (IPR) and trans-zeatin-riboside (tZR) are two crucial CKs, which can regulate and control plant growth and differentiation although the amount is low. In the background of grafting, we can clearly find that the content of IPR and tZR is both enriched in stem, which may be related to its physiological effect. Moreover, the trend of content change is identical in IPR and tZR. In the scion part, the purity of IPR and tZR is elevated among A. annua leaf and stem, whereas in rootstock, IPR and tZR show a downward trend.

After the graft stage, the content of ABA is mainly down-regulated except N. benthamiana stem, which is clearly higher than before. Besides, the relative content of ABA is highest in N. benthamiana leaf and stem. JA plays an important role in defending plant resistance, plant pathology and pest disaster, thus it is clearly show that the content of JA derivative, named JA-LIE, represent substantial increases in each part after grafting.

Heterografting always have identical effects on finding coding RNA through different constructed systems. However, there are individual differences among each species. Before the systematic understanding of mechanisms involved in the grafting process, the changing patterns of individual species should be analyzed. Although there was no various epigenetic characteristic change, endogenous physiology and components may already be occurring and changing. We speculated that the artemisinin may create a mobile behavior forwardly or passively from scion to rootstock during grafting. This phenomenon showed that the heterografted plant tends to grow less than normal, and need more nutrients to keep its physiological activity.

Nicotiana benthamiana is common model plant used in genetic modification, subcellular localization and grafting. β-1,4-glucanases^[29] produced by N. benthamiana could facilitate reconstruction of cell walls in graft union, thus using N. benthamiana as rootstock may greatly improve graft compatibility. ATP binding cassette (ABC) was an enormous family and widely existed in plants, like tobacco, Arabidopsis, rice^[46−48]. It can play an important role in response biotic stress, especially ABCG subfamily. In tobacco, NtPDR1, the members of ABCG subfamily, could transport diterpenes, thus promoting the ability of plant defense. In our grafted model, the expression level showed no obvious changes, which may be related to the existing good grafting conditions for both.

Artemisia annua is an important traditional Chinese medicine, and it is famous for its effect of healing malaria. The content of artemisinin has been the focus of research. Through the grafted model, we preliminary understood the genes involved in resistance of A. annua. By using the genes related to plant growth and stress resistance, we may find an effective way to research biological properties in heterografted plants.

After 30 DAG, most of the grafted samples could still survive, with some deaths due to the influence of the external environment or original sample situation. However, it is not possible to determine whether A. annua and N. benthamiana exhibit graft compatibility due to the short co-growth time. Through the analysis of the transcriptome of scion and rootstock, we can obtain a large amount of information. Most of DEGs belong to plant resistance and growth, which is consistent with what we suspected. We found mobile mRNAs moved from rootstock to scion or scion to rootstock. These genes may play a key role in the grafting process. Among the mobile genes in A. annua scion, AA476810 caught our attention owing to the function of cutin and wax biosynthesis. In A. annua, the biosynthesis of cutin and wax could directly influence the development of secretory trichome in leaf, thus affecting the content of artemisinin. It is not sure if the movement of AA476810 caused the lower artemisinin level in grafted scion. In addition, it is not known whether the movement of this gene will result in impaired growth of new A. annua leaves or affect the growth of tobacco leaves. In 50 mobile N. benthamiana genes, an ABC transporter named Niben101Scf01719g08010 was brought forward for research. The ABC transporter family could participate in various physiological processes, especially the transport of terpenoids^[49]. So, we guess that the mobile of ABC transporter may influence the biosynthesis of artemisinin derivative. The transport function of the N. benthamiana gene was validated by constructing into the PDR196 vector, whereas lacking good repeatability, we failed to obtain clear results.

The grafting process may be closely related to the activity of miRNA

Until a few years ago, owing to technical limitations and inadequate knowledge of miRNA, the important situation of small RNA has not received enough attention. However, with the successful blowout of related research, the importance of miRNA is becoming increasingly known. In plants, miRNA always influences morphological development of leaves, and could participate in plant hormone signal transduction and cope with environmental stress. We analyzed the different expression of miRNAs in scion and rootstock, finding that the number of up-regulated miRNAs is more than the down-regulated. We then chose the transcription factor belonging to high expression abundance and significant differences in miRNAs for miRNA-gene-GO analysis. In A. annua scion, the screened TFs can be separated into three classes, including plant development, phytohormone transduction and abiotic stress responses. The 15 target genes may play an important role in the growth and development of A. annua. However, only one type of TF was present in N. benthamiana rootstock, NF-Y, which is the key point for miR169 binding. The class of NF-Y is able to take part in a variety of plant physiological activities, such as hypocotyl elongation, flowering process and responding to abiotic stress, thus we speculated that miR169 and NF-Y played an indispensable role in the grafting process.

Firstly, miR159, miR396, miR166 and miR159, miR6149, miR166 were screened owing to the high abundance in Aa scion and Nb rootstock, respectively. The different expression miRNAs like miR159 and miR166 have been reported in other grafted models, which indicated the conservative and importance of these miRNAs. It is well known that miRNAs could participate in various plant physical activity and responding to environmental change^[18]. Owing to its short nucleotide sequence, miRNA can transfer from root to shoot to exercise the role, vice versa^[16,17]. In Aa/Nb grafted model, miR159 and miR166 are the same high expression abundance. It has been marked that the two microRNA are vital network hubs in plant stress response, indicating it may participate in callus formation and plant physiological activity. Besides, miR159 and miR166 also have the function of plant fertility^[50,51] and parthenocarpy, respectively. Additionally, miR396^[52] and miR6149^[53] are both ancient non-coding small RNAs, which can play crucial roles in various activities and stress response.

To date, miR169 has been regarded as a ubiquitous regulator to responding various abiotic stress by obtaining new promoter^[54]. And the target gene of miR169 in plants is nuclear transcription factor Y. Nuclear transcription factor Y is widely present in eukaryotes and can specifically bind CCAAT-box^[55]. As a conserved transcription factor, NF-Y plays an important role in plant growth, tissue development and responding to stress. It has been reported that 3'UTR of NF-YA (subunit A) has a cleavage target for miR169, thus miR169 could regulate the expression level of NF-Y via transcript cleavage. The tobacco rootstock not only inputs nutrients from roots into scion, but endures grafting and external stress during the graft process.

As is well known, MYB TF is one of the largest plant families and can clearly influence the trichome formation and content of artemisinin^[56]. Additionally, GAMYB could encode R2R3MYB domain TF in plants, which is the target gene of miR159. Owing to the conservatism and universality of miR159, GAMYB plays a key role in inhibition of plant development. Besides, MADS-box, MYC2 and GRAS TFs are key nodes of ABA, JA and SA, JA, respectively. Thus, we can speculate that due to the different abundance of miRNAs, the expression of target genes may show huge changes, which further influence hormone levels and growth of Aa. Additionally, SBP-box and IIS transcription factor mainly play roles in defending abiotic stress and improve tolerance in adversity.

Combined with existing reports and the results of this experiment, we put forward that the six miRNAs can act as potential biomarkers in successful heterografted plants. We conferred that there also existed mobile miRNAs during the grafting process, however, because of the short length and conservative between different samples of miRNAs, we cannot clearly identify the movement trajectory. The mobile miRNAs may perform an important function in grafted plants.

Hormonal responses in grafted plants

Hormones are also important for proper growth and development in plants, especially in promoting blood vessel growth and responding stress signal. We analyzed seven types of phytohormone transduction pathway and three types of hormone levels to realize its changing tendency during the grafting process.

As a hormone that can regulate a variety of life activities, auxin could canalize the pathway to rebuild vascular tissue in graft union^[11]. Additionally, it can regulate expansion of cortex cells and cell proliferation of vascular tissue during graft tissue remodeling, including differentiation of xylem and phloem. TIR1 possesses a F-box protein, part of formation of the SCF^TIR1 ubiquitin-ligase complex^[57], which interacts with Aux/IAA transcriptional repressor and mediates its degradation in the presence of auxin, thus controlling the initiation of adventitious root. In A. annua scion, the expression of TIR1 in the whole pathway showed an obvious up-regulation trend, whereas there is no evidence of significant genetic differences in tobacco during grafting. After auxin treatment, AUX/IAA, most of the GH3 genes and SAUR are all immediately expressed, this is the reason why these three were considered as primary-response genes^[58]. The effect of GH3 family genes on plant growth and development is visible in plants. For instance, the mutation of GH3.9 in Arabidopsis thaliana leads to taproot length increase, suggesting members of the GH3 family may be involved in root development. Furthermore, GH3 can also exhibit the ability to alter leaf phenotypes and response to biotic or abiotic stress. ARF, short for auxin response factors, are considered as a new family formed during the evolution of land plants^[59]. ARF genes are included in the control of many physiological effects, such as cell division and leaf morphology. ARF act as both transcriptional activators and repressors.

It has been noticed that CKs (cytokinin) could regulate shoot regeneration, plant morphogenesis, biotic and abiotic stress. In addition, the role of CKs during the graft process is the focus point. Through regulatory networks of CKs, it can influence cell division and promote graft interface healing, thus leading to successful grafting^[60]. A-ARR and B-ARR are mainly response regulators of cytokinin transduction. Type A-ARRs member, play a crucial role in negatively regulating the transduction of cytokinin, thus influencing multiple plant phenotypes and developmental process^[61]. Type B-ARRs are positive regulators of cytokinin response, it can control the expression of downstream genes. Besides, B-ARR could also influence shoot development by constructing a hormone transcriptional network.

According to our transcriptome data, the expression pattern of genes on the ABA transduction pathway appears to be an up-regulated stream, which is consistent with ABA physiological function. However, phytohormones do not always exhibit their role alone, especially ABA and auxin CKs, thus gene expression should be studied from a more comprehensive and systematic perspective^[58]. As is well known, ABA (abscisic acid) refers to a kind of plant hormone that causes bud dormancy, leaf shedding and inhibit cell growth. In addition, many researchers have demonstrated that ABA, auxin, CKs have antagonistic action^[42]. PYR/PYLs is a family of START proteins, which can directly influence pyrabactin and ABA signal. Besides, PYR/PYL act as receptors of ABA, which can in turn influence the expression of PP2C. SnRK2 (SNF1-related protein kinase2), mediates positive effects of ABA signal transduction, and then regulates root growth and seed germination^[62]. ABFs are a series of transcription factors that could act in coordination on the ABA signal pathway.

Gibberellin and ethylene have been identified to collaboratively trigger a variety of plant development processes, such as response to abiotic stress and regulation of ovule death^[63]. GID1, short by GA receptor protein, and GID2 (a F-box protein), were successively proved to participate in GA signal transduction, which protein complex bind to the key factor DELLA protein.

Through the synergistic or antagonistic action of JA and SA, many plant pathogen diseases could be defended. Besides, JA and SA also play an important role in the synthesis of metabolites respectively and adaptation to biotic and abiotic stress. MYC2, a bHLH transcription factor, has been recognized as the most essential role in JA signal transduction. During plant development and defence against pathogens, MYC2 could integrate multifaceted hormone signals to balance physical activity.

However, in the analysis of transcriptome, owing to the positive and negative feedback adjustment among plants, genes in each signaling pathway did not show a consistent trend. At the level of gene expression, we can only focus on the key points reported, such as DELLA in gibberellin, MYC2 in JA. Additionally, through the pathway of hormone transduction, it is clearly that the regulation process of hormones is in a dynamic change law, which is closely linked to plant physiological state and changes in the external environment. The change trend of different phytohormones is not well understood, it is often coordinated by a variety of physiological and environmental factors. However, through the quantitative experiment we built, we can make a preliminary judgment that the changes of hormones in scion part and rootstock part were identical, and then the up-regulated content in scion did not show the same trend in rootstock. Additionally, phytohormone is always in a dynamic process in terms of content, the sampling of 30DAG just represents the hormonal equilibrium of grafted plant after a period of symbiosis.

The content of endogenous phytohormone is always low, to clearly understand the changing pattern among grafted plants, five parts were divided. As shown in the bar chart of hormone content, hormone levels vary in each part of the plant (Fig. 6). ABA could promote leaf abscission and inhibit cell elongation. In Aa scion, Nb leaf and root, there was no significant gap in control group and treated group. However, it is clear that in Nb stem, the content of ABA is substantially increased after grafting, which may be related to the fact that plants perceive changes in the outside world as a factor that inhibits their own growth as rootstocks.

Similarly, the content of iPR and tZR in Aa stem and Nb stem also showed huge dynamic change. Surprisingly, the content of JA was a substantial improvement in five parts, indicating JA may be involved in various physiological development among whole plants.

Heterografting is an important method to understand regulatory mechanisms of plant stress response, including long-distance transport small RNA, coding RNA, phytohormones and other metabolites. In the present study, we aimed to find out more new potential mechanisms during the grafting process of A. annua and N. benthamiana. We constructed an Aa/Nb grafting model to analyze its change in mRNA, miRNA and phytohormone level. A month of co-growth will help us understand quickly the molecular and metabolic changes occurring in plants. 7,794 DEGs (different expression genes) and 8,214 DEGs were identified in Aa scion and Nb rootstock, respectively, which mainly belong to defense response and signal transduction in scion and substance metabolism in rootstock. During the grafting process, 50 Nb genes and 20 Aa genes were identified as potentially active genes. miR159 and miR166 were considered as biomarkers of successful grafting plants owing to its conservation and physiology. Besides, R2R3-MYB, bHLH, GRAS, GAMYB, SBP-box, MADS-box, IIS in scion and NF-Y were regarded as key genes involved in growth and development of grafted plants. The target genes screened may be important for improving plant stress resistance and regulating metabolite production in vivo, which could be used for directional breeding.