Phenotypic changes of the cut positions in the grafted and separated treatments
We grew two tomato cultivars, ‘JiaHong No. 4’ (‘JH4’) and ‘JiuLv 787’ (‘JL787’) to 4-leaf stage (Fig. 1a-b). ‘JL787’ is a common rootstock cultivar in China, while ‘JH4’ is usually used as the scion. In this experiment, eight groups of plant internodes were prepared as follows (Fig. 1c-h); (1) separated ‘JH4’ shoot parts, (2) separated ‘JH4’ root parts, (3) separated ‘JL787’ shoot parts, (4) separated ‘JL787’ root parts, (5) grafted ‘JH4’ shoot parts, (6) grafted ‘JH4’ root parts, (7) grafted ‘JL787’ shoot parts, and (8) grafted ‘JL787’ root parts.
Based on our previous experience, reunion of the scion and rootstock parts of tomato seedlings could be complete < 6 days after grafting (DAG) [27]. In the autografts of both cultivars, the adhesion force between graft interfaces could be detected at 3 DAG. We therefore performed microscopic observations of internode sections near the cut position at this stage to survey phenotypic changes at the cellular level. Differences were found between the upper and lower interfaces of the separated samples (Fig. 1i-l). More parenchyma cells were observed at the upper interfaces, especially near the vascular bundles, while none were found at the lower interface. In addition, the cortex cells were unusually swollen, which could be caused by discontinuous water transport from the root to the shoot. In the grafted samples, the shoot and root parts were connected by parenchyma cells (Fig. 1m-n), which subsequently become specialized into vascular bundles.
RNA sequencing and analysis of gene co-expression networks
We performed high-throughput RNA sequencing to compare the transcriptome profiles that underlie the phenotypic changes. Because initial adhesion happened within 3 days after grafting, two time points of 48 h and 72 h were chosen with regard to the previous eight groups of internode sections. Each collecting with two biological repetitions. Therefore, a total of 32 RNA-seq libraries were constructed and named using the convention of Cultivar (C for ‘JH4’ and T for ‘JL787’) | Treatment (S for separated and G for grafted) | Position (U for above and D for below the cut site) | Coltime (collection time; 1 for 48 h and 2 for 72 h) | Biorep (biological replicate; R1 for replicate #1 and R2 for the replicate #2). The internode sections collected are illustrated in Fig. 2a. After sequencing and quality filtering of the raw reads, normalized FPKM values were obtained using the Tophat and Cufflinks analysis pipeline (see Methods). Mapping information is shown in Additional file 1: Table S1. Genes with low abundance and low variability were filtered out in order to reduce noise. Log2-transformed expression values of the remaining 3409 genes (Additional file 1: Table S2) were imported into the WGCNA package for co-expression network analysis. The overall relationships among the different samples is shown by a cluster dendrogram (Fig. 2b). The two levels of Biorep show good correlation. In addition, samples were clustered regularly based on all experimental traits, especially for Position and Treatment.
As a result of weighted co-expression network analysis, the expression profiles of the 3409 genes were grouped into 11 modules (MEs), reflecting 10 different co-expression networks (ME1-ME10) and the outliers that do not belong to any cluster (ME0). ME1–10 and ME0 are distinguished by different colors (Fig. 3a). The network structure was visualized using multi-dimensional scaling (MDS) scatter plots (Fig. 3b). Genes in the different modules were distributed in three dimensional space like fingers. Modules positioned close to one another are expected to have similar expression patterns.
In order to measure the relationships between different modules further, the eigengenes (the first principal component) of the modules were calculated in each sample set (Fig. 3c). Some module pairs indicated by eigengenes are significantly positively or negatively correlated. For example, ME1 (turquoise) is positively correlated with ME3 (brown) and ME4 (yellow), while ME1 is negatively correlated with ME2 (blue), ME5 (green), ME8 (pink), and ME10 (purple). ME9 (magenta) showed a significant positive correlation with ME2, ME5, and ME6 (red). In addition, ME7 (black) showed no correlation with any of the other colored modules.
Certain co-expression networks show obvious correlations with certain experimental traits
When associating the co-expression networks with certain experimental traits, diverse correlation patterns indicated by module eigengenes were obvious (Fig. 4a). As expected, the biological replicates did not display correlations with any module. Besides, there were no significant correlations in the Coltime column, indicating that there were no considerable changes in the gene expression profiles at 48 h and 72 h after treatment. Obviously, most of the modules showed strong correlations with Position. There were positive correlations between Position and ME1, ME3, and ME4, implying the presence of up-regulated eigengenes in these modules in the section above the cut site. In contrast, there were negative correlations between Position and ME2, ME5, ME8, and ME10, and ME9 had a lower coefficient value but the correlation was still significant. We found that no modules were negatively correlated with the most important factor, Treatment, while ME3 and ME6 showed marked positive correlations, suggesting that graft treatment may induce increased expression levels in groups of genes. Therefore, by eigengenes, we found integral features of gene co-expression networks that were related to every experimental trait.
At the same time, we also focused on the expression profiles of individual genes. We made use of Gene Significance (GS) measures, the coefficients of the Pearson correlations between individual genes and trait factors. The higher the absolute GS value, the more biologically significant is the gene. GS provides detailed information, which is especially important in the unsigned network. As shown in Fig. 4b, the GSs exhibited diverse patterns when placed in the context of the modules. The color intensities in the matrix are generally consistent with the module eigengene-trait image. Also, no GS was detected for Biorep and Coltime. The GSs for Cultivar were found mainly in ME7, with very few in ME0. As for Position, ME1, ME3, and ME4 had more positive GSs, while ME2, ME5, ME8, ME9, and ME10 had more negative GSs. In addition, the GSs for Treatment clustered mainly in ME3 and ME6, with fewer in ME9, implying biological significance for the graft treatment for the genes in these modules.
As WGCNA is a method to get a broad view of the gene expression levels in the network scale, we also performed pair-wise differential expression tests in parallel to co-expression network analysis. Differentially expressed genes were shown in Fig. 4c, corresponding to the arrangement order of individual genes in Fig. 4b. It was obvious that the GSs from co-expression network analysis and the DEGs from pair-wise tests had a very similar profile. Therefore, GSs were reliable to reflect expression patterns of particular genes.
Venn diagrams were used to visualize the overlaps between GSs of the Cultivar, Position, and Treatment traits. As shown in Fig. 4d, GS for Position comprised the largest circle, implying that the largest number of genes were differentially expressed between internode sections above and below the cut site. With the intersection areas weighted by gene numbers, it was clear that the GSs for Position and Treatment overlapped considerably with each other.
Significant biological functions are enriched in the different co-expression networks
Based on the previous results, the expression profiles of genes in modules and focused traits were displayed. Because all trait factors had two levels, genes in each module were separated into two parts, often with opposite modality (Fig. 5). For example, because ME1 was significantly related to the factor Position, genes with average values in the U samples that were larger than the average values in the D samples were labeled as ME1-U, and the remaining were labeled as ME1-D. We considered that both parts have biological meaning. Accordingly, GO annotation and enrichment analyses were performed (Fig. 5, Additional file 1: Table S3, Additional file 2: Figure. S1). Tomato genes were first used as queries. When the tomato database was not adequate, Arabidopsis ortholog genes were then used as queries. GO terms derived from the two annotations are given in Additional file 1: Table S3.
As can be seen from ME1-U and -D, the expression levels of the U and D samples diverged dramatically when they were separated, and the divergence was abbreviated from both sides when grafted. The most enriched functions in ME1-U included sterol biosynthesis, water transport, and mitotic cell cycle. And in ME1-D, except for cell redox homeostasis and xyloglucan metobolism, more items were found in annotations of Arabidopsis orthologs, including response to chitin, photomorphogenesis, response to ethylene, response to jasmonic acid, response to wounding, and organic anion transport.
In ME2-U, the expression levels of all the D samples were consistently low, and the U sample levels in the separated treatments were high, but were low in the grafted treatments. Biological processes like water deprivation, inorganic anion transport, plant organ senescence, lipid storage, and starch metabolic process were included. In ME2-D, the expression levels of all the D samples were consistently high. Genes that participate in energy and defense functions, such as photosynthesis, response to other organism, ROS metabolism, and toxin metabolism were enriched, and they showed consistently high expression in the D samples and became closer to the U samples when grafted.
In addition, expression levels in ME4-U remained high in the U samples. Few transcripts were in the separated D samples, but the number became closer to the U samples when grafted. Biological processes related to auxin transport and signaling were found to be enriched in this group of genes. In the separated samples in ME4-D, hexose biosynthesis and serine family amino acid metabolism were enriched.
Gene expression profiles in ME5, ME8, and ME10 were similar to ME1, with minor differences. Of the first two modules, sterol metabolism was the main term in the U samples, and photosynthesis was the main term in the D samples. In the ME10 module, acetyl-CoA metabolism and positive regulation of cell communication were enriched in the U and D samples, respectively.
ME3, ME6, and ME9 were the modules related to Treatment. Although ME3 and ME9 were related to the Position trait, grouping by Position and by Treatment resulted in the same two sets of genes. Uniformly, ME3-D/S, ME6-S and ME9-U/S, with S (separated) characteristics, were found to be enriched in either no or very few functions. In contrast, G samples (grafted) were found to be enriched in some specific functions in the different modules. As in ME3-U/G, genes with functions associated with cell wall biosynthesis were significant in the U samples when grafted. In ME6-G, except for L-phenylalanine or lignin metabolic processes, the most significant biological processes were related to detoxification, such as response to toxic substance, hydrogen peroxide catabolism, and response to oxidative stress. The situation was similar in ME9-D/G. These results suggest that the biological processes cell wall biosynthesis and detoxification are important in internode sections of the grafted samples.
ME7 was the only module associated with the variable Cultivar. Genes with higher expression levels in T samples (from ‘JL787’) were enriched in L-phenylalanine metabolic processes and the production of siRNA involved in RNA interference compared to the C samples (from ‘JH4’). This may imply that the rootstock cultivar ‘JL787’ has stronger biotic or abiotic resistance than does the scion cultivar ‘JH4’. Except for this, the gene expression patterns of the two cultivars were similar in all of the other co-expression clusters. However, the number of differentially expressed genes used to measure the cultivar related-effect was limited even for the two selected cultivars, because some germline-specific reads were ignored when mapping to the tomato reference genome in this work. Therefore, a more detailed effort should be made to investigate the cultivar-related effect.
In order to view the detailed expression patterns of pivotal individual genes, we assigned hub genes in each module. The top 20 genes with the largest hubness were selected to represent the modules (Fig. 6). The biological functions in which the hub genes participate were generally in agreement with the enriched GO terms, but provided more detailed information. For example, ME4 was enriched in the auxin signaling and efflux pathways, and included the candidate hub genes PIN1 (Solyc03g118740.2), GH9C2 (Solyc12g055970.1) and IAA9 (Solyc04g076850.2).