Yu Wei ,Lingjing Co ,Xucheng Hung ,Xun Wng ,Hun Wng ,Yongcho Song ,Qing He ,Mingjie Lyu ,Xinwen Hu ,Jun Liu ,*

a National Key Facility for Crop Gene Resources and Genetic Improvement,Institute of Crop Sciences,Chinese Academy of Agricultural Sciences,Beijing 100081,China

b Institute of Tropical Agriculture and Forestry,Hainan University,Haikou 570228,Hainan,China

c Department of Biological Sciences,Xi’an Jiaotong-Liverpool University,Suzhou 215123,Jiangsu,China

d Biotechnology Research Institute,Chinese Academy of Agricultural Sciences,Beijing 100081,China

ABSTRACT The diploid wild goat grass Aegilops tauschii (Ae.tauschii,2n=14;DD),as the D-sub genome of common wheat,provides rich germplasm resources for many aspects of wheat breeding.Abscisic acid (ABA) is an essential phytohormone that plays a pivotal role in plant adaptation to abiotic stresses.However,the gene regulation network of Ae.tauschii in response to ABA stress remains unclear.Here,we conducted a time-course strand-specific RNA-sequencing study to globally profile the transcriptome that responded to ABA treatment in Ae.tauschii.We identified 4818 differentially expressed transcription units/genes with time-point-specific induction/repression patterns.Using functional annotation,one-to-one ortholog and comparative transcriptome profiling analyses,we identified 319 ABA-responsive Ae.tauschii orthologs that were also induced/repressed under ABA treatment in hexaploid wheat.On the quantitative trait loci(QTL)used in wheat marker-assisted breeding,we found that the ABA-responsive expression patterns of eight Ae.tauschii orthologs were associated with drought stress tolerance,flowering process and/or grain quality.Of them,the ABA-responsive gene encoding sucrose:sucrose 1-fructosyltransferase in fructan and glucose metabolism pathways showed the most significant association with wheat drought tolerance.The characterization of ABA early-responsive genes in this study provides valuable information for exploring the molecular functions of the regulatory genes and will assist in wheat breeding.

Keywords:Wheat Aegilops tauschii Transcriptome Abscisic acid Fructan

1.Introduction

Approximately 8000 years ago,the diploid wild goat grassAegilops tauschii(Ae.tauschii,2n=14;DD) hybridized with the cultivated tetraploid wheatTriticum turgidum(2n=4x=28;AABB)and generated the allohexaploid wheat (Triticum aestivumL.,AABBDD) [1].These species carrying the three sub-genomes,AA,BB and DD,were suggested to originate from a common ancestorTriticeaeapproximately 2.5 to 4.5 million years ago [2].The most distinct feature of the allohexaploid wheat gained fromAe.tauschiiare its improved grain quality,which enables the allohexaploid wheat flour as a versatile material to produce various foods,such as bread,naan,and buns.Moreover,Ae.tauschiiis a species derived from arid and semi-arid areas that increased the abiotic stress tolerance of the allohexaploid wheat.Thus,the D sub-genome plays important roles in grain quality and abiotic resistance in wheat,even though it contains fewer genes than those of the A and B genomes[3,4].The recent high-resolution integrated genetic map and improved assembly of the referenceAe.tauschiigenome [5] will highly accelerate the understanding of wheat allopolyploidization.

Abscisic acid (ABA),as a plant hormone,can mediate gene expression by regulating multiple key developmental processes,such as hydraulic conductivity and seed germination [6,7].The ABA signalling pathway is involved in multiple stress responses,including drought,salinity,cold and/or heat [8].In the ABA signalling network,ABA molecules are sensed and bound to the ABA receptor complex[9],resulting in the inactivation of the enzymatic activity of 2C protein phosphatases (PP2Cs) [10].The inhibition of PP2Cs leads to the activation of sucrose-non-fermenting kinase 1-related protein kinase 2 (SnRK2) [11],which affects ABAresponsive transcription processes by interacting with the bZIP transcription factor(TF)ABI5 and B3 domain-containing transcription factor ABI3 [8],triggering downstream ABA responses in plants,such as promoting stomatal closure and repressing root elongation [12].ABI3,ABI4,and ABI5 are ABA-responsive binding factors(ABFs)that bind to promoter regions and induce or repress downstream genes [13].

