Biological characteristics of aging in human acute myeloid leukemia cells: the possible importance of aldehyde dehydrogenase, the cytoskeleton and altered transcriptional regulation

Patients included in the study

Based on the AML cell samples collected at the first time of diagnosis and the predefined genetic criteria we selected 18 low-risk patients with favorable genetic prognosis and 15 high-risk patients with adverse genetic prognosis (see Materials and Methods for group definitions; Tables 1, 2 and Supplementary Tables 1, 2). In our present context the terms high- and low-risk refer to the risk of having a chemoresistant relapse even after the most intensive antileukemic therapy. The two groups differed significantly with regard to age, cause of their leukemia and morphological signs of differentiation. Low-risk patients were generally younger, none of them had secondary AML and half of them had normal karyotype (Tables 1, 2). The frequency of patients with monocytic AML cell differentiation (i.e. FAB M4/5) was higher in the low-risk group. However, the expression of the CD34 stem cell marker did not differ significantly between the two groups. Most high-risk patients were older and had complex cytogenetic abnormalities (i.e. ≥3 abnormalities).

In order to study the impact of morphological signs of differentiation in patient groups that differed significantly with regard to age, we carried out proliferation assays with primary AML cells derived from an alternative cohort of consecutive patients in the presence of hematopoietic growth factors (Supplementary Analysis Tables 1, 2). Our results showed that the proliferative responsiveness of patients with or without morphological signs of differentiation and with age above or below 65 years did not differ significantly. Thus, based on this analysis that included consecutive patients, i.e. not only high-/low-risk karyotypes but also normal and intermediate risk karyotype, we could not find any evidence for a general association between proliferative capacity and patient age/AML cell differentiation.

The 33 AML patient samples were used to perform MS-based proteomics and phosphoproteomics analyses as follows. Firstly, we have carried out an age-dependent analysis in low-risk patients and compared the proteomic and phosphoproteomic profiles between nine elderly low-risk and nine younger low-risk (see Materials and Methods for subgroup definitions). Secondly, we have investigated whether the identified age-dependent differences were also observed when comparing high-risk and all the low-risk patients that differed in age and cytogenetics-related prognosis.

Effects of aging on the AML cell proteome; a proteomic comparison between elderly and younger low-risk patients

We obtained the proteome profiles of AML cells derived from nine elderly low-risk and nine younger low-risk patients using our FASP-based workflow (Figure 1). We quantified 5966 proteins, of which 4369 had a quantitative value in at least five patients in each group (Figure 2A). The proteome profiles from the elderly low-risk and younger low-risk patients were compared using t-test based statistical analysis and resulted in a small set of only 29 differentially expressed proteins, with 18 of them being upregulated and 11 downregulated for the elderly low-risk patients (Figure 2A, Supplementary File 1, Supplementary Table 3 upper part, Supplementary Table 4). Gene ontology (GO) enrichment analyses showed that regulation of T cell mediated immunity (i.e. syntaxin-7, STX7, and galectin-10, CLC) and oxidoreductase activity (i.e. aspartyl/asparaginyl beta-hydroxylase, ASPH, and ribosomal oxygenase 2, RIOX2) were more abundant GO terms in elderly low-risk patients (Figure 2B, top plot). The proteins enriched in this group were primarily located in organelles such as mitochondria and endoplasmic reticulum. Kyoto Encyclopedia for Genes and Genomes (KEGG) pathways analysis demonstrated that the upregulated elderly low-risk proteome was enriched with histidine, ascorbate and aldarate metabolism pathways (i.e. mitochondrial aldehyde dehydrogenase, ALDH2). In contrast, proteins involved in tRNA aminoacylation for protein translation (e.g. mitochondrial tryptophan-tRNA ligase, WARS2, and mitochondrial methionine-tRNA ligase, MARS2) and 5’-3’ exodeoxyribonuclease activity (e.g. aprataxin, APTX) were less plenteous in elderly low-risk patients (Figure 2B, bottom plot).

Figure 1. Overview of the high-risk and low-risk AML patient cohort and the liquid chromatography tandem mass spectrometry (LC-MS/MS) workflow for the proteome and phosphoproteome analysis. The study included AML cell samples from 15 high-risk and 18 low-risk patients collected at the time of first diagnosis. Patients were classified after cytogenetic and molecular genetic analyses. Low-risk patients were further split into elderly low-risk and younger low-risk patients. AML sample preparation steps for proteome and phosphoproteome analysis included AML cell enrichment by density gradient separation, genetic analyses with classification of patients, cell lysis, addition of the super-SILAC (stable isotope labeling with amino acids in cell culture) mix, filter-aided sample preparation (FASP)-based protein digestion and additional immobilized metal affinity chromatography (IMAC) enrichment of phosphopeptides before data-dependent acquisition (DDA) on the mass spectrometer. *Median age of each patient group or subgroup.

