Molecular characterization of wheat TaSKs

A cDNA fragment encoding a protein with high identity to the mammalian Glycogen synthase kinase 3 (GSK-3) and to the Drosophila serine/threonine kinase SHAGGY (SGG) was isolated in the screen of an embryonic cDNA library constructed by means of a suppression subtractive hybridization (SSH) approach [28]. Using this fragment, two new cDNA sequences were obtained by SMART RACE cDNA amplification and named Triticum aestivum Shaggy-like Kinase 1 and 2 (TaSK1 and TaSK2). As a part of the 5′ of TaSK1 could not be cloned by means of the latter technique, additional cloning was performed. Thus cloning followed by alignment of 23 cDNA and 56 genomic clones of TaSK1 as well as 18 cDNA and 21 genomic clones of TaSK2 provided evidence for the occurrence of three expressed gene copies of TaSK1 and TaSK2 named TaSK1-A,B,C and TaSK2-A,B,C. A manual approach and the algorithms/programs CLUSTALX, MAFFT, MUSCLE, Figtree and Quicktree were utilized to assemble, align and subgroup these cloned sequences [29–32]. Genomic and cDNA consensus sequences of TaSK1-A,-B and -C, and of TaSK2-A,-B and -C were extracted from the alignments.

Consensus genomic sequences of TaSK1-A,B,C had a size of respectively 4436, 4422 and 4195 bps. Intron 1 of TaSK1-C could not be cloned. The sizes of the consensus genomic sequences of TaSK2-A,B,C were 3825, 3999, and 3824 bps respectively. Their genomic structure is similar to the one reported for Arabidopsis ASKs [12] namely 12 exons interrupted by 11 introns.

TaSK1-A,C and TaSK1-B predicted proteins considering the longest open reading frame (ORF) contained respectively 400 and 401 amino acids. Their calculated molecular weight was respectively 44.9 and 45.0 kDa. The three TaSK2 consensus cDNAs encoded predicted proteins of 402 amino acids (longest ORF) with a molecular weight of 45.2 kDa. Identity among TaSK1-A,B,C was ranging from 98.8 to 99% while identities among TaSK2-A,B,C were ranging from 99.3 to 99.5% (Table 1C). In comparison, TaSK1 and TaSK2 displayed 88.3 to 88.8% identity (Table 1C). Consequently, identities among the three TaSK1 and among the three TaSK2 were higher as the ones between TaSK1 and TaSK2.

Copy number and chromosome localization of TaSKgenes in the hexaploid wheat genome

Global alignment provided strong evidence for the presence of three expressed gene copies of TaSK1 as well as three expressed copies of TaSK2. The complexity of the hexaploid wheat genome gives rise to the question whether these copies were homoeolog gene copies and/or paralog genes.

A PCR-product was obtained with specific primers for TaSK1-A in all lines tested except in the line N3B-T3D providing strong evidence for the localization of TaSK1-A on chromosome 3B (Figure 1A). Amplification was obtained with primers specific for TaSK1B for all lines except N3D-T3A strongly suggesting a localization of TaSK1-B on chromosome 3D (Figure 1A). Similar approaches led to the conclusion that TaSK1-C most probably is located on chromosome 3A as the only line without amplification is N3A-T3B (Figure 1A).

Similarly, TaSK2-A probably localizes on chromosome1B while TaSK2-B is located on chromosome 1A (Figure 1B). Unfortunately such a Polymerase Chain Reaction (PCR) approach based on sequence specific primers was not suitable for TaSK2-C as this sequence did not have an appropriate specific nucleotide insertion or deletion that distinguished it from the other ones. However exon 4 of TaSK2-C carried a conserved nucleotide exchange that created an Rsa1 restriction site absent in the other TaSK2 and in TaSK1 sequences. Specific primers for TaSK2 sequences were designed upstream and downstream of this restriction site (Table 2). Amplification followed by Rsa1 digestion gave rise to a digestion product in all lines tested (arrows) except in line N1D-T1B suggesting that the digestion site was absent in this line (Figure 1B). We therefore conclude that TaSK2-C is located on chromosome 1D.

In summary, TaSK1-A,B,C were located on the homoeologous chromosomes 3B, 3D, 3A while TaSK2-A,B,C were located on the homoeologous chromosomes 1B, 1A, 1D. Most probably each gene copy was present only on one chromosome. These data strengthen the results of the sequence alignment analysis and provided strong evidence that the three TaSK1 copies on one hand and the three TaSK2 copies on the other hand were homoeolog gene copies.

