Cell-of-origin Discovery in Infant Leukemia through Integration of 3D Models and Patient Transcriptomic Data

Summary

This protocol aims to explore the cellular composition and temporal placement of candidate cell-of-origin for leukemias that arise in utero by integrating single-cell and/or bulk RNA sequencing from hemogenic gastruloids with patient data.

Abstract

Pediatric hematological malignancies remain challenging to investigate and model due to the age group-specificity of certain genetic abnormalities. In utero origin has been demonstrated for a subset of pediatric leukemias, placing their respective cell of origin (CoO) during embryonic development. We recently reported a 3D hemogenic gastruloid (haemGx) model of embryonic blood formation derived from mouse embryonic stem cells, resolving the spatio-temporal complexity of developmental hematopoiesis. Importantly, it allows genetic engineering to introduce disease-relevant mutations. Using haemGx, we modeled the most common acute myeloid leukemia exclusive to infants (infAML), subtype t(7;12)(q36;p13), which arises in utero and is characterized by MNX1 overexpression. Here, we detail a method to define susceptibility to specific mutations that integrate phenotypic and transcriptional changes in the haemGx system and compares them with patient data. By proxy of our MNX1-overexpression haemGx, we show a pipeline from cell engineering to downstream analyses of leukemogenic potential. In particular, we focus on the clinical relevance of the model by integrating single-cell and/or bulk RNA sequencing from the haemGx platform with patient data to extract cellular composition and temporal placement of the putative CoO. This method is adaptable to the introduction of other oncogenic mutations, chromosomal rearrangements, or epigenetic modifications, as well as to chemical perturbations, including drug vulnerability and growth factor dependence. This flexibility allows for broad application across diverse disease contexts, enabling mechanistic dissection of how specific alterations disrupt early developmental trajectories with clinical relevance.

Introduction

Pediatric leukemias can exhibit age-specific genetic abnormalities that distinguish them from those in older patients. Age-specific features configure distinct biological properties of the lineages from which the malignancies arise—their cell of origin (CoO)¹. In particular, identifying a CoO for infant leukemias (infAML) remains challenging. Leukemia initiation in utero^2,3,4, is confounded by the spatio-temporal complexity of developmental hematopoiesis, which utilizes yolk sac (YS), aorta-gonad mesonephros (AGM), and fetal liver (FL) niches in a time-dependent manner for unique cell type specification (YS, AGM) and expansion/maturation (FL)⁵.

The identification of CoO has relied on the detection of leukemia-associated abnormalities (LAA), for example, cytogenetic aberrations or fusion genes in different hematopoietic compartments by fluorescence in situ hybridization (FISH) or polymerase-chain reaction (PCR)-based methods⁶. Functionally, the introduction of LAA in mice via transplantation of transduced hematopoietic cells or via germline genetic manipulation can confirm their leukemogenic potential by expansion of specific populations, albeit not always with complete recapitulation of clinical features⁷. The increasing availability of next-generation sequencing data has improved the characterization of transcriptional profiles and cellular compositions in both experimental models and patient samples^8,9,10, allowing tracing potential CoO and understanding their trajectories. Advanced tools such as patient-derived organoids and induced pluripotent stem cell (iPSC) technology have allowed the investigation of LAA in physiologically relevant conditions to their origin, such as appropriate cellular backgrounds and/or supporting microenvironment^11,12.

In infant forms, the identification of CoO is constrained by the availability of models that recapitulate the fetal environment in space and time, where CoO is likely to be found. Several pediatric abnormalities have been mapped to embryonic windows^13,14, with direct evidence for t(8;21)/RUNX1-RUNX1T1 and t(7;12)/MNX1-ETV6 arising in utero^15,16. The myeloproliferative disorder juvenile myelomonocytic leukemia (JMML) has been shown to arise from YS-specified erythro-myeloid progenitors (EMP) prior to FL colonization¹⁷. Similarly, the CoO for infant ALL harboring t(4;11)/KMT2A-AF4 was pinpointed at the FL lympho-myeloid primed progenitor (LMPP)^18,19. Nevertheless, CoO discovery experiments are often performed by transplantation or ex vivo cultures, limiting the ability to simultaneously capture dynamic changes and supporting structures.

We recently used hemogenic gastruloids (haemGx) to model the rare form of infAML carrying t(7;12)(q36;p13)²⁰, which results in ectopic MNX1 overexpression²¹. HaemGx is a scalable 3D model of developmental hematopoiesis derived from mouse embryonic stem cells (mESC), which achieves stepwise recapitulation of mesoderm formation, hemogenic endothelium (HE) specification, endothelial-to-hematopoietic transition (EHT), and hematopoietic progenitor emergence, in time-congruent YS-like and AGM-like niches²⁰. The unique association of t(7;12) in infancy^21,22, and the recent discovery of its antenatal origin¹⁶ are indicative of a development-stage-specific cell type underlying the leukemic effects of the translocation. In fact, MNX1 overexpression can transform FL but not adult hematopoietic cells^23,24. Using haemGx, we placed its putative CoO at the HE-to-EMP transition, closely resembling transcriptional profiles observed in t(7;12) patient samples²⁰.

Here, we describe an integrative approach to explore transcriptional profiles of LAA within an embryonic context using haemGx, with the overall goal of CoO discovery (Figure 1), based on our previous work on modeling t(7;12) AML in haemGx²⁰. Protocol section 1 describes the use of engineered LAA in haemGx and downstream analyses to assess leukemogenic features, while protocol section 2 details an in silico method to infer the temporal placement of candidate CoO, comparing RNA sequencing from haemGx and patient data. This method is most suitable for the discovery of cell types involved in hematological malignancies with an embryonic component and has been optimized for use with mESC to ensure full compatibility with transcriptomic data.

