Coagulation FXIII-A Protein Expression Defines Three Novel Sub-populations in Pediatric B-Cell Progenitor Acute Lymphoblastic Leukemia Characterized by Distinct Gene Expression Signatures

Background: Leukemic B-cell precursor (BCP) lymphoblasts were identified as a novel expression site for coagulation factor XIII subunit A (FXIII-A). Flow cytometry (FC) revealed three distinct expression patterns, i.e., FXIII-A negative, FXIII-A dim, and FXIII-A bright subgroups. The FXIII-A negative subgroup was significantly associated with the “B-other” genetic category and had an unfavorable disease outcome. Methods: RNA was extracted from bone marrow lymphoblasts of 42 pediatric patients with BCP-acute lymphoblastic leukemia (ALL). FXIII-A expression was determined by multiparameter FC. Genetic diagnosis was based on conventional cytogenetic method and fluorescence in situ hybridization. Affymetrix GeneChip Human Primeview array was used to analyze global expression pattern of 28,869 well-annotated genes. Microarray data were analyzed by Genespring GX14.9.1 software. Gene Ontology analysis was performed using Cytoscape 3.4.0 software with ClueGO application. Selected differentially expressed genes were validated by RT-Q-PCR. Results: We demonstrated, for the first time, the general expression of F13A1 gene in pediatric BCP-ALL samples. The intensity of F13A1 expression corresponded to the FXIII-A protein expression subgroups which defined three characteristic and distinct gene expression signatures detected by Affymetrix oligonucleotide microarrays. Relative gene expression intensity of ANGPTL2, EHMT1 FOXO1, HAP1, NUCKS1, NUP43, PIK3CG, RAPGEF5, SEMA6A, SPIN1, TRH, and WASF2 followed the pattern of change in the intensity of the expression of the F13A1 gene. Common enhancer elements of these genes revealed by in silico analysis suggest that common transcription factors may regulate the expression of these genes in a similar fashion. PLAC8 was downregulated in the FXIII-A bright subgroup. Gene expression signature of the FXIII-A negative subgroup showed an overlap with the signature of “B-other” samples. DFFA, GIGYF1, GIGYF2, and INTS3 were upregulated and CD3G was downregulated in the “B-other” subgroup. Validated genes proved biologically and clinically relevant. We described differential expression of genes not shown previously to be associated with pediatric BCP-ALL. Conclusions: Gene expression signature according to FXIII-A protein expression status defined three novel subgroups of pediatric BCP-ALL. Multiparameter FC appears to be an easy-to-use and affordable method to help in selecting FXIII-A negative patients who require a more elaborate and expensive molecular genetic investigation to design precision treatment.