Figure 2. The regulated proteome in the study of elderly low-risk vs younger low-risk patients. (A) Overview of proteomic data analysis from elderly low-risk and younger low-risk patient samples. Volcano plot analysis of proteins quantified in at least five patients per group. Points (in magenta) above the non-axial horizontal grey line represent proteins with significantly different abundances (P <0.05). Heatmap of the 29 differentially expressed proteins in the elderly low-risk and younger low-risk groups. The log₂ of fold change (FC) of protein levels in the elderly low-risk relative to the younger low-risk group is displayed on the right of the heatmap. (B) Gene Ontology (GO; BP, biological processes, with lilac bars; CC, cellular compartments, with light grey bars; MF, molecular functions, with yellow bars) and KEGG pathways (orange bars) analyses of upregulated and downregulated proteins in the elderly low-risk group. The various enriched GO terms and KEGG pathways are displayed on the y-axis while the corresponding –log₁₀P values are shown on the x-axis. The number of genes associated to a specific GO term or KEGG pathway is shown on the right side of the corresponding bar. Abbreviations were used in cases of long GO term (Immunoglob. for Immunoglobulin; Pos. for Positive and Neg. for Negative).

It can be seen from Supplementary Table 3 that a major part of the proteins with differential expression is involved in either transcriptional regulation (five proteins), protein homeostasis/modulation (10 proteins) or mitochondrial functions/metabolism (10 proteins). Altered epigenetic/transcriptional regulation is regarded as a characteristic feature of aging in hematopoietic stem cells [7, 22, 23] together with alterations in metabolism [24–26] and/or mitochondrial functions [27–30], protein homeostasis [13, 22] and DNA repair/genomic instability [7, 23, 31]. Thus, even though the comparative analysis of the proteome from elderly and younger low-risk AML patients showed only 29 differentially expressed proteins, most of these proteins are involved in the regulation of cellular processes that are known to be altered in aging cells. Among them, NUMB (protein numb homolog) is of particular interest because it is involved in the regulation of asymmetrical cell division [32] and altered cell polarity is a hallmark of aging [31].

The effect of aging on the AML cell phosphoproteome; a comparison of elderly and younger low-risk patients

We constructed a dataset comprising 14,574 identified phosphopeptides, from which 11,962 class I protein phosphorylation sites were quantified on 2818 proteins of nine elderly low-risk and nine younger low-risk patients. We identified 105 differentially regulated phosphorylation sites based on statistical analysis of 4767 phosphosites that could be quantified in at least five patients in each of these two groups (Figure 3A and Supplementary File 2). A cluster including 43 significantly upregulated and another cluster including 62 downregulated phosphosites in the elderly relative to the younger low-risk group were subjected to GO and KEGG pathways enrichment analysis (Supplementary Tables 3, 5, 6).

Figure 3. The regulated phosphoproteome in the study of elderly low-risk vs younger low-risk patients. (A) Overview of phosphoproteomic data analysis from elderly low-risk and younger low-risk patient samples. Volcano plot analysis of phosphosites quantified in at least five patients per group. Points (in magenta) above the non-axial horizontal grey line represent phosphosites with significantly different phosphorylation levels (P <0.05). (B) GO (BP terms with lilac bars; CC terms with light grey bars; MF terms with yellow bars) and KEGG pathways (orange bars) analyses of proteins with increased and decreased phosphorylation in elderly low-risk patients. The various enriched GO terms and KEGG pathways are displayed on the y-axis while the corresponding –log₁₀P values are shown on the x-axis. The number of genes associated to a specific GO term or KEGG pathway is shown on the right side of the corresponding bar. (C) Networks of protein-protein interactions (PPI) based on STRING database and visualized in Cytoscape after ClusterONE analysis. The significance of networks with high cohesiveness is shown with the P value of a one-sided Mann-Whitney U test. The differentially regulated phosphorylation sites are shown in yellow boxes next to each protein. FC of phosphorylation are color-coded; red-colored proteins showed a higher phosphorylation in elderly low-risk patients and blue-colored proteins showed a higher phosphorylation in the younger low-risk group. (D) Sequence motif analysis of the ± eight amino acids flanking the differentially regulated phosphorylation sites for each patient group and kinase-substrate enrichment analysis (KSEA) of differentially regulated and unregulated phosphorylation sites. The kinase z-score (x-axis) is the normalized score for each kinase (y-axis), weighted by the number of identified substrates indicated on the right side of the plot. Only significant predicted kinases with false discovery rate (FDR) values <0.05 were shown.

Regulation of filopodium assembly, vesicle cytoskeletal trafficking and RNA polymerase II transcription cofactor activity (Figure 3B, upper plot) represented GO terms with higher protein phosphorylation in elderly low-risk patients. These terms include phosphoproteins such as formin-binding protein 1-like (FNBP1L), protein kinase C-binding protein 1 (ZMYND8), Rho family-interacting cell polarization regulator 2 (FAM65B), dysbindin (DTNBP1), TP53-binding protein 1 (TP53BP1) and arginine-glutamic acid dipeptide repeats protein (RERE). The phosphoproteins enriched in elderly low-risk patients were primarily located in nuclear chromosomes, WASH complex, transverse tubules and endosome membranes. ARID1A (AT-rich interactive domain-containing protein 1A) may be of particular importance because this protein is a regulator of CDC42 (cell division control protein 42 homolog) which plays an essential role in the regulation of actin and tubulin organization and thereby regulation of cellular polarity that can be lost in aging hematopoietic cells [13, 22]. KEGG pathways analysis showed similarities with the results from the GO term enrichment, i.e. the upregulated elderly low-risk phosphoproteome was enriched for protein export pathways. Several phosphoproteins involved in nitrogen compound metabolic process (e.g. histone-lysine N-methyltransferase 2A, KMT2A, and eukaryotic translation initiation factor 3 subunit F, EIF3F), regulation of cell cycle arrest (e.g. cyclin-dependent kinase 1 and 2, CDK1 and CDK2; nucleophosmin, NPM1) and transcriptional activator activity (e.g. ribosomal RNA processing protein 1 homolog B, RRP1B) were less phosphorylated in elderly low-risk patients (Figure 3B, bottom plot). Proteins less phosphorylated in this group showed a significant enrichment of DNA replication and p53 signaling KEGG pathways.

