Comparative Analysis of the Transcriptional Response of Tolerant and Sensitive Wheat Genotypes to Drought Stress in Field Conditions

Drought stress is one of the most adverse environmental limiting factors for wheat (Triticum aestivum L.) productivity worldwide. For better understanding of the molecular mechanism of wheat in response to drought, a comparative transcriptome approach was applied to investigate the gene expression change of two wheat cultivars, Jimai No. 47 (drought-tolerant) and Yanzhan No. 4110 (drought-sensitive) in the field under irrigated and drought-stressed conditions. A total of 3754 and 2325 differential expressed genes (DEGs) were found in Jimai No. 47 and Yanzhan No. 4110, respectively, of which 377 genes were overlapped, which could be considered to be the potential drought-responsive genes. GO (Gene Ontology) analysis showed that these DEGs of tolerant genotype were significantly enriched in signaling transduction and MAP (mitogen-activated protein) kinase activity, while that of sensitive genotype was involved in photosynthesis, membrane protein complex, and guard cell differentiation. Furthermore, 32 and 2 RNA editing sites were identified in drought-tolerant and sensitive genotypes under drought compared to irrigation, demonstrating that RNA editing also plays an important role in response to drought in wheat. This study investigated the gene expression pattern and RNA editing sites of two wheat cultivars with contrasting tolerance in field condition, which will contribute to a better understanding of the molecular mechanism of drought tolerance in wheat and beyond.


Introduction
Drought is one of the most hazardous environmental stresses limiting plant growth and development, which gradually becomes a major threat to the world's agricultural production nowadays and leads to huge yield losses of major crops annually [1][2][3].According to statistics, only 20% of cropland worldwide is available for irrigation and it provides approximately 40% of global food production, whereas rain-fed agriculture provides the remaining 60% [4].Understanding of the molecular mechanism of drought response in crops is crucial for genetic improvement and breeding for drought tolerance, which could meet the challenge of population boom and food security in the 21st century [5].To date, extensive studies have been conducted to uncover the complexity of mechanisms of plants in response to drought at the morphological, physiological, and molecular level [6][7][8][9].Generally, plants rapidly close their stomata when subjected to drought stress to decrease water losses from leaves, and then a series of downstream response processes are triggered [10][11][12].On a cellular level, an osmotic adjustment is first activated and the osmolytes such as proline, glutamate, and mannitol as well as sorbitol and trehalose are accumulated, which could prevent the plant cell from dehydration by increasing the osmotic stress to keep cell membrane integration and enzyme function under drought stress [13,14].In addition, these substances have been used as the physiological indictors to assess the drought tolerance in plants.On molecular level, several genes and proteins have been reported to be induced in response to drought tolerance, such as dehydration-responsive element binding protein (DREB), C-repeat-binding factor (CBF) and myeloblastosis oncogene (MYB) [15][16][17].These drought-responsive genes are mostly transcription factors, which play the hub role in drought-signaling transduction and regulation pathway, such as ROS (reactive oxygen species) [18,19].It is well known that drought is a complex quantitative trait, which is affected by various factors including environmental condition, genotype, developmental stage, drought severity, and duration [16,[20][21][22].Thus, an increasing number of studies are further needed to uncover the complex regulatory mechanism of drought tolerance.
Wheat is one of the most important cereal crops all over the world, occupying 17% of cultivated lands and providing the main food source for 30% of the global population [23,24].Furthermore, wheat is widely grown in a large range of lands under both irrigated and rain-fed conditions.Due to global warming, drought has become the most serious environmental constraint to wheat production and has caused about 5.5% average loss annually [25,26].Therefore, mining and using drought-tolerant genes to improve wheat varieties with enhanced drought tolerance is urgently needed to meet the challenge of global climate change and food security [27].Recently, great progress has been made in revealing the molecular mechanism of drought response and many drought-responsive genes have been identified in wheat [28][29][30][31].It demonstrated that over-expression of TaNAC69 could enhance the drought tolerance in bread wheat [32], and TaSAP5 could alter the drought stress responses by promoting the degradation of DRIP (DREB interacting protein) proteins [33].In addition, the wheat MYB gene TaMYBsdu1 is found to be up-regulated under drought stress but showed differential expression between tolerant and sensitive genotypes, suggesting it plays a crucial role in regulating drought tolerance [34].With the advent of high-throughput sequencing technology, RNA deep sequencing has been widely used to investigate the gene differential expression profiles involved in drought response in wheat at the transcriptome level, which provided a direct and effective method to identify drought-inducible genes and also contribute to better understanding of drought-signaling pathways [35,36].It has been demonstrated that different wheat cultivars had diverse molecular basis for drought response and adaptation and the performance of wheat under controlled conditions showed less correlation with field performance [20,37,38].However, most current studies were conducted using the limited genotypes (drought-tolerant or sensitive) under the controlled conditions.The transcriptional difference of drought-tolerant and drought-sensitive wheat varieties under irrigated and drought-stressed field conditions are not well understood up to now.Here, we investigated the gene expression profiles of two elite wheat cultivars with contrasting drought tolerance, namely Jimai No. 47 (drought-tolerant) and Yanzhan No. 4110 (drought-sensitive) in the field under irrigated and drought-stressed conditions to provide more information for better understanding of the molecular mechanism of drought tolerance in wheat and beyond.

