Root transcriptomic responses of grafted grapevines to heterogeneous nitrogen availability depend on rootstock genotype.

In many fruit species, including grapevine, grafting is used to improve scion productivity and quality and to adapt the plant to environmental conditions. However, the mechanisms underlying the rootstock control of scion development are still poorly understood. The ability of rootstocks to regulate nitrogen uptake and assimilation may contribute to this control. A split-root system was used to grow heterografted grapevines and to investigate the molecular responses to changes in nitrate availability of two rootstocks known to affect scion growth differently. Transcriptome profiling by RNA sequencing was performed on root samples collected 3 and 24 h after nitrogen supply. The results demonstrated a common response involving nitrogen-related genes, as well as a more pronounced transcriptomic reprogramming in the genotype conferring the lower scion growth. A weighted gene co-expression network analysis allowed the identification of co-regulated gene modules, suggesting a role for nitrate transporter 2 family genes and some transcription factors as main actors controlling this genotype-dependent response to heterogeneous nitrogen supply. The relationship between nitrate, ethylene, and strigolactone hormonal pathways was found to differ between the two genotypes. These findings indicated that the genotypes responded differently to heterogeneous nitrogen availability, and this may contribute to their contrasting effect on scion growth.

Introduction

Grafting is increasingly used in fruit species and vegetable crops to improve crop yield and to reduce disease susceptibility (Goldschmidt, 2014; Warschefsky ). In viticulture, this practice is widely used in Europe, allowing grape growth (Vitis vinifera L.) in soils infected by phylloxera since the 19th century. In addition to their contribution to improving root traits such as resistance to many soil-borne pests, rootstocks also affect scion growth, biomass allocation, and mineral nutrition, resulting in differences in yield and fruit quality (n class="Chemical">Csikasz-Krizsics and Diofasi, 2008; Dalbo ; Ollat ). The influence of rootstocks on the development and productivity of the scion is an important trait to consider in the selection of rootstock genotypes in agriculture. Water and/or nutrient uptake and transport, hormonal regulation, and long-distance signalling have been investigated to explain the control of vegetative development of the scion by the rootstock in various species, including grapevine (Albacete ; Zhang ). In grapevine, among the different mechanisms involved, it has been shown that the heterografting process induces profound transcriptomic changes in the scion (Cookson and Ollat, 2013) and that the exchange of mRNA between the scion and rootstock could be influenced by the environment and the genotype (Yang ). The rootstock genotype affects the nitrogen (N) status of the grafted plant, and in some cases this effect was correlated with the ability of the rootstock to control scion vigour (Nikolaou ; Zerihun and Treeby, 2002). A recent study demonstrated that rootstock control of scion growth was dependent on N supply and correlated with distinct root and leaf ionome profiles (Lecourt ).

N nutrition is one of the major factors that influence plant growth. In aerobic soils, n class="Chemical">nitrate is the main source of N for most plants and represents an important regulator of biomass allocation, acting not only as a nutrient but also as a signal (Crawford and Glass, 1998). As the concentration of N in the soil is heterogeneous, plants have evolved adaptive mechanisms of N sensing and signalling, involving local and systemic signals, to utilize this resource efficiently (O’Brien ). In Arabidopsis thaliana, the primary nitrate response involves the transcriptional regulation of some genes, which occurs within minutes after a change in local nitrate supply (Medici and Krouk, 2014). As part of this response, the nitrate transceptor NPF6.3 controls the transcriptional regulation of nitrate-responsive genes such as NRT2.1 (Gojon ). This response is modulated by phosphorylation/dephosphorylation events involving the CBL-interacting protein kinases CIPK8 and CIPK23 (Bouguyon ) and protein phosphatase 2C family members, such as ABI2 (Léran ). Several nitrate-responsive gene regulators (LBD37/38/39, NLP7, SPL9 or NRG2) have also been identified (Vidal ; Xu ).

Nitrate uptake and assimilation are affected by not only local but also systemic signals coming from other parts of the plant, such as distant roots or shoots, to integrate the distribution of n class="Chemical">nitrate content in the soil with the nutritional and energy status of the plant (Nacry ; Krapp ; Krapp, 2015; O’Brien ). The plant’s response to nitrate content heterogeneity includes morphological and physiological modifications that allow adaptive root foraging to ensure optimal N uptake. Several studies on the model plant Arabidopsis based on a split-root system have shown that the differential of root growth, that is, growth increase in N-rich zones and growth limitation in N-poor zones, is controlled by both local and systemic signals integrating systemic N demand and N supply signals (Ruffel ; Alvarez ; Bouguyon ). Some transcription factors involved in this regulation have been identified, such as ANR1, a MADS-box transcription factor (Gan ; Mounier ). NPF6.3, which acts upstream of ANR1, also plays an important role in the control of lateral root growth through the regulation of auxin export from lateral root primordia (Gojon ). In addition to auxin, cytokinins, abscisic acid, and ethylene hormonal pathways, as well as peptide hormones such as C-TERMINALLY ENCODED PEPTIDE, are also controlled by nitrate concentration, and in turn regulate morphological and physiological responses to nitrate (Kiba ; Tabata ; Krouk, 2016; Ohkubo ).

Recently, the molecular networks that control N responses in relation to homogeneous or heterogeneous n class="Chemical">nitrate availability have been explored, mostly in herbaceous species (Vidal ; O’Brien ). Systems biology approaches performed in Arabidopsis have led to the identification of some master regulators that integrate systemic signals (Vidal , 2013; Alvarez ; Canales ; Li ). In contrast, there is very little information available about the global transcriptomic changes induced by N supply in woody species, and even less in grafted perennial plants.

The aim of the present study was to determine whether two grapevine rootstocks [Vitis riparia cv. Riparia Gloire de Montpellier (RGM) and Vitis berlandieri × Vitis rupestris hybrid cv. 1103 Paulsen (1103P)], known to confer low and high scion vigour, respectively (Lecourt ; Zhang ), differ in their n class="Chemical">nitrate perception and signalling in response to a heterogeneous nitrate supply. These rootstocks were grafted with a scion of Vitis vinifera cv. Cabernet Sauvignon (CS). Grafted plants were grown in a split-root system under low N conditions in both root compartments and subjected to 24 h of a localized N supply (in one compartment). A global transcriptomic analysis by RNA sequencing (RNA-seq) was conducted on roots collected 3 and 24 h post-treatment (hpt). The results showed that for both rootstocks, the N-related gene expression pattern was extensively regulated in response to a local N supply. Interestingly, CS/RGM presented many more differentially expressed genes (DEGs) than CS/1103P according to the N treatment, and a weighted gene co-expression network analysis (WGCNA) approach allowed the identification of N-related modules specific to the CS/RGM combination.

Materials and methods

Plant materials and split-root experiment