Figure 1: Overview of methodologies to be used for cell-of-origin discovery in leukemia by integrating haemGx phenotypic data with patient transcriptomics​. This approach allows the engineering of leukemia-associated abnormalities in mESC to be investigated in a hemogenic gastruloid (haemGx) model via downstream molecular and bioinformatics analyses. Abbreviations: LAA = leukemia-associated abnormalities; mESC = mouse embryonic stem cells; haemGx = hemogenic gastruloid model; GSEA = gene set enrichment analysis. Please click here to view a larger version of this figure.

Protocol

1. Use of engineered leukemia-associated abnormalities (LAA) in haemGx and downstream analyses to assess leukemogenic features

1. Culturing mESC with LAA for haemGx aggregation
   NOTE: mESC can be edited to introduce the LAA of interest by overexpression systems (here, using lentiviral transduction, MNX1 is overexpressed using a pWPT-LSSmOrange-PQR under the EF-1α promoter²⁵), or gene editing approaches (e.g., CRISPR/Cas9 to introduce mutations). A control line (e.g., transduction by empty vector or Cas9-only control) is also generated. Various reporter lines can be used, for instance, Kdr(Flk1)-GFP²⁶, Sox17-GFP²⁷, or T/Bra-GFP²⁸. It is strongly recommended to validate the LAA before proceeding; this can be achieved by qPCR and/or flow cytometry to assess overexpression (see example in Supplemental File 1 Supplemental Figure S1AB) or DNA sequencing for other genomic edits.
   1. Coat tissue culture vessels with 0.1% gelatin in PBS for 10 min at room temperature (5 mL for a 25 cm² flask).
   2. Prior to the addition of cells, discard all gelatin and add ES-LIF medium (500 mL of Glasgow MEM BHK-21, 50 mL of fetal bovine serum, 5 mL of glutamine supplement, 5 mL of MEM Non-Essential Amino Acids, 5 mL of sodium pyruvate solution (100 mM), and 1 mL of 2-Mercaptoethanol 50 mM). Murine Leukemia Inhibitory Factor (LIF) is added at 1,000 U/mL.
   3. Culture mESC at a density of 300,000 cells per T25 flask. Medium must be changed daily. When passaging, reseed the cells at the same density. Allow the cells to expand for at least two passages before using them for aggregation of haemGx. Once the culture reaches 60–80% confluency, prepare for aggregation.
      ​NOTE: Appropriate maintenance and health of the mESC culture is crucial for haemGx assembly. Inspect cell culture under an inverted microscope daily to assess morphology and confluency. The culture should appear with compact, well-defined islands with smooth colony edges (Supplemental Figure S1C).
2. Generation of haemGx from mESCs
   1. Culture haemGx in N2B27 medium (refer to Table 1 for recipe).
   2. Prewarm PBS containing Ca²⁺ and Mg²⁺, ES-LIF medium, N2B27, and Trypsin–EDTA in a 37 °C water bath before use.
   3. Remove the culture medium from the tissue culture flask and gently rinse the cells 2x with 5 mL of prewarmed PBS. Aspirate the PBS, then add 1 mL of prewarmed 0.25% Trypsin–EDTA to detach the cells. Incubate the flask at 37 °C for approximately 5 min or until the cells have completely detached from the surface. Gently pipette the suspension up and down with a 5 mL pipette to break up any remaining clumps.
   4. Add 5 mL of ES-LIF medium to neutralize the Trypsin–EDTA and gently rinse the culture surface to improve cell recovery. Transfer the cell suspension to a 15 mL centrifuge tube and spin at 400 × g for 5 min to pellet the cells.
   5. Carefully remove the supernatant and resuspend the pellet in 1 mL of ES-LIF medium. Mix thoroughly by pipetting to obtain a single cell suspension, then count the cells using Trypan blue and a hemocytometer.
      NOTE: Cell viability must be at least 80%.
   6. Prepare a working suspension at 10 cells/µL, which corresponds to 400 cells per well in a 40 µL droplet. Calculate the total number of cells required using this formula:
      Total cells = number of wells × cells per well × 1.1
      Where 1.1 accounts for a 10% excess to compensate for pipetting variability.
      NOTE: This ensures sufficient volume and cell number for all wells while minimizing variation between aggregates. For 60 wells with 10% excess, this equals 2.64 × 10⁴ cells in 2.64 mL of N2B27 medium.
   7. Transfer 2.64 × 10⁴ cells into 5 mL of prewarmed PBS in a clean 15 mL centrifuge tube and centrifuge at 400 × g for 5 min. Carefully aspirate the PBS without disturbing the cell pellet, then gently add another 5 mL of prewarmed PBS for a second wash. Disperse the pellet between the washes and centrifuge again at 400 × g for 5 min.
   8. Gently remove as much PBS as possible without disrupting the pellet. Resuspend the cells in 1 mL of warm N2B27 medium using a P1000 pipette to obtain a single cell suspension. Finally, add additional N2B27 to reach the required final volume (e.g., add 1.64 mL).
   9. Dispense 40 µL of the prepared cell suspension into each of the 60 inner wells of an ultra-low adherence U-bottom 96-well plate to initiate aggregation. To minimize evaporation and maintain consistent humidity across the plate, add approximately 150 µL of PBS to each of the outer wells. Cover the plate with the lid and verify the presence and even distribution of cells under a microscope. Following confirmation, incubate the plate for 48 h in a humidified incubator at 37 °C with 5% CO₂ to allow the cells to aggregate.
      NOTE: The plate can be centrifuged briefly at this stage using a plate centrifuge to aid aggregation (400 × g for 3 min).
