Differential Expression Analysis of Phytohormone-Related Genes of Korean Wheat (Triticum aestivum) in Response to Preharvest Sprouting and Abscisic Acid (ABA)

Preharvest sprouting (PHS) is a key global issue in production and end-use quality of cereals, particularly in regions where the rainfall season overlaps the harvest. To investigate transcriptomic changes in genes affected by PHS-induction and ABA-treatment, RNA-seq analysis was performed in two wheat cultivars that differ in PHS tolerance. A total of 123 unigenes related to hormone metabolism and signaling for abscisic acid (ABA), gibberellic acid (GA), indole-3-acetic acid (IAA), and cytokinin were identified and 1862 of differentially expressed genes were identified and divided into 8 groups by transcriptomic analysis. DEG analysis showed the majority of genes were categorized in sugar related processes, which interact with ABA signaling in PHS tolerant cultivar under PHS-induction. Thus, genes related to ABA are key regulators of dormancy and germination. Our results give insight into global changes in expression of plant hormone related genes in response to PHS.


Introduction
Preharvest sprouting (PHS) is a key agronomic issue because it can reduce grain yield and crop quality, especially of milling and baking quality wheat [1,2], causing economic losses [3]. In Korea, the highest rainfall is in summer, with continuous rain during the harvest period, which is a major concern for wheat cultivation. Thus, controlling dormancy and germination in cereal grains and improving PHS tolerance has become a major aim of Korean wheat breeding.
Seed dormancy is a trait that blocks germination of seeds under apparently favorable conditions. Many cereal crops, including wheat, have been selected for rapid and uniform germination, and therefore show weak dormancy, making them susceptible to PHS [4,5]. In wheat, quantitative trait loci (QTLs) for dormancy and PHS tolerance have been studied by many researchers. Ogbonnaya et al. [6] utilized recombinant inbred lines (RILs) progeny populations from a cross between PHS-tolerant and -susceptible parents to study the genetic basis of dormancy and PHS tolerance, and identify the chromosomal position of QTLs for PHS tolerance. The relationships between seed color and PHS tolerance QTLs located on group 3, 4, 5, and 6 chromosomes have been studied in hexaploid wheat in Canada in breeding of white wheat [7].
Hormonal regulation is a primary mechanism of seed dormancy and germination in seed plants. Previous studies indicated that germination is regulated by phytohormones including abscisic acid (ABA) and gibberellic acid (GA). ABA and GA are key regulators of dormancy maintenance and release, respectively, with antagonistic functions [8][9][10].
Liu et al. [11] investigated the expressions of ABA, GA, jasmonic acid (JA), and IAA (indole-3-acetic acid) related genes in after-ripened wheat seeds using microarray. Seed ABA levels are regulated by balanced biosynthetic and catabolic pathways [12,13]. Genes encoding key enzymes in ABA metabolism and catabolism have been intensively studied in dormant and dormancy-released seeds [14,15]. GA breaks dormancy and promotes germination. Typically, GA levels increase and ABA levels decrease in imbibed non-dormant seeds. Light and GA promote the degradation of ABA and the expression of the GA hormone metabolism gene increase [16]. Other phytohormones, such as auxin and cytokinin, e also associated with seed dormancy and germination. Auxin, which plays an important role in regulating plant growth and development [17,18], regulates seed dormancy via an interaction with ABA signaling and ABI3 activation [19]. Crosstalk with cytokinin, a phytohormone that affects various plant developmental processes, has been thoroughly reviewed by El-Showk et al. [20].
To identify the putative key genes involved in responses to PHS and ABA treatment, RNA-seq analyses were performed. In this study, differentially expressed genes (DEGs) related to ABA and GA hormone metabolism and signaling pathways in two Korean wheat cultivars with different responses to PHS (tolerant/susceptible) were identified. Furthermore, DEGs related to auxin and cytokinin hormone metabolism and signaling were analyzed for their interactions with ABA on seed dormancy regulation. Here, we suggest the putative key genes related in phytohormone metabolism and signaling pathways for regulation of seed dormancy using comparative transcriptome analysis.