The two rootstock genotypes, RGM and 1103P, were both grafted with the scion genotype CS. A double-grafting system (Tandonnet ) was used to obtain two-roots–one-shoot plants. Grafted plants resembled an inverted ‘Y’ (Fig. 1A). The combinations n class="Chemical">CS/RGM and CS/1103P were analysed. After callusing and rooting, grafted plants were cultivated in two sand-filled pots of 3 l capacity (one per root system) in a greenhouse and irrigated with a full nutrient solution (Cookson ) for 24 d of acclimation (Fig. 1A). A low-nitrate-content nutrient solution [LN: 0.8 mM KNO3, 0.57 mM K2HPO4, 0.69 mM MgSO4, 1.39 mM CaCl2, 0.8 mM K2SO4, 0.3 mM CaSO4, and micronutrients as described by Lecourt ] was then applied for 2 weeks. At 0 hpt, one side of the root system was irrigated with 1 l of 5 mM high-nitrate (HN) solution (‘HN roots’) and the other side was irrigated with LN solution (‘LN roots’). The HN solution was produced from LN solution by adding KNO3 up to 5 mM, and adjusting the K2SO4 and CaSO4 concentrations to equilibrate the K+ between LN and HN solutions. Root tips (the first 2 cm from the root cap) were harvested at 0, 3, and 24 h after application of the HN solution, immediately frozen in liquid nitrogen, and stored at –80 °C. Three plants (biological replicates) of each combination were sampled at each time point, with two root samples per plant (LN and HN roots) for 3 and 24 hpt (Fig. 1B, C).

Fig. 1.

Split-root experiment in the greenhouse. (A) A scion was grafted on to two rootstock cuttings of the same genotype. Both sides of the root system were grown in two separated pots (left). The arrow indicates the grafting point. Right panel, photograph of a root system harvested during the experiment. (B) Diagram outlining the experimental protocol. The whole root system was supplied with LN solution (0.8 mM nitrate) for 2 weeks. Root samples from three plants per combination were then harvested (0 hpt samples). At the same time, one side of the root system of 6 other plants was supplied with 1 l of HN solution (5 mM nitrate). Root tip samples from each part of the root system (two root samples per plant: HN roots and LN roots) were harvested 3 or 24 h after HN supply (3 and 24 hpt samples). (C) Nomenclature of harvested samples. Three individual plants (biological replicates) were harvested per condition. CS, Cabernet Sauvignon; RGM, Riparia Gloire de Montpellier; 1103P, 1103 Paulsen. (This figure is available in colour at JXB online.)

Split-root experiment in the greenhouse. (A) A scion was grafted on to two rootstock cuttings of the same genotype. Both sides of the root system were grown in two separated pots (left). The arrow indicates the grafting point. Right panel, photograph of a root system harvested during the experiment. (B) Diagram outlining the experimental protocol. The whole root system was supplied with LN solution (0.8 mM n class="Chemical">nitrate) for 2 weeks. Root samples from three plants per combination were then harvested (0 hpt samples). At the same time, one side of the root system of 6 other plants was supplied with 1 l of HN solution (5 mM nitrate). Root tip samples from each part of the root system (two root samples per plant: HN roots and LN roots) were harvested 3 or 24 h after HN supply (3 and 24 hpt samples). (C) Nomenclature of harvested samples. Three individual plants (biological replicates) were harvested per condition. CS, Cabernet Sauvignon; RGM, Riparia Gloire de Montpellier; 1103P, 1103 Paulsen. (This figure is available in colour at JXB online.)

RNA extraction, library preparation, and sequencing

