Genomic and Transcriptomic Dissection of the Large-Effect Loci Controlling Drought-Responsive Agronomic Traits in Wheat

Abstract

Drought tolerance is one of the most important targets for wheat breeding. Previous population genetics studies have uncovered 20 large-effect quantitative trait loci (QTLs) that contribute to stress-responsive agronomic traits. Here, we identified 19,035,814 single nucleotide polymorphisms and 719,049 insertion/deletion variations in the genomes of two popular winter wheat cultivars, Lu-Mai 14 and Han-Xuan 10, using a whole-genome re-sequencing assay. There were 4972 loss-of-function mutations carried by protein-coding genes, such as CCA1/LHY, AGO1, ABI3/VP1, EIN3, TPP, and ARFs. We carried out a time-course abscisic acid (ABA)-treatment experiment and profiled 61,251 expressed genes in the roots using a strand-specific RNA sequencing approach. A large number of genes showed time-point specific and/or cultivar-preferential responsive expression patterns. Gene ontology enrichment analysis revealed that ABA-responsive genes were associated with stress-related functions. Among the 20 QTLs, we uncovered 306 expressed genes with high- and/or moderate-effect variations and 472 differentially expressed genes. Detailed analysis and verification of the homozygous genomic variations in the candidate genes encoding sulfotransferase, proteinase, kinase, nitrate transporter, and transcription factors suggested previously unexpected pathways associated with abiotic stress responses in wheat.

1. Introduction

Wheat (Triticum aestivum L., AABBDD, 2n = 42) is one of the major cereal crops, providing approximately 19% of global dietary energy [1]. Exposure to drought or other abiotic stresses during the wheat life cycle seriously restricts crop productivity [2,3,4]. With the increasing global demand for food, breeding high-yielding wheat cultivars tolerant to abiotic stresses is of primary importance [5,6]. Based on current advances in technologies and theories of molecular biology [7,8,9,10], marker-assisted selection utilizing genetic linkage between phenotypes and agronomic traits and target genes or quantitative trait loci (QTLs) is highly effective for crop improvement. Systemic profiling of trait-linked markers associated with agronomic traits that respond to drought and abiotic stresses will accelerate the breeding process of wheat cultivars with extensive adaptability, high yield, and stress endurance [11,12,13].

The hexaploid winter wheat cultivars Lu-Mai 14 (LM14), Han-Xuan 10 (HX10), and their derived varieties have been widely cultivated in North China [12]. These two cultivars represent distinct types of environmental adaptations in breeding selection. LM14 is a high-yielding cultivar with low drought tolerance and is mainly planted in well-irrigated regions. In contrast, HX10 has a moderate-level yield, high drought tolerance, and is widely cultivated in arid and semi-arid regions. In our previous study [12], 20 large-effect QTLs associated with drought and/or heat stress were identified using linkage-based QTL mapping analysis of 150 doubled haploid lines and a genome-wide association study (GWAS) of 277 wheat accessions phenotyped across 30 environments.

Abscisic acid (ABA) is a phytohormone that can increase plant tolerance to multiple abiotic stresses by directing the stimulation of the stomatal aperture, accumulating osmocompatible solutes, and regulating plant development [14,15]. ABA stimulation can trigger reactive oxygen species (ROS) production, which reduces the activities of ABI1, ABI2, and MAP kinases [16,17]. In wheat, the application of exogenous ABA can enhance dormancy, stress endurance, and resistance to pre-harvest sprouting [18,19,20]. In higher plants, ABA molecules are sensed by ABA receptor proteins [21,22], leading to the inactivation of type 2C protein phosphatases (PP2Cs) and activation of sucrose-non-fermenting 1-related protein kinase 2 (SnRK2) [23]. Subsequently, SnRK2 binds to transcription factors, such as ABI3 and ABI5, and interacts with ABA-responsive element (ABRE) motifs on the promoter regions of downstream genes to trigger ABA-responsive transcription processes, mediate root elongation repression, and promote stomatal closure [15,24,25,26,27,28,29,30].

Single nucleotide polymorphism (SNP) and insertion/deletion (InDel) variations are high-density genomic variations, which can lead to phenotypic variations; therefore, both are often used to develop molecular markers [31,32]. SNP molecular markers are important auxiliary breeding markers in crops because of their abundant loci, high sensitivity, reproducibility, accuracy, strong specificity, and ease of use. SNP molecular markers are used in applications such as high-density genetic linkage mapping, genotyping, and interspecific genetic relationship analysis [32,33,34]. InDel markers are moderate polymorphism differences markers, with a simple design and stable genotyping. InDel markers are successfully used in the genetic research of wheat, using regular primers, basic PCR instruments for amplification, polyacrylamide gel (PAGE) and high-resolution melting (HRM) technology for gene analysis [35,36,37]. Strand-specific RNA sequencing (ssRNA-seq) can eliminate the expression interference from complementary strand transcripts and make gene quantification more accurate [38,39]. It plays an important role in new transcript discovery, genome annotation, and expression profile analysis [40].

