Identification of RHOBTB2 aberration as an independent prognostic indicator in acute myeloid leukemia

Transcriptional levels of RHOBTB genes in patients with AML

The RHOBTB genes comprise three members, RHOBTB1, RHOBTB2, and RHOBTB3, in mammalian cells. We compared the transcriptional levels of RHOBTB genes in the bone marrow of AML patients (TCGA-LAML, n = 173) with those in normal samples (GTEx, n = 70) through Gene Expression Profiling Interactive Analysis 2 (GEPIA2) (Figure 1). Gene expression analysis using box plots indicated that the transcriptional levels of RHOBTB1 and RHOBTB3 were decreased (P < 0.05, Figure 1A and 1C), while that of RHOBTB2 was significantly increased (P < 0.05, Figure 1B). Consistently, three datasets indicated increased RHOBTB2 expression; six datasets and five datasets showed reduced RHOBTB1 and RHOBTB3 expression, respectively, in leukemia compared to normal samples in the ONCOMINE database (Supplementary Figure 1). Downregulation of RHOBTB1 and RHOBTB3 has been reported in various types of tumors, and this pattern was confirmed in AML. RHOBTB2 showed notably divergent expression patterns between AML and other tumor types, which led us to further explore the underlying clinical significance.

Prognostic values of RHOBTB genes in patients with AML

Intriguingly, a high expression level of RHOBTB2 was associated with poor overall survival (OS) (Hazard ratio (HR) (high) = 2.9; P = 0.00041, Figure 2B) while a low expression level of RHOBTB3 was associated with poor OS for patients with AML (HR (high) = 0.44; P = 0.0045, Figure 2C) in GEPIA2. However, there was no difference for RHOBTB1 (HR (high) = 1.1; P = 0.85, Figure 2A). The prognostic evaluation capacity of RHOBTB2 and RHOBTB3 was validated with the UALCAN database (Supplementary Figure 2) and Kaplan-Meier (KM) survival analysis (log-rank P-value = 0.000239 for RHOBTB2, log-rank P-value = 0.00024 for RHOBTB3) in the TCGA-LAML cohort (n = 151) obtained from https://portal.gdc.cancer.gov/ in January 2020 (Figure 3A–3B). RHOBTB2 overexpression was recognized as a risk factor for OS with HR = 1.672 (95% confidence interval (CI), 1.285–2.176) (Figure 3A) via Cox proportional hazards analysis, while high expression of RHOBTB3 as a protective factor with HR = 0.444 (95% CI, 0.288–0.685) (Figure 3B). The median survival time of the RHOBTB2 high-expression group was 0.8 years and that of the low-expression group was 2.3 years. In comparison, the median survival times of the RHOBTB3 high and low-expression groups were 2.3 years and 0.7 years, respectively (Figure 3A–3B). Time-dependent receiver operating characteristic (ROC) analysis of RHOBTB2 and RHOBTB3 was performed to compare each gene's predictive accuracy. RHOBTB2 had a larger area under the curve (AUC) than RHOBTB3, especially for 3- and 5-year survival (3-year AUC = 0.732 vs. 0.673, 5-year AUC = 0.802 vs. 0.720) (Figure 3C–3D). Therefore, the results above suggest that both RHOBTB2 and RHOBTB3 may be potential prognostic factors for patients with leukemia, and RHOBTB2 showed better prognostic performance.

Identification of RHOBTB2 as an independent prognostic indicator in AML

We performed univariate and multivariate Cox regression analyses to determine whether RHOBTB2 and RHOBTB3 are robust AML OS-related genes that can be used for prognosis prediction.

Multiple clinical factors, such as age, WBC count, blast cell percentage, and cytogenetic abnormalities, impact the prognosis of AML. Some individual genes, including four-and-a-half LIM domain 1(FHL1), HOPX and FAM124B have been recently identified as candidate prognostic factors through a genome-wide Cox regression screening project [17]. Thus, RHOBTB2 and RHOBTB3, combined with risk factors including age (≥60 years old), WBC count (≥30 × 10⁹/L), blast cell percentage, cytogenetic abnormalities, therapeutic agent target (FLT3, DNMT3A, and TP53 mutations, etc.), and de novo prognostic indicators (FHL1, HOPX, and FAM124B) (Supplementary Table 1) were used for multivariate Cox regression analysis of the TCGA-LAML cohort (n = 151). As shown in Figure 3E, the forest plots indicated that high RHOBTB2 expression, but not high RHOBTB3 expression, was strongly predictive of poor outcome in AML patients (HR = 1.581; 95% CI, 1.102–2.270; P = 0.012), independent of clinical features including age, WBC count, blast cell percentage and gene mutation status (Figure 3E). Compared to RHOBTB3 and the three potential prognostic indicators (FHL1, HOPX and FAM124B), only RHOBTB2 displayed prognostic value (P-value = 0.003, vs. P-value = 0.057 for RHOBTB3, P-value = 0.896 for FHL1, P-value = 0.087 for FAM124B, and P-value = 0.184 for HOPX), and a higher HR (HR = 1.685; 95% CI, 1.188–2.388) in the Cox model (Figure 3F, right panel), even though all five genes were statistically significant in univariate Cox regression analysis (Figure 3F, left panel). The results above verify that RHOBTB2 can be used as an independent and effective predictor of the OS for AML patients.