3. Medium change and cytokine supplementation
   NOTE: The medium is to be changed daily at the same time. Figure 2 depicts a schematic representation of the protocol. A detailed method for medium change for gastruloids in 96-well plates has been previously described²⁹.
   1. After 48 h, gently add 150 µL of fresh prewarmed N2B27 medium containing 100 ng/mL Activin A and 3 µM GSK3 inhibitor CHIR99021 to each well (this combination of Activin A + CHIR99021 is hereafter known as the AC pulse). Ensure the pipette tip contacts only the well wall, without disturbing the gastruloids.
   2. At 72 h, remove 150 µL of medium from each well to terminate the AC pulse. To avoid disturbing gastruloids, hold the pipette at approximately 30° near the bottom of the well and aspirate gently at ~50 µL/s. 
   3. After removing the AC pulse, gently add 150 µL of fresh prewarmed N2B27 medium supplemented with VEGF and FGF2, each at a final concentration of 5 ng/mL, and incubate for an additional 24 h.
   4. At 96 h, remove 150 µL of medium from each well and replace it with 100 µL of fresh prewarmed N2B27 containing VEGF and FGF2 at the same concentrations as above.
   5. Continue the medium exchange process by removing 100 µL and adding 100 µL of N2B27 with VEGF and FGF2 per well every 24 h until 144 h.
   6. At 144 h, remove 100 µL of medium from each well and replace it with 100 µL of fresh N2B27 medium containing VEGF (5 ng/mL), FGF2 (5 ng/mL), and mSHH at a final concentration of 20 ng/mL.
   7. At 168 h, remove the mSHH-containing medium (100 µL per well) and replace it with 100 µL per well of prewarmed N2B27 supplemented with VEGF (5 ng/mL) and a cytokine cocktail consisting of mSCF (100 ng/mL), mTPO (20 ng/mL), and mFLT3L (100 ng/mL). Incubate for an additional 24 h, then repeat the same medium replacement at 192 h using the same cytokine cocktail.
   8. At 216 h, the protocol is complete; collect gastruloids for downstream analyses.
4. Collection and dissociation of haemGx
   NOTE: To enable detailed cellular and molecular analyses, gastruloids can be dissociated from their 3D structures into single-cell suspensions for downstream characterization at any time point. It is strongly recommended to confirm the persistence of the desired LAA until endpoint (e.g., by qPCR as shown in Supplemental Figure S1D).
   1. Gently remove 100 µL of medium from each well using a single-channel pipette. Be sure not to disturb the gastruloids at the bottom of the wells.
   2. Using a P1000 pipette, transfer the remaining 100 µL of medium along with the gastruloid into a labelled 1.5 mL microcentrifuge tube. Allow the gastruloids to settle naturally at the bottom of the tube (1–2 min). Once settled, carefully remove the supernatant.
      NOTE: Multiple gastruloids can be collected per tube.
   3. Wash the sample by adding 1 mL of PBS (without Ca²⁺ or Mg²⁺) to remove residual medium.
   4. Centrifuge the tubes at 400 × g for 5 min at room temperature. Gently discard the PBS supernatant, leaving the gastruloid pellet untouched.
   5. Add 200 µL of prewarmed recombinant enzyme-based cell dissociation reagent. Ensure gastruloids are fully immersed in the solution. Place the tubes in a 37 °C water or pebble bath for 5 min to facilitate dissociation.
   6. After incubation, use a P200 pipette to vigorously pipette the sample up and down at least 10x to mechanically aid dissociation. Observe a small aliquot under the microscope to check for single-cell suspension. If large clumps remain, repeat the incubation step once more for 2–3 min and pipette again.
   7. Add 300 µL of ES-LIF medium to each tube to neutralize the cell dissociation reagent. Pipette up and down several times to ensure all cells are detached and there are no clumps.
   8. Centrifuge again at 400 × g for 5 min. Carefully discard the supernatant to remove residual enzyme and debris.
   9. Gently resuspend the resulting cell pellet in 200–1,000 µL of ES-LIF medium, depending on the pellet size and intended downstream applications. 