Resequencing and ssRNA-seq are high-throughput genomic variation identification methods. To systemically construct the molecular genomic map of the two widely cultivated wheat cultivars, LM14 and HX10, we used a re-sequencing assay to identify SNPs and InDels on a genome-wide basis. Since numerous variant loci can be obtained through re-sequencing a joint analysis of re-sequencing and large-effect QTLs was carried out to uncover molecular markers of importance, find causal variants or genes with important mutations, and narrow the range of candidate genes. Moreover, we performed a time-course ABA treatment experiment and ssRNA-seq analysis to identify ABA-responsive genes and monitor their expression dynamics during ABA processing in the two wheat cultivars. The genes of the QTL regions were analyzed for differential expression in the transcriptome, to verify the relationship between ABA-responsive genes and drought resistance. Finally, we assessed the drought-related candidate genes in terms of both gene structure variation and gene expression. This study provides a new data source for drought resistance molecular markers and drought-related genes in wheat. It provides a theoretical and experimental basis for the identification of drought-resistant wheat varieties and their breeding.

2. Materials and Methods

2.1. Plant Materials and Growth Conditions

Seeds of the wheat cultivars LM14 and HX10 were placed in Petri dishes (diameter 15 cm) lined with wet paper and germinated for 3 days in an incubator at 25 °C; three biological replicates were prepared. All seeds were pre-sterilized to prevent possible microbial contamination. The seeds were surface-sterilized with 70% ethanol for 3 min, treated with 2% sodium hypochlorite (NaClO) for 5 min, and rinsed thrice in sterile distilled water [41]. After germination, the seedlings were cultivated in Hoagland solution in a chamber under a 16/8 h light/darkness cycle and 22/18 °C, respectively, until the three-leaf stage. For ABA and polyethylene glycol (PEG) treatments, seedling roots were immersed in Hoagland solutions containing either 200 μM ABA or 20% PEG6000, respectively. Tissues were harvested at four time points (0, 6, 12, and 24 h post treatment (hpt)), immediately frozen in liquid nitrogen, and stored at −80 °C prior to RNA extraction.

Following the protocols described previously [42,43], we determined the relative water content (RWC) of the third leaf from the top and the primary roots. The same tissues were used to determine ion leakage using a previously reported method [44,45].

2.2. DNA and RNA Extraction and Library Construction

Genomic DNA was extracted using the Hi-DNAsecure Plant Kit (Tiangen Biotech, Beijing, China) and quantified using a Qubit 4 (Thermo Fischer Scientific, Waltham, MA, USA). The TruSeq™ DNA Kit (Illumina, San Diego, CA, USA) was used to construct the DNA sequencing library using genomic DNA that had been sheared with a Bioruptor Pico instrument (Diagenode, Seraing, Belgium) according to the manufacturer’s instructions.

Total RNA was isolated using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer’s instructions. RNA concentration was measured using a NanoDrop 2000 (Thermo Fischer Scientific, Waltham, MA, USA), and RNA integrity was evaluated using the RNA 6000 Nano kit on a 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). RNA sequencing libraries were prepared using the UltraTM Directional RNA Library Prep kit (New England Biolabs Inc., Ipswich, MA, USA) according to the manufacturer’s instructions.