Plant Samples and RNA Isolation
Two wheat elite varieties, Jimai No. 47 and Yanzhan No. 4110, were used in this study, of which Jimai No. 47 is a widely grown in the arid area of northern China with excellent drought tolerance and Yanzhan No. 4110 is a high-yield but drought-sensitive variety.These two varieties were grown in field plots of Luoyang A&F institute (Luoyang, Henan, China) with the same plant density in the 2017-2018 crop season.Each plot was 6 m 2 with 2 m in width and 3 m in length.The normal field management was used but differed in water condition.The irrigated treatment was applied three times in November 2017, January 2018, and March 2018 while the drought-stressed treatment was rain-fed with 280 mm of rain in overall crop season.Six replication and randomized block design were used.At the grain-filling stage, the flag leave of 10 wheat plants in each plot were randomly collected for RNA isolation.
Total RNA of the above prepared samples was firstly isolated using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) and treated with RNase-free DNase I to remove any contaminating genomic DNA according to the manufacturer's instructions.Then, the quality of RNA was checked by agarose gel electrophoresis and the quantity was measured by NanoDropND-1000 Spectrophotometer (NanoDrop Technologies, Waltham, MA, USA).Finally, the equal quantity RNA of three replications for the same treatment was pooled and used for RNA sequencing.

RNA Deep Sequencing and Data Analysis
The four pooled RNA samples were used to construct the RNA sequencing libraries following Illumina's standard pipeline (Illumina, San Diego, CA, USA).High-throughput sequencing was performed on an Illumina HiSeq3000 platform following the standard protocol at Sangong Bio-Technology Corporation (Shanghai, China).The RNA-seq data have been deposited into the genome sequence archive (GSA) database in BIG Data Center, Beijing Institute of Genomics (BIG), Chinese Academy of Sciences, with the accession number CRA001148, and are publicly accessible at http://bigd.big.ac.cn/gsa [39].
The quality of raw data obtained for sequencing was tested by FastQC and then quality filtered by FASTX-toolkit.The adaptor contamination, low-quality reads (quality scores < 20), reads with ambiguous "N" bases more than 10% bases, as well as reads less than 20 bases were removed to obtain the clean data.Then, all the clean data was aligned to wheat reference genome The International Wheat Genome Sequencing Consortium (IWGSC) RefSeq v1.0 [40] by HISAT2 with the default parameters (version 2.0.5)[41].Then, gene expression levels were calculated based on the FPKM (fragments per kilobaseof exon per million fragments mapped) method and only the high-confidence gene models annotated in wheat RefSeq v1.0 were used.Pearson correlations between biological replicates were also calculated based on the FPKM values of all expressed genes, to assess the reliability.

Identification of Differentially Expressed Genes (DEGs)
Differential expression genes were identified using DESeq2.0[42].The adjusted p-value < 0.05 and fold-changes > 2 (Log 2 (treatment/control) | ≥ 1) were used as thresholds for differentially expressed analysis.Then, a K-means clustering was used to extract the fundamental patterns of gene expression [43].GO terms that are significantly over-represented in each cluster were determined by the AgriGO [44].Singular Enrichment Analysis (SEA) in AgriGO was used to detect over-represented GO categories in each cluster compared to the whole genes.GO terms with corrected FDR of less than 0.05 were taken as significant ones.KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways were obtained by searching against KEGG database [45].

