Clonally propagated plants from the GF677 rootstock (Prunus dulcis Miller ×P. persica), which was selected for its high level of drought tolerance, were acquired from a commercial nursery (Agromillora Iberia, S.L., Barcelona, Spain). The rootstocks were grown for 2 weeks in 300 cm³ pots containing a peat substrate, and then they were micrografted with P. persica var. Catherina.

Subsequently, 30 representative plants were transplanted into 15 L containers with TKS-1, a 1:1 ratio of sand to peat substrate (Floragard, Oldenburg, Germany) and 2 g kg^-1 Osmocote 14-13-13 (The Scotts Company LLC, Maryville, OH, USA). The plants were grown in an experimental greenhouse in Zaragoza, Spain (41°43′ N, 0°48′ W) under controlled environmental conditions (23°C day/18°C night, 14 h light/10 h dark photoperiod) for 21 days before the start of the experiment. During this period (April 2011 to May 2011) the plants were watered daily until runoff was visible.

The drought stress experiment started on May 14 and continued for 16 days. The 30 plants were randomly separated into two groups: well-watered plants (the control plants) and water-deprived plants (the drought-stressed plants). The control plants were watered daily to field capacity while the stressed plants were watered with 80% of the quantity of water that had evapotranspired the previous day (Marcelis et al., 2007). The soil water content was measured using time domain reflectometry (TDR), with 20 cm-long probes inserted vertically into the containers, as described in Moret-Fernández et al. (2012). The soil water content and physiological parameters, namely, stem water potential (Ψs), stomatal conductance (g_s), photosynthetic rate (A_N), and intercellular CO₂ concentration (Ci), were recorded for both groups. These measurements were taken on days 0, 7, 12, and 16 after the start of the experiment, on clear days between 10:00 and 12:00, as reported in Jiménez et al. (2013). On day 16, leaf and root samples of three randomly selected biological replicates were collected from both the control and drought-stressed plants (12 samples in total). The plant tissues were then frozen in liquid nitrogen and stored at -80°C.

The total RNA from the three biological replicates was extracted from the root and leaf samples of the control and water-deprived plants according to the method described by Meisel et al. (2005), which was adapted to mini-preparations (Jiménez et al., 2013). Subsequently, samples were treated with DNase I (Thermo Fisher Scientific, Waltham, MA, USA) to remove the contaminating genomic DNA. RNA integrity and purity were assessed using the NanoDrop 2000 Spectrophotometer (Thermo Fisher Scientific) and the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Only RNA samples with A260/A280 ratios from 1.9 to 2.1, A260/A230 ratios ≥ 2, and RNA integrity values > 8 were used in the subsequent experiments.