Triticum aestivumTaSKs display a GSK3/SGG signature

The catalytic domains of TaSK1-A,B,C and TaSK2-A,B,C shared high sequence identity with the catalytic domains of Arabidopsis thaliana BIN2 (91-90%), Drosophila melanogaster SHAGGY (67-68%) and Homo sapiens GSK-3β (70-71%). As comparison, human GSK-3β and Drosophila SHAGGY showed 85% identity in their catalytic domain. The mentioned percentages were pairwise alignment scores obtained by means of ClustalW2.

All TaSKs contained the highly conserved motifs CDFGSAK and SYICSR (AYICSR in the case of TaSK2s) that distinguish the GSK-3 subfamily among serine/threonine protein kinases (Figure 2) [23, 34].

Within the latter motif, all TaSKs had a tyrosine (Tyr) residue in equivalent position to the Tyr 216 of GSK-3β (Figure 2). Phosphorylation of this residue is implicated in the modulation of kinase activity in human and in Arabidopsis[20, 35–37].

Residues in equivalent position to Arg 96, Arg 180, Lys 205 of GSK-3β were present in TaSKs (Figure 2). In the case of GSK-3β, these residues define a pocket for binding of primed substrates [35, 36, 38]. Pre-phosphorylated (primed) substrate by another kinase binds to this pocket and is thereby correctly positioned for a phosphorylation by GSK-3β [35, 36, 38]. However, although Arabidopsis BIN2 contains this pocket, its phosphorylation activity is not based on priming phosphorylation but rather requires a direct interaction with BRZ1 [39].

Inhibition of GSK-3β in the insulin signaling pathway relies on N terminal residue serine 9 (Ser9) phosphorylation [38]. Ser9 was absent in TaSKs as it is the case for BIN2 (Figure 2) [7].

TaSKs were classified in group II due to the presence of the SIDIW box characteristic for plant group II GSKs (Figure 2) [23]. Like Arabidopsis BIN2, TaSKs contained the TREE motif (Figure 2). Almost all Bin2.1 gain of function mutations localize to this motif [13, 15, 19]. Bin2.1 protein was shown to be more stable than its wild type form [19]. TaSKs contained also the motif MEYV that contains key residues for docking of Arabidopsis GSK inhibitor Bikinin (Figure 2) [40]. This inhibitor poorly inhibits human GSK-3β whose motif in equivalent position is LDYV [40].

Protein sequence analysis provides strong evidence for a classification of TaSKs in the GSKs subfamily of protein kinases.

TaSK1 and TaSK2 are functional kinases

An in vitro kinase activity assay was performed to evaluate whether TaSK1 and TaSK2 had a kinase activity. Full length TaSK1 and TaSK2 (longest ORF), Arabidopsis BIN2, OsGSK7 (TaSKs homolog in Oryza sativa) and wheat TaGSK1 were overexpressed in E. coli as GST fusion proteins and affinity purified in native conditions. Transphosphorylation activities towards a bovine myelin basic protein fragment (MBP) that contains several consensus phosphorylation sites for kinase proteins was assayed using radiolabeled ATP (Figure 3). TaSK1, TaSK2, OsGSK7, BIN2 and TaGSK1 expressed as fusion proteins were phosphorylating the MBP fragment while GST alone was not able to phosphorylate MBP (Figure 3, top panel). Furthermore, a long exposure of SDS-PAGE gel to the X-ray film showed a clear autophosphorylation activity for TaSK1, BIN2, OsGSK7 and TaGSK1 (Figure 3, lower panel) while a weaker signal was observed for GST-TaSK2.

These data indicate that cloned TaSK1 and TaSK2 were functionally active kinases.

TaSK1 and TaSK2 belong to clade II of plant GSKs

The genome of the core eudicotyledonous Brassicaceae Arabidopsis thaliana contains 10 different GSKs. These ASKs have been grouped based on their sequences into either three [12] or four [41] major clades. Little is known about Liliopsida, resp. Poaceae GSKs and their phylogenetic relationship in comparison to Arabidopsis. Therefore, the phylogenetic relationship of GSKs of selected Poaceae to Arabidospsis ASKs and to GSKs belonging to other selected eudicotyledons was investigated and analysed.