Several protein-protein interactions (PPI) networks of significant cohesiveness were found after ClusterONE analysis based on STRING interactions of differentially phosphorylated proteins (Figure 3C). The most significant network (cluster 1) consisted of three phosphoproteins with higher phosphorylation in elderly low-risk patients and being involved in plasma membrane organization and endocytosis (i.e. protein export found in the KEGG analysis). The second significant cluster involved mRNA stabilization and splicing proteins. Serine/arginine repetitive matrix protein 1 and 2 (SRRM1 and SRRM2) showed higher phosphorylation on multiple sites in younger low-risk relative to elderly low-risk patients. These splicing factors are phosphorylated on multiple serine and threonine residues by dual specificity tyrosine-phosphorylation-regulated kinase 3 (DYRK3) during the G2-to-M transition, after the nuclear-envelope breakdown [33] (i.e. regulation of transcription and cell cycle, see the KEGG analyses).

As shown in Supplementary Table 3, a majority of the total set of the differentially regulated phosphorylated sites were located in proteins involved in RNA synthesis/function (i.e. transcription factors, epigenetic regulation, histone modulation, RNA splicing and ribosomal regulation). Such altered transcriptional regulation is regarded as a hallmark of aging [31]; it is also a characteristic of hematopoietic stem cell aging [19, 20] and thus an important characteristic of AML in elderly patients [7].

Differences in AML cell kinase activity when comparing elderly low-risk and younger low-risk patients

To identify protein kinases differentially activated in elderly low-risk and younger low-risk patients we performed phosphorylation site motif analysis with the WebLogo tool [34] (Figure 3D, left plots). We found higher activity of several kinases such as protein kinase A and C (PRKACA, PRKCD; basophilic motif upstream to the differentially phosphorylated site), calmodulin-dependent protein kinase II (CaM kinase II; basophilic motif upstream to the differentially phosphorylated site), serine/threonine protein kinase PAK2 (basophilic motif upstream to the differentially phosphorylated site), casein kinase 2 (CSK2; acidic amino acid-based motif downstream to the differentially phosphorylated site) and extracellular signal-regulated kinases (ERK1/2; proline-directed motif downstream to the differentially phosphorylated site) in elderly low-risk patients. In contrast, cyclin-dependent kinases (CDKs; proline-directed motif downstream to the differentially phosphorylated site) substrates appeared more phosphorylated in the younger low-risk group.

The kinase-substrate enrichment analysis (KSEA) [35, 36], which is based on phosphorylation fold changes (FCs) to estimate kinase’s activity, confirmed the higher activity of PAK2 and mitogen-activated protein kinases (MAPKs) in elderly low-risk patients and the higher activity of CDKs in the younger low-risk group (Figure 3D, right plot). KSEA analysis revealed a high number of CDK1 and CDK2 substrates from proteins such as ARID1A, LIG1 (DNA ligase 1), FOXK2 (forkhead box protein K2), NPM1 and retinoblastoma-associated protein (RB1).

We found five phosphosites in three different protein kinases in this data set using the activation loop analysis tool (see Materials and Methods). All of them (CDK1 T14, CDK2 T14, CDK1 Y15, CDK2 Y15 and GSG2 S147) can be phosphorylated during mitosis and were upregulated in younger low-risk relative to elderly low-risk patients.

Proteomic comparison of high-risk and all low-risk patients; different expression of neutrophil degranulation, platelet degranulation and cytoskeleton proteins

We compared the proteome profiles of AML cells derived from 15 high-risk and 18 low-risk patients. We quantified 6569 proteins, of which 5009 had a quantitative value in at least five patients in each group. The t-test based statistical analysis resulted in 205 differentially expressed proteins, 82 proteins were upregulated and 123 were downregulated in high-risk relative to low-risk patients (Supplementary Figure 1A and Supplementary File 3). Hierarchical clustering based on the 205 regulated proteins identified two main patient clusters (Figure 4A, left part), which corresponded to the high-risk and low-risk samples, although a distinct separation was not obtained. GO enrichment analysis showed that cytoskeleton organization, actin binding and integrin binding as biological process and molecular function GO terms were over-represented in high-risk patients (Figure 4A, right part). These terms include several tubulin (TUBB) chains, protein-tyrosine kinase 2-beta (PTK2B), cytoplasmic linker-associated protein 1 (CLASP1), actin-related protein 2/3 complex subunit 1A (ARPC1A), NCK-interacting protein with SH3 domain (NCKIPSD), TNF receptor-associated factor 2 (TRAF2) and NCK-interacting protein kinase (TNIK), alpha-actinin-1 (ACTN1) and integrin-linked protein kinase (ILK). Moreover, KEGG pathways analysis confirmed that the focal adhesion and pathogenic E. coli infection pathways (i.e. several tubulin chains and actin-related protein 2/3 complex subunits) were enriched in the upregulated high-risk proteome.

