Identification of predictive genetic signatures of Cytarabine responsiveness using a 3D acute myeloid leukaemia model

Abstract This study reports the establishment of a bone marrow mononuclear cell (BMMC) 3D culture model and the application of this model to define sensitivity and resistance biomarkers of acute myeloid leukaemia (AML) patient bone marrow samples in response to Cytarabine (Ara‐C) treatment. By mimicking physiological bone marrow microenvironment, the growth conditions were optimized by using frozen BMMCs derived from healthy donors. Healthy BMMCs are capable of differentiating into major hematopoietic lineages and various types of stromal cells in this platform. Cryopreserved BMMC samples from 49 AML patients were characterized for ex vivo growth and sensitivity to Ara‐C. RNA sequencing was performed for 3D and 2D cultures to determine differential gene expression patterns. Specific genetic mutations and/or gene expression signatures associated with the ability of the ex vivo expansion and response to Ara‐C were elucidated by whole‐exome and RNA sequencing. Data analysis identified unique gene expression signatures and novel genetic mutations associated with sensitivity to Ara‐C treatment of proliferating AML specimens and can be used as predictive therapeutic biomarkers to determine the optimal treatment regimens. Furthermore, these data demonstrate the translational value of this ex vivo platform which should be widely applicable to evaluate other therapies in AML.

different partners and alter the expression of genes involved in the development of AML. 3 Acute myeloid leukaemia classifications by recent WHO criteria are focused on significant cytogenetic and molecular genetic subgroups. 5 Chemotherapy has remained as a main treatment for AML for the past several decades. Breakthrough new AML therapies of two targeted drugs, Midostaurin and Enasidenib, were approved by the US Food and Drug Administration (FDA) in 2017. Midostaurin targets AML harbouring fms like tyrosine kinase 3 (FLT3) mutations, whereas enasidenib treats AML with an isocitrate dehydrogenase (IDH) 2 mutations. Most recently Venetoclax, a Bcl-2 inhibitor, has also been approved by the FDA to treat AML. Selection of appropriate patients for clinical trials has been a challenging task, and researchers have been constantly searching for genetic biomarkers and gene signatures that can be used to predict drug responsiveness by whole-exome sequencing (WES) and RNAseq technologies. 4,6 Drug responsiveness was usually evaluated in 2D ex vivo culture during short-term treatment of 2-4 days using viability assay, which is considerably shorter than the duration of chemotherapy in the clinic that lasts ten days per cycle. 1 The 2D culture has been recognized for a long time as an insufficient system to accurately predict drug responsiveness due to the lack of physiologically relevant microenvironment. 7,8 To overcome this limitation, 3D culture models have been developed in recent years for various cell types and tissues. 7,8 The 3D ex vivo models can recapitulate the complex physiology and maintain functional in vivo responses that are not observed in routine 2D cultures. In the case of 2D culture of leukaemic BMMCs, lack of bone marrow microenvironment prevents long-term culture of leukaemic cells without supplementing maintenance media with growth factors or co-culture with stromal cells. However, co-culture with stromal cells further complicates drug testing because specific killing of AML needs to be carefully determined, and often stromal cell growth can become dominant unless irradiated stromal cells are used. In recent years, several 3D models have been reported for culturing bone marrow samples, including matrigel, hydrogel, synthetic polymers and human amniotic membrane with or without co-culture of stromal cells. [9][10][11][12] The stiffness of matrix also plays an important role in AML cell growth and responsiveness to drug treatment. Shin