ABA can increase plant tolerance to abiotic stresses by stimulating stomatal aperture,accumulating osmo-compatible solutes,and regulating plant development [14].ABA stimuli induce reactive oxygen species (ROS) production,which reduces the activities of ABI1,ABI2,and MAP kinases forming a feedback loop [15].An exogenous application of ABA in maize can enhance the activities of several antioxidative enzymes,such as superoxide dismutase and ascorbate peroxidase [16].In wheat,the application of exogenous ABA can enhance the dormancy or resistance to pre-harvest sprouting [17].ABA concentration has been found to be closely associated with the activity of fructan-and sucrose-metabolizing enzyme 1-SST [18].Because the transcriptional changes in many ABA-responsive genes occur at an early stage,approximately 30 min to 2 h with a peak at 1–2 h[19,20],investigatingAe.tauschiitranscriptome alterations under a time course of ABA stress will help elucidate how ABA affects gene expression during crop development.

In recent years,efforts have been made to detect the transcriptome changes in wheat under different environmental stresses,such as freezing [21],salt [22],and drought [23].Here,we conducted a time-course study and profiled the transcriptome of ABA-responsive genes inAe.tauschiiusing strand-specific RNA sequencing (ssRNA-seq).These genes were enriched in stress response-related functions.Ortholog and quantitative trait locus(QTL)analyses uncovered a group ofAe.tauschiiorthologous genes controlling grain quality,rust resistance,and drought stress tolerance,such asLOX,LR67,SR36,andSST,in hexaploid wheat.Our study provides transcriptome information to explore regulatory genes inAe.tauschiiand identifies a group of loci that may encode functional genes to assist wheat breeding.

2.Materials and methods

2.1.Plant materials for RNA-seq

Ae.tauschii(2n=14;DD) AL8/78 was used in this study.The seeds were germinated on wet filter paper in Petri dishes with three biological replicates and incubated at 20 °C in a growth chamber in the dark.The seedlings were grown in a growth chamber under a long day condition (24 °C/18 °C,16-h light/8-h dark photoperiod cycles) with 50% humidity.The 14-day-old seedlings were subjected to ABA treatment.The plants in 18 pots were sprayed with 200 μmol L-1ABA solution,and the plants in the other 18 pots were sprayed with the same volume of doubledistilled H2O (ddH2O) and served as mock treatment.The tissues were separately harvested from the shoots at five time points(0 h,1 h,6 h,12 h,and 24 h).To reduce the circadian effects,we harvested the tissues at 10:00 AM on each sampling day.The samples were frozen immediately in liquid nitrogen and stored at-80°C for further use.Two PP2C genes(PP2C30 andPP2C68)were used as marker genes to evaluate the response to ABA inAe.tauschii.

2.2.Total mRNA extraction and qRT–PCR analysis

Total RNA extraction of individual plants at each time point and under each treatment was conducted using an RNeasy Plant Mini Kit (Qiagen,Hilden,Germany) according to the manufacturer’s instructions.The RNA concentration was initially determined using a Nanodrop 2100 spectrophotometer.The integrity of the RNA samples was determined using an Agilent Bioanalyser with all samples;all the samples had an RIN score>8.5.

The total RNA was reverse transcribed using a PrimeScript RT Reagent Kit (Takara,Dalian,China) and gDNA Eraser (Perfect Real Time) (code No.RR047A)to synthesize the first-strand DNA.Realtime PCR was performed in optical 96-well plates using SYBR Premix Ex Taq 582 II (Takara,Dalian,China) and the ABI 7300 Real-Time PCR system(Applied Biosystems,Foster City,USA).Ae.tauschiiGAPDH(glyceraldehyde-3-phosphate dehydrogenase)was used as an endogenous reference gene (forward primer:5′-TTAGACTTGC GAAGCCAGCA-3′,reverse primer:5′-AAATGCCCTTGAGGTTTCCC-3′)for sample normalization[24].Three technical and three biological replicates were performed to produce the relative expression to the reference gene using the 2-ΔΔCTmethod[25].

2.3.RNA library construction and Illumina sequencing

RNA sequence libraries were prepared using NEBNext Ultra Directional RNA Library Prep kit for Illumina according to the manufacturer’s instructions.The quality and size of the cDNA libraries were checked using the Agilent 2200 TapeStation system (Agilent Technologies,Waldbronn,Germany).The RNA libraries were sequenced using the HiSeq X10 sequencing system (Illumina,CA,USA) with a 150-cycle paired-end sequencing protocol [26].

2.4.Analysis of RNA-seq datasets

