To identify the TaPP2C genes in hexaploid wheat, amino acid sequences of all PP2Cs in Arabidopsis and rice were downloaded from The Arabidopsis Information Resource (TAIR¹) and Rice Genome Annotation Project (RGAP²) databases with the National Center for Biotechnology Information (NCBI³) as a complementary database. These sequences were then used as queries to perform BLASTp and tBLASTn searches (threshold e-value < 1e-10) against the T. aestivum reference sequences in the Ensembl Plants database⁴ supported by the IWGSC database⁵. Then, the candidate TaPP2C genes were used as queries to perform BLASTn searches of the wheat genome to obtain more potential genes. The hmmsearch program of the HMMER software⁶ (version 3.2.1) was also applied to the identification of TaPP2Cs using protein phosphatase 2C domain (PF00481, PF07830, and PF13672) in Pfam 32.0 database⁷. The above obtained protein sequences were further screened for conserved domains using SMART⁸ and NCBI Conserved Domains with automatic mode (threshold = 0.01; maximum number of hits = 500)⁹, and proteins without a typical PP2C catalytic domain were removed. Multiple sequence alignment of TaPP2C amino acid sequences was performed using ClustalX 2.1 then a phylogenetic tree was generated using MEGA 6.0 based on the neighbor-joining (NJ) method with 1000 bootstrap replicates (Larkin et al., 2007; Tamura et al., 2011). Nonsynonymous (Ka) and synonymous (Ks) substitution rates were calculated by Ka/Ks Calculator 2.0¹⁰ using the Nei and Gojobori (NG) method (Wang et al., 2010).

To map all TaPP2C genes to wheat chromosomes, the genome annotation file IWGSC RefSeq v1.0 was downloaded from the IWGSC database. Multiple sequence alignment of TaPP2Cs was performed to analyze gene duplication events among the three subgenomes (A, B, and D) of hexaploid wheat. Subsequently, synteny blocks of TaPP2Cs were calculated using MCScanX with e-value ≤ 1e-10 (Wang et al., 2012). Chromosomal locations and syntenic relationships were illustrated using Circos-0.67. To analyze the chromosomal translocation events in the wheat genome, data provided by Ma et al. (2013) and Clavijo et al. (2016) were used. The coding sequences and genome sequences of TaPP2Cs were used to determine the exon-intron structures by Gene Structure Display Server¹¹ (Hu B. et al., 2015). To identify conserved motifs within TaPP2Cs, the MEME motif search tool (Bailey et al., 2009) was applied with an optimum motif width of 6–50 and each motif having 2–600 sites. The results were rearranged by TB tools.

To analyze the expression patterns of TaPP2Cs, RNA-seq data of the project choulet_URGI (Ramírez-González et al., 2018), DRP000768 (Oono et al., 2013), SRP041017 (Zhang et al., 2014), SRP043554 (Li et al., 2015), and SRP045409 (Liu et al., 2015) were downloaded from the expVIP platform¹². Heatmaps were generated from log 2 based transcripts per million (TPM) values using pheatmap package of R project¹³. Expression data of those TaPP2Cs represented in all the three homoeologous subgenomes were chosen to analyze the homoeologous expression patterns using SigmaPlot¹⁴.

The hexaploid wheat (T. aestivum L. cv. Chinese Spring) seeds were surface sterilized, then were soaked in distilled water in a greenhouse (16 h light/8 h dark cycle at 22°C). After 2 weeks, young seedlings were steeped in and sprayed with 200 mM NaCl, 20% (w/v) polyethylene glycol (PEG) 6000, 100 μM ABA, 100 μM gibberellin (GA) and 100 μM methyl jasmonate (MeJA) for 24 h, respectively. The seedlings treated with distilled water were used as the controls. The leaf tissues from seedlings were harvested at six different time points (0, 1, 3, 6, 12, and 24 h) after treatments. All leaf samples, including treated and control samples, were collected with three biological replicates at each time point, and were stored at -80°C till the extraction of total RNA.