Plant Materials and Treatments
The research materials used were 'Keumgang' (Korea RDA accession no.IT213100), which is a vulnerable cultivar to PHS, and 'Woori' (Korea RDA accession no.IT175538) known as PHS tolerance cultivar. Seeds of both varieties were subjected to vernalization for 4 weeks at 4 • C, and then transplanted into plastic pots with a diameter of 11 cm to grow until 35 days after fertilization (DAF 35). The DAF 35 spikes were harvested as a control and stored −70 • C until RNA extraction. The sandbury method was conducted to confirm the dormant reaction of the seed by performing PHS treatments under artificial conditions, and was conducted based on the Bayer [21] experimental method. DAF 35 spikes were designated as an experimental group, and one-third of spikes of each cultivar were exposed to the ground in a plastic box of 590 mm × 380 mm containing wet sand. PHS of spikes was induced by mist treatment for 12 h for a total of 7 days.
Abscisic acid (ABA) treatment was injected with 0.1 mM ABA every 2 days to the harvested DAF 35 spikes for 7 days. After 7 days treatment of PHS or ABA, the spikes were harvested as a control and stored −70 • C until RNA extraction 'Keumgang' (KC),'Woori' (WC) spikes harvested on DAF 35 were used as controls. 'Keumgang' (KP),'Woori' (WP) spikes treated with PHS and 'Keumgang'(KA), 'Woori' (WA) spikes treated with ABA were considered as the experimental group.

RNA Isolation and Illumina Sequencing
Total RNA extraction was performed on harvested spikes using the TRIzol method (Invitrogen, MA, USA) according to the manufacturer's instructions. Libraries for RNA-seq analysis were prepared using the TruSeq TM RNA sample preparation kit and 100-bp pairedend sequencing was performed on the Illumina HiSeq TM -2500 platform (Illumina, CA, USA). The DynamicTrim and LengthSort provided in SolexaQA [22]. All clean readings were assembled using Oases [23] and Velvet assembler [24], and the contacts created by de novo assembly of readings are called unigenes.

Bioinformatic Analysis of Rna-Seq Data
For gene annotation, de novo assembled unigenes were compared against databases including the Phytozome, Gene Ontology (GO), UniProt, eukaryotic Ortholog Groups (LOG), Kyoto Encyclopedia of Genes and Genomes (KEGG), and The Arabidopsis Information Resource (TAIR). Unigene sequences were also mapped to IWGSC MIPS v2.2 wheat transcriptome sequences (IWGSC 2014) using BLAST with an e-value cutoff of 1 × 10 −10 ; gene IDs with >95% identity and <2 bp mismatch were selected [25].
Paired-end clean reads were mapped back to de novo assembled unigenes to obtain read counts for expression analysis using Bowtie 2. Gene expression levels from read counts were normalized to Reads Per Kilobase per Million (RPKM) with the edgeR package of Bioconductor in the R software environment. Unigenes with > 500 raw read counts from at least one library were selected and the log2 value of the relative RPKM ratio was applied to determine relative expression levels. Genes with an absolute log2 value of the relative RPKM ratio > |1| (2-fold change) were defined as differentially expressed genes (DEGs). Heatmaps were generated by using the heatmap. 2 function in the gplots package in the R software environment. To analyze gene ontology (GO) enrichment and network analysis of DEGs, the lists of Arabidopsis thaliana orthologous were uploaded to ClueGO [26] plug-in of Cytoscape [27] and the hypergeometric test was applied (p < 0.05).

Validation of RNA-Seq by qRT-PCR
RNA-seq results were validated by qRT-PCR. Total RNA extraction was performed using RNA purification (Ribospin TM Seed/Fruit, GeneAll ® Biotechnology, Seoul, Korea) and extracted according to the protocol suggested by the manufacturer. After adjusting the total RNA concentration to be used for cDNA synthesis to 1000 ng/µL, first strand cDNA synthesis was performed using the Power cDNA Synthesis Kit (iNtRON Biotechnology, Seoul, Korea) as suggested by the manufacturer's instructions. Quantitative Real-time PCR was conducted using QIAGEN's Rotor-Gene Q (QIAGEN, Hilden, Germany) Real Time-PCR machine with Rotor-Gene SYBR Green PCR kit (QIAGEN, Hilden, Germany). Primer set for qRT-PCR were designed using Primer3 program and listed in Table S13 [28]. The relative expression level was normalized by DAF 35 Keumgang spike using the 2 −∆∆Ct [29] method.