The Poaceae family includes 12 subfamilies among them the Pooideae, Ehrhartoideae and the Panicoideae that provide the bulk of human nutrition. Besides Triticum aestivum (wheat), the Pooideae contain species like Hordeum vulgare (barley) and Brachypodium distachyon the latter being proposed recently as new grass model system [42]. Oryza sativa (rice) and Zea mays (maize) were selected as representatives of the Ehrhartoideae and Panicoideae, respectively.

Arabidopsis thaliana is a paleopolyploid that has been subjected to two additional whole genome duplications (WGDs) called α and β events after the so-called γ event, the latter being a triplication event resulting in the hexaploid common ancestor of many or most angiosperms [43–45]. Therefore, the GSK sequences of two core eudicotyledons, the Brassicale Carica papaya and the basal Rosid Vitis vinifera, both not exposed to the more recent α and β WGD events specific to Brassicaceae or to any other events after the γ duplication, were included in this study [46–48]. In addition, GSK sequences of Aquilegia coerulae a member of the basal-most or stem eudicotyledons (Ranunculales) and the moss Physcomitrella patens as representative of non-seed plants were added to this phylogenetic analysis. The genomes of all these plants are fully sequenced.

Except for published rice GSK and Arabidopsis ASK gene sequences [7, 41], the sequences of the other GSKs were identified in different databases by means of annotation mining and BLAST (Basic Local Alignment Search Tool) searches (Additional file 1).

Besides TaSK1-A,B,C and TaSK2-A,B,C only 4 other wheat GSK sequences have been identified in the databases. Considering the complexity of the wheat genome due to its size (16,000 Mb) and polyploidy, this number appears low. Probably more sequences will be identified once full genome sequencing data will be available in open access databases. Twenty-eight different maize GSKs were found in the databases (Additional file 1). However, among them 10 subgroups were distinguished in which identities between the predicted proteins were ranging from 97 to 99%. Maize is diploid and the tissues used to generate the cDNAs were derived from different hybrids [49, 50]. Therefore, we hypothesized that the predicted proteins within each group may be the same and that the difference observed may be due to strain polymorphism or might be caused by sequencing artifacts inherent to high throughput sequencing approaches, respectively difference in gene structure prediction. As a consequence, one accession within each group was chosen as representative for this group (Additional file 1).

As already mentioned, previous phylogenetic analyses led to the identification of different plant GSK clades named I to IV [12, 41]. This clade nomenclature was kept here based on the Arabidopsis and rice GSK clade members. With regard to these previously defined subfamilies, the Neighbour-Joining (NJ) tree supported clade II (in blue), clade IV (green) and clade I (yellow) as monophyletic with bootstrap support 89–99 (Additional file 2), while the Maximum Likelihood (ML) tree supported clade III (purple), clade IV (green) and clade I (yellow) with 59–72 quartet support (Additional file 3) and the Bayesian Inference (BI) tree supported the clade III (purple), clade II (blue) and clade IV (green) with posterior probabilities of 0.84-0.99 (Figure 4). Thus, clade IV (green) is supported by all three methods, while the blue (II), yellow (I) and purple (III) clades were supported as monophyletic by two out of three methods each.

In all three trees, TaSK1-A,B,C and TaSK2-A,B,C sequences clustered reliably in respectively two closely related subclades together with related sequences from other grasses (Figure 4, Additional file 2 and 3). In the BI and NJ trees, these two subclades containing TaSK1-A,B,C and TaSK2-A,B,C were embedded into a monophyletic clade II (blue) (Figure 4, Additional file 2). Interestingly clade II includes ASKeta/BIN2, and its close relatives ASKdzeta and ASKiota (Figure 4), all three being involved in brassinosteroid signaling [51]. None of the wheat GSKs identified so far in the databases was classified in clade III that includes ASKtheta shown to be also involved in brassinosteroid signaling [22]. Nevertheless all the other Poaceae selected had at least one GSK representative in this clade (Figure 4).

The presence of the Physcomitrella patens paralogs exclusively at the base of the green clade, supported by all three analyses, suggests that the clade IV (green) was the ancestral one (present already in the earliest land plants), from which all other family members have subsequently evolved (Figure 4, Additional file 2 and 3). Interestingly, lineage-specific gain and retention of paralogs have occurred in Physcomitrella patens, bringing the extant number of paralogs to seven, comparable to the numbers in seed plants (see below). Three paralog pairs are probably derived from the WGD known to have occurred in this lineage [52].