Association of RHOBTB2 expression with AML classification and clinical characteristics

The French–American–British (FAB) classification, devised in the 1970s and 1980s, recognizes eight subtypes of AML (FAB M0-M7) based mostly on morphology and cytochemistry [18]. The more precise World Health Organization (WHO) classification has been used in recent years, which considers clinical features, morphology, immunophenotyping, cytogenetics, and molecular genetics [19]. A set of recurring chromosomal and genetic lesions related to oncogenes, tumor suppressor genes, and other regulatory elements that control vital cell functions, such as t(15;17) (q22;q12), PML-RARα, t(8;21) (q22;q22), RUNX1-RUNX1T1, FLT3-ITD mutation and NPM1 mutation, have been introduced in many patients [19–22].

We further clarified the expression profile of RHOBTB2 in different subtypes of AML using UALCAN. The box plots showed that RHOBTB2 was remarkably increased in FAB subtypes M0, M1, M2, M4, and M5 compared to that in M3 (also known as acute promyelocytic leukemia (APL)), which has a higher degree of differentiation (P < 0.001 for M0 vs. M3, M1 vs. M3, M2 vs. M3, M4 vs. M3 and M5 vs. M3, Figure 4A). FAB subtypes M6 and M7 were not considered because of the small sample size for each. It seems like progenitor cells have higher RHOBTB2 expression levels than their progeny.

Next, the correlation between RHOBTB2 expression and clinical features, such as age, gender, FLT3 mutation, PML/RAR-fusion, and RAS activation status in AML patients, was analyzed through UALCAN (Figure 4B–4F). Box plots showed that the expression level of RHOBTB2 was higher in the 61–80-year-old group than in the 21–40-year-old group (P = 0.015, Figure 4B) and higher in the male group than in the female group (P < 0.001, Figure 4C). We did not analyze the 80–100-year-old group, as it had only five samples. In addition, RHOBTB2 expression was increased in AML patients without FLT3 mutation (P = 0.024, Figure 4D) and patients without PML-RAR fusion (P = 0.0038, Figure 4E). The presence of FLT3 mutations in AML enabled the recent approval of targeted drugs that can help patients achieve prolonged remission. PML-RAR is a fusion gene that is associated with the specific subtype of leukemia APL.

Furthermore, the correlation between RHOBTB2 and disease-related genes was analyzed. The statistical scatter plots from the GEPIA2 database showed that RHOBTB2 expression had a positive association with the expression of FLT3 (Pearson’s correlation = 0.37, P = 5.56E-07), FHL1 (Pearson’s correlation = 0.37, P = 3.8E-07), RUNX1 (Pearson’s correlation = 0.17, P = 0.023), RUNX2 (Pearson’s correlation = 0.3, P = 4.9E-05), and RUNX3 (Pearson’s correlation = 0.42, P = 6.2E-09) (Figure 4G). FLT3 is overexpressed on CD34⁺ blast cells in approximately 93% of AML cases and might share the same expression pattern with RHOBTB2 in bone marrow samples. High expression of FHL1 and RUNX1-3 is related to poor prognosis. NPM1, the high expression of which predicts a superior prognosis, was negatively correlated with RHOBTB2, but the correlation was not significant.

Functional enrichment analyses of RHOBTB2 and co-expressed genes in AML

To further explore the potential function and molecular pathways of the RHOBTB2 gene in AML, we utilized the LinkedOmics database [23] to identify co-expressed genes of RHOBTB2 in data of 173 patients from TCGA. A total of 8,198 genes related to RHOBTB2 were altered, which reflects the considerable impact of the core gene RHOBTB2 on AML pathogenesis. The 2,309 gene clusters of these related genes that were positively related to RHOBTB2 are displayed as red dots, whereas the 2,028 gene clusters that were negatively associated with RHOBTB2 are represented by green dots in the volcano plot (P < 0.01 and FDR < 0.01, Figure 5A). The top 20 significant gene sets positively and negatively associated with RHOBTB2 are presented in Table 1.