RNA-Seq Data Analysis
A total of six libraries (control, PHS-induced and ABA treated samples of both cultivars, KC, KP, KA, WC, WP, and WA) were generated by paired-end sequencing on the Illumina HiSeq TM -2500 platform. After quality control and sequence pre-processing, clean reads with an average length of 82.22 bp were obtained and mapped to IWGSC MIPS v2.2 wheat transcriptome sequences. However, the average mapping rate for all samples was 45% and especially 10% for the WC library. The mapping results were not suitable for expression analysis and thus de novo assembly was performed. The Oases/Velvet assemblers generated 104,542 and 98,607 contigs for 'Keumgang' and 'Woori,' respectively; the 'Keumgang' assembled contigs, with an average length of 733 bp, were used as mapping reference sequences and regarded as unigene. A summary of the RNA-seq data is shown in Table 1. Among 104,542, a total of 16,727 unigenes according to our selection criteria (>500 raw read counts from at least one library) were selected for analysis.

Gene Ontology (GO) Network Analysis of DEGs
To identify the genes that interacted with phytohormones metabolisms or signaling, the DEGs were divided into 8 groups based on their comparative expression levels. DEGs included in phytohormone metabolisms and signaling pathways were excluded for the GO analysis. The comparative expression levels and number of DEGs are illustrated in Table 2. Group-I to group-IV represents the DEGs that showed higher expression in 'Woori' compared to 'Keumgang', and group-V to group-VIII represents the DEGs that showed higher expression in 'Keumgang' compared to 'Woori'. Group-I, II, V, and VI represents the DEGs that showed higher expression in ABA treatment, including the genes with differential expression between PHS treatment (group-I and V), or including the genes that are not the DEGs (group-II and VI). As we are interested in the DEGs expressed in the PHS tolerant cultivar, 'Woori', the DEGs that were known to interact with phytohormones were screened in the group-I to IV. GO analysis of group-V to VIII are represented in Figure S1 and Tables S5-S8. To understand the biological relevance of the DEGs of each analysis groups, we performed functional enrichment analysis using ClueGO with ontology sources of cellular component (CC), biological process (BP), and molecular function (MF). Among 456 DEGs included in group-I, 230 genes were matched to Arabidopsis orthologous and were mainly involved in the flavonoid biosynthetic process, benzene-containing compound metabolic process, and alpha-amino acid catabolic process ( Figure 1A, Table S1). Group-II have DEGs of four categories of L-phenylalanine metabolic processes, sugar:proton symporter activity, drug transmembrane transport, and positive regulation of translation ( Figure 1B, Table S2). In group-III, 85 Arabidopsis orthologous were identified and functionally categorized. A total of 30 GO terms were categorized into 12 functional groups including filiform apparatus and trehalose biosynthetic processes ( Figure 1C, Table S3). Group-IV has 7 categories including the trehalose biosynthetic process, potassium/proton antiporter activity, and GMP synthase (glutamine-hydrolyzing) activity ( Figure 1D, Table S4). Group-V has 10 categories including cysteine-type peptidase activity, cytoplasmic side of plasma membrane, and DNA conformation change ( Figure S1A, Table S5). Group-VI DEGs, which showed comparatively higher expression in 'Keumgang' under ABA-treatment but showed no significantly different expression between two cultivars under PHS-induction, were divided into 5 categories (nucleosome, regulation of anthocyanin metabolic process, cellulose biosynthetic process, GA mediated signaling pathway, and regulation of ABA-activated signaling pathway) ( Figure S1B, Table S6). Among 428 DEGs included in group-VII, 308 genes were matched to Arabidopsis orthologous (TAIR gene IDs) and were mainly involved in chloroplast thylakoid, including photosynthetic processes ( Figure S1C, Table S7). In group-VIII, there were 32 Arabidopsis orthologous genes among 43 DEGs divided into 7 functional categories including abscisic aldehyde oxidase activity, linoleate 9S-lipoxygenase activity, and the chitin catabolic process ( Figure S1D, Table S8).

