This was a prospective observational cohort study that included consecutive de novo paediatric (≤18 years) patients with AML registered from July 2016 to December 2019 at medical oncology outpatient clinic of our cancer centre. The workflow of the study is depicted in Figure 1 . Study was ethically approved by institute ethics committee and informed consent was taken from care givers and assent was obtained from all participants (≥8 years). Patients with granulocytic sarcoma without marrow involvement, acute promyelocytic leukaemia (AML M3), and mixed phenotypic acute leukaemia were excluded. Fifty age-matched patients of solid malignancies without marrow involvement were also enrolled as controls. Baseline clinical details, European LeukemiaNet (ELN) risk stratification (2), were recorded and all patients were treated uniformly as per institutional protocol (Methods S1 and S2) (13). Remission status and survival outcomes were noted.

Bone marrow mononuclear cells were isolated by Ficoll-Hypaque (Sigma diagnostics, USA) density gradient centrifugation followed by isolation of total RNA and DNA (Methods S3). All the samples were assessed for mtDNA copy number as per previously described protocol (Methods S4) and classified into three separate groups based on relative mtDNA copy number (12). Patients were categorized into: AMLCN_H (mtDNA copy number ≥ 75^th percentile), AMLCN_I (mtDNA copy number 50th to 75th percentile) and AMLCN_L (mtDNA copy number< 50th percentile) groups. A subset of samples was randomly selected from each of the three sub-groups and controls with RNA integrity score above 7 and a total of 15 samples (12 patients including 3 from AMLCN_H group, 4 from AMLCN_I group, 5 from AMLCN_L group and 3 controls) were sent for whole transcriptome profiling for the identification of DEGs compared to controls (Methods S5 and S6; Figure S1A ).Absolute fold change value ≥ 2 (a ≥ two-fold change in expression, either upregulated or downregulated) and adjusted p value (q ≤ 0.05) threshold compared to controls was considered as differentially expressed genes (DEGs). The sequencing raw data was submitted to NCBI SRA (Sequence Read Archive) and available at PRJNA778747.

Out of all identified DEGs from transcriptome sequencing, mitochondria-related genes were filtered using Cytoscape compartment mitochondrion score (0 being minimum and 5 being highest) (14). DEGs with topmost mitochondrial compartment score were selected for validation in a cohort of paediatric patients with AML. Along with this, Hub genes as well as maximum interactive genes were identified using CytoHubba and molecular complex detection (MCODE) clustering algorithm respectively (15, 16). The genes of MCODE cluster 1 and Hub genes were assessed for their mitochondrial localization as above and genes in each group with highest mitochondrial compartment score were selected for validation (Methods S7). Based on these selection strategies, a total of 20 mitochondria-related DEGs were selected for validation.

Real time PCR was performed to validate the selected mitochondria-related genes using specific primers ( Table S1 ) and the gene expressions were quantified per previously described protocol (12).