For ABA treatment, seedling roots were immersed in Hoagland solutions containing 200 μM ABA. The root tissue RNA extracted from the two cultivars was repeated thrice at four time points (0, 6, 12, and 24 hpt). Moreover, the DNA of two cultivars (LM14 and HX10) was extracted directly. In total, 24 RNA strand-specific RNA sequencing libraries and 2 DNA re-sequencing libraries were constructed. The quality and quantity of each library were assessed using an ABI 7300 (Applied Biosystems, Foster City, CA, USA) and Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) according to the manufacturer’s instructions. DNA and RNA libraries were sequenced using a HiSeqX10 sequencing system (Illumina, San Diego, CA, USA) with a 150-cycle paired-end sequencing protocol. The fastq-formatted RNA-seq datasets are publicly available under accession number PRJCA002233 in the Genome Sequence Archive database (http://gsa.big.ac.cn/ (accessed on 30 June 2020)) of the BIG Data Center.

2.3. Bioinformatics Analysis for RNA-seq Experiment

The wheat genome sequence and annotation files were downloaded from the Ensembl Plants database under the release IWGSC Ref1.0 and IWGSC/INSDC Assembly GCA_900519105.1 [7,46]. After removing low-quality reads and adapter contaminants, we aligned the reads to the reference genome using HISAT2 [47,48,49]. Using a pipeline composed of Cufflinks, SAMtools, HTseq-count, and DESeq2 with strand-specific parameters [50,51,52,53,54], we assembled the transcription units and calculated the uniquely mapped read counts and transcripts per million reads (TPM) values of the annotated genes. The expressed genes were filtered using the criterion of TPM > 1 in at least one sample.

2.4. Genotyping Analysis

The sequencing reads were mapped to the wheat reference genome using BWA software (v.0.7.17) with default options [55]. Then, duplicate reads were removed using the GATK (v.4.0.10.0) MarkDuplicates function [56]. SAMtools (v.1.8) mpileup was used to generate raw variants (SNPs and InDels) with -ugf option [57]. BCFtools (v.1.9) was then used for variant calling with the -vm option [58], and VCFtools (v.0.1.13) was used for downstream filtration [59]. The effects of these SNPs/InDels were predicted using SnpEff (v.4.3) [60].

2.5. Experimental Verification by Sanger Sequencing Assay

Using the aforementioned genomic DNA as a template, PCR products were purified using the Axygen gel extraction kit (Axygen, Union City, CA, USA). The PCR products were subcloned using the pEASY^®-Blunt Zero vector (Transgen, Beijing, China) for Sanger sequencing to identify the SNPs in the two cultivars. Primer sequences are listed in Table S1.

2.6. Statistical Analysis

Data were analyzed using Systat version 13.2 (Systat software, Inc., Chicago, IL, USA). Student’s t-test was performed with SPSS software version 25.0 for windows (SPSS Inc., Chicago, IL, USA).

3. Results

3.1. Whole-Genome Re-Sequencing of LM14 and HX10 Discovered a Large Number of SNP/InDel Variations in Wheat

In total, 801,126,119 and 714,694,329 sequencing reads were obtained in LM14 and HX10, respectively. Of those, 798,608,429 (99.7%) and 709,880,835 (99.3%) reads could be aligned to the reference wheat genome assembly IWGSC RefSeq 1.0, covering approximately 8.23× and 7.32× average depths, respectively. We identified 19,035,814 SNPs and 719,049 InDels. The A- and B-subgenomes of wheat contained 7,448,150 and 9,562,457 SNPs/InDels, respectively, whereas the D-subgenome had 2,485,967 SNPs/InDels (Figure 1A). The A- and B-subgenomes contained more SNPs/InDels than the D-genome, which is consistent with previous studies [61,62]. Of these, 19,216,255 SNPs/InDels were identified in the intergenic regions. The remaining 538,608 were located in intragenic regions and comprised 237,899 exons and 300,709 introns. According to the VCF annotation [59], the SNPs/InDels leading to serious structural destruction of the open reading frame (ORF) and rare amino acid changes were classified as high-effect variations; those causing UTR elongations/truncations and missense variations were classified as moderate-effect variations; and the remaining ones without structural changes in the ORF or exon, such as synonymous amino acid change, were considered as low-effect variations (Table S2). In the LM14 and HX10 genomes, we identified 4972 high-effect, 129,002 moderate-effect, and 121,433 low-effect SNP/InDel variations (Figure 1B).

The 4972 high-effect SNPs/InDels causing serious structural destruction of the ORF could be considered as loss-of-function variations of their hosting genes. By orthologous search of Arabidopsis and rice, we identified a group of regulatory genes with known biological functions carrying high-effect variations in the two wheat cultivars. There were 307 receptor kinases and 258 NBS-LRR proteins encoded by genes carrying high-effect SNPs/InDels, which is consistent with the fact that the exon structures of defense genes are highly variable [63]. Additionally, some gene-encoding proteins controlling multiple pathways were found in the list, such as Circadian Clock Associated 1 (CCA1/LHY, TraesCS7B01G188000, TraesCS1B01G227400), Time For Coffee (TraesCS2A01G197800), Argonaute-1 (TraesCS2A01G258100, TraesCS2A01G258200, TraesCS6D01G403900), Clavata-3 (TraesCS1A01G393200), the homeobox protein Knotted-1 (TraesCS7A01G511800), Ethylene Insensitive 3 (EIN3, TraesCS3B01G428000), Abscisic Acid Insensitive 3 (ABI3/VP1, TraesCS6D01G384000, TraesCS6B01G440000, TraesCS5B01G073600, TraesCS5A01G363700, TraesCS4A01G457200, TraesCS2A01G009100, TraesCS1A01G075800), Root Hair Defective 6 (TraesCS7A01G229900), the histone-lysine N-methyltransferase ASH2 (TraesCS2B01G015500), Trehalose-6-phosphate Synthases (TPSs, TraesCS5A01G203500, TraesCS1D01G065300, TraesCS1A01G064800), Trehalose 6-phosphate Phosphatase (TPP, TraesCS6A01G301800), EID1-like F-box protein 3 (TraesCS3A01G331200), the auxin transport protein BIG (TraesCS2A01G038500), the phosphate transporter PHO1 (TraesCS6B01G418600), and the potassium transporters (TraesCS7A01G209900, TraesCS5A01G233200, TraesCS2D01G118400, TraesCS1D01G048200, TraesCS1B01G061200) (Table S3). The biological functions of these genes in wheat are unknown.

Of the 20 QTLs associated with abiotic stress tolerance, we identified 79 high-effect, 1733 moderate-effect, and 1325 low-effect SNP/InDel variations (Figure 1B). The chromosomal distributions of gene density, SNP/InDel density, and moderate/high-impact SNP/InDel density are shown in Figure 1C. These results show a large number of genetic variations in the two wheat cultivars. QTL analysis of the doubled haploid population could effectively narrow down the genetic variations associated with the agronomic traits that respond to abiotic stresses.

3.2. Identification of the Expressed Genes in LM14 and HX10 Roots Using RNA-seq

To measure the morphological features associated with the drought stress response of the two cultivars in the greenhouse, we grew the seedlings of the LM14 and HX10 at the 3-leaf stage under water culture conditions and treated the roots with 20% PEG-6000 and 200 μM ABA according to the method used in previous studies [64,65]. After a 5-day treatment with ABA and PEG, LM14 showed yellow leaf phenotypes under both treatments, whereas HX10 presented a relatively slightly stressed phenotype (Figure 2A). The relative water content of HX10 leaves and roots was significantly higher than that of LM14 under ABA and PEG treatments, whereas there was no significant difference between the two cultivars under the control treatments (Figure 2B). We also measured the stress status of the samples by determining the ion leakage. Ion leakage levels were significantly lower in the leaves and roots of HX10 than in those of LM14 (Figure 2C). These results confirmed the different phenotypes between LM14 and HX10 under greenhouse conditions. Moreover, it revealed that the drought sensitivities of the two cultivars are associated with ABA signaling sensitivities.

To monitor the gene expression network underlying the ABA response process of the wheat cultivars with different stress sensitivities, we treated the roots of the 3-leaf LM14 and HX10 seedlings with 200 μM ABA and isolated the root samples at 4 time points (6 hpt,12 hpt, 24 hpt, and 0 hpt (untreated)) with 3 replicates. We constructed strand-specific RNA sequencing libraries of the 24 samples and sequenced them using a 150-cycle paired-end sequencing protocol.

After removing the low quality reads and adapter contaminants, 1,223,375,840 reads (85.8% average mapping rate and 61.2% average unique mapping rate) were aligned to the wheat reference genome release IWGSC RefSeq v1.0 (Table S4). We calculated the TPM values of the annotated genes using uniquely mapped reads. The gene expression levels were highly correlated among the biological replicates and principal component analysis results indicated high reproducibility of the RNA-seq datasets (Figure 3A,B). We identified 61,251 expressed genes (TPM > 1) in at least one sample (Figure 3C; Table S5). During analysis of the 41,509 homeologous genes from the A-, B-, and D-subgenome, 37,158 (90.4%) genes were expressed in all three subgenomes, 3092 genes were expressed in two of the three subgenomes, and the remaining 1259 genes were expressed in A-, B- or D-subgenomes (Figure 3D).

3.3. Time-Point Specific Expression Pattern of ABA-Responsive Genes in LM14 and HX10 Roots

We performed differential expression analysis to profile the expression patterns of the ABA-responsive genes. To identify the time-point specific responsive genes, we compared the expression levels of the genes in 6-, 12-, and 24-hpt samples with those in untreated samples. We identified 40,010 and 38,149 significantly differentially expressed genes (DEGs) in LM14 and HX10, respectively. A total of 30,506 (64.0%) genes were identified in both cultivars (Figure 4A). Next, we searched for cultivar-preferential responsive genes by comparing the expression levels of genes between LM14 and HX10 at 0, 6, 12, and 24 hpt with ABA. In total, we identified 19,696 significantly differentially expressed genes between LM14 and HX10, including 10,298 genes that also showed time-point specific responsive expression patterns in the two cultivars. Moreover, 5769 differentially expressed genes were identified in the untreated samples between LM14 and HX10, suggesting that, in addition to the background expression difference, the cultivars may have a group of DEGs endogenously associated with ABA signaling and downstream pathways. The statistics of the time-point specific and cultivar-preferential responsive genes are shown in Figure 4B.

We further analyzed the expression patterns of DEGs by clustering their normalized TPM values. The 28,134 genes (TPM > 5 in at least one sample) were classified into 120 clusters. Figure 4C–H shows several representative expression patterns: sustained induction in both cultivars (Figure 4C), sustained induction in LM14 (Figure 4D), sustained induction in HX10 (Figure 4E), background expression difference in the two cultivars (Figure 4F), sustained repression in HX10 (Figure 4G), and sustained repression in LM14 (Figure 4H). The expression patterns are shown in Figure S1.

3.4. Functional Categorization Suggesting the Specific Biological Processes of ABA-Responsive Genes with Specific Expression Patterns

The time-course experiment with two cultivars enabled us to classify four representative types of genes with specific ABA-induced/repressed expression patterns: (1) 20,208 ABA-responsive genes in both cultivars without a cultivar-preferential expression pattern, which were considered as the general ABA-responsive genes (Figure 5A); (2) 10,298 genes showing differential expression levels not only between ABA treated and untreated samples but also between cultivars, which were classified as cultivar-preferential ABA-responsive genes (Figure 5B); (3) 11,166 DEGs between cultivars at 0 hpt, reflecting the endogenous expression difference between cultivars (Figure 5C); (4) 8530 genes responsive to exogenous ABA application and expressed at similar levels between cultivars at 0 hpt but differentially expressed between strains under ABA treatment (Figure 5D). The biological functions of these genes were explored by gene ontology (GO) enrichment analysis. All the GO enrichment results were listed in Table S6 (p-value < 0.05). The 20 representative GO terms identified in each group are shown in Figure 5. The GO terms associated with stress-response, such as “response to stress”, and “oxidoreductase activity”, were found in all four groups, indicating that the molecular profiling analysis was consistent with the morphological and physiological results of ABA time-course treatment experiments. In the general ABA-responsive genes, the significantly enriched GO terms were “microtubule”, “microtubule-based process”, “DNA replication”, and “glycolytic process”. In the cultivar-preferential ABA-responsive genes, significantly enriched GO terms were associated with “nicotianamine biosynthetic process”, “transaminase activity”, “apoplast”, and “cellulose synthase activity”. The endogenous cultivar-preferential genes were enriched in “nicotianamine biosynthetic process”, “ammonium transport”, “cell growth”, “chitin binding”, “cell wall”, and “calmodulin binding”. The genes responsive to exogenous ABA application were associated with multiple processes, such as “sodium ion transport”, “nitrate assimilation”, “potassium ion transport”, “glutamate metabolic process”, “sucrose synthase activity”, “trehalose biosynthetic process”, and “isoprenoid biosynthetic process”. These genes are likely to act as downstream regulators of multiple biological processes directed by ABA signaling pathways in wheat. The four groups of genes with distinct expression patterns were enriched in the specific GO terms, suggesting time-point specific and/or cultivar-preferential downstream pathways involved in the ABA and abiotic stress responses in wheat. Detailed clarification of the pathways warrants further study.

3.5. Screening Candidate Genes Associated with Abiotic Stress Responses in Wheat

The genome annotation located 2469 genes in the previously identified 20 large-effect QTLs associated with drought and/or heat stresses. Variations in genomic structures or expression differences of several master genes at these loci may play causal roles in abiotic stress responses. Using whole-genome sequencing and time-course RNA-seq analyses, we identified 1540 expressed genes, 411 carrying SNPs/InDels in their coding sequence (CDS) region, 306 genes with high- and/or moderate-effect variations, and 472 DEGs between the two cultivars on the loci (Figure 6).

On the basis of their expression patterns, SNP/InDel variations, and annotated molecular functions, we selected 12 candidate genes. TraesCS2B01G200100 encodes a sulfotransferase that can transfer a sulfonate group from 3-phosphoadenosine-5′-phosphosulfate to a hydroxyl/amino group on substrates such as flavonoids, brassinosteroids, and salicylic acid, and is known to be a functional regulator of plant root development and ABA-mediated abiotic responses [66,67,68]. In our dataset, we found that the gene was highly induced by ABA treatment, confirming that it is a stress-responsive gene in wheat (Tables S5 and S7). The re-sequencing analysis identified a homozygous single nucleotide guanine deletion in LM14, leading to a frameshift loss-of-function mutation in its ORF. The mutation and its adjacent SNP were confirmed by Sanger sequencing in LM14 and HX10 (Figure 7A). TraesCS2B01G205300 encodes a serine proteinase that is highly induced by ABA. We identified a homozygous C-to-G stop codon gain SNP in LM14 that caused an ORF truncation and a loss-of-function mutation (Figure 7B). TraesCS2B01G216800 encodes a serine/threonine-protein kinase, whose expression is repressed by ABA treatment. We found a homozygous T-to-G missense SNP and a homozygous AGC in-frame insertion in LM14 (Figure 7C). TraesCS2A01G538500, encoding a helix-loop-helix (bHLH) transcription factor, was highly induced by ABA treatment. We identified eight homozygous SNPs on the gene body in LM14, and the C-to-T, A-to-C, and G-to-T SNPs led to missense mutations in its ORF (Figure 7D). TraesCS5A01G430700 encoding peroxidase was repressed by ABA treatment. The gene carrying the A-to-G SNP resulted in a missense mutation in HX10 (Figure 7E). TraesCS6A01G217100 encodes an expansin protein associated with cell growth. Its expression levels were high in LM14 and HX10 roots and were repressed by ABA treatment. We identified a thymine insertion that caused a loss-of-function mutation with frameshift and stop codon gain (Figure 7F). TraesCS6A01G209900, which encodes a high-affinity nitrate transporter (NRT) protein, is an important regulator of plant development and biomass accumulation. The expression of NRT was repressed by ABA treatment. In HX10, NRT carried a homozygous C-to-G missense SNP (Figure 7G). TraesCS2B01G167900 is a bZIP transcription factor-encoding gene that is highly induced by ABA treatment. We identified a homozygous C-to-A missense SNP, which was confirmed using Sanger sequencing (Figure 7H). In addition to the SNP and/or InDel variations, we identified 472 genes that were differentially expressed between LM14 and HX10. For example, the expression levels of the microsomal glutathione S-transferase (TraesCS4D01G042000) and the sulfite reductase hemoprotein beta-component (TraesCS1A01G323400) were specifically induced under ABA treatment in HX10 compared with those in LM14, whereas the FAD-binding berberine gene (TraesCS2B01G172300) and the IgA FC receptor (TraesCS2B01G205100) were specifically induced by ABA in LM14, but not in HX10 (Figure 7I–L). Our analyses systemically profiled the genomic variations and expression patterns of genes on the large-effect loci associated with ABA-mediated abiotic responses in wheat.

4. Discussion

4.1. Genomic Features Characterized by Re-Sequencing in Wheat

As an allohexaploid Triticeae plant, wheat has a large genome with abundant genomic variations [7]. Previous studies have identified and screened SNP and InDel variations in landraces, germplasm resources, modern cultivars, and purified inbred lines using genotyping by sequencing (GBS) [69,70], kompetitive allele-specific PCR (KASP) [32,71], SNP microarray [72,73], and simple sequence repeat (SSR) detection assays [74,75]. With the emergence of high-quality assembly of the wheat genome, whole-genome re-sequencing assays are becoming one of the most robust high-throughput approaches to identify high-density variations at the single-nucleotide level [61,76]. In our study, we identified a large number of SNPs and InDels using re-sequencing with 7–8× coverage in two widely cultivated cultivars with distinct environmental adaptations. The A- and B-subgenomes contained more variations than the D-subgenome. Of the 41,509 A-, B-, and D-subgenome homeologous genes, 37,158 (90.4%) genes were expressed in all three subgenomes. These genomic features are consistent with those of previous studies on wheat [77].

4.2. The Phenotypic and Physiological Differences in Wheat Based on Differences in ABA Sensitivity

ABA, as a very important signal in drought stress, has a significant impact on the drought resistance of plants. Consistent with the research of Zhao et al., the yellowing of leaves treated with ABA and PEG could be attributed to the chlorophyll degradation caused by the accumulation of ABA [78,79]. The decrease in the relative water content in the leaves and roots indicated increased water loss caused by transpiration. The increase in ion leakage indicated that continuous stress led to the collapse of the antioxidant system, and subsequently, plasma membrane damage. The response to PEG and ABA was similar between the two wheat cultivars. However, phenotypic and physiological differences indicated that there are remarkable differences in the degree of damage caused by the differences in the variety. These results are in line with earlier reports [80,81,82]. Chlorophyll degradation and transpiration water loss inhibit photosynthesis. These changes in phenotype and physiology indicate that exogenous ABA and PEG treatments induce various stress-related signals that are transported long distances from the root to the leaves. There are significant differences in the accumulation of ABA in different genotypes of wheat under water stress. The stronger the drought resistance, the lesser the ABA accumulated [83,84,85]. Drought resistant varieties accumulate less ABA, while drought sensitive varieties accumulate more ABA [83,84]. The amount of ABA accumulated affects the growth of plants, which is consistent with the results of this study. The capacity to accumulate ABA and ABA sensitivity is heritable genotypic variations; therefore, it is important to elucidate the reasons underlying these genotypic differences [83,86].

4.3. The Gene Expression Patterns Responsed to ABA Treatment in Wheat

Our time-course transcriptome analysis uncovered a large number of ABA-responsive genes and their expression patterns. We used the K-means clustering algorithm to divide the DEGs into 120 patterns (Figure S1, Table S5). These gene expressions and clustering information provide valuable basic data for studying ABA and drought response in wheat.

GO enrichment analysis found the ABA-responsive genes associated with stress-response functions, which were consistent with the morphological and physiological features of ABA response processes [79,82]. Moreover, we also identified a group of cultivar-preferential genes that were significantly associated with apoplast and nicotianamine biosynthesis. It is known that stress-induced ABA and ABA-glucose as its conjugated form can be secreted into the apoplast [87]. Nicotianamine can promote the transportation activity of apoplast to load biosynthesis products for secretion, such as metals, leading to the accumulation of apoplastic ROS (including OH, O²⁻, H₂O₂, ¹O₂) in cells [88]. These results suggested that the transportation activity of apoplast and the deposition processes of ABA in apoplast may be a distinguished feature of a drought-resistant wheat cultivar. In addition, the cultivar-preferential genes responsive to exogenous ABA application showed the enriched function terms in multiple biological processes, such as transportation of nitrate, potassium, and sodium, sucrose synthase activity, and isoprenoid biosynthesis, suggesting that nutrient balance, sucrose synthesis, and secondary metabolites also play important roles in stress tolerance in the wheat cultivar. It is highly warranted to make these hypotheses and their values in marker-assisted breeding clear.

The candidate genes obtained from the combined analysis of transcriptome and large-effect QTLs are shown in Figure 7 and Table S3. The functions of some genes remain unclear; however, some of them are clearly associated with drought resistance. The expression of thioredoxin (TraesCS5D01G454800) and basic region/leucine zipper motif (bZIP) transcription factor (TraesCS6D01G160200) (Figure S2B,F) did not differ significantly between the two cultivars at 0 hpt (mock) and 24 hpt; however, both were strongly induced at 6 hpt and their expression decreased subsequently. The induction intensity of HX10 was significantly higher than that of LM14, indicating that these genes belong to the group of fast response genes and could be located upstream of the signal transduction. Thioredoxin is widely believed to protect the activity of cell membrane, form disulfide through oxidation, reduce lipid peroxidation and play a role in the tolerance of plants to oxidative stress [89,90,91]. This implies that HX10 functions through thioredoxin to regulate the activity of antioxidant enzymes, inhibit the oxidation of cell membranes, and reduce their sensitivity to abiotic stress [92,93,94,95]. bZIP, as a transcription factor, is frequently reported to participate in the ABA and drought resistance pathways. The analysis of the bZIP genes in the wheat gene family revealed that TabZIP (Traes_7AL_25850F96F.1) can improve drought, salt, and heat tolerance, reduce the peroxide accumulated by ROS in plants under stress, and induce the high expression of antioxidant enzymes related to ROS pathways, such as APX and CAT [96]. wlip19 was induced by both ABA and drought, however, it was confined to the leaves and not observed in the roots. It is speculated that the bZIP genes in wheat may have the differentiation of gene function that can improve the drought tolerance of plants [97]. TabZIP174 and TabZIP60 were induced by PEG and ABA, and the transgenic Arabidopsis of both genes showed significantly increased drought resistance and soluble sugar content in drought, consistent with the sugar synthesis related pathway in the Figure 5D, suggesting that the increase of the transcription factor probably affected the downstream regulation of sugar metabolism [98,99]. The above two genes were subjected to the same pattern of induced responses, differing only in degree, and it is speculated that such genes may have a dosage effect.

The heat shock transcription factor HSF (TraesCS2B01G232000) (Figure S2C) exhibits heat, drought, salt, and cadmium tolerance in Arabidopsis, maize, rice, and other plants [100,101,102,103,104]. There were 78 HSF genes in the wheat genome, which can be classified into three classes, A, B, and C, and which were able to regulate the expression of the heat shock proteins. In recent years, HSFs have been gradually studied in wheat, and it has been reported that TaHSFA2e-5D increases heat and drought resistance with an ABRE-binding cis-acting element upstream, implying that it is regulated by ABA signaling [105]. The expression patterns following PEG and ABA treatments are almost identical. The expression reaches a maximum at 12 hpt and decreases at 24 hpt [106]; however, the gene expression increased continuously in our study. This difference could be attributed to the differences in the varieties due to the use of the heat- and drought-resistant cultivar TAM107 in the previous study [106]. The time taken to reach the peak expression could be different for different cultivars of wheat; however, the overall trend remained the same. Assessing the expression of this gene for longer periods is necessary for future experiments. The two genes HSF gene(TraesCS2B01G232000) and basic helix-loop-helix (bHLH) DNA binding superfamily gene (TraesCS2A01G538500) (Figure S2C,G) had similar expression patterns and were almost unchanged in LM14, while they were strongly induced in HX10. The expression differences following induction were obvious, indicating that two genes were induced only in the drought-resistant cultivar and belonged to cultivar-specific drought-responsive genes.

AsA acts as a reactive oxygen scavenger, scavenging large amounts of reactive oxygen species produced during growth and stress, affecting photosynthetic efficiency and plant stress resistance [107]. GDP-L-galactose phosphorylase 1 (TraesCS4B01G252800) (Figure S2H) is the rate-limiting enzyme of the L-galactose pathway. OsGGP knockout resulted in an 80% reduction of AsA and affected AsA levels as a key gene [108,109]. Previous studies have found that drought decreases AsA levels in wheat leaves, and we provide a perspective on root AsA anabolism through the transcriptome, hypothesizing that HX10 is relatively less damaged due to the rapid response of the GGP gene, consistent with the physiological results [110]. The two genes in Figure S2C,G were both similar and different from those in Figure S2D (Ran GTPase-activating protein 2 (TraesCS4B01G252900), H. Two genes in Figure S2A (Dehydrin (TraesCS6B01G383600)) and E (Lipase (TraesCS2B01G160100)) were significantly more expressed after induction in LM14 than in HX10. In many past studies, dehydrin improved drought tolerance in plants, and the more resistant varieties had higher dehydrin expression, in contrast to our results [111,112]. This could be attributed to the analysis at the protein level in the earlier reports, while we analyzed at the transcriptional level. The function of dehydrin protein is affected by the phosphorylation level and post-transcriptional modifications; therefore, further studies will be needed.

4.4. Identification of Loss-of-Function Variations of Master Genes in Wheat

Genetic variations included high-effect (4972) and moderate-effect (129,002) SNPs and InDels and resulted in elongations/truncations of ORFs and exons and missense changes of amino acids in the two cultivars. These include the known master genes regulating development and stress responses, such as CCA1/LHY, AGO1, ASH2, ABI3/VP1, EIN3, TPS, TPP, NBS-LRRs, and ARFs. The phenotypes caused by structural and missense variations in these genes may be weakened in wheat by the dosage compensation effects of their homologous genes encoded by subgenomes. The potential functions of these genes can be further elucidated using overexpression transformation and CRISPR-Cas genome editing experiments [113,114,115,116].

4.5. Developing Molecular Markers for Functional Verification and Wheat Breeding

Recent improvements in high-throughput platforms and bioinformatics analysis tools, such as GWAS [11,12,117], BSA [118,119] and QTL analyses [120,121], have facilitated the advancement of genetic approaches for assigning biological functions to a certain genomic region(s) in wheat. However, performing a fine mapping experiment for identification of the causative gene(s) in large segregating populations is often labor-intensive and inefficient in wheat because of the nature of its large and allohexaploid genome. In our study, we identified structural variations and expression patterns of the genes on the large-effect QTLs by re-sequencing and time-course transcriptome profiling at the specific response process and identified the candidate genes associated with the abiotic stress response. With the improved genetic transformation system of wheat [113,114], the biological functions of the variations in these genes can be clarified in follow-up studies. In addition, the identified variations should be compared with the haplotypes identified from wheat breeding and germplasm resources using GBS, re-sequencing, and SNP microarray analyses to develop optimal molecular markers. With these markers, more wheat accessions and their derived lines can be tested and screened for stress-tolerance enhancement. The high-density genomic variations identified in our study could serve as a rich resource for further genetic haplotyping analyses, functional verification of molecular biology, and development of molecular markers for breeding selection.

5. Conclusions

We systemically identified 19,035,814 SNPs and 719,049 InDels in the genomes of two popular winter wheat cultivars, LM14 and HX10, using a re-sequencing assay with 7–8× average coverage. We carried out a time-course experiment of ABA treatment, and profiled and identified 61,251 expressed genes using a strand-specific RNA sequencing approach. A large number of genes showed time-point specific and/or cultivar-preferential responsive expression patterns. GO enrichment analysis revealed that the ABA-responsive genes were associated with stress-related functions and multiple biological processes, such as nicotianamine biosynthesis, nitrate and potassium transport, sucrose and trehalose biosynthesis, and cell growth. On the 20 large-effect QTLs, we uncovered 411 expressed genes carrying SNPs/InDels in their CDS regions, 306 genes with high- and/or moderate-effect variations, and 472 cultivar-preferentially expressed genes. Detailed analysis and verification of the homozygous genomic variations in the candidate genes encoding sulfotransferase, proteinase, kinase, nitrate transporter, and transcription factors suggested previously unexpected pathways associated with abiotic stress responses in wheat.