Phenotypic difference between vegetative hybrids and self-rooted grapes
The 20 vegetative hybrids of Vitis vinifera / Schisandra chinensis were grown normally for 25 years (Fig. 1A). The leaves of Vs hybrids above the graft junction were significantly different from those of self-rooted of V. vinifera ‘Hongruibao’ (Fig. 1B and 2C, top right): the leaf apex of the hybrids was more obtuse, and the reverse side was covered with shorter and whiter hairs (Fig. 1B and 2C, top left), while below the graft union, the Vs hybrids produced five dehiscent leaves, and the veins were radial rather than reticulate, and every vein was clearly visible (Fig. 1B and 2C, centre down). The highly significant characteristic of Vs hybrids was the reverse side of the leaves above the graft union which were white colour (Fig. 1C, top left), but the leaves of Vv plants and Vs hybrids below the graft union were deep green (Fig. 1C).
The internode of self-rooted V. vinifera ‘Hongruibao’ was the longest (Fig. 1D) at 14.7 cm, and was covered by a lot of hairs. The internode of Vs hybrids above the graft union was of medium length at 7.4 cm or approximately half that of V. vinifera ‘Hongruibao’, with a hairless surface (Fig. 1E). The length of the Vs hybrids internode below the graft union was the shortest at 5.3 cm which was only approximately one third as long as that V. vinifera ‘Hongruibao’ (Fig. 1F). Vv, VS-above and Vs-below showed significant differences in internode length (Fig. 1M). In total, the internode length of vegetative hybrids showed a tendency that was between rootstocks and scions.
The fruits of self-rooted V. vinifera ‘Hongruibao’ were oval and lighter green (Fig. 1G), but the fruits of Vs hybrids were round and more purple-redder than those V. vinifera ‘Hongruibao’ (Fig. 1H) with a sourer teste. Both types of fruits were covered with thick wax, but Vs hybrids was thicker than self-rooted. The difference in terms of the length (30.07 ± 1.14 mm) and weight (10.46 ± 1.15 g) of mature fruits of Vv plants differed significantly (P ≤ 0.01) from those of mature fruits of Vs hybrids (23.06 ± 1.21 mm and 6.86 ± 0.82 g, respectively), but there were no significant differences in the widths of the fruits (Fig. 1N).
The graft union was swollen (Fig. 1I), and the anatomical structure indicated the formation of numerous bridges for the vascular bundle (Fig. 1J), but those of the self-rooted V. vinifera ‘Hongruibao’ were not obvious (Fig. 1K). At 200 μm scale, the numberous brides of vascular bundle from hybrids Vs was 3 paris, but self-rooted was only 1 pari (Fig. 1I-K). In total, the area above and below the graft union was balanced, indicating that both growths were concordant and compatible.
We reproduced the Vs hybrids using twig cuttings above the graft junction in a garden in June 2018, and found that the traits of the seedlings, including leaf type, colour and hairs on the reverse side of leaves. The radial veins were stable with the adult tree (Fig. 1L) and the result showed the new traits were stable inheritance. However, the rate of survival in these Vs hybrids was 15%–20% which was much lower than that V. vinifera ‘Hongruibao’ (80%–85%), suggesting that the unclear barrier materials was produced in graft hybrids.
Transcriptome sequencing and functional annotation of DEGs
Since there are many morphological differences between heterografted V. Vinifera/S. chinensis (Vs) and self-grafted grapes (Vv), we hypothesized that heterografted may resulted in significant changes in genetic information and performed transcriptome sequencing on the scions of Vv and Vs.
A total of 43.68 GB of sequencing data contained 291,145,166 raw reads. Clean data were obtained by filtering out unqualified sequences using Cutadapt [32], after which it still contained 288,679,444 reads. Hisat [22] was used for the reference genome comparison of clean data and produced 30,609 genes from 131,224,950 unique mapped reads. The summary of the Illumina Hiseq4000 transcriptome sequencing of the two groups is shown in Table 1. The screening criteria for differentially expressed genes were |log2fc| > 1 and P < 0.005. There were 222 DEGs; the up- and down- regulated DEGs are shown in Table S1. The GO annotation of each type was sequenced from high to low, according to the number of target genes annotated, and are presented in Table S1. The five most abundant transcripts involved with biological processe being those associated with the oxidation-reduction process, metabolic process, protein phosphorylation, regulation of transcription and transport. Among the 15 cellar components, the three most highly represented were the membrane, and integral components of the membrane and nucleus. Of the 10 molecular function categories, the two most significant categories were ATP binding and oxidoreductase activity.
To better understand the biological roles of these genes, we analysed them with KEGG annotation. All DEGs and their annotations are listed in Table S1. KEGG enrichment analysis results showed that the degradome target significant enrichment KEGG pathways were mainly associated with carbon metabolism, ribosomes, biosynthesis of amino acids, protein processing in the endoplasmic reticulum, glyoxylate and dicarboxylate metabolism, photosynthesis, and glycolysis and gluconeogenesis.
We took the materials of leaves, fruits, and phloem of the scions for RT-qPCR experiments to verify the sequencing results. The RT-qPCR results were basically consistent with the sequencing results. While some genes were up-regulated or down-regulated in all three parts, while some genes were up-regulated or down-regulated in only one or two organs (Fig. 2A).
sRNA-seq and identification of miRNAs
The expression of many genes is regulated by miRNAs, so we conducted miRNA sequencing to search for miRNAs that may play key functions in heterologous grafting. The small RNA library was assembled from mixed miRNA pools following the same procedures used for transcriptome sequencing. Identification of V. vinifera miRNAs and their targets using high-throughput sequencing and degradome analysis were conducted using methods previously reported [33-36]. Comparison of results of the Rfam database [37], The length distribution of unique sRNA is summarized in Fig. 3, and most sRNAs are 24 nt in length.
The unique sequences mapped to species-specific mature miRNAs in hairpin arms were identified as known miRNAs. The unique sequences mapping to the other arm of known species-specific precursor hairpins opposite to the annotated mature miRNA-containing arm were considered to be novel 5p- or 3p-derived miRNA candidates. The remaining sequences were mapped to other selected species precursors (with the exclusion of specific species) in miRBase 21.0 using a BLAST search. The mapped pre-miRNAs were further BLASTed against the specific species genomes to determine their genomic locations. The sequences could each be mapped to species-specific miRNAs we defined as known miRNAs. Unmapped sequences were BLASTed against the specific genomes, and RNAs containing hairpin structures were predicted from the flank 120 nt sequences using RNAfold software (http://rna.tbi.univie.ac. at/cgi-bin/RNAfold.cgi). Those that met the criteria for secondary structure prediction were potentially novel miRNA candidates. In total, 1,195 miRNAs were found in the RNA library and categorised into six groups based on their match situations in miRbase. Group 1a (reads were mapped to miRNAs/pre-miRNAs of V. vinifera in miRbase and the pre-miRNAs were further mapped to the genome and EST) and group 4 (reads were not mapped to pre-miRNAs of selected species in miRbase, but were mapped to the genome and extended genome sequences from genome that may form hairpins) displayed the highest redundancy. RT-qPCR was used to verify the sequencing results, and the expression trend of most miRNAs was consistent with the sequencing results (Fig. 2B).
Differential expression analysis and function annotation of miRNAs
To identify the effects of miRNAs in vegetative hybrids, the differential expression of miRNAs in six libraries of two groups was compared using the read counts generated from the high-throughput sequencing (Table S2). In total, 79 miRNAs (P ≤ 0.05) containing 27 known miRNAs and 52 new miRNAs showed differential expression patterns. Among these, 43 miRNAs were up-regulated in Vs hybirds as compared to Vv plants. Furthermore, there were 42 differentially expressed miRNAs with mean expression abundance > 10, including 22 up-regulated and 20 down-regulated miRNAs. When the screening condition was relaxed to P ≤ 0.1, there were 83 up-regulated miRNAs and 76 down-regulated miRNAs. All detected miRNAs belonged to 46 gene families. Most miRNAs of the same gene family had the same expression pattern. For example, the six miRNAs in the miR162 family (ptc-MIR162b-p5, vvi-MIR162-p5, vvi-miR162, tcc-miR162_R+1, hpe-miR162a_R+2, and stu-miR162a-5p_R+3_1ss4GT) were down-regulated in Vs hybrids compared to Vv plants, and the five members in the miR530 family (ptc-MIR530a-p3_2ss12CA21GA, tcc-miR530b_L+2R-1_1, tcc-miR530b_R+1_1ss9CT, and tcc-miR530b_L+2R-1_2, stu-MIR530-p5) were up-regulated in Vs hybrids compared to Vv plants. We also found some variable expression patterns between miRNAs in same family. There were 20 miRNAs in the miR339 family; nine had low expression levels, the other nine (vvi-MIR399a-p5, vvi-miR399a, vvi-miR399e, vvi-MIR399h-p5, vvi-MIR399d-p5_1ss20TC, vvi-miR399d_L+2R-2, vvi-MIR399i-p5, vvi-miR399i, and tcc-miR399a_1ss21GT) were down-regulated in Vs hybrids compared to Vv plants, while two (vvi-MIR399g-p5 and vvi-miR399) were up-regulated. The cluster analysis data are shown in Table S3.
psRNATarget has been used to predict the target genes with significantly different expressed miRNAs [38]. Results of differential miRNA target gene prediction showed the target gene information corresponding to miRNA and provided the GO and KEGG annotation information of the target genes (data shown in Tables S4 and S5). GO Enrichment showed that target genes were mainly classified as being associated with the defence response, biosynthetic process, salicylic acid biosynthetic process, signal transducer activity, and drug transmembrane transporter activity.
Target prediction and category statistics of miRNAs by degradome sequencing and targets functional enrichment
Because miRNAs mainly perform biological functions by degrading target genes, we performed degradome sequencing. Through degradome sequencing, 22,945,508 raw reads containing 5,922,957 unique reads were obtained from the mixed degradome pools. There were 14,163,908 (61.73% of all reads) transcript mapped reads containing 4,607,411 unique transcript mapped reads. The number of input transcripts was 52,365 and 45,728 transcripts were covered. All potential cleaved transcripts were divided into five categories according to the signature abundance at each occupied transcript position [28]. All resulting reads (t-signature) were reverse complemented and aligned to the miRNA identified in our study. The targets were selected and categorised as 0, 1, 2, 3, or 4.
In total, 1,442 targets of 221 known miRNAs and 117 targets of 32 novel miRNAs were identified from the degradome sequencing (target sites were detected by both target gene prediction and experimental sequencing results). Table S6 shows the sequencing results of the target genes predicted by miRNA that were detected by degradation and were statistically analysed are shown in Table S6. Novel miRNAs and their targets are shown in Table S11. Most of the miRNAs were found to cleave two or more different transcripts, 73 targets were cleaved by ath-miR5658_R-3_1ss13AT, and 31 targets were cleaved by PC-5p-284876_17; they respectively represent the highest amounts of transcripts cleaved by the same miRNAs in the known and novel miRNAs in this study. Only 35 known miRNAs and 11 novel miRNAs cleaved a single transcript target (Table S7). Furthermore, the prediction by psRNA Target combined with the mRNA generated in the density file of the degradation group allowed the number of miRNAs that interacted with the transcript to be calculated. It showed that most of the transcripts (606 of 840) were cleaved by a particular miRNA and one transcript target could be cleaved by no more than 11 miRNAs (Table S8). Targets distributed among five categories revealed that 372 targets (approximately 24%) of conserved miRNAs belonged to category 0, 21 targets (approximately 1.35%) of conserved miRNAs belonged to category 1, 675 targets (approximately 43%) of conserved miRNAs belonged to category 2, 22 targets of conserved miRNAs belonged to category 3 (approximately 1.65%), and 469 targets (about 30%) of conserved miRNAs belonged to category 4 (Table S7).
We analysed target transcripts homologous to other species in our degradation library using the BLASTX algorithm. To further understand the functions of the identified target genes, a GO functional classification was conducted. Of the biological processes, the most abundant processes were regulation of transcription and DNA-templates, transcription, DNA-templates, and auxin-activated signalling pathways. Among the cellar components, the three most highly represented were the nucleus, cytosol, and chloroplast stroma. Of the molecular function categories, the four most highly represented categories were DNA binding, transcription factor activity, and sequence-specific DNA binding, structural constituents of cell walls, and electron carrier activity. The up- and down-regulated genes and major DEGs and DEMs are shown in Fig. 4A. And the important DEGs and DEMs were mainly associated with defence responses, fruit development, signal transduction, photosynthesis, and metabolism pathways are shown in Fig. 4B.
Measurement of metabolite variation and identification-annotation through metabolomic analysis
Metabolically related DEGs and DEMs play a key role in grafting plantsand many small molecular metabolites can also be transported through the grafting site, thereby affecting the growth and development of organisms. In this study, we performed metabolomic sequencing analysis on mixed samples of fruit, phloem, and leaves, respectively, for Vv plants, Vs hybrids, and Sc plants. The LC-MS system was used to collect all samples according to the manufacturer’s instructions. The samples were scanned with negative and positive mass spectrometry. To evaluate the stability of the LC-MS during the entire acquisition, a quality control (QC) sample (pool of all samples) was acquired after every 10 samples. All detected significant metabolites are listed in Table2. Based on non-targeted metabolomics, principal component analysis (PCA) was performed to reveal the variation in the groups of metabolites. The results showed that the first two principal components accounted for 40.19% and 10.41% of the metabolic variance for all samples during negative mass spectrometry, and with positive mass spectrometry these values were 43.19% and 10.84%, respectively (Fig. 5A). According to the PCA score plot, the three groups were sharply divided in positive mass spectrometry, with Sc plants being the farthest. In negative mass spectrometry, Sc plants were clearly separate, and Vv plants and Vs hybrids had some degree of overlap, although they could still be distinguished.
Because some compounds have isomers, the first-order mass spectra could not firmly identify them, and each substance could be shattered in the mass spectrometer. These fragments could be detected by mass spectrometry to generate second-order mass spectra, which could be matched and scored with the second-order mass spectra of the standard substance in the in-house database. Finally, the second-order identification results of the metabolites were obtained. All metabolites, as well as comparative analysis and KEGG annotations of Vv-Sc, Vv-Vs, and Sc-Vs, are listed in Table S8-11. To screen for differential metabolites, univariate analysis of variance multiple (fold-change) and t tests were used to obtain a q-value by BH correction, combined with the multivariate statistical analysis of the VIP value of PLS-DA. Different ions simultaneously satisfied the following criteria: ratio > 2 or ratio < 1/2; q-value < 0.05; VIP > 1. The heatmap of different metabolites is shown in Fig. 5B.
Through metabolomic analysis, we found that one term (Table S10, M553T311) for the enrichment of gomisin was up-regulated in Vs hybrids. This component is a main medicinal ingredient in S. chinensis and have not been found in any DEGs associated with these two substances; thus, the difference in this content may be caused by direct product transport. There were seven terms related to flavonoid metabolism with significant differences in Vv plants compared to Vs hybrids in the positive mass, and three in the negative mass, which were related to fruit colour. As described above, there were also significant differences in fruits colour in Vs hybrids and Vv plants. For all significantly different metabolites, the screening software MBRole (http://csbg.cnb.csic.es/mbrole/) was used for function enrichment analysis, and the positive mass spectrometry results of Vv-Sc, Vv-Vs, and Sc-Vs are listed in Tables S8-11. These results demonstrate that differences between Vv plants and Vs hybrids were minimal at approximately 30 items, with the three most frequent categories being purine metabolism, arginine and proline metabolism, and ABC transporters. The Sc-Vs and Sc-Vv, respectively, had 68 and 65 different entries. Among the 30 significant categories in Vv-Vs, 27 had marked differences in Sc-Vv and 25 had equal or greater KEGG compounds in Sc-Vv.