Total RNA was extracted from each sample using a Plant Total RNA extraction Kit (Zomanbio, Beijing, China), according to the manufacturer’s instruction, and stored at -80°C. First-strand cDNA was synthesized from total RNA (50 ng–2 μg) in a 20 μl volume using FastKing RT Kit (Tiangen, Beijing, China), according to the instructions. The concentration of total RNA varied from 0.01 to 1 μg/μl. Next, qRT-PCR was performed on a real-time PCR instrument (CFX96; Bio-Rad, Hercules, CA, United States) using AceQ qPCR SYBR Green Master Mix (Vazyme, Nanjing, China). To identify cis-regulatory elements in gene promoters, approximately 2 kb upstream sequences of genes were analyzed via PlantCARE search tool¹⁵.

A total of six TaPP2Cs in group A and ten TaSnRK2s were amplified from the wheat cDNA. The TaPP2C and TaSnRK2 genes were cloned into pGADT7 and pGBKT7 vectors, respectively. Primers for the amplification of TaSnRK2s were obtained from Zhang et al. (2016). Yeast two-hybrid assay was performed according to the manufacturer’s protocol (Clontech, United States) using yeast strain AH109. Positive transformants picked from SD medium lacking leucine and tryptophan (SD/-Leu/-Trp) were subsequently transferred to auxotrophic medium for further selection.

The pSN1301-TaPP2C135 plasmid and pSN1301 empty vector were transformed into Arabidopsis using the floral-dip method with Agrobacterium tumefaciens strain EHA105 (Clough and Bent, 1998). Seeds of transgenic Arabidopsis were selected using Murashige and Skoog (MS) medium (pH 5.8) supplemented with 20 mg/L hygromycin B. Homozygous lines of T₃ and T₄ generations were used for germination analysis. For germination assay, approximately 50–60 seeds were sown on MS plates containing various concentrations of ABA. After stratification for 4 days, the germination greening ratio was scored daily for consecutive 7 days.

After the genome-wide searching and characterization of PP2C catalytic domain, a total of 257 PP2C homoeologs in wheat were identified (Supplementary Table S1). These TaPP2C genes were renamed based on the order of wheat subgenomes (A, followed by B and D), chromosomes (1–7), and positions on each chromosome. Phylogenetic analysis of TaPP2Cs and OsPP2Cs was performed to analyze the evolutionary relationships (Supplementary Figure S1). An individual phylogenetic tree of the TaPP2C proteins was also made to separately check their phylogenetic relationships (Figure 1). The result indicated that these TaPP2C proteins could be divided into 13 groups (A-M) with 11 ungrouped proteins, which was consistent with the PP2C groups found in rice and Arabidopsis.

Hexaploid wheat contains three (A, B, and D) homoeologous subgenomes. Consequently, every hexaploid wheat gene potentially owns three homoeologs from three homoeologous chromosomes as a triad. Blast searches against the wheat genome revealed that the 257 TaPP2C genes represented 257 homoeologs of 95 genes: 76 genes had all the three homoeologs, ten genes had two homoeologs, nine genes had only one homoeolog (Table 1 and Supplementary Table S2). The distribution of TaPP2Cs in each group was similar to that of AtPP2Cs and OsPP2Cs, except for groups A, H, and K: these groups contained more members in wheat than in Arabidopsis and rice (Table 2).