Figure 4. Proteomic differences between high-risk and low-risk patients; the importance of the cytoskeleton reflected in levels of neutrophil degranulation, platelet degranulation and endomembrane trafficking proteins. (A) Hierarchical clustering of 33 patients was based on the expression (SILAC log₂ ratio) of 205 proteins with significantly different regulation in AML cells from high-risk (dark grey squares) and low-risk patients (light grey squares). Two vertical main clusters were observed, one dominated by proteins with higher abundance in mostly high-risk patients (upper cluster) and the other by proteins with higher abundance in low-risk patients (lower cluster). GO and KEGG pathways analyses of the two protein clusters were performed to reveal enriched BP, CC and MF terms for the high-risk and low-risk patients. The various enriched GO terms and KEGG pathways are displayed in the scatter plot. The number of genes associated to a specific GO term or KEGG pathway (count) and the corresponding –log₁₀P values are shown on the x-axis and y-axis, respectively. Abbreviations were used in cases of long GO term or KEGG pathway name (Pathogen. for Pathogenic; fil. for filament). (B) Reactome term enrichment was performed using the STRING app (1.5.1) in Cytoscape. The five Reactome pathways with highest significance are shown with the corresponding FDR values. The protein nodes are colored according to their high-risk/low-risk FC, i.e. orange indicates increased abundance in the high-risk group and blue increased abundance in the low-risk group.

In the low-risk proteome, proteins involved in neutrophil degranulation (i.e. extracellular secretion) and lipid metabolic processes such as complement factor D (CFD), myeloperoxidase (MPO), lysosome-associated membrane glycoprotein 1 and 2 (LAMP1/2), mitochondrial enoyl-CoA delta isomerase 1 (ECI1), 3-ketoacyl-CoA thilase, mitochondrial (ACAA2) and lysophosphatidylcholine acyltransferase 2 (LPCAT2) were higher expressed (Figure 4A, right part). Arylsulfatase and procollagen galactosyltransferase activities were the molecular functions more abundant in the low-risk group. As much as 43% of the proteins with higher abundance in low-risk patients were annotated to the organelle membrane in the cellular component analysis.

Analysis of Reactome terms showed that neutrophil degranulation (Reactome pathways I), platelet degranulation (Reactome pathways II) -both these granulation processes reflecting the extracellular secretion in myeloid cells-, post-translational protein modification (Reactome pathways III), COPI-dependent Golgi-to-ER retrograde traffic (Reactome pathways IV) and cell-extracellular matrix interactions (Reactome pathways V) pathways were enriched in the set of 205 regulated proteins (Figure 4B and Supplementary Table 7). We observed several clusters of upregulated proteins in the high-risk patient group in Reactome pathways III, IV and V including cytoskeleton proteins such as tubulin chains, ACTN1, ILK, beta-parvin (PARVB) and Ras suppressor protein 1 (RSU1). Clusters of upregulated proteins in the low-risk patient group were found in other Reactome pathways, such as the one comprised of arylsulfatase A (ARSA), arylsulfatase B (ARSB), MPO, N(4)-(beta-N-acetylglucosaminyl)-L-asparaginase (AGA), beta-glucuronidase (GUSB) and serine/threonine-protein kinase 11-interacting protein (STK11IP) in Reactome pathways I; ARSA, ARSB and inactive C-alpha-formylglycine-generating enzyme 2 (SUMF2) in Reactome pathways III; and syntaxin-18 (STX18), neuroblastoma-amplified sequence (NBAS) and coatomer subunit gamma-2 (COPG2) in Reactome pathways IV.

Few proteins in Reactome pathways I and II show an expression mainly or limited to myeloid cells, whereas a majority of them are expressed by various tissues/organs/cells and are important for membrane/organellar functions or protein metabolism/modulation. Thus, these networks seem to mainly reflect differences in fundamental and common cellular processes more than differences in the differentiation status of the cells.

Only ALDH2 had a similar differential expression both when comparing elderly low-risk vs younger low-risk and high-risk vs all low-risk patients