and Mooney
reported that matrix mechanics influenced both cell proliferation and sensitivity to chemotherapeutic agents of several AML cell lines ex vivo and in vivo using a hydrogel system. 13 However, it remains to be investigated whether any of these 3D systems can be used to identify gene signatures and mutations predictive of responsiveness to drug treatment of BMMCs from AML patients.
In the current study, a 3D platform was established and evaluated to culture BMMCs from both normal donors and AML patients. Under optimized conditions, 3D bone marrow cultures properly maintained growth and showed differentiation of major hematopoietic cell lineages as well as stromal cells. Gene expression signatures of 2D and 3D cultures were evaluated by comparing to the cryopreserved BMMCs using RNAseq. This platform was then used to identify AML responders and non-responders to Ara-C, a first-line chemotherapy drug for treating AML. 1 Eight Ara-C responders and five non-responders were assessed for differential gene expression patterns in vehicle vs Ara-C treatment.
Furthermore, novel genetic mutations associated with Ara-C responsiveness were revealed through WES. These results indicate that 3D ex vivo translational platforms can assist in identifying prospective biomarkers, accelerate discovery of novel treatments as well as to define AML treatment regimens in the clinical setting prior to initiation of therapy.

| Tartrate-resistant acid phosphatase (TRAP) staining
Cells were washed three times with PBS and fixed with 4% paraformaldehyde at room temperature without permeabilization.
Following three more washes with PBS, chromogenic substrate was added according to manufacturer's instruction. Cells were incubated for 20-60 minutes at room temperature and visualized at 10× magnification to determine best colour development timing for imaging.

| Alkaline phosphatase (AP) staining
Cells were washed three times with PBS and fixed for 1 minute with 4% formaldehyde. Following another three washes with PBS /0.05% tween-20, 20 μL of substrate (one BCIP/NBT tablet

| Adipocyte staining
Cells were washed three times with PBS, followed by fixation with 4% formaldehyde and permeabilization in 4% formaldehyde/0.4% Triton X-100. After three washes with PBS, cells were stained with 5 μg/mL Hoechst and 2.5 μg/mL deep red cell mask dyes for 1 hour at room temperature. Cells were washed three more times with PBS, and lipid tox was added to the plate. After shaking at room temperature for 30 minutes, cells were imaged at 10× magnification.

| Mineralization staining
Cells were washed three times with PBS, followed by fixation with 4% formaldehyde for 30 minutes. After three washes with distilled water, cells were stained with Alizarin Red Stain solution for 45 minutes in dark with gentle shaking. After additional four washes with distilled water, cells were stored in PBS and imaged at 10× magnification.     Figure S1).  To further investigate the differences between 2D and 3D cultures, gene expression signatures were compared by RNAseq using BMMCs from AML donor 3 after 10-day culture. 2D and 3D signatures were also compared to uncultured cryopreserved samples from the same donor. Results of this study indicate that, the gene expression patterns of 3D and 2D cultures are more similar to each other than to the uncultured frozen samples ( Figure 2B). A total of 4284 and 3203 genes were differentially regulated in 3D vs 2D system, respectively, when using uncultured samples as a reference ( Figure 2C). When looking at the union signature (5090 genes), the majority of differentially expressed genes as compared to the cryopreserved sample are similar in both 2D and 3D systems ( Figure 2D). However, a cluster of 461 genes was significantly distinct in 3D platform as compared to the 2D system, indicating existence of additional functionalities in the 3D environment ( Figure 2E). Several Gene Ontology (GO) terms were over-

| Determination of the Ara-C response in 49 AML patients
The newly defined 3D platform was then used to culture BMMCs from 49 AML donors (Table 1) Table 2. IC50 values ranged from 17 nmol/L to higher than 10 μmol/L. Ara-C responsiveness in 2D vs 3D culture was also compared using nine donors, and IC50 value shifts of more than 10-fold were observed for three donors (Table 2) with the relative IC50 ranking being different in 2D vs 3D platform, suggesting that microenvironment affects drug response. Figure 4 represents bright field light microscope images of a representative responder treated with either vehicle ( Figure 4A) or 10 μmol/L Ara-C ( Figure 4B) and IC50 curves of the same responder ( Figure 4C) and a non-responder ( Figure 4D).

