Genome-wide identification of alcohol dehydrogenase (ADH) gene family under waterlogging stress in wheat (Triticum aestivum)

Identification of ADH gene family in wheat

The wheat genome data was downloaded from the wheat genome database (https://urgi.versailles.inra.fr/download/iwgsc/IWGSC_RefSeq_Assemblies/v1.0/) (Appels et al., 2018). We first used 26 melon ADH family protein sequences as search inquiry sequences and searched melon ADH homologous genes using the local BLAST program (amino acid identity > 70%, E value < 10⁻¹⁰). All melon ADH proteins were derived from the literature of Jin (Jin et al., 2016). Secondly, the ADH domain (PF00107.26, PF08240.12 and PF13602.6) was obtained from the Pfam database (http://pfam.xfam.org/) (Finn et al., 2010). We first used the HMMER3 software package to create the Hidden Markov Model file (Eddy, 2011), and then used HMMsearch with default parameters to search the wheat protein database. All the above candidate proteins were verified by Batch-CDD (https://www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb.cgi) and SMART (http://smart.embl-heidelberg.de/) (Letunic, Doerks & Bork, 2012), and all proteins without GroES-like domain and zinc-binding domain were deleted.

The number of amino acids,protein molecular weight and isoelectric point of the candidate proteins were calculated by the ExPASy website (https://web.expasy.org/compute_pi/). Plant-mPLoc (http://www.csbio.sjtu.edu.cn/bioinf/plant-multi) (Fang et al., 2020) was used to predict the subcellular location of the TaADH genes.

Evolutionary analysis of TaADH genes

To clarify the evolutionary relationship of ADH protein in wheat and the evolutionary relationship of ADH protein between wheat and several other species, all the ADH protein sequences were aligned by the ClustalW program and then used to constructed the phylogenetic tree by MEGA7.0 (Kumar, Stecher & Tamura, 2016). All the phylogenetic trees were constructed by the neighbour joining (NJ) method, 1000 replicates of bootstrap values, the pairwise deletion option. ADH protein sequences from Arabidopsis thaliana, Cucumis melo, Cucumis sativus, Glycine max, Hordeum vulgare, Lycopersicon esculentum, Oryza sativa and Vitis vinifera were downloaded from the NCBI (National Center for Biotechnology Information) database by using their gene IDs from the Pyrus bretschneideri reference (Jin et al., 2016).

Analysis of structural characteristics of TaADHs

To clarify the structural characteristics of the ADH gene, the intron-exon distribution map of TaADH genes was generated by the Gene Structure Display Server 2.0 (GSDS2.0, http://gsds.gao-lab.org/) (Hu et al., 2015). The generation of the intron-exon distribution map depends on the cDNA sequence and the corresponding genomic DNA sequence of the wheat ADH gene.

To clarify the domain of wheat ADH protein, we used DNAMAN software to carry out the multiple sequence alignment of all TaADH protein sequences. Besides, we used Multiple EM for Motif Elicitation (MEME, http://meme-suite.org/tools/meme) (Bailey et al., 2006) to analyze the motif of TaADH proteins. The program was set according to the following parameters: the optimal width of motif is 650 amino acid residues, and the maximum number of motifs is 15. Finally, the motifs map of TaADH protein were presented by TBtools (https://github.com/CJ-Chen/TBtools) (Chen et al., 2020).

Distribution of TaADHs on chromosomes and gene duplication events

To determine the position of TaADHs on the chromosome, the starting position of TaADHs was extracted from the Chinese Spring wheat genome (https://urgi.versailles.inra.fr/download/iwgsc/IWGSC_RefSeq_Assemblies/v1.0/) and finally presented by TBtools (Chen et al., 2020).

To identify the TaADH gene duplication events, the open reading frame of all TaADHs were BLAST compared with each other by the local BLASTN program (Identity > 80%, e-value < 1e⁻¹⁰). Gene alignment coverage was then acquired using the previously calculated method: gene alignment coverage = (alignment length-mismatch length)/the length of larger genes (Wang et al., 2010). When the gene alignment coverage was more than 0.75, they were considered to be a duplicated gene (Wang et al., 2010). Besides, in the 100 kb region, two duplicated genes separated by other genes were considered as tandem duplicated genes. When the distance between two duplicated genes was more than 100 kb or the duplicated genes were distributed on different chromosomes, they were named as fragment duplicated genes. Non-synonymous substitution rate (Ka), Synonymous substitution rate (Ks) and Ka/Ks were calculated by DnaSP software (http://www.ub.edu/dnasp/) (Rozas et al., 2003). The formula: T = Ks/2λ × 10⁻⁶ Mya was used to calculate divergence time (T), λ = 6.5 × 10⁻⁹ represented the rate of divergence of synonymous substitutions per site per year, the unit of evolution time was millions of years (Mya) (Emanuelsson et al., 2000).

Analysis of tissue expression pattern of TaADHs

To analyze the expression pattern of TaADH genes in root, stem, leaf, spike and grain, we downloaded the RNA-seq reads data from Chinese Spring wheat through expVIP (an expression visualization and integration platform) platform (http://www.wheat-expression.com/). Raw data was deposited as DRP000768 and SRP028357 (https://urgi.versailles.inra.fr/files/RNASeqWheat/) on NCBI and were derived from the Chinese Spring (123 non-stressed samples, 15 tissues) (Borrill, Ramirezgonzalez & Uauy, 2016). TPM (transcripts per kilobase of exon model per million mapped reads) values were used as expression units on expVIP software. After obtaining the average TPM expression level of all TaADH genes, the heatmap was drawn with TBtools software (https://github.com/CJ-Chen/TBtools) (Chen et al., 2020).

Cis-acting element analysis of TaADH genes

To analyze the cis-acting element of TaADH genes, the promoter sequences (2,000 bp before the start codon) of all TaADH genes were extracted from the Chinese Spring wheat database. The cis-acting elements of these genes were predicted by online tool PLANTCARE (https://bioinformatics.psb.ugent.be/webtools/plantcare/html/) (Lescot et al., 2002) and visualized by TBtools (https://github.com/CJ-Chen/TBtools) (Chen et al., 2020).

The functional enrichment analysis was performed using g:Profiler (version e102_eg49_p15_7a9b4d6) (https://biit.cs.ut.ee/gprofiler/gost) with g:SCS multiple testing correction method applying significance threshold of 0.05 (Raudvere et al., 2019). The 22 TaADH gene IDs were uploaded into the program, and ‘Triticum aestivum’ was chosen as the reference organism to analyze molecular function, cellular components, and biological processes.

Expression analysis of TaADHs under waterlogging stress

To study the waterlogging tolerance of TaADHs, ‘Zhoumai 22’ (ZM22) and ‘Bainong 607’ (BN607) were used as materials in this study. The materials were cultivated and provided by professor Ou’s team of the School of Life Science and Technology of Henan Institute of Science and Technology.

Experimental group: 50 seeds of each variety were submerged in the 14-cm-diameter glass Petri dish and filled with 200 mL of sterilized deionized water (pH 6.8; electrical conductivity 1.5 μS cm⁻¹), and enveloped in aluminium foil to minimize gas exchange at 20 °C in the dark for 3 days. The seeds were germinated for 24 and 72 h (cultured in a Petri dish containing a layer of filter paper and 10 mL of aseptic deionized water) after waterlogging treatment. Control group: the seeds of each variety germinated normally without waterlogging treatment. All the Petri dishes containing seeds were placed in a growth chamber with 25 °C, 75% relative humidity, 16 h light/8 h dark cycle. The experimental group and the control group were repeated five times. The characteristics of ZM22 and BN607 were identified by phenotypic analysis and bud length measurement between the experimental group and the control group.

To determine the response of TaADH genes to waterlogging, we used the qRT-PCR to determine the expression level of ADHs. The seeds were collected at 24 and 72 h after waterlogging stress (the buds and roots of the seeds were removed), and the RNA was extracted by RNAsimple Total RNAKit (Tiangen, Beijing, China). To avoid the contamination of genomic DNA, RNase-free DNase I (Takara, Tokyo, Japan) was used to remove DNA from the total RNA. The first-strand cDNAs were synthesized using the PrimeScript RT Reagent Kit with gDNA Eraser (Takara, Dalian, China) according to the manufacturer’s protocol. The qRT-PCR assays were performed with the Primer Script RT Reagent Kit (Takara, Dalian, China) and 18S (AJ272181.1) was used as a reference gene. The quantitative primers of all TaADH genes were designed by Primer 6.0 software, and BLAST searched against the wheat database to determine the specificity of primers (Table S1). The PCR conditions were as follows: 95 °C for 10 s and 40 cycles of 95 °C for 5 s and 60 °C for 30 s. The qRT-PCR was performed using an ABI Step One Plus. All the experiments were performed with three biological replicates. The relative expression was calculated using the 2^−ΔΔCt method (Livak & Schmittgen, 2001).

Genome-wide identification and physicochemical characteristics of ADH genes in Triticum aestivum

A total of 22 ADH genes were identified in the wheat genome based on the BLAST program. According to the location distribution of these genes on chromosomes (from Ta1A, Ta1B, Ta1D to Ta7A, Ta7B, Ta7D, from top to bottom), they were named TaADH1~TaADH22 (Table 1, Table S2). The length of predicted coding sequences of 22 TaADH genes ranged from 1,044 to 1,244 bp, the number of corresponding amino acids ranged from 347 to 415 aa. The theoretical isoelectric point (PI) ranged from 5.68 to 8.2, and the molecular weight ranged from 34.4 to 44.2 kDa. Through the subcellular localization prediction of TaADH genes, it was found that they were all localized to the cytoplasm.

Classification and conserved domain analysis of 22 TaADHs

To analyze the evolutionary relationship of 22 TaADH proteins, the phylogenetic tree was constructed by MEGA7.0. According to the amino acid sequence identity, 22 TaADH proteins were divided into two subfamilies (subfamily I, subfamily II) (Fig. 1A). Subfamily I contained the largest number of (10) TaADH proteins, and the remaining three TaADH proteins (TaADH20~22) belonged to subfamily II.

To analyze the conserved motif of TaADH protein, we used MEME online software to analyze the 22 TaADH protein sequences. A total of 14 motifs were detected in the wheat ADH family (Fig. 1B, Fig. S1). Motif 1, 3, 4 and 5 existed in all TaADH proteins, indicating that these TaADH proteins were highly conserved (Fig. 1B). Motif 2, 8, 9 and 11 existed only in subfamily I, while motif 12, 13 only existed in subfamily II (Fig. 1B), indicating structural differences and specificity among different subfamilies. Motif 7 was only missed in TaADH1, motif 11 was only missed in TaADH13, TaADH9 and subfamily II, and motif 10 was only missed in TaADH14 and subfamily II.

Comparing the sequences of 22 TaADH proteins, all of these protein sequences contained conserved GroES-like domain and Zinc-binding domain (Fig. S2). GroES-like domain contained 35~164 amino acids, and the Zinc-binding domain contained 206~340 amino acids, which correspond to the structural characteristics of the ADH gene family.

Exon-intron analysis of 22 TaADHs

Through the analysis of intron-exon structure, we found that 22 TaADH genes contained 8~10 exons (Fig. 2). In subfamily I, three genes (TaADH5, TaADH8 and TaADH11) contained six exons, two genes (TaADH13 and TaADH16) contained eight exons, and the remaining TaADH genes contained nine exons (Fig. 2). In subfamily II, all TaADH genes contained eight exons (Fig. 2). As a whole, we found that the same branch contained similar structural patterns.

Chromosomal distribution and gene duplication of 22 TaADHs

According to the initial position of 22 TaADHs on the chromosome, they were visualized by TBtools. Twenty-two TaADHs were unevenly distributed on 15 of 21 chromosomes (they do not exist on chromosomes Ta2A-D and Ta3A-D) (Fig. 3A). There were three genes on chromosomes Ta4A-D, respectively; two genes on chromosome Ta1A; and one gene on the other chromosomes, respectively (Fig. 3).

By analyzing the duplication events of 22 TaADHs, we found 64 pairs of duplicated genes (Fig. 3, Table S3). Most of these duplicated genes were on different chromosomes, TaADH1-TaADH2, TaADH1-TaADH3, TaADH6-TaADH7 and TaADH9-TaADH10 were on the same chromosome, respectively. However, we found that only the distances of TaADH1-TaADH2 and TaADH9-TaADH10 were less than 100 kb, indicating that the two duplicated gene pairs had tandem duplication events, and the remaining duplicated gene pairs had undergone fragment duplication events (Table 1, Table S3).

To reveal the selection pressure of TaADH family genes in the process of evolution, the non-synonymous substitution rate (Ka), synonymous substitution rate (Ks) and Ka/Ks for 64 duplicated pairs were calculated (Fig. 3B, Table S3). All Ka/Ks of these duplicated pairs were less than 1, which tended to a pure selection, indicating that the sequence similarity of TaADH genes was very high and relatively conservative in the process of evolution. The evolution time of the duplicated events of TaADH genes can be divided into three evolution periods (Fig. 3C, Table S3). The first period was 11.19~16.42 million years ago (Mya), with 30 duplicated gene pairs. Twelve TaADH duplicated gene pairs occurred at 7.56~9.56 Mya. The remaining 22 TaADH duplicated gene pairs occurred at 0.9~5.53 Mya. Although these gene sequences were conserved, they were different in evolutionary time.

Phylogenetic relationship of ADHs in Triticum aestivum and A. thaliana

To study the evolutionary relationship between Triticum aestivum (22) Arabidopsis thaliana (7), Cucumis melo (13), Cucumis sativus (12), Glycine max (3), Hordeum vulgare (1), Lycopersicon esculentum (7), Oryza sativa (1) and Vitis vinifera (8), a phylogenetic tree was constructed. According to the amino acid sequence identity, all ADH proteins were divided into three subfamilies (subfamily I, II, III) (Fig. 4), which was consistent with the classification of TaADH in Fig. 1A. TaADHs existed only in subfamily I and subfamily II. Subfamily I contained the most (59) ADH proteins, subfamily III contained the least (5) ADH proteins, and the remaining proteins (10) belonged to subfamily II. Besides, according to the number of amino acid residues, all ADH proteins in subfamily I belonged to medium-chain ADH proteins, 3 ADH proteins from subfamily II belonged to long-chain-ADH, the remaining proteins from subfamily II belonged to medium-chain ADH proteins (Fig. 4). All TaADHs belonged to medium-chain ADH protein (Fig. 4).

The same subfamily usually contained ADH proteins of several species. For instance, subfamily II contained ADH proteins from Arabidopsis thaliana, Cucumis melo, Cucumis sativus, Glycine max, Triticum aestivum and Vitis vinifera, which indicated that these species came from the same ancestor a long time ago.

Tissue expression patterns of TaADH genes

According to the transcriptome data of different wheat tissues, we analyzed the tissue expression pattern of 22 TaADH genes in roots, leaf, stem, spike and grain (Fig. 5). Nine genes (TaADH8, TaADH11, TaADH18, TaADH17, TaADH19, TaADH5, TaADH22, TaADH20 and TaADH21) were highly expressed in all tissues. Six genes (TaADH12, TaADH13, TaADH1, TaADH9, TaADH2 and TaADH15) were not expressed or low expressed in all tissues; the remaining genes were only highly expressed in specific tissues, such as TaADH4 and TaADH6 were highly expressed in grain, but not or low expressed in other tissues.

For the duplicated genes, we found that they had the same tissue expression patterns, such as TaADH8_TaADH11, TaADH17_TaADH18, TaADH20_TaADH21, and so on. There were also differences in the expression patterns of some other duplicated genes (Fig. 5). For example, in the duplicated gene pair TaADH9_TaADH16, TaADH9 was only expressed in grain and root. In contrast, TaADH16 was expressed in all tissues, and the expression profiles in grain and root were lower than that in other tissues (Fig. 5), so it was speculated that duplicated genes had been diversified in the process of evolution.

The cis-regulatory element analysis and GO annotation of TaADH genes

To analyze the potential function of 22 TaADH genes, we analyzed the cis-acting elements of TaADH gene promoters. Through the analysis of 2,000 bp before the start codon, we found that a total of 563 cis-acting elements, which responded to 11 stresses (biotic and abiotic), including hormone response, anaerobic response, defence, and stress response, drought induction, light response, low-temperature response, etc. (Fig. 6, Table S4). Seventeen genes (77%) respond to ABA, of which the TaADH4 gene contained eight abscisic acid responsiveness elements (ABREs), which suggested that the gene may play a key role in the abscisic acid response. The TGACG-motif and CGTCA-motif were methyl jasmonate response elements (Chen & Qiu, 2020). TaADH3 each contained seven TGACG-motifs and seven CGTCA-motifs, so it was speculated that this gene played a key role in methyl jasmonate response. Anaerobic response element (ARE) contained a conserved AAACCA sequence, mainly involved in anaerobic induction. In this study, we found that except for TaADH13 did not contain ARE, all other TaADH genes contained ARE elements. For example, TaADH6 and TaADH9 contained six AREs, respectively, which suggested that TaADH family genes may play key roles under anaerobic stress.

By analyzing the functional characteristics of TaADH genes, namely ‘molecular function’, ‘biological process’, and ‘cellular components’, it can help us understanding the function of proteins at the molecular level. For the category of molecular function (Fig. S3, Table S5), in addition to TaADH1 (TraesCS1A02G370100.1), other 21 TaADH genes were involved 15 GO categories, and they were all associated with “zinc ion binding” (GO:0008270), “oxidoreductase activity” (GO:0016491), “metal ion binding” (GO:0046872), “cation binding” (GO:0043169), “ion binding” (GO:0043167) and “catalytic activity” (GO:0003824). TaADH3 (TraesCS1B02G389200.1) ~ TaADH8 (TraesCS4B02G106300.1) and TaADH11 (TraesCS4D02G103000.1) were involved in “alcohol dehydrogenase activity” (GO:0004022) and “oxidoreductase activity” (GO:0016491). From the perspective of ‘biological process’ categories, all TaADHs participated in 30 GO categories, and they were associated with “oxidation-reduction process (GO:0008152)” and “metabolic process (GO:0055114)”. TaADH4 (TraesCS1D02G376300.1) ~ TaADH8 and TaADH11 genes participated in 24 GO categories. Moreover, TaADH17 (TraesCS6A02G386600.1) ~ TaADH19 (TraesCS6D02G371200.1) were involved with “ethanol metabolic process (GO:0006067)”, “ethanol oxidation (GO:0006069)”, “primary alcohol metabolic process (GO:0034308)” and “alcohol metabolic process (GO:0006066)”. The ‘cell components’ annotation predicted the TaADH4 ~ TaADH8 and TaADH11 existed in cytosol (GO:0005829).

Phenotypic analysis of ZM22 and BN607 under waterlogging stress

ADH gene plays a key role in plants under anaerobic conditions (Bailey-Serres & Voesenek, 2008). To study the expression of the ADH gene (Fig. 7, Table S6) in wheat seeds at 24 and 72 h after waterlogging for 3 days, we selected the seeds of two wheat varieties (Zhoumai 22 (ZM22) and Bainong 2008 (BN607)) for analysis.

At 24 h after the waterlogging treatment, the bud length of ZM22 and BN607 had no significant difference with the control seeds (no waterlogging treatment). However, all the seeds of ZM22 and BN607 germinated at 72 h after the waterlogging treatment. The bud length of ZM22 was significantly lower than that of the control seeds, while the bud length of BN607 was not significantly different from that of the control seeds (Fig. 7, Table S6), which further indicated that BN607 was more tolerant to waterlogged stress than ZM22.

Expression analysis of TaADH genes under waterlogging treatment

The expression profiles of 22 TaADH genes in ZM22 and BN607 seeds at 24 and 72 h after waterlogging and control treatment were analyzed (Fig. 8, Table S6). The results showed that in ZM22, the expression of TaADH1/2, TaADH13, TaADH17, TaADH18, TaADH19 and TaADH20 was significantly up-regulated at 24 h of germination. However, among these genes, only the TaADH13 gene was significantly induced at 72 h after waterlogging. There was no significant difference in the relative expression level of TaADH1/2, TaADH17, TaADH18, TaADH19 and TaADH20 between the waterlogging and control treatment at 72 h after waterlogging.

In BN607, the expression of TaADH1/2, TaADH3-6, TaADH8-13, TaADH19 and TaADH20 was significantly up-regulated at 24 h after waterlogging. Besides, TaADH1/2, TaADH3 and TaADH9 genes were also significantly induced at 72 h after waterlogging. Moreover, the relative expression levels of TaADH5, TaADH6, TaADH14 and TaADH16 were significantly lower than that of the control at 72 h after waterlogging (Fig. 8, Table S6). Based on the above analysis, we found that the response of waterlogging-sensitive ZM22 and waterlogging-tolerant BN607 to waterlogging stress at different times mainly depended on the time-specific expression of some key ADH genes in wheat.