SBCC073 is a Spanish landrace-derived line, developed by three generations of head-row multiplication starting from a single seed of the original landrace, with no phenotypic selection. It is fully homozygous (http://www.eead.csic.es/EEAD/barley/core.php?var=73). Scarlett is a German spring malting cultivar. Seed was provided by Saatzucht Breun GmbH. SBCC073 is an intermediate growth habit type, according to Szucs et al. (2007), and it needs just 15 days of vernalization to promote timely flowering. Seeds of SBCC073 and Scarlett were sown and seedlings were allowed to grow for 1 week. Then, they were vernalized for 24 days, in order to synchronize the development of the two genotypes. At the end of the vernalization period, plants at the 3-leaf stage were transferred to 28.0 × 20.8 cm (height × diameter) black plastic pots (one seedling per pot) with standard substrate made of peat, fine sand and perlite Europerl B-10 (Europerlita Española SA, Barcelona, Spain), from a mix with 46 kg, 150 kg, and 50 L, respectively. Two series of pots were placed in a greenhouse (natural photoperiod, controlled maximum temperature 28°C, average daily temperature 25 ± 2°C during the day and 21 ± 3°C at night) and in a growth chamber (16 h light/8 h dark, 21°C daytime/18°C night temperature). Additional pots filled only with substrate were used to estimate dry weight and field capacity (FC). Soluble fertilizer was provided with irrigation. Plants were treated with fungicide (Triadimenol 25%) to prevent powdery mildew build-up.

Plants were subjected to a control and two different stress treatments: severe drought, and mild drought combined with heat. Control conditions and the severe drought treatment were imposed in a growth chamber. There, plants remained well-watered (70% FC) for 30 days. Then, water application was gradually reduced, to resemble a slow drying soil, based on weight of each pot relative to the estimated FC. Once the target fraction of FC was reached, the pots were watered to keep such weight constant. Treatment levels in the growth chamber were 70% (control) and 20% FC (severe drought). At the sampling date, 30 days after the beginning of the drought treatment, all plants had been at the target fraction of FC for at least 14 days. Plants in the greenhouse were subjected to heat stress right after transplant, for 61 days, from Zadoks growth stage (GS) 20 (Zadoks et al., 1974) up to GS 73, on average. Out of these, the first 30 days the plants were well-watered (80% FC). The last 31 days, the drought treatment was applied gradually, as in the growth chamber, targeting an intermediate drought stress level of 50% FC. During all 61 days, day temperature was kept under 30° with a cooling system, whereas temperature was not controlled at night. Average temperature over all period was 23.6°C (daily range 20.9°–26.3°C), average maximum temperature was 27.6°C, and average minimum temperature was 18.8°C. The pace of accumulation of growing degree days at this treatment clearly exceeded rates of accumulation found in typical field winter and spring sowings, considering the same period of crop development (Figure S1). During the experiment, average temperature exceeded 25°C half of the days. During the treatments, all pots were weighted, watered, rotated and their positions swapped every 2 days.

Several traits were recorded 60 days after transplant. Leaf water potential (LWP) in leaves was measured at noon using a Scholander chamber (SF-PRES-70, Solfranc Tecnologías SL, Vila-Seca, Spain). Stomatal conductance (SCo) was measured, starting at 9 am, using a leaf porometer (Decagon Devices, Pullman, WA, USA). Relative water content (RWC) was also estimated, as described in Talame et al. (2007). For each plant, three independent measurements were taken for LWP, SCo, and RWC. In addition, tiller number per plant (TN) was counted. All measurements were taken at two biological replicates.

Two tissues, young inflorescences and leaves were sampled at 60 days after transplant. Flag leaves were collected. In some cases, flag leaves of main tillers were too small to ensure having enough tissue for RNA extraction. In these cases, last fully expanded leaves of secondary tillers were also harvested, as a backup. Fresh material was harvested and frozen in liquid N₂ before RNA extraction with the NucleoSpin® RNA Plant kit (Macherey-Nagel, Düren, Germany). RNA quality was assessed with a NanoDrop 2000 spectrophotometer (Thermo Scientific, Wilmington, DE, USA) and with Bioanalyzer 2100 hardware (Agilent, Santa Clara, CA, USA; average RIN: 6.7 for leaves, 8.1 for flowers). Barcoded cDNA libraries were prepared at CNAG (Barcelona, Spain) following Illumina TruSeq standard procedures, and eventually sequenced in an Illumina HiSeq2000 sequencer, using a full flow-cell, 4 samples per lane, to produce 2 × 101 bp paired-end reads. The whole dataset consisted of 2 biological replicates from greenhouse plants (2 tissues × 2 replicates × 2 genotypes), 2 biological replicates of developing inflorescences and 3 biological replicates of leaves from plants subjected to drought and well irrigated plants in the growth chamber (5 × 2 genotypes × 2 treatments). Leaf samples used for sequencing finally came from flag leaves and last fully expanded leaves in equal numbers. The results of both types of leaves were indistinguishable, and thus were treated as the same tissue.

Raw reads were sequentially processed with FASTQC v0.10.0 (Andrews, 2010) and Trimmomatic v0.22 (Bolger et al., 2014), discarding stretches of mean Phred score <28 and cropping the first nucleotides to ensure a per-position A, C, G, T frequency near 0.25. Only reads of length ≥80 nucleotides were kept for further analysis. Surviving reads were error-corrected with Musket v1.0.6 (Liu et al., 2013) and default parameters. Then, reads were assembled following two different procedures, de novo and reference-guided.

De novo assemblies were obtained using Trinity r2013-02-25 recommended procedures (Haas et al., 2013). First, reads from sample replicates were pooled together and in silico normalized, to a maximum coverage of 30. This procedure was repeated with the resulting read sets to obtain, for each genotype, a final set of normalized reads. These were used for de novo assembly of SBCC073 and Scarlett transcriptomes.

A reference-guided assembly (RGA) was generated with the Tuxedo pipeline (Trapnell et al., 2012). First, clean reads were mapped to the IBGSC cv. Morex assembly (Mayer et al., 2012) with Tophat2 (v2.0.9; –b2-very-sensitive, –b2-scor-min C, − 28,0 –read-mismatches 4–read-gap-length 12 –read-edit-dist 12-G 21Aug12_Transcript_and_CDS_structure.gff). This mapping procedure was performed in two steps, a first one to exclude reads with multiple mappings to the whole reference assembly (-M, -g 1, –no-discordant) and a second one to identify reads mapping unambiguously to gene coding loci (-g 2, –no-discordant, –no-mixed). Mappings were used as input for Cufflinks (v2.2.1). Individual assemblies were merged with the reference Morex assembly with Cuffmerge.

Clean reads were mapped back to the de novo transcriptomes using Trinity script alignReads.pl with Bowtie (Langmead et al., 2009). In addition, the newly assembled isoforms were mapped to Morex, Bowman, Barke WGS (Whole Genome Shotgun) assemblies (Mayer et al., 2012) and Haruna Nijo flcDNAs (Matsumoto et al., 2011) with the script bmaux_align_fasta from the Barleymap package (Cantalapiedra et al., 2015) (hierarchical = yes query-mode = cdna thres-id = 98 thres-cov = 10), keeping together sequences matching the same reference sequence. Sequences in each of these groups were clustered with WCD-express v0.6.3 (Hazelhurst and Liptak, 2011) using threshold = 24, which is equivalent to a 98% identity cut-off.

Presence of these isoforms in existing references was further confirmed by aligning them iteratively to additional sequence repositories. These were the Haruna Nijo genome assembly (Sato et al., 2016), genome contigs of Chinese Spring wheat (Mayer et al., 2014), barley ESTs from HarvEST assembly 36 (Close et al., 2007), the MIPS repeat database (Nussbaumer et al., 2013), and sequences from Hordeum, Brachypodium, Triticum, Oryza or Aegilops in the nt NCBI database (ftp.ncbi.nlm.nih.gov/blast/db). Alignment to Morex, Bowman and Barke WGS assemblies, and to Haruna Nijo genome and flcDNAs was repeated with a more stringent coverage threshold (thres-cov = 80). Finally, transcripts were scanned for the presence of sequencing vectors by comparison with the EMVec database (ftp://ftp.ebi.ac.uk/pub/databases/emvec/) and as a result 64 sequences were removed.

Gene annotation of assembled contigs was performed with the script transcripts2cdsCPP.pl (-n 50) from GET_HOMOLOGUES-EST (v 04052016, https://github.com/eead-csic-compbio/get_homologues, Contreras-Moreira et al., 2017), which uses Transdecoder (https://transdecoder.github.io/) and blastx alignments to SwissProt proteins to define CDS sequences. Clusters obtained with GET_HOMOLOGUES-EST (get_homologs-est.pl -t 0 -M -S 96 -A –L), requiring percentage sequence identity >96, were used to obtain reciprocal correspondences between transcripts from SBCC073 and Scarlett assemblies. PFAM domains in translated CDS sequences were also annotated (get_homologs-est.pl –D).

Differential expression contrasts were performed for each genotype, tissue and treatment; both for isoforms and genes. For this purpose, we compared three different pipelines.

For the first one, estimation of expression levels of isoforms and genes was done with RSEM v.1.2.11 (Li and Dewey, 2011), using Bowtie2 (Langmead and Salzberg, 2012) and otherwise default parameters. RSEM “expected counts” were used as input for differential expression analyses with the “glm” functions of the R (R Development Core Team, 2008) Bioconductor package edgeR v3.8.6 (Robinson et al., 2010) (false discovery rate function “BH” set to 0.001). A minimum CPM (counts per million) of 0.4, equivalent to around 10 RSEM “expected counts” based on a linear regression (R-square = 1, intercept ~ 0, slope = 25), was required in at least half of the samples to include an isoform or a gene in the analysis.

A second method relied on kallisto v0.42.5 (Bray et al., 2016) to obtain “expected counts” and to generate 100 bootstrap samples for each replicate, followed by test for differential expression with sleuth v.0.28.0 Wald test (Pimentel et al., 2016), using the previously generated bootstrap samples.

For the third method, Cuffquant and Cuffdiff v.2.2.1 (Trapnell et al., 2013) were used to test differential expression, with FDR 0.05, on the RGA transcripts.

Principal component analyses (PCA) of the resulting expression estimates from kallisto were done with the function PCA from R package FactoMineR 1.29 (Lê et al., 2008). Correlation analysis was performed using the R package corrplot 0.73 (Wei and Simko, 2014).

Reference genes for calculating relative expression were either searched in the literature or selected from our RNAseq data. The latter were those with the smallest coefficient of variation of expression values across samples, among isoforms not reported as differentially expressed (DE) by edgeR. DE isoforms to be checked with RT-qPCR were chosen randomly from bins covering the range of edgeR logFC. All the selected DE isoforms had TPM (transcripts per million) >1. Primers for both reference genes and DE isoforms were designed with Primer Express 3.0 (Applied Biosystems). Conservation of the target sequences was checked in both SBCC73 and Scarlett isoforms. Whenever possible, one of the primers of the pair was set over an exon-exon junction and toward the 3′ end.

The same DNase I-treated RNA samples used for RNAseq were utilized for the RT-qPCR assays. First strand cDNA synthesis was made from 2 μg of total RNA to a final volume of 40 μl containing oligo(dT)20 for priming and SuperScript III Reverse Transcriptase (Invitrogen, Cat.No. 18080-044). All the RT-qPCR reactions were performed in an ABI7500 (Applied Biosystems, Foster City, CA, USA) with the following PCR profile: 95°C 10 min pre-denaturation step; 95°C 15 s denaturation and 60°C 50 s annealing (40 cycles), followed by a melting curve 60°–95°C default ramp rate. The efficiency of primers was obtained from calibration curves with 1:5 dilution series and at least 4 points fitted in a linear regression with R-square over 0.99. We used NormFinder (Andersen et al., 2004) to analyze the stability value of the reference genes. Relative change of expression was calculated according to Pfaffl (2001), but using the geometric mean of three reference genes as normalization factor (Vandesompele et al., 2002).

Software CPC (Kong et al., 2007) was used to tag DE isoforms as coding or non-coding, and to obtain Uniref90 best hits. In addition, contained CDS sequences were deduced and PFAM protein domains annotated, as explained earlier for all the isoforms of each transcriptome. GO terms for each DE isoform were obtained with in-house script barleyGO (http://www.eead.csic.es/compbio/soft/barleyGO.tgz). Enrichment tests for PFAM domains and GO terms were performed in R using the Fisher exact test (p < 0.01). For the GO terms, we used the R package topGO (Alexa and Rahnenfuhrer, 2016).

DE isoforms were searched in metabolic pathways databases, including KEGG (Kanehisa et al., 2016), PlantReactome (Tello-Ruiz et al., 2016), and PlantCyc (Plant Metabolic Network, 2016). For KEGG, we obtained the list of genes of Oryza sativa (“osa”), from which we retrieved Orthology identifiers and pathways. DE isoforms were aligned to those genes with blastn (-perc_identity 75 –num_alignments 1), discarding hits with low query coverage in the alignment (“qcovs” < 70). PlantReactome (file “gene_ids_by_pathway_and_species.tab”) was explored with Morex gene identifiers to obtain the pathways involved in differential expression. The gene identifiers were derived from mappings of de novo assemblies to the Morex reference genome from the validation step using the Barleymap package, as explained above. In the case of PlantCyc, we obtained the blast set “plantcyc.fasta” and enzymes annotation (“PMN11_June2016/plantcyc_pathways.20160601”), and used a custom script to match annotated enzymes with blastx (e-0.00001–num_alignments 1), filtering hits with percentage identity ≥75. Enzymes and pathways were grouped in broader categories manually, by merging their textual descriptions in KEGG and PlantCyc.

The literature was surveyed to obtain protein and transcript sequences which had been previously associated with response to water deprivation in barley. These drought-related sequences were aligned with Blast[p|x] to genes from the Haruna Nijo genome assembly, which allowed mapping them to their corresponding DE isoforms from this study.

DE isoforms were clustered based on their TPM values (from kallisto). Distance between each pair of isoforms was calculated with Pearson correlation. This metric was weighted with Euclidean distance, under the hypothesis that isoforms sharing their expression pattern, but differing in magnitude, might have promoters which could be overlooked when clustered together with Pearson correlation only. These distances were used to perform hierarchical clustering (R package hclust, method = “complete”). To declare the final number of clusters, the dendrogram was pruned when 95% of clusters had an internal average distance below 0.001% of the initial average distance of all DE isoforms.

The following procedure was used to recover promoter sequences corresponding to the genes present in the expression clusters. DE isoforms from each cluster were mapped to transcripts from the Morex WGS assembly (Blastn -perc_identity 98). For each cluster containing 10 or more genes, repeat-masked promoter sequences (−1,000, +200 nucleotides around TSS) were retrieved from the RSAT::Plants server (http://plants.rsat.eu, version Hordeum_vulgare.082214v1.29) (Medina-Rivera et al., 2015). As negative controls, promoter sequences were retrieved from randomly generated gene clusters of the same size. Enrichment in GO terms and motif discovery with oligo-analysis and dyad-analysis were performed following the protocol of (Contreras-Moreira et al., 2016). Motif scores within upstream regions of co-expressed genes and their orthologous genes in Brachypodium distachyon reference (v1.0.29) (International Brachypodium Initiative, 2010), were obtained with the program matrix-scan from RSAT::Plants. These scores were also calculated for motifs generated by permutation of the bases of each discovered motif. Therefore, two types of evidences were used to assess the reliability of discovered motifs: (i) their statistical significance compared to the negative controls, and (ii) their matrix-scan scores compared to the scores of permuted motifs. Discovered motifs were annotated by comparison to plant regulatory motifs in the footprintDB repository (Sebastian and Contreras-Moreira, 2014). The highest scoring motif, in terms of footprintDB “Ncor” score, was selected as the best hit. The full report on the promoter analysis, including source code, is available at http://floresta.eead.csic.es/rsat/data/barley_drought_clusters.

Finally, deduced peptide sequences of DE isoforms annotated as transcription factors with iTAK (http://bioinfo.bti.cornell.edu/cgi-bin/itak/index.cgi), were used to predict their putative DNA-binding motifs with footprintDB.

Two different experiments were set up, in which plants were placed in a growth chamber or in a greenhouse. The growth chamber was kept at strictly controlled environmental conditions, whereas the greenhouse underwent a natural photoperiod (August–September, 2012, starting with 14 h 23 min and ending with 11 h 46 min daylight, http://www.fomento.gob.es/salidapuestasol/2012/Zaragoza-2012.txt) and controlled temperature, surpassing typical temperatures expected under field conditions (Figure S1), but more variable irradiation and humidity. Both daytime and night temperatures in the greenhouse were higher than in the growth chamber, whereas relative humidity was similar on average (Figure S2). In both settings, water stress was imposed after initiation of the stem elongation stage. Growth chamber plants were watered in order to stay at 70% field capacity (FC) (controls, C), or instead subjected to reduced irrigation until reaching 20% FC (drought, D). Plants in the greenhouse were irrigated to an intermediate 50% FC (mild drought and heat, MDH). These experiments are outlined in Figure 1.

Daily loss of water, based on the weights of pots, was largest in C plants, intermediate under MDH and lowest under D (Figure S3). The same trend was observed for leaf water potential (LWP), summarized in Figure 1. LWP was proportional to the three imposed water regimes, with plants subjected to drought (D and MDH) showing larger absolute LWP that those well-watered. The largest absolute value corresponded to Scarlett plants under D, in which SBCC073 plants had values comparable to those of both SBCC073 and Scarlett plants under MDH. Likewise, minimum values for stomatal conductance (SCo) were recorded for plants under D (Table 1). However, the largest SCo was found under MDH. Relative water content (RWC) was lowest for plants under D, in both genotypes, whereas under MDH, it was closer to that of C plants in SBCC073, and closer to that of plants under D in Scarlett. Tiller number (TN) was also affected by water deprivation, being larger in C than under D, both in SBCC073 and Scarlett. Under MDH, similarly to the RWC observations, TN was less affected in SBCC073 than in Scarlett.

Sequencing of cDNA libraries, derived from leaf (LF) and young inflorescence (YI) transcripts, yielded 1.18 billion paired-end reads. From this dataset, we assembled separate de novo transcriptomes for Scarlett and SBCC073, as well as a reference-guided assembly (RGA) (Figure S4).

The de novo assemblies yielded similar numbers and lengths of isoforms for both genotypes (Table 2). These sets, with 103,623 genes in SBCC073 and 113,962 in Scarlett, were comparable but larger than the annotated gene sets for the Morex cultivar (Mayer et al., 2012), with 75,258 high and low confidence genes, and with the results from the RGA (75,204 genes). Validation and correction of the de novo isoforms was performed in three stages. First, the clean reads were mapped back to the assembled transcripts, to compute the fraction of well aligned pairs of reads (both reads mapped, correct orientation and insert size), which was near 83% for both cultivars. Second, de novo subcomponents were revised for re-clustering. This requires some explanation. Whereas RGA contigs are isoforms associated to known genes from the reference, de novo assembly generates contigs which are isoforms clustered in so called subcomponents. In some cases, these subcomponents accumulate closely related sequences, for instance from paralogous genes or expressed pseudogenes, which should be separated. Therefore, this second step consisted in re-clustering isoforms from subcomponents to genes, by alignment to annotated references (see Methods), and assigning them to different loci when appropriate. The final number of genes in the de novo assemblies was 112,923 in SBCC073 and 123,582 in Scarlett. Third, the isoforms were matched to a variety of genomic and transcriptomic sequence repositories of barley, wheat and other grasses. In total, 93% of SBCC073 and 87% of Scarlett genes could be confirmed. These sequence comparisons are further illustrated in Figure 2. Note that at least 10% alignment coverage was required in all cases. Further, the alignment against Morex, Barke, Bowman and Haruna Nijo was repeated, with a more strict minimum coverage of 80%. This test confirmed that 88,293 (78% of SBCC073) and 92,713 (75% of Scarlett) genes map with high confidence to previously reported barley sequences.

Clean paired-end reads were mapped back to SBCC073 and Scarlett assemblies, to estimate expression counts for each transcript. These estimates were subsequently used to identify DE tags (genes and isoforms) between stressed treatments and C, for each tissue and genotype. For this purpose, we compared three different pipelines, which rely on different software for each of the two steps: RSEM-edgeR, kallisto-sleuth and Cuffquant-Cuffdiff. In addition, a set of isoforms from YI were randomly chosen to test their expression by RT-qPCR, using genes selected from the literature and from our RNAseq expression data as references (Table S1).

The results of differential expression computed with kallisto-sleuth had the best agreement with those of RT-qPCR (Figure 3). The outcome of the RSEM-edgeR pipeline was comparable to kallisto-sleuth after discarding a few outliers. Moreover, PCA and clustering of samples, using expression estimates from kallisto, showed good correlation between replicates (Figures S5, S6). When the expression estimates, obtained with the three methods, were directly compared, RSEM-edge and kallisto-sleuth showed the best agreement (Figures S7-S9 and Table S2). In order to reduce false positives, final DE tags were obtained from the intersection between those two methods.

Overall, the response differed between genotypes in YI, and between treatments in LF (Figure 4). Under D, we found almost no response in SBCC073, either in YI or LF samples, whereas in Scarlett, YI samples had many DE tags. On the contrary, abundant changes in gene expression were observed under MDH, with the exception of YI from SBCC073, which remained mostly unaltered. Regarding the proportion of up-regulated tags over total DE tags, in LF under MDH it was close to 50%, in both genotypes, whereas in YI from Scarlett plants it increased under D (62.6% in isoforms, 61.4% in genes) and decreased dramatically under MDH. There was high agreement between DE genes and DE isoforms in all contrasts, although some DE genes were different to those found when analyzing isoforms (Table S3). On the other hand, common DE tags between different contrasts were negligible, with the exception of LF under MDH, in which Scarlett and SBCC073 shared a low but sizable fraction (Figure S10).

We explored the cases in which there were several isoforms associated to a single gene differentially affected by the treatments. First, we studied those genes which had more than one isoform in the same contrast (Data Sheet 2). For example, in LF from Scarlett under MDH there were 54 DE isoforms corresponding to 21 genes, ranging from 2 to 6 isoforms per gene. Out of these, just 2 isoforms from a single gene showed opposite directions of change in gene expression. As in that example, for all the other contrasts, the direction of change was the same for most groups of isoforms. The outstanding exception was YI from Scarlett under D, with more than one third of isoforms showing different induction/repression. This suggests that the transcriptional landscape of YI from Scarlett was dramatically affected by D. Second, we assessed the function of isoforms from a single gene, which were differentially expressed in different contrasts (Data Sheet 3). In SBCC073, we found only 4 isoforms, from a single gene annotated as CCA1/LHY (see below), which were all of them induced. Of these, 2 of them occurred in YI and LF under MDH, and 2 in YI only. All 4 isoforms were annotated as coding, putative splice variants. In Scarlett, we found 61 isoforms, associated to 24 genes, 9 of them with isoforms with opposite directions of change in gene expression. For example, a gene annotated as alpha-glucan water dikinase had 9 isoforms of which 3 were induced in YI and 6 were repressed in LF, under MDH. These two groups differed in functional annotation. Moreover, a gene annotated as beta-glucosidase showed 4 DE isoforms under MDH, one induced in LF and 3 repressed in YI. Also, we observed a transcript likely encoding a protein kinase for which one isoform was induced in YI under D and another one was induced in LF under MDH.

Overall gene expression changes (number of DE tags and cumulative logFC from each contrast) were compared with the physiological measurements. Some large correlations were obtained (Table S4), although these results must be considered with care due to the small sample size. For LWP, we found a positive correlation with YI overall logFC of isoforms (r 0.97, p-value 0.03) and number of DE tags (r 0.99, p-value 0.01). SCo exhibited strong positive correlation with gene expression changes in LF (ranges: r 0.95–0.98, p-values 0.05–0.02).

DE isoforms were annotated combining different strategies, as described in Materials and Methods. The main annotation results are detailed in the following sections, whereas the complete results are provided in Data Sheet 4.

As explained in the previous section, just a few isoforms were DE in LF under D. In both genotypes, we found an up-regulated isoform encoding a polyamine oxidase, involved in spermine and spermidine degradation. In addition, an isoform corresponding to a chlorophyll apoprotein from photosystem II was down-regulated in Scarlett. However, this change was not observed in SBCC073, which instead showed induction of transcripts of three proteins, namely ABA/WDS (abscisic acid/water deficit stress) induced protein, ribonuclease T2 and calcineurin-like phosphoesterase. Other DE isoforms were annotated as non-coding or of unknown function.

There were more DE tags in LF under MDH, and involved a more diverse array of gene functions than under D. The same polyamine oxidase induced in LF under D was also observed up-regulated in Scarlett under MDH. Intriguingly, in SBCC073 we found up-regulated a transcript encoding a spermidine synthase.

Some GO terms were enriched in both genotypes, including “phosphorelay signal transduction system,” “pyrimidine-containing compound biosynthesis process,” “response to temperature stimulus,” “response to water deprivation” and “thiamine biosynthetic process” (Data Sheet 5). Other pathways and cellular processes involved in the responses of both genotypes were starch phosphorylation, chorismate biosynthesis, L-ascorbate biosynthesis and recycling, DMNT biosynthesis (a volatile homoterpene), and other proteins involved in protein folding, proteolysis and defense response (Figure 5). We also found in both genotypes up-regulation of isoforms annotated as CCA1/LHY MYB-related TF (Table S5). Moreover, we found another DE gene annotated as MYB-related TF in both genotypes, which is similar to Arabidopsis thaliana TCL2, and an additional uncharacterized MYB-related TF in SBCC073 only. At the same time, down-regulation of other genes related with circadian rhythm was detected, like adagio-like protein 3 and a PRR1 (HvTOC1) transcription regulator. In SBCC073, we found also down-regulation of another circadian clock related gene, annotated as APRR3. Another gene up-regulated in both genotypes was annotated as protein kinase CIPK9. Regarding transporters, repressed transcripts encoding aquaporins were noticed in both genotypes. There were a few other protein domains regulated in both genotypes, most of them repressed.

Differences between genotypes were also seen among DE transcripts in LF under MDH. For instance, in SBCC073 there was enrichment of terms such as “actin filament-based movement,” “ammonium ion metabolic process” and “defense response by cell wall thickening,” while in Scarlett a greater variety of response-related terms were obtained, such as “response to abscisic acid,” “response to bacterium,” “response to ethylene,” “response to hydrogen peroxide” or “response to wounding” (Data Sheet 5). Also, DE isoforms related to glycine betaine biosynthesis and to abscisic acid (ABA) biosynthesis were seen in SBCC073, whereas trehalose biosynthesis was involved in the response of Scarlett LF to MDH (Figure 5). Moreover, isoforms involved in cell wall, epidermis (wax esters) and membrane lipids (glycerophospholipids, ceramide) metabolism were up-regulated in Scarlett but not present among SBCC073 DE isoforms. This was also the case of some defense response metabolic pathways (benzoxazinoids and dhurrin biosynthesis), xanthophylls metabolism, several antioxidation related proteins (like baicalein peroxidase or glutathione S-transferase) or sulfur metabolism related proteins. We also found differences among TFs and protein kinases (PKs) (Table S5). For instance, CIPK17 and a C2C2-Dof TF, whose best SwissProt hit is Arabidopsis protein CDF2, were up-regulated, and an AP2/ERF-AP2 TF (related to Brassica napus BBM2) down-regulated, all in SBCC073. Instead, repression of a TUB TF, similar to O. sativa subsp. japonica TULP7, and induction of both a bZIP TF and a jasmonate ZIM TIFY TF, the latter related to O. sativa subsp. japonica TIFY6B, was noticed in Scarlett. Besides aquaporins, already mentioned, DE isoforms related to transport processes were different between genotypes, being more abundant in Scarlett. These included lipid transfer proteins, phosphate, potassium, triose-phosphate, adenine, vacuolar amino acid and ABC transporters, and a repressed NUCLEAR FUSION DEFECTIVE 4 (NFD4) protein.

In YI, the transcriptional responses were markedly different between genotypes, with only minor responses in plants of genotype SBCC073 under both treatments. Indeed, a single down-regulated transcript was identified in SBCC073 under D, annotated as Pollen Ole e 1 allergen/extension. Under MDH, a repressed isoform was annotated as “non-coding,” whereas four up-regulated isoforms corresponded to CCA1/LHY.

In contrast with what was seen in SBCC073, YI from Scarlett showed abundant gene expression changes. Enriched GO terms found both under D and under MDH were scarce (Table 3), including cell wall-related processes “beta-glucan biosynthetic process,” “lignin metabolic process,” “phenylpropanoid metabolic process,” and “cell wall organization or biogenesis,” and others like “response to carbon dioxide” and “sucrose metabolic process.” Other shared DE tags included isoforms involved in tetrahydrofolate biosynthesis and a subtilase serine protease (Figure 6). Among DE TFs in YI, we found B3-ARF isoforms (Auxin response factors with B3 and PB1 domains) induced under both treatments (Table S6). However, reciprocal alignment revealed that they belong to different genes (blastn, alignment coverage 48% and percentage of identity 63%). The most similar protein of the isoform in the D treatment was ARF21, also known as OsARF7b, whereas the closest homolog of the isoform found under MDH was ARF11.

Besides B3-ARF TFs, only a few other isoforms were up-regulated in Scarlett YI under MDH, corresponding to an elongation factor EF-1, a DNA topoisomerase, a kinesin motor domain, CCA1/LHY, and a condensing complex subunit protein. All the others were down-regulated, whose enriched GO terms included “cellulose biosynthetic process,” “xylan biosynthetic process,” “flavonoid biosynthetic process,” “mitotic chromosome condensation,” “plasmodesmata-mediated intercellular transport” and “mucilage extrusion from seed coat” (Table 3). Other differences with respect to the D treatment were the involvement of enzymes from thiamine biosynthesis, triglyceride catabolism, epoxidation, berberine alkaloid biosynthesis or auxin biosynthesis (Figure 6). Among repressed isoforms related with transporters, we found sugar and lysine-histidine transporters, a PRA1 family protein B2 (a protein family related to regulation of vesicle trafficking, (Kamei et al., 2008), and several ABC transporters (Table S6; Data Sheet 5). Other proteins (and protein domains) which were found DE only under MDH included an expansin-B3, a putative cell wall protein, a PMR5/Cas1p, and several germin-like proteins.

Under D, Scarlett YI showed almost twice as many induced than repressed isoforms. The number of enriched GO terms was greater than for all the other contrasts (Data Sheet 5), including numerous enriched processes (Table 3) and metabolic pathways (Figure 6), related with responses to abiotic stress (cell wall thickening, biosynthesis of wax, triglyceride mobilization, expansin-A7), development (seed, embryo and root development), central metabolism (starch, glucose, pyruvate, many amino acids, fatty acids biosynthesis, activation and beta-oxidation), hormones (ethylene, jasmonate), energy (ATP and NADP metabolism related proteins, F and V-type H+-transporting ATPases), nucleic acids and proteins metabolism, antioxidation, proteolysis, protein folding, numerous proteins involved in transport and vesicle trafficking (Table S6), tRNA synthetases, an up-regulated MADS-MIKC TF whose best hit in SwissProt is O. sativa subsp. japonica MADS6, several PKs (like CIPK30) and phosphatases (like phosphoinositide phosphatase SAC7), proteins involved in interactions and signal transduction (SNF2, ASPR1 topless-related protein 1, 14-3-3 protein epsilon, CypP450), cytoskeleton proteins (tubulin, myosin, fimbrin and villin domains), and even processes related with photosynthetic tissues, like biosynthesis of chlorophyll a or tetrapyrrole, or induction of a Rubisco activase.

All these evidences indicate that responses to D and MDH of Scarlett YI were different, and that reproductive tissues were undergoing large gene expression changes, especially under D.

We surveyed the literature reporting genes and proteins expressed in barley in response to water deprivation. The goal was to compare those sequences to the DE transcripts identified in this work. The studies listed in Table 4 include 5 microarray experiments, 8 based on proteomics, 1 RNAseq study, 1 QTL work, 1 surveying expression QTL and 1 meta-analysis. Most of them focused on barley plants under drought, with a few exceptions. The work “matsumoto2014” surveyed responses to desiccation, salt stress and ABA. In addition, both “ashoub2015” and “rollins2013” combined drought and heat stress. The meta-analysis “shaar-moshe2015” compared drought related genes from different plant species. Although many of these works (10) sampled leaves, other tissues were also analyzed in some of them (mainly shoots, roots, spikes and grain).

Out of 4389 DE tags (proteins, genes and transcripts) reported overall in the studies above, more than half (2730) were barley genes included in the meta-analysis “shaar-moshe2015” and, indeed, that study matches the largest number of DE tags of the current work. However, in relative terms, the most similar were those of “ashoub2013,” “ashoub2015,” “vitamvas2015,” “wang2015,” “kausar2013,” and “rollins2013,” in decreasing order, whose DE tags were also found in the present study in proportions ranging from 52 to 32% (Figure 7). Interestingly, these are all proteomics studies. The next most similar study, “chmiewlewska2012,” and also “wendelboe2012,” were both proteomics studies. However, in contrast with the previous ones, focused on leaf tissue, these also included roots. Note that when considering only DE proteins from leaf samples, the agreement with “chmielewska2012” increased approximately from 22% to around 27%. On the other hand, DE transcripts from Scarlett YI under D matched the largest percentage of DE tags from the surveyed studies.

We also recorded the number of DE tags found in individual contrasts of our study, which had already been identified in previous studies. These figures for the four main contrasts of our study, Scarlett YI under D, Scarlett YI under MDH, SBCC073 LF under MDH and Scarlett LF under MDH, were 49, 30, 58, and 55%, respectively. The largest figures found for the leaf contrasts likely reflect the prevalence of studies which sampled LF tissues.

A total 466 DE isoforms were not found in previous studies, whereas 153 were in just one study and 54 in two. Only 23 DE isoforms were in common in three or more studies. These DE isoforms included several 70 and 90 kDa heat shock proteins, a S-methyltransferase and a S-adenosylmethionine synthase 2 from S-adenosyl-L-methionine cycle, an N-methyltransferase involved in choline biosynthesis, transcripts related with photosynthesis and carbon fixation, a sucrose synthase, a phosphoglycerate mutase and a triose-phosphate isomerase, a glutathione peroxidase, a ferredoxin-NADP+ reductase, a phenylalanine/tyrosine ammonia-lyase, a ClpC subunit from an ATP-dependent Clp protease, an ATP synthase and a V-type H+-transporting ATPase subunit, an aspartate kinase, a protein with Potato inhibitor I family domain and a spermidine synthase (Table 5).

DE isoforms were clustered based on their expression patterns across samples (Figure S11), with the aim of identifying shared regulatory motifs in their upstream genomic regions. We obtained 23 clusters, 14 of them with more than 10 isoforms (Table S7). Several clusters contained mostly isoforms from a given contrast while others had mixed DE tags from different treatments (Figures S12, S13).

In order to validate the expression-based gene clusters, they were tested for Gene Ontology (GO) enrichment. Moreover, to test the hypothesis that co-expressed genes might share cis-regulatory sequences, their upstream sequences were subjected to motif discovery algorithms and the DNA motifs found were annotated. Finally, the resulting regulatory motifs were compared to the binding predictions of DE expressed TFs identified in this work, trying to link these TFs to clusters of DE tags.

The results are summarized in Figure 8. Upstream sequences of genes from cluster 1, with functional annotations related to the metabolism of carbohydrates, contain a wtATAAAAGw site, which is similar to motifs of TATA-binding proteins and Dof TFs (Yanagisawa, 2002). We observed a C2C2-Dof TF up-regulated in SBCC073 LF under MDH (see previous sections), although we were not able to identify DNA-binding domains associated to it. Therefore, we cannot confirm whether or not C2C2-Dof protein binds to this motif to regulate genes in cluster 1, but the possibility deserves further investigation. Promoter sequences of genes in clusters 9 and 10, which group mostly transcripts down-regulated in LF under MDH, contain sites identical to the consensus of CCA1/LHY, which belongs to the MYB/SANT family (Green and Tobin, 1999). These sites were independently predicted by oligo-analysis (AAAATATCTy) and dyad-analysis (aAAAkaTCTw), indicating that they are high-confidence predictions. Genes of this cluster are annotated as components of thiamine biosynthesis in the chloroplast. Accordingly, CCA1/LHY, which was up-regulated in SBCC073 and Scarlett samples under MDH, binds to the same motif (aAAATATCTkY). Cluster 12 had predicted yaCGTACGtr cis-elements. Genes in this cluster were induced in LF under MDH, and are annotated as heat shock proteins. Finally, genes in cluster 14 are annotated as components of salinity response, and share cis-elements of consensus smACACTbm.

Out of 11 DE TFs, 7 were associated with DNA-binding domains (Table 6), including CCA1/LHY (see above), the MYB-related TF of unknown function DE in SBCC073 LF under MDH, the MADS-MIKC up-regulated in Scarlett YI under D (AwRGaAAaww), the B3-ARF TFs induced in Scarlett YI either under D or MDH (yTTGTCtC), the bZIP up-regulated in Scarlett LF under MDH (cayrACACGTgkt) and the AP2/ERF-AP2 down-regulated in SBCC073 LF under MDH (CACrrwTCCCrAkG). It is possible that these genes were in part regulating the changes in gene expression in response to the treatments. However, these could not be linked to the motifs identified in promoters.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.