INTRODUCTION
According to current knowledge, acute lymphoblastic leukemia (ALL) can be best characterized by an integrated set of clinical, pathological, morphologic, immunophenotypic, and genetic properties. Detailed characterization of the neoplastic clones of children with ALL allowed the development of highly effective risk-tailored therapies. Introduction of advanced diagnostic tools helped in identifying an increasing number of molecular targets and may contribute to the application of personalized treatment offering cure for each individual patient. Recurrent genetic aberrations represent the basis of the classification of ALL according to the 4th edition of the WHO Classification of Tumors of Hematopoietic and Lymphoid Tissues in 2008, and its 2016 revision (1,2). Gene expression profiling by oligonucleotide microarray was shown to contribute to conventional and molecular cytogenetics by improving diagnostic accuracy and prognostic relevance as well as by defining new entities. The "Bother" genetic subgroup, i.e., BCP-ALL cases without established recurrent genetic abnormalities exhibits alterations to be revealed by advanced genomic technologies. Within "B-other" ALL, Philadelphia-like (Ph-like) or BCR-ABL1-like ALL, a provisional entity of the 2016 revision of the WHO classification, has been defined based on gene expression signature similar to Phpositive/BCR-ABL1-positive ALL and characteristically distinct from the rest of BCP-ALL cases ("non-B-other") (3)(4)(5)(6). Recently, fourteen B-cell precursor (BCP)-ALL subgroups were defined by analyzing a large sample set of an international study using RNAsequencing. The study revealed six new subgroups in addition to previously identified ones demonstrating the diagnostic utility and prognostic power of gene expression studies in pediatric BCP-ALL (7).
Our group identified BCP-ALL blasts as a new expression site for coagulation factor XIII subunit A (FXIII-A). In contrast to the BCP-ALL blasts, neither normal bone marrow B-cell progenitors, nor mature B-cells, both normal and leukemic once, nor Tlymphoblasts/lymphocytes express FXIII-A (8). We observed three characteristic expression patterns by flow cytometry (FC): FXIII-A negative, FXIII-A dim, and FXIII-A bright lymphoblasts with about two thirds of pediatric BCP-ALL patients representing FXIII-A positive (FXIII-A dim and bright) cases. According to our retrospective clinical investigation, patients with FXIII-A negative lymphoblasts had significantly worse event-free and overall survival than patients with FXIII-A positive lymphoblasts. Moreover, the "B-other" genetic subtype was significantly more frequent within the FXIII-A negative than in the FXIII-A positive subgroup (9). The prognostic importance of FXIII-A expression by FC was confirmed and further specified in a recently concluded prospective study suggesting a favorable prognostic effect of the FXIII-A dim expression pattern within the FXIII-A positive subgroup (unpublished results). Here we present results, for the first time, on gene expression signatures associated with the three characteristic FXIII-A protein expression patterns. In addition, we investigated overlaps between the gene expression profile of the FXIII-A negative protein expression subgroup with the gene expression profile of the "B-other" subgroup.  (9,22) ] were excluded, since they were not treated according to BFM ALL-IC 2009 protocol. Quantitatively and qualitatively suitable RNA was separated from 42 patients with FXIII-A negative (14), FXIII-A dim (21), and FXIII-A bright (7) BCP-ALL.

Patients
Excess bone marrow (BM) samples were collected after signed informed consent obtained from study participants and their legal caregivers.

RNA Preparation
Excess BM samples from BCP-ALL patients were collected into PAXgene Blood RNA Plastic Tube (PreAnaltyX, Hombrechtikon, Switzerland). Total cellular RNA was isolated using PAXgene Blood miRNA Kit (PreAnalityX) according to conventional protocol. Qualitative and quantitative analyses of RNA were performed using Agilent Bioanalyzer (Agilent Technologies, La Jolla, CA, USA). Samples with an RNA Integrity number above 8.0 were used in MicroArray analysis.