5. Downstream analysis: Colony forming and replating assay
   NOTE: Optimal seeding concentrations need to be determined empirically, as haemGx cells may have different colony-forming capacity at the first plating. However, untransformed cells cannot replate indefinitely, and differences between conditions can be discerned at later replating (P3–5).
   1. Plating on methylcellulose
      1. Thaw a prealiquoted 3 mL tube of methylcellulose (Mouse Methylcellulose Complete Media) on ice. Once fully thawed, vortex thoroughly to ensure the medium is well mixed, given the high viscosity of methylcellulose. Use the top speed setting for 1–3 periods of 5 s; adjust timings and repeats empirically. As vortexing introduces air bubbles, allow the methylcellulose to rest at room temperature until the bubbles dissipate completely.
         NOTE: Ensure that the entire volume of the aliquot is being mixed; angle the tube while vortexing to check that the methylcellulose medium at the bottom/tip of the tube is also being mixed. Thorough mixing is critical to achieve uniform growth factor concentrations.
      2. Meanwhile, prepare the gastruloid single-cell suspension for the colony-forming assay. Ensure that a total of 100,000 cells (or empirically determined optimal concentration) is resuspended in 300 µL of N2B27 medium.
      3. Once the methylcellulose is free of bubbles, carefully add the 300 µL cell suspension to the 3 mL of methylcellulose using a P1000 pipette. Dispense slowly into the methylcellulose medium immediately below the meniscus to minimize additional bubble formation; do not pipette to mix. Vortex the mixture thoroughly to achieve a uniform distribution of cells within the methylcellulose. 
      4. Allow the methylcellulose medium to rest on ice for at least 5 min to eliminate any remaining bubbles. After debubbling, plate the medium into two wells of a 6-well plate by dispensing half the volume, typically 1.3–1.5 mL per well (corresponding to 50,000 cells per well), using a 5 mL serological pipette and pipettor.
         NOTE: Make sure that the methylcellulose is aspirated and dispensed slowly using the low-speed setting on the pipettor and a continuous flow.
      5. Add sterile PBS to the surrounding empty wells to maintain humidity and prevent drying. Incubate the plate at 37 °C with 5% CO₂ for 7–10 days before colony scoring.
      6. After 7–10 days of incubation, image the plates and score the colonies. To calculate colony-formation capacity (frequency), calculate the number of colonies/number of cells seeded.
         NOTE: haemGx produce heterogeneous colonies encompassing hematopoietic and non-hematopoietic cells. Standard scoring of hematopoietic colonies can be used at late time points (192–216 h); however, earlier progenitor-like cells may produce undefined morphologies³⁰.
      7. Score total colony numbers in preliminary analyses and focus on the morphology of late replating colonies at later stages. Adjust colony scoring criteria empirically, but define a colony as having at least 50 cells.
         NOTE: It is important for individual users to keep qualitative and quantitative scoring criteria constant. At different steps of the protocol, colonies are scored independently and blindly by a second observer to ensure consistency and reproducibility of results. Qualitative colony scoring is checked against cell morphology in cytospins (see below).
   2. Colony replating
      1. Prewarm PBS to room temperature and add 1.5 mL to each well. Gently resuspend the methylcellulose by slowly pipetting up and down using a P1000 pipette set to 500 µL. Avoid foaming up the medium.
      2. Once the methylcellulose has been fully dispersed with no visible clumps remaining, transfer the entire contents of each well into a 15 mL conical tube. Rinse the same well with an additional 1.5 mL of PBS to collect any remaining cells and pipette slowly to resuspend the remaining methylcellulose; add this suspension to the same tube. Perform a final rinse with 500 µL of PBS to recover residual cells and pool with the previous washes.
      3. Centrifuge the collected suspension at 400 × g for 5 min at room temperature to pellet the cells. Carefully aspirate the supernatant, then wash the pellet with 5 mL of PBS to remove any remaining methylcellulose.
      4. Centrifuge again under the same conditions (400 × g for 5 min). Discard the supernatant and resuspend the final cell pellet in 200–1000 µL of N2B27 culture medium (volume adjusted based on pellet size).
      5. Count viable cells using a hemocytometer.
      6. Replate the cells by following the same colony-forming assay preparation protocol described in Steps 1.5.1.1–1.5.1.5. Continue serial replating as required to assess colony-forming potential over successive generations.
      7. Dissociate colonies into cell suspensions at the desired endpoint as described in step 1.5.1 and subject them to downstream characterization, including cytospin preparations and Giemsa Wright staining for detailed cellular morphology (described in detail elsewhere³¹).
         NOTE: In summary, section 1 describes the generation of haemGx that harbor LAA. Once the resulting phenotype is assessed in vitro, the transcriptomes of these cells are analyzed using a computational approach. Gene expression signatures derived from the haemGx model can be compared with bulk transcriptomic datasets from leukemia patients to understand phenotypic similarities. Using Gene Set Enrichment Analysis (GSEA), these signatures are mapped with known populations across differentiation time points to infer a temporal window of susceptibility to the LAA (section 2).

Figure 2: Timeline of haemGx protocol. Schematic representation of the generation of haemGx production from mESC cells over a 216 h protocol, highlighting the addition of an A/C pulse at 48 h, followed by the specific chemical cues of cytokines from 144 h to 216 h to promote the specification for hemato-endothelium. Abbreviations: mESC = mouse embryonic stem cells; haemGx = hemogenic gastruloid model. Please click here to view a larger version of this figure.

2. In silico method to infer the temporal placement of candidate CoO, comparing RNA-seq from haemGx and patient data