For external validation of mitochondrial related DEGs, the RNA-sequencing data (Illumina HiSeq 2000) of TCGA adult AML(LAML) dataset was chosen, which is one of the largest datasets of transcriptomic profile in AML with recorded clinical outcome(https://www.cbioportal.org/study/summary?id=laml_tcga). The adult dataset was specifically chosen to see the impact of prognostic impact of the validated mitochondria-related age group in a different age group as well. The expression of validated DEGs was compared with LAML data set using online available GEPIA2 (Gene Expression Profiling Interactive Analysis) web server (http://gepia2.cancer-pku.cn/#index) (17).

The role of demographic and clinical features, including gender, age, haemoglobin, hyperleukocytosis (≥50000/µl), platelet count, presence of chloroma and ELN risk stratification (2)on survival outcome was analysed using the Cox regression. Factors with p<0.1 in univariable analysis were included for multivariable Cox regression in a forward stepwise manner using log likelihood change. Clinico-demographic factors which were significant in multivariable analysis were included in a multivariable Cox regression model along with gene signature risk score to explore the independent predictive value of gene signature. The timed AUC using 12-months survival and 18-months survival as the outcome and Harrel’s C-index of the clinical prognostic model and combined clinical and gene signature prognostic model were compared for identifying the additional prognostic benefits of gene signature over clinical parameters. The impact of mtDNA copy number on survival outcome was also analysed similarly.

The prognostic impact of our gene signature risk score on OS was done in TCGA LAML (n=179) dataset by Cox regression analysis. Patients were similarly sub-grouped into high-risk and low-risk category based on median value of the gene signature; the survival outcomes of the patients were compared between high-risk score vs low-risk score patients using Kaplan Meier analysis and timed AUC of 12-month and 18-month survival was evaluated. Based on available karyotyping data, patients of the TCGA dataset were grouped into poor-risk karyotype and others (including good and intermediate-risk karyotype). The association of risk score with clinical features such as age, sex, and karyotype were also evaluated in TCGA dataset. The karyotype category and mitochondrial gene risk score were assessed for their impact on OS by a multivariable Cox regression model to explore the independent predictive value of the gene signature in the external cohort as well.

Total 170 patients were enrolled, out of which 27 patients (5 patients were AML M3, 4 had granulocytic sarcoma without marrow involvement, and 18 patients had insufficient samples) were excluded. The baseline demographic and clinical characteristics of final 143 patients are summarized in Table S2 . Median age was 10 years (range: 0.8-18 years) and 50% of the patients were classified as ELN good risk. Total 104 patients (72.7%) achieved complete remission (CR) after induction therapy. At median follow-up of 36 months (32.67-39.33 months), the median OS was 21.93 months (13.54–30.31months). The clinical characteristics of the TCGA LAML dataset are summarized in Table S3 .

We identified 898, 769, and 953 significantly dysregulated transcripts in AMLCN_H, AMLCN_I and AMLCN_L groups respectively by whole transcriptome sequencing as represented in volcano plots ( Figures S1B–D ). Majority of genes were found significantly downregulated in all three groups whereas the number of dysregulated genes were higher in AMLCN_H group compared to other two groups. A total of 351 DEGs (59 upregulated and 292 downregulated) were identified in AMLCN_H. Similarly, AMLCN_I and AMLCN_L groups had 290 (66 upregulated and 224 downregulated) and 332 (47 upregulated and 285 downregulated) DEGs respectively as compared to controls.

Out of all DEGs, 78, 58, and 71 mitochondria-related DEGs were identified in AMLCN_H, AMLCN_I and AMLCN_L groups respectively. Among them, 35 genes were common in all three subgroups, whereas 18, 12 and 14 mitochondria-related DEGs were exclusively present in AMLCN_H, AMLCN_I and AMLCN_L groups respectively ( Figure S1E ). In AMLCN_H, AMLCN_I and AMLCN_L groups, we identified 17, 18 and 17 hub genes respectively using CytoHubba analysis, of which eight were common among all subgroups ( Table S4 , Figures S1F–H ). Furthermore, using MCODE analysis, clusters with maximum scores were generated and seed gene was determined in the three groups ( Figures S1I–K ). MMP9 was identified as seed node with maximum MCODE score in both AMLCN_H and AMLCN_I group ( Table S5 ). Based on the mitochondrial compartment score, CytoHubba and MCODE analyses, a total of 20 DEGs were selected for further validation ( Table S6 ). The expression pattern of these selected DEGs in RNA sequencing data were represented in the heatmaps ( Figures 2A, C ).

In the validation cohort of 143 AML patients, the expression of SLC25A3, SDHC, RACK1/GNB2L1, FASTKD1, ATP5J, CLIC1, GLUD1, and SLC25A29 were found to be significantly upregulated ( Figure 2B , Table 1 ) while FASLG, HRK, ALAS2, SLC25A21, CYP1B1, SNCA, MMP9, and OLFM4 were significantly downregulated ( Figure 2D , Table 1 ) compared to controls. Two selected genes, LIG1 and MRPL51 did not show significant dysregulation while LONP1 had a reverse expression in the validation compared to transcriptomic expression profile. Upon comparison with TCGA dataset of adult AML patients, similar dysregulation was observed for ALAS2, SLC25A21 and SLC25A29 genes while a reverse expression pattern was observed for ATP5J and CLIC1 genes; none of the other genes showed significant dysregulation in the TCGA dataset ( Table 1 ).

On univariable analysis, increased mtDNA copy number was significantly associated with poor event free survival (HR= 2.14; 95%CI (1.39-3.29); p=0.001) and overall survival (HR= 2.77; 95% CI (1.70-4.59); p<0.001) ( Figures S2A, B ). The timed AUC of mtDNA copy number for predicting 12 months and 18 months survival was 0.66 and 0.68 respectively ( Figures S2C, D ). In patients with increased mtDNA copy number, expression of SLC25A3, SDHC, RACK1/GNB2L1 and FASTKD1, were significantly higher compared to those with low mtDNA copy number ( Figures S2E–H ). Exclusive elevated expression of these 4 genes were also observed in transcriptome of samples with high/intermediate mtDNA copy number(AMLCN_H/AMLCN_I) compared to low mtDNA copy number (AMLCN_L). On correlation analysis, these 4 genes along with 2 other genes CLIC1 and ATP5J showed significant positive correlation with mtDNA copy number ( Table S7 ).

On multivariable analysis, upregulated expression of 2 genes, SDHC (HR 1.29; 95% CI (1.14-1.41); p<0.001) andCLIC1(HR 1.20; 95% CI (1.04-1.38); p=0.013), and downregulation of SLC25A29(HR 0.88; 95% CI (0.83-0.93); p<0.001) were found to be independently predictive of worse OS ( Table 2 ) and they were included for the development of a prognostic gene signature model. All these 3 genes (SDHC, CLIC1, SLC25A29) satisfied internal validation by bootstrapping ( Table S8 ), and were finally selected for prognostic model building. Beta coefficient of each of the variables were used for calculation of risk score as follows:

The formula was used to calculate risk score of all the patients. Risk score median value (10.382) was taken as the cut-off for subgrouping patients into high-risk and low risk group. Patients with high-risk scores (≥10.382) had inferior OS (HR 1.010; 95% CI (1.007-1.014): p<0.001) compared to those with low-risk score (<10.382) ( Figure 3A ). Harrel’s C-index of the prognostic model was 0.675. The timed AUC of the risk score for 12 months and 18 months survival was 0.747 and 0.736 respectively ( Figures 3B, C ).

On multivariable Cox regression analysis, upregulation of SDHC (HR 1.225; 95% CI (1.100-1.363); p<0.001) and downregulation of SLC25A29 (HR 0.905; 95% CI (0.860-0.952); p<0.001) were also predictive of worse EFS. We also found that patients with high-risk score had significantly lower EFS as compared to low-risk score patients (HR 1.008; 95% CI (1.001-1.012); p<0.001) ( Table 2 ). Harrel’s C-index of prognostic model was 0.626. The timed AUC of the risk score for 12 months and 18 months EFS was 0.617 and 0.612 respectively.

On univariable Cox regression analysis of clinical variables, ELN intermediate/poor risk and absence of chloroma were significantly associated with inferior OS and only ELN category came out to be an independent prognostic factor in multivariable analysis ( Table S9 , Figure 3D ). Furthermore, on multivariable analysis, both the ELN risk category (p=0.040) and risk score (p<0.001) were found to be independent prognostic factors for OS. We also performed multivariable analysis including mtDNA copy number and observed that all three factors i.e. risk score (p<0.001), ENA risk categories(p=0.012) and mtDNA copy number(p=0.012) were independent prognostic factors for OS ( Table S10 ).

To compare the predictive ability of our gene signature risk score and ELN risk stratification on OS of AML patients, a time dependent AUC was constructed. Harrel’s C-index of the ELN risk stratification was 0.59 and the timed AUC of ELN risk category on 12 months and 18 months survival was 0.60 and0.64 respectively ( Figures 3E, F ). We combined the ELN risk strategy with our risk score and calculated the predictive ability of the model. The Harrel’s C-index of the model was 0.688 and the timed AUC of combining ELN risk strategy with gene signature risk score for 12 months and 18 months was 0.761 and 0.765 respectively ( Figures 3H, I ).

We found that a high-risk score was significantly associated with poor risk cytogenetics(p=0.021), absence of RUNX1-RUNX1T1 translocation (p=0.027) and ELN intermediate/poor risk group (p=0.016). Furthermore, the proportion of patients achieving CR was significantly higher in the low-risk group as compared to the high-risk group(p=0.017) ( Table 3 ). On subgroup analysis, it was observed that the mitochondria-related gene signature risk score category was significantly predictive of survival outcome across all clinically relevant subgroups except in those with intermediate/poor-risk karyotype ( Figure S3 ).

The predictive ability of gene signature score along with ELN risk stratification on survival outcome of paediatric AML patients was also assessed. Patients with low gene signature score (low risk) belonging to ELN good risk category had significantly better survival outcome (Median OS: Not reached) and predicted 12-months (80% ± 6%), as well as 18-months (75% ± 7%) survival. Similarly, patients with high gene signature score (high risk) belonging to ELN intermediate/poor risk category had significantly inferior outcome (4.67 months (0-3.71)) with 12-months and 18-months predicted survival of 33% ± 8% and 25% ± 7% respectively. On the other hand, patients belonging to other groups (ELN intermediate/poor risk and low-risk; ELN good risk and high-risk score) had intermediate survival outcome (median survival of 27.77-22.90 months respectively) between the two other groups ( Figure 3G ; Table 4 ).

Using our risk calculation model, we calculated the risk score in TCGA dataset (n=179) and similarly, patients were further sub-grouped as high-risk score (higher than median) and low risk score (lower than median) based on the median value (43.434). Kaplan Meier analysis showed that patients with a high-risk score (≥43.434) had inferior OS (HR 1.01;95% CI (1.00-1.02); p<0.019) compared to those with a low-risk score (<43.434) ( Figure 3J ; Table 2 ). Along with this, poor risk karyotype patients had worse overall survival (HR 1.89; 95% CI (2.95-1.20); p=0.004) compared to patients with good risk or intermediate risk karyotype (7.03 vs 18.96months). On multivariable analysis karyotype (poor vs good risk/intermediate risk) and risk score were found to be independently predictive of (p=0.002; p=0.025 respectively) for worse OS. The timed AUC of risk score for 12-months and 18-months survival in the TCGA dataset were 0.64 and 0.63 respectively ( Figures 3K, L ) and Harrel’s C-index of the prognostic model was 0.600. In addition to this, high risk score was also found to be associated with adverse clinical feature of intermediate/poor risk cytogenetics in TCGA dataset as well ( Table S11 ).