Microarray Analysis
Affymetrix GeneChip Human Primeview array (Affymetrix, Santa Clara, CA, USA) was used to analyze global expression pattern of 28,869 well-annotated genes. 3'IVT Expression Kit (Affymetrix) and GeneChip WT Terminal Labeling and Control Kit (Affymetrix) were used for amplifying and labeling 250 ng of RNA samples. Samples were hybridized at 45 • C for 16 h and then standard washing protocol was performed using GeneChip Fluidics Station 450 (Affymetrix) and the arrays were scanned on GeneChip Scanner 7G (Affymetrix) procedure. Data of this study have been deposited in NCBI's Gene Expression Omnibus and are accessible through GEO Series accession number GSE134480 (https://www.ncbi.nlm.nih.gov/geo/query/ acc.cgi?acc=GSE134480).

Gene Ontology Analysis
Gene Ontology (GO) analysis was performed using Cytoscape 3.4.0 software (cytoscape.org) with ClueGO application. The settings were the following: GO biological process and GO immune system process. For statistical analysis two-sided hypergeometric test and Benjamini-Hochberg FDR were used. Significantly enriched GO categories were considered to p-value <0.05 and κ score <0.4.

Real Time Quantitative Polymerase Chain Reaction (RT-Q-PCR) Validation of Microarray Data
For RT-Q-PCR validation of microarray data pre-designed, factory-loaded 384-well TaqMan low-density array (ThermoFisher Scientific, Waltham, MA, USA) was used to determine the level of expression of selected genes using technical duplicates. Genes for validation were selected based on fold-changes (with fold-changes >2.0 cut-off) between any of the pre-defined subgroups and based on putative biological and clinical relevance of gene functions as defined by GO analysis data ( Table 1). RT-Q-PCR expression levels of target genes were normalized to the mean of B2M, GAPDH, and GUSB reference genes. Normalized gene expression values were calculated based on the C t method, where relative expression equals 2 − Ct , where C t represents the threshold cycle (C t ) of the target minus that of the mean of reference genes.

In silico Investigation of Validated DE Genes
Interactions of validated genes and F13A1 gene were investigated using STRING v11. (12) and GeneHancer (13) databases. STRING v11 database contains putative protein-protein interactions predicted on a well-defined score system. GeneHancer portrays 285 000 integrated candidate enhancers and subsequently links enhancers to genes.

Statistical Analysis
Microarray data were analyzed by Genespring GX14.9.1 software (Agilent Technologies, La Jolla, CA, USA). To identify statistically significant genes, we used volcano plot analysis. The resulting scatterplot showed statistical significance (p-value) vs. magnitude of change (fold-change).  Affymetrix data files were imported using the Robust Multi-Array Average algorithm and median normalization was performed. To identify differentially expressed genes between pre-defined data sets, statistical analysis was performed using ANOVA with Tukey post-hoc test (14) and moderated Ttest, Benjamini-Hochberg False Discovery Rate was used for multiple testing corrections, p-value <0.05 was considered as significant difference.

Characterization of Samples of Children
With BCP-ALL RNA samples obtained from 42 pediatric patients with BCP-ALL were used for gene expression profiling investigations ( Table 2). Cytoplasmic expression pattern of FXIII-A by FC successfully divided patients into three categories: patients with

FXIII-A negative, FXIII-A dim and FXIII-A bright BCP-ALL.
Representative FC dot plots and histograms are shown in Figure 1. Regarding genetic categories, 27 patients had recurrent genetic abnormalities and 15 patients were assigned to the "Bother" subgroup: 7/12 FXIII-A negative patients, 6/21 FXIII-A dim patients, and 2/7 FXIII-A bright patients.  (Figure 2). Comparing the gene expression signature of the "B-other" subgroup with the rest of patient samples, the so called "non-B-other" group, 142 DE genes were found after filtering to 1.5-fold-change. Heat map analysis clearly showed two distinct sub-populations: gene expression signature of the "B-other" subgroup differed characteristically from the gene expression signature of the "non-B-other" subgroup ( Figure 3). Then we investigated how DE genes, characteristic for "B-other" status, were related to DE genes of subgroups according to FXIII-A expression pattern. We found 32 DE genes expressing exclusively in the FXIII-A negative and "B-other" group vs. FXIII-A negative and the "non-B-other" group. Heat map analysis confirmed that, with the exception of one outlier, gene expression profile characterizing FXIII-A negative samples overlapped with gene expression profile characterizing samples of the "B-other" subgroup (Figure 4).

Functional Characterization of Differentially Expressed Genes in the BCP-ALL Subgroups
Identification of enriched functional categories according to the FXIII-A protein expression pattern, DE genes were categorized into 156 GO processes. Nevertheless, most of them, in particular those with the strongest statistical p-values, were related to epigenetic and/or gene expression regulatory processes such as histone modification, chromatin organization, RNA destabilization, post-transcriptional regulation of gene expression, etc., or other regulatory and cellular processes, such  Table 1). In addition, we identified biological processes resulting in peptidyllysine modification, which are in relation with the known physiological function of FXIII-A catalyzing the formation of γglutamyl-ε-lysyl amide crosslinks between fibrin monomers to form an insoluble clot (Supplementary Table 1).

