Based on a survey of 156 elite commercial wheat varieties in the North China Plain, we selected Luyuan 502 for its drought tolerance (Li, 2014). Seeds of Luyuan 502 were field sown in early October 2014 at the Experimental Station of the Hebei Agricultural University in Gaocheng County, China (114.84° E 38.03° N, 56 m above sea level). Plants were arranged in three blocks (A, C, and D), each being 12 m long and 8.5 m wide, with a rain shelter to keep the plants from rain and snow (Figure 1). To simulate local farming practices, a two-time irrigation strategy was adopted. Block D (control block) was watered twice, whereas block A and block C was watered once. On 1 April 2015, when the first internodes of most plants (>80%) were visible (Feekes' growth stage 6), the plants in blocks C and D were watered. On 1 May 2015, when most of the plants (>80%) had begun to flower (Feekes' growth stage 10.5.1) (Large, 1954), plants in blocks A and D were watered. About 574 kg of water (equal to 60 mm) was used for each irrigation. During wheat growth, the relative soil water content (RSWC) was recorded at nine time points (T1 to T9) to monitor the level of drought stress in the three blocks (Figure 2). The RSWC was calculated as RSWC = [(current pot weight − soil dry weight)/weight of soil watered to field capacity] × 100.

Plant samples were taken in blocks A and D at T4 (7 April; 7 days after first irrigation), T5 (14 April; 14 days after first irrigation) and T6 (21 April; 21 days after first irrigation), and in blocks C and D at T8 (7 May; 7 days after second irrigation), and T9 (14 May; 14 days after second irrigation) (Figure 2). At T4, T5, and T6, the upper parts of main tillers were collected, while at T8 and T9, the developing spikes were collected after removing awns and rachis. The developmental stages of plants at these five time points were checked by microscope or the naked eye. Plant materials from five individual plants in each block were pooled as a biological replicate for RNA extraction. Three biological replicates were collected at each time point, resulting in 30 samples for cDNA library construction and sequencing (three biological replicates × two blocks × five time points = 30 samples). The collected samples were frozen immediately in liquid nitrogen and stored at −80°C until required. To examine the effects of drought on plant growth and development in each block, total dry biomass before flowering, total dry biomass at maturity, spike and spikelet numbers, number of fertile and infertile spikelets per spike, thousand grain weight and grain yield were recorded. The differences between grain setting rate and grain yield difference to the control were also calculated.

Total RNA was extracted using TRIzol reagent (Invitrogen) according to the manufacturer's instructions. DNA was removed by digestion with RNase-free DNase (Qiagen), and RNA was purified and concentrated using an RNeasy column (Qiagen). RNA quality was evaluated by gel electrophoresis, spectrophotometer analysis and an Agilent 2100 bioanalyzer. RNA samples were sent to Novogene Bioinformatics Technology Co. Ltd., Beijing, China (http://www.novogene.cn) for further processing and sequencing. The sequencing libraries were sequenced on Illumina HiSeq 2500 (Illumina, San Diego, USA) and 125/150 bp paired-end reads were generated. The RNA sequences have been deposited at the National Center for Biotechnology Information (NCBI) with the accession number of SRP102636. The raw reads were filtered by removing adapter sequences, ambiguous nucleotides (if the proportion of “Ns” exceeded 10%), and low-quality sequences (when the percentage of bases with low Phred quality score ≤20 is greater than 50% in a read). At the same time, Q20, Q30, and GC contents of the clean data were calculated. All downstream analyses were based on clean data.