A main difference between our first comparison of elderly low-risk vs younger low-risk patients and this second comparison of high-risk vs all low-risk patients is the higher number of differentially expressed proteins in the latter analysis. Furthermore, only three proteins, ALDH2, Ufm1-specific protease 2 (UFSP2) and BAG family molecular chaperone regulator 2 (BAG2) were differentially regulated in both studies (Supplementary Figure 2A). Only ALDH2, a poor prognosis predictor in AML as well as in urothelial cancer ([37]; http://www.proteinatlas.org) showed a similar upregulation in both comparisons. Western blot analyses using lysates from nine high-risk and nine low-risk (five elderly low-risk and 4 younger low-risk) patient cells showed higher ALDH2 expression in high-risk and elderly low-risk when compared to all low-risk and younger low-risk patients, respectively (Supplementary Figure 3), although differences between groups were not statistically significant according to the Mann-Whitney test.

These observations are consistent with the hypothesis that ALDH2 expression reflects an adverse prognostic impact of age in AML.

Phosphoproteomic comparison of high-risk vs all low-risk patients; detection of differences in mitotic cell cycle regulation

We identified and quantified 14,990 class I protein phosphorylation sites from 3279 proteins when comparing 15 high-risk and 18 low-risk patients. We found 239 differentially regulated phosphorylated sites based on statistical analysis of 6682 phosphosites, which were quantified in at least five patients in each group (Supplementary Figure 1B and Supplementary File 4). Hierarchical clustering using these 239 phosphosites clearly distinguished the phosphoproteome of the two patient groups (Figure 5A, left part). Two clusters, one containing 124 phosphosites and another with 115 phosphosites, were upregulated in the high-risk and in the low-risk patient group, respectively.

Figure 5. Phosphoproteomic differences between high-risk and low-risk patients; the importance of the cytoskeleton, mitotic cell cycle regulation and CDK activities. (A) Hierarchical clustering of the 33 patients based on the phosphorylation level (SILAC log₂ ratio) of 239 phosphosites with significant differences between high-risk and low-risk patient samples. Two vertical main clusters were observed, one dominated by phosphosites with higher phosphorylation in low-risk patients (upper cluster) and the other by phosphosites with higher phosphorylation in high-risk patients (lower cluster). GO and KEGG pathways analyses of the two corresponding phosphoprotein clusters were performed to reveal enriched BP, CC and MF terms in the high-risk and low-risk patient samples. The various enriched GO terms and KEGG pathways are displayed in the scatter plot. The number of genes associated to a specific GO term or KEGG pathways (count) and the corresponding –log₁₀P values are shown on the x-axis and y-axis, respectively. Abbreviations were used in cases of long GO term or KEGG pathway name (reg. for regulation). (B) Visualization of hit Reactome pathways was performed using the ReactomeFIViz app (7.2.3) in Cytoscape. Two significant Reactome networks (FDR <0.05) that mapped phosphoproteins with differential phosphorylation in our dataset are shown. The protein nodes are colored according to their high-risk vs low-risk log₂ phosphorylation FC, i.e. orange indicates increased phosphorylation in the high-risk group and blue increased phosphorylation in the low-risk group. (C) Networks of PPI based on STRING database and visualized in Cytoscape after ClusterONE analysis. The significance of networks with high cohesiveness is shown with the P value of a one-sided Mann-Whitney U test. The differentially regulated phosphorylation sites are shown in light green boxes next to each protein. FC of phosphorylation are color-coded; orange-colored proteins showed a higher phosphorylation in the high-risk group and blue-colored proteins showed a higher phosphorylation in the low-risk group. (D) Causal relationships between phosphoproteins with differentially regulated phosphorylation sites in the high-risk vs low-risk phosphoproteome set was studied with SIGNOR. The analysis showed the pivotal role of CDKs in the control of cell cycle, cytoskeleton and translation phosphoproteins. Nodes and types of relationships are displayed as indicated on the bottom part of the design.

Cellular component GO analysis revealed an enrichment of upregulated cytoplasmic, cytoskeleton and membrane phosphoproteins for the high-risk patients, whereas organelle, organelle envelope and nuclear envelope structures were enriched in low-risk patients (Figure 5A, right part). While cell adhesion molecule binding, positive regulation intracellular signal transduction and small molecule metabolic process were the biological process and molecular function GO terms enriched in the high-risk group, positive regulation chromosome segregation, cell cycle checkpoint and Rho guanyl-nucleotide exchange factor activity were enriched in the low-risk patient group (Figure 5A, right part).

Reactome pathways for apoptotic execution phase and for mitotic prophase were found significantly enriched with proteins mainly of higher phosphorylation in the high-risk and in the low-risk patient group, respectively (Figure 5B and Supplementary Table 7). Two PPI networks of significant cohesiveness were found after ClusterONE analysis based on STRING interactions of differentially phosphorylated proteins (Figure 5C and Supplementary Table 7). The most significant network (cluster 2) consisted of six phosphoproteins involved in RNA processing and mRNA splicing, most of them with higher phosphorylation in low-risk patients. The other significant cluster included phosphoproteins of the mitotic cell cycle process (cluster 1). All phosphoproteins, except nuclear mitotic apparatus protein 1 (NUMA1), were significantly more phosphorylated in high-risk patients. A sequence logo analysis of the amino acids surrounding the phosphosites in cluster 1 suggested PRKACA and kinases of the PRKC family involved in the phosphorylation of the six cell cycle proteins (Supplementary Figure 4).

A final analysis of the regulated high-risk vs low-risk phosphosite set to study the interactions of their corresponding phosphoproteins in signal transduction with SIGNOR [38] (Figure 5D) confirmed the relevant roles of CDKs in the regulation of the cell cycle, RNA processing, translation and cytoskeleton function that are observed in AML patients with different risk-related cytogenetics abnormalities.

Several protein kinases are differentially activated in AML cells derived from high-risk and all low-risk patients

To identify protein kinases differentially activated in the two groups we performed phosphorylation site motif analysis with IceLogo [39]. We found a basophilic motif upstream to the differentially phosphorylated site in high-risk patients when compared to low-risk patients, suggesting an activation of PRKACA, PRKCA and PRKCD (Figure 6A). Furthermore, KSEA confirmed the higher activity of PRKACA and predicted several other serin/threonine protein kinases (PRKG1, PRKD1 and PRKCD), MAPKs and RAC-alpha serine/threonine protein kinase (AKT) isoforms activated in the high-risk patient group (Figure 6B). Although PRKACA and AKT1 phosphorylated a large number of substrates (52 and 51, respectively), CDK1 (significantly predicted with unadjusted P= 0.022 in the KSEA; therefore not displayed in Figure 6B which shows predicted kinases at FDR <0.05) phosphorylated 140 substrates in this group. Serine/threonine protein kinase MAK, which phosphorylates proteins involved in protein ubiquitation and in the transcriptional coactivation of androgen receptor, was activated in the low-risk patient group.

Figure 6. Kinase prediction analysis of the high-risk vs low-risk phosphoproteome. (A) Sequence motif analysis of the ± six amino acids flanking the differentially regulated phosphorylation sites for either group. (B) KSEA of differentially regulated and unregulated phosphorylation sites. The kinase z-score (x-axis) is the normalized score for each kinase (y-axis), weighted by the number of identified substrates indicated on the right side of the plot. Significant predicted kinases with FDR <0.05 are shown.

We found eight phosphosites on five different protein kinases in the data set of 239 differentially regulated phosphorylation sites. Three phosphosites on serine/threonine-protein kinase D2 (PRKD2; S197, S198 and S706; the latter located in the activation loop of the kinase), involved in the regulation of cell proliferation via ERK1/2 signaling, in Golgi membrane trafficking and cell adhesion, were significantly more phosphorylated in the high-risk patient group. PRKD2 S706 is probably phosphorylated by the members of the PRKC family [40]. SRC S17 and CDK1 Y15 were also significantly more phosphorylated in the high-risk patient group. Phosphorylation of CDK1 Y15 by Wee-1 like protein kinase 1/2 (WEE1/2) inhibits the protein kinase activity and acts as a negative regulator of entry into mitosis (G2 to M transition) whilst phosphorylation by PRKCD activates the G2-M DNA damage checkpoint after UV irradiation [41]. Dephosphorylation at CDK1 Y15 by active M-phase inducer phosphatase 1/2 (CDC25A/B) leads to CDK1 activation at the G2-M transition [42]. The higher phosphorylation of CDK1 (and CDK2) Y15, along with T14 phosphorylation, was further confirmed by a separate MS-based immune-affinity enrichment (Supplementary Table 8 and Supplementary File 5). The significantly higher phosphorylation of CDK1/2 T14 alone in high-risk patients was also observed in the latter analysis (not significant in the general phosphopeptide enrichment analysis with P= 0.078). Finally, tyrosine-protein kinase Fes/Fps (FES), a regulator of the actin cytoskeleton and microtubule assembly, S408; and STK26 (serine/threonine-protein kinase 26), a mediator of cell growth, T327 and T328 were significantly more phosphorylated in the low-risk patient group.

The results from Western blotting did not show any significant difference in protein phosphorylation of several kinase candidates between the high-risk and the low-risk groups (Supplementary Figure 5). However, despite the different spectrum of sensitivity of both methods, the increased activity of CDK1, PRKCD and PRKACA in the high-risk group detected by the MS-based data was also observed in the Western blots (Supplementary Figure 5). We also noticed a distinct CDK1 T161-Y15 phosphorylation pattern for each of the patient groups. While these two phosphosites were similarly phosphorylated in low-risk patients, CDK1 T161 phosphorylation was different of that observed in CDK1 Y15 in high-risk patients suggesting a different CDK1 regulation for each of the groups. Moreover, the phosphorylation levels of PRKCD S645 and PRKACA T197 were parallel for each patient showing a joint activation of both kinases.

Only LSP1 had a similar differential phosphorylation both when comparing elderly low-risk vs younger low-risk and high-risk vs all low-risk patients

Five phosphosites on proteins lymphocyte-specific protein 1 (LSP1), lamina-associated polypeptide 2, isoform alpha (TMPO), CDK1, CDK2 and N-acetyltransferase ESCO2 (ESCO2) were quantified in both our comparisons, i.e. elderly low-risk vs younger low-risk and high-risk vs all low-risk patients (Supplementary Figure 2B). Only LSP1_S177, an actin-binding cytoskeleton protein involved in cell migration and possibly phosphorylated at that residue by MAPKAPK2 and/or protein kinase C [43], showed significant associations with high age/increased risk of relapse in both our comparisons. However, CDK1/2 Y15 showed increased phosphorylation in high-risk and in younger low-risk patients. These apparently conflicting observations on CDK’s activity illustrate that the impact of certain age-associated proteomic differences may differ between patient subsets and depend on the biological/genetic context.

Selection of patients and preparation of patient cells

All cell samples in the present study were derived at the first time of diagnosis before start of any antileukemic treatment. Our institution is responsible for diagnosis and treatment of patients with AML in a defined geographical area, and the patients included in our present study represent all patients from this area during a defined time period that fulfilled defined criteria established before start of the study. First, to ensure a high quality of the analyzed samples with at least 95% AML cells we included only patients with high levels of peripheral blood AML cells, i.e. at least 10 x 10⁹ total leukocytes and at least 80% of these circulating leukocytes being leukemic cells. Second, among these consecutive patients with high levels of circulating AML cells, we selected all patients who fulfilled predefined genetic criteria for adverse and favorable prognosis as can been seen from Supplementary Table 1. These criteria are based on the European Leukemia Net (ELN) guidelines [1]. The ELN classification of AML is based on karyotyping together with molecular-genetic analyses, and patients are classified as having favorable, intermediate or adverse prognosis. This risk stratification does not take into account other pretreatment risk factors like peripheral blood blast count at the time of diagnosis, previous hematological malignancy (i.e. myelodysplastic syndrome, chronic myelomonocytic leukemia or chronic myeloproliferative neoplasia) and previous exposure to cytotoxic therapy for other diseases [1]. The classification does not take into account response to first induction cycle and minimal residual disease (MRD) after remission induction either [1, 45, 106, 107].

The patients with adverse prognosis included ten patients with complex karyotype and five patients with either monosomal karyotype, del 5 or -7. These adverse prognosis patients represent all patients who fulfilled the cytogenetic criteria for high-risk disease. Patients were screened for FLT3 and NPM1 mutations but additional mutational analyses for identification of additional adverse prognosis patients were not available during the defined time period, and for this reason these adverse prognosis patients were selected based on cytogenetic criteria alone. These patients will be referred to as high-risk patients. The median age of this patient group is 74 years.

The patients with favorable prognosis included all patients with the cytogenetic abnormalities inv(16), t(16;16) and t(8;21). We also included all patients with normal karyotype, wild-type FLT3 and NPM1 insertion as well as one patient with normal karyotype, low ratio of internal tandem duplication mutation of FLT3, FLT3-ITD, and NPM1 insertion (Supplementary Table 1, patient F17) and another patient with CEBPA mutation (Supplementary Table 1, patient F18). Analysis of the FLT3-ITD ratio was not available as a routine analysis during this time period and only one patient could be classified as favorable based on these criteria (Supplementary Table 1). These patients will be referred to as low-risk patients. The median age of this patient group is 64.5 years.

More characteristics of the high-risk and low-risk patient groups are summarized in Tables 1, 2, and more detailed comparisons are presented in Supplementary Tables 1, 2. All patients were Caucasians.

As we included patients with high levels of circulating AML blasts at the time of first diagnosis, we could prepare highly enriched AML cell populations (>95% purity) by a standardized method based on density gradient separation alone (for a detailed discussion see references [53, 54, 65]). All samples were cryopreserved by using the same standardized method. The patient samples included in our study did not differ significantly with regard to peripheral blood blast count or storage time in liquid nitrogen.

Patient grouping for MS-based proteomic and phosphoproteomic analysis

In order to compare patient groups with different cytogenetics-based prognosis, we selected 15 high-risk and 18 low-risk patients for subsequent proteomics analyses. The difference of the median age between the two groups was statistically significant (P= 0.005, Table 1). In order to study the influence of age in low-risk patients, we further divided the 18 patients in two subgroups of nine patients each, according to an age threshold of 65 years at the time of diagnosis, i.e. the elderly low-risk (median age of 68 years) and the younger low-risk (median age of 47 years). This threshold was chosen because it corresponds to the median age of patients at the first diagnosis of AML and it has also been used to distinguish between young/middle aged and elderly patients (i.e. not fit for the most intensive chemotherapy) in previous clinical studies [1, 108, 109]. The difference of the median age between the two subgroups was statistically significant (P <0.0001). This and other patient characteristics of the low-risk subgroups are shown in Table 2.

Peptide preparation

Our standard sample preparation of patient cell lysate in 4% sodium dodecyl sulfate (SDS)/0.1 M Tris-HCl (pH 7.6), the filter-aided sample preparation (FASP) procedure with AML patient samples and the immobilized metal affinity chromatography (IMAC) for phosphopeptide enrichment have been described previously [65, 110]. In short, 20 μg of each patient lysate was mixed with 10 μg of an AML-specific super-SILAC mix [111] for proteomic analyses, and digested according to the standard FASP protocol [110, 112]. The super-SILAC spiked peptide samples were fractionated using styrenedivinylbenzene-reversed phase sulfonate (SDB-RPS) plugs (Empore Discs, 3M) for proteome analysis [113]. The phosphoproteomics samples (226-2506 μg) were mixed with the super-SILAC mix at 1:2 ratio (w:w; super-SILAC mix:AML patient sample), FASP processed and enriched using the IMAC procedure. Extra labeled peptides samples (600-1479 μg) were also FASP prepared for phosphotyrosine immunoprecipitation (IP). IP was performed using PTMScan pTyr antibody beads (p-Tyr-1000; Cell Signaling Technology) according to the manufacturer’s protocol. All the peptides samples were brought to equal volumes of binding buffer before antibody incubation. Eluted Tyr-phosphopeptides were cleaned up with 3-C18-disks-stage tips before MS analysis.

LC-MS/MS measurements

Peptide sample preparation prior to proteome and phosphoproteome analysis and settings of the LC-MS/MS runs on a Q Exactive HF Orbitrap mass spectrometer coupled to an Ultimate 3000 Rapid Separation LC system (Thermo Scientific) were carried out as described earlier [65]. Tyr-phosphopeptides were pre-concentrated on a 2 cm × 75 μm ID Acclaim PepMap 100 trapping column and separated on a 25 cm × 75 μm ID EASY-spray PepMap RSLC analytical column (Thermo Scientific). The Tyr-phosphopeptides were eluted during a 105 min binary gradient with solvent A (0.1% formic acid) and solvent B (0.1% formic acid/acetonitrile). The gradient started at 5% B from 0–5 min and increased to 7% B from 5–5.5 min, then to 22% B from 5.5–65 min, to 35% B from 65–87 min, and to 90% B from 87–92 min. Hold at 90% from 92-102 min, then ramped to 5% B from 102–105 min. The Q Exactive HF mass spectrometer was operated in data dependent acquisition (DDA) mode. Full MS scans (scan range 375–1500 m/z) were acquired in profile mode with a resolution R = 60 000 using an AGC target value of 3 × 10⁶ charges. MS/MS scans were acquired in profile mode for the top 10 precursors. The AGC target value was set to 1 × 10⁵ charges with a maximum injection time of 110 ms and a resolution R = 60 000. The normalized collision energy was 28 and the isolation window was 1.2 m/z with null m/z offset. The dynamic exclusion lasted for 20s.

The super-SILAC proteomic, phosphoproteomic and Tyr-phosphoproteomic samples were analyzed as three separate experiments in a controlled randomized order (i.e. samples from each patient group were distributed equally over the analysis sequence) with a LC-MS quality control (HeLa protein digest) run approximately every 10 patient samples.

Data analysis

LC-MS/MS raw files were processed with MaxQuant software version 1.5.2.8 [114, 115]. The spectra were searched against the concatenated forward and reversed-decoy Swiss-Prot Homo sapiens database version 2018_02 using the Andromeda search engine [116]. The Perseus 1.6.1.1 platform was used to analyze and visualize protein groups and phosphosites [117]. MaxQuant-normalized SILAC ratios were inverted, log₂ transformed and normalized again using width adjustment. Hierarchical clustering of significantly differential proteins and phosphosites was done with Perseus using the Pearson correlation function and complete linkage. Volcano plots were done with Prism8 (GraphPad). GO and KEGG pathways analysis was performed using a GO tool [118]. The most significantly over-represented GO and KEGG pathways terms with P <0.05 were displayed in bar or scatter plots in Prism8. The amino acid distribution surrounding the phosphosites was analyzed using iceLogo (P= 0.05) with the sequence windows obtained in the MaxQuant-generated phosphosite output file [39]. Sequence logo analyses from a small number of phosphopeptide sequences were generated with WebLogo [34]. Kinase activity estimates were inferred by the KSEA App. [35, 36]. Regulated and unregulated phosphosites were analyzed with the PhosphoSitePlus [119] and NetworKin [120] databases using a substrate count and a NetworKin score cutoff of 5. Kinase activation loop analysis was performed with the tools for phosphoproteomics data analysis at http://phomics.jensenlab.org. PPI networks were obtained by using the STRING database version 11.0 with interactions derived from experiments and databases at a high confidence score of 0.7 [121]. Networks were visualized using the Cytoscape platform version 3.3.0 [122]. The ClusterONE plugin was used to identify protein groups of high cohesiveness [123]. Reactome term enrichment and visualization of hit pathways were performed using the STRING app (1.5.1) [124] and the ReactomeFIViz app (7.2.3) [125, 126], respectively. Causal relationships between phosphoproteins were studied with the SIGnaling Network Open Resource (SIGNOR) 2.0 [38]. Venn diagrams were made with Biovenn [127].

Western blotting

Western blotting from nine high-risk and nine low-risk patient cells were performed. Twenty μg of each SDS-based cell lysate was loaded on a NuPAGE 4-12% Bis-Tris protein gel (ThemoFisher Scientific) and transferred onto a nitrocellulose membrane (Amersham Protran, GE Healthcare Life Sciences) in an XCell II Blot Module (ThermoFisher Scientific). Antibodies were purchased from Cell Signaling Technology and Abcam. They were used according to manufacturer’s guidelines. Chemiluminescence was developed with SuperSignal West kits (Thermo Scientific) and measured on a LAS-3000 imager (Fujifilm). Band intensities for each protein were determined by densitometry software Image J [128]. Band intensities of a protein spotted at approximately 62 kDa on Ponceau-stained membranes were used for normalization.

Statistical analysis

Proteins and phosphosites (localization probability >0.75) with at least five individual SILAC ratios in each patient group were selected for two-sample unequal variance t-test and Z-statistics [129] to find significantly different FC for proteins and phosphosites between the different patient groups. Tyrosine phosphoproteomic data were tested for significance using linear models for microarray data (LIMMA) using the VSClust app [130]. Data from Western blotting bands were expressed as the median ± 95% confidence interval. Statistical analysis was performed using the Mann-Whitney test in Prism8.

Ethics statement

Primary AML cells were collected from AML patients after written informed consent in accordance with the Declaration of Helsinki. The storage of cells in our biobanks (REK 1759/2015) and the use of cells in the present project (REK 305/2017) were approved by the Regional Ethics Committee.

Data availability

The LC-MS/MS raw files and MaxQuant output files have been deposited to the ProteomeXchange consortium via the PRIDE partner repository [131, 132] with dataset identifier PXD019785.