The general theme as revealed by the phylogeny is that each of the eudicotyledon species is represented with one (green, purple), two (blue), or two/three (yellow) sequences per subfamily (Figure 4). In the green (ancestral) clade no A. coerulea sequence is present. Most likely, this sequence has been subject to a secondary loss or there is no adequate gene model present for this sequence. The purple clade lacks a C. papaya sequence, however, this sequence has been removed due to a truncated gene model (cf. Methods). The yellow clade harbours three V. vinifera genes and two each of the other three species, the tree topology suggesting secondary loss of the latter two (Figure 4). In summary, these results point out the presence of 6 (possibly 7) GSKs in the ancestral eudicotyledonous genome after the gamma event, resp. after the separation from the Liliopsida.

The phylogeny of the Liliopsida, considering fully sequenced genomes (all except H. vulgare and T. aestivum), indicate that each of these species has in general one (purple, green), three (yellow) or three-/four (blue) GSK sequences per clade. Clade I (yellow) does not include a third Zea mays GSK. This sequence has been removed because it was incomplete (cf. Methods). A fourth rice GSK was identified in clade II. Taken together these data indicate that 8 GSKs were present in the genome of the Liliopsida-ancestor at the time it diverged from the eudicotyledons.

Cloning and sequencing of TaSKcDNA and genomic clones

A cDNA fragment of TaSK1 was originally isolated by screening of an embryonic cDNA library constructed by means of suppression subtractive hybridization method (SSH). For the SSH, total RNA was isolated from embryo material of Triticum aestivum cv Sonora using TRIzol® Reagent (Invitrogen). Dynabeads Oligo (dT)₂₅ from Dynal were used to purify mRNA. SSH was performed using the PCR-select cDNA substraction kit (Clontech) according to the suppliers instructions. 5′ and 3′ ends of the gene fragment were generated by means of BD SMART™ RACE cDNA Amplification Kit (Clontech) as recommended by the manufacturer.

Several cloning approaches were used to obtain genomic and cDNA sequences of TaSK1-A,B,C and TaSK2-A,B,C.

Total RNA and genomic DNA were extracted from Triticum aestivum cv Sonora tissues using respectively the RNeasy® Mini Kit and the DNeasy® Plant Mini Kit (QIAGEN) as recommended by the manufacturer.

Reverse transcription was performed using either the SuperScript™ II Reverse Transcriptase (Invitrogen), the RevertAid™ H Minus First Strand cDNA Synthesis Kit (Fermentas) or the BcaBEST™ RNA PCR Kit (Takara) as recommended by the manufacturers. Standard PCR procedure or PCR using BcaBEST™ RNA PCR Kit (Takara) were used to amplify the cDNAs.

Genomic DNA was amplified either by means of conventional PCR that may include the use of BcaBEST™ RNA PCR Kit (replacing cDNA by genomic DNA), inverse PCR [69] or thermal asymmetric interlaced PCR [70, 71].

PCR products were cloned in the pCR®II-TOPO®, pCR®2.1-TOPO® and pCR®-Blunt® vectors (Invitrogen). Sequencing was outsourced to GATC Biotech AG, Agowa GmbH, Eurofins MWG GmbH, and Sequence Laboratories Göttingen GmbH.

Sequenced TaSK1 / TaSK2 cDNA and genomic clones were assembled, aligned, subgrouped manually and by means of CLUSTALX2.0.12, Multiple Alignment using Fast Fourier Transform (MAFFT) v6.717b, and MUSCLEv3.8 algorithms. The phylogenetic tree programs Figtree v1.3.1 and Quicktree-SD were used in addition to subgroup the cloned genomic and cDNA sequences.

Accession numbers for genomic and cDNA sequences deposited at the GenBank database are as follows: TaSK1A [GenBank: JX307288, GenBank:JX294419], TaSK1B [GenBank:JX307289, GenBank:JX294420], TaSK1C [GenBank:JX307290, GenBank:JX307292], TaSK2A [GenBank:JX307291, GenBank:JX307293], TaSK2B [GenBank:JX312689, GenBank:JX312688], TaSK2C [GenBank: JX312691, GenBank:JX312690].

Chromosome localization of TaSKs