Total RNAs were isolated from 30 samples (three biological replicates per condition; the conditions are illustrated in Fig. 1C), as described by Cookson . Total RNAs were sent to IGBMC Microarray and Sequencing Platform (Strasbourg, France) and their quality was further determined using an Agilent 2100 Bioanalyzer (Agilent Technologies). The preparation of the libraries was performed by the IGBMC Microarray and Sequencing Platform using the TruSeq RNA Sample Prep Kit (Illumina Inc.). Libraries were sequenced with an Illumina HiSeq 2000 sequencer (Illumina Inc.) and 50 bp single reads were generated. The data have been deposited in the National Center for Biotechnology Information (NCBI) Sequence Read Archive (http://www.ncbi.nlm.nih.gov/sra, last accessed 19 June 2017) under accession number PRJNA342391.

Pre-processing of RNA-seq data

Quality checks were done with fastQC on the Galaxy website hosted on the Abims Bioinformatics Platform (Roscoff, France) (http://abims.sb-roscoff.fr/resources/galaxy, last accessed 19 June 2017). The raw reads were then filtered by the Genotoul Bioinformatin class="Chemical">cs Platform (Toulouse, France) with trimming of the low-quality reads with an average Phred score <28 (see Supplementary Table S1 at JXB online).

De novo assembly

De novo assembly was performed by the Genotoul Bioinformatics Platform using DRAP v1.3 (Cabau ). This meta-assembler first performed an assembly with Trinity (Grabherr ) or Oases (Schulz ) using the multi-kmers approach. After assembly compaction, poorly supported and misassembled transcripts were filtered out using a reads per kilobase per million mapped reads (RPKM) threshold of 1, 3, or 10 (Supplementary Table S2). Considering the results of the different assemblies (Supplementary Table S2), the assembly with Oases multi-kmers and an RPKM threshold of 1 was selected.

Functional annotation of de novo assembled transcriptome

Functional annotation was performed by the Genotoul Bioinformatics Platform with series of BLASTX against different protein databases from Ensembl hosted on the Genotoul server. Contigs were post-filtered on read alignment and the best-hit annotation was selected for each of them. Gene Ontology (GO) terms were predicted using Interproscan 5. All the data were available on a web interface RNAbrowse (http://ngspipelines2.toulouse.inra.fr:9005/ngspipelines/#!/NGSpipelines/Vitis%20rupestris%20-%20RGMand1103P, last accessed 19 June 2017), which is described as an RNA-seq de novo assembly results browser (Mariette ). The best annotation file was downloaded and a BIN code assignment was performed using Mercator (Lohse ). Enrichments of functional categories of the Mercator annotation in the significantly differentially expressed gene sets were tested for significance by Fisher’s exact tests with a Bonferroni correction for multiple tests, using Mefisto version 0.23beta (http://www.usadellab.org, last accessed 19 June 2017).

Quality evaluation

To assess the reliability of the de novo assembled transcriptome, raw reads were aligned against it using the aligner BWA MEM (Li and Durbin, 2009). The completeness of the assembly was determined through the use of the Core Eukaryotic Genes Mapping Approach (CEGMA) (http://korflab.ucdavis.edu/Datasets/cegma/, last accessed 20 July 2016) (Parra ), which allows comparison between the transcriptome and a dataset containing highly conserved annotated core proteins.

Analysis of differential expression

Raw counts were determined using BWA MEM and SAMTools (Li ), generating an overall counts matrix. Then, the R package edgeR (Robinson ) was used to identify differentially expressed genes using a stringent threshold: absolute value of Log Fold Change (LFC) >1 and False Discovery Rate (FDR) <0.01.

Network analysis

A co-expression gene network was constructed using the WGCNA software package (v1.51) in R (Langfelder and Horvath, 2008) using all the libraries except those corresponding to time 0 hpt (24 libraries). In order to remove low-expressed contigs (reflecting noise), a filter was applied to keep only contigs that had at least 20 counts in 70% of the libraries. A total of 45358 contigs satisfying the above threshold were obtained, and their counts data were transformed using the function varianceStabilizingTransformation of the package DESeq2 (Love ). The resulting set of counts was used for network construction and module detection using the function blockwiseModules. Briefly, an adjacency matrix was created by calculating the biweight mid-correlation raised to a power β of 8 (soft-threshold estimated with the pickSoftThreshold function) and the maxPoutliers parameter set to 0.05. The subsequent Topological Overlap Matrix (TOM) was used for module detection using the DynamicTreecut algorithm with a minimal module size of 30 and a branch merge cut height of 0.25. The module eigengenes were used to evaluate the association between the 26 resulting modules and traits (genotype, treatment, and time).

Validation of RNA-seq analysis by quantitative real-time PCR

Total RNAs were reverse transcribed into cDNA using the SuperScript III First-Strand Synthesis System for RT-PCR (Invitrogen). Quantitative real-time PCR (qPCR) reactions were performed using SYBR Green on an iCycler iQH (Bio-Rad), according to the procedure described by the supplier. The relative expression of the genes was calculated using the 2–ΔΔCt method (Livak and Schmittgen, 2001), and using the reference genes EF1γ and GAPDH for normalization. Primer sequences are listed in Supplementary Table S3.

Nitrate measurements

Total nitrate concentration was measured according to the modified method of Cataldo . Samples (150 mg) of whole root were ground and mixed in 400 µl 7.5% n class="Chemical">trichloroacetic acid, and centrifuged at 16000 g for 5 min. Then, 200 µl 5% salicylic acid in concentrated H2SO4 was added to 62.5 µl of the supernatant, and after 20 min at room temperature, 4.75 ml 2M NaOH was added. The absorbance at 410 nm of each sample was measured using an Epoch 2 microplate spectrophotometer (Bio Tek Instruments). A standard curve was constructed using a 0.25 g l–1 NO3– N solution. For each sample, a negative control was done where the salicylic acid solution was replaced with 0.2 ml concentrated H2SO4.

Results

De novo assembly and annotation of a root transcriptome from grafted Vitis plants

To obtain a comprehensive overview of the transcriptome changes in response to rootstock genotype and/or N supply, 30 cDNA libraries were constructed from root tips of two scion/rootstock combinations, CS/1103P and CS/RGM, grown in a split-root system (Fig. 1). The sequencing was performed using an Illumina Hiseq 2000 sequencer and approximately 40 million 50 bp single reads were obtained for each sample (Supplementary Table S1). More than 97% of the raw reads passed the quality evaluation based on a Phred score threshold.

In grapevine, the reference genome is that of PN40024, a n class="Species">V. vinifera accession (Jaillon ). Since the two rootstocks used in this study belong to different species of the genus Vitis, relying on the reference genome may underestimate the variability among the genotypes and may not allow the identification of specific genes or isoforms. Therefore, a de novo assembly strategy was chosen in order to build the most comprehensive transcriptome allowing differential expression analysis between the two studied genotypes. Briefly, a total of 1.19 billion cleaned reads from the 30 cDNA libraries were assembled into contigs using the Oases assembler, with an RPKM threshold of 1, which allowed a higher percentage of contigs showing protein assignments (Supplementary Table S2). A total of 51250 contigs (47519774 bp) were generated, with an N50 at 1481 (Supplementary Table S2). The assembly generated a high number of transcripts enriched for smaller-sized contigs (51.56% of contigs were in the size range 201–599 bp; Supplementary Fig. S1).

To assign putative gene functions and annotations to all the assembled contigs, BLASTX alignment was performed against protein databases such as Swissprot, Refseq, and several Genotoul bioinformatics platform internal databases of different species, including n class="Species">Medicago truncatula, Theobroma cacao, Populus trichocarpa, A. thaliana, Prunus persica, and V. vinifera, using a cut-off E-value of 1e-5. The best hit for each contig was selected and they were further categorized by GO functional annotation (Fig. 2). Among the 51250 contigs, 17683 were associated to at least one GO category, generating 48257 GO annotations. Furthermore, the BLASTX series provided insights into the taxonomic distribution of the transcripts, with, as expected, more than 75% of contigs having top hits to sequences from V. vinifera (Supplementary Fig. S2). This species distribution suggested that the assembly and annotation of this Vitis root transcriptome were correct and reliable. To complete the annotation process, Mercator, another classification tool for next-generation sequencing data (Lohse ), was used to assign BIN codes. In total, 56268 annotations were generated, corresponding to 51012 contigs covering almost all the transcriptome (Fig. 3).

Fig. 2.

Gene ontology (GO) distribution of the de novo merged transcriptome The distribution of the contigs among level 1 GO categories of biological process, cellular component, and molecular function is shown. Molecular function was the category that contained the highest number of contigs of the assembled transcriptome.

Fig. 3.

Classification of the de novo merged transcriptome into BIN code classes. The percentage of contigs annotated in each BIN code compared with the total number of annotations generated by Mercator is presented. Category number 35, ‘not assigned’, is not shown, and contained 44.7% of the assignments.

Classification of the de novo merged transcriptome into BIN code classes. The percentage of contigs annotated in each BIn class="Chemical">N code compared with the total number of annotations generated by Mercator is presented. Category number 35, ‘not assigned’, is not shown, and contained 44.7% of the assignments.

In order to examine the representation of the transcriptome of the RNA-seq reads, reads were mapped back to the de novo assembly, and showed more than 90% of mapping for all the samples (data not shown). Finally, the completeness of the de novo assembled transcriptome was evaluated by comparison of its 51250 contigs with a set of highly conserved core proteins using the CEGMA pipeline (Parra ), allowing the identification of 98.39% of the 248 core eukaryotic genes contained in this core dataset.

Heterogeneous nitrate supply induces a larger modification in the root transcriptome of CS/RGM compared with CS/1103P

Representation of each sample in a multidimensional scaling plot view (Supplementary Fig. S3) showed a clear clustering of the samples depending on the rootstock genotype and the N treatment. Moreover, the separation between Ln class="Chemical">N and HN root transcriptomes was more pronounced for the combination CS/RGM than for CS/1103P. DEGs were then identified by comparing LN and HN roots at 3 or 24 hpt for the two combinations, with the threshold of |LFC|>1 and FDR<0.01 (Supplementary Table S4). As shown in Fig. 4, the number of DEGs in response to high local N treatment was very different in the two combinations. CS/RGM exhibited a much higher number of DEGs (1369 DEGs in total; Fig. 4B) than CS/1103P (212 DEGs) (Fig. 4A). In CS/1103P, the majority of DEGs (205) was found at 3 hpt, and 132 of them were up-regulated in HN roots compared with LN roots. In the combination CS/RGM, most of the DEGs at 3 hpt (81%) were up-regulated in HN roots, while the majority of the DEGs found at 24 hpt (76%) were down-regulated.

Fig. 4.

Analysis of the differentially expressed genes (DEGs) in roots in response to heterogeneous N availability in the two scion/rootstock combinations. Venn diagrams show the number of transcripts that were up-regulated (bold) and down-regulated (grey) in HN roots compared with LN roots at 3 and 24 hpt in (A) CS/1103P plants and (B) CS/RGM plants. In (B), in addition to the 1369 contigs that were differentially expressed, three contigs were duplicated since they followed a different regulation pattern (i.e. up- or down-regulated) between 3 and 24 hpt.

Analysis of the differentially expressed genes (DEGs) in roots in response to heterogeneous N availability in the two scion/rootstock combinations. Venn diagrams show the number of transcripts that were up-regulated (bold) and down-regulated (grey) in n class="Chemical">HN roots compared with LN roots at 3 and 24 hpt in (A) CS/1103P plants and (B) CS/RGM plants. In (B), in addition to the 1369 contigs that were differentially expressed, three contigs were duplicated since they followed a different regulation pattern (i.e. up- or down-regulated) between 3 and 24 hpt.

Comparing the lists of the DEGs allowed the identification of 172 common DEGs between CS/1103P and n class="Chemical">CS/RGM (Supplementary Table S4), among which 115 transcripts exhibited the same response to N treatment in both combinations (Fig. 5A). At 24 hpt, only four transcripts shared a common pattern, one being up-regulated and three down-regulated. The majority of the DEGs (54%) found in CS/1103P were common to those in CS/RGM. However, CS/RGM possessed a higher number of genotype-specific DEGs (Fig. 4), and only 8% of the total number of DEGs found in CS/RGM was shared with CS/1103P.

Fig. 5.

A core set of N-related genes was shared between CS/1103P and CS/RGM in response to heterogeneous N availability. (A) Venn diagram of the common DEGs following the same expression pattern at 3 and 24 hpt in both rootstock genotypes (up-regulation is represented in orange text and down-regulation in purple). (B) BINs enriched in the list of common DEGs (n=172). The Contingency column shows the number of genes (i) in the BIN in the input list, (ii) in the background, (iii) not in the BIN in the input list, and (iv) not in the background. P-values were adjusted with a Bonferroni correction. Values were filtered with an adjusted P-value threshold <0.01 and an enrichment >1. (C) Hierarchical clustering of the transcripts and the different conditions (i.e. genotype × treatment × harvesting time) in a heatmap presenting the expression pattern of each DEG (rows) within the different conditions (columns). The harvesting time 0 hpt was excluded during the hierarchical clustering process. Expression values are RPKM log2-transformed with up-regulation to down-regulation varying from orange to purple. Right panel, transcript clusters were extracted using the gene hierarchical clustering tree. The x-axis of each plot represents the harvesting time (in hpt); the y-axis represents the mean-centred RPKM log2-transformed values. The HN condition is indicated by dark-coloured solid lines and the LN condition by light-coloured dashed lines. For each cluster by genotype, mean±SE values are represented in black.

A core set of N-related genes was shared between n class="Chemical">CS/1103P and CS/RGM in response to heterogeneous N availability. (A) Venn diagram of the common DEGs following the same expression pattern at 3 and 24 hpt in both rootstock genotypes (up-regulation is represented in orange text and down-regulation in purple). (B) BINs enriched in the list of common DEGs (n=172). The Contingency column shows the number of genes (i) in the BIN in the input list, (ii) in the background, (iii) not in the BIN in the input list, and (iv) not in the background. P-values were adjusted with a Bonferroni correction. Values were filtered with an adjusted P-value threshold <0.01 and an enrichment >1. (C) Hierarchical clustering of the transcripts and the different conditions (i.e. genotype × treatment × harvesting time) in a heatmap presenting the expression pattern of each DEG (rows) within the different conditions (columns). The harvesting time 0 hpt was excluded during the hierarchical clustering process. Expression values are RPKM log2-transformed with up-regulation to down-regulation varying from orange to purple. Right panel, transcript clusters were extracted using the gene hierarchical clustering tree. The x-axis of each plot represents the harvesting time (in hpt); the y-axis represents the mean-centred RPKM log2-transformed values. The HN condition is indicated by dark-coloured solid lines and the LN condition by light-coloured dashed lines. For each cluster by genotype, mean±SE values are represented in black.

Both CS/RGM and CS/1103P share a core set of genes differentially regulated in response to N variation

The 172 common DEGs shared between both rootstock genotypes were mainly associated with hormone metabolism and N-related functional categories (Fig. 5B). In both combinations, the n class="Chemical">HN treatment altered the expression of genes linked to N metabolism and transport, oxidative pentose phosphate (OPP), tetrapyrrole synthesis, and RNA regulation, as shown by the enrichment analysis (Fig. 5B). A double hierarchical clustering performed on these common DEGs highlighted two major clusters for samples and three major clusters for gene expression profiles. The two clusters for samples were mainly separated by the N conditions for both genotypes (Fig. 5C). With regard to the transcript expression profiles, most of the genes associated with N metabolism functional categories were present in cluster 1 (Fig. 5C), showing an up-regulation in HN roots compared with LN roots at 3 hpt. The 15 transcripts of cluster 2 were also induced in HN roots only at 3 hpt in both genotypes. However, these transcripts were more highly expressed in CS/1103P compared with CS/RGM. Lastly, cluster 3 contained 63 transcripts, most of which encoded putative proteins associated with ethylene hormone metabolism, a functional category that showed significant enrichment (Fig. 5B). Interestingly, this cluster contained 57 transcripts that showed a different regulation pattern in response to N treatment between the two combinations (i.e. they are not presented in the Venn diagram in Fig. 5A). In CS/1103P, these transcripts were down-regulated only at 3 hpt after the HN treatment, while they were differentially regulated specifically at 24 hpt in CS/RGM between LN and HN roots.

N treatment regulates gene expression in a rootstock-dependent manner

In addition to the 172 common DEGs, each scion/rootstock combination possessed a specific set of DEGs, since 212 and 1369 transcripts were found to be differentially expressed in CS/1103P and n class="Chemical">CS/RGM, respectively (Fig. 4A, B). The enriched functional categories present in the total list of DEGs in response to the N treatment at 3 and 24 hpt in CS/1103P (Supplementary Fig. S4A) were almost the same as those found for the common DEGs (Fig. 5B). Among the 40 DEGs specific to CS/1103P, 25 were annotated as ‘not assigned’, and the others were present in clusters 1 and 2 (8 and 7 DEGs, respectively). They were associated to different functional categories, such as cell wall, miscellaneous, and development. Two ERF/AP2 transcription factors were identified and these belonged to cluster 1 (Supplementary Table S4; Supplementary Fig. S4B).

In CS/RGM, the enrichment analysis performed with all the DEGs showed functional categories that were not significantly enriched either in the common list of DEGs or in the n class="Chemical">CS/1103P-specific DEGs (Fig. 6A). Categories related to abiotic stress, fermentation, secondary metabolism (flavonoids), and gibberellin metabolism were among them (Fig. 6A; Supplementary Fig. S5A, B). The hierarchical clustering of the 1369 DEGs allowed the identification of nine clusters (Fig. 6B). Genes related to fermentation and secondary metabolism were present in cluster 4, and followed the same pattern as N-related genes that were up-regulated at 3 hpt in HN roots. Cluster 1 transcripts showed decreased expression at 3 hpt in the LN roots and were mainly linked to regulation of transcription. Cluster 6, which contained 40% of the DEGs, was related to hormone metabolism (mainly ethylene and gibberellin), RNA regulation of transcription, and signalling. Several gibberellin-metabolism-associated genes were also included in cluster 8. The remaining 14% of the CS/RGM DEGs were separated into five clusters, and most of them were not assigned.

Fig. 6.

In response to higher N availability, six-fold more genes were differentially expressed in CS/RGM compared with CS/1103P. (A) Enrichment analysis of the DEGs from CS/RGM (n=1369). The Contingency column shows the number of genes (i) in the BIN in the input list, (ii) in the background, (iii) not in the BIN in the input list, and (iv) not in the background. P-values were adjusted with a Bonferroni correction. Values were filtered with an adjusted P-value threshold <0.01 and an enrichment >1. (B) Hierarchical clustering of the transcripts and the different conditions (i.e. treatment × harvesting time) in a heatmap presenting the expression pattern of each DEG (rows) within the different samples (columns). The harvesting time 0 hpt was excluded during the hierarchical clustering process. Expression values are RPKM log2-transformed with up-regulation to down-regulation varying from orange to purple. Right panel, transcript clusters were extracted using the gene hierarchical clustering tree. The x-axis of each plot represents the harvesting time (in hpt); the y-axis represents the mean-centred RPKM log2-transformed values. The HN condition is indicated by dark-coloured solid lines and the LN condition by light-coloured dashed lines. Mean±SE values are represented in black.

In response to higher N availability, six-fold more genes were differentially expressed in n class="Chemical">CS/RGM compared with CS/1103P. (A) Enrichment analysis of the DEGs from CS/RGM (n=1369). The Contingency column shows the number of genes (i) in the BIN in the input list, (ii) in the background, (iii) not in the BIN in the input list, and (iv) not in the background. P-values were adjusted with a Bonferroni correction. Values were filtered with an adjusted P-value threshold <0.01 and an enrichment >1. (B) Hierarchical clustering of the transcripts and the different conditions (i.e. treatment × harvesting time) in a heatmap presenting the expression pattern of each DEG (rows) within the different samples (columns). The harvesting time 0 hpt was excluded during the hierarchical clustering process. Expression values are RPKM log2-transformed with up-regulation to down-regulation varying from orange to purple. Right panel, transcript clusters were extracted using the gene hierarchical clustering tree. The x-axis of each plot represents the harvesting time (in hpt); the y-axis represents the mean-centred RPKM log2-transformed values. The HN condition is indicated by dark-coloured solid lines and the LN condition by light-coloured dashed lines. Mean±SE values are represented in black.

Both rootstock genotypes showed regulation of hormone metabolism genes in response to N treatment, but with some differences. The distribution of all the functional categories within the set of DEGs at 3 hpt in each combination (Supplementary Fig. S5A, B) showed that this category was over-represented only for n class="Chemical">CS/1103P, and contained 14% of the DEGs in that combination, while it was represented by only 2.7% of the DEGs in CS/RGM. In both combinations, most of the genes related to this category were down-regulated in HN compared with LN roots.

Taken together, these results highlighted that RGM differs from 1103P in response to a heterogeneous N supply with a higher number of genes involved and a specific pattern of expression for genes involved in hormonal or secondary metabolism pathways.

Genes encoding nitrate transporters and enzymes involved in N metabolism are highly regulated in CS/RGM grafted plants

In order to understand the differences between the two rootstock genotypes in their response to changes in N availability, we focused on the expression pattern of genes involved in n class="Chemical">N nutrition (absorption, metabolism, transport, and signalling). Among the genes known to be regulated by nitrate in several species, transcripts corresponding to 18 genes were found to be significantly up- or down-regulated in response to heterogeneous nitrate supply in at least one condition (Table 1). Genes encoding nitrate assimilation-related enzymes such as nitrate and nitrite reductases (NR and NIR, respectively), glutamine synthetase (GS2), and uroporphyrin-III C-methylase (UPM1), and the nitrate transporters NRT2.4 (NRT2.4a and NRT2.4b) and NRT3 were significantly up-regulated in HN compared with LN roots. For each of the transcripts corresponding to these genes, the LFC value was higher in CS/RGM than in CS/1103P (Supplementary Table S4). CS/RGM was more affected by the N availability modification, as the transcripts corresponding to these 18 genes were all significantly differentially expressed only in this combination, while only seven of them exhibited an |LFC|>1 (and an FDR<0.01) in CS/1103P. CS/RGM-specific responsive transcripts encoded the nitrate transporter NRT2.5 and NPF transporters such as NPF6.3, the cytosolic isoform of the glutamine synthetase (GSR), and 6-phosphogluconate dehydrogenase (6PGDH). Some genes were down-regulated in HN roots compared with LN roots specifically in CS/RGM (AMT, NPF2.13, NPF3.1, and NPF4.5).

Table 1.

List of genes associated with N metabolism and transport that were differentially expressed between the LN and HN conditions in CS/1103P and/or CS/RGM at 3 and 24 h post-treatment (hpt)

Name	Function	Contig name	CRIBI accession v1	CS/1103P		CS/RGM
				HN vs. LN		HN vs. LN
				3 hpt	24 hpt	3 hpt	24 hpt
6PGDH	6-phosphogluconate dehydrogenase	mix_LOC100241717	VIT_02s0025g00900	0.91	0.5	1.85	1.04
AMT3.1	Ammonium transporter 3.1	mix_LOC100252515	VIT_07s0031g02950	–0.97	0	–1.3	–2.31
AMT3.3	Ammonium transporter 3.3	mix_LOC100248822.2.2	VIT_08s0058g00140	–0.58	0	–1.1	–1.46
NR	Nitrate reductase	mix_LOC100264320	VIT_18s0001g03910	1.21	0	2.02	0.84
NIR1	Nitrite reductase 1	mix_contig_09056	VIT_03s0063g00370	1.24	0	2.39	1.42
		mix_contig_10888	VIT_03s0063g00370	1.37	0.66	2.5	1.45
GLT1	Glutamate synthase	mix_LOC100246868	VIT_16s0098g00290	0.69	0.84	1.48	1.32
GS2	Glutamine synthetase	mix_contig_00751	VIT_05s0020g02480	1.09	0	1.83	0.93
		mix_LOC100261413.1.2	VIT_05s0020g02480	1.05	0	1.82	0.91
		mix_LOC100261413.2.2	VIT_05s0020g02480	1.03	0	1.96	0.92
GSR	Glutamate synthase	mix_contig_00892	VIT_14s0006g00350	0.72	0.63	1.2	0.55
		mix_contig_01347	VIT_17s0000g01910	0.69	0.58	1.06	0
		mix_contig_02263	VIT_14s0006g00350	0.72	0.64	1.21	0.57
		mix_contig_02304	VIT_14s0006g00350	0.7	0	1.28	0.6
		mix_contig_06858	VIT_14s0006g00350	0.7	0.6	1.14	0
		mix_GLNA2.1.7	VIT_14s0006g00350	0.84	0	1.4	0
		mix_GLNA2.2.7	VIT_14s0006g00350	0.69	0.59	1.26	0.63
		mix_GLNA2.3.7	VIT_14s0006g00350	0.69	0.61	1.3	0.6
		mix_GLNA2.4.7	VIT_14s0006g00350	0.72	0.61	1.2	0
		mix_GLNA2.5.7	VIT_15s0024g01530	0.69	0.57	1.17	0.56
		mix_GLNA2.6.7	VIT_14s0006g00350	0.78	0.64	1.16	0
		mix_GLNA2.7.7	VIT_14s0006g00350	0.74	0.61	1.26	0.58
LBD39	LOB domain-containing protein 39	mix_LOC100261250.1.3	VIT_07s0129g00330	0.94	0.64	1.51	1.05
		mix_LOC100261250.2.3	VIT_07s0129g00330	0.62	0	1.13	0.65
		mix_LOC100261250.3.3	VIT_07s0129g00330	0.74	0	1.27	0.79
NPF2.13	NRT1/PTR FAMILY 2.13	mix_LOC100250071	VIT_01s0026g01490	–0.58	0	–1.77	–3.03
NPF3.1	NRT1/PTR FAMILY 3.1	mix_LOC100250961	VIT_01s0011g03400	0	0	–1.07	0
NPF4.5	NRT1/PTR FAMILY 4.5	mix_contig_05016	VIT_18s0001g11280	0	–0.54	0	–1.39
NPF6.3	NRT1/PTR FAMILY 6.3	mix_contig_00176	VIT_02s0154g00260	0.63	0.78	1.05	1.45
		mix_contig_09300	VIT_02s0154g00260	0.84	0	1.41	1.51
NRT2.4a	Nitrate transporter 2.4a	mix_contig_00726	VIT_06s0061g00320	1.29	0.97	2.23	0.95
		mix_contig_09409	VIT_06s0061g00320	1.34	0.95	2.28	1
		mix_LOC100241340	VIT_06s0061g00320	1.31	0.95	2.25	0.93
NRT2.4b	Nitrate transporter 2.4b	mix_LOC100263699.1.2	VIT_08s0040g01500	1.4	0.87	2.27	0.89
		mix_LOC100263699.2.2	VIT_08s0040g01500	1.66	0	2.18	0
NRT2.5	Nitrate transporter 2.5	mix_LOC100260250	VIT_01s0127g00070	0.5	–0.51	1.08	–0.69
NRT3.1	Nitrate transporter 3.1	mix_contig_00610	VIT_17s0000g09470	1.32	0.88	1.97	1.02
		mix_LOC100258771.1.2	VIT_17s0000g09470	1.12	0.75	1.8	0.99
		mix_LOC100258771.2.2	VIT_17s0000g09470	1.3	0.96	1.95	1.04
		mix_NAR21	VIT_17s0000g09470	1.33	0.94	1.94	1.03
UPM1	Uroporphyrin methylase 1	mix_CICLE_v10031826mg	VIT_13s0064g01470	1.46	0.71	2.35	0.92
		mix_LOC100852901.1.4	VIT_13s0064g01470	1.52	0.69	2.39	0.99
		mix_LOC100852901.2.4	VIT_13s0064g01470	1.49	0	2.36	0.92
		mix_LOC100852901.3.4	VIT_13s0064g01470	1.33	0	2.31	0.84
		mix_LOC100852901.4.4	VIT_13s0064g01470	1.54	0	2.37	0.9

The gene names have been associated to each contig according to the CRIBI annotation v1. Log Fold Change (LFC) values are indicated for each contig and condition. When the differential expression between HN and LN roots at a given time post-treatment was significant [|LFC|>1 and False Discovery Rate (FDR)<0.01], the numbers are highlighted in grey. Bold numbers indicate when genes were found to be up-regulated in the HN root side compared with LN side. Italicized numbers indicate when genes were down-regulated.

The gene names have been associated to each contig according to the CRIBI annotation v1. Log Fold Change (LFC) values are indicated for each contig and condition. When the differential expression between HN and Ln class="Chemical">N roots at a given time post-treatment was significant [|LFC|>1 and False Discovery Rate (FDR)<0.01], the numbers are highlighted in grey. Bold numbers indicate when genes were found to be up-regulated in the HN root side compared with LN side. Italicized numbers indicate when genes were down-regulated.

The expression profiles of seven genes with contrasting patterns were verified by qPCR. NR, n class="Gene">NIR, NPF6.3, NRT2.4a, NRT2.4b, NRT3, and GS2 were differentially expressed in both combinations, following the same pattern of expression ratio between the HN and LN conditions (Fig. 7). The ratio of transcript abundance between HN and LN roots was 0.82- to 6.73-fold depending on the gene, and was higher in CS/RGM (average ratio=3.04) than in CS/1103P (average ratio=2.12). Pairwise Pearson’s correlations between LFC values from edgeR analysis on RNA-seq data and those from qPCR data showed a strong correlation, indicating that the two methods gave very similar profiles (Fig. 7). Interestingly, the pronounced response observed for CS/RGM was found to be correlated with a higher difference of nitrate content between each side of the root system. Nitrate quantification performed on root samples showed that the nitrate content was significantly higher in HN roots than LN roots at 24 hpt in both combinations. Moreover, the ratio between each N condition was higher for CS/RGM than for CS/1103P (Fig. 8).

Fig. 7.

Validation of expression profiles of N-related genes by qPCR. Left panel, LFC values of normalized gene expression quantified using qPCR. Transcript levels are normalized to the reference genes EF1γ (VIT_12s0035g01130) and GAPDH (VIT_17s0000g10430) and relative to the control condition LN at 0 hpt for each combination. CS/1103P is represented in black and CS/RGM in grey. The HN condition is indicated by solid lines and the LN condition by dashed lines. Data are presented as mean±SE (n=3 biological replicates). Significant differences between conditions at each time point are indicated as *P<0.05 and **P<0.01 (Student’s t-test). Right panel. Pearson’s correlations of the LFC values obtained in qPCR (x-axis) and in RNA-seq with edgeR (y-axis) relative to the control condition LN at 0 hpt for each combination. The correlation coefficient and P-value (Bonferroni adjusted) are presented in the grey boxes. (This figure is available in colour at JXB online.)

Fig. 8.

Nitrate concentrations in root samples of CS/1103P and CS/RGM at 24 hpt. Data are presented as mean±SE (n=3 biological replicates). Statistical analyses of the rootstock genotype (genotype), nitrogen supply (condition), and their interaction (genotype × condition) effects were performed by analysis of variance: *P<0.05, **P<0.01, ***P<0.001. Bars are shaded depending on condition, with black corresponding to HN roots and grey to LN roots.

Validation of expression profiles of N-related genes by qPCR. Left panel, LFC values of normalized gene expression quantified using qPCR. Transcript levels are normalized to the reference genes EF1γ (VIT_12s0035g01130) and GAPDH (VIT_17s0000g10430) and relative to the control condition Ln class="Chemical">N at 0 hpt for each combination. CS/1103P is represented in black and CS/RGM in grey. The HN condition is indicated by solid lines and the LN condition by dashed lines. Data are presented as mean±SE (n=3 biological replicates). Significant differences between conditions at each time point are indicated as *P<0.05 and **P<0.01 (Student’s t-test). Right panel. Pearson’s correlations of the LFC values obtained in qPCR (x-axis) and in RNA-seq with edgeR (y-axis) relative to the control condition LN at 0 hpt for each combination. The correlation coefficient and P-value (Bonferroni adjusted) are presented in the grey boxes. (This figure is available in colour at JXB online.)

Nitrate concentrations in root samples of n class="Chemical">CS/1103P and CS/RGM at 24 hpt. Data are presented as mean±SE (n=3 biological replicates). Statistical analyses of the rootstock genotype (genotype), nitrogen supply (condition), and their interaction (genotype × condition) effects were performed by analysis of variance: *P<0.05, **P<0.01, ***P<0.001. Bars are shaded depending on condition, with black corresponding to HN roots and grey to LN roots.

Weighted gene co-expression network analysis to identify hub genes in response to N availability

A gene co-expression network analysis was conducted using the WGCNA R package (Langfelder and Horvath, 2008). This approach was used to define clusters of highly correlated genes (modules). Including the grey module (which contained the unconnected genes), 26 modules were found and analysed for their association with each experimental trait (correlation with genotype, n class="Chemical">N condition, or time) (Fig. 9). Then, a kME value (module eigengene-based connectivity) was calculated for each transcript to every module (Supplementary Table S5). For each module eigengene, highly correlated transcripts were filtered with a correlation coefficient >0.80 and a P-value <0.01. The first 100 transcripts were selected to perform an enrichment analysis using Mefisto; the resulting enriched categories are summarized in Supplementary Table S6. The genes corresponding to the transcripts that showed the highest correlation coefficient (>0.9) with the module eigengene were considered as potential hub genes.

Fig. 9.

Module–trait relationships. Experimental traits correspond to each column and their association with each module eigengene (rows) is represented by a correlation coefficient and P-value within parentheses. The colour of the cell indicates the correlation coefficient between the traits: orange indicates a high positive correlation and turquoise a high negative correlation. In the left panel, the number of contigs included in each module is presented in parentheses.

With regard to N treatment, the modules ‘greenyellow’ and ‘purple’ showed the highest correlation value with the n class="Chemical">HN condition and the lowest P-value (Fig. 9). The module ‘greenyellow’ included transcripts up-regulated in the HN roots in both combinations (Fig. 10A) and functional categories were found to be linked to the N response, including nitrate metabolism- and transport-related categories, but also OPP electron transfer and tetrapyrrole synthesis. Most of the nitrate-related genes were found in this module (6PGDH, GS2, NR, NIR, and LBD39) and some of the corresponding transcripts showed a high connectivity with the eigengene (kME >0.8; Supplementary Table S5). TGA1, TGA4, and two-component response regulator ARR9-encoding genes were identified with a kME value >0.9 in this module.

Fig. 10.

Eigengene average expression. For the selected modules (A) ‘greenyellow’ and (B) ‘purple’, samples are represented in columns. Values used for the eigengene average expression in the barplot (upper panel) come from the top 100 associated contigs according to the module membership value (kME). These values are those presented in rows on the heatmaps (orange indicates high expression and purple low expression).

The module ‘purple’ correlated positively with the HN condition and n class="Chemical">CS/RGM (Fig. 10B). The most enriched categories were linked to nitrate transport and ammonia metabolism (Supplementary Table S5). N-related transcripts such as GLT1, NRT2.4a, and NPF6.3 were present in this module. The gene with the highest kME value encoded a BTB/POZ and TAZ domain-containing protein 1-like. Genes encoding TGA1, ERF, and LBD38/39 transcription factors, as well as CBL-interacting protein kinases, were also found as potential hub genes in this module.

Three modules, ‘lightyellow’, ‘pink’, and ‘lightgreen’, were found to be associated with the LN treatment. The first two of these shared the n class="Chemical">ethylene functional category, the module ‘lightyellow’ contained also cell-wall-associated genes, and the module ‘pink’ contained some genes linked to RNA regulation. The module ‘lightgreen’ contained genes related to a hormone class described in the BIN code as ‘abscisic acid synthesis-degradation’, which corresponds to the strigolactone pathway. These genes belong to the CS/RGM cluster 6 (Fig. 6B), indicating that they showed a significant increase in transcript level in the LN roots compared with the HN roots only in this combination.

Modules ‘blue’ and ‘turquoise’ were associated with RGM and 1103P, respectively. They contained genes whose expression pattern was correlated to one genotype regardless of the treatment or the time. In both of these modules, many highly connected genes were related to biotic stress. An LHY-encoding transcript was found to be positively connected to the module ‘blue’ (kME value 0.98) and negatively connected to the module ‘turquoise’.

Conversely, an ELF4-like gene was highly positively connected to module ‘turquoise’ and negatively to module ‘blue’. n class="Gene">LHY and ELF4 are both required for circadian clock function in Arabidopsis (Kikis ). These results suggested that some genes involved in the circadian clock regulation of plant growth might be differentially regulated in the two rootstock genotypes.

Discussion

To understand the molecular mechanisms involved in rootstock control of scion growth, the present study focused on two rootstock genotypes, RGM and 1103P, known to induce contrasting scion growth for the grape cultivar CS. Grafted plants were placed in a heterogeneous n class="Chemical">N availability condition by means of a split-root experiment and root transcriptomic analyses were performed. As this experimental design tended to mimic natural conditions in which N availability is fluctuating, the plants sensed both low and high N content, and gene expression was regulated in response to local and systemic signals (Li ). The transcriptomic responses were investigated at 3 hpt in order to detect the changes due to rapid regulation in response to a local change of N supply, and at 24 hpt, as the response to systemic N signalling integrating whole-plant functioning appears at later time points (Ruffel ; Li ).

To compare the different combinations, the transcriptome was de novo assembled, including RNA-seq reads from both combinations. The quality of the transcriptome was validated and showed a good reliability. As expected, a majority of transcripts corresponded to the n class="Species">V. vinifera reference genome, reinforcing confidence in this de novo assembly approach, even if other works that also focused on grapevine, such as Corso , provided satisfying results through a guided assembly strategy.

Heterogeneous nitrate supply modulates the root transcriptome in the two scion/rootstock combinations in different ways

The differential analysis of transcript levels between HN and Ln class="Chemical">N roots at 3 and 24 hpt highlighted that gene expression was more profoundly impacted by the difference in nitrate availability in CS/RGM than in CS/1103P. Nitrate measurements confirmed these results and suggested a close regulation of nitrate uptake, depending on its availability, for CS/RGM. In addition, a WGCNA analysis defined modules containing genes for which the expression profiles were associated to the different experimental traits and identified genes that could act as a hub in more than one module.

A common set of genes exhibiting the same transcriptional regulation in both genotypes, that is, up- or down-regulation in HN roots compared with Ln class="Chemical">N roots, was identified. These N-responsive genes, which are involved in N uptake and assimilation (e.g. NRT2.4a, NR, and NIR), are specific key actors of the systemic N signalling in A. thaliana and M. truncatula (Li ).

Other known actors responded to the N availability in both genotypes and were found to be included in the ‘greenyellow’ module, which correlated positively with n class="Chemical">HN treatment. The OPP pathway represented the most significantly enriched category. In this category, the genes encoding glucose-6-phosphate dehydrogenase (G6PDH) and 6-phosphogluconate dehydrogenase (6PGDH) were previously shown to be induced by high nitrate levels (Wang ; Lejay ). The OPP electron-transfer-related genes correspond to the ferredoxin precursors essential for NIR. Furthermore, the nitrite reduction catalysed by NIR needs a siroheme as a binding site for the nitrite. UPM1 is involved in the production of siroheme (Tanaka and Tanaka, 2007) via the tetrapyrrole biosynthesis pathway, which corresponds to another category that was significantly enriched in the ‘greenyellow’ module. The genes encoding TGA1 and TGA4 were highly connected to this module eigengene, suggesting that they might be key regulators of the core N response in grapevine roots. In Arabidopsis, the TGA1/TGA4 transcription factors may function in the same nitrate-signalling pathway of NRT2.1/NRT2.2 to regulate lateral root density (Alvarez ).

Interestingly, the module ‘purple’ correlated positively with HN treatment and RGM rootstock genotype. n class="Chemical">N-related genes such as GLT1, NRT2.4a, and NPF6.3 were also present in this module. These genes were significantly up-regulated in RGM in HN roots compared with LN roots at 3 hpt. They were also up-regulated in 1103P, although the level of transcripts (evidenced in both RNA-seq and qPCR data) and the LFC values were lower in CS/1103P than in CS/RGM. Thus, this specific RGM response could be attributed either to genes not responding or genes showing a weaker induction in 1103P.

Although NRT2.4a was up-regulated in n class="Chemical">HN roots in both 1103P and RGM, it was found in this RGM-associated module. In grapevine, among the four genes belonging to the NRT2 family, two encode proteins showing a high similarity to the Arabidopsis AtNRT2.4/AtNRT2.2/AtNRT2.1 protein clade. These nitrate transporters function as two-component high-affinity nitrate influx transporters in Arabidopsis, and form a tetrameric protein complex with AtNAR2.1 (also named AtNRT3.1 or WR3) (Kotur ). The corresponding grapevine genes have been named NRT2.4a and NRT2.4b, and their expression has been previously shown to be affected by nitrate supply (Pii ; Tomasi ). Interestingly, NRT2.4a responded differently depending on the root system side and exhibited the same expression profile as NRT3.1. In Arabidopsis, NRT2.1 is subject to both local and systemic regulation and has been shown to regulate lateral root development (Little ; Remans ), while AtNRT2.4 is induced by nitrate starvation and is involved in the uptake of nitrate in low N conditions (Okamoto ; Kiba ). Our results suggest that grapevine NRT2.4a may be an orthologue of AtNRT2.1 and may play an important role in the response to heterogeneous nitrate availability.

The top hub gene of the module ‘purple’ encodes a BTB/POZ and TAZ domain-containing protein 1-like. Using a systems biology approach, Araus recently identified BT2, a BTB/POZ and TAZ domain protein-encoding gene, as the most central and connected gene in an n class="Species">Arabidopsis N use efficiency (NUE) network. They suggested that BT gene family members act as negative regulators of nitrate uptake and NUE in plants. In our study, this gene appeared as a hub gene in a module correlated with the rootstock genotype that confers the lower scion growth (i.e. RGM). The present results cannot demonstrate a role for BT genes in grapevine in response to N availability; however, these genes are interesting candidates for future research on NUE control by the rootstock.

The transcription factor gene TCP20 plays a key role in the systemic signalling pathway that directs root foraging in heterogeneous n class="Chemical">N availability (Guan ). Interestingly, a gene encoding a homologue of AtTCP20 was identified in the module ‘blue’. This module was uncorrelated to nitrate treatment, which is consistent with the fact that in Arabidopsis this gene is not nitrate inducible (Wang ; Wang ). However, the module ‘blue’ was strictly positively correlated with RGM, suggesting that this transcription factor represents an important actor for further characterization of the influence of rootstock genotype on systemic signalling in response to N availability.

Transcriptional regulation of the ethylene and strigolactone pathway genes is rootstock dependent

Genes belonging to the ethylene-related functional category were significantly over-represented in the DEGs between Ln class="Chemical">N and HN roots in both genotypes, but with a time-dependent expression pattern. Most of these genes encode members of the apetala2/ethylene response factor (AP2/ERF-TF) super family. A recent study showed that the expression of these genes was regulated by N availability (Zhao ). The interaction between ethylene and N affects several physiological processes, including root architecture (Khan ). In addition, in Arabidopsis, ethylene regulates some nitrate transporters, such as NPF6.3 and NRT2.1 (Tian ). The WGCNA approach confirmed that ERF genes were found in modules associated with the LN treatment. Interestingly, this approach also highlighted several genes putatively involved in strigolactone biosynthesis in a module that was positively associated with LN treatment. These genes were significantly up-regulated in LN roots compared with HN roots only in CS/RGM. Strigolactones belong to a recently identified group of plant hormones known to regulate plant development and architecture in response to the environment, particularly phosphorus and N availability (Xie ; Waldie ; Lopez-Obando ; Pandey ). Interestingly, their expression was modulated in the rootstock that is known to confer the lower scion vigour (RGM), particularly when plants are facing limiting N availability (Lecourt ). These results highlighted a contrasting N response between the rootstocks, implying different regulation between hormonal pathways and N content.

Various circadian-clock-related genes are highly connected to modules associated with each rootstock genotype

The circadian clock regulates many aspects of plant biology, including primary metabolism, hormone signalling, responses to biotic and abiotic stresses, and plant development. The mechanisms of the plant circadian clock involve multiple interlocking transcriptional feedback loops (Greenham and McClung, 2015). n class="Gene">CCA1 and LHY encode MYB transcription factors, which are components of a negative feedback loop at the centre of the Arabidopsis circadian clock. CCA1/LHY may be involved in a morning loop, and ELF3 and ELF4 are members of an evening loop (McClung, 2014). It has been suggested that in roots the circadian clock orchestrates diurnal carbon allocation and growth (Yazdanbakhsh ). Here, we showed that some components of the circadian transcriptional feedback loops followed different expression patterns in the roots of the two rootstocks, and that two of these genes (ELF4 and LHY) may be hub genes conversely connected to modules associated with each genotype. This result suggests that circadian-clock-related genes, and thus control of metabolism and development, are differentially regulated in the two genotypes.

Conclusion

This comparative global transcriptomic analysis in roots of grafted grapevines showed that the two studied rootstocks responded differently to heterogeneous N availability. These results highlighted that, in addition to a core n class="Chemical">N response common to both genotypes, the transcriptomic response was, unexpectedly, enhanced in the rootstock conferring the lower scion vigour (RGM). The difference between the two genotypes was even more pronounced at 24 hpt and involved genes related to hormonal pathways such as the ethylene and strigolactone pathways. These results suggest the involvement of N nutrition in the control of grafted scion vigour, since two rootstocks differing in this developmental trait presented a genotype-specific transcriptomic signature. They also showed a difference in the timeframe at which a systemic response occurred in response to N heterogeneity; this difference may originate from the integration of different metabolic signals from the whole plant. In addition to the differential gene expression analysis, a WGCNA analysis provided insights into the gene networks and the hub genes connected to genotype-specific modules. These hub genes are interesting candidates for further investigations of the control of scion vigour by the rootstock and of the control of NUE in grapevine.

Supplementary data

Supplementary data are available at JXB online.

Fig. S1. Length distribution of de novo assembled contigs.

Fig. S2. Top-hit species distribution of the merged final transcriptome.

Fig. S3. Multidimensional scaling plot of count data from the different RNA-seq samples

Fig. S4. The few DEGs responding to N availability in CS/1103P involve mainly N-related genes.

Fig. S5. Functional categories distribution within DEGs.

Table S1. Read number before (raw) and after trimming (cleaned).

Table S2. Summary table of the de novo assemblies.

Table S3. List of the primers used for qPCR experiments.

Table S4. Pairwise comparisons resulting from differential expression analysis performed using edgeR.

Table S5. Gene module membership.

Table S6. BIN code enrichment for each WGCNA module.

Click here for additional data file.

Click here for additional data file.