Expression Analysis of Genes Related to ABA Hormone Metabolism and Signaling
A total of 8 unigenes were identified as encoding ABA metabolic enzymes, one for zeaxanthin epoxidase (ZEP), one for 9-cis-epoxycarotenoid dioxygenase (NCED), 4 for ABA-aldehyde oxidase (AAO), and two for ABA 8 -hydroxylases (CYP707A) (Figure 2A). The genes encoding ZEP were up-regulated differentially (RPKM ratio > 1) in PHS tolerant cultivar 'Woori' under both PHS and ABA treatment (WP and WA, respectively) and not differentially (RPKM ratio < 1) in PHS susceptible cultivar 'Keumgang' under both PHS and ABA treatment (KP and KA, respectively). Along with ZEP, NCED is involved in the early phase of ABA hormone metabolism (conversion of violaxanthin to xanthoxin). The gene encoding NCED was up-regulated only in 'Woori (WP and WA)'. Most AAO genes involved in ABA metabolic pathway were differentially up-regulated in all samples, and especially highly up-regulated in WA. Two ABA catabolic genes CYP707A were identified in our dataset and showed contrasting expression patterns. Based on our RNA-seq results, 28 unigenes (7 PYR/PYL/RCAR, 4 PP2C, 7 SnRK, 6 ABIs, and 4 ABF/AREBs) were annotated as ABA signaling related genes ( Figure 2B). Among 7 genes encoding ABA receptors PYR/PYL/RCAR, 5 genes showed same expression patterns (differentially up-regulated in KP, WP, and WA, and down-regulated in KA). Most signaling pathways downstream of PYR/PYL/RCAR were differentially up-regulated, but ABI5s and one AREB/ABF were differentially down-regulated in KP. The detailed information and comparative expression value of ABA metabolic and signaling genes are shown in Table S9.

GA Hormone Metabolism and Signaling Genes
Based on our RNA-seq dataset, a total of 12 unigenes were identified as encoding GA metabolic pathways. Among the GA metabolic genes, we identified 4 copalyl diphosphate synthase (CPS), 5 ent-kaurene synthase (KS), 2 ent-kaurene oxidase (KO), and one gibberellin 2-beta dioxygenase (GA2ox) ( Figure 3A). Starting with geranylgeranyl pyrophosphate (GGPP), CPS converts GGPP to ent-Copalyl-PP. One CPS (Traes_7BL_02F3E2C78.1) was differentially down-regulated in KP and up-regulated in WP, whereas other 3 CPS (Traes_2AS_A4BAAEE05.1, Traes_2DS_7917F0E44.1, and Traes_2BL_80AC86A1A.1) were highly up-regulated under KP and WP. The former is located on chromosome 7 (7BL), and the latter are located on chromosome 2 (2AS, 2DS, and 2BL). The five KS unigenes proved to be located on chromosome group 2. Four of these genes were highly up-regulated in KP and WP (Traes_2AS_2610FD515.1, Traes_2AS_D3212BACF.1, Traes_2BL_212C1CE5C.1, and Traes_2BL_B06E350C2.1) and identified to be located on 2AS and 2BL. The other KS was located on 2DL (Traes_2DL_BEDDAC0D4.1). GA2ox, which catalyzes biologically active gibberellin G 1 to in-active GA 8 , was differentially up-regulated in KP, WP, and WA, and not differentially down-regulated in KA. In the GA signaling pathway, 3 unigenes were identified as encoding GA receptors, GID1, 7 unigenes as encoding SCF complex (2 CUL1, 4 SKP1, and 1 RBX1), and one DELLA protein (RGA1) ( Figure 3B). Although most of the GA receptors and SCF complex genes were up-regulated, gene encoding RGA1, which repress the GA response, is also up-regulated in all treated samples and showed higher expression in 'Woori'. The detailed information and comparative expression value of GA metabolic and signaling genes are shown in Table S10.