TheAe.tauschiigenome sequence (Aet v.4.0) and annotation files were downloaded from the ENSEMBL Plants database [5].HISAT2-build was used to build index files for all the scaffold sequences.We used HISAT2 with ‘‘--rna-strandness RF”to align RNA-seq reads against the index files.Gene structures were assembled using Cufflinks and Cuffcompare with the ‘‘--library-type frfirststrand”parameter.Read counts were obtained using HTseqcount with the‘‘-s reverse”parameter[27].Transcripts per million(TPM) were calculated and normalized using HTseq-count and DESeq2 with global normalization parameters [28,29].

2.5.Differential gene expression analysis

Differential expression analysis was performed using DESeq2[24,29].The read count on each gene was aligned and calculated using HTseq-count based on the GTF-formatted file produced by Cufflinks and Cuffcompare.Subsequently,we applied DESeq2 to normalize the gene expression levels and perform differential expression analysis based on the negative binomial distribution[30].Genes with a normalized expression fold-change greater than 2,a significantP-value less than 0.01 and a Benjamini-Hochberg adjustedP-value/false discovery rate(FDR)less than 0.1 were considered to be differentially expressed genes.A custom Perl script was used to integrate and summarize the results produced by HTseq-count and DESeq2.

2.6.Gene ontology(GO)and kyoto encyclopedia of genes and genomes(KEGG) analyses

GO enrichment analysis was performed using agriGO,GOEAST and ggplot2 based on the hyper geometric test [31–33].KEGG pathway analysis was conducted using clusterProfiler [34].

2.7.Identification of orthologs in wheat

Reciprocal best hit blast was performed to identify orthologs ofAe.tauschiiin bread wheat (Triticum aestivumL.,Chinese Spring).High-confidence protein sequences ofAe.tauschiiwere downloaded from http://aegilops.wheat.ucdavis.edu/ATGSP/annotation/[5],and high-confidence protein sequences of wheat were downloaded from https://wheat-urgi.versailles.inra.fr/Seq-Repository/Assemblies [35].The criteria of max_target_seqs=1,evalue<10-5,5060 were employed to filter the target genes [36–38].

In addition to the wheat annotation,we also searched forAe.tauschiiorthologous genes inArabidopsisand rice.Matches to the transcription factor annotation ofArabidopsiswere downloaded from AGRIS and DATF databases [39,40].We used a custom Perl script to match the annotation information ofArabidopsisand/or rice genes to that of wheat genes.

2.8.Ae.tauschii transcriptome annotation

