Abstract

Currently, Dual Specificity YAK1-Related Kinases (MNB/DYRK) were found in slime molds, protista, fungi, and animals, but the existence of plant homologues is still unclear. In the present study, we have identified 14 potential plant homologues with the previously unknown functions, based on the strong sequence similarity. The results of bioinformatics analysis revealed their correspondence to DYRK1A, DYRK1B, DYRK3, and DYRK4. For two plant homologues of animal DYRK1A from Physcomitrella patens and Arabidopsis thaliana spatial structures of catalytic domains were predicted, as well as their complexes with ADP and selective inhibitor d15. Comparative analysis of 3D-structures of the human DYRK1A and plant homologues, their complexes with the specific inhibitors, and results of molecular dynamics confirm their structural and functional similarity with high probability. Preliminary data indicate the presence of potential MNB/DYRK specific phosphorylation sites in such proteins associated with plant cytoskeleton as plant microtubule-associated proteins WVD2 and WDL1, and FH5 and SCAR2 involved in the organization and polarity of the actin cytoskeleton and some kinesin-like microtubule motor proteins.

1. Introduction

From the time of Yak1p kinase discovery in budding yeast more than two decades ago, DYRK (MNB/DYRK) family members have been identified in all eukaryotes. At present, different members of this family are regarded as key players in a number of cellular processes [1, 2]. The analysis of the genomic structure and the conservation degree within the kinase domain of the mammalian DYRK subfamily members (DYRK1A, DYRK1B (Mirk), DYRK2, DYRK3 (REDK), and DYRK4) reveal that all these kinases originated via gene duplication during the late periods in metazoan evolution [3, 4]. Mammalian Dual Specificity YAK1-Related Kinases (DYRKs) comprise a group of tyrosine-regulated kinases within the CMGC protein kinase family (MNB/DYRK subfamily) [5], homological yeast Yak1 [1, 6, 7], and Drosophila minibrain kinases [7]. Animal DYRK kinases participate in several signaling pathways involved in the development and support of cell homeostasis [3, 8]. One main feature of DYRKs is their dual ability to autophosphorylate tyrosine residues and to phosphorylate their substrates on serine and threonine residues within an RPx(S/T)P consensus sequence [9, 10]. At the same time, like many other protein kinases, mammalian DYRK1A prefers serine over threonine in substrate peptides and does not detectably phosphorylate tyrosine [10]. It was established that once Y-autophosphorylated, DYRKs lose their tyrosine kinase activity and further function as serine/threonine kinases [2].

DYRK family members phosphorylate many substrates and play a key role in the signaling pathways regulating cell cycle, nuclear functions in cell proliferation, differentiation, and brain development [3, 8, 11]. It is well known that phosphorylation of tau proteins (TAU_HUMAN, P10636: 523-PGSRSRTPSLPTP-535) by DYRK1A is involved in the brain development [1214], and its gene mutations and overexpression are strongly associated with cytoskeleton malfunction and Down syndrome [5, 15]. However, many aspects of MNB/DYRK kinases role remain unclear. The existing data regarding the role of MNB/DYRK members in the regulation of the dynamics of cytoskeleton in general and microtubule (MT) in particular are not complete, as their ability to phosphorylate tubulin in animals and fungi remains unknown [16, 17]. At the same time, Murakami et al. (2012) have reported that DYRK1A can efficiently phosphorylate both MT-associated proteins MAP1A and MAP2 within MTs [18]. Mass spectrometric analysis revealed that both mammalian MAP2 and MAP1A contain the consensus sequence for the substrate specificity of DYRK1A-R(x)xxS/T(P/V) (RPxSP and RxSP for MAP2 and MAP1A) [10, 1820]. MAP1A and MAP2 phosphorylation by DYRK1A was also confirmed in vitro by coimmunoprecipitation [18]. Thus, such modification of microtubule-associated proteins is the only known pathway for animal DYRKs to modulate MT dynamics [18].

The recent data suggest that the DYRK family members are found in 4 of 5 main taxa of slime molds, protista, fungi, and animals, though the number of these protein kinases genes varies considerably at different kingdoms. Thus, Yak1—the first known DYRK—turned out to be the only member of MNB/DYRK family in Saccharomyces cerevisiae. However, all known DYRK proteins share common structural, biochemical, and functional properties with their yeast ancestors [11, 22].

Several years ago the evidence for the existence of DYRK-like kinases in plants appeared [3, 23]. It is believed that almost all plants have one to three YAK1-related genes, but none of the DYRK genes in plants have been characterized yet. Also, plants do not have apparent HIPKs (Homeodomain-Interacting Protein Kinases) that regulate transcription by interacting with homeodomain proteins in animal cell [24] but do have expanded PRP4 kinase subfamilies [23]. The functions of these genes also remain unclear, as well as the existence of DYRK1 homologues in plants. Therefore, the plant MNB/DYRK kinases and their putative functional role are to be identified.