1. Cell type mapping of bulk patient transcriptomic data using GSEA
   NOTE: This analysis allows the mapping of cell types that are overrepresented in bulk RNA-seq data of specific patient subtypes compared to other classes using Gene Set Enrichment Analysis (GSEA)³². GSEA requires specific normalized units for RNA-sequencing data: https://docs.gsea-msigdb.org/#GSEA/GSEA_and_RNA-Seq/.
   1. Download patient data from the appropriate repository, for example, TARGET data for pediatric acute myeloid leukemia (AML).
      NOTE: TARGET data are accessible from multiple sources: UCSC Xena Browser: https://xenabrowser.net/datapages/³³, cBioPortal https://www.cbioportal.org/datasets³⁴, GDC Data Portal https://portal.gdc.cancer.gov/³⁵
   2. Download clinical data and RNA-sequencing data in TPM or normalized counts to be compatible with GSEA (Supplemental Figure S2). For TARGET-AML, download TPM data and clinical data to identify the desired cohort. Select the patients to be included in the analysis making use of phenotype and clinical annotations (e.g., see example for cBioPortal, Supplemental Figure S3).
   3. Download cell type gene sets. Gene sets in a compatible format for GSEA are available from EnrichR³⁶: https://maayanlab.cloud/Enrichr/#libraries. Download PanglaoAugmented 2021³⁷ for cell type gene sets, also available from PanglaoDB: https://panglaodb.se/index.htmL.
   4. Download and install GSEA software: https://www.gsea-msigdb.org/gsea/downloads.jsp.
   5. Prepare files for GSEA software using a custom gene set.
      NOTE: This analysis can also be run using GSEA’s inbuilt signature collection MSigDB (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp) for further biological insights.
   6. Prepare the files required for GSEA according to the template (Supplemental Figure S4), populated in spreadsheets and saved as .txt files. Change the file extension to .gct, .cls, or .gmt accordingly:
      .gct file: expression file with samples, genes, and expression values
      .cls file: phenotype labels matching the .gct file; contains the classes/phenotypes to be compared.
      .gmt file: gene set file containing gene set names and corresponding genes defining each set.
      NOTE: Different gene nomenclatures can be used. Gene symbols often contain duplicated values that can affect results. Ensembl identifiers are preferred.
   7. Upload all files into GSEA. Select Load Data | Method 1 | Browse for Files…
   8. Run GSEA to compare two phenotype classes.
   9. Select Run GSEA. Populate the required fields.
      • Expression dataset: select the file containing counts for desired samples.
      • Gene set database: select the gene sets to compute (e.g., PanglaoAugmented 2021).
      • Number of permutations: 1,000 is the default (range: 1,000–10,000 for higher statistical resolution).
      • Phenotype labels: select the .cls file containing phenotype labels for the categories to compare.
      • Collapse/Remap to gene symbols: Depending on the format of the .gct file, select whether to remap gene identifiers to gene symbols to match gene set file (.gmt). Choose from three options:
        • Collapse: if multiple gene identifiers map to the same gene symbol, GSEA will collapse them into a single gene by its maximum expression value. To be used for microarray datasets with probe IDs.
        • Remap-only: converts gene identifiers without collapsing duplicates (i.e., gene identifiers in .gct do not match identifiers in .gmt).
        • No collapse: does not convert gene identifiers (i.e., gene identifiers match identifiers in .gmt).
      • Permutation type: phenotype (shuffles phenotype labels to compute gene statistics) or gene set (randomly reassigns genes in each set to compute probability of gene membership). Gene set is recommended for RNA-seq data, if fewer replicates are available (≤4) and to minimize phenotype variability.
      • Chip platform: select the ID type of genes in the .gct file.
   10. Populate the basic fields.
       Metric for ranking genes: for categorical classification (not continuous), choose Signal2Noise, tTest, (log2)Ratio_of_classes, or Diff_of_classes depending on data types. For RNA-seq data, choose Ratio_of_classes (absolute fold change between classes) or log2_Ratio_of_classes (relative fold change between classes).
       NOTE: Parameters are adjustable depending on specific needs and data type. Refer to GSEA website for details.
   11. Select Run. The status of the analysis is reported in the GSEA reports box. Once complete, click on COMPLETE to directly open the summary of the analysis. The full analysis is saved locally in the specified folder.
   12. To retrieve results and analysis, open the index file (index.html) to view a summary of results. Individual enrichments of gene sets with enrichment curves are also generated. To extract the leading edge (or core enriched genes that account for the enriched set), select all genes marked by “Yes” in the “core enriched” column for each gene set (access via index.html | Detailed enrichment results in html format | GS Details | Table: GSEA details)
       NOTE: Additional downstream analysis of enriched genes in the leading edge can be conducted by gene ontology / term enrichment analysis (e.g., via EnrichR).
   13. Visualization
       Find gsea_report_*.tsv files for each phenotype in the folder to export results for graphical analysis.
2. Inference of differential cell type representation by timepoint in haemGx differentiation
   NOTE: This analysis enables mapping bulk RNA-seq data from disease-modeled haemGx or patient samples onto a temporally resolved atlas of haemGx differentiation, with the aim of pinpointing differentially enriched cellular populations at specific developmental time points.
   1. Retrieve bulk RNA-seq data from public databases (see 2.1.1 for TARGET AML) or RNA-seq sequenced at desired timepoints in (e.g. haemGx-MNX1 vs haemGx-EV at endpoint 216 h). Units must be compatible with GSEA (see 2.1.1).
   2. Retrieve cluster classifiers for normal haemGx differentiation protocol²⁰ at specific timepoints from https://github.com/deniseragusa/haemGx-CoO/releases/tag/v1to be used as custom gene sets (.gmt).
      NOTE: The haemGx classifier list was derived from single-cell RNA-seq of normal haemGx differentiation protocol and represent time-point resolved populations in haemGx. Ortholog mapping was applied using BioMart to be compatible in human / mouse comparisons. However, exact species matching is still limited by inherent differences in developmental processes between mouse and human.
   3. Run GSEA and retrieve results as step 2.1.1.
   4. Visualization. A time-resolved UMAP image is available at https://github.com/deniseragusa/haemGx-CoO/releases/tag/v1 for visualization of cluster mapping.

Results

We used haemGx to model the most common infAML subtype, t(7;12)(q36;p13), via MNX1 overexpression as a proxy, and to infer the developmental window of susceptibility and its clinical relevance to patient transcriptomics using GSEA.

To introduce LAA, we used lentiviral transduction to introduce MNX1 overexpression with the pWPT-LSSmOrange-MNX1-OE-PQR vector to overexpress MNX1 (mESC-MNX1) (Supplemental File 1 Supplemental Figure S1AB) and used an empty vector as a control (mESC-EV). MNX1-overexpressing mESC were cultured to form haemGx (haemGx-MNX1), alongside the empty vector control (haemGx-EV). Appropriate haemGx assembly can be assessed by inverted microscopy for successful aggregation—cells form a circular, compact structure with defined borders (Figure 3A). If using a Flk1-GFP marker line, the appearance of a polarized GFP signal at 96 h confirms the induction of patterned endothelium (Figure 3A).

