Ethical approval, banking and processing of human samples
Primary human samples (PB or BM; described in Supplementary Table 1) were analyzed with approvals from the Inserm Institutional Review Board Ethical Committee (project C19-73, agreement 21-794, CODECOH DC-2020-4324) and from the INForMeD Study (REC: 199833, 26 July 2016, University of Oxford). Patients and normal donors provided written informed consent in accordance with the Declaration of Helsinki for sample collection and use in research. For secondary AML patients, we specifically selected samples from patients with known TP53 mutation.
Cells were subjected to Ficoll gradient centrifugation and for some samples, CD34 enrichment was performed using immunomagnetic beads (Miltenyi). Total mononuclear cells (MNCs) or CD34+ cells were frozen in FBS supplemented with 10% dimethyl sulfoxide for further analysis.
Targeted bulk sequencing
Bulk genomic DNA from patient samples’ mononuclear or CD34+ cells was isolated using DNeasy Blood & Tissue Kit (Qiagen) or QIAamp DNA Mini Kit (Qiagen) as per the manufacturer’s instructions. Targeted sequencing was performed using a TruSeq Custom Amplicon panel (Illumina) or a Haloplex Target Enrichment System (Agilent Technologies) with amplicons designed around 32, 44 or 77 genes53. Targets were chosen based on the genes/exons most frequently mutated and/or likely to alter clinical practice (diagnostic, prognostic, predictive or monitoring capacity) across a range of myeloid malignancies (for example, MDS/AML/MPN). Targets covered in all panels include ASXL1, CALR, CBL, CEBPA, CSF3R, DNMT3A, EZH2, FLT3, HRAS, IDH1, IDH2, JAK2, KIT, KRAS, MPL, NPM1, NRAS, PHF6, RUNX1, SETBP1, SF3B1, SRSF2, TET2, TP53, U2AF1, WT1 and ZRSR2. Sequencing was performed with a MiSeq sequencer (Illumina), according to the manufacturer’s protocols. Raw sequence data in FASTQ format were analyzed using the following variant callers and as previously described16,53: BWA v-0.7.12 (read alignment); Picard-tools (marking duplicates); samtools v-1.2; v-1.139 (BAM file creation); GATK HaplotypeCaller v-3.4-46 GRVC v-1.1; snpEff v-4.0 (variant calling). Run quality control included %DP_100X (>95%), %DP_200X (>90%), number of reads per sample and % reads q30 forward and reverse (>85%), read quality mean (>30) and percentage of mapped reads (>75%). A minimum of ten reads was required for variant calling. Results were analyzed after alignment of the reads using two dedicated pipelines, SOPHiA DDM (Sophia Genetics) and an in-house software GRIO-Dx. All pathogenic variants were manually checked using Integrative Genomics Viewer software. The analysis is presented in Extended Data Figs. 1a and 8a.
Pathogenic scores for each TP53 variant (Extended Data Fig. 8e) were derived from the Catalog of Somatic Mutations in Cancer using the FATHMM-MKL algorithm. The FATHMM-MKL algorithm integrates functional annotations from ENCODE with nucleotide-based hidden Markov models to predict whether a somatic mutation is likely to have functional, molecular and phenotypic consequences. Scores greater than 0.7 indicate that a somatic mutation is likely pathogenic, while scores less than 0.5 indicate a neutral classification.
The type and location of TP53 mutations from this study, de novo AML patients and CHIP individuals represented in Extended Data Fig. 8f were generated using Pecan Portal54. De novo AML TP53 mutations were downloaded from ref. 55 and ref. 27; CHIP-associated TP53 mutations were obtained from refs. 56,57,58.
Sanger sequencing of patient-associated mutations in PDX models
Genomic DNA from PDX sorted populations (LMPP: hCD45+Lin–CD34+CD38−CD45RA+CD90− and GMP: hCD45+Lin–CD34+CD38+CD45RA+CD123+) was extracted using QIAamp DNA Mini Kit (Qiagen). Sanger sequencing was performed with forward or reverse primers (Supplementary Table 6a) targeting mutations identified by targeted bulk sequencing in the corresponding primary samples using Mix2seq kit (Eurofins Genomics) and sequences were analyzed with the ApE editor.
SNP array sample preparation, copy number variant and loss of heterozygosity analysis
Bulk genomic DNA from patients’ MNCs was isolated using DNeasy Blood & Tissue Kit (Qiagen) as per the manufacturer’s instructions. 250 ng of gDNA was used for hybridization on an Illumina Infinium OmniExpress v1.3 BeadChips platform.
To call mosaic copy number events in primary patient samples, genotyping intensity data generated were analyzed using the Illumina Infinium OmniExpress v1.3 BeadChips platform. Haplotype phasing, calculation of log R ratio (LRR) and B-allele frequency (BAF), and calling of mosaic events were performed using MoChA WDL pipeline v2021-01-20 (MoChA: a BCFtools extension to call mosaic chromosomal alterations starting from phased VCF files with either BAF and LRR or allelic depth) as previously described59,60. In brief, MoChA comprises the following steps: (1) filtering of constitutional duplications; (2) use of a parameterized hidden Markov model to evaluate the phased BAF for variants on a per-chromosome basis; (3) deploying a likelihood ratio test to call events; (4) defining event boundaries; (5) calling copy number and (6) estimating the cell fraction of mosaic events. A series of stringent filtering steps were applied to reduce the rate of false positive calls. To eliminate possible constitutional and germline duplications, we excluded calls with lod_baf_phase <10, those with length <500 kbp and rel_cov >2.5, and any gains with estimated cell fraction >80%, log(R) > 0.5 or length <24 Mb. Given that interstitial LOH are rare and likely artefactual, all LOH events <8 Mb were filtered59. Events on genomic regions reported to be prone to recurrent artifact59 (chr6 < 58 Mb, chr7 > 61 Mb and chr2 > 50 Mb) were also filtered, and those where manual inspection demonstrated noise or sparsity in the array.
To find common genomic lesions on a focal and arm level, Infinium OmniExpress arrays were initially processed with Illumina Genome Studio v2.0.4. Following this, LRR data were extracted for all probes and array annotation was obtained from Illumina (InfiniumOmniExpress-24v1-3_A1). LRR data were then smoothed and segmentation called using the CBS algorithm from the DNACopy61,62 v1.60.0 package in R. A minimum number of five probes was required to call a segment, and segments were analyzed using GenomicRanges63 v1.38.0. Definitions of amplification, gain, loss and deletion events were as outlined in ref. 64. Segmentation data were then analyzed in GISTIC65 v2.023.
For PDX models, genomic DNA from sorted populations (LMPP: hCD45+Lin–CD34+CD38−CD45RA+CD90− and GMP: hCD45+Lin–CD34+CD38+CD45RA+CD123+) was extracted using QIAamp DNA Mini Kit (Qiagen). SNP-CGH array hybridization was performed using the Affymetrix Cytoscan HD (Thermo Fisher Scientific) according to the manufacturer’s recommendations. DNA amplification was checked using BioSpec-nano spectrophotometer (Shimadzu) with expected concentrations between 2,500 ng μl−1 and 3,400 ng μl−1. DNA length distribution post fragmentation was checked using D1000 ScreenTapes on Tapestation 4200 instrument (Agilent Technologies). Cytoscan HD array includes 2.6 million markers combining SNP and nonpolymorphic probes for copy number evaluation. Raw data CEL files were analyzed using the Chromosome Analysis Suite software package (v4.1, Affymetrix) with genome version GRCh37 (hg19) only if achieving the manufacturer’s quality cut-offs. Only CNAs >10 kb were reported in the analysis presented in Extended Data Fig. 3k,l.
Single-molecule cloning and sequencing of patient-derived cDNA
To experimentally verify the bi-allelic nature of TP53 mutations in TP53-sAML patients, cDNA from a selected patient with putative TP53 bi-allelic status (patient ID GR004) was PCR-amplified using cDNA-specific primers spanning both TP53 mutations (fwd: 5′-GACCCTTTTTGGACTTCAGGTG-3′ and rev: 5′-CCATGAGCGCTGCTCAGATAG-3′). PCR amplification was performed with KAPA 2X Ready Mix (Roche), a Taq-derived enzyme with A-tailing activity, for direct cloning into a TA vector (pCR2.1 TOPO vector, TOPO TA Cloning Kit, Invitrogen) as per the manufacturer’s instructions. Sanger sequencing for 10 different colonies was performed using M13 forward and reverse primers; a representative example is shown in Extended Data Fig. 1h.
Fluorescence-activated cell sorting (FACS) and single-cell isolation
Single-cell FACS-sorting was performed as previously described16, using BD Fusion I and BD Fusion II instruments (Becton Dickinson) for 96-well plate experiments or bulk sorting experiments, and SH800S or MA900 (SONY) for 384-well plate experiments. Experiments involving the isolation of human HSPCs included single color stained controls (CompBeads, BD Biosciences) and Fluorescence Minus One controls (FMOs). Antibodies used for HSPC staining are detailed in Supplementary Table 7a (combinations indicated as Panel A or B).
Briefly, single cells were directly sorted into 384-well plates containing 2.07 μl of TARGET-seq lysis buffer66. Lineage-CD34+ cells were indexed for CD38, CD90, CD45RA, CD123 and CD117 markers, which allowed us to record the fluorescence levels of each marker for each single cell. The 7-aminoactinomycin D (7-AAD) was used for dead cell exclusion. Flow cytometry profiles of the human HSPC compartment (Extended Data Figs. 2 and 9) were analyzed using FlowJo software (version 10.1, BD Biosciences).
Single-cell TARGET-seq cDNA synthesis
Reverse transcription (RT) and PCR steps were performed as previously described66. Briefly, SMARTScribe (Takara, 639537) retrotranscriptase, RNAse inhibitor (Takara, 2313A) and a template-switching oligo were added to the cell lysate to perform the retrotranscription step. Immediately after, a PCR mix comprised of SeqAMP (Takara, 638509) and ISPCR primer (binding to a common adapter sequence in all cDNA molecules) was used for the PCR step with 24 cycles of amplification. Target-specific primers spanning patient-specific mutations were also added to RT and PCR steps (Supplementary Table 6a). After cDNA synthesis, cDNA from up to 384 single-cell libraries was pooled, purified using Ampure XP Beads (0.6:1 beads to cDNA ratio; Beckman Coulter) and resuspended in a final volume of 50 μl of EB buffer (Qiagen). The quality of cDNA traces was checked using a high-sensitivity DNA kit in a Bioanalyzer instrument (Agilent Technologies).
Whole transcriptome library preparation and sequencing
Pooled and bead-purified cDNA libraries were diluted to 0.2 ng μl−1 and used for tagmentation-based library preparation using a custom P5 primer and 14 cycles of PCR amplification66. Each indexed library was purified twice with Ampure XP beads (0.7:1 beads to cDNA ratio), quantified using Qubit dsDNA HS Assay Kit (Invitrogen, Q32854) and diluted to 4 nM. Libraries were sequenced on a HiSeq4000, HiSeqX or NextSeq instrument using a custom sequencing primer for read1 (P5_seq: GCCTGTCCGCGGAAGCAGT GGTATCAACGCAGAGTTGC*T, PAGE purified) with the following sequencing configuration: 15 bp R1; 8 bp index read; 69 bp R2 (NextSeq) or 150 bp R1; 8 bp index read; 150 bp R2 (HiSeq).
TARGET-seq single-cell genotyping
After RT-PCR, cDNA + amplicon mix was diluted 1:2 by adding 6.25 μl of DNAse/RNAse free water to each well of 384-well plate. Subsequently, a 1.5 μl aliquot from each single-cell derived library was used as input to generate a targeted and Illumina-compatible library for single-cell genotyping66. In the first PCR step, target-specific primers containing a plate-specific barcode (Supplementary Table 6b) were used to amplify the target regions of interest. In a subsequent PCR step, Illumina compatible adapters (PE1/PE2) containing single-direction indexes (Access Array Barcode Library for Illumina Sequencers-384, Single Direction, Fluidigm) were attached to pre-amplified amplicons, generating single-cell barcoded libraries. Amplicons from up to 3,072 libraries were pooled and purified with Ampure XP beads (0.8:1 ratio beads to product; Beckman Coulter). These steps were performed using Biomek FxP (Beckman Coulter), Mosquito (TTP Labtech) and VIAFLO 96/384 (INTEGRA Biosciences) liquid handling platforms. Purified pools were quantified using Qubit dsDNA HS Assay Kit (Invitrogen, Q32854) and diluted to a final concentration of 4 nM. Libraries were sequenced on a MiSeq or NextSeq instrument using custom sequencing primers as previously described66 with the following sequencing configuration: 150 bp R1; 10 bp index read; 150 bp R2.
Targeted single-cell genotyping analysis
Data preprocessing
For each cell, the FASTQ file containing both targeted gDNA and cDNA-derived sequencing reads was aligned to the human reference genome (GRCh37/hg19) using Burrows–Wheeler Aligner (BWA v0.7.17) and STAR67 (v2.6.1d). Custom perl scripts were used to demultiplex the gDNA and mRNA reads in the BAM file into separate SAM files based on targeted-sequencing primer coordinates (https://github.com/albarmeira/TARGET-seq). Next, Samtools68 (v1.9) was used to concatenate the BAM header to the resulting SAM files before reconverting the SAM file to BAM format, which was subsequently sorted by genomic coordinates and indexed. Both gDNA and mRNA reads were tagged with the cell’s unique identifier using Picard (v2.3.0) ‘AddOrReplaceReadGroups’ and duplicate reads were subsequently marked using Picard ‘MarkDuplicates’. The sequencing reads overhanging into intronic regions in the mRNA reads were additionally hard-clipped using GATK (v4.1.2.0) SplitNCigarReads69,70.
Variant calling
Variants were called from the processed BAM files using GATK Mutect2 with the options (–tumor-lod-to-emit 2.0 –disable-read-filter NotDuplicateReadFilter –max-reads-per-alignment-start) to increase the sensitivity of detecting low-frequency variants. The frequency of each nucleotide (A, C, G, T) and indels at each predefined variant site were also called using a Samtools mpileup as previously described16. Lastly, the coverage at each predefined variant site was computed using Bedtools71 (v2.27.1).
To determine the coverage threshold of detection for each variant site, the coverage for ‘blank’ controls (empty wells) was first tabulated. A cut-off coverage outlier value was computed as having a coverage exceeding 1.5 times the length of the interquartile range from the 75th percentile. Next, a value of 30 was added to this outlier value to yield the final coverage threshold to be used for genotype assignment.
Genotype assignment
For each predefined variant site, the number of reads representing the reference and alternative (variant) alleles for indels (insertion and deletions) and single nucleotide variants (SNVs) were tabulated from the outputs of GATK Mutect2 and Samtools mpileup, respectively.
Here a genotype scoring system was introduced to assign each variant site into one of the following three possible genotypes: WT, heterozygous or homozygous mutant. Chi-square (\({{\rm{\chi }}}^{2}\)) test was first used to compare the observed frequency of reference and alternative alleles against the expected fraction of reference and alternative alleles corresponding to the three genotypes. The expected fraction of the reference alleles was 0.999, 0.5 and 0.001, and the expected fraction of the alternative alleles was 0.001, 0.5 and 0.999 for WT, heterozygous and homozygous mutant genotype, respectively. The \({{\rm{\chi }}}^{2}\) statistics were then tabulated for each fitted model and converted to genotype scores using the following formula:
$${\rm{Score}}_{\rm{genotype}}=\,\frac{1}{{\rm{log}}_{10}({\chi }^{2}+1)}$$
The genotype assigned to the variant site was based on the genotype model with the highest score.
Next, the variant (alternative) allele frequency (VAF) was computed and variant sites with 2 < VAF < 4 and 96 < VAF < 98 were reassigned as ‘ambiguous’. For cells with no variants detected at the specific variant sites by the mutation callers (either due to the absence of the variants, that is WT genotype, or that such variants were present below the detection limit), a ‘WT’ genotype was assigned to those cells with coverage above the specific threshold and ‘low coverage’ to those cells with coverage below such threshold.
Taken together, each variant site was assigned one of the five following genotypes: WT, heterozygous, homozygous mutant, ambiguous or low coverage. Variants with ambiguous or low coverage assignments for a particular cell were excluded from the analysis.
Computational reconstruction of clonal hierarchies
Genotypes for each single cell were recoded for input to SCITE in a manner inspired by ref. 72; each mutation in each gene was coded as two loci, representing two different alleles. In the first recorded loci, all homozygous calls from each mutation where coded as heterozygous genotype calls. In the second recorded loci, all heterozygous and homozygous genotype calls in the original mutation matrix were coded as homozygous reference (that is, WT) and heterozygous, respectively. For example, if for a certain mutation 0 represents WT status, 1 encodes heterozygous and 2 refers to homozygous status, these would be encoded as (0,0), (1,0) and (1,1), respectively, where the first term in the parenthesis corresponds to the first loci and the subsequent to the second loci.
Then, SCITE was used (git revision 2016b31, downloaded from https://github.com/cbg-ethz/SCITE.git; ref. 73) to sample 1,000 mutation trees from the posterior for every single-cell genotype matrix corresponding to a particular patient, where all possible mutation trees are equally likely a priori. For patients in which several disease time points were available, all time points were merged for SCITE analysis. As parameters for every SCITE run ‘-fd 0.01’ (corresponding to the allelic dropout (ADO) rate of reference allele in our adapted SCITE model), ‘-ad 0.01’ (corresponding to the ADO of the alternate allele), a chain length (-l) of 1e6 and a thinning interval of 1 while marginalizing out cell attachments (-p 1 -s) were used.
To summarize the posterior tree sample distribution, the number of times a particular sample matched each tree was computed. For each patient, the most common tree topology in the posterior tree samples is reported (Extended Data Figs. 2b–o and 9e–m), where ‘pp’ is the proportion of samples that match this tree. For each clade in the most common posterior tree, clade probabilities were estimated as the proportion of trees in the posterior that contained the clade. These are indicated in each square for each mutation in Extended Data Figs. 2b–o and 9e–m.
Clone assignment
For every patient’s most common posterior tree, we assigned every cell to the tree node that matches the genotype of that particular cell. If an exact match was not found, then for every tree node, the loss of assigning a cell to that node was calculated using the following loss function:
$$\begin{array}{l}l(M)={\rm{log }}\left(\text{ADO}\right)\left(M\left[1,2\right]+M\left[3,2\right]\right)+{\rm{log }}\left({FD}\right)\left(M\left[2,1\right]+M\left[2,3\right]\right) \\ +{\rm{log }}({A\text{DO}}^{2}\text{FD})(M\left[1,3\right]+M\left[3,1\right])\end{array}$$
where M is a confusion matrix generated across all loci of a cell in which the first index represents the genotype that was measured for that particular cell (1 = homozygous reference, 2 = heterozygous, 3 = homozygous alternate), and the second index represents the genotype implied by the tree node. ADO = 0.01 and FD = 0.001 were used. Every cell was assigned to the node with the lowest loss \(l\). For the trees presented in Extended Data Figs. 2b–o and 9e–m, only the numbers of cells with exact genotype matches were reported.
Testing for evidence of homozygous genotypes
Due to the nature of our loci-specific mutation encoding (each gene is encoded as two loci), homozygous mutations are placed in the clonal hierarchy independently of their accuracy. Therefore, for every patient and at every locus with observed homozygous alternate genotype calls, the tested null hypothesis was that all homozygous alternate genotype calls are due to ADO at a level not exceeding 0.05 using a one-tailed binomial test. The total number of draws for the test is the number of heterozygous and homozygous alternate genotype calls at the locus, the number of successful draws is the number of homozygous alternate calls and the success rate is 0.05. Only homozygous alternate genotype calls below this 0.05 cut-off were reported in Extended Data Figs. 2b–o and 9e–m; the results of the binomial test are reported for each patient and mutation in Supplementary Table 8.
Computational validation of TP53 bi-allelic status from single-cell targeted genotyping datasets
To further validate the bi-allelic status of TP53 mutations in our dataset, the patterns of ADO in TARGET-seq single-cell genotyping data from patients carrying at least two different TP53 mutations were investigated (n = 6; Extended Data Fig. 1j).
To test the hypothesis that the observed TP53-WT/TP53-homozygous (TP53-WT/HOM; or (0,2)) cells are the result of a chromosomal loss (and therefore, in different alleles), the following null hypothesis (H0) was formulated: observed TP53-WT/HOM cells are double ADO events. Under H0, every TP53-WT/HOM cell (0,2), TP53-HOM/WT cell (2,0), TP53-HOM/HOM (2,2) as well as an unknown number of TP53-WT/WT (0,0) are the result of a TP53-HET/HET (1,1) cell undergoing ADO at both sites. The following assumptions were made: (1) ADO is unbiased toward HOM or WT and (2) ADO events at each TP53 site are independent. The null hypothesis was then tested with a binomial test, where the number of (2,2) events should be half the sum of (0,2) + (2,0) events (Extended Data Fig. 1j). (0,0) events were disregarded.
If TP53 mutations are bi-allelic, the expected number of WT/HOM and HOM/WT would be higher than HOM/HOM cells considering TARGET-seq expected ADO rates (1–5%).
Single-cell 3′-biased RNA-sequencing data preprocessing
FASTQ files for each single cell were generated using bcl2fastq (version 2.20) with default parameters and the following read configuration: Y8N*, I8, Y63N*. Read 1 corresponds to a cell-specific barcode, index read corresponds to an i7 index sequence from each cDNA pool and read 2 corresponds to the cDNA molecule. PolyA tails were trimmed from demultiplexed FASTQ files with TrimGalore (version 0.4.1). Reads were then aligned to the human genome (hg19) using STAR (version 2.4.2a), and counts for each gene were obtained with FeatureCounts (version 1.4.5-p1; options–primary). Counts were then normalized by dividing each gene count by the total library size of each cell and multiplying this value by the median library size of all cells processed, as implemented in the ‘normalize_UMIs’ function from the SingCellaR package74 (version 1.2.1; https://github.com/supatt-lab/SingCellaR). A summary of the preprocessing pipeline can be found at https://github.com/albarmeira/TARGET-seq-WTA.
Quality control was performed using the following parameters: number of genes detected >500, percentage of ERCC-derived reads <35%, percentage of mitochondrial reads <0.25% and percentage of unmapped reads <75%. Cells with less than 2,000 reads in batch1, 5,000 reads in batch2 and 20,000 reads in batch3 were further excluded. This QC step was performed independently for each sequencing batch owing to differences in sequencing depth (mean library size: 42,949 batch 1, 93,580 batch 2 and 171,393 batch 3). After these QC steps, 7,123 cells passed QC for batch 1, 5,779 for batch 2 and 6,319 for batch 3 (79.3%, 68.9% and 80.3% of cells processed, respectively). Then, 2,734 cells from a previously published study16 corresponding to 8 MF patients and 2 normal donor controls were further integrated, encompassing a final dataset of 21,955 cells in total.
Identification of highly variable genes
Highly variable genes above technical noise were identified by fitting a gamma generalized linear model (GLM) of the log2(mean expression level) and coefficient of variation for each gene, using the ‘get_variable_genes_by_fitting_GLM_model’ from SingCellaR package and the following options: mean_expr_cutoff = 1, disp_zscore_cutoff = 0.1, quantile_genes_expr_for_fitting = 0.6 and quantile_genes_cv2_for_fitting = 0.2. Those genes with a coefficient of variation above the fitted model and expression cut-off were selected for further analysis, excluding those annotated as ribosomal or mitochondrial genes.
CNA inference from single-cell transcriptomes
InferCNV (v1.0.4) was used to identify CNAs in single-cell transcriptomes75 (https://github.com/broadinstitute/inferCNV/wiki). Briefly, inferCNV creates genomic bins from gene expression matrices and computes the average level of expression for each of these bins. The expression across each bin is then compared to a set of normal control cells, and CNAs are predicted using a hidden Markov model. For each patient, inferCNV was performed with the following parameters: ‘cutoff = 0.1, denoise = T, HMM = T’, compared to the same set of normal donor control cells (n = 992). To identify CNA subclones, inferCNV in analysis_mode = ‘subclusters’ was used. CNAs identified by inferCNV were manually curated by removing those with size <10 kb, merging adjacent CNA calls with identical CNA status into larger CNA intervals and comparing them to SNP-Array bulk CNA calls. Finally, to generate combined TARGET-seq single-cell genotyping and CNA-based clonal hierarchies, the CNA status from each inferCNV cluster was assigned to its predominant genotype.
Dimensionality reduction, data integration and clustering
PCA was performed using ‘runPCA’ function from the SingCellaR R package, and Force-directed graph analysis was subsequently performed using the ‘runFA2_ForceDirectedGraph’ with the top 30 PCA dimensions to generate the plots in Extended Data Fig. 4a.
For the analysis of patient IF0131 presented in Extended Data Fig. 3m, PCA was performed using ‘runPCA’ function from the SingCellaR R package and then UMAP was performed using the ‘runUMAP’ function with the top ten PCA dimensions and the following options: n.neighbors = 20, uwot.metric = ‘correlation’, uwot.min.dist = 0.30, n.seed = 1.
Integration of TARGET-seq single-cell transcriptomes from 10,459 cells corresponding to 14 TP53-sAML samples was performed using ‘runHarmony’ function implemented in the SingCellaR package, using the patient ID as covariate and the following options: n.dims.use = 20, harmony.theta = 1, n.seed = 1. Diffusion map analysis was performed using ‘runDiffusionMap’ with the integrative Harmony embeddings and the following parameters: n.dims.use = 20, n.neighbors = 5, distance = ‘euclidean’. Signature scores were calculated using ‘plot_diffusionmap_label_by_gene_set’ to generate the plots in Figs. 2a and 3a. Only cells with assigned genotypes ‘TP53 multihit’ and ‘TP53-WT’ are shown.
Pseudotime trajectory analysis
Monocle3 (ref. 76; https://cole-trapnell-lab.github.io/monocle3/) was used to infer differentiation trajectories from single-cell transcriptomes. Raw UMI count matrix and clustering annotations were extracted from the SingCellaR object to build a Monocle3 ‘cds’ object. Next, we retrieved the first two components of the diffusion map (DC1 and DC2), and the ‘learn_graph’ function was used to calculate the trajectory on the two-dimensional diffusion map, using TP53-WT preleukemic cell cluster as the root node. Pseudotime was calculated using ‘order_cells’ function and overlayed on the diffusion map embeddings to generate the plot in Fig. 2b.
Differential expression analysis
Differentially expressed genes from TARGET-seq datasets were identified using a combination of nonparametric Wilcoxon test, to compare the expression values for each group, and Fisher’s exact test, to compare the frequency of expression for each group, as previously described17. Logged normalized counts were used as input for this comparison, including genes expressed in at least two cells. Combined P values were calculated using Fisher’s method and adjusted P values were derived using Benjamini–Hochberg procedure. Significance level was set at P-adjusted < 0.05. For the analysis presented in Extended Data Fig. 4b and Supplementary Table 2, the top 100 differentially expressed genes with log2(FC) > 0.3 and at least 20% expressing cells are shown. Differentially expressed genes identified between TP53-multihit versus TP53-WT cells were further assessed for the enrichment of known p53 target genes (337 curated p53 target genes from ref. 77) for the analysis presented in Extended Data Fig. 4c. We assessed the extent of overlap of these gene lists using the R package GeneOverlap. The overlapping genes were further assessed for the enrichment of p53-related pathways using the R package clusterProfiler.
For the analysis presented in Fig. 2k,l, only genes overexpressed in TP53 multihit cells and log2(FC) > 0.75 were included; for Fig. 4d, only those with log2(FC) > 1 were considered. Violin plots (Fig. 4e and Extended Data Fig. 9n) from selected differentially expressed genes were generated using ‘ggplot2’ package in R.
Gene-set enrichment analysis
For analysis involving <600 cells (Fig. 4c and Supplementary Table 5), GSEA was performed using GSEA software version 4.0.3 (www.gsea-msigdb.org/gsea/index.jsp) with default parameters and 1,000 permutations on the phenotype, using log2(normalized counts).
For analysis involving >600 cells per group (Fig. 3k and Extended Data Figs. 4d and 9o), GSEA was performed with ‘identifyGSEAPrerankedGene’ function from SingCellaR R package with default options. Briefly, differential expression analysis was performed between two cell populations using the Wilcoxon rank sum test, and the resulting P values were adjusted for multiple testing using the Benjamini–Hochberg approach. Before the differential expression analysis, down-sampling was performed so that both cell populations had the same number of cells. Next, −log10(P value) transformation was performed and the resulting P values were multiplied by +1 or −1 if the corresponding log2(FC) was >0.1 or <−0.1, respectively. The gene list was ranked using this statistic in ascending order and used as input for GSEA analysis using ‘fgsea’ function from the fgsea R package with default options.
MSigDB HALLMARK v7.4 50-gene sets or previously published signatures (https://www.gsea-msigdb.org/gsea/msigdb/cards/GENTLES_LEUKEMIC_STEM_CELL_UP) were used for all analysis. Normalized enrichment scores were displayed in a heatmap using pheatmap R package. Gene sets with false discovery rate (FDR) q value lower than 0.25 were considered significant.
Projection of single-cell transcriptomes
A previously published human hematopoietic atlas was downloaded from https://github.com/GreenleafLab/MPAL-Single-Cell-2019 and used as a normal hematopoietic reference to project TP53-sAML and de novo AML transcriptions using Latent Semantic Index Projection24. Common genes to all datasets were selected, and then TP53-sAML or previously published de novo AML cells25 were projected using ‘projectLSI’ function for the analysis presented in Fig. 2c,d. A previously published human MF atlas78 was used as a reference to project TP53-sAML multihit cells in the analysis presented in Extended Data Fig. 5d,e, using previously defined force-directed graph embeddings.
Velocyto analysis
Loom files were generated for each single cell using velocyto (v0.17.13) with options -c and -U, to indicate that each BAM represents an independent cell and reads are counted instead of molecules (UMIs), respectively79. The individual loom files were subsequently merged using the combine function from the loompy Python module.
HDs with at least 300 cells with RNA-sequencing data and patients with at least 300 cells consisting of >50 preleukemic (TP53 WT) cells and >50 TP53 multihit cells were included for analysis. For each individual, the Seurat object was created from the merged loom file and processed for downstream RNA-velocity analysis80. Specifically, for each patient, the spliced RNA counts were normalized using regularized negative binomial regression with the SCTransform function81. Next, linear dimension reduction was performed using RunPCA function and the first 30 principal components were further used to perform nonlinear dimension reduction using the RunUMAP function. Ninety-six multiple rate kinetics (MURK) genes previously shown to possess coordinated step-change in transcription and hence violate the assumptions behind scVelo were removed82. The processed and MURK gene-filtered Seurat object was then saved in h5Seurat format using the SaveH5Seurat function and finally converted to h5ad format using the ‘Convert’ function.
AnnData object was created from the h5ad file using the scvelo python module for RNA velocity analysis83. Highly variable genes were identified and the corresponding spliced and unspliced RNA counts were normalized and log2-transformed using the scvelo.pp.filter_and_normalize function. Next, the first- and second-order moments were computed for velocity estimation using the scvelo.pp.moments function. The velocities (directionalities) were computed based on the stochastic model as defined in the scvelo.t1.velocity function, and the velocities were subsequently projected on the UMAP embeddings generated from Seurat above. Finally, the UMAP embeddings were annotated using the HSPC and erythroid lineage signature scores74 and TP53 mutation status. For each cell, the cell lineage signature score was computed using the average SCTransform expression values of the individual cell lineage genes.
Analysis of bulk BeatAML and TCGA gene expression datasets
Data retrieval and preprocessing
Two publicly available AML cohorts with genetic mutation and RNA-sequencing data available were used to validate findings from our single-cell analysis, namely BeatAML26 and TCGA27. Gene expression values in FPKM (fragments per kilobase of transcript per million mapped reads) were retrieved from the National Cancer Institute (NIH) Genomic Data Commons (GDC)84. Gene expression values were then offset by 1 and log2 transformed. TP53 point mutation status was retrieved from the cBio Cancer Genomics Portal (cBioPortal)85. Clinical data including survival data for BeatAML and TCGA were retrieved from the BeatAML data viewer (Vizome) and NIH GDC, respectively.
We selected samples from the BeatAML cohort with an AML diagnosis (540 de novo AML and 96 secondary AML) collected within 1 month of the patient’s enrollment in the study, and with both TP53 mutation status and RNA-sequencing data available. For patients for whom multiple samples were available, samples were collapsed to obtain patient-level data. Specifically, the mean gene expression value for each gene from multiple samples was used to represent patient-level gene expression value. Furthermore, patients with at least one sample with a TP53 mutation were considered TP53-mutant. Analysis of TP53 VAF and reported karyotypic abnormalities indicated that the vast majority of patients could be classified as ‘multi-hit’, and therefore patients were classified as TP53-mutant or WT without taking into account TP53 allelic status. In total, 360 patients with TP53 mutation status (329 TP53 WT and 31 TP53 mutant) and RNA-sequencing data available were included for analysis. Of these, 322 patients had concomitant survival data available (294 TP53 WT and 28 TP53 mutant).
The TCGA cohort consisted of 200 de novo AML patients represented by one sample each, of which 151 patients had TP53 mutation status (140 TP53 WT and 11 TP53 mutant) and RNA-sequencing data available, and were included for analysis. Of these, 132 patients had concomitant survival data available (124 TP53 WT and 8 TP53 mutant).
Cell lineage gene signature scores
For each sample, a given cell lineage gene signature score was computed as the mean expression values of the individual genes belonging to the cell lineage gene signature. Here the gene signature scores for two cell lineages were computed, namely myeloid and erythroid populations. Two gene sets for each cell lineage were compiled. The first gene set was based on cell lineage markers previously reported in the literature, whereas the second gene set was based on cell lineage markers derived from analyzing a published single-cell dataset24. Genes from each score are described in Supplementary Table 3.
For the former approach, six erythroid genes (KLF1, GATA1, ZFPM1, GATA2, GYPA and TFRC; Fig. 2e and Extended Data Figs. 5k,m and seven myeloid genes (FLI1, SPI1, CEBPA, CEBPB, CD33, MPO and IRF8; Fig. 2f) were identified. For the latter approach, the expression values of erythroid and myeloid cell clusters were first compared separately against all other cell clusters using Wilcoxon ranked sum test. The erythroid cluster consisted of the early and late erythroid populations, while the myeloid cluster consisted of granulocyte, monocyte and dendritic cell populations. Erythroid and myeloid-specific gene signatures were defined as genes having FDR values <0.05 and log2(FC) > 0.5 in ≥20 and 17 comparisons, respectively. In total, 100 erythroid genes and 135 myeloid genes were identified from this single-cell dataset (Supplementary Table 3) and were used to compute the scores presented in Extended Data Fig. 5g–j.
TP53 target gene score
Genes downregulated in TP53-multihit compared to TP53-WT cells (defined as per ‘differential expression analysis’ section above; related to Extended Data Fig. 4b) and p53 targets positively regulated from ref. 77 were used to compute a TP53-target gene score presented in Extended Data Fig. 5k.
Prognostic signatures and Cox-regression survival models
LSC signature score
The 17-gene LSC17 gene set was retrieved from ref. 31. For each sample, the LSC17 score was defined as the linear combination of gene expression values weighted by their respective regression coefficients.
To identify TP53-sAML LSC signatures from our TARGET single-cell dataset, two different approaches were used. First, differentially expressed genes were identified as overexpressed in all Lin–CD34+ TP53-multihit cells regardless of their transcriptional classification (p53-all-cells) versus MF, HD and TP53-WT preleukemic cells; this gene set consists of 29 genes (Supplementary Table 4a). For the second approach, the same analysis was performed, but TP53-multihit cells transcriptionally defined as LSCs (falling in the LSC-like cluster; Fig. 2a, middle) were specifically selected; this gene-set is comprised of 51 genes (p53LSC; Supplementary Table 4a).
Next, lasso cox regression with tenfold cross-validation implemented in the glmnet R package (version 4.1-1) was used to identify p53-all-cells and p53-LSC genes that were associated with survival and to estimate their respective regression coefficients86. Specifically, Harrel’s concordance measure (C-index) was used to assess the performance of each fitted model during cross-validation. The best model was defined as the fitted model with the highest C-index. Subsequently, the coefficient for each gene estimated using the best model was used to compute the gene signature scores. Only genes with nonzero coefficient values were included in the final gene set. In total, 9 and 44 genes were retained from the p53-all-cells and p53-LSC gene sets, respectively. For each sample, the gene signature score for each gene set was defined as the linear combination of gene expression values weighted by their respective regression coefficients31,86. The list of p53-LSC and p53-all-cells gene signatures is provided in Supplementary Table 4b.
Survival analysis
For each gene expression signature, patients were first split using the median gene expression signature score. This resulted in two groups of patients, namely patients with high expression scores (greater than or equal to the median) and patients with low expression scores (lower than the median), exemplified in Extended Data Fig. 6a,b.
The Cox proportional hazards regression model implemented by the survival R package (version 3.5–5) was fitted to estimate the hazard ratio associated with each feature. The log-rank test was used to test the differences between survival curves. The features analyzed here were LSC17, p53-all-cells and p53-LSC signatures. Patients with low gene expression signature scores (below median) and patients with TP53 WT status were specified as the reference groups in the model. Kaplan–Meier curves were plotted using the survminer R package (version 0.4.9) to visualize the probability of survival and sample size at a respective time interval.
In vitro assays
Short-term liquid culture experiments
For short-term liquid culture differentiation experiments (Fig. 3j and Extended Data Fig. 7h,i), single cells from different Lineage-CD34+ HSPC populations (HSC: CD34+CD38–CD45RA–CD90+, MPP: CD34+CD38−CD45RA−CD90−, LMPP: CD34+CD38−CD45RA+ and more committed progenitors CD34+CD38+) were directly sorted into a 96-well tissue culture plate containing 100 μl of differentiation media: StemSpan (StemCell Technologies, 09650), 1% penicillin+streptomycin, 20% BIT9500 (StemCell Technologies, 9500), 10 ng ml−1 SCF (Peprotech, 300-07), 10 ng ml−1 FLT3L (Peprotech, 300-19), 10 ng ml−1 TPO (Peprotech, 300-18-10), 5 ng ml−1 IL3 (Peprotech, 200-03), 10 ng ml−1 G-CSF (Peprotech, 300-23), 10 ng ml−1 GM-CSF (Peprotech, 300-03), 1 IU ml−1 EPO (Janssen, erythropoietin alpha, clinical grade) and 10 ng ml−1 IL6 (Peprotech, 200-06).
For all liquid culture experiments, 50 μl of fresh 1× differentiation media was added on day 4. Readout was performed by flow cytometry after 12 d of culture using the antibodies detailed in Supplementary Table 7c (combination indicated as Panel D).
LTC-IC assay
Fifty cells from each Lin–CD34+ population (HSC; MPP; LMPP and CD38+) and donor type (HDs, MF and TP53-sAML) were sorted in triplicate. Cells were resuspended in 100 μl of MyeloCult (StemCell Technologies, H5150) supplemented with hydrocortisone (10−6 M; StemCell Technologies, 74142) and plated into an irradiated supportive stromal cell layer (5,000 SI/SI cells and 5,000 M2-10B4 cells per well) in a 96-well tissue-culture plate coated with Collagen type I (Corning, 354236).
The medium was changed weekly and after 6 weeks of culture, cells were washed in IMDM + 20% fetal calf serum (FCS) and plated into 1.4 ml of cytokine-rich methylcellulose (StemCell Technologies, MethoCult H4435). Colonies were scored 14 days later under an inverted microscope, and each colony was classified according to its morphology as BFU-E (Burst-forming unit erythroid), CFU-G (Colony Forming Unit granulocyte), CFU-GM (granulocyte-macrophage), CFU-M (macrophage) or CFU-GEMM (granulocyte, erythrocyte, macrophage and megakaryocyte). Selected colonies were used for cytospin and genotyping as outlined below.
LTC-IC colony genotyping
LTC-IC colonies were picked from methylcellulose media, washed, resuspended in 10 μl of PBS and transferred to individual wells in a 96-well PCR plate. 15 μl of lysis buffer (Triton X-100 0.4%, Qiagen Protease 0.1 AU ml−1) was added to each well, and samples were incubated at 56 °C for 10 min and 72 °C for 20 min. A 3 μl aliquot from each lysate was used as input to generate a targeted and Illumina-compatible library for colony genotyping. The preparation of single-cell genotyping libraries involves three PCR steps. In the first PCR step, target-specific primers spanning each mutation of interest are used for amplification (Supplementary Table 6a); in the second PCR step, nested target-specific primers (Supplementary Table 6b) attached to universal CS1/CS2 adapters (forward adapter—CS1: ACACTGACGACATGGTTCTACA; reverse adapter—CS2: TACGGTAGCAGAGACTTGGTCT) further enrich for target regions; and in the third PCR step, Illumina-compatible adapters containing sample-specific barcodes are used to generate sequencing libraries.
TP53 knockdown and differentiation of human CD34+ cells
shRNA sequence for p53 knockdown has been previously cloned into the lentiviral vector pRRLsin-PGK-eGFP-WPRE and validated87. Primary human CD34+ cells from patients with MPN (Supplementary Table 1) were infected twice with scramble (shCTL) or shTP53 with a multiplicity of infection of 15 and sorted 48 h later on CD34 and GFP expression. Cells were cultured in serum-free medium with a cocktail of human recombinant cytokines containing EPO (1 IU ml−1, Amgen), FLT3-L (10 ng ml−1, Celldex Therapeutics), G-CSF (20 ng ml−1, Pfizer), IL-6 (10 ng ml−1, Miltenyi), GM-CSF (5 ng ml−1, Peprotech), IL-3 (10 ng ml−1, Miltenyi), TPO (10 ng ml−1, Kirin Brewery) and SCF (25 ng ml−1, Biovitrum AB).
On day 12 of the culture, cells were stained with the antibodies detailed in Supplementary Table 7c (combination indicated as Panel C). DAPI was used for dead cell exclusion before acquisition on a FACSCanto II (BD Biosciences) instrument and on a BD FACS Diva software (version 8.0.2). Analysis of FACS data was performed using Kaluza (version 2.1, Beckman Coulter) software.
Quantitative real-time PCR in shRNA experiments
In TP53 knockdown experiments, RNA from either CD34+ cells sorted after transduction or bulk cells at day 12 of culture was extracted using Direct-Zol RNA MicroPrep Kit (Zymo Research) and reverse transcription was performed with SuperScript Vilo cDNA Synthesis Kit (Invitrogen). Quantitative RT–PCR was performed on a 7,500 real-time PCR Machine using SYBR-Green PCR Master Mix (Applied Biosystems). Expression levels were normalized to PPIA (housekeeping gene). Primers used are listed in Supplementary Table 6c.
Xenotransplantation
Purified CD34+ cells from AML patients were transplanted via retroorbital vein injection in sublethally irradiated (1.5 Gy) NOD.CB17-Prkdcscid IL2rgtm1/Bcgen mice (B-NDG, Envigo) (female, 8 weeks old, n = 1 for IF0131, n = 3 for GR001). All experiments were approved by the French National Ethical Committee on Animal Care (2020-007-23589). Blood cell counts were performed monthly by submandibular sampling of mice with blood chimerism assessed by flow cytometry using hCD34, hCD45 and mCD45 antibodies (Supplementary Table 7b, PDX PB panel). The following endpoints were applied: >50% of human blast cells in the blood, abnormalities of blood cell count (hemoglobin <7 g dl−1, platelets <150 × 109 l−1 or white blood cells >60 × 109 l−1), altered general conditions or >15% of weight loss. At sacrifice, BM was stained with the antibodies listed in Supplementary Table 7b (PDX BM panel) and HSPC fractions were sorted on an Influx Cell sorter (BD Biosciences).
Evaluation of cell morphology
Cell morphology from PDX models (Extended Data Fig. 3d) and in vitro LTC-IC cultures (Extended Data Fig. 7f) was assessed after cytospin of 50–100,000 cells onto a glass slide (5 min at 500 r.p.m.) and May–Grünwald Giemsa staining, according to standard protocols. Images were obtained using an AxioPhot microscope (Zeiss).
Mouse bone marrow chimeras and ethical approval
Trp53tm2Tyj Commd10Tg(Vav1-icre)A2Kio or Trp53tm2Tyj Tg(Tal1-cre/ERT)42-056Jrg C57/BL6 mice (hereafter referred to as Vav-iCre Trp53R172H/+ or SCL-CreERT Trp53R172H/+, respectively) and C57/BL6 WT mice used for BM chimera experiments were bred and maintained in accordance to UK and France Home Office regulations. All experiments carried out were performed under Project License P2FF90EE8 approved by the University of Oxford Animal Welfare and Ethical Review Body or under Project License no. 2020-007-23589, approved by the French National Ethical Committee on Animal Care. Trp53tm2Tyj (ref. 88), Commd10Tg(Vav1-icre)A2Kio (ref. 89; Jackson Laboratory, 008610) and Tg(Tal1-cre/ERT)42-056Jrg (ref. 90) have been previously described.
For in vivo experiments, two different chimera settings were used. For the first setting (Fig. 5a), 1 million BM cells from Vav-iCre Trp53R172H/+ CD45.1 mice and 1 million BM CD45.2 WT cells from competitor mice were transplanted intravenously into lethally irradiated (10 Gy total body irradiation, split dose) congenic CD45.2 mice. Male and female recipient CD45.2 mice were used and were 6–8 weeks old at transplantation, while male and female CD45.1 experimental BM donors were 5–6 weeks old at the time of BM collection. For the second setting (Fig. 5h), 0.9 million BM cells from Trp53LSL–R172H/+ CD45.2 mice (two males for two independent experiments, 8 and 13 weeks old) and 2.1 million BM CD45.1 WT competitor mice (four males for two independent experiments, 11 and 17 weeks old) were transplanted intravenously into lethally irradiated (9.5 Gy total body irradiation) congenic CD45.2 mice (females, 8 weeks old) and Trp53 mutation was induced 4 weeks after transplantation by tamoxifen (gavage 200 mg kg−1, Sigma-Aldrich) during 4 d, followed by tamoxifen feeding during 2 weeks (Ssniff Diet). In each cohort, a selection of mice was injected intraperitoneally with three rounds of six injections, each of 200 μg poly(I:C) (first setting) or 100 μg poly(I:C) (second setting; GE Healthcare, 27-4732-01) or placebo (PBS1X). Alternatively, Vav-iCre Trp53R172H/+ mice were injected with three rounds of eight injections, each of 35 μg LPS from Escherichia Coli O111:B4 (LPS; L4391-1MG and L5293-2ML; Sigma-Aldrich).
Poly(I:C) and LPS were administered during weeks 6–8, 10–12 and 14–16 (setting 1), or during weeks 7–8, 11–12 and 15–16 (setting 2) post-transplantation. Within each round, injections were spaced one or two days apart. Blood cell counts and analysis of PB chimerism along with mature lymphoid and myeloid populations were performed every 2–4 weeks by submandibular sampling of mice, while BM chimerism and HSPC populations were analyzed 18–20 weeks after transplantation. The antibodies used are detailed in Supplementary Table 7d. For dead cell exclusion, 7-AAD (Sigma-Aldrich) or DAPI (BD Biosciences) were used. FACS analyses were carried out on BD Fortessa or BD Fortessa X20 (BD Biosciences) and profiles were later analyzed using FlowJo (version 10.1, BD Biosciences) or Kaluza (version 2.1, Beckman Coulter) software.
LSK apoptosis and cell cycle
BM LSK cells (setting 2) were stained with Annexin-V and DAPI in Annexin V binding buffer 1X (BD Biosciences) for apoptosis analysis. BM LSK cell cycle was assessed by flow cytometry using Ki-67 and DAPI staining, after fixation and permeabilization (BD Cytofix/Cytoperm and Permeabilization Buffer Plus, BD Biosciences).
M-FISH
Fifty CD45.1 (Trp53R172H/+) or CD45.2 (WT) LSK (Lin–Sca1+c-Kit+) cells from poly(I:C)-treated and control recipient mice were sorted and cultured for 1 week into Complete X-vivo15 media (BE-04-418Q, Lonza) supplemented with 10% FCS (Sigma-Aldrich, F9665), 0.1 mM 2-mercaptoethanol (Gibco, 21985023), 1% penicillin-streptomycin (PAA laboratories), 2 ng ml−1 mouse stem cell factor (mSCF; PeproTech, 250-03), 10 ng ml−1 mouse granulocyte–monocyte colony-stimulating factor (mGM–CSF; Immunex), 5 ng ml−1 human thrombopoietin (hTPO; PeproTech, 300-18-10), 10 ng ml−1 human granulocyte colony-stimulating factor (hG-CSF; Neopogen) 5 ng ml−1 human FLT3 ligand (hFL; Immunex, 300-19), 5 ng ml−1 mouse interleukin 3 (mIL-3; PeproTech, 213-13). Cells were cultured at 37 °C 5% CO2. On day 7 of culture, metaphase spreads were collected following synchronization with Colcemid (KaryoMAX; Thermo Fisher Scientific, 11519876) 50 ng ml−1, for 3 h at 37 °C. The cells were then incubated with KCl 75 mM for 15 min at 37 °C and spun down. Following this, the cells were fixed in a methanol-acetic acid and then dropped onto glass slides.
M-FISH was performed with the 21XMouse-Multicolor FISH probe kit (Metasystem Probes, D-0425-060-DI), following the manufacturer’s instructions. For microscopy analysis, slides were mounted in Vectashield Antifade Mounting Medium with DAPI (2BScientific, H-1200). Images were acquired and analyzed using Leica Cytovision software (v7.3.1), on an Olympus BX-51 epifluorescence microscope equipped with a JAI CVM4+ progressive-scan 24 fps B&W fluorescence CCD camera. All cells were karyotyped, excluding metaphases severely damaged for technical reasons.
The analysis of the M-FISH hybridized cells was blinded. The cells on each slide were scored for the presence of structural aberrations (translocations, and/or derivative chromosomes and fragments) and/or numerical abnormalities. The presence of more than 40 chromosomes per cell was considered a numerical abnormality, except for cases where it could clearly be attributed to the presence of adjacent metaphases. Chromosome counts lower than 40 were not scored as numerical abnormalities for the impossibility to rule out technical issues (that is, metaphases bursting at the hypotonic step). We scored as follows: translocations and presence of one chromosome plus one or more extra chromosomal fragment(s)/derivative(s) as ‘structural abnormalities’ (except for sex chromosomes); presence of two chromosomes (or one in case of sex chromosomes) plus one or more extra chromosomal fragment(s)/derivative(s) as ‘partial chromosome gains’; two chromosomes (or one in case of sex chromosomes) plus one or more extra chromosomes as ‘whole chromosome gains’; two chromosomes plus two chromosomes with at least five different chromosomes present in number = 4n as ‘tetraploidy or sub-tetraploidy’. Counts of numbers of karyotypic aberrations per cell were performed scoring every type of event occurring on one chromosome as a single event (that is, presence of four chromosomes is counted as one aberration).
IFNγ ELISA assay
WT mice were injected intraperitoneally with a single dose of 200 μg poly(I:C) and spleens were collected from injected mice and nontreated controls 4 h after injection. Spleens were processed into a single-cell suspension in 200 μl PBS, spun down at 500 g for 5 min and supernatant was collected and used as spleen serum. IFNγ levels were assessed using mouse IFNγ Quantikine ELISA assay (R&D Systems, MIF00) following the manufacturer’s instructions. Optical densities of 450 nm and 540 nm were determined using Clariostar microplate reader (BMG Labtech).
Statistical analysis
Statistical analyses are detailed in figure legends (Figs. 2–6 and Extended Data Figs. 4–10) and performed using GraphPad Prism software (7 or later version) or R (version 3.6.1 and 4.0.5) software. The number of independent experiments, donors and replicates for each experiment are detailed in figure legends.
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.