Here we report that a few plant proteins with the previously unknown functions have the strong sequential and structural homology to mammalian DYRK and yeast YAK1 protein kinases. It is also supposed that MNB/DYRK specific phosphorylation sites might be present in some plant proteins associated with cytoskeleton as well.

2. Methodology

2.1. Database Search and Sequence Analysis
2.1.1. Database Search

The reviewed (annotated) sequences of DYRK kinases were taken from the Protein Knowledgebase UniProtKB http://www.uniprot.org/ [25]. The search for plant DYRK homologues was performed via BLASTp scanning of the UniProtKB (Swiss_Prot and TrEMBL) database within “Viridiplantae” group using human DYRK1A (Q13627) catalytic domain as reference sequence.

2.1.2. Sequence Alignments

The multiple alignments of the amino acid sequences were constructed using ClustalX (v. 2.0.5) program (http://www.clustal.org/) [26]. In turn, the amino acid sequences were aligned using a set of BLOSSUM substitution matrices [26].

2.1.3. Domain Architecture Analyses

The domain architecture of proteins was analyzed using the SMART tool http://smart.embl-heidelberg.de/ [27, 28] based on the ExPASy Proteomics Server information (http://www.expasy.org/).

2.1.4. Protein-Protein Interactions

The potential function of plant homologs was also analyzed by protein-protein interactions modeling against the signaling cascades of Homo sapiens using the STRING 9.05 tool (http://string-db.org/) [29].

2.1.5. Identification of Consensus Site for DYRK1A Phosphorylation

For identification of consensus site (motif) for MNB/DYRK phosphorylation a query pattern against a UniProtKB database were searched. R-x(1,3)-[ST]-[PVL], R-P-x-S-P and R-x-S-P patterns were constructed based on published data [10, 1820]. For database screening we used Pattern Search tool from Protein Information Resource (PIR—http://pir.georgetown.edu/), PROSITE (http://prosite.expasy.org/) and BLASTp (http://blast.ncbi.nlm.nih.gov/ for RPxSP).

2.1.6. Cladistic Analysis

The cladistic analysis was performed based on the catalytic kinase domains multiple alignments [30, 31] using the Neighbor-Joining method [32, 33]. The dendrograms were visualized and analyzed with the help of the MEGA4 (http://www.megasoftware.net/) program [34].

2.2. Spatial Structure Modeling, Docking, and Molecular Dynamics Simulations
2.2.1. D-Structures Visualisation, Alignment, and Analysis

PyMol 1.4 package was used for molecular visualization, structural alignment and analysis (http://www.pymol.org/).

2.2.2. Molecular Modeling Studies: Template Selection and Homology Modeling

Amino acid sequences of proteins from Physcomitrella patens and Arabidopsis thaliana were obtained from UniProt knowledge base (A9U0E5 and Q8RWH3) [25]. Homology protein structure modeling was performed in Modeller 9v8 (http://salilab.org/modeller/) [35] based on template X-ray RCSB Protein Data Bank structures: 3ANQ (2.6 Å) [36], 2VX3 (2.4 Å) [37, 38], and 1Z57 (1.7 Å) [39] specified with PDB-BLAST (http://www.rcsb.org/pdb/search/advSearch.do?st=SequenceQuery) [39, 40].

Several methods for helix prediction were applied in a direct manner using their websites to allocate helix segments of DYK1A including LOBO (http://protein.bio.unipd.it/lobo/) [41], ArchPred (http://manaslu.aecom.yu.edu/loopred/) [42], and SooperLooper [43]. The accuracy and stereochemical quality of the models were assessed by the PROCHECK program (http://www.ebi.ac.uk/thornton-srv/software/PROCHECK/) [44] and online server MolProbity (http://molprobity.biochem.duke.edu/) [45]. A profile of the resulting model was checked using Verify3D (Structure Analysis and Verification Server v. 4 (SAVES v4)—http://services.mbi.ucla.edu/SAVES/) [46].

All root mean square deviations (RMSDs) of the structures were calculated using MODELLER (http://salilab.org/modeller/) [35] and PyMol functions.

For ligand alignment, protein active-sites comparison and ligand binding investigation we used LigAlign plaguing integrated in PyMol package (http://compbio.cs.toronto.edu/ligalign/) [47]. The role of phosphate binding pocket conservative residues was determined by alignment of the several Protein Data Bank (PDB ID: 2VX3, 3KVW, 3ANQ, 1Q8Y) DYRK1A structures in complexes with EHB ((1z)-1-(3-ethyl-5-hydroxy-1,3-benzothiazol-2(3H)-ylidene)propan-2-one) [36], d15 (N-(5-{[(2S)-4-amino-2-(3-chlorophenyl)butanoyl]amino}-1H-indazol-3-yl)benzamide) [48], IRB ((2Z,3E)-7′-bromo-3-(hydroxyimino)-2′-oxo-1,1′,2′,3-tetrahydro-2,3′-biindole-5-carboxylic acid) [38, 49] and ADP (adenosine 5′-diphosphate) [50] was performed.

2.2.3. Molecular Docking

Docking procedures were done in CCDC Gold 5.1 program (http://www.ccdc.cam.ac.uk/) [51]. 12 Å (the radius of the active site), 100 (the number of genetic algorithm (GA), GoldScore and Astex Statistical Potential (ASP) (scoring and rescoring functions resp.), 4.0 and 2.5 (Van der Waals forces (VdW) and hydrogen bonding (h-bonding) were set as the docking protocol parameters; the parameters for genetic algorithm were set automatically. The additional function implemented during the docking search was flexibility of amino acids. Grounding on the results of the pharmacophore building we assigned flexibility to the list of flexible side chains.

The ligand models were built using the Marvin Sketch (http://www.chemaxon.com/) editor. Topology and parameter files were generated for the ligands using the SwissParam web-server (http://www.swissparam.ch/) [52].

2.2.4. Analysis of Ligands Pharmacophore Properties

Ligand-enzyme interactions were investigated in the LigandScout (v. 3.0) program (http://www.inteligand.com/) [53].

2.2.5. Molecular Dynamic Simulations

To obtain a relaxed protein structure we provided a short molecular dynamic simulation (1 ns) using a G53a6 force field with the GROMACS 4.5.3 package (http://www.gromacs.org/) [54] without the ligand. All free MD simulations were carried out in GROMACS (in vacuum and then in the water) with Charmm27 (http://www.charmm.org/) force field. All computations were performed in a high performance IFBG cluster (VO CSLabGrid, Ukrainian National Grid—http://ung.in.ua/) [54, 55].

Each system was inserted in a water box where the layer of water molecules was equal to 10 Å. The systems were subjected to steepest descent energy minimization for 20,000 steps. Then protein backbone was frozen, and solvent molecules with counter ions were allowed to move during a 100 ps position restrained MD run.

All simulations were run under periodic boundary conditions with constant NPT ensemble (number of particles (N)/system pressure (P)/temperature (T), http://www.gromacs.org/) by using Berendsen’s coupling algorithm for the constant temperature (310 K) and the pressure (1 bar) maintainance. The time step for the simulations was 2 fs. The electrostatic interactions were calculated using the Particle-mesh Ewald (PME) algorithm with an interpolation order of 4 and a grid spacing of 0.12 nm. Van der Waals forces were studied by using a cutoff of 10 Å and the coordinates were stored every 100 ps. Free dynamics simulation was performed at 310 K for 26 ns for each complex. Ligand-binding site interactions were evaluated based on average RMSD and short range (SR) Coulomb interaction energy.

3. Results

3.1. The Sequence Homology Search of Plant DYRK-Like Kinases

A UniProtKB http://www.uniprot.org/ database search, revealed presence of 18 “Reviewed” sequences of animal DYRKs from Homo sapiens (DYRK1A, DYRK1B, DYRK2, DYRK3, DYRK4), Macaca fascicularis (DYRK3), Mus musculus (DYRK1A, DYRK1B, DYRK2, DYRK3, DYRK4), Rattus norvegicus (DYRK1A, DYRK3), Xenopus laevis (DYRK1A), X. tropicalis (DYRK1A), Gallus gallus (DYRK2), Drosophila melanogaster (DYRK2 (smi35A), DYRK3), and 2 (DYRK1, DYRK2) from Dictyostelium discoideum (Mycetozoa).

In order to find potential plant DYRK1A homologues, we implemented SIB BLAST-search within “Viridiplantae” group, using human DYRK1A (UniProt: Q13627) catalytic domain as reference sequence. The primary group of plant homologues comes to 90 potential protein kinases having entire catalytic domain (assigned based on analysis in SMART—http://smart.embl-heidelberg.de/). In compare with human DYRK1A, the group exhibit 35–51% of sequence identity, 55–70% of sequence similarity and a gap percentage of 1–10%. All plant sequences were represented by proteins of unknown function, and had “Unreviewed” UniProt status. The Neighbor-Joining clustering of plant homologues and “Reviewed” DYRK catalytic domains from Animalia, Fungi, Protista and slime molds resulted in common clade. In the N-J tree, catalytic domains of   “Reviewed” DYRK protein kinases formed common clade with 34 potential plant homologs. Out of them, only 14 plant proteins have been deposited (August 2014) in UniProt knowledgebase as the full-length sequences (Table 1, Figure 1).

In animals, DYRK members, such as DYRK1A, DYRK1B, DYRK2, DYRK3, and DYRK4, in spite of sequence similarity of their catalytic domains, in general, exhibit differences in their complete sequences, functionally important motives, 3D-structure and substrate proteins specialization [1, 5, 8]. They have various partner proteins and specific differences in the range of protein-protein interactions. We apply “search by protein sequence” of “protein-mode” uses quantitative sequence similarity search and often distributes a given interaction fractionally among several protein pairs of the target organism. Thus, the STRING v.9 (http://string-db.org/) search with complete sequences request and using scan against a well-studied human protein-protein interactions database (selecting “organism: Homo sapiens”) allowed us to perform in silico analysis of plant proteins with unknown function and to confirm them as potential Dual Specificity YAK1-Related Kinases.

Thus, plant homologues, deposited in the UniProt knowledgebase as a group of proteins of unknown function, have been identified bioinformatically, and DYRK1A homologs from P. patens subsp. patens (UniProt: A9U0E5) and A. thaliana (UniProt: Q8RWH3), were selected for further spatial structure and function prediction. The choice of these proteins was not only due to their high sequence identity with mammalian DYRK1A, but also by the fact that P. patens (NCBI-Taxonomy ID: 3218) and A. thaliana (NCBI-Taxonomy ID: 3702) are model organisms and intensively investigated in the framework of genomic and proteomic projects.

3.2. Plant DYRK1A Homologues Spatial Structure Modeling and Ligand-Binding Sites Prediction

Spatial structure modeling was performed for two proteins: A9U0E5 from P. patens and Q8RWH3 from A. thaliana. Suitable template was selected through a BLAST search of the RCSB Protein Data Bank in order to develop reliable models, and human DYRK1A (PDB: 3ANQ) was chosen as the best template for A9U0E5 and Q8RWH3 3D-structure modeling. After the phases of energy minimization in vacuum and in water, Tip3p model, using Charmm27 force field (Figure 2), the quality of geometry of human DYRK1A and two plant models were confirmed with PROCHECK and MolProbity servers. Finally, the models were checked with structure evaluation server Verify3D. Structural superimposition of 3D-models of catalytic domains from P. patens (UniProt: A9U0E5) and A. thaliana (UniProt: Q8RWH3) plant homologues with template 3ANQ and another (control) PDB structure 2VX3 (chain A) of human DYRK1A demonstrated sufficient global and local structural similarity (RMSD < 1) (Figure 3).

Besides the global similarity, the identity of some functionally important motifs and residues was demonstrated by structural superimposition. It is known that DYRK family members are distinguished from other homologous kinases by their activation mechanism [2, 56], since DYRK kinase active state depends on tyrosine autophosphorylation at a conserved YxY site (Y319xY321 in human DYRK1A) in the catalytic domain [2, 3]. This activation domain located directly in front of Subdomain VIII was shown to be intramolecularly phosphorylated only during translation, while the mature (Y-phosphorylated) DYRK family members have only serine/threonine kinase activity [21, 57]. Therefore, Tyr321 (Y321) phosphorylation (pY321) is necessary for DYRK1A (PDB: 2VX3) activation [4, 38]. Structural superimposition demonstrates the identity of site (YxxxS-motif) in the activation loops of human DYRKA1 and plant homologues, and is important for DYRKs activation through tyrosine residues (Y321 in human DYRK1A, Y233 in Physcomitrella and Y284 in Arabidopsis homologues) autophosphorylation [3] (see Figures 3 and 4). It should be noted that we identified YxY-motif in activation loops of all studied plant homologues except Q019C0 (Ot05g02500) from O. tauri containing sFtiqSr motif in position identical to that in canonical yeast Yak1-like protein kinase (acYak1, Q5A3L7_CANAL) from Candida albicans (Figure 4).

Also, we have identified some other functionally important canonical MNB/DYRK residues, highly specific for mammalian DYRK, that never occur in other CMGC group members: Cys289, Cys312, Tyr319, and Gln323 [58]. It is established that Gln323 in human DYRK1A activation loop is involved into the unique activation mechanism of DYRK-related kinases [56]. Moreover, some conserved Gln residues are identical in all studied plant homologues (Figures 4 and 5). Their similar positions in spatial structure were Q235 in A9U0E5 from P. patens and Q286 in Q8RWH3 from A. thaliana (Figures 4 and 5). Tyr321 are highly conserved in all known protein kinases of the DYRK family, including those from yeasts, slime molds, protozoans, and plant homologues (except Q019C0_OSTTA from O. tauri and acYak1 from C. albicans). These residues are extremely important, as far as it is known that Gln323 and Tyr321 mutations critically reduce catalytic activity of the mammalian DYRK1A [56].

Besides that, it is known that first residue (glutamine) of helixG is a highly conserved in DYRK family and rarely found in other kinases [56]. In human DYRK1A Glu366 (E366) in helixG was identical to Glu329 (E329) in A. thaliana homolog (Figure 4). In A9U0E5 from moss P. patens in this position was cysteine Cys278 (D278), identical to the first residue in helixG of lmDYRK from L. mansonii (Figure 4) [56].

As it was already mentioned, Tyr321 (Y321) phosphorylation (pY321) is necessary for DYRK1A (PDB: 2VX3) activation [4, 11, 38]. In order to justify inferences from sequence similarity on phosphorylated (active) and non-phosphorylated (inactive) states of animal DYRK and plant homologues, a structural comparison of activation loops regions was performed. The models of nonphosphorylated (Y321 and corresponding residues) structures with “active” (pY321) DYRK1A (2VX3) and DYRK2 (3K2L and 3KVW) phosphorylated structures from Protein Data Bank were compared. Nonphosphorylated models have undergone the stage of molecular dynamics simulation, and then, YxYxxS regions of PDB structures of activated human DYRKs (Figure 6, see 1, 2, 3) and models of inactive human DYRK1A (Figure 6, see 4) were superimposed with reconstructed plant homologues from P. patens (Figure 6, see 5) and A. thaliana (Figure 6, see 6). The conformation of plant tyrosine residues which fit human DYRK1A Y321 was similar to that in inactive state of human DYRK1A. Therefore, we assume that the initial self-activation of plant homologues probably occurs in the same manner by autophosphorylation of conserved tyrosine residue corresponding Y321 of human DYRK1A. It should be noted that canonical YxxxS motif has been identified in activation loops of all studied plant homologues except Q019C0 (Ot05g02500) from O. tauri.

Since most common experimental evidence for a role of protein kinases is based on the use of various inhibitors, the next logical step was to analyze the possible mechanisms of binding of known specific inhibitors of animal DYRK1A with appropriate plant homologues [59, 60]. All protein kinases catalyze the same type of reaction and inhibitory agents interact with the relatively conserved ATP-binding pocket [5961]. Though, numerous variable and sometime quite unique features surrounding this site in the diverse kinases [62] give rise to the specific pockets and/or docking sites for the putative ATP-competitive inhibitors. Thus, the majority of the existing inhibitors operate in ATP-competitive manner and their selectivity depends on amino acid composition and structural similarities/differences in their ATP-binding pocket and its environment [59]. Crystallography revealed several ATP-competitive inhibitors binding mechanisms for animal DYRK1 [17, 36, 48, 63, 64]. Hence, the alignment of the several Protein Data Bank 2VX3, 3KVW, 3ANQ, 1Q8Y DYRK1A structures in complexes with EHB, d15, IRB and ADP was performed to study the concerning functional (catalytic) similarity of human DYRK1A and its plant homologues, and to understanding the role of phosphate binding pocket conservative residues.

Subsequent alignment of ligand-protein complexes identified some common features in heterocyclic core of studied compounds (Figure 7). Such procedure allowed to determine reactive atoms in the ligands for further experiments, and, in particular, determined the choice of ADP and d15 as candidate structures for pharmacophore specification and docking. Respectively, a spatial structure models for human DYRK1A with ADP and d15 buried in binding site were created (Figure 8).

As a result, we have identified a list of amino acid residues, potentially critical for ligand binding site formation in human DYRK1A: Gly166, Gly168, Phe170, Gly171, Gln172, Val173, Ala186, Val222, Phe238, Glu239, Leu241, Asn244, Leu245, Lys289, Glu291, Asn292, Leu294, and Asp307. On the basis of the comparison of human DYRKA1 3D-models structural alignment with its A. thaliana and P. patens homologues, the identity of amino acids forming the surface of the binding site was confirmed. In turn, the ligand position analysis assuming the comparison of some interaction features (h-bonding, hydrophobic properties, center of negative/positive charge) allowed us to discard EHB and IRB from further studies.

Binding sites comparison in yeast (S. cerevisiae) SR protein kinase Sky1p (PDB: 1Q8Y) and human DYRK1A (PDB: 2VX3) allowed us to definitively identify residues that form phosphate binding cavity and those directly involved in ADP (Lys187, Phe169, Ser170, Glu247, Leu249, Asn252, Glu298, and Asp550) and d15 (Glu239, Leu241, and Asn292) binding (Figure 8).

Moreover, the comparison of ligands positions confirms d15 from 2VX3 PDB-structure as the best object for subsequent docking and MD operations of human DYRKA1 and two plant homologues.

3.3. Reconstruction and Comparison of Human DYRK1A and Plant Homologues Complexes with Specific ATP-Competitive Inhibitor d15

As it was shown earlier, Y321 tyrosine residue in YxxxS motif is extremely important for DYRK1A activation. Taking into account this fact, we have reconstructed similar conserved pY-residues in Physcomitrella and Arabidopsis homologues using the Marvin Sketch editor. d15 molecular docking was performed using CCDC Gold genetic algorithm. As a result, models of complexes with d15 were built for human DYRKA1, Q8RWH3 from A. thaliana and A9U0E5 from P. patens. All molecular docking and dynamics operations were performed in conformity with protocols described in Section 2.2.

Subsequent MD calculations demonstrate that average RMSD of ligand and SR Coulomb interaction energy between binding site and d15 were around 0.1 nm and −221.4 kJ/mol, respectively. The slight differences for average number of H-bonds between active site and the substrate proteins were observed and exhibit 4.1 for human and plant complexes. This confirms structural similarity of ligand binding in human DYRK1A and plant homologues and give evidence for further investigation of their catalytic activity (Figures 9(a), 9(b), and 9 (c)). However, it should be noted that after 2 ns of MD it appears that ligand binding state is not so stable and it is evident from the RMSD graphics, especially for protein from A. thaliana. However, the values of both energy and RMSD do not exceed estimated scores for protein from P. patens (Figure 10).

At the same time, molecular docking and MD calculations confirm our assumption concerning structural similarity of human Dyrk1A and plant homologues from A. thaliana (UniProt: Q8RWH3) and P. patens (UniProt: A9U0E5).

3.4. Presence of MNB/DYRK Specific Phosphorylation Sites in the Proteins Associated with Plant Cytoskeleton

The profile search for MNB/DYRK specific phosphorylation site motifs (see Section 2.1.5) in A. thaliana proteome has revealed 370 (67 UniProt Reviewed) potential target proteins. It is interesting that among these hits we identified were proteins associated with plant cytoskeleton. In particular, MNB/DYRK specific phosphorylation sites were identified in plant MAPs-WVD2 (163-PLTRPKSPKLN-173) and WDL1 (202-PLTRPKSPKLI-212) that modulate helical organ growth and anisotropic cell expansion in Arabidopsis [65]. WDL1 (At3g04630, Protein WVD2-like 1, Q8GYX9) and WVD2 (At5g28646, protein WAVE-DAMPENED 2, Q84ZT9) are known as MAPs that regulate orientation of interphase cortical microtubules. Consensus MNB/DYRK specific phosphorylation site motifs were identified in FH5 (163-PPTRPKSPPPR-173) and SCAR2 (1298-RLPRPRSPLVD-1309)-A. thaliana proteins involved in actin cytoskeleton organization and polarity.

In addition, potential MNB/DYRK sites of specific phosphorylation were identified in four “Unreviewed” proteins, presumably from the kinesin superfamily: At5g06670 (Kinesin heavy chain-like protein, UniProt: Q9FG03; 71-VRFRPLSPREI-81), At5g42490 (ATP binding microtubule motor family protein, UniProt: Q9FIG8; 835-MEIRPESPADS-845), At3g12020 (Kinesin heavy chain-like protein; UniProt: Q9LHL9; 74-VRFRPLSPREI-84) and At4g39050 (Kinesin like protein F19H22.150, UniProt: Q9SVI8; 9-RSSRPPSPASS-19).

4. Discussion

Our primary in silico search for the putative plant MT- and cell cycle-related serine-threonine kinases [6668] was the background for the current study focused on the comparative analysis of the specific sequential and structural features of the selected plant homologues with known DYRK family members. It is proved that mammalian DYRK isoforms regulates cytoskeletal dynamics, and DYRK1A phosphorylates human microtubule-associated proteins tau (at Thr212) and MAP1B (at Ser1392) [10, 18–20]. The proteins, evolutionarily close to DYRK kinases, were revealed in plant cell as the results of NJ clustering of catalytic kinase domains suggest. Phylogenetic mapping of the conserved protein kinase catalytic domains is a basic step of the functional characterization of new family members [69]. After that, we focused on two plant homologues from the model plants: A9U0E5 from P. patens and Q8RWH3 from A. thaliana and found that A9U0E5 demonstrate the higher total weight of the alignment with human DYRK1A. The detailed study of both proteins revealed the strong similarities not only in whole catalytic domain sequences and spatial structures, but also in the mechanisms of kinase domains activation (autophosphorylation) that is recognized as the hallmark of animal and fungi MNB/DYRK kinases [14, 7073]. Sequence comparison of the activation loop (P+1 site) and helixG regions of the known DYRK family kinases and the potential plant homologues indicates strong similarity of these important regions with the mammalian DYRK, dmMNB from D. melanogaster, YAK1 from S. cerevisiae, ddYakA from D. discoideum, and some other known DYRK members.

Structural superposition of the complete atomic models of A9U0E5 from P. patens, Q8RWH3 from A. thaliana and DYRK1A from H. sapiens demonstrates identity of YxY motif (YxxxS) topology with the activation loops of human DYRK1A and plant homologues as well as with the important for DYRKs activation tyrosine residues corresponding to Y321 in human DYRK1A (See Figures 4 and 6) [3]. At the same time, total RMSD < 1 at superimposed plant models and experimentally confirmed PDB-structure 2VX3 of human DYRKA1 indicate close spatial structure relationship of their catalytic domains [40]. It should be mentioned that all experimentally proven structures of MNB/DYRK members deposited in RCSB Protein Data Bank, initially contain the second phosphorylated tyrosine in YxpxxS motif of activation loop. Undoubtedly, the study of conformational changes after Y321 modification is still far from its completion. The detailed up-to-date research of the molecular dynamics reckoning with the topology of the modified residues would be at the focus of our further experiments.

One of the key instruments for the protein kinase functional studies remains the selective inhibitors of individual or groups of kinases [59]. Now, most DYRK1A inhibitors are ATP-competitive [38]. In case of animal DYRK1A, the most studied inhibitor is indazol compound d15 [74]. Molecular docking, MD calculations, and comparison of human and plant complexes with d15 demonstrate similarity of amino acid composition in binding site and in mechanisms of ligand binding by DYRK1A and studied plant homologues from A. thaliana and P. patens. Thus, it not only once again demonstrates similarity of DYRK1A and plant homologues, but also confirms the legitimacy of d15 application for further in vitro research, including the study of the role of plant DYRK homologues in phosphorylation of cytoskeletal proteins and regulation of cell division.

However, the role of plant MNB/DYRK homologues in cytoskeleton regulation remains unclear. Plant homologues for mammalian MAP2 have not been identified earlier, and their divergence probably occurred at early stages of evolution [75]. Existence of plant homologues for MAP1A is also controversial. Blast scanning (against human MAP1A) of Arabidopsis proteome has revealed some homologues with low similarity (identity ≤ 21). In turn, the At2g22795 was the most similar: Identities = 21%, Positives = 43%; Gaps = 3% [75], but its interplay with the cytoskeleton is questionable. Furthermore, none of the revealed proteins contained RPxSP motif that was found only in other uncharacterized hypothetical protein, At2g12875 (Identities = 19%, Positives = 38%, Gaps = 8%) with the unknown relation with cytoskeleton. These motifs were also absent in the previously discovered plant MAPs from RP/EB (EB1) family and in particular, EB1A (Q7XJ60), EB1B (Q9FJJ5), and EB1C (Q9FGQ6) from A. thaliana [7577]. Taken altogether, our data suggest that the putative plant DYRKs input in microtubule dynamics regulation in higher plants via MAP1A, MAP2, and EB1 phosphorylation (A, B, and C) is unlikely.

Nevertheless, the search for MNB/DYRK specific phosphorylation sites in proteome of A. thaliana has revealed 370 (67 UniProt Reviewed) potential hits [65]. It has to be noted that among these hits we also identify plant proteins associated with microtubules. In particular, MNB/DYRK specific phosphorylation siteswere identified in plant MAPs-WVD2 (163-PLTRPKSPKLN-173) and WDL1 (202-PLTRPKSPKLI-212) modulating helical organ growth and anisotropic cell expansion in Arabidopsis [65]. WDL1 (At3g04630, Protein WVD2-like 1, Q8GYX9) and WVD2 (At5g28646, Protein WAVE-DAMPENED 2, Q84ZT9) are known as MAPs that regulate interphase cortical microtubules orientation. WVD2 and WDL1 modulate both rotational polarity and anisotropic cell expansion during organ growth and promote clockwise root and etiolated hypocotyls coiling, clockwise leaf curling, and left-handed petiole twisting [65, 78].

Consensus MNB/DYRK specific phosphorylation site motifs were identified in A. thaliana FH5 and SCAR2 proteins involved in the organization and polarity of actin cytoskeleton. It is known that FH5 (At5g54650, Formin-like protein 5, Q94B77) interacts with the barbed end of actin filaments and nucleates actin-filament polymerization in vitro. FH5 seems to participate in cytokinesis and its expression occurs in endosperm and localizes in a cell plate, a plant-specific membranous component that is assembled at the plane of cell division [79]. In turn, SCAR2 (At2g38440, protein SCAR2, Q5XPJ9) involved in the regulation of actin and microtubule organization. It is known that SCAR2 is a part of a WAVE complex that activates the Arp2/3 complex [80, 81] able to regulate trichome branch positioning and expansion in A. thaliana [10, 31, 82]. In addition, potential MNB/DYRK sites of specific phosphorylation were identified in four “Unreviewed” proteins, presumably belonging to the kinesin superfamily [83]: At5g06670, At5g42490, At3g12020, and AT4g39050.

Therefore, plant protein kinases MNB/DYRK as well as their mammalian counterparts input into cytoskeleton regulation. However, the function of circa a half of A. thaliana proteins with MNB/DYRK phosphorylation site is still unknown, though in the nearest future new data about MNB/DYRK molecular targets and new players in the regulation of plant cytoskeleton might be identified.

5. Conclusions

Summing up the study, our results reveal initial assumption about the existence of MNB/DYRK protein kinase family (CMGC group) in higher plants. In particular, it is shown by strong global sequence and structure similarity in catalytic domains with human DYRK1A, and unique DYRK family activation loop (P+1 site) and helixG motifs. Also, superimposition of the active site residues of human DYRK1A and plant homologues from P. patens and A. thaliana, as well as the results of docking and MD-simulations confirmed considerable similarity in their ATP-biding pockets. This allows the proceeding of subsequent experimental confirmation of plant homologues using known mammalian DYRK1A inhibitors. In particular, it confirms possibility of application of mammalian DYRK1A ATP-competitive inhibitor d15, for further experimental study of plant DYRK kinases. Also, based on results of profile search of DYRK-specific phosphorylation site motifs, we presuppose that possible involvement of plant DYRK kinases regulate phosphorylation of some proteins that are involved in the organization and polarity of the MT- and actin-cytoskeleton in A. thaliana.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work was supported by the Ukrainian National target scientific and technical program “The development and implementation of grid-technologies in 2009–2013”— http://ung.in.ua/. Molecular dynamics simulations were carried out using computing facilities of the Ukrainian National Grid (UNG)—http://infrastructure.kiev.ua/. The authors also appreciate Dr. Yuliya Krasylenko suggestions for the improvement of the paper’s language.

Supplementary Materials

S1 - Complete FASTA sequences of all proteins used in the study (for details please see Table 1, Figure 1):

(A) 18 "Reviewed" UniProtKB sequences of animal DYRKs from Homo sapiens (DYR1A (Q13627), DYR1B (Q9Y463), DYRK2 (Q92630), DYRK3 (O43781), DYRK4 (Q9NR20)), Macaca fascicularis (DYRK3 (Q4R6S5)), Mus musculus (DYR1A (Q61214), DYR1B (Q9Z188), DYRK2 (Q5U4C9), DYRK3 (Q922Y0), DYRK4 (Q8BI55)), Rattus norvegicus (DYR1A (Q63470), DYRK3 (Q4V8A3)), Xenopus laevis (DYR1A (Q2TAE3)), X. tropicalis (DYR1A (Q0IJ08)), Gallus gallus (DYRK2 (Q5ZIU3)), Drosophila melanogaster (DYRK2 (Q9V3D5, smi35A), DYRK3 (P83102)), and 2 (DYRK1 (Q76NV1), DYRK2 (Q54BC9)) from Dictyostelium discoideum (Mycetozoa).

(B) 14 plant DYRK-like kinases selected based on the Blastp results of UniProt knowledgebase scanning: A9U0E5-Physcomitrella patens subsp. patens, A9RDR5-Physcomitrella patens subsp. patens, A9TA59-Physcomitrella patens subsp. patens, C0PR87-Picea sitchensis, D7MHY5-Arabidopsis lyratalyrata, Q8RWH3-Arabidopsis thaliana, B9SDF7-Ricinus communis, C1MRI0- Micromonas pusilla, Q00ZG2-Ostreococcus tauri, A4RXI5-Ostreococcus lucimarinus, Q019C0-Ostreococcus tauri, D8RDP5-Selaginella moellendorffii, A9TKQ3-Physcomitrella patens subsp. patens, D8RRY9-Selaginella moellendorffii.

S2 - Kinase domain sequences (*.fasta) of all protein kinases used in the study (for details please see Table 1, Figure 1, S1).

S3 - ClustalW2 (http://www.clustal.org/clustal2/) multiple sequence alignment of kinase domains of all protein kinases used in the study (for details please see Table 1, Figure 1, S1).

S4 - Supplement for Figure 1. *.dnd file of Neighbour-Joining tree built based on catalytic domain sequences clustering (S1). (For visualization please use MEGA5: http://www.megasoftware.net/)

S5 - Supplement for Figure 3. PyMOL session file (*.psf): Models of spatial structure of DYRK1A homologues (green) catalytic domain from P. patens (a) and A. thaliana (b) and their structural superimposition with 2VX3 (chain A) Protein Data Bank structure of human DYRKA1 (pink). YxY site - YxYxxS motif in the activation loop, also known as T-loop. (For visualization please use PyMOL: http://www.pymol.org/)

S6 - 3D-model of catalytic domain from Arabidopsis thaliana (UniProtKB: ARATH_Q8RWH3). Please use any program for *.pdb file visualization.

S7 - 3D-model of catalytic domain from Physcomitrella patens (UniProtKB: PHYPA_A9U0E5). Please use any program for *.pdb file visualization.

S8 - Supplement for Figure 7. Structural alignment of DYRK1A inhibitors (EHB, d15, IRB) and ADP. (For visualization please use PyMOL: http://www.pymol.org/)

S9 - (a, b, c). PyMOL session file (*.psf) demonstrating final frames from molecular dynamics simulations for human and plant d15-protein complexes displaying final position of ligand in the ATP-binding site: (a) human DYRKA1; (b) Q8RWH3 from Arabidopsis thaliana; (c) A9U0E5 from Physcomitrella patens. (For visualization please use PyMOL: http://www.pymol.org/)

S10 - Supplement for Figure 9 (d). PyMOL session file (*.psf) demonstrating superimposition of d15 in binding site in human DYRK1A (red) and the similar regions in plant homologues from Arabidopsis thaliana (green) and Physcomitrella patens (light green). (For visualization please use PyMOL: http://www.pymol.org/)

S11 - Supplement for Figures 2 and 10. Microsoft Excel file (*.xls) with results of molecular dynamic simulations carried out during the study: Sheet 1 (Figure 2) - The RMSD of DYRK1A proteins for 26 ns MD simulation: human DYRK1A (1), Q8RWH3 from Arabidopsis thaliana (2), A9U0E5 from Physcomitrella patens (3); Sheet 2 (Figure 10(a)) / Sheet 3 (Figure 10(a)) - The RMSD (Sheet 2) of a ligand and short range columbic interaction energy (Sheet 3) of DYRK1A proteins in complex with d15 during 26 ns MD simulation: human DYRK1A (1), Q8RWH3 from Arabidopsis thaliana (2), and A9U0E5 from Physcomitrella patens (3).

  1. Supplementary Material