Differences in morphology or reporter expression can be tracked by imaging and quantified for differential phenotypes, for example, by size difference between control and MNX1-haemGx (Figure 3B). Downstream analyses following haemGx disassembly can include flow cytometry staining for specific markers to identify timepoints with differential effects (e.g., an expansion of Kit+ cells at 144h by MNX1, Figure 3C). Leukemogenic potential can be assessed by colony formation and replating capacity in methylcellulose, quantifying the number of colonies generated at each replating under different conditions, with the expectation that non-transformed cells will extinguish colony formation at early platings (Figure 3D). Figure 3E shows representative colony outputs from haemGx at different replatings, highlighting colony diversity and quantity under different conditions; in this case, haemGx-EV fails to generate colonies after replating 2 P2 (Figure 3E). Colony scoring criteria can be adapted for specific experimental designs, including by total colony number or morphology. Here, potential colonies from haemGx-MNX1 are shown, with scorable colonies pointed by the full arrows (Figure 3F). Smaller aggregations of cells or sparse individual cells are not scored as colonies (dashed arrows in Figure 3F). Colony quantifications from haemGx disassembled at different timepoints can be plotted in Figure 3G to present colony frequency and serial replating capacity as in vitro measures of transformation. Additional assays for characterizing colony morphology include microscopy or Giemsa-Wright staining of cytospin preparations. Disassembled haemGxs can then be subjected to virtually any omics-based analyses following standard protocols.

Figure 3: Representative results using haemGx to model t(7;12) infant AML. (A) Images of individual haemGx of MNX1 vs EV conditions acquired using a microscopy plate reader. Scale bar = 300 µm for all panels. The green signal indicates polarization of the FLK1-GFP marker throughout the 216 h protocol. (B) Time-course size differences of haemGx-MNX1 vs haemGx-EV calculated on cross diameter using images acquired in A. Mean ± SD of three replicate experiments; two-tailed t-test, p<0.05 (*), 0.001 (**), 0.0001 (***), and 0.00001 (****). (C) Representative example of flow cytometry analysis of C-Kit⁺ marker between haemGx-MNX1 and haem-Gx EV disassembled at each time point between 120 h and 216 h. Mean ± standard deviation of 3–7 independent experiments; two-way ANOVA and Sidak’s multiple comparison test significant at p<0.05 for C-Kit⁺ cells only (construct contribution to variance p=0.0191; 144 h comparison *p=0.0190). (D) Schematic representation of serial replating of colony-forming assays to assess leukemogenic potential. Serially replating colonies (bottom) indicates a transformed phenotype. P=plating. (E) Representative photographs of serial replating of colony-forming cell assays initiated from haemGx-MNX1 and haemGxEV at different platings. Plate stitching was performed at 4X magnification on a plate reader. (F) Representative image of potential morphologically distinct colony outputs from haemGx-MNX1. An example of scoring criteria is shown by arrows: full arrows indicate a scorable colony, while dashed arrows show small cell depositions not to be counted as colonies. (G) Quantification of colony-replating efficiency of EV and MNX1 144 h and 216 h-haemGx cells. Mean + SD of n>3 replicates; Kruskal-Wallis with Dunn’s multiple comparison testing at significant q<0.05. Images in panels A, B, C, and G are taken from Ragusa et al.²⁰. Please click here to view a larger version of this figure.

GSEA can be used to deconvolute bulk RNA-seq data to infer cell type enrichment and clinical similarity to patients when single-cell approaches are not available. Here, patient data was used from the TARGET database to map cell type composition of t(7;12) patients against other AML forms (we combined inv(16), KMT2A (or MLL/KMT2A), normal karyotype, and t(8;21)), by revealing over-represented transcriptional programs of cell identity by proxy of NES values by GSEA (Figure 4). The output of GSEA analysis is reported in tabular format in Figure 4A and as representative enrichment plots for each gene set in Figure 4B. Correct file formatting and phenotype labelling is critical to complete GSEA successfully (see protocol step 2.1) (Supplemental Figure S2, S3, and S4). NES values reflect the enrichment in respective phenotype categories (here we compared t(7;12) vs other AML): a positive NES value indicates an overrepresentation in the first phenotype label (here, t(7;12) AML); a negative NES value indicates an overrepresentation in the second phenotype label (here, other AML). Statistically significant enrichments (p val < 0.05, FDR & FWER < 0.05) can be visualized by plotting NES values and FDR (in –log10) for each phenotype. Cell types can be grouped in families to help visualization (Figure 4C). These results indicate which cell types are enriched and/or over-represented in the t(7;12) cohort compared to other AML subtypes.