| Determination of gene signatures that predict AML patient responsiveness to Ara-C
The 23 tested AML donors were classified as 8 responders (IC50 < 100 nmol/L), 5 non-responders (IC50 ≥ 1.7 μmol/L) and 10 moderate responders (100 nmol/L < IC50 < 1.7 μmol/L) as indicated in Figure 4E. To investigate whether certain gene signatures can be used to predict the responsiveness of AML donors to Ara-C, a RNAseq experiment was designed as shown in Figure 4F. Donors Upon Ara-C treatment at IC90 values, a cluster of 248 genes was differentially expressed compared to vehicle-treated samples, when responder and non-responder samples were combined ( Figure 6A).
Interestingly, the transcriptional regulation by Ara-C treatment was more robust in the non-responder samples (365 genes, Figure 6B) than in the responder samples (six genes, Figure 6C). Genes involved in several pathways were regulated by Ara-C either exclusively or more profoundly in non-responder samples, such as the p53 pathway (for example TP53I3, GADD45A, CDKN1A and MDM2).

| Determination of gene fusions and mutations in AML samples
Abnormal gene fusions are frequently observed in AML patients due to chromosomal abnormalities. 3 In our cohort of samples, a total of 14 gene fusions were identified in 13 AML donors by RNAseq after merging the results from FusionCatcher and EricScript Methods. 14,15 Table 3 shows the 14 gene fusions and the samples that they occur.
Eight out of the 14 gene fusions have been previously reported, 4,28-30 leaving 6 novel ones. To confirm these novel gene fusions, the gene fusion supporting reads were aligned to the fusion sequences to check the break points (supplemental Figure 2).
For analysis of somatic mutations by WES, the tumour only pipeline was selected because of the presence of more than 80% of malignant blasts in the majority of BMMCs. Fisher exact test was applied to identify mutated genes to differentiate bone marrow F I G U R E 6 Gene expression analysis following Ara-C treatment. A, The heat map contains the 248 genes that are significantly differentially expressed (±1.5-fold change with FDR_BH < 0.01) following Ara-C treatment compared to vehicle treatment in all 8 donor samples. The four responders and four non-responders were combined as one group. Depicted are the fold change values of each individual Ara-C-treated sample vs the pool of vehicle-treated samples (eight donors). B, The heat map contains the 365 genes that are significantly differentially expressed (± 1.5-fold change with FDR_BH < 0.01) following Ara-C treatment compared to vehicle treatment in the four non-responders. Depicted are the fold change values of each individual Ara-C-treated sample vs the pool of vehicle-treated samples (eight donors). C, The corresponding heat map of the Ara-C signature in responder donor samples (six genes) sample growth status. Different variants within a gene were merged.
If a gene has at least one somatic variant, this gene was assigned as a 'mutant'; otherwise, it was a 'non-mutant'. Only genes mutated in at least two samples were included. Due to the small sample size, these associations are not statistically significant after correcting the P-values with false discovery rate. However, interesting trends were seen from the 12 genes with P-value < .05 (Table 4A). If only the non-synonymous somatic mutations were involved into the association analysis, only one gene differentiates bone marrow sample growth status with P-value < .05 ( the Ara-C H50 values and could be used to predict Ara-C responsiveness (Table 5A). If only the non-synonymous somatic mutations were counted, five genes (HDAC8, CRYBG3, PRSS3, MYH7 and ZAN) were found to associate with Ara-C responsiveness with Pvalue < .05 (Table 5B). in our top 100 gene list for predicting Ara-C responsiveness and was reported to be one of the potential candidate pathway genes of relevance for pharmacogenetic studies on ara-C and other nucleoside analogs. 33 A most recent gene profiling work performed with uncultured and untreated samples presented a large functional genomic data set of primary AML bone marrow mononuclear cells and revealed new markers and mechanisms of drug sensitivity and resistance. 6 Additionally, the ex vivo drug testing of 122 small molecule inhibitors was done in 2D culture for a short duration of 4 days and did not include Ara-C. There is about 20% overlap in mutated genes identified from this study compared to ours, indicating genetic markers predicting common drug responsiveness.

AUTH O R CO NTR I B UTI O N S
All authors reviewed the manuscript and approved its content.

DATA AVA I L A B I L I T Y S TAT E M E N T
This manuscript does not contain sharable data.