Analysis of RNA Editing Sites
To identify the RNA editing sites of these two wheat varieties between irrigated and drought-stressed, all the clean reads of RNA-seq were mapped to the wheat genome IWGSC RefSeq v1.0 using SPRINT (SnP-free RNA editing IdeNtification Toolkit) software with default parameters [46].To avoid any errors, the independencies mapping software BWA (http://biobwa.sourceforge.net)[47] and SNP nomenclature tools Samtools (http://www.htslib.org)[48] were integrated to identify the SNP between RNA reads and reference genomic DNA.The overlapping results obtained by both methods were retained for further analysis.Then, using the irrigation samples as background, the SNPs at the same position in drought-stressed sample which were different from the background were considered as the potential RNA editing sites.Finally, the editing sites were further filtered using the following parameters: (1) the edited sites having more than 5 mapped reads; (2) the ratio of editing reads and total mapped reads more than 50%; (3) the editing site identified in both biological replications.

RNA-seq Analysis
To enrich the knowledge of drought tolerance mechanism in wheat, we investigated the gene expression profiles of wheat cultivar Jimai No. 47 (drought-tolerant, JM) and Yanzhan No. 4110 (drought-sensitive, YZ) in the field under irrigated (I) and rain-fed (R) conditions by RNA-seq technology.A total of 8 pooled RNA samples were sequenced, and about 61.43 Gb of raw data was obtained, with an average of about 53 million pair-end reads with 150 bp in size for each sample (Table 1).After a quality filter, 55.81 Gb of clean data remained, representing 90.8% of the raw data.In addition, the clean reads of these samples ranged from 41,653,874 to 51,913,178, with an average of 48,402,376 reads representing 6.97 Gb.Then, all the clean reads were mapped to the wheat genotype Chinese Spring reference genome IWGSC_V1 (accessed from URGI (Unité de Recherche Génomique Info) database on 6 June 2018).Results showed that on average, about 80% of clean reads could be mapped to the wheat reference genome, and approximately 77% were uniquely mapped, which was similar to previous studies [9,35].The unmapped 20% reads might be due to genotype specifics or incompletion of the reference genome.Furthermore, most of RNA-seq reads were mapped onto the exon region of the reference genome, although they were also mapped to intergenic and intron regions with very low frequency.Based on the mapping result, the expression level of the annotated genes in these mapped regions were obtained.Then, we used the gene expression level to calculate the correlating coefficient of two biological replicates of all samples to examine the repeatability and reliability.Results showed the biological replicates of all samples showed strong correlation relationships with the coefficient (R 2 ) of more than 0.90.Furthermore, we compared the abundance of expressed gene in these samples.Out of 110,790 high-confidence genes annotated in wheat genome, 58,498 (52.8%) genes were detected to be expressed in these samples, of which 45,258 genes were found to be expressed in all tested samples, representing the core gene set of wheat (Figure 1).A total of 50,981, 51,140, 52,622 and 52,961 genes were expressed in JM_I, JM_R and YZ_I and YZ_R, respectively.This means that compared to irrigation condition, more genes were induced to express by drought stress in both drought-tolerant and sensitive genotypes, which is consistent with previous studies [20,36].From the point of view of the genotype, the variety YZ has more abundant expressed genes than the JM variety.Additionally, 1167, 903, 1189 and 1223 genes were found to be specifically expressed in JM_I, JM_R and YZ_I and YZ_R respectively (Figure 1).

Identification of the Differentially Expressed Genes
To identify the drought-responsive genes, the differentially expressed gene (DEGs) were analyzed between different treatments and different genotypes using padj < 0.05 and |log2Ratio| ≥1 as thresholds.In drought-tolerant genotypes (JM), 2754 DEGs were identified between irrigation and drought condition, of which 1152 gene were up-regulated, and 1612 gene were down-regulated (Figure 2).In drought-sensitive genotype (YZ), there were 2325 DEGs, and 1075 and 1250 were up-regulated and down-regulated, respectively (Figure 2).GO analysis found that the DEGs of drought-tolerant genotype JM were enriched into signal transducer activity (GO:0005057), MAP kinase activity (GO:0004707), intracellular signal transduction (GO:0035556) and cellular response to abiotic stimulus (GO:0071214) while those of drought-sensitive genotype YZ were mainly involved in photosynthesis (GO:0015979), membrane protein complex (GO:0098796), ER membrane protein complex (GO:0072546), guard cell differentiation (GO:0010052), and positive regulation of response to oxidative stress (GO:1902884) (Figure 3).These results indicated that the drought-tolerant genotype activated a series of signaling pathways in response to drought stress and made adaptive adjustment, so the cell process and membrane activity were not negatively affected by drought stress, while the drought-sensitive genotype did not rapidly activate signaling transduction but activated photosynthesis and cellular process to cope with drought stress.The enriched KEGG pathway of the DEGs also showed the drought-tolerant genotype as mainly involved in signaling pathway while the drought-sensitive genotype was involved in photosynthesis and cellular activity (Figure 4).The different mechanisms of tolerant and sensitive genotypes responding to drought will contribute to develop an effective method to cope with drought stress and to assess the ability of stress tolerance.Furthermore, it is found that 377 DEGs were overlapped in drought-tolerant and sensitive genotypes, which could play fundamental roles in the regulation of drought response in wheat, including AP2/ERF (APETALA2/ethylene responsive factor), MYB and WRKY transcription factors.Additionally, 2377 and 1948 DEGs were specifically identified in genotypes JM and YZ respectively, which might be involved in genotype-specific regulatory pathway in response to drought.

Identification of the Differentially Expressed Genes
To identify the drought-responsive genes, the differentially expressed gene (DEGs) were analyzed between different treatments and different genotypes using padj < 0.05 and |log2Ratio| ≥1 as thresholds.In drought-tolerant genotypes (JM), 2754 DEGs were identified between irrigation and drought condition, of which 1152 gene were up-regulated, and 1612 gene were down-regulated (Figure 2).In drought-sensitive genotype (YZ), there were 2325 DEGs, and 1075 and 1250 were up-regulated and down-regulated, respectively (Figure 2).GO analysis found that the DEGs of drought-tolerant genotype JM were enriched into signal transducer activity (GO:0005057), MAP kinase activity (GO:0004707), intracellular signal transduction (GO:0035556) and cellular response to abiotic stimulus (GO:0071214) while those of drought-sensitive genotype YZ were mainly involved in photosynthesis (GO:0015979), membrane protein complex (GO:0098796), ER membrane protein complex (GO:0072546), guard cell differentiation (GO:0010052), and positive regulation of response to oxidative stress (GO:1902884) (Figure 3).These results indicated that the drought-tolerant genotype activated a series of signaling pathways in response to drought stress and made adaptive adjustment, so the cell process and membrane activity were not negatively affected by drought stress, while the drought-sensitive genotype did not rapidly activate signaling transduction but activated photosynthesis and cellular process to cope with drought stress.The enriched KEGG pathway of the DEGs also showed the drought-tolerant genotype as mainly involved in signaling pathway while the drought-sensitive genotype was involved in photosynthesis and cellular activity (Figure 4).The different mechanisms of tolerant and sensitive genotypes responding to drought will contribute to develop an effective method to cope with drought stress and to assess the ability of stress tolerance.Furthermore, it is found that 377 DEGs were overlapped in drought-tolerant and sensitive genotypes, which could play fundamental roles in the regulation of drought response in wheat, including AP2/ERF (APETALA2/ethylene responsive factor), MYB and WRKY transcription factors.Additionally, 2377 and ,1948 DEGs were specifically identified in genotypes JM and YZ respectively, which might be involved in genotype-specific regulatory pathway in response to drought.The expression patterns of genes will provide the crucial information for understanding their biological function [49].Thus, all the identified DEGs were further used to define clusters based on their specific expression patterns in each sample.A K-mean clustering method was conducted with the squared Euclidean distance measure, and all DEGs could be classified into 10 categories (Figure 5).The cluster I-V showed the relatively higher expression in sensitive genotype than tolerant genotype under both well-watered and drought conditions.Among them, Cluster I comprised 84 genes showing high expression in YZ_R but showed almost similar expression in other samples, suggesting they might be specifically induced by drought in sensitive genotypes.Cluster 2 with 647 genes showed low expression in JM_I and showed high expression in other samples.These genes might be involved in response to drought and GO enrichment analysis found they mainly functioned in cellular response to stress (GO:0033554), plant organ senescence (GO:0090693) and ARF (Auxin response factor) protein signal transduction (GO:0032011) and regulation of ARF protein signal transduction (GO:0032012).The cluster VI-X showed relatively high expression in the tolerant genotype compared to the sensitive genotype.Cluster 9 with 169 genes showed specifically high expression in JM_R, suggesting these genes play a crucial role in regulating tolerance to The expression patterns of genes will provide the crucial information for understanding their biological function [49].Thus, all the identified DEGs were further used to define clusters based on their specific expression patterns in each sample.A K-mean clustering method was conducted with the squared Euclidean distance measure, and all DEGs could be classified into 10 categories (Figure 5).The cluster I-V showed the relatively higher expression in sensitive genotype than tolerant genotype under both well-watered and drought conditions.Among them, Cluster I comprised 84 genes showing high expression in YZ_R but showed almost similar expression in other samples, suggesting they might be specifically induced by drought in sensitive genotypes.Cluster 2 with 647 genes showed low expression in JM_I and showed high expression in other samples.These genes might be involved in response to drought and GO enrichment analysis found they mainly functioned in cellular response to stress (GO:0033554), plant organ senescence (GO:0090693) and ARF (Auxin response factor) protein signal transduction (GO:0032011) and regulation of ARF protein signal transduction (GO:0032012).The cluster VI-X showed relatively high expression in the tolerant genotype compared to the sensitive genotype.Cluster 9 with 169 genes showed specifically high expression in JM_R, suggesting these genes play a crucial role in regulating tolerance to drought.Then, GO analysis of these 169 genes found that they mainly enriched the positive regulation of cellular response to oxidative stress (GO:1900409), positive regulation of response to ROS (GO:1901033), stomatal movement (GO:0010118) as well as auxin-activated signaling pathway (GO:0009734), which proved their roles in regulating wheat defense systems to tolerate drought stress.
drought.Then, GO analysis of these 169 genes found that they mainly enriched the positive regulation of cellular response to oxidative stress (GO:1900409), positive regulation of response to ROS (GO:1901033), stomatal movement (GO:0010118) as well as auxin-activated signaling pathway (GO:0009734), which proved their roles in regulating wheat defense systems to tolerate drought stress.drought.Then, GO analysis of these 169 genes found that they mainly enriched the positive regulation of cellular response to oxidative stress (GO:1900409), positive regulation of response to ROS (GO:1901033), stomatal movement (GO:0010118) as well as auxin-activated signaling pathway (GO:0009734), which proved their roles in regulating wheat defense systems to tolerate drought stress.

Analysis of RNA Editing Sites
RNA editing is a process that occurs in the RNA molecular base change or modification when transcribed, which is one of most important mechanisms regulating gene expression and enriching genetic information at the post-transcription level [50].A larger number of studies have reported that RNA editing not only controls plant organ formation, growth, and development, but also plays an indispensable role in the response to diverse stresses [51,52].The RNA-seq data provides a resource to identify the RNA editing sites at the whole transcriptome level.Using methods described in the Material and Methods section, we detected the RNA editing sites between water and drought condition to identify the drought-induced RNA editing sites (Table 2).In total, 32 drought-responsive RNA editing sites in 22 genes were found, of which 30 were found in drought-tolerant genotype JM and 2 were in drought-sensitive genotype YZ.TraesCS6B01G079400.1 have 4 editing sties, of which 3 were in genotype JM and 1 in YZ, followed by TraesCS4B01G293600.3 showing 3 editing sites and another 6 genes with 2 editing sites.The remaining 13 genes owned one sites.Function annotation of these genes with RNA editing sites found that they included the transcription factor gene, such as MYB (TraesCS6B01G012800.1) and bHLH (TraesCS5A01G279200.1), kinase proteins as well as histone and plasma membrane, suggesting RNA editing might function as the key regulator in activating the processes and pathways of drought tolerance.Furthermore, 13 sites were found in UTR regions and the remaining 19 sites were in coding regions, of which 9 sites were edited at the third position of the codon, 8 at the second position and 2 at the first position.The 10 edited sites occurred at the first and second position of the codon caused the alteration of amino acid, which could be considered as candidates for further functional study.Finally, a total of 7 types of base change were introduced by RNA editing in these 32 sites, of which C to T mutation is the most abundant type with the value of 10 times, followed by T to C and G to A with the number of 5, C to A with the number of 4, G to T and G to C with the number of 3 as well as A to G with the number of 2. This result showed that the transition (22) was significantly higher than transversion (10) in these identified RNA editing sites in wheat, which was consistent with the previous reports in the plastomes of einkorn wheat and Aegilops tauschii.L. [51,52].
Note: a JM = the drought-tolerant genotype Jimai No. 47; YZ = the drought-sensitive genotype Yanzhan No. 4110; b the position represents the position of the edited site at the mRNA molecular; c I1 and I2 = the two samples under irrigation condition whereas R1 and R2 = the two samples under rain-fed condition.The capital = the abbreviation of nucleotide; the number = the number of mapped reads; d "-" = the editing sites located in UTR region.

Table 1 .
Summary of RNA-seq and mapping data.

Table 2 .
The drought-responsive RNA editing sites identified in wheat.