Peripheral blood (PB) and bone marrow (BM) samples were collected from 41 adult patients with AML-M1/M2 and from 20 adult healthy volunteers (HV). Each individual provided signed informed consent for treatment and for their participation in this study. Appropriate approval was also obtained from the Bioethical Commission of the Karol Marcinkowski Poznan University of Medical Sciences, Poznan, Poland. The patients were diagnosed and treated at the Department of Hematology and Bone Marrow Transplantation at the Poznan University Hospital of the Lord’s Transfiguration of the University of Medical Sciences in Poznan, Poland. Standard AML therapy using cytosine arabinoside plus daunorubicin (3′+7′) was administered to all patients to induce complete remission (CR), which was defined according to the European Leukemia Net guidelines (2). When available, the samples were collected at the 3 following time-points: when AML was diagnosed, prior to first therapy (T0), when CR was established [between day +21 and +28 after the start of induction therapy (T1)], and when the disease relapsed (T2). However, from the T1 samples, the number of cells was much lower than that of cells from the T0 and T2 samples, and often, there was not sufficient material to perform replicate experiments. Moreover, not all micro-array images met the required quality criteria and had to be filtered out. The material from the T1 time-point was the most heterogeneous, and preliminary microarray analysis revealed that the T1 samples did not cluster together and did not exhibit clear common characteristics. The T2 samples were generally similar to the T0 samples and we were not able to distinguish them by unsupervised hierarchical clustering. Therefore, we decided to limit our analyses to the samples collected at the T0 time-point, with the use of available information about the history of treatment. The advantage of this selection was that the T0 samples were the most homogenous (the fractions of blood and bone marrow-derived mononuclear cells isolated from patients with AML-M1 contained ~90% leukemic cells, whereas from those from the patients with AML-M2 contained ~70% myeloblasts). Table I presents the summarized information on the patient and HV samples.

Mononuclear cells form the peripheral blood (termed PBMCs) and bone marrow (termed BMMCs) were separated through density gradient centrifugation (Gradisol L; Aqua-Medica, Bogucin, Poland) and washed 3 times with 1X phosphate-buffered saline (PBS; Ca and Mg-free; IBSS Biomed S.A., Warsaw, Poland). The cell pellet was suspended in lysis buffer from a mirVana miRNA Isolation kit (Ambion/Thermo Fisher Scientific, Waltham, MA, USA) and immediately frozen at −80°C.

Total RNA was extracted from the PBMCs and BMMCs using a mirVana miRNA Isolation kit and the DNase-treated (TURBO DNA-free kit) (both from Ambion/Thermo Fisher Scientific). RNA integrity was evaluated with the use of a Bioanalyzer 2100 and Total RNA Nano assay (Agilent Technologies, Santa Clara, CA, USA). Only RNAs with an RNA integrity number (RIN) ≥7.5 were used for the gene expression analyses.

Screening of the literature enabled us to identify the following groups of genes: i) Proven and postulated acute leukemia biomarkers; ii) general oncogenes; and iii) predicted to be specifically involved in leukemic transformation. In addition, we selected a set of control genes. Human housekeeping genes, not associated with oncogenesis, were used as the positive controls, whereas plant- and bacterial-specific genes served as the negative controls. We ordered 903 chemically synthesized, amino-modified microarray probes from 2 companies: Ocimum Biosolutions Ltd. (Hyderabad, India; 783 probes 50-nt-long) and Operon Biotechnologies GmbH (Cologne, Germany; 120 probes 70-nt-long). As an additional control, we used 4 short DNA probes (24-28-nt-long; Institute of Biochemistry and Biophysics, Warsaw, Poland), 8 probes complementary to external spike RNAs (AM1781; Ambion/Thermo Fisher Scientific) and one random oligo probe (Operon Biotechnologies GmbH). In total, the microarray contained 916 unique probes: 896 probes complementary to 838 human genes, 8 probes specific for 5 bacterial genes, 3 probes specific for 3 plant genes and 9 artificial control probes. The probes were diluted in Pronto Epoxide Spotting Solution to a 20 µM concentration and spotted in triplicates onto Epoxide-Coated Slides (both from Corning Inc., New York, NY, USA) with a SpotArray 24 instrument (Perkin-Elmer, Waltham, MA, USA). The design of our AML-array has been deposited in the ArrayExpress database (http://www.ebi.ac.uk/arrayexpress) under the accession no. A-MEXP-2220.

A quantity of 10 µg of each total RNA sample was reverse transcribed using anchoredoligo(dT)₂₀, aminoallyl-modified dNTPs and SuperScript III reverse transcriptase from a SuperScript Plus Indirect cDNA Labeling System (Invitrogen, Carlsbad, CA, USA). Amino-modified cDNAs from the patients with AML and HV were labeled with Alexa Fluor 647 and the reference cDNA from the HL60 cell line sample (obtained from Dr Marcin Schmidt from Poznan University of Life Sciences) with AlexaFluor 555. The labeled cDNA was purified (MinElute Reaction Cleanup kit; Qiagen, Hilden, Germany), dissolved in hybridization buffer (5X SSC, 0.1% SDS and 0.1 mg BSA/ml) and kept up to 30 min at 50°C prior to hybridization in Corning microarray hybridization chambers in a HybArray12 (Perkin-Elmer) or a water bath by using a step-down hybridization protocol (5 h/50°C, 5 h/45°C and 5 h/40°C). Three subsequent wash steps were applied: i) 2X SSC and 0.1% SDS at 40°C for 5 min; ii) 2X SSC at room temperature for 5 min; and iii) 0.2X SSC at room temperature for 5 min. The slides were dried through centrifugation with the Microarray High-Speed Centrifuge (Arrayit Corp., Sunnyvale, CA, USA; 5 sec, 2,000 × g, room temperature) and scanned with a ScanArrayExpress system (Perkin-Elmer) at 5-µm resolution.