Chromosomal locations as well as syntenic relationships among the A, B, and D subgenomes of TaPP2C genes are illustrated in Figure 2. All TaPP2C homoeologs were mapped to 21 wheat chromosomes, which were highlighted in the middle circle in Figure 2. The TaPP2C genes in different groups (indicated in different colors in Figure 2) showed an uneven distribution across the A, B, and D subgenomes and unbiased distribution among the seven chromosomes of each subgenome. TaPP2C homoeologs involved in chromosome translocation and pericentromeric inversion events were identified (Table 3), and these crosslinks were also represented in the inner circle of Figure 2. Six triads (18 homoeologs) were involved in pericentromeric inversions between the long and short arms of chromosome 4A; three triads (nine homoeologs) were involved in reciprocal translocations between the long arms of chromosomes 4A and 5A; and one triad (three homoeologs) was involved in translocation of 7BS and 4AL. Additionally, it was worth mentioning that synteny analysis on TaPP2C194/207/217 and TaPP2C195/205/218 suggested pericentromeric inversion between the long and short arms of chromosome 6B.

To determine the mode of selection of duplicated TaPP2C genes in groups A-I and K, Ka/Ks ratios were calculated for each gene-pair (Supplementary Table S3). All of the computed gene pairs showed a Ka/Ks ratio < 1, suggesting that TaPP2C genes in these groups underwent purification or negative selection. The Ka/Ks value of three gene pairs (TaPP2C48/63, TaPP2C202/213, and TaPP2C42/57) was zero, indicating strong purifying selection. The average Ka/Ks ratios of gene pairs in different groups ranged from 0.0808 (group G) to 0.3258 (group I), while the Ka/Ks ratios of genes in group F varied from 0 to 0.756.

Exon–intron structural diversity within a gene family is an important clue for the evolutionary and functional analyses of gene family members. To examine the structural features of TaPP2C genes, one homoeolog of each TaPP2C gene was selected and exon-intron structure was analyzed (Figure 3A). The results revealed that genes in the same group shared a similar number of exons but with different exon and intron lengths. A few exceptions were noted in most groups. In group D, while most genes harbored four or five exons, TaPP2C-d8 and TaPP2C-d10 contained only three and two exons, respectively. Genes in group K exhibited wide variation in exon number ranging from one to 12, and TaPP2C-k1, -k2, and -k3 genes contained only one exon. We also examined protein domains and conserved motifs in amino acid sequences of TaPP2Cs. Protein domain analysis showed that most members contained typical PP2C catalytic domains, whereas all six members of group G contained PLN03145, which also belonged to PP2C protein family (Figure 3B).

To identify common motifs among different groups of TaPP2C proteins, we used the MEME motif search tool. Eleven conserved motifs were identified (Table 4); the distribution of these motifs in TaPP2C proteins is illustrated in Figure 3C. The motifs identified in TaPP2C proteins resembled those identified in OsPP2Cs and AtPP2Cs, except motifs 9 and 10. Proteins in the same group exhibited similar motif distribution patterns (Figure 3C). Some motifs (1, 3, 4, 6, and 7) were present in most groups, whereas other motifs were specific to one or two groups; for example motifs 10 and 11 were observed only in groups E and D, respectively. Motif 4 containing [DG]X₂[G](D, aspartic acid; G, glycine; X, any amino acid), is a PPM-type phosphatase signature.

To clarify the biological roles of TaPP2C genes in wheat, expression patterns of all TaPP2C homoeologs were analyzed. RNA-seq data across 15 tissues (five tissues at three different developmental stages) of the wheat cultivar Chinese Spring under non-stress conditions were used to analyze the spatial and temporal expression patterns of TaPP2C genes (Figure 4). There were four homoeologs with no detectable expression, which might merely express at other specific tissues or under special conditions. Most TaPP2C genes exhibited a broad range of expression in stem, spike, root, leaf and grain tissues of wheat plants at different developmental stages.