RHOBTB2 and its top 200 associated gene clusters were subjected to gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analyses to identify enriched categories and signaling pathways in the TCGA-LAML cohort. The bubble diagrams showed that the gene clusters were located in cell morphogenesis involved in differentiation, small GTPase mediated signal transduction, and positive regulation of catalytic activity (http://amigo.geneontology.org/amigo, GO:0000902, GO:0000904, GO:0007264, GO:0043085, etc) in the aspect of biological processes (Figure 5B). The functions of these categories above may be related to cell division, cell cycle, proliferation, and differentiation of blast cells in AML. In the aspect of cellular components, Chang et al. found that RhoBTB2 is distributed in a vesicular pattern [24], and coincidentally these genes are putative structural constituents of phagocytic vesicle (GO:0045335), endocytic vesicle (GO:0030139) and cytoplasmic vesicle membrane (GO:0030659) (Figure 5C). When it comes to the aspect of molecular functions, they are localized in RHOBTB2-related utilities such as Rho GTPase binding (GO:0017048), Rac GTPase binding (GO:0048365) and GTPase activator activity (GO:0005096) (Figure 5D). The top 10 pathways related to RHOBTB2 and co-expressed gene clusters were defined by KEGG analysis. These pathways such as cytotoxicity (https://www.kegg.jp/kegg/, hsa04650), phagosome (hsa04145), leukocyte transendothelial migration (hsa04670), and apoptosis (hsa04210), have been intriguingly implicated in vital pathological processes including cell fate determination and migration of leukocytes (Figure 5E). Given our results and the architecture, localization and biological functions of RHOBTB2 postulated in previous research, we hypothesized that overexpression of RHOBTB2 in leukemic blast cells could regulate cell differentiation, cell cycle, proliferation, apoptosis, migration and vesicle transport.

GEPIA2

GEPIA2 (http://gepia2.cancer-pku.cn/, Beijing, China) is an interactive web-based tool for analyzing cancer-related RNA sequencing data provided by TCGA and GTEx projects [25]. General gene expression profiles, survival analysis, and correlation analysis were conducted through the “Expression Analysis” module with the TCGA- LAML cohort (n = 173) and normal tissues (n = 70), the data of which are available in the panel “dataset sources”. Student’s t-test was used to perform expression analysis. The survival results were displayed by Kaplan-Meier curves with HRs and P values from a log-rank test. A P-value = 0.05 was used as the threshold of statistical significance.

Univariate and multivariate Cox regression analyses

Univariate and multivariate Cox regression analyses were performed to identify candidate prognostic genes in the LAML cohort from TCGA. The data (available through https://portal.gdc.cancer.gov/) were updated in Jan 2020, and the cohort contains 151 AML patients with high-throughput sequencing (RNA-Seq) data and detailed clinical information [26].

A forest plot with the P-value, HR and 95% CI of each variable was built through “survival”, “survminer” and “forestplot” R packages in RStudio 4.0.3. Gene expression levels were dichotomized based on the median expression level in the cohort as the cutoff value.

LinkedOmics

LinkedOmics (http://www.linkedomics.org) provides a unique portal to analyze cancer multi-omics data and clinical data for 32 cancer types and 11 158 patients from TCGA project [23]. Genes associated with RHOBTB2 were identified in the TCGA-LAML cohort (n = 173) and are presented in volcano plots. Pearson’s correlation test was used to evaluate the statistical relationship.

UALCAN

UALCAN (Birmingham, AL, USA, http://ualcan.path.uab.edu) serves as a platform for validating specific genes and screening tumor candidate biomarkers [27]. RHOBTB gene expression in AML subgroups based on various clinicopathologic features and survival outcomes was investigated via “Expression Analysis” and “Survival Analysis” modules, respectively. The processed RNA-sequencing data and survival profiles of the AML cohort (n = 163) were obtained using TCGA assembler (http://www.compgenome.org/TCGA-Assembler/). A P-value of 0.05 was used as the threshold for significance.

GO and KEGG pathway enrichment analyses

The WebGestalt database (http://www.webgestalt.org/option.php) for deriving biological insights from gene lists was exploited to perform GO and KEGG pathway enrichment analyses for RHOBTB2 and the top 200 co-expressed genes. The built-in reference human protein-coding genome was selected as the background parameter. Bubble plots with (-log10) P-values, FDRs, and enrichment ratios were generated through “ggplot” and “dplyr” R packages. A (-log10) P-value > 1.3 was considered to indicate enrichment of a meaningful pathway.