Validation of DEG Data by qRT-PCR
To confirm the expression patterns, qRT-PCR analysis was performed using randomly selected unigenes from hormone metabolism, signaling related DEGs. The selected DEGs for qRT-PCR were listed in Table S13. As shown in Figure 6, the RNA-seq results fold change values were agreed with the qRT-PCR expression patterns. The comparison of RNA-seq data and qRT-PCR data suggested accuracy of sequencing methods.

Discussion
In this study, transcriptome analysis was performed to identify DEGs for PHS-induced and ABA-treated two Korean wheat cultivars with different PHS tolerance. The analysis was focused on the expression patterns of the genes related to hormones metabolism and signaling and differentially expressed genes between treatment and cultivars.
In ABA signaling pathway, ABI and AREB/ABF lead the ABA responses and are activated by SnRK2. Soon et al. [32] reported that activities of SnRK2.6 are controlled by ABA and four SnRK2.6 (Traes_2AS_64E0EF151.1, Traes_2DL_6CD064E13.1, Traes_5AL_4EE2511FA.1, and Traes_4AL_44E99BA01.1) were identified in our RNA-seq results (Table S9). The degree of up-regulation of SnRK2 and downstream genes ABIs and AREB/ABFs were higher in 'Woori' than 'Keumgang' in both PHS-induced and ABA-treated samples. Figure 3 showed the GA hormone metabolism and the expression patterns of related genes in our RNA-seq data. Huang et al. [33] reported the identification and chromosomal location of genes related to GA hormone metabolism in wheat using 'Chinese Spring' nulli-tetrasomic lines. They isolated three TaCPS genes localized on chromosomes 7A, 7B, and 7D and three TaKS genes localized on chromosomes 2A, 2B, and 2D. DNA gel blotting using HvCPSL1 probe showed strong hybridization bands on chromosome group 7 and weak hybridization bands on chromosome group 2 of wheat [34]. Gibberellins have several forms of bio-active and bio-inactive (or precursors) such as GA 1 and GA 8 , respectively. In our raw data there were GA20ox and GA3ox that convert bio-inactive GA 19 to GA 20 and bio-active GA 1 from bio-inactive GA 20 , respectively, but these were not selected in analysis due to the gene selection criteria, unfortunately. GA20ox was down-regulated in 'Keumgang' (KP and KA) and up-regulated in WA and the contents of GA 20 were lower in 'Keumgang' than 'Woori' (Figure 3). GA3ox showed no different expression between cultivars except PHS treated samples (higher in KP with RPKM ratio 2.362).
After the first identification of GID1 in GA signaling mutants of rice [35], GID1 orthologous were identified in Arabidopsis [36,37]. The ubiquitin-proteasome system (UPS), via SCF complex, is utilized in signaling via many plant hormones, including IAA, jasmonic acid (JA), GA, strigolactone, ethylene, and ABA [38], and SCF complex components were recently studied in developing wheat spikes [33,39]. Although all the genes encoding GA receptors and SCF complexes were up-regulated, DELLA gene was also up-regulated in 'Woori' (RPKM ratio 4.35 in WP and 4.72 in WA).
The other hormones, IAA and cytokinin, are known to regulate seed dormancy and germination via interactions with ABA. IAA interacts with ABA and ABI3. With high IAA content, ARF10/16 is activated and ABI3 is expressed to maintain seed dormancy [19,40], and PHS-tolerant cultivar 'Woori' showed up-regulation in all genes encoding ABI3 and ARF10/16, whereas 60% of the genes were up-regulated in 'Keumgang'. Wang et al. [41] proposed that cytokinin regulates germination via interactions between cytokinin signaling and ABA, so that the Type A ARR represses ABI5, promoting germination. Although expression of Type A ARR showed higher up-regulation in 'Woori' than in 'Keumgang,' expression of all three ABI5 was up-regulated in WP and down-regulated in KP. All the transcriptome analysis results of phytohormone hormone metabolism and signaling suggest that GA, IAA, and cytokinin may play a subordinate role to ABA in terms of dormancy and germination and ABA is the primary regulator of dormancy and inhibition of germination.
The DEG analysis also supports the role of ABA in seed dormancy and germination inhibition. DEG analysis group-II, III, and IV, which showed significant higher expression of DEGs in PHS-tolerant cultivar 'Woori' than 'Keumgang', had GO categories related to sugar biosynthetic processing or transport (trehalose biosynthetic process and sugar/proton symporter activity). In particular, 'Woori' showed higher expression of genes related in the trehalose biosynthetic process when PHS induced (group-III and IV). Trehalose is a disaccharide synthesized from two glucose molecules (UDP-Glc + Glc-6-P) [42,43] and widely used economically in areas such as food industry and medicine [44,45]. Dekkers et al. [46] reported, via the effect of sugars in 10 different Arabidopsis accessions, that glucose induced the delay of germination. They also investigated the germination delay effect other sugars (galactose, sucrose, maltose and trehalose) and the interaction with phytohormones. GA or ethylene was independent from the sugars, and ABA contents led the glucose response although ABA seemed to be independent on ABI2, 4, and 5 signaling. The sugar-induced delay of germination is involved in ABI3-mediated ABA signaling, and RGL2 (RGA-like 2) and SPY (SPINDLY) in GA signaling. The effects of exogenous glucose on germination inhibition is studied using abi3, rgl2, and spy5 Arabidopsis mutants [47]. RNA-seq data revealed two RGL2 genes and one SPY gene which showed higher up-regulation in 'Woori' than in 'Keumgang', but these were not selected for investigation due to the gene selection criteria. There were DEGs categorized into flavonoid biosynthetic process and L-phenylalanine metabolic process in group-I and II, respectively, which showed comparatively higher expression in WA than KA whether genes were expressed differentially or not in KP and WP along with drug membrane transport (including ABCG40s in our raw dataset). A total of 6 DEGs were identified as encoding phenylalanine ammonia-lyase 1 (PAL1) in group-I and II (each group had 3 PAL1 genes). PAL1 catalyze the phenylpropanoid pathway and regulate the content of flavonoids. Trehalose is also known to activate PAL to enhance powdery mildew resistance [48,49]. Therefore, the majority of categorized DEGs that showed comparatively higher expression in 'Woori' (DEG analysis group-I to IV) were involved in sugar-ABA interactions and consequently the ABI3-mediated ABA signaling.