Equal amounts of total RNA of each tissue from each experimental group (see Figure 1) were used to construct 12 RNA libraries. Total RNA was submitted to Otogenetics Corporation (Atlanta, GA USA) for RNA-Seq assays. Briefly, 1–2 μg of cDNA was generated using Clontech Smart cDNA kit (Clontech Laboratories, Inc., Mountain View, CA, USA, catalog# 634925) from 100 ng of total RNA. cDNA was fragmented using Covaris (Covaris, Inc., Woburn, MA, USA), profiled using Agilent Bioanalyzer, and subjected to Illumina library preparation using NEBNext reagents (New England Biolabs, Ipswich, MA, USA, catalog# E6040). The quality and quantity and the size distribution of the Illumina libraries were determined using an Agilent Bioanalyzer 2100 (Agilent Technologies). The libraries were then submitted for Illumina HiSeq2000 sequencing according to the standard operation. Paired-end 90–100 nucleotide (nt) reads were generated and checked for data quality using FASTQC (Babraham Institute, Cambridge, UK). FASTQ file was sent to customer for downstream analysis. All the raw reads data were deposited with the European Nucleotide Archive (ENA) as part of project PRJEB12334.

The quality of the raw paired-end reads was monitored using FastQC v0.10.0. Low quality segments, erroneous base calls, or adapter fragments were trimmed using ConDeTri Perl script v5.8.9, which discarded data with poor quality scores (Q < 25) or read lengths < 35 bp (Smeds and Kunstner, 2011). Post-processed reads were mapped to the P. persica Whole Genome Assembly v1.0 using TopHat software v2.0.3. This software was selected as a mapping tool because it can generate a database of splice junctions based on the gene model annotation (Trapnell et al., 2012). The TopHat parameters were set at < 3 mismatches when mapping reads and a maximum of 20 multiple hits for each library. The resulting aligned sequences were used with the Cufflinks suite v1.3.0. In order to generate assembled transcripts for each tissue type and each experimental group. These assemblies were then merged together using the Cuffmerge tool to provide a uniform basis for calculating transcript expression. All the assembled transcripts were deposited with the ENA under accession numbers HADJ01000001–HADJ01060423.

The changes in the relative abundance of transcripts in drought versus control conditions were quantified and normalized to the number of reads per kilobase of transcripts per million mapped reads (RPKM). The gene expression levels were then estimated using the Cuffdiff program from the Cufflinks suite. The statistical significance of the level of differential expression of each gene in the roots and leaves was determined by initially setting the false discovery rate (FDR)-adjusted P-values at 0.05 (which are known as Q-values), using the Benjamini and Hochberg method (Benjamini and Hochberg, 1995). However, in order to reduce the likelihood of false positives, only DEGs with a Q-value < 0.01 were considered for further analysis. The complete workflow of the RNA-seq analysis is provided in Figure 1.

Significant DEGs were annotated using BLASTX by scanning four standard resources: the NCBI nr database², Phytozome11³, the Genome Database of Rosaceae (GDR)⁴, and the Kyoto Encyclopedia of Genes and Genomes (KEGG). The E-value cut-off was set to 1E-5 and a 70% query coverage threshold was used to discard partial/single-domain protein matches. A gene ontology (GO) analysis was performed using standalone Blast2GO v3.2 with the same E-value cut-off. This software assigned GO terms to each DEG to allow their putative functions to be predicted in terms of molecular functions (MFs), biological processes (BPs), and cellular components (CCs) (Gotz et al., 2008). These annotations were enhanced by merging them with InterProScan-assigned GO terms, and then running the annotation augmentation module (Annex). The resulting GO terms were plotted and visualized with the Web Gene Ontology Annotation Plot (WEGO) tool⁵ (Ye et al., 2006).

Gene ontology enrichment was carried out for the DEGs associated with both tissues using the singular enrichment analysis (SEA) function of the web-based tool AgriGO⁶. The input list consisted of the whole set of DEGs, while the annotation of Peach Genome (v1.0 Joint Genome Institute) was used as a pre-computed background. Overrepresented terms in the three main categories (BP, MC, and CC) were filtered using Fisher’s exact test and the Benjamini-Hochberg multiple testing correction (Q-value < 0.05).

In order to confirm the reliability and accuracy of the RNA-seq analysis, RT-qPCR was performed on a set of 34 root and leaf genes selected for their putative drought-related functions, including 16 upregulated genes (Log₂FC > 2), seven unchanged genes (|Log₂FC| < 2), and 10 downregulated genes (Log₂FC < -2). In these expressions, FC is the fold change ratio between the drought-stressed and control group RPKM expression (Supplementary Table S1).

The total RNA samples (1.5 μg) were treated with DNase I to remove the contaminating genomic DNA. Subsequently, the samples were reverse transcribed using oligo (dT)₁₈ as a primer with RevertAid H Minus First Strand cDNA Synthesis Kit (Thermo Scientific). The RT-qPCR reactions were performed with the 7500 Fast Real-Time PCR System v2.0.1 (Applied Biosystem by Life Technologies, Grand Island, NY, USA) using bi-technical replicates and tetra-biological replicates for each tissue-experimental group (two of the biological replicates were from the same plants used in the RNA-seq analysis). The reactions were performed using 10 μl of SYBR Green Master Mix (Kapa Biosystems, Boston, MA, USA), 1 μl of each primer (making a total of 4 μM), and 5 μl of diluted cDNA in a final volume of 20 μl. Control cDNA and control primer were included for each run. The primers were designed using NCBI primer-BLAST software⁷ according to the following criteria: primer size of 18–22 bp, GC content between 40 and 60%, amplicon size of 90–160 bp, and annealing temperatures from 57 to 62°C. Moreover, the primers were aligned to the target gene sequence using BioEdit software v7.2 to ensure specific annealing. The likely secondary structures were also assessed to avoid hairpins and primer dimers. A BLASTN scan of the theoretical amplicon was carried out to test the homology to the target genes. Finally, each of the products underwent gel electrophoresis to confirm the presence of a single amplicon of the expected length. The primer sequences and features are listed in Supplementary Table S1. The efficiencies and quantification cycle (Cq) for each gene were calculated using the LinRegPCR program (Ruijter et al., 2009). Gene expression measurements were determined using the gene expression difference (GED) formula (Schefe et al., 2006). Actin 2 and AGL-26 were used as reference genes for data normalization. The relative expression was calculated with respect to the GF677 rootstock control group. A correlation analysis between the levels of gene expression according to the RNA-seq and RT-qPCR analyses was performed using SPSS v23.0.

After a period of drought stress of 16 days, the main visible effects on the plants were wilting and slight defoliation due to the decreased turgor pressure and the shrinkage of the leaf cells. In order to verify that the visual symptoms were indicative of exposure to water deficit conditions, the soil water content and standard physiological parameters of drought-induced effects on leaves were measured (Table 1). The soil water content dropped remarkably from 26.63% in the control plants to 10.69% in the drought-stressed plants, which indicated the presence of decreased turgor pressure and therefore explained the wilting. A decrease was also observed for the stem water potential (Ψs) in response to the reduction in soil water content, further confirming that the plants experienced drought stress. In addition, as leaf water status is considered to be a reliable indicator of plant water balance, stomatal conductance (g_s), intercellular CO₂ concentration (Ci), and the net photosynthetic rate (A_N) were measured. The results revealed that the drought-stressed plants exhibited a lower stomatal conductance compared to the control plants. This reduction led to a significant decline in intracellular CO₂ concentrations, which decreased from 277.30 μmol CO₂ mol^-1 to 261.27 μmol CO₂ mol^-1. Taken together, these results explain the observed slowdown of the photosynthetic machinery (A_N) as a result of the drought conditions.

As roots are the first organs to be exposed to drought, and leaves are the first to sense water loss, both tissues were sampled from the control and the drought-stressed plants and used for transcriptome analysis to obtain an overview of the responses of P. persica during water deprivation. Three biological replicates were processed in order to construct 12 RNA libraries (Figure 1). A total of 315 M paired 100 bp reads were generated, ranging from 20.81 to 61.30 M raw reads per library. Among those, more than 188 M (59.85%) high-quality sequences were retained after pre-processing and filtering out reads containing adaptors, short reads (<35 bp), and reads with low quality scores (Q < 25). The remaining 136 M single and paired reads (72.34% of all the high-quality sequences) were mapped to the Prunus reference genome (from 10.18 to 14.78 M reads per library), which showed that the quality of these mapped genes was sufficient to conduct the subsequent analysis. A summary of the raw data generated, and the trimmed and mapped reads, is summarized in Table 2.

Consequently, the aligned sequence reads were used for reference-guided assembly and thereafter merged using the Cufflinks–Cuffmerge workflow. Overall, 34,559 and 26,062 transcript isoforms were obtained from roots and leaves, respectively, which correspond to 22,079 and 17,854 genes, respectively (Figure 1).

In order to explore the transcriptional response to drought stress, genes in both roots and leaves were tested for differential expression between the control and drought conditions. Expression levels for each gene were calculated and normalized to RPKM values. Initially, the multiple testing correction involved a Q-value < 0.05 and a total of 1,171 genes were found to exhibit differential expression (813 in the roots and 358 in the leaves). Subsequently, a more stringent Q-value was applied (< 0.01) in order to identify the most reliable DEGs. In this analysis, 500 DEGs were identified in the roots and 236 in the leaves. The distribution of these genes is provided in Figure 2. As illustrated, in drought stress conditions, the number of downregulated genes was slightly higher than the number of upregulated genes. Furthermore, the analysis revealed that there were approximately twice the number of DEGs in the roots than in the leaves, confirming the expectations that (i) the root is the first organ that senses and is affected by drought stress and (ii) roots respond faster to stress than leaves, undergoing more complex gene regulation during water deprivation. These results further highlight the key role of rootstocks as stress buffers.

In order to assign putative functions to the DEGs, all identified transcripts with a Q-value < 0.01 were subjected to Blast2GO annotation. After enhancing the annotation by running InterProScan and Annex, 1,904 root annotations and 815 leaf annotations were obtained. Out of the 500 DEGs in the roots, 361 (72.2%) were successfully annotated; however, 52 (10.4%) DEGs did not match any of the sequences characterized in the databases, which may indicate the presence of novel genes. Among the root sequences, 53 (14.7%) were either hypothetical or uncharacterized proteins. Regarding the leaves, 169 out of 236 sequences (71.6%) were annotated, including 19 (11.2%) with hypothetical or uncharacterized functions, while 34 (14.4%) sequences had no significant hits. A summary is presented in Supplementary Figure S1. Regarding the matched species, the majority of the highest-scoring hits were from P. persica (70.49% in the root and 79.69% in the leaves) and P. mume (24.63% in the roots and 16.75% in the leaves), which belongs to the family Rosaceae. These results confirm the quality of our data and the assembly process. A graph displaying the species distribution and the top BLASTX hits is provided in Supplementary Figure S2.

The entire set of DEGs was subjected to GO analysis in order to achieve a broader functional characterization. As a result, 500 DEGs in the roots and 236 in the leaves were classified into 56 subcategories within three main categories (BP, MF, and CC). In total, 283 DEGs in the roots and 137 in the leaves were associated with BP terms, 298 root DEGs and 139 leaf DEGs were associated with MF terms, and 211 root and 108 leaf DEGs were annotated with CC terms (Figure 3). Note that in many cases the same sequence can be assigned to more than one category. In both tissues, the most represented BP subcategories were “metabolic process,” followed by “cellular process” and “response to stimulus” (Supplementary Figure S3). As for MF, the major subcategories were “binding” and “catalytic activity”. Finally, among the CC terms, “membrane” was the most dominant subcategory for the roots, followed by “cell” and “cell part”, which were associated with the leaves. These results, with a comprehensive list of GO subcategories, are plotted in Figure 3.

The annotations within the BP category were the most informative, as they are easier to interpret in the context of drought responses. The list of DEGs annotated with “metabolic process” is an important resource for the identification of novel genes involved in drought acclimation. In addition, Figure 3 shows that “locomotion” (GO: 0040011), “hormone metabolic process” (GO: 0042445), “detection of stimulus” (GO: 0051606), and “cell killing” (GO: 0001906) were exclusively associated with root DEGs, thus highlighting the essential role of roots in plants’ responses to drought. Moreover, DEGs involved in “response to stimulus” (GO: 0006950) seem to play a pivotal role in drought sensing and the responses.

The GO enrichment analysis was performed using a Q-value < 0.05. The results revealed that, in the analysis of the root tissue, there were significant differences in the 26 GO terms between the DEGs and the genome reference, while only three molecular GO terms were enriched in the analysis of the leaves (“heme”, “iron”, and “tetrapyrrole binding” data not shown). In the analysis of the roots, the enriched GO terms were related to BP and MF (12 and 14 GO terms, respectively, as shown in Supplementary Figure S4). The most significant BP terms, such as “responses to stimulus” (GO: 0050896), “responses to stress” (GO: 0006950), and “biotic stimulus” (GO: 0009607), are illustrated in Supplementary Figure S5.

A heat map of DEGs involved in the “response to stimulus” is shown in Figure 4 (see list in Supplementary Table S2). These genes were clustered into five clades according to their expression patterns. The genes in clusters C1, C4, and C5 had higher levels of expression in the control plants than in the drought-stressed plants; they mainly encode peroxidases and proteins related to “responses to biotic stimulus”, such as major allergen proteins, which might indicate that drought turns off this machinery. The remaining clusters (C2 and C3) comprised genes with higher levels of expression in the drought-stressed plants than in the control plants and they included kinases, transcription factors (TFs), and genes related to phosphate starvation.

After annotation, the DEGs in both tissues were classified into two major categories: regulatory genes (genes implicated in signaling and transcriptional regulation) and functional genes (genes that encode proteins that are directly involved in cell protection and damage repair). These genes are described in the next sections and they are further detailed in Supplementary Tables S3A,B.

In order to look into the pathways that the DEGs were involved in, KEGG analysis was carried out. The major classes identified in roots and leaves are included in Supplementary Figure S6, which were Metabolism of amino acid in roots (31 DEGs, Supplementary Figure S6A) and Xenobiotics biodegradation and metabolism in both tissues (22 DEGs in the roots and 25 DEGs in the leaves). The complete set of matched pathways is summarized in Supplementary Table S4. Of the 500 DEGs in the roots, 152 had significant matches in the KEGG database (128 enzymes) and they were classified into 63 pathways (Supplementary Table S4A). In the leaves, 88 genes were assigned to 36 KEGG pathways and associated with 50 enzymes (Supplementary Table S4B). The major pathways identified in the roots were phenylpropanoid biosynthesis (nine DEGs and three enzymes) and aminobenzoate degradation (nine DEGs and four enzymes). In the leaves, glutathione metabolism (10 DEGs and four enzymes), aminobenzoate degradation (10 DEGs and two enzymes), and drug metabolism associated with cytochrome P450 (six DEGs and one enzyme) had the highest levels of differential expression.

In order to further confirm the accuracy of the RNA-seq expression estimates, a total of 34 candidate genes were selected for RT-qPCR validation according to their RPKM transcript abundance and Log₂FC. As illustrated in Figure 6, the expression values of the selected DEGs in both tissues significantly correlated with the RPKM values, with the exception of the chloroplastic ribulose bisphosphate carboxylase small chain (RBCS) gene, which may be a result of the unstable expression of this chloroplastic gene. The correlation between the RNA-seq and RT-qPCR measurements was evaluated using linear regression, based on the following equation: RT-qPCR value = b (RNA-Seq value) + a (Figure 7). Interestingly, the linear regression analysis indicated a highly significant correlation between the methods, indicating a general agreement regarding the transcript abundance determined by both methodologies (r = 0.89 and r = 0.95 for root and leaf DEGs, respectively). In conclusion, the obtained results confirm the reliability of the transcriptomic profiling data estimated from RNA-seq data.

As sessile organisms, plants are unable to escape when environmental conditions become unfavorable. Nevertheless, they can successfully deploy complex physiological and molecular strategies to cope with environmental stresses. In this study, the physiological measurements confirmed that the plants were effectively subjected to a water deficit, and the plants elicited physiological responses to combat it. In fact, when water availability is limited, plants change their biochemistry in order to be able to retain as much water as possible and increase their chances of survival. One of the earliest responses to minimize water loss is the reduction of stomatal conductance, which leads to a reduction in CO₂ diffusion through the stomata pores. The results concurred with the findings of several other studies that showed that water scarcity significantly reduces the rate of photosynthesis (Jiménez et al., 2013) by affecting the CO₂ balance and stomatal status (Rahmati et al., 2015).

In order to investigate the dynamic changes in gene expression in the roots of GF677 rootstock and the leaves of Catherina cultivar budded together, RNA-seq was employed using the Illumina platform. Surprisingly, after quality control, the total number of clean reads generated from the RNA libraries of the drought-stressed plants was lower than the number for the control plants (Table 2). This is in contrast with previous studies, which have reported an activated plant transcriptome in response to drought (Tang et al., 2013; Garg et al., 2016). We thus hypothesize that a drought period of 16 days may be too short to drive the full upregulation of the P. persica genome, especially as the rootstock used was selected for its tolerance to drought (Jiménez et al., 2013).

The most frequent BLASTX top hits in our sequence homology searches were from P. persica and P. mume. These similarities highlight the quality of the assembly process. The annotation of the DEGs revealed a considerable number of hypothetical or uncharacterized protein functions, and some of these genes had large changes in expression. These genes could provide a good starting point for further experimental characterization. Generally, proteins with unknown functions are widespread across species, even in model plants. Indeed, in Arabidopsis thaliana, 13% of the genes encode proteins with unknown functions (Luhua et al., 2013). In spite of this, such genes are potentially interesting as they may encode proteins that would be valuable to breeders. Elucidating their biological roles in Prunus spp. is thus an important challenge, which we aim to achieve in future studies.

After GO annotation, the DEGs were labeled with 56 terms within three main categories (BP, MF, and CC). The most dominant terms, illustrated in Supplementary Figure S3, concur with the findings of previous research, confirming their universal involvement in the response to drought stress conditions (Dong et al., 2014; Li et al., 2016). In addition, some GO terms such as “locomotion”, “hormone metabolic process”, “detection of stimulus”, and “cell killing” were exclusively associated with root DEGs (Figure 3), indicating that despite being sessile, P. persica roots exhibit dynamic changes in architecture in response to water scarcity. An induced phospholipase D (PLD), associated with the regulation of cell migration and root development, is likely to be involved in locomotion. This observation likely reveals that, under conditions of drought stress, roots move downward in order to find water and escape from the harmful external factors. In fact, PLD hydrolyzes lipids, which results in the formation of phosphatidic acid (PA), a compound that is responsible of inducing cell proliferation and primary root growth (McLoughlin and Testerink, 2013). Functions related to the detection of stimuli were found to be increased in the roots during the drought responses of P. persica, confirming that these organs are responsible for sensing water deprivation. According to previous studies, the main steps for handling any type of abiotic stress are signal perception, signal transduction, and expression of stress-inducible genes. Thus, we propose that root cells first perceive drought through sensors located in the cell wall and/or membrane and then they convey the signals to other organs. Indeed, we have found that “membrane” and “cell” were the most dominant GO terms in the CC category, particularly regarding root DEGs. The “hormone metabolic processes” (GO: 0042445) associated with the root DEGs will be further discussed in Section “Signaling and Regulatory Proteins.”

We observed that the initiation of drought stress triggered a wide range of responses, which implies that there are many genes and mechanisms involved in drought tolerance in P. persica. According to their associated proteins, we classified the DEGs as signaling and regulatory DEGs or functional DEGs.