Figure 4: GSEA analysis of cell type enrichment in patient transcriptomes. (A) Example of tabular output from GSEA analysis showing enrichment scores in the “t(7;12)” phenotype. The table reports gene set information (“GS”, “GS DETAILS”, “SIZE”), enrichment score metrics (“ES” and “NES”), statistical significance (“NOM p-val", “FDR q-val", “FWER p-val"), and core enrichment statistics (“RANK AT MAX”, “LEADING EDGE”). Statistically significant gene sets can filtered by FDR or FWER scores < 0.05. (B) Enrichment plots generated by GSEA software for each gene set computed against the chosen phenotype classes. (C) Representative visualization of enrichment scores for cell type analysis of t(7;12) vs other AML. NES and FDR values were extracted from the tabular format. Cell types were grouped by manual categorization. Darker color indicates higher NES values and size represents statistical significance in –log10(FDR). Please click here to view a larger version of this figure.

Timepoint-resolved single-cell RNA-seq from the normal haemGx differentiation can be used as a reference atlas of developmentally accurate cellular populations, which we had previously characterized with reference to mouse datasets of hematopoietic specification¹⁸. In this analysis, we make use of cluster-defining genes of hemato-endothelial and hematopoietic populations, as well as supporting stromal and epithelial clusters (see “identity” classification in Figure 5A), corresponding to specific timepoints of the haemGx protocol (see “time point” in Figure 5A). Using the UMAP plot of this dataset as a map makes cluster identification and timepoint contributions easier. This framework allows the placement of specific transcriptomics onto temporal windows of the haemGx protocol. Using the cluster-defining genes as custom gene sets, GSEA can be run to compare bulk RNA-seq from haemGx with MNX1 overexpression with EV at the endpoint (216 h), with the aim of mapping transcriptional programs that are overrepresented (by cell proliferation or by a block in differentiation) in MNX1-overexpressing haemGx. GSEA results (see Figure 4A,B for output format) can be visualized by plotting enrichments by NES against clusters/timepoints (Figure 5A), indicating temporal windows of susceptibility. Interrogating patient RNA-seq with the same approach, for example, t(7;12) vs other AML subtypes, allows a visual alignment of populations captured by haemGx compared to patient profiles (Figure 5A). As inferred from the cell type mapping in patients supporting enhanced hemato-endothelial specification in t(7;12) patients (Figure 4C), Figure 5A shows that haemGx-MNX1 are enriched for hemato-endothelial clusters 5 and 0 (144 and 192h, respectively), as well as of EMP-like cluster 4 (144h) and MLP-like cluster 8 (192/216h). The 144h/192h timepoint and hemato-endothelial populations match the enrichment seen in t(7;12) patients vs other AML (Figure 5A). We also compared RNA-seq data from t(9;11)/KMT2A-MLLT3 patients, showing a significant enrichment of the MLP-like cluster 8 (Figure 5A), consistent with the myeloid progenitor origin of KMT2A-rearranged AML, which differs from t(7;12). These enrichments can be visualized by projections onto the UMAP of normal haemGx differentiation (Figure 5B).

Figure 5: Temporal mapping of bulk RNA-seq data from disease-modeled haemGx or patient data against normal haemGx differentiation. (A) Bubble plot representation of GSEA NES values and statistical significance by –log10(FDR) of enrichments in specific clusters of haemGx differentiation, comparing overrepresented populations in haemGx-MNX1 compared to haemGx-EV (top), t(7;12) AML samples compared to other AML (middle), and MLL AML (0-2 years old) against other AML (bottom). On the x-axis, cluster numbers corresponding to normal haemGx differentiation are aligned to their corresponding timepoint and cluster identity. (B) UMAP of time-resolved global clustering of single-cell RNA-seq of normal haemGx. Boxes represent statistically significant GSEA results from A and are superimposed to the corresponding cluster number, showing the enriched clusters in haemGx-MNX1 (blue boxes), t(7;12) AML (red boxes), and MLL (green boxes). Images in panels A, B, C, and G are taken from Ragusa et al.²⁰. Please click here to view a larger version of this figure.

This analysis can be expanded to map similarities in cellular composition as a means to identifying conditions / timepoints with the highest alignment with patient features. Here, we compared the enrichment in haemGx or replating CFC colonies in MNX1 vs EV at 144 h and 216 h timepoints, showing the most extensive cell type similarity in replating colonies from 144 h (Figure 6A). For further dissection of gene sets contributing to the enrichment in specific phenotypes, the “leading edge” list of genes can be retrieved for each set, for example, in cluster 2 at 144h in haemGx—MNX1 (Figure 6A) and be subjected to additional ontology—based enrichment analyses. Here, we corroborated the endothelial contribution of this cluster to the haemGx—MNX1 phenotype by cell-type enrichment and by exploring pathway enrichments using the KEGG database (Figure 6B). More complex visualization approaches can be employed to map cell type enrichment with individual genes contributing to the separation (Figure 6C).

Figure 6: Representative visualization strategies for cell type analysis. (A) Bubble plot of enrichment scores (NES) for cell type analysis allowing visualization of multiple comparisons, including t(7;12) vs other AML, haemGx and colony-forming assays from MNX1 vs EV at 144 h and 216 h. Darker color indicates higher NES values and size represents statistical significance in –log10(FDR). (B) Example of leading-edge extraction from tabular results of GSEA. From each gene set, the leading edge corresponds to genes marked by “Yes” in the “CORE ENRICHED” column. The leading edge list can be subjected to additional analyses by gene ontology, for example by EnrichR probing the PanglaoDB repository of cell types, or KEGG repository of biological pathways. Results can be plotted as bar charts showing the statistical significance in –log10(FDR). (C) Example of combinatorial plot integrating cell type analysis from Panglao DB database enriched to leading edge genes in haemGX-MNX1 (blue dots) and t(7;12) AML (red dots) by mapping genes and corresponding time points in the normal haemGx protocol. Please click here to view a larger version of this figure.

Name of Material	Volume	Description
N2 Supplement (100x)	500 µL	Medium Supplement
B27 Supplement (50x)	1 mL	Medium Supplement
2—Merceptoethanol (50Mm)	100 µL	Reducing agent
Glutamax	1 mL	Medium Supplement
Neurobasal Medium	48.7 mL	Media
DMEM/F-12, with Glutamax	48.7 mL	Media

Table 1: Materials and recipe of N2B27 prepared in-house (for 100 mL)

Supplemental File 1: Quality control, data acquisition, and preprocessing workflow supporting haemGx analysis. Supplemental figures illustrate (S1) validation of mESC cultures prior to haemGx, including transduction efficiency, MNX1 overexpression, and assessment of morphology and confluency; (S2) retrieval of RNA-sequencing datasets from the UCSC Xena Browser; (S3) access to corresponding clinical data from cBioPortal; and (S4) preparation of input file formats (.gct, .cls, .gmt) required for gene set enrichment analysis.Please click here to download this file.

Discussion

This protocol is amenable to model a variety of LAA that can be investigated in the context of embryonic hematopoietic development, with the advantage of spatio-temporal resolution and compatibility with established downstream molecular, functional, and biochemical analyses. Here, we focused on gastruloids that recapitulate hemato-endothelial specification to YS-like EMP and AGM-like HSPC emergence; however, other gastruloid / developmental organoid models that recapitulate specification of different tissues and organs can, in principle, be adapted for other disease contexts. Unlike alternative models, such as engineered cell lines or murine transplantation systems, haemGx capture the spatio-temporal dynamics of embryonic hematopoiesis within a physiologically relevant 3D environment, allowing observation of when and where leukemogenic transformation may occur.

Successful generation of haemGx depends on appropriate mESC culture, daily medium supplementation, and microscopy-based culture morphology inspection. Troubleshooting of gastruloid aggregation steps has been described before²⁹, with aggregation failure relating to incorrect cell counting or medium supplementation. Critical steps in the haemGx protocol lie in the initial cell culture conditions for mESC cells and in the medium supplementation at the specified timepoints, to ensure appropriate dynamics and comparisons between control and LAA lines. It is also recommended to validate the presence of the desired modification in the mESC line and in haemGx until endpoint, for example, by confirming successful integration of a vector by flow cytometry detection of a fluorescent reporter, or presence of the oncogene or the mutation by qPCR or sequencing. The effects on haemGx differentiation will differ depending on the abnormality of interest, requiring exploratory experiments to identify relevant phenotypes. Recommended markers to test the appropriate dynamics of key populations in the haemGx protocol include Flk1/KDR (replaceable by VE-cadherin) from 96 h onwards, CD41 and c-Kit at 144 h, and CD45 at 216 h; phenotype refinement by combination of CD41 and/or CD45 with Kit, CD31 and CD34 extend characterization of EHT and HSPC-like cells²⁰. Failure to appropriately differentiate and express the relevant markers is due to errors in supplementation dosage and daily consistency; it is recommended to store supplements as aliquots to avoid thaw-freezing cycles.

Common difficulties in the collection of haemGx for downstream analyses include failure to dissociate in the dissociation reagent. The volume of dissociation reagent can be increased at 50 µl increments and frequent pipetting facilitates the dissociation into homogeneous solution. The capacity to perform molecular and biochemical assays is directly proportional to the number of gastruloids generated, requiring up- or downscaling depending on the application. On average, one gastruloid at endpoint (216 h) yields 50,000–70,000 cells. Colony-forming assays, in particular, are sensitive to cell numbers plated, as the resolution of colony formation detection can be affected by overcrowding and the inability to discern individual colonies. Common issues with colony formation include failure to form colonies or overconfluency, requiring empirical determination of suitable cell numbers for plating, especially in LAA conditions of unknown clonogenic potential. A limitation of this approach is that it has only been applied to the identification of clear hematopoietic and hematopoietic-like colonies²⁰. However, other morphologies have not been systematically characterized by colony-formation and rather rely on single-cell transcriptomics of haemGx populations instead.

The advantage of a 3D system encompassing different cell types is that it enables exploration of a wide range of ontogeny, which can also be studied through high-throughput sequencing methods. While single-cell approaches would be preferred to capture the cellular composition at high resolution, GSEA³² provides a simple tool to resolve heterogeneity in bulk transcriptomes and identify overrepresented programs that match cell types. This is particularly useful to maximize the value of clinical data, where single-cell data are not routinely or universally available. The PanglaoDB repository encompasses an extensive atlas of gene sets for cell types from integration of mouse and human single-cell RNA-seq datasets, which also capture embryonic markers, allowing a comprehensive coverage of tissue identities and lineages³⁷. Additional gene sets with more detailed cell types (in-house markers or from other repositories^36,37) can also be incorporated in this protocol for targeted analyses. The species differences in genes between the mouse haemGx system and human patient data may be a limitation, which can be attenuated by filtering gene matches for ortholog identification (for example, https://www.informatics.jax.org/homology.shtmL). Implementation of a human version of the haemGx model will overcome this limitation. Successful implementation of GSEA is dependent on the correct formatting of input files, which may cause the software to flag errors (e.g. “ERRORS #:1 Parsing trouble… " "Bad format expect ncols: [x] but found: [y] on line >”). Other common points of failure include incompatibility of gene nomenclature between GSEA .gct and .gmt files, raising the error “The collapsed dataset was empty when used with chip:…”. A full troubleshooting guide is available for the GSEA software (https://docs.gsea-msigdb.org/#GSEA/GSEA_FAQ/).

In conclusion, this protocol allows exploratory mapping of developmental windows of susceptibility to specific LAA, by integration of phenotypic observations in the haemGx model of embryonic hematopoiesis and transcriptomics data.

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
Activin A Plus	Qkine	Cat. #QK005	Peptide, recombinant protein
B-27 Supplement (50x), serum free	Thermo Fisher Scientific	Cat. #17504044	Medium supplement
CHIR99021 (Chiron)	BioGems	Cat. #2520691	Peptide, recombinant protein
Gibco 2-Mercaptoethanol (50 mM)	Fisher Scientific	Cat. #11528926	Reducing agent
Gibco DMEM/F-12, with GlutaMAX	Fisher Scientific	Cat. #10565018	Medium
Gibco Glasgow's MEM	Fisher Scientific	Cat. #11570576	Medium
Gibco Glutamax	Fisher Scientific	Cat. #35050038	Medium Supplement
Gibco Neurobasal Medium	Thermo Fisher Scientific	Cat. #21103049	Medium
Mouse Methylcellulose Complete Medium	R&D Systems	Cat. #HSC007	Medium
Murine FGF-basic	PeproTech	Cat. #450-33	Peptide, recombinant protein
Murine Flt3-Ligand	PeproTech	Cat. #250-31L	Peptide, recombinant protein
Murine LIF	PeproTech	Cat. #250-02	Peptide, recombinant protein
Murine SCF	PeproTech	Cat. #250-03	Peptide, recombinant protein
Murine Sonic Hedgehog (Shh)	PeproTech	Cat. #315-22	Peptide, recombinant protein
Murine TPO	PeproTech	Cat. #315-14	Peptide, recombinant protein
Murine VEGF165	PeproTech	Cat. #450-32	Peptide, recombinant protein
N-2 Supplement (100x)	Thermo Fisher Scientific	Cat. #17502048	Medium supplement