Row clustering was applied, and as a result, the 257 homoeologous TaPP2C genes fell into five groups (I-V) based on expression patterns. Genes in groups I and IV, which accounted for nearly two-thirds of the TaPP2C genes, maintained similar expression levels in all 15 tissues. Additionally, expression levels of genes in groups I and IV were higher or lower, respectively, than those in other groups. Most genes in groups II and III exhibited moderate expression levels, although the expression level varied at different developmental stages. Expression levels of group II genes including TaPP2C82, 93, and 106 were dramatically higher in leaf and stem tissues than in other tissues at all developmental stages. Notably, TaPP2C82, 93, and 106 identified as homoeologs of the same gene in group K, displayed parallel expression pattern. However, not all homoeologs of the same gene showed similar expression patterns, TaPP2C86, 98, and 111, homoeologs of gene TaPP2C-m1, showed diverse expression patterns in groups II and IV respectively, suggesting functional diversification of homoeologous genes. Expression levels of genes in group V were relative low in most tissues with the exception of particular one or two tissues. Group V contained 20 homoeologs, all of which showed highly preferential expression in grain; TaPP2C41, 56, and 69 showed grain-specific expression especially at 30 days post anthesis. However, no obvious comparable expression pattern was observed within different groups of TaPP2C genes.

To further investigate the potential responses of TaPP2C genes to different stresses, RNA-seq data of four abiotic (heat, drought, cold, and phosphate starvation) and two biotic (stripe rust and powdery mildew) treatments were acquired. These expression profiles were clustered according to the groups of TaPP2C homoeologs to identify the potential biological roles of each group (Figure 5). Several homoeologs with missing expression data under all treatments were observed; these are displayed as blank cells in Figure 5. Comparison of data shown in Figure 5 with those shown in Figure 4 revealed that homoeologs displayed as blank cells either maintained relatively low expression levels in all 15 tissues or preferentially expressed in particular tissues; most RNA-seq data were obtained from leaves of wheat seedlings, which may explain the presence of blank cells (represented by TaPP2C78, 155, and 223). However, blank cells could also be the result of expression induced under a specific condition.

All TaPP2C genes showed varied expression patterns after cold treatment except those representing blank cells (Figure 5). Most of the genes in group A were significantly up-regulated after drought and cold treatments but not in response to heat treatment, with the exception of TaPP2C24 (homoeolog of TaPP2C-a2) and TaPP2C245 (homoeolog of TaPP2C-a11), which distinctly responded to heat stress but to not drought stress. Interestingly, the expression pattern of TaPP2C24 was contrary to that of TaPP2C13 and TaPP2C38 (the other two homoeologs of TaPP2C-a2), probably because of differences in the upstream regulatory regions of these genes. Additionally, most genes in of group A were up- and down-regulated in response to powdery mildew infection and phosphate (Pi) starvation, respectively, indicating their roles in fungal pathogen defense and Pi metabolism. Thus, group A PP2Cs are involved in several plant processes, beyond the ABA signaling pathway. All genes in group B responded to drought stress and pathogen infection but showed no remarkable changes under Pi-stress. Several genes in group C were up-regulated under heat and drought stresses, such as all three homoeologs of TaPP2C-c2. The expression patterns of genes in groups E, H, and I under these stress conditions were similar to those of genes in group D, and 50% of the genes in these groups were significantly up-regulated by heat treatment. Genes in groups F and K exhibited diverse expression patterns under heat and drought treatments, but all of these genes were up-regulated upon fungal disease infection.

To better understand the roles of homoeologs in TaPP2C gene family under different stresses, a comprehensive analysis of expression patterns was conducted (Figure 6). Firstly, genes in possession of triad homoeologs were selected from all TaPP2C genes. Then heat or drought stress treatments with the lowest number of blank cells were selected. Finally, 159 homoeologs (belonging to 53 TaPP2C genes) were put into analysis. Relative expression abundance of homoeologs within triads before and after heat and drought treatments is shown in Figure 6A. A simple impractical assumption was that each homoeolog contributed equally to the total amount of mRNA of a gene (Leach et al., 2014). In fact, nearly one third of TaPP2C homoeologs displayed expression bias within triads under normal condition, and this bias was more evident under heat and drought stresses, suggesting unbalanced functions of homoeologs in stress responses. According to their expression abundance, the unbalanced homoeologous expression patterns were divided into six categories: homoeolog-dominant or homoeolog-suppressed (Figure 6B). Overall, homoeologs from the D subgenome within triads had slightly higher abundance than those from B and A subgenomes, with more D-homoeolog dominance and less D-homoeolog suppression. Moreover, D-homoeolog dominance within triads increased obviously after abiotic stress treatments, especially after heat treatment. Additionally, B-homoeolog dominance within triads was notably enhanced under stresses, particularly drought stress, while A-homoeolog dominance was maintained at the lowest level consistently throughout all the treatments. These results showed significant variation in homoeologous expressions among TaPP2C genes, and stresses either increased or decreased the expression abundance of homoeologs of the same gene.