The clean reads were aligned to wheat reference genome sequences released by the International Wheat Genome Sequencing Consortium (IWGSC) (ftp://ftp.ensemblgenomes.org/pub/release-23/plants/fasta/triticum_aestivum/dna/) using TopHat 2.0.12 (Trapnell et al., 2010). Mismatches of no more than two bases were allowed in the alignment. After alignment, HTSeq v0.6.1 was used to count the read numbers that were mapped to each gene. The expected number of Fragments Per Kilobase of transcript sequence per Million base pairs sequence (FPKM) of each gene was calculated to determine the expression values of this gene. Pearson's correlation coefficients between the three biological replicates were calculated in the R software package to quantify the correlation between biological replicates.

A consensus FPKM across the three biological replicates was calculated and used in the DEG analysis (Trapnell et al., 2010). For each gene, the relative homologous rice gene, homologous rice gene annotation, and NCBI non-redundant protein database annotation were investigated. DEGs were identified using DESeq 1.10.1 (Anders and Huber, 2010). To assess the variability among samples, we performed principal component analysis (PCA) for the wheat genes identified at T4, T5, and T6 using the prcomp command with default parameters in the R software package (Robinson et al., 2010).

Five pair-wise comparisons between the three blocks were made to identify the drought-responsive genes. At T4, T5, and T6, the comparison was made between blocks A and D while at T8 and T9, the comparison was made between blocks C and D. A corrected p-value < 0.05 was set as the criteria for determining DEGs. We created heatmaps to demonstrate the gene expression data of T4 using Java TreeView (Robinson et al., 2010). The subcellular locations of proteins for the DEGs at T4 were analyzed using PSORT Prediction software (http://psort.hgc.jp/form.html).

To study the biological significance of the DEGs, an enrichment analysis of gene ontology (GO) terms was conducted with GO seq R package. Statistically significant over-representation of GO categories (corrected p-value < 0.05) in response to the drought treatment was determined separately for each time point. Similarly, KOBAS software was used to test the statistically enriched pathways associated with the DEGs at each time point in the KEGG (Kyoto Encyclopedia of Genes and Genomes) database.

To confirm the RNA-seq results, among the identified DEGs in Luyuan 502, 21 genes were randomly selected and assessed using RT-qPCR. RT-qPCRs were conducted on a StepOne Plus Real-Time PCR system (Applied Biosystems) using SYBR Green I for the detection of PCR products. Each reaction was performed in a final volume of 16 μL, containing 8 μL SYBR Green PCR Master Mix (Applied Biosystems), 250 nM of each primer and 50 ng cDNA template. The thermal cycling conditions were 94°C for 10 min, followed by 40 cycles of 94°C for 15 s, 55°C for 30 s, and 60°C for 1 min, with fluorescence detection at the end of each cycle. The amplification of a single product per reaction was confirmed by melting curve analysis. All reactions were performed in three technical triplicates. Wheat α-tubulin (forward: ATCTGTGCCTTGACCGTATCAGG; reverse: GACATCAACATTCAGGACACCATC) was used as an internal reference gene to normalize Ct values of each reaction (Chen et al., 2016), which were determined using the CFX96 software with default settings. The sequences of the 21 genes and primers used in the RT-qPCR analysis are listed in Table S5.

A drought-tolerant variety Cangmai 6001 and a drought-sensitive variety Hanmai 9 were subjected to normal watering as a control, and 5 days of non-watering as a drought treatment until the RSWC reached the threshold (35%) at the tilling stage in the glasshouse. To confirm their different sensitivities to drought, four drought-related physiological indexes—ascorbate peroxidases (APX) activity, catalase (CAT) activity, H₂O₂ content and MDA content—were measured in both varieties before sampling.

Wheat leaves subjected to drought stress or mock treatments were sampled for physiological index determination. The MDA (malondialdehyde) level was estimated according to Li et al. (2014). H₂O₂ (hydrogen peroxide) accumulation was assessed using commercial kits (Jiancheng Biotech Inc., Nanjing, China) according to Yang et al. (2016). Each sample was homogenized in pre-cooled phosphate-buffered saline (PBS) using 1 mL of buffer per 0.1 g of fresh tissue. The homogenate was centrifuged at 10,000 g for 10 min at 4°C. Freshly isolated supernatant fractions were used immediately for measuring H₂O₂ content. Adduct formation was measured spectrophotometrically at 405 nm using Thermo Scientific Multiskan FC (Shanghai, China) in strict accordance with the manufacturer's instructions. Protein contents were determined using an Enhanced BCA Protein Assay Kit (Beyotime, Shanghai, China). The activities of antioxidant enzymes, including catalase (CAT) and ascorbate peroxidase (APX), were measured as described previously (Rao et al., 1996; Li et al., 2011; Tian et al., 2014). The assay of enzyme activity was carried out using a spectrophotometer at 25°C.

As we mainly focused on genes differentially expressed during early reproductive stages, we selected five DEGs during early reproductive stages for RT-qPCR validation: Traes_5DS_CCCDA48421 (T4), Traes_5BS_9584239E51 (T4), Traes_2DL_77F25CE27 (T4 & T5), Traes_3DL_304C8DD67 (T4 & T5) and Traes_7DS_1D74598FD (T4 & T6). These five DEGs were associated with floral organ development, stomatal movement or photosynthesis activity etc. RNA of Hanmai 9 and Cangmai 6001 was extracted using TRIzol reagent (Invitrogen). Three trials were conducted for the measurement of four drought-related physiological indexes and RT-qPCR analysis. A general mean across each trial was calculated and used. Two-tailed unpaired Student's t-tests was used to determine if the differences between the two varieties were significant or not. Wheat α-tubulin was used as the reference gene (Chen et al., 2016). RT-qPCR primer sequences are listed in Table S5.

Initially, RSWC in blocks A, C, and D was similar at T1 and T2. The first irrigation on 1 April in blocks C and D, increased RSWC in these two blocks to 62.2 and 63.0% at T3, then decreased to 52.8 and 53.5% at T4, 46.5 and 46.9% at T5, and 35.5% at T6. In contrast, RSWC in the unirrigated block A changed little (36.7–30.8%) during the same period. The second irrigation on 1 May in blocks A and D increased RSWC in these two blocks to a much higher level than in block C at T7, T8 and T9 (Figure 2).

Total plant biomass was a direct growth parameter to measure the effects of drought on wheat growth and development. As expected, the plants in block A which suffered drought stress during early reproductive stages had the lowest total dry biomass before flowering. At maturity, due to second irrigation, the biomass of plants in block A had increased to a similar level as that in block C, but was still much lower than block D (control). A similar trend was observed for grain setting rate and yield. Blocks A and C had lower values (−6.6 and −6.1%; −4.9 and −2.4%) on these two important traits compared with block D. Detailed information on the traits measured in the three blocks is shown in Table 1.

Microscopic inspection determined that plants at T4, T5, and T6 were at the pistil and stamen differentiation stage, anther differentiation stage and tetrad stage, respectively. Plants at T8 and T9 were at the early flowering and grain formation stages, respectively (Figure 1).

The quality of total RNA was good (Figure S1). Of the 1.7 billion raw reads generated from the 30 cDNA libraries, approximately 1.6 billion clean reads were identified, ranging from 47 million to 69 million reads per library. Clean reads were mapped to the wheat genome released by IWGSC using TopHat 2.0.11. The mapping rates of each library ranged from 98.74 to 99.14%. A total of 99,411 genes from IWGSC and 16,244 novel isoforms of known genes were identified (Table 2). The genes with FPKM (transcript abundance of the gene) above 1.0 accounted for around 40%, and those above 3.0 accounted for around 25%. The biological replicates were highly correlated. The Pearson's correlation coefficients among the three replicates for each time point ranged from 0.82 to 0.98 (Figure S2), indicating that the three replicates were consistant. The PCA analysis indicated that most of the variation in gene expression among different ellipses was a consequence of the developmental stage. Furthermore, two distinct groups formed within each ellipse; one comprised plants in block A and the other comprised plants in block D, indicating that the transcriptomes of drought and irrigated plants clearly differed from each other (Figure 3).

Drought induced 309 DEGs across different time points: 226, 44, 8, 6, and 25 DEGs at T4, T5, T6, T8, and T9, respectively (Figure 5, Table S1). The GO terms and KEGG pathways associated with those DEGs are listed in Tables S2, S3, respectively.

At the pistil and stamen differentiation stage (T4), 193 up-regulated DEGs and 33 down-regulated DEGs were identified. Their expression patterns are shown in Figure 4. The genes in clusters 1–4 were up-regulated genes by drought, whereas those in cluster 5 were down-regulated genes. Furthermore, subcellular location analysis of 100 proteins of the DEGs revealed that they were mainly located in the nucleus (18), peroxisomes (18), mitochondria (13), plasma membrane (13), and chloroplasts (9) (Figure 4 and Table S4). These results suggest that these organelles play critical roles in drought stress, and that the peroxisome, nucleus and mitochondria may be more sensitive to drought stress damage in wheat at early developmental stages. GO analysis showed that the “protein disulfide oxidoreductase activity” and “cell redox homeostasis” were enriched (Table 3). Four KEGG pathways-“galactose metabolism,” “circadian rhythm–plant,” “starch and sucrose metabolism” and “flavonoid biosynthesis”-were significantly impacted by drought stress (Table 4).

At the anther differentiation stage (T5), 37 up-regulated DEGs and seven down-regulated DEGs were identified. Three GO terms “magnesium ion binding,” “carbon fixation,” and “ribulose-bisphosphate carboxylase activity” were over-represented (Table 3). The “carbon fixation in photosynthetic organisms,” “glyoxylate and dicarboxylate metabolism” and “carbon metabolism” pathways were significantly enriched (Table 4). Interestingly, five DEGs were shared between the pistil and stamen differentiation stage (T4) and the anther differentiation stage (T5) (Figure 5). Of these, the drought up-regulated DEG Traes_3DL_304C8DD67 was annotated to the vacuolar-processing enzyme (VPE), which has been associated with growth inhibition, cell death and stomatal movement in Arabidopsis (Hara-Nishimura et al., 2005; Yamada et al., 2005; Albertini et al., 2014).

At the tetrad stage (T6), two up-regulated DEGs and six down-regulated DEGs were identified. The pathways of metabolism of porphyrin and chlorophyll were significantly affected (Table 4). This enriched pathway only contains one DEG, Traes_2DS_BB50DEEF8, which was homologous to the rice gene LOC_Os03g20700 encoding magnesium chelatase H subunit (CHLH). In our study, the expression of Traes_2DS_BB50DEEF8 in the water-stressed block A was 4.53-fold higher than that in the control block D, suggesting a positive response of this gene to drought stress. Traes_2DS_BB50DEEF8 was one of the five genes shared between T4 and T6 (Figure 5). The CHLH subunit affected plant hormone abscisic acid (ABA) signaling in stomatal guard cells, and overexpression of the CHLH gene in guard cells improved plant drought tolerance by promoting stomatal closure (Tsuzuki et al., 2013).

Compared with earlier reproductive periods, the number of DEGs detected during flowering declined. During early flowering (T8), six DEGs were identified, including one up-regulated gene and five down-regulated genes. The down-regulated gene TRAES3BF077200010CFD_g encoded LTPL38 which may contribute to mineral accumulation in wheat grain (Singh et al., 2014). At the grain formation stage (T9), 25 DEGs were identified (two up-regulated and 23 down-regulated). The GO terms “alpha-amylase activity” and “calcium ion binding” were enriched (Table 3). The “starch and sucrose metabolism” pathway was enriched (Table 4).

To verify the RNA-seq results, RT-qPCR analysis was conducted against 21 randomly selected genes. The expression patterns of these 21 genes assessed by RT-qPCR correlated well (R² = 0.9345) with those obtained from the RNA-seq analysis (Figure 6). These results confirm the accuracy of the transcriptome changes obtained by RNA-seq in this study.

The two varieties used for RT-qPCR validation responded differently to drought in terms of drought-related physiological indexes. Cangmai 6001 had higher APX and CAT activities and lower MDA and H₂O₂ contents than Hanmai 9 in response to drought stress, thus confirming the better drought tolerance of Cangmai 6001 (Figure 7).

The RT-qPCR confirmed the involvement of five selected DEGs in drought tolerance under different genetic backgrounds. For example, the up-regulated DEG Traes_3DL_304C8DD67 (encoding VPE) shared by T4 and T5 increased in Cangmai 6001 and Hanmai 9, but more so in the drought-sensitive Hanmai 9. The two down-regulated CC-type GRXs genes (Traes_5DS_CCCDA48421 and Traes_5BS_9584239E51) declined less in Cangmai 6001 than Hanmai 9, while Traes_2DL_77F25CE27 and Traes_7DS_1D74598FD increased more in Cangmai 6001 than Hanmai 9 under drought, which were in line with RNA-seq results (Figure 8).

The total dry biomass of plants in blocks A and C declined, indicating a negative impact of drought on wheat growth and development, compared with the control block D. Drought had the greatest effect on grain number related traits of plants in block A, with reductions of 6.45% in spike number, 11.02% in spikelet number per spike and 6.64% in grain setting rate, compared with the control block D. Plants in block C had the lowest thousand grain weight, declining by 6.24% compared with control block D (Table 1). It can be concluded from these observations that drought during early reproductive periods (before flowering) mainly impacted grain number related traits whereas drought during flowering affected grain weight. This is consistent with other studies on wheat, where drought during stem extension caused floret and whole spikelet death and drought during grain filling reduced grain size and weight (Oosterhuis and Cartwright, 1983; Dorion et al., 1996; Ji et al., 2010; Pradhan et al., 2012).

A critical and direct parameter for measuring drought damage is grain yield. It is well established that drought during flowering and grain filling can severely reduce wheat yield (Farooq et al., 2014). Some studies have focused on the drought-responsive genes during these periods (Aprile et al., 2009, 2013). In contrast, genes responsive to drought at early reproductive periods were mostly unknown. A study conducted by Ugarte et al. (2007) indicated that the high temperature during stem elongation reduced wheat yield by 46%, compared with that from booting to anthesis (27%) and heading to anthesis (15%) (Ugarte et al., 2007) Our results with drought stress followed a similar trend, where grain yield in block A (drought stress applied during early reproductive stages before flowering) declined by 6.1% compared with 2.4% in block C (drought stress applied after flowering periods) (Table 1). Therefore, like heat, drought stress during early reproductive periods may affect wheat yields more than drought during flowering periods. Investigation of genes responsive to drought during the early reproductive stage is needed.

In contrast to previous wheat transcriptome studies where hundreds or even thousands of DEGs have been identified (Aprile et al., 2009, 2013; Li et al., 2012; Liu et al., 2015), only 309 DEGs were identified in our study. The likely reason for this is our experimental design (irrigation strategy). The small difference in the level of drought stress between control and test samples can reduce the number of DEGs identified in transcriptome analyses (Aprile et al., 2009). We adopted a two-times irrigation schedule in the control block as this schedule is widely adopted in the North China Plain, where about two-thirds of China's total wheat output is produced (Lu and Fan, 2013). This two-times irrigation strategy is sufficient to maintain wheat production under drought stress conditions (Li et al., 2005; Yang et al., 2006). By simulating this strategy, we aimed to study the molecular mechanisms to maintain wheat production in local areas and to provide useful genetic resources for drought tolerance improvement.

In our study, the most DEGs were identified at T4 (pistil and stamen differentiation stage). The chalcone synthase gene (CHS) and flavonoid 3′-hydroxylase gene (F3H) in the enriched flavonoid biosynthesis pathway, and fructan biosynthesis genes (1-FFT, 6-SFT, and 1-SST) in the enriched starch and sucrose metabolism pathways were known to be associated with dorught tolerance in wheat (Kawakami and Yoshida, 2002; Ithal and Reddy, 2004; Ji et al., 2010; Ma et al., 2014; Wang et al., 2016). The two enriched GO terms at this stage (protein disulfide oxidoreductase activity and cell redox homeostasis) were associated with the same group of DEGs, including one up-regulated and seven down-regulated DEGs. These DEGs were annotated to CC-type glutaredoxins (GRXs), which are small (10–15 kDa) oxidoreductases that catalyze the reduction of disulfide bonds of their substrate proteins in the presence of glutathione (GSH). The loss-of-function mutants of CC-type GRXs genes in Arabidopsis, maize and rice had fewer petal primordia, disrupted anther lobe differentiation, and disabled meiotic entry of sporogenous cell progenies (leading to male sterility; Xing and Zachgo, 2008; Wang et al., 2009; Hong et al., 2012; Kelliher and Walbot, 2012). Therefore, down-regulation of CC-type GRXs genes in wheat may explain the spikelet death caused by drought during the early reproductive stage (Oosterhuis and Cartwright, 1983), but this needs further investigation. Meanwhile, a previous study reported increased sensitivity to oxidative stress in Arabidopsis mutants with knock down GRX genes (Li et al., 2014). Under drought stress, the two CC-type GRXs genes (Traes_5DS_CCCDA48421 and Traes_5BS_9584239E51) declined less in Cangmai 6001 than in Hanmai 9, which was in line with its drought tolerance (Figures 7, 8).

The subcellular location of the proteins of the DEGs at T4 indicated that the peroxisome, nucleus and mitochondria might be sensitive to drought stress damage. This finding differed from those of Hippophae rhamnoides seedlings, where the proteins responsive to drought stress were mainly located in the chloroplast, mitochondria and secretory pathway (He et al., 2016). Plant oxidative stress can be caused by excess H₂O₂ accumulation under drought stress, which severely damages biomolecules due to the elevated and non-metabolized cellular H₂O₂ (Sofo et al., 2015). In this study, Hanmai 9 accumulated significantly more H₂O₂ under drought stress than Cangmai 6001, and significantly less H₂O₂-metabolizing enzymes such as CAT and APX. Peroxisomes are involved in the response of plants to biotic and abiotic stresses (Kaur et al., 2009; Hu et al., 2012; Smith and Aitchison, 2013). Importantly, plant peroxisomes are involved in conserved functions (e.g., detoxification of reactive oxygen species, or ROS) and plant-specific functions (e.g., photorespiration and metabolism of hormones) (Kaur et al., 2009; Hu et al., 2012; Smith and Aitchison, 2013). The DEGs encoding proteins in the peroxisome may be involved directly or indirectly in ROS scavenging in wheat at early developmental stages under drought.

At the anther differentiation stage (T5), the genes involved in photosynthetic activities were significantly affected by drought. The enriched GO terms (ribulose-bisphosphate carboxylase activity and carbon fixation) contained three up-regulated DEGs encoding ribulose-1,5-bisphosphate carboxylase/oxygenase (Rubisco). Rubisco catalyzes the first step in net photosynthetic CO₂ assimilation and is a central component of photosynthesis—a new “green revolution” (Whitney et al., 2011; Long et al., 2015). Studies have shown that Rubisco activase acclimates in situ to high temperature when the stress is imposed at a slow rate in the field. Altered expression of Rubisco activase might be crucial for continued CO₂ fixation under drought stress, protecting plant photosynthetic capacity (Crafts et al., 1998; Law and Crafts-Brandner, 2001). In addition, Pelloux et al. (2001) and Fu et al. (2007) reported that the expression pattern of some Rubisco genes under drought stress had a temporal manner in Aleppo pine and rice, which increased at an early stage of growth but then decreased at a later stage. Similarly, of the three wheat Rubisco genes in our study, one (Traes_7DS_1D74598FD) was up-regulated by drought at T5 and later down-regulated at T8. Further experiments are required to investigate the mechanism for different expression of Rubisco genes during wheat growth.

At tetrad stage (T6), eight DEGs were identified. Most of these DEGs had unknown functions in wheat. However, the stomatal movement related gene Traes_2DS_BB50DEEF8 (encoding CHLH) which was shared between T6 and T4 arose our interests. As plants lose over 95% of their water via transpiration through stomata, the engineering of stomatal activity is believed to be a promising approach to reduce the water requirement of crops and to enhance productivity under stress conditions (Schroeder et al., 2001).

Compared with early reproductive stages, drought stress during flowering had less effect on wheat gene expression in the developing grains because fewer DEGs were detected during early flowering (T8; 6 DEGs) and grain filling (T9; 25 DEGs). We found that drought stress induced several grain development-related genes. For example, during early flowering, the DEG TRAES3BF077200010CFD_g was annotated to the lipid-transfer protein LTPL38, which may contribute to mineral accumulation in wheat grains (Singh et al., 2014). At grain filling (T9), four DEGs in the enriched KEGG pathway “starch and sucrose metabolism” stood out: three were annotated to the wheat a-Amy3 gene that encodes a-amylase, and one was annotated to a rice gene encoding glycoside hydrolase family 31 proteins, a-glucosidase (ONG2) (Nakai et al., 2007). Both a-amylase and ONG2 can hydrolyze raw starch granules (Sissons and Macgregor, 1994). a-Amy3 mainly provides energy for early developing wheat grains in the spikes (Zanetti et al., 2000). Rice ONG2 and its mRNA is only produced during ripening and remains an active and key enzyme in starch degradation metabolism during the initial stage of germination (Sissons and Macgregor, 1994; Nakai et al., 2007). The change in these genes under drought may affect wheat grain development; we will investigate the functions of these genes in future experiments.

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.