Nullisomic-tetrasomic lines (cv Chinese Spring) originally established by Sears et al. (1966) [33] were provided by the National small grains Germplasm Research facility of the United States Department of Agriculture (USDA).

Genomic DNA was extracted using the standard CTAB-DNA isolation [72].

Specific primers for amplification of TaSK1-A, TaSK1-B and TaSK1-C sequences were respectively SF97/SR96, SF61/SR98 and SF99/SR101 (Table 2). Specific primers for amplification of TaSK2-A, TaSK2-B and TaSK2-C were respectively SF104/SR105, SF102/SR102, dCAPS-T2C-F/dCAPS-T2C-R (Table 2). TaSK2-C amplicons were subsequently digested with RsaI endonuclease.

DNAs extracted from lines of cv Sonora and Bobwhite were used as control as the primers were designed based on cv. Sonora TaSK sequences.

In vitrokinase activity assays

The N-Terminus of TaSK1, TaSK2, BIN2, OsGSK7 and TaGSK1 full length proteins were cloned in frame with a Gluthatione-S-Transferase (GST) tag.

GST-TaSK1, GST-TaSK2, GST-BIN2, GST-OsGSK7, and GST-TaGSK1 fusion proteins were overexpressed in E. Coli. and affinity purified in native conditions on Gluthatione Sepharose 4B resin. In vitro kinase reactions were performed by adding to 10 μl of the purified GST-fusion protein, 20 μl of the kinase activity buffer (20 mM HEPES, pH 7.4, 15 mM MgCl_2, 5 mM EGTA, 1 mM DTT) containing ATP γ³²P and the bovine myelin basic protein (MBP, fragment) as described by Jonak et al., (2000) [27]. The reaction was incubated at room temperature for 45 minutes and subsequently stopped by adding 10 μl of SDS-Page loading buffer. After denaturation at 95°C for 1 minute, protein phosphorylation was analysed by autoradiography after migration on a 12% SDS/PAGE gel.

Phylogenetic analysis

After the initial selection of GSK homologs by annotation mining and BLAST searches, homologs were detected using BLAST (cutoff 30% alignment identity and 80 amino acids alignment length) in selected genomes. The GSK sequences of selected lineages were obtained from EnsemblPlants, The Arabidopsis Information Resources (TAIR), Plant Genome DataBase (PlantGDB), Phytozome, GenBank, Hawaii Papaya Genome Project ASGPB, and the Rice genome annotation project (RGAP) databases. The current dataset contains only genes confirmed from completely sequenced genomes, except for Hordeum vulgare and Triticum aestivum. Identified sequences were analyzed for the presence at the protein or predicted protein level of the relevant GSK motifs and residues. Only Poaceae sequences with evidence at transcript level were included in this study. Curation of initial phylogenetic trees led to discarding duplicated sequences. Sequences not belonging to annotated shaggy kinases were separated by a long branch in the initial trees; all but two of the sequences were discarded, the remaining serving as outgroup. Protein accessions or locus name are listed in Additional file 1.

The multiple alignment of the full length amino acid sequences was generated using M.A.F.F.T. v6.717b [32] in the ‘auto’ mode and was subsequently manually curated (removing regions of poor alignment quality) using Jalview v2.7 [73], resulting in 501 columns that were used for phylogenetic inference. Two Zea mays [EnsemblPlants:GRMZM2G075992_P01, GRMZM2G332798_P01] and one Carica papaya [ASGPB: CARPA_18.208] sequences were not included in the phylogenetic analyses because the gene models appeared incomplete, covering less than 50% of the alignment, and the sequences thus were placed on very long branches in the initial phylogenies. Final topologies were inferred by Neighbour-Joining (NJ) using QuicktreeSD [30, 74] with 1,000 bootstrap samples, by Maximum Likelihood (ML) using TreePuzzle v5.2 [75] with quartet puzzling and eight gamma distributed rates and by Bayesian Inference (BI) using MrBayes v3.1.2 with eight gamma distributed rates and two hot / two cold chains for two million generations. ProtTest v1.3 [76] was used to determine the model best suited to the dataset and turned out to be JTT with gamma distributed rates and invariant sites, which was hence applied for ML and BI inference. The trees were outgroup-rooted on the branch leading to two related, non-shaggy type A. thaliana kinases [TAIR: AT1G73690.1, AT1G67580.1].