Validation of Global Transcriptomics Data
From the oligonucleotide microarray results of DE genes, either according to FXIII-A expression status or according to "B-other" genetic status we selected 45 genes for validation by RT-Q-PCR. Selection of 13/45 genes was based on fold change results, whereas an additional 32/45 genes were selected according to enriched functional categories of potential interest as defined by the GO analysis ( Table 1). We were not able to detect transcripts of RORA by RT-Q-PCR which might have a technical reason.

FXIII-A Expression-Based Results
Expression of F13A1 gene was detected and readily validated by RT-Q-PCR in every sample. Intensity of gene expression; however, was characteristically different among samples of the three different FXIII-A protein expression subgroups with an increasing intensity in terms of relative fold-changes measured by RT-Q-PCR from the FXIII-A negative, through dim to bright subgroups ( Figure 5).
Similarly, most of the genes (8/13 p < 0.05, 9/13 p < 0.10) from the group selected on the basis of highest fold changes between any two groups among FXIII-A negative, dim and bright groups according to the microarray results were validated by RT-Q-PCR. Of the 25 genes selected on basis of functional significance according to GO analysis, four genes could be validated. Within the GO group of translation regulation there were 2/6 with p < 0.05, and 3/6 with p < 0.10 ( Table 3) genes that could be validated, providing a considerably better ratio than it was found in the subgroups of apoptosis regulation (0/3), chromatin modification (1/5), and nuclear export (0/5) (Tables 1, 3). Generally, foldchanges were enhanced in RT-Q-PCR compared to microarray data (Supplementary Table 3). PLAC8 and HAP1 represent impressive examples. DE of PLAC8, a placenta-specific gene of unrevealed function, was not found significant by microarray but turned out significant upon validation. HAP1, a member of the translation regulation GO group, exhibited a very slight overexpression by microarray and it proved much stronger by RT-Q-PCR.
Considering individual normalized gene expression values, in most of the cases, there was a clear trend of a continuous increase from FXIII-A negative through dim to bright subgroups that was endogenously validated by the F13A1 relative expression. ANGPTL2, NUCKS1, RAPGEF5, and SEMA6A followed this trend (Figure 5). Relative fold-changes were similar whether determined by microarray measurements or RT-Q-PCR. In case of FOXO1, HAP1, and TRH, separation of the bright subgroup from the other two subgroups seemed more prominent and showed a relatively higher fold-change value as determined by RT-Q-PCR than compared to the microarray results. With the exception of PLAC8, validated DE genes were downregulated in the FXIII-A negative subgroup compared to the two other subgroups (Supplementary Table 3).
Potential interactions could not be revealed between the protein products of validated genes and FXIII-A using STRING v11 functional protein association networks database. Nevertheless, we could identify a FXIII-A dependent expression of FOXO1. Using GeneHancer database we could in silico identify enhancers that might have roles in the parallel upregulation of F13A1 and validated genes. Transcription factor binding sites of ATF7, POLR2A, RAD21, SMARCA5 gene products were identified for all of the 14 genes shown in Table 3.

"B-OTHER" STATUS-BASED RESULTS
For the "B-other" status-based validation, DE genes related to biological processes (macrophage migration, lymphocyte apoptotic process, T cell apoptotic process, and their regulation), were selected for validation, as GO analysis results were less diverse than in the case of FXIII-A expression-based GO annotations (Supplementary Table 1). We were able to validate the differential expression of five genes by RT-Q-PCR. DFFA, GIGYF1, GIGYF2, and INTS3 were overexpressed in "B-other" vs. "non-B-other" samples and in turn, CD3G exhibited a relatively lower expression in the "B-other" vs. "non-B-other" samples ( Figure 6). RORA, IL7R, CCL5, PLAC8, CX3CR genes, selected for validation based on their functions and fold-changes detected by microarray, could not be validated.