The microarray images were processed using GenePix Pro version 6.0 software (Molecular Devices, LLC, Sunnyvale, CA, USA). Spots that did not meet a set of criteria (SNR >5; Dia., 110–250 nm; F CV ≤100; SD ≥1000; satur. ≤1%) were filtered out. The raw data files (deposited in the ArrayExpress database under the accession no. E-MTAB-5434) were loaded into R Bioconductor version 3.2.0 (R Development Core Team 2009) and processed and analyzed with the limma package (36). Briefly, the foreground and background median signal values (excluding those with weights of 0) were background-corrected with the subtraction method and normalized with the global loess method. Probed replicates were averaged, and the replicate arrays from the same patient were merged. Differential gene expression was tested using a linear model suitable for experiments with common reference design. The final p-values were adjusted with FDR correction (37). Data clustering was conducted and visualized with the Heatplus R package.

DNA-free RNA (3 µg per sample) was reverse transcribed using SuperScript III RT and oligo(dT) (Invitrogen). The reaction mixtures (20 µl vol) were incubated for 2 h at 50°C, and additional enzyme (1 µl) was added after the first hour. Following reverse transcription, the samples were incubated for 20 min at 70°C with 10 µl of 1 M NaOH. Subsequently, 10 µl of HCl 1 M were added for neutralization, and the cDNA was precipitated overnight at −20°C with 100 µl (2.5 vol) of 96% ethanol and 4 µl (1/10 vol) of 3 M sodium acetate, pH 5.2. The centrifuged pellet was washed twice with 70% ethanol and dissolved in with 60 µl of DEPC-H₂O. A total of 1 µl of each cDNA, diluted 3-fold, served as a template for real-time PCR in reactions, including MESA Green qRT-PCR MasterMix Plus (Eurogentec, Seraing, Belgium) and primers specific for the STMN1, NPM1, S100A8 and S100A9 transcripts. From 4 reference genes [ACTB, glyceraldehyde 3-phosphate dehydrogenase (GAPDH), PGK1 and PPIA], the last 2 were selected as the most stable using geNorm (https://genorm.cmgg.be). The sequences of all the primers are presented in Table II. Each gene analysis was performed in duplicate. The standards were generated by PCR amplification of the HeLa cDNA template, obtained as a result of reverse transcription of control HeLa RNA with a SuperScript Plus Indirect cDNA Labeling system (Invitrogen). Real-time analysis was performed in a Rotor-Gene Q thermocycler (Qiagen), using the following program: initial denaturation (94°C, 2 min), 40 cycles of denaturation (94°C, 1 min), annealing (60°C, 1 min) and elongation (72°C, 30 sec), final elongation (72°C, 10 min) and a final hold (4°C). The melting curve was generated in the range of 60–95°C. The reaction volume was 20 µl. Each reaction was performed in triplicate. The real-time PCR products were verified by melting curve analysis. The product quantities were estimated using the two standard curves’ method. The values obtained for the triplicates were averaged. Subsequently, the results obtained for the 2 biological replicates (separate reverse transcription reactions and separate real-time analyses) were averaged and divided by the average of the replicates of the 2 reference genes.

DNA-free RNA (2.5 µg per sample) was reverse transcribed and purified, as described above for real-time PCR, with the exception that the RT reaction lasted 1.5 h and only 1 µl of enzyme was used. The precipitated cDNA, dissolved in 50 µl of DEPC-H₂O, served as a template for quantitative PCR with the use of a QX200 ddPCR system, QX200 EvaGreen ddPCR Supermix (Bio-Rad, Hercules, CA, USA), and primers specific for 4 candidate genes (ABL1, ANXA3, CAT and WT1) and one reference gene (PGK1). The reaction volume was 20 µl, with a primer concentration of 250 nM, and the volume of cDNA (1–3 µl) was optimized for each candidate gene and sample. The PCR conditions were as follows: initial denaturation (95°C, 5 min), 40 cycles of denaturation (95°C, 30 sec), annealing (60°C, 30 sec) and elongation (72°C, 45 sec), cooling (4°C, 5 min), final denaturation (90°C, 5 min) and final hold (12°C). The temperature ramping rate was 2°C/sec. Following PCR, the plates were directly analyzed with a QX200 Droplet Reader. The data were processed using Quanta Soft version 1.5.38.1118 software (Bio-Rad). The number of droplets for each candidate gene was divided by the number of droplets obtained for the reference gene analyzed in the same run (for ABL1, ANXA3 and WT1) or in the same multiplex reaction (CAT). Each reaction was performed in duplicate. The sequences of all the primers are presented in Table II.

Statistical analyses and bar plots were made in R version 3.2.0 R and R Studio version 0.98.1102. Welch two sample t-tests (unpaired) were applied for pairwise comparisons and ANOVA was applied for comparisons across 2 or more groups. For two- or more-way ANOVA, Tukey honestly significant difference (HSD) tests based on multiple comparisons of means were applied to determine which pairwise comparisons were statistically significant. To test the correlations between the expression values of 2 genes or between a particular gene’s expression and clinical data (WBC, age), Pearson’s correlations were calculated. The threshold p-value was always set as 0.05.

The following open-source software tools and databases were explored: NCBI databases (www.ncbi.nlm.nih.gov/), including PubMed, Gene and Gene Expression Omnibus (GEO) and GDSbrowser, Expression Atlas (www.ebi-ac-uk/gxa/) and DAVID (david.abcc.ncifcrf.gov/).

To study gene expression in AML, we designed a small boutique AML-array, as deposited in the ArrayExpress database (A-MEXP-2220). The array contained >900 long oligonucleotide probes that are complementary to human genes implicated, inter alia, in hematopoietic cell differentiation and maturation, proliferation, apoptosis and leukemic transformation. Gene functional analysis with DAVID (38) revealed that our microarray was significantly enriched in genes linked with immune, infection and cancer classes of disease. Out of OMIM diseases, only AML was found to be statistically significant. The genes covered by the AML-array encode proteins from the following KEGG pathways: Cancer (111 genes), cytokine-cytokine receptor interaction (102 genes), Jak-STAT signaling (64 genes), chemokine signaling (59 genes), MAPK signaling (54 genes), focal adhesion (49 genes), hematopoietic cell lineage (42 genes) and regulation of actin cytoskeleton (41 genes).

The AML-array was used to estimate relative gene expression in PBMCs and/or BMMCs obtained from 33 de novo AML patients, prior to therapy, and 15 HV (Table I). The AML samples were classified as M1 or M2 FAB subtypes, which were characterized on the basis of the early blast cell differentiation arrest and either no myeloid cell maturation (M1) or partial maturation (M2). The microarray experiment was based on two-color hybridization, wherein each studied sample was matched with a common reference (HL60 cell line). Data from 118 hybridizations, as deposited in the ArrayExpress database under the accession no. E-MTAB-5434, were background-corrected and normalized, and technical replicates (of probes and arrays) were then merged and analyzed using the R Bioconductor limma package.

First, we assessed whether there were any differences between the transcriptomes of the BMMCs and PBMCs. Similar to the proteomics analysis results of our previous study (33), there were no statistically significant differences between the AML-PBMCs and AML-BMMCs (differential expression analysis with a limma linear model, adjusted p-value >0.11) (Table III). Therefore, in the following transcriptome analyses, we did not divide the samples into BMMC- and PBMC-derived groups.

Our next aim was to identify the differences between the AML and HV samples. Assuming a significance level α 0.05 as the threshold, we selected 163 differentially expressed genes (DEGs), of which 78 were underexpressed and 85 were overexpressed in AML (Table III). However, after sample clustering with these DEGs, 2 samples were still misclassified (one AML sample was included in the HV cluster, and one HV was in the AML cluster). The application of a more restrictive threshold (α 0.01) limited the list of DEGs to 83 (42 underexpressed and 41 overexpressed) (Table IV). In total, 36 genes from the list (43%, up- and downregulated in equal proportions) encoded proteins involved in signal transduction and transcriptional regulation, affecting processes, such as cell growth, proliferation, adhesion and apoptosis. Of note, 10 genes (12%, 4 up- and 6 downregulated) were associated with cytoskeletal organization and functioning, and 12 genes (14,5%, exclusively underexpressed in AML) were engaged in immune and inflammatory response. With the use of these 83 DEGs, we were able to perfectly distinguish the AML from the HV samples (Fig. 1A); the entire set of arrays was clearly divided into 2 clusters: 15 control HV samples constituted a separate cluster, whereas all 33 AML samples belonged to a second cluster.

Considering another criterion, the level of expression, we selected from the list of 83 DEGs, 23 genes with the greatest change in expression (log₂FC values >1.0 or <−1.0, equivalent to a >2-fold change in expression in both directions): 12 up- and 11 downregulated (Fig. 1B). The most significantly overexpressed gene in AML (log₂FC close to 3, equal to an 8-fold enrichment) was STMN1, which encodes a protein involved in the regulation of the microtubule filament system, stathmin 1. STMN1 was also the top statistically significant gene in the comparison between the AML and HV samples (adjusted p-value 1.67×10⁻¹⁰) (Tables III and IV). Other upregulated genes encoded proteins associated with myeloid leukemia and cancer development: Oncogenes, cell cycle and apoptosis regulators [CCAAT/enhancer binding protein α (CEBPA), Kirsten rat sarcoma viral oncogene homolog (KRAS), v-myc myelocytomatosis viral oncogene homolog, avian (MYC), cyclin-dependent kinase 6 (CDK6), minichromosome maintenance complex component 5 (MCM5), SET nuclear proto-oncogene (SET), proto-oncogene c-Kit ligand (KITLG) and serglycin (SRGN)], proteins involved in angiogenesis [angiopoietin 1 (ANGPT1)], cell growth, glucose metabolism and the transcriptional regulator, enolase 1 (ENO1) and ribosomal protein RPLP0.

Among the downregulated genes, predominantly those encoding defense/immunity proteins, were the following: Cell surface receptors CD74 and Fc fragment of IgE (FCER1G), ficolin 1 (FCN1) typically expressed in the peripheral blood leukocytes, antimicrobial agent lysozyme (LYZ), tumor necrosis factor (TNF) family member lymphotoxin β (LTB), interferon induced transmembrane protein 1 (IFITM1) and calcium binding protein S100A9. Moreover, decreased gene expression in AML was observed for the cell cycle regulator protein tyrosine phosphatase (PTPRE), apoptosis suppressor baculoviral IAP repeat-containing 3 (BIRC3), regulator of actin polymerization, cell proliferation, migration and differentiation thymosin β4 (TMSB4X) and membrane skeletal protein adducin 3 (ADD3).

To determine whether AML FAB subtypes could be distinguished by using our boutique array, we compared the gene expression profiles obtained for 11 patients diagnosed as AML-M1 and for 22 patients classified as AML-M2. However, we did not detect any statistically significant differences (adjusted p-value >0.35) (Table III). A separate analysis of AML-M1 vs. HV and AML-M2 vs. HV revealed 54 and 149 DEGs, respectively (α=0.05). Among the 54 DEGs from the comparison between AML-M1 and HV, 50 were shared with the list of 149 DEGs obtained from the comparison between AML-M2 and HV, thus suggesting that the AML-M1 subtype is a subset of AML-M2 and the latter shows some unique features. The list of genes shared between the M1 and M2 FAB subtypes included all the genes shown in Fig. 1B apart from SRGN, CD74, PTPRE and ADD3, which were unique for the comparison between AML-M2 and HV. Other genes shared between M1 and M2 were ERCC2, ITGB4, HSPA8, EEF1G, HOXB6, IGFBP7, KIT, PKM2, RASSF5, MN1, ARHGEF12, HSP90AB1, MYB, HSPA9 and GSN (overexpressed), and CCL19, MAP4K2, PF4, CCL5, B2M, HCK, CAPN10, HOXA9, TRADD, PGK1, PLBD1, IFNA1, BTG1, NFKBIA and PDGFRB (underexpressed). To 4 genes unique for AML-M1 belonged 2 underexpressed (ACTB and CST7) and 2 overexpressed (INPP5D and RHOH) genes. The list of overexpressed genes unique for AML-M2 included inter alia CTSG, CTSD, GAL, MYH11, GAPDH, ADRA2C, ATP6V0C, MAP7, ARF1, METTL3, BCR, TUBB, ABL1, NPM1, CD34, JUNB and MPO, encoding myeloperoxidase, marker protein of AML with maturation.

As regards the effects of mutations, we analyzed the gene expression profiles from the following AML subgroups: RUNX1/RUNX1T1⁺ [5 patients with translocation t(8;21), typical of the AML-M2 FAB subtype], NPM1⁺ (8 patients with an NPM1 mutation) and FLT3⁺ (4 patients with FLT3-ITD). Owing to the co-existence of the FLT3-ITD and NPM1 mutation in some patients, we were also able to distinguish 3 additional subgroups: NPM1⁺/FLT3⁺ (2 patients), NPM1⁺/FLT3⁻ (6 patients) and NPM1⁻/FLT3⁺ (2 patients). However, we were not able to identify any distinct clusters or detect significant differences in any of the comparisons (Table III).

The last result worth mentioning regarding the AML analysis using the AML-array is the classification of the patients according to their response to therapy. The association of gene expression with the prognosis and treatment outcome is a clinically important aspect of leukemia studies. Although all the studied samples were collected at the time of diagnosis, further monitoring of patients enabled to divide them into 3 general subgroups corresponding to their responses to therapy: Patients who reached complete remission (CR, 15 patients), patients resistant to therapy (RES, 9 patients) and patients who succumbed to the disease during therapy (X, 9 patients). The CR subgroup was further divided into 2 smaller subgroups: Patients with complete remission lasting up to 1 year (CR-short, 12 patients) and patients with complete remission lasting >1 year (CR-long, 3 patients). We did not detect any statistically significant differences in gene expression between any of the above subgroups (Table III).

In all analyses performed with the AML-array, significant differences were detected only when an AML subgroup was compared with HV. Notably, the obtained lists of DEGs, despite being of variable length, consisted a high number of shared genes, and STMN1 was always at the top of the ranking (Table III). To eliminate the HV impact, we repeated all the analyses, focusing solely on the AML samples. The result was the same: No statistically significant differences were detected between AML subgroups apart from one: The overexpression of the CPA3 gene, encoding a carboxypeptidase, in NPM1⁺/FLT3⁺ AML vs. remaining AML samples (p=0.043) (Table III). Unsupervised hierarchical clustering of the data also indicated that the patients with AML were generally a homogenous group, despite individual variations (Fig. 2).

To verify the results of the microarray analysis, we selected 4 genes, STMN1, NPM1, S100A9 and S100A8, and quantified their expression with the use of real-time PCR and a Rotor-Gene Q instrument. In this case, our sample set contained 14 HV and 21 AML samples. The STMN1 gene was selected as the most upregulated gene in AML compared with HV and NPM1 was selected as upregulated only in some AML subgroups (AML-M2, AML-NPM1⁺ and AML X), whereas S100A9 was selected as one of the most downregulated genes in AML and different AML subgroups compared with HV (Table III). S100A8 was underexpressed with statistical significance only in AML and in the AML-CR subgroup (Table III); however, we were interested whether 2 genes belonging to the same gene family present similar expression patterns.

First, owing to the variable expression of many so-called ‘housekeeping genes’, we tested 4 commonly used reference genes: ACTB, GAPDH, PGK1 and PPIA. The expression stability of each gene across the samples was estimated with geNorm. PGK1 and PPIA were indicated as the most stable, whereas GAPDH was the most variable, showing a trend toward upregulation in the AML samples (data not shown). This finding was consistent with the results of microarray analysis, in which GAPDH was in the list of DEGs between the AML and HV samples. Finally, 2 reference genes, PGK1 and PPIA, were used as normalization factors. Gene expression was analyzed in the context of AML subtypes, mutation status and response to therapy (Figs. 3 and 4).

The most consistent results between the microarray and real-time PCR experiments were obtained for the STMN1 gene, which was significantly induced in AML compared with HV (p=0.00044), although the fold change was somewhat lower (~4.5-fold compared with an almost 8-fold change in the microarray data analysis) (Fig. 3A). Similar to the microarray data analysis results, STMN1 upregulation was visible in each AML subgroup compared with HV (t-test p=0.038–0.000019) (Figs. 3 and 4). There were no significant differences in STMN1 expression between particular AML subgroups with the following two exceptions: NPM1⁻/FLT3⁻ vs. NPM1⁻/FLT3⁺ (p=0.015) and NPM1⁻/FLT3⁺ vs. NPM1⁺/FLT3⁻ (p=0.029). The highest STMN1 expression was observed in the NPM1⁻/FLT3⁺ subgroup (Fig. 4D).

As regards NPM1 gene expression, we observed a trend toward an upregulation in AML (Fig. 3A); however, the increase in NPM1 expression was statistically significant in only one comparison (RES vs. HV; t-test p=0.021) (Fig. 3C and D).

Real-time PCR-based analysis of S100A8 and S100A9 genes confirmed their general suppression trend in AML or at least some AML subgroups compared to HV. The expression levels of these 2 relative genes, coding for calcium-binding proteins, highly correlated (Pearson’s correlation, 1; p-value <2.2e⁻16); however, a statistically significant decrease was observed only for S100A9: In the AML-M1 subtype compared with HV (t-test p=0.0012), as well as with AML-M2 (t-test p=0.038) (Fig. 3B), in NPM1-mutated AML vs. HV (t-test p=0.013) (Fig. 4B), irrespective of the FLT3 mutation status (t-test p=0.048 for HV vs. NPM1⁺/FLT3⁻ and t-test p=0.01 for HV vs. NPM1⁺/FLT3⁺) (Fig. 4D), and in the RES subgroup vs. HV (t-test p=0.018) (Fig. 3C and D).

As the same group of AML samples was also previously analyzed in a proteomic analysis in our laboratory (33,34), in this study, we attempted to compare the results of the AML transcriptome and proteome analyses. However, owing to the limitations of both approaches, the comparison could be made for only a small subset of genes/proteins. The AML-array consisted of probes that were complementary for >900 genes, representing ~4% of all human genes. However, the 2-DE-MS analysis of AML proteomes enabled the detection of <200 proteins (33), representing much <1% of all the human proteins. In detail, from the 184 unique proteins detected in the proteomic study, only 52 (28%) had corresponding gene-specific probes on the AML-array. From the 184 identified proteins, 40 were determined to be differentially accumulated in the studied samples, of which 20 were represented on the AML-array: 13 by gene-specific probes and 7 by the probes complementary to the related genes (from the same gene family). Half of the results of the proteomics and transcriptomic analyses were concordant or partially concordant, but the other half was discordant (Table V).

To verify some of the discrepancies between the transcriptomic and proteomic results, we applied another quantitative PCR method, ddPCR, which is one of the most precise, sensitive and accurate PCR types. We expanded our research group to 20 HV samples and 41 AML samples. For the analysis, performed with a QX200 Droplet Digital PCR system (Bio-Rad), we selected 4 genes. Two of these, ANXA3, encoding Annexin 3, and CAT, encoding catalase, were determined as being differentially accumulated in the proteomics analyses. None of these genes were identified with the AML-array as being differentially expressed, and the expression trend was concordant for CAT, but was discordant for ANXA3. As a type of positive control, we used the ABL1 oncogene, which was slightly, but significantly, increased in AML-M2 and AML-RES compared with HV in the microarray experiment. Additionally, we analyzed Wilms tumor 1 gene (WT1), which has recently been described as overexpressed and prognostically relevant in most acute leukemias (39,40), but was not included in our AML-array probe set. As a reference gene, we selected PGK1, on the basis of our previous real-time PCR analysis. As regards the CAT gene, we were able to multiplex reactions with a reference gene, due to the differences in PCR product sizes (Table II). For the ABL1, ANXA3 and WT1 genes, a reference gene was analyzed in a separate reaction, but always within the same run. The data from duplicate experiments were separately processed in Quanta Soft version 1.5.38.1118 software (Bio-Rad). The relative expression value of each gene was calculated as the ratio of droplet numbers for the studied gene vs. the reference gene. The results of the statistical analysis, performed in R Bioconductor, are presented in Figs. 5 and 6.

Our positive control gene, ABL1, revealed a trend towards an upregulation in AML, statistically significant in FLT3- and NPM1-mutated AML vs. HV (HSD Tukey test p=0.02 for HV vs. FLT3⁺; HSD Tukey test p=0.0003 for HV vs. NPM1⁺/FLT3⁺) (Fig. 6C and D). ABL1 gene expression was also significantly upregulated in the FLT3⁺ samples when compared with the FLT3⁻ samples (HSD Tukey test p=0.02) (Fig. 6C) and significantly different in the following comparisons: NPM1⁻/FLT3⁻ vs. NPM1⁺/FLT3⁺ (HSD Tukey test p=0.0003), NPM1⁺/FLT3⁺ vs. NPM1⁺/FLT3⁻ (HSD Tukey test p=0.002) and NPM1⁺/FLT3⁺ vs. NPM1⁻/FLT3⁺ (HSD Tukey test p=0.0097). The ANXA3 gene was generally expressed at a very low level, with a trend towards overexpression in the AML vs. HV; however, the greatest increase in ANXA3 expression was observed in AML-M2 (Fig. 5B). The only statistically significant difference in the expression of this gene was observed between AML-M1 and AML-M2 (t-test p=0.023) (Fig. 5B), what definitively confirmed the proteomics results. Partially concordant results were obtained for the CAT gene. In agreement with the proteomics results, we observed a significant increase in CAT expression in AML and some AML subgroups compared with HV (Figs. 5 and 6, t-test p=0.0096 for HV vs. AML-M1; t-test p=0.022 for HV vs. AML-M2; t-test p=0.00085 for HV vs. RUNX1/RUNX1T1⁻; p=0.006 for HV vs. NPM1⁻; p=0.02 for HV vs. NPM1⁺; p=0.00085 for HV vs. FLT3⁻; t-test p=0.04 for HV vs. FLT3⁺; t-test p=0.008 for HV vs. NPM1⁻/FLT3⁻; t-test p=0.03 for HV vs. NPM1⁺/FLT3⁻). We also detected a significant increase in CAT expression in the FLT3⁻ compared to the FLT3⁺ samples (t-test p=0.023). However, the differences in the CAT gene expression levels between the M1 and M2 AML subtypes were subtle and not significant (t-test p=0.52), in contrast to the proteomics findings (Fig. 5B). As regards the response to therapy, the only statistically significant difference, and the most elevated CAT expression was observed in the CR and CR-long AML subgroups (t-test p=0.0044 for HV vs. CR and t-test p=0.013 for HV vs. CR-long) (Fig. 5C and D). As regards the WT1 gene, the ddPCR experiment revealed an evident and statistically significant upregulation of this gene in all AML subgroups, whereas in the HV samples, its expression was scarcely detectable (t-test p=0.046–0.000066) (Figs. 5 and 6). An interesting observation was the increasing trend in WT1 expression in the therapy-related AML subgroups: The higher the WT1 level, the less responsive was AML to therapy (Fig. 5C and D). High levels of WT1 were also detected in the AML-M1 (Fig. 5B) and AML samples without the RUNX1/RUNX1T1 mutation (Fig. 6A), but harboring NPM1, FLT3 or both of these mutations (Fig. 6B–D). The differences across the particular AML subgroups were not significant, apart from 2 comparisons: AML-M1 vs. AML-M2, where WT1 was upregulated in AML-M1 (HSD Tukey’s test p= 0.039) (Fig. 5B) and RUNX1/RUNX1T1⁺ vs. RUNX1/RUNX1T1⁻, where WT1 was upregulated in AML without translocation t(8;21) (t-test p=0.0023) (Fig. 6A).

In order to identify links between gene expression measured by RT-qPCR, the mutation status and the clinical characteristics of the patients with AML (WBC, age and sex), we performed a series of statistical tests (t-test, one- and two-way ANOVA, Tukey multiple comparison of means and Pearson’s correlation). The strongest correlations were observed between the WBC counts and the AML FAB subtypes (one-way ANOVA, p=0.0101, WBC counts were higher in AML-M1), WBC counts and NPM1 mutation status (one-way ANOVA, p=0.0229), WBC counts and FLT3 mutation status (one-way ANOVA, p=0.0308), and WBC counts and both FLT3 and NPM1 mutation status (two-way ANOVA, p=0.0160), and the highest WBC count was observed in mutated AML samples.

As regards the effects of age, we observed a moderately positive correlation between patient age and NPM1 gene expression (Pearson’s correlation, 0.348; p=0.1222), a moderately negative correlation between age and the expression of the 2 genes, S100A9 (Pearson’s correlation, −0.380; p=0.0899) and ANXA3 (Pearson’s correlation, −0.33; p=0.046) and a correlation between age and the mutation status of NPM1 (one-way ANOVA, p=0.0411) or NPM1 and FLT3 together (two-way ANOVA, p=0.0446). There were no correlations between the RT-qPCR-determined expression of any other gene and age (Pearson’s correlations equal to −0.23, 0.019, 0.15, 0.07 and 0.2 for S100A8, STMN1, ABL1, CAT and WT1, respectively). We also did not find any correlations between the gene expression levels and sex (one-way ANOVA, p-values equal to 0.45, 0.659, 0.286, 0.555, 0.633, 0.302, 0.628 and 0.263 for NPM1, S100A8, S100A9, STMN1, ABL1, ANXA3, CAT and WT1, respectively) or WBC count (Pearson’s correlations equal to −0.02, −0.25, −0.27, −0.02, 0.07, −0.19, −0.18 and 0.28 for NPM1, S100A8, S100A9, STMN1, ABL1, ANXA3, CAT and WT1, respectively).

Of note, there was no correlation between NPM1 gene expression and the NPM1 mutation status (one-way ANOVA; p=0.634). The presence of the translocation t(8;21) and fusion gene RUNX1/RUNX1T1 did not correlate with any molecular or clinical characteristic.

To compare our results with other AML gene expression profiling results, we explored the Gene Expression Omnibus (GEO) database (www.ncbi.nlm.nih.gov/sites/GDSbrowser) to identify the most relevant datasets. We selected 7 datasets, including 6 derived from Affymetrix Gene Chips, the most commonly used microarray platform, and one from Agilent microarrays. Among these datasets, only 3 (GDS3057, GDS1059 and GDS2251) contained healthy control samples. These 3 datasets were screened with 238 names of genes, which we had identified as differentially expressed between AML (or particular AML subgroups) and HV by using the AML-array, RT-qPCR and proteomics. Table VI shows the results of comparative analysis. The percentage of genes with an expression trend concordant in our data and any GEO dataset, was in the range 47.8–51.9% (mean, 49.2%) or in the range 51.8–78.3% (mean, 65.3%), when including partially concordant genes, for which expression trends were concordant for some (but not all) conditions/cell lines/microarray probes. Interestingly, the level of coincidence between GEO datasets was much lower: The percentage of concordant genes remained in the range 35.3–43.4% (mean, 39.3%) or in the range 47.4–64.2% (mean, 56.6%), when including partially concordant genes. In reference to discordant genes, an opposite effect could be observed, although the difference was smaller: The percentage of discordant genes in comparative analysis between our data and any GEO data set was in the range 15.5–31.3% (mean, 22.8%) while within GEO datasets in the range 10.0–27.1% (mean, 19.3%). Analyzing 238 gene expression trends, we were able to select 45 strong candidate genes with expression trends concordant in our dataset and all 3 GEO datasets. In total, 26 genes (ABCF1, ABL1, AK2, ANGPT1, BCR, CD34, CDK6, CEBPA, HOXA10, HSP90AB1, HSPA8, HSPA9, IGFBP7, MAP7, MCM5, MT1H, MYB, NPM1, RBBP4, RPLP0, SET, SLC25A1, SMAD3, STMN1, TRAF4 and WT1) were consequently increased in AML when compared to the normal control, whereas 19 were decreased in AML (BTG1, CEBPB, CSF3R, CX3CR1, CYBB, FCER1G, FCN1, HCK, IGHM, LCP1, LGALS3, LYZ, PF4, S100A8, S100A9, STAT5B, STK4, TMSB4X and TRADD). The above-mentioned list contained 14 (10 increased and 4 decreased) from the 23 genes presented in Fig. 1B, identified as the genes with the highest expression fold change in AML-array-based analysis, and one gene (RBBP4) identified in our data only by proteomic approach. Not all of these genes were up-to-date described as AML biomarkers. It is also worth mentioning that 2 members of the CCAAT/enhancer binding protein family of transcription factors, CEBPA and CEBPB, exhibited opposite expression trends in all studied datasets.

As regards FAB subtypes, we analyzed the expression trends of 6 genes, including 5 genes encoding proteins determined as differentially accumulated between AML-M1 and AML-M2 by proteomic approach (ANXA3, CAT, LCP1, PGD and PRDX6) (33) and one gene identified with ddPCR (WT1). In comparison to the GDS1064 dataset, referring to AML in different stages of maturation, we observed concordant results for only 2 genes (33%), WT1 and ANXA3, the only 1 gene out of 5 selected by proteomic approach.