Conclusions
In conclusion, a total of 104,542 unigenes were identified from RNA-seq which were used as mapping reference of six samples of artificial PHS-induced and ABA-treated 'Keumgang' and 'Woori'. Among them, 123 unigenes were identified to be related to hormone metabolism and signaling pathways of phytohormones, ABA, GA, IAA, and cytokinin and 1862 of differentially expressed genes. The transcriptome analysis of phytohormone metabolism and signaling related genes and GO categorization of DEGs suggest that the main factor affecting PHS is ABA and ABI-mediated signaling (Figure 7). These results provided the data of transcriptional changes during PHS and consequently provide the clues for understanding processes of PHS and seed dormancy. Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/app11083562/s1, Figure S1: Network analysis of DEGs group-V to VIII, Table S1: Categorization of group-I differentially expressed genes (DEGs), Table S2: Categorization of group-II differentially expressed genes (DEGs), Table S3: Categorization of group-III differentially expressed genes (DEGs), Table S4: Categorization of group-IV differentially expressed genes (DEGs), Table  S5: Categorization of group-V differentially expressed genes (DEGs), Table S6: Categorization of group-VI differentially expressed genes (DEGs), Table S7: Categorization of group-VII differentially expressed genes (DEGs), Table S8: Categorization of group-VIII differentially expressed genes (DEGs), Table S9: Expression of unigenes related in ABA hormone metabolism and signaling, Table S10: Expression of unigenes related in GA hormone metabolism and signaling, Table S11: Expression of unigenes related in IAA hormone metabolism and signaling, Table S12: Expression of unigenes related in cytokinin hormone metabolism and signaling, Table S13: List of selected ABA, GA, IAA, Cytokinin pathway genes and gene identification. Constructed pathway gene specific primers are listed as forward (F) and reverse (R) primers used in qRT-PCR analysis.