The assembled transcripts were annotated by Trinotate v3.1.1(https://trinotate.github.io/).The Transdecoder was employed to detect peptide coding regions.BLASTp and BLASTx v2.7.1 (release 2017-10,e-value cutoff of 10-5) were performed to determine the sequence homology to UniProt/SwissProt.HMMR v3.1b2(http://hmmer.org/) and Pfam database v31.0 (http://pfam.xfam.org) were used to identify conserved protein domains.The transmembrane regions were predicted using TMHMM 2.0 (http://www.cbs.dtu.dk/services/TMHMM-2.0/) [41],and potential signal peptides were generated using SignalP 4.1 (http://www.cbs.dtu.dk/services/SignalP).The RNAmmer 1.2 server was used to identify the rRNAs[42].Additionally,EggNOG v4.5.1(Evolutionary genealogy of genes:Non-supervised Orthologous Groups)was applied to predict the GO functional classification of unigenes based on known orthologous genes [43],and Enzyme Commission (EC)terms and metabolic pathways were generated by the KEGG database [44].Moreover,GO terms were inferred from BLAST hits and Pfam domains with blast2GO [45].

2.9.Verification of gene expression by qRT–PCR

A total of 1 to 2 μg of DNase I (RNeasy Plant Mini Kit)-treated RNA was reverse transcribed using SuperScript III (Invitrogen,Carlsbad,CA,USA) and oligo(dT) primer.cDNA was analysed by quantitative PCR using SYBR Premix Ex Taq(Takara,Dalian,China)and the Biorad CFX96 real-time PCR system.The normalized cDNA was amplified by 40 cycles at an annealing temperature of 60 °C using ABI StepOnePlus (Applied Biosystems,Foster City,CA,USA)and Roche LightCycler 480 II (Roche,Basel,Switzerland).We used three technical replicates to generate the average expression levels relative to the reference genes using the 2-ΔΔCTmethod [25].The primers are listed in Table S1 online.

2.10.Accession numbers

The experimental information of the RNA-seq datasets and fastq-formatted raw sequences are available in the Genome Sequence Archive database (http://bigd.big.ac.cn) under accession CRA001085.

3.Results

3.1.Genome-wide identification of expressed genes in Ae.tauschii

The genome sequences ofAe.tauschiiaccession AL8/78 was recently assembled [5].To explore the gene regulation network underlying ABA treatment,AL8/78 seedling shoots with leaves were isolated after 1 h,6 h,12 h,and 24 h of ABA and mock treatments,respectively (Fig.1A).Before ABA treatment,we observed that approximately 24% of the stomata were closed.After 3 min of treatment,approximately 97% of them were closed (Fig.1B),indicating a physiological impact of ABA treatment.Previous studies in model plants have uncovered that the early transcriptional response to ABA occurs approximately 30 min to 2 h with a peak at 1–2 h,including the first-round responses of the regulatory genes [19,20].To further confirm the time-point-specific response ofAe.tauschiiplants under ABA treatment,we used twoPP2Cgenes,AtPP2C30 andAtPP2C68,as marker genes (Fig.1C).Realtime quantitative reverse transcription–polymerase chain reaction(qRT-PCR)showed that the two genes were specifically induced at 1-h ABA treatment,whereas the genes under mock treatment were not highly induced.This result confirmed the existence of an early transcriptional response to ABA inAe.tauschii.Therefore,we used the 1-h treatment as the first time point and constructed ssRNAseq libraries of nine samples and sequenced using the Illumina platform.

After removing the low-quality reads and adapter contaminants,850,425,685 reads were aligned to theAe.tauschiireference genome(Aet v.4.0).Approximately 92%of the RNA-seq reads could be mapped to the reference genome (Table S2).Using a bioinformatics pipeline with strand-specific parameters (methods),we assembled the expressed genes and calculated the read counts uniquely mapped to the genome.In total,we identified 35,548 expressed transcription units (TUs)/genes (normalized uniquely mapped reads with counts >3 in at least one sample).Among them,32,931 genes were overlapped with annotated genes,while the remaining 617 genes were unannotated compared with the referenceAe.tauschiigene annotation.

3.2.Expression pattern of ABA-responsive genes

Using pair-wised differential expression analysis based on TPM values,we identified 4818 significantly differentially expressed genes (fold-change of TPMs between at least two samples >2;P-value of DESeq2 <0.05),including 4688 annotated genes and 130 genes encoding predicted protein-coding ORFs.Subsequently,we carried out differential expression analysis between ABA-treated and mock-treated samples.In total,2017 genes were significantly induced or reduced following ABA treatment (the log2fold-change between the treated and mock samples at the same time point was more than 1).Among them,1930 genes were identified to be overlapped with annotated genes,87 genes were unannotated and significantly induced or repressed under ABA treatment (Fig.2A).

We further carried out expression pattern analysis of the differentially expressed genes by clustering the normalized expression levels.The expression levels of 4818 genes were classified into 100 clusters,which could be categorized into five major expression patterns according to their expression specificity,including 1-h to 6-h short-term induction (Fig.2B;Group 100),1-h to 6-h shortterm reduction(Fig.2B;Group 62),12-h to 24-h long-term induction(Fig.2B;Group 68),12-h to 24-h long-term reduction(Fig.2B;Group 8),and 1-h to 24-h fluctuating expression (Fig.2B;Group 6 and Group 19).As shown in Fig.2B,the Group 68 showed stable gene expression following 1-h to 12-h ABA treatment but increased expression levels of genes under 24-h long-term ABA treatment;thus,we considered this group of genes to contain 12-h to 24-h long-term induction genes.The genes of Group 62 showed a reduced expression profile following 6-h ABA treatment,while the gene expression of Group 100 was promoted following 6-h ABA treatment by more than one fold;we referred these two groups as containing 6-h short-term reduction or induction genes.

Fig.1.The early ABA-response in Ae.tauschii.(A) Schematic model of the ABA and mock treatments.The 14-day-old Ae.tauschii seedlings were sprayed with 200 μmol L-1 ABA solution and the seedling sprayed with ddH2O were taken as mock treatment.Seedling shoots with leaves were isolated at three time points(0 h,1 h,and 6 h).(B)Guard cell responses to the ABA signal.The stomatal complex of Ae.tauschii is formed by two dumbbell-shaped guard cells flanked by two subsidiary cells.The opening stomatal pores under mock treatment and closed stomatal pores under ABA treatment are indicated by arrows.(C) Induced expression of the PP2C genes under ABA treatment.The expression levels of PP2C30 and PP2C68 was significantly upregulated under 1-h ABA treatment.The bars indicate the standard error of the values (n=3).

The gene expression specificity of each group illustrated their specific functional roles in plants following ABA treatment.At each time point of ABA treatment,more downregulated genes were observed than upregulated genes (Fig.3C),suggesting that a genome-wide transcriptional repression event plays a major role in the response to ABA stress.Additionally,the repression event mainly occurred at the middle stage of ABA treatment (6–12 h;Fig.3C),implying that many genes responded negatively to ABA treatment inAe.tauschii,which possibly due to ABA playing vital roles in dormancy in plant organs [46].This result was similar to that reported in sorghum and strawberry [47,48].

3.3.Identification of functional ABA-responsive genes

Our clustering analysis characterized a subset of genes showing time-point specific induction/repression expression patterns at different time points (Fig.3B).To investigate the biological functions of these genes,we searched for their orthologs in wheat,rice,andArabidopsis.The representative genes,including those encoding TFs,were listed in Table S3.

Furthermore,we investigated the expression trends of top expressed TFs at different time points.The expression levels of 58 TF genes belonged to the upregulated group under 1-h ABA treatment and those of 60 ones were downregulated at 6-h treatment (Fig.3D).WRKY,ERF,bHLH,NAC,andMYBwere the top five highly expressed TF families.A similar number of TFs in the same family across different time points suggests a relatively consistent expression pattern at different time points (Fig.3D).

In the early ABA-responsive gene groups at 1 h and/or 6 h,we identified genes encoding TFs in the ABA signalling pathway(Table S4),such asWRKY24,WRKY28,WRKY50,WRKY57,ERF1,andMYB15.Additionally,we identified a group of downstream genes belonging to the ABA signalling pathway,such as the genes encoding peroxidase 2,probable 2-oxoglutarate-dependent dioxygenase ANS,cationic peroxidase SPC4,and proline dehydrogenase 2.The genes encoding ABA receptors PYL5 and PYL8 and their interacting regulators PP2C9,PP2C29,PP2C30,PP2C32,PP2C34,and PP2C68 showed a similar time-point-specific induction pattern under 1-h ABA treatment.PP2Cs are direct regulators of SnRK2,which controls the ABA signal transduction pathway in plants [11].Their expression patterns profiled in this study indicated the early response following a 1-h ABA signal.

Several ATP-,NAD(P)H-,photosystem-and cytochromeassociated genes were upregulated under 1-h to 6-h ABA treatment but were downregulated under 12-h to 24-h ABA treatment(Table S4,column V).A group of genes encoding photosystem proteins,such as the ones encoding photosystem I assembly protein Ycf4,photosystem II protein D1 (PSBA),photosystem II reaction centre protein H (PSBH),and photosystem I P700 chlorophyllaapoprotein A1 (PSAA),were successively induced by ABA treatment from 1 h to 6 h,and then the expression was reduced under 12 h and/or 24 h ABA treatment(Fig.2B;Group 100).Furthermore,some cytochrome-encoding genes,such ascytochrome P450 andcytochrome f,presented the same induced expression pattern.Similarly,the expression of most ATP synthase-encoding genes,includingATP synthase subunit alpha(ATPA),ATP synthase subunit beta(ATPB),ATP synthase subunit a(ATPI),andATP synthase subunit b(ATPF),showed similar expression patterns as the genes encoding photosystem proteins.The proteins encoded by these genes can catalyze light-driven ATP synthesis,thus supplying the energy for photosynthesis[49].The expression levels ofNAD(P)H-quinone oxi-doreductase subunit I,NADPH-dependent aldo–keto reductase,andNAD(P)H-quinone oxidoreductase subunit1 (NU1C) were also strongly enhanced by ABA treatment from 1 h to 6 h,and these genes are functionally related to the energy-converting process[50].These genes were successively induced by 1-h to 6-h ABA treatments,reflecting a continuous transcriptome response to the ABA signalling pathway at an early stage.In addition to the genes responding to 1-h and/or 6-h ABA stress,some genes encoding stress-related proteins were induced by successive ABA treatments at 12 h and/or 24 h,such asWRKY71,NAC22,PER2,andDehydrin(DHN3).

Fig.2.Genome-wide identification of ABA-responsive genes in Ae.tauschii.(A) Schematic diagrams showing the expressed TUs/gene in mock-treated and ABA-treated samples.The number represents the TUs/gene in different samples.The numbers of expressed TUs,differentially expressed TUs between two samples,and differentially expressed TUs between mock-treated and ABA-treated samples at the same time point are given from left to right,respectively.(B)Clustering analysis of DEGs based on the expression patterns.Each cluster represents a set of differentially expressed genes.The light grey lines in each cluster indicate the ratios of gene expression levels under mock and ABA treatments relative to those of the 0-h samples,and bold pink lines represent the average expression ratios of the genes at each sampling point.The y-axis represents the log2 Fold-change.TU,transcript unit;DETU,differentially expressed transcript unit.

To further confirm our RNA-seq analysis and validate the expression patterns of ABA-responsive genes at different time points,eight differently expressed genes,PER2 (AET2Gv20218200),WAK3 (AET5Gv20612600),ZFP1 (AET4Gv20165000),WRKY76(AET5Gv20544800),DHN3 (AET6Gv20865700),TIFY10c(AET5Gv20519300),HOX24 (AET6Gv20626600) andMYB20(AET7Gv20519300),were randomly selected to verify their expression levels using qRT–PCR assay(Fig.4).The expression of most of these genes was consistent with that detected by ssRNA-seq,indicating a high level of reliability of our ssRNA-seq analysis.

Next,we compared our transcriptome profiling results with those of ABA-responsive genes identified in common wheat after 24-h ABA treatment in a recent study [51].Among 319 ABAresponsive wheat orthologs,we identified 154 with consistent expression profiles under 24-h ABA treatment,including 43 upregulated and 111 downregulated genes (Table S5).One hundred forty-four orthologous genes with functional similarities were located on the wheat D genome.The allohexaploid hybridization events ofTriticum turgidumwithAe.tauschiioccurred approximately 8000 years ago.The sequence similarity between theAe.tauschiigenome and D-subgenome is much higher than that between the A-and B-subgenomes.Thus,it is not surprising that most of theAe.tauschiiorthologous genes are located on the D-subgenome.Moreover,the expression specificities of these orthologs in two studies suggest their functional associations,highlighting the functional importance ofAe.tauschiigenes in wheat breeding and improvement.

To elucidate the functions of the induced/repressed differentially expressed genes,we conducted GO classification analysis.We identified 167 significantly enriched GO terms comprising 2350 DEGs.The top 20 enriched GO terms are shown in Fig.5A.These genes were significantly enriched in molecular functions and biological process (false discovery rate <0.05),such as ‘‘catalytic activity”,‘‘response to stimulus”,‘‘response to hormone”,and‘‘response to abscisic acid”,suggesting that the genes classified in these terms may be involved in the ABA signalling pathways inAe.tauschii.Furthermore,many genes were designated to cellular component terms,such as‘‘intrinsic component of membrane”and‘‘integral component of membrane”,indicating that the membrane and cell junction could function in signal transduction.

Fig.3.The expression patterns of ABA-responsive genes in Ae.tauschii.(A)The Venn diagram illustrates the number of annotated genes in the top five databases.(B)Heatmap depicts the fold-change expression of upregulated and downregulated genes involved in the ABA response.(C)Number of upregulated genes(blue)and downregulated DEGs(red) following ABA treatment.(D) The bar plot demonstrates the number of identified TFs following ABA treatment.WRKY,transcription factors in which N-terminal ends contain the WRKYGQR amino acid sequence and a zinc-finger motif;ERF,ethylene responsive element binding factor;bHLH,basic helix-loop-helix transcription factors;NAC,no apical meristem (NAM), Arabidopsis thaliana transcription activation factor (ATAF1/2),and cup-shaped cotyledon (CUC2);MYB,myeloblastosis viral oncogene homolog.(E) Upset plot of the overlapped DEGs in different groups.up,upregulated;down,downregulated.

Fig.4.Relative expression levels of eight randomly selected genes as detected by ssRNA-seq and qRT–PCR,respectively.The bars indicate standard error (n=3).Dark blue and red,RNA-seq analysis;Light blue and yellow,qRT–PCR analysis;TPM,transcripts per million.

Fig.5.GO(A)and KEGG(B)pathway enrichment of differentially expressed genes between ABA-and mock-treated samples in Ae.tauschii.

To better elucidate the biological functions of the DEGs,we conducted KEGG analysis.In total,505 DEGs were significantly enriched in 32 pathways,and the top 20 most highly enriched pathways are shown in Fig.5B.Consistent with the observation of GO analysis,the most abundant KEGG pathways were located in ‘‘circadian rhythm-plant”,‘‘peroxisome”,and ‘‘plant hormone signal transduction”.These results imply that most of the GO and KEGG assignments of the identified genes were stress and hormone signal-responsive genes associated withAe.tauschii.

3.4.Ortholog analysis of Ae.tauschii genes

To systemically explore the biological functions of the assembled genes inAe.tauschii,we aligned the peptide sequences encoded by the predicted ORFs to identify homologous similarities with known protein sequences,domains,motifs,and/or secondary structures (Fig.3A).We used BLASTX and BLASTP to align the sequences against the NCBI Non-Redundant protein database and identified 32,660 and 27,744 significantly matched sequences,respectively[52].We also aligned the sequences to the Pfam database with HEMMER/PFAM and found 29,529 genes containing significant matches.Moreover,we predicted the signal peptides encoded by 6450 genes using SignalP as well as transmembrane regions of proteins encoded by 11,417 genes using TmHMM.Finally,we annotated the GO terms of the assembledAe.tauschiigenes by aligning their protein sequences to the NCBI Non-Redundant protein database and Pfam database using blast2GO,and significantly matched GO terms were found on 30,933 and 19,058 genes,respectively.In total,329,910 annotations were identified for theAe.tauschiiassembled genes.These results and gene annotation information are presented in Table S6,respectively.

Next,we compared the protein sequences ofAe.tauschiiwith those of the hexaploid wheat cv.Chinese Spring to identify orthologs[53].In total,we detected 21,842Ae.tauschiigene orthologs in the wheat D-subgenome,including 10,688 that displayed one-toone 100%identical matches and were referred to as representative orthologs(Table S7).Among them,10,333Ae.tauschiigenes owned orthologs encoded by the wheat D-subgenome.The remaining 355 orthologs were identified in wheat A-and B-subgenomes,or were unknown,suggesting that these genes were diversified during evolution,a finding consistent with that in the previous study [5].

3.5.ABA-responsive genes associated with important wheat QTL

We explored the potential functions of the ABA-responsive genes in wheat breeding by comparing these genes with 13 widely used wheat genes/QTL markers associated with drought tolerance and/or pre-harvesting sprouting [54,55].We found that oneAe.tauschiiortholog(AET7Gv20017900,SST)of two genes in wheat(sucrose:sucrose1-fructosyltransferaselocus,TaSST-4DandTaSST-4A),a key enzyme in fructan biosynthesis,responded to ABA treatment.SST can regulate the carbon flow of photosynthesis products between different tissues by transferring fructose from one sucrose to another and generating glucose and 1-kestose [56].SST is associated with drought,cold and osmotic stresses and plays an important role in counteracting oxidative stress [57].In this study,we found that the expression level of theTaSSTsortholog was increased at 1 h and/or 6 h under mock treatment;whereas this rapid induction was arrested under ABA condition at 6–24 h(Fig.6).

We also compared the ABA-responsive orthologs with 77 wheat genes/QTL associated with photoperiod sensitivity,disease,flowering,and grain quality [54,55,58–65].We found the orthologs of seven gene loci,Ppd-D1,Lox-B1,PRR-B1,Lr67,Pm21,Lr34,andSr36,in 77 agronomically important genes/QTL (Fig.6).According to the QTL study,TaLox-B1 co-segregated with two lipoxygenase(LOX) gene markers and correlated with the simple sequence repeat (SSR)locus on wheat chromosome 4B [66].TaLr67 was significantly associated with stripe rust resistance-related markers on chromosome 4D [65,67,68].

4.Discussion

4.1.Transcriptome profiling of Ae.tauschii and differential gene expression under ABA treatment

Ae.tauschii,a wild relative of wheat,is the D-genome progenitor of common wheat [1,69].By transcriptome analysis,we identified 33,548 expressed genes,including 617 unannotated genes thus far.Their uncharacterized coding/noncoding biological functions remain to be uncovered in the future.These results demonstrate that most of the genes have been characterized inAe.tauschii,and this good annotation will facilitate a deeper understanding of A/B/D homologous gene function and evolution in wheat.

4.2.DEGs involved in ABA signalling pathway

Our analysis identified several major classes of ABA response regulators following 1-h and/or 6-h ABA treatment,including genes encoding the ABA receptor (PYL8),ABA regulators (PP2Cs),and ROS-related genes,as well as two light signalling related genes(PIF3andPIF1-like helicase).These genes showed a time-pointspecific induction pattern at 1 h,demonstrating their early response following ABA signalling.

Fig.6.ABA-responsive genes associated with important wheat QTL detected by RNA-seq.TPM,transcripts per million.Bars indicate standard error(n=3).

The orthologous genes ofNAC,MYB,andbZIPfamilies have been previously characterized to participate in the stress response in plants such asArabidopsis,rice,and wheat[70–73].WRKYsare also essential components in the ABA signalling pathway through mediating the transcription process of their targeted genes,including well-known ABA-responsive genes such asABF2,ABF4,ABI4,andABI5 [74].ThebZIPsubfamily genes could be enhanced by exogenous ABA and other stress treatments in vegetative tissues[75].Several studies have shown that ABA induces leaf senescence and NACs play crucial roles in ABA-induced leaf senescence signalling [76].In summary,the identification of these ABAresponsive TF genes revealed that the transcriptional control of stress-responsive gene expression is a pivotal part of plants in response to various abiotic stresses.

In addition,we observed that ethylene,auxin,and ABA signalling coordinately modulate specific genes associated with plant senescence,such asRAP24,EAR,andDRMH3.Moreover,some of the typical stress-inducible genes,such asNACs,WRKYs,MYBs,HSFs,andEIN2 detected in our datasets,were responsive to other hormones [77–81].Together,these lines of evidence demonstrate that different phytohormone signalling pathways possibly share common regulators in plants,as well as a complex cascade of gene expression in response to osmotic stresses,especially the crosstalk in ABA-dependent and ABA-independent signalling pathways,remains to be unveiled.

The identification of photosystem-,cytochrome-and energyassociated genes at 1-h or/and 6-h ABA treatment implies that the ABA signalling may regulate the photosystem at an early stage by influencing the carbon fixation-and cytochrome-related genes.In addition to the genes induced by 1-h and/or 6-h ABA treatment,a group of genes showed induction/repression expression dynamics following 12-h and/or 24-h ABA treatment,comprisingPER2,DHN3and a set of transcription factor-encoding genes,includingWRKY51,WRKY71,ERF2,NAC22,andNAC29,demonstrating the downstream transcriptional response of these genes in the ABA pathway.Generally,during plant adaptation to osmotic stress,various signalling proteins,such as TFs,protein kinases,and proteins that control photosynthesis and that perform functions ranging from stress signal perception to signal transduction,were strongly enhanced/reduced [82].

4.3.ABA-responsive genes on the important QTL

SST is an important fructan-synthesizing enzyme,and its decreased transcript levels with 6-h to 24-h ABA treatment suggests the accompanied inhibition of fructan synthesis inAe.tauschiileaves.This finding is consistent with supplemental ABA significantly repressing SST activity in the wheat stem [18].Because the SST found in our study was identified as an ortholog of two published wheat QTL markers,TaSST-4DandTaSST-4A,its transcript changes in response to ABA inAe.tauschiiimply its potential function in resistance to ABA stress in wheat.Consequently,our results raise the hypothesis that fructan metabolism and transport pathways play crucial roles in wheat stress-resistance breeding,which remain to be clarified in the future.Additionally,the orthologs of seven agronomically important genes/QTL,including those ofTaLox-B1,TaLr67,TaPm21,TaPpd-D1,TaPRR-B1,TaLr34,andTaSr36,were identified in this study.The finding highlights that the identification of ABA-responsive genes could contribute to molecular breeding and disease resistance programmes in wheat.A recent study demonstrated that all 107 examined accessions ofAe.tauschiicontained the same HS2-type allele ofTaSG-D1as that in common wheat,where the geneTaSG-D1that regulates the formation of round grains of wheat is encoded by this allele.This evidence indicates that understanding the gene reservoir ofAe.tauschiicould assist the genetic engineering of grain shape and optimize the flour extraction rate in modern wheat cultivars [83].

CRediT authorship contribution statement

Xucheng Huang,Jun Liu,Yu Wei,and Liangjing Cao:designed the study.Liangjing Cao and Yongchao Song:sampled the biological materials,extracted the RNA,and sequenced the libraries.Xucheng Huang,Xuan Wang,Huan Wang,Qiang He,Mingjie Lyu,and Yu Wei:analysed the datasets.Xucheng Huang,Yu Wei,and Jun Liu:prepared the figures and wrote the manuscript.All the authors reviewed and approved the final version of the manuscript.

Declaration of Competing Interest

The authors declare no conflict of interest.

Acknowledgments

This research was supported by the National Key Research and Development Program of China (2016YFD0101001) and the Agricultural Science and Technology Innovation Program of CAAS.We thank Long Mao,Aili Li,and Ruilian Jing at CAAS for providing technical assistance.

Appendix A.Supplementary data

Supplementary data for this article can be found online at https://doi.org/10.1016/j.cj.2020.10.008.