DISCUSSION
FXIII-A is a useful marker of acute myeloid leukemia (AML) of the monocyte and megakaryocyte lineages (31)(32)(33)(34). However, its expression in BCP-ALL blasts was an unexpected finding, since, in contrast to the normal myeloid counterparts of AML blasts, normal B-lymphocytes and their precursors do not express FXIII-A (8). Intracellular expression of FXIII-A protein can be consistently demonstrated in about two thirds of pediatric BCP-ALL samples by FC. Interestingly, the FXIII-A negative status was shown to be significantly associated with the "B-other" genetic subtype and patients with FXIII-A negative BCP-ALL had a significantly worse disease outcome than patients with FXIII-A positive lymphoblasts (9). These facts together suggested that pathological expression of FXIII-A in leukemic BCP blasts may define one or more sub-populations according to the FXIII-A expression pattern by FC. There was a continuous increase in normalized gene expression levels from FXIII-A negative through dim to bright subgroups that was endogenously validated by the F13A1 differential expression within the three FXIII-A protein expression groups. ANGPTL2, NUCKS1, RAPGEF5, and SEMA6A followed this trend. Based on the intensity of the differential expression, separation of FOXO1, HAP1, and TRH genes of the FXIII-A bright subgroup were more prominent. PLAC8 expression was most intensive in the FXIII-A dim subgroup.
Little is known on the transcriptional regulation of the F13A1 gene. Molecular investigations of myeloid leukemia cell lines revealed the presence of binding sites for ubiquitous (NF-1 and SP-1) and myeloid enriched (MZF-1-like protein, GATA-1, and Ets-1) transcription factors in the promoter region of the gene (35). Interestingly, transcription factor binding sites for SP-1, GATA-1, and Ets-1 could not be confirmed by the ENCODE ChIP-seq datasets (data not shown) (36). A more recent study demonstrated that FXIII-A was synergistically regulated by IL-4 and dexamethasone in alternatively activated macrophages and these regulatory molecules are relevant for both normal and leukemic BCPs as well (37). However, transcriptional regulation of the F13A1 gene in leukemic lymphoblasts has not yet been studied. Therefore, we decided to investigate the gene expression profile of BCP lymphoblasts according to their FXIII-A expression status. Preliminary results of an ongoing prospective, multi-centric clinical study performed by our group suggested the clinical relevance of the sub-populations characterized by the FXIII-A expression pattern, since patients with FXIII-A dim but not bright lymphoblasts had a better disease outcome than patients with FXIII-A negative lymphoblasts (unpublished preliminary results). Accordingly, we investigated the global gene expression signature of three sample sets.
We have revealed that the three groups according to FXIII-A expression pattern had unique gene expression signatures which were characteristically different from each other. GO analysis of data resulted in a number of biologically relevant enriched functional categories. Identification of biological processes resulting in peptidyl-lysine modification supported the clinical relevance of our findings. FXIII-A is an enzyme catalyzing the formation of γ-glutamyl-ε-lysyl amide crosslinks. In addition to its well-known role in the formation of stable fibrin clot, FXIII-A participates also in much less characterized intracellular regulatory processes including differentiation of monocyte/macrophages, osteoblast, and osteoclast (38). In platelets, FXIII-A has been shown to crosslink cytoskeletal proteins (39). Should FXIII-A participate in similar processes in leukemic BCP lymphoblast, intracellular FXIII-A activity might contribute to autophagy and apoptosis of FXIII-A expressing lymphoblast. Moreover, the product of four genes validated in our cohort, EHMT1, ING5, MDM2, NUP43 were shown to methylate and acetylate lysyl groups of intranuclear and intracellular proteins, among others, p53, and that way they can regulate survival of neoplastic cells (25).
In addition, we showed that gene expression pattern of the FXIII-A negative subgroup overlapped with the gene expression profile of the "B-other" genetic subgroup. However, our results indicated that the FXIII-A expression status was a more powerful determinant of gene expression profile than the "B-other" status.
Fourteen genes were validated according to the FXIII-A expression status. Importantly, we were able to detect the  Table 3). NUP43 has not yet been shown to be associated with any forms of human cancer in contrast to other members of the NUP gene family (40). PLAC8 which is a trophoblast lineage marker physiologically, was most intensively expressed in the FXIII-A dim subgroup. This gene has been shown to be aberrantly activated in various types of cancer arising in mammals and mammalian cancer cell lines, but not in any subtype of human ALL (24,25). Based on in silico investigations we were not able to reveal a direct link between DE genes as related to the three different FXIII-A expression groups and the regulation of F13A1 gene. Common enhancer elements of the validated DE genes and the F13A1 make likely that common transcription factors may regulate the expression of these genes in a similar fashion. Validated DE genes of the "B-other" subgroup overlapping with the gene expression profile of the FXIII-A negative subgroup suggested the definition of a different sub-population from the BCR-ABL1-like/Ph-like B-ALL subgroup discovered originally by oligonucleotide microarrays (4,5). Evaluating the microarray results, we identified 14 DE genes in our "Bother" subgroup (ANXA5, ARL4C, CD34, IFNGR1, MAGT1, MAPKBP1, MAP3K2, ME3, OSBPL8, PAPOLA, PTPRC, SUZ12, TACC1, TEMPO) overlapping with DE genes in the "B-other" subgroup defined by Den Boer et al. (5). However, all the DE genes within our "B-other" group differed from class-defining genes of the so called "novel" ALL subtype introduced by Yeoh et al. (3) implying that they specify a different subtype. Each validated gene in this category, i.e., CD3G, DFFA, GIGYF1, GIGYF2, and INTS3 were shown to play a role in neoplastic processes, including leukemia (41)(42)(43)(44)(45). However, none of these genes were published previously to be associated with pediatric BCP-ALL.
The small sample size and the low number of validated genes due to restricted budget represent an important limitation of the present study. The small number of patients involved in this investigation would make a clinical outcome analysis unreliable, even if statistically significant. Nevertheless, there was a substantially larger ratio of children who died due to their disease in the FXIII-A negative group: 5/14 children (36%), than in the FXIII-A dim group: 5/21 (24%) and in the FXIII-A bright group: 1/7 (14%). To investigate the clinical significance of FXIII-A expression in children with BCP-ALL we have started a pilot study of the ongoing BFM ALL-IC 2009 clinical trial with the participation of the Hungarian, Polish and Slovak national study groups. Preliminary analysis of this pilot study confirmed the unfavorable outcome of patients with FXIII-A negative BCP-ALL (results to be published upon completion of the pilot study). Similarly, an exact relationship between gene expression signature of FXIII-A negative samples and "Bother" samples, including BCR-ABL1-like/Ph-like signature, can be defined by investigating a larger number of samples from children with BCP-ALL.
In conclusion, we were the first to demonstrate the general expression of F13A1 gene in pediatric BCP-ALL samples. The intensity of F13A1 expression corresponded to the expression of FXIII-A protein, determined by FC in these samples. Three well-defined categories of FXIII-A protein expression: FXIII-A negative, FXIII-A dim, and FXIII-A bright subgroups defined characteristic and distinct gene expression signatures detected by Affymetrix oligonucleotide microarrays. Gene expression signature of the FXIII-A negative subgroup showed an overlap with that of the subgroup with "B-other" genetics. Validated genes proved biologically and clinically relevant. We described differential expression of genes not shown previously to be associated with pediatric BCP-ALL. Protein products of the newly identified genes may offer therapeutic targets for precision treatment of BCP-ALL. Multiparameter FC appears to be an easy-to-use and affordable method to assist in selecting pediatric patients with FXIII-A negative BCP-ALL who require a more elaborate and expensive molecular genetic investigation to design individualized therapeutic protocols.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in the NCBI's Gene Expression Omnibus and are accessible through GEO Series accession number GSE134480 (https://www.ncbi.nlm.nih. gov/geo/query/acc.cgi?acc=GSE134480).

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Scientific Research Ethical Committee (TUKEB) of the Medical Research Council of Hungary: No. 43033-1/2014/EUK(423/2014). Written informed consent to participate in this study was provided by the participants' legal guardian/next of kin.