Several members of group A PP2Cs have been shown to function as negative regulators of ABA signaling pathway in Arabidopsis. To evaluate functions of TaPP2C genes in group A, the expression of these genes was examined by qRT-PCR under ABA, salt and drought stress treatments (Figure 7). Eight genes from group A were randomly selected for this analysis, and specific primers matching one or two homoeologs of the same gene were designed based on the expression patterns of TaPP2C homoeologs described above. Primers used for qRT-PCR were listed in Supplementary Table S4. Under ABA stress, the expression of all eight genes was up-regulated by more than ten-fold, which was consistent with the findings in rice and Arabidopsis. Seven genes were up-regulated under both drought and salt treatments, whereas TaPP2C-a1 showed only subtle changes under both stresses. Expression levels of TaPP2C-a8, -a9, and -a10 genes were dramatically increased after ABA, salt, and drought treatments; thus, these genes are excellent candidates for functional characterization in future studies.

Since group A PP2C genes mediated the crosstalk between ABA and other phytohormones in Arabidopsis (Manohar et al., 2017), we investigated the responses of group A TaPP2C genes to GA and MeJA (Figure 7). The results showed that TaPP2C-a5, -a7, -a9, and -a10 genes were up-regulated to different extents after GA treatment, and TaPP2C-a3, -a7, -a8, and -a9 were down-regulated after MeJA treatment.

To systematically assess the interactions between group A TaPP2Cs and TaSnRK2s, we performed yeast two-hybrid assays. Ten TaSnRK2 members have been isolated previously (Zhang et al., 2016). However, in this study, we found that the originally identified TaSnRK2.7 gene was a homoeolog of TaSnRK2.6 gene, while another gene identified as TaSnRK2 gene was renamed TaSnRK2.7. The re-identified TaSnRK2s were presented in Supplementary Table S5. According to the result of qRT-PCR analysis of group A TaPP2C genes, six corresponding TaPP2C homoeologs were successfully cloned. Primers used were listed in Supplementary Table S4. The results showed that all six TaPP2Cs interacted with one or two members of subclass III TaSnRK2s (Figure 8). However, weak interaction was detected between TaPP2Cs and TaSnRK2.1 and TaSnRK2.2, which belonged to subclasses II SnRK2s.

TaPP2C-a9 was greatly induced by ABA treatment (Figure 7). Therefore, transgenic Arabidopsis plants expressing TaPP2C135 (B homoeolog of TaPP2C-a9) were generated, and lines OE4, OE5 and OE6 were randomly selected for further analysis (Figure 9). Transgenic Arabidopsis of pSN1301 vacant vector (VC) was used as control. Expression levels of TaPP2C135 in the transgenic lines and the wild type were verified by RT-PCR (Figure 9A). The seed germination greening ratios of transgenic lines and the wild type on MS medium with or without ABA treatment were calculated (Figure 9B,C). While on MS medium, no difference was found among TaPP2C135 transgenic lines, VC transgenic lines and the wild type. With increasing ABA concentrations in the medium, the seed germination greening ratio significantly decreased. However, the seed germination greening ratios of TaPP2C135 transgenic lines were higher than those of the wild type and VC controls, especially at higher ABA concentration, indicating that TaPP2C135 transgenic lines were more tolerant to ABA than the wild type and VC controls.