Identification of the Thyrotropin-Releasing Hormone (TRH) as a Novel Biomarker in the Prognosis for Acute Myeloid Leukemia

Acute myeloid leukemia (AML) is a biologically and genetically heterogeneous hematological malignance with an unsatisfactory risk stratification system. Recently, through the novel single-cell RNA sequencing technology, we revealed heterogeneous leukemia myeloblasts in RUNX1-RUNX1T1 AML. Thyrotropin-releasing hormone (TRH), as biomarkers of CD34+CD117bri myeloblasts, were found to be prognostic in RUNX1-RUNX1T1 AML. However, the clinical and genetic features of TRH in AML patients are poorly understood. Here, with data from TCGA AML, TRH was found to be downregulated in patients older than 60 years old, with DNMT3A and NPM1 mutations, while overexpressed in patients with KIT mutations. This was further validated in three other cohorts of primary AML including Beat AML (n = 223), GSE6891 (n = 461), and GSE17855 (n = 237). Furthermore, we demonstrated that the expression of TRH in AML could be used to improve the ELN 2017 risk stratification system. In conclusion, our preliminary analysis revealed that TRH, a novel biomarker for AML patients, could be used to evaluate the survival of AML.


Introduction
Acute myeloid leukemia (AML) is a genetically and biologically heterogeneous hematological malignance with impaired hematopoiesis and bone marrow failure [1]. To induce a complete remission (CR), the backbone of intensive chemotherapy is still a combination of cytarabine and anthracyclines, including doxorubicin and daunorubicin, which reached a more than 65% CR rate in AML patients under 65 years old. However, relapse remains the leading cause of treatment failure [2,3]. It is estimated that the median overall survival (OS) of AML is no more than 8.5 months, and 5-year OS is only 24.0% [4].
Clinical risk stratification for AML is complex due to patient-associated factors and disease-related factors [3]. Especially for disease-related factors, the evaluation of leukemic cell genetic changes is still under exploration. Recently, through whole-exome sequencing (WES) and whole-genome sequencing (WGS), the Cancer Genome Atlas (TCGA) identified potential driver mutations and generated the genomic landscape of AML [5]. Moreover, the European Leukemia Net (ELN) guidelines incorporated gene mutations including RUNX1 and TP53 for AML stratification, which underlined the importance of genetic mutations in disease prognosis [6]. Gene expression based on RNA-seq was also demonstrated to be efficient in the risk stratification of AML [7].

Analysis of Gene Expression and Activated Pathways in RUNX1-RUNX1T1 AML
Differentiation expression analysis was performed by using DESeq2 [20]. The expression level of genes was evaluated by fragments per kilobase million (FPKM). Salmon [21] was used to determine and quantify the RUNX1-RUNX1T1 transcript in t (8;21) AML as previously described [8]. Gene set enrichment analysis (GSEA) was performed with the GSEA software using the Kyoto Encyclopedia of Genes and Genomes (KEGG) gene set as described before [8,9].

Immune Infiltration Analysis
CIBERSORT was used to explore the proportion of immune cells in the bone marrow microenvironment as previously described [22]. Specifically, CIBERSORT was performed using the LM22 signature gene file with 100 permutations as the minimum. For each sample, a value was generated and only when p < 0.05 was viewed as significant and further analyzed. Correlation between TRH expression and immune cells inferred by CIBERSORT was determined by Spearman correlation.

Statistical Analysis
The Kaplan-Meier method was used to estimate the probability of overall survival (OS). Log-rank tests were used to compare the P value. Statistical analyses were performed with R software (version 4.0.2, https://www.r-project.org/ (accessed on 20 August 2020)).

Pan-Cancer Analysis Revealed the Unique TRH Expression Pattern in AML
The expression of TRH in pan-cancer is rarely reported. Here, GEPIA2 [23] was utilized to analyze the RNA-seq data of pan-cancers from TCGA and the GTEx projects. We found that the expression level of TRH was actually low in most cancers and matched normal tissues ( Figure 1). For hematological diseases, the expression levels of TRH in patients with t (8;21), inv (16), and t (15;17) were higher than those in patients with acute lymphoblastic leukemia (ALL), chronic lymphocytic leukemia (CLL), chronic myeloid leukemia (CML), and myelodysplastic syndrome (MDS) ( Figure S1a). With single-cell RNA-seq data from Human Cell Landscape database [29], we found that the expression of TRH was rather low in most tissues ( Figure S1b). In addition to the high expression level of TRH in AML, only ovarian serous cystadenocarcinoma (OV) had a significantly differential expression of TRH compared with the normal ovaries (p < 0.05) ( Figure 1). Further survival analysis of TRH expression in the OV showed no prognostic significance, both for disease-free survival (p = 0.07) and overall survival (p = 0.23) ( Figure S2a,b). Thus, the expression and clinical analysis in pan-cancers showed that TRH was a unique biomarker for AML.
Given that AML patients had a higher expression of TRH compared to healthy populations (Figure 1), and the Kaplan-Meier analysis of the TRH-low group indicated an inferior outcome for AML patients [9], we further explored the prognostic value of TRH expression in different AML FAB subgroups. For FAB-M2 and M4 AML patients, the TRH-high group presented a significantly better outcome than the TRH-low group ( Figure S2e,f). For M0 and M1 AML patients, we observed a trend of better prognosis for the TRH-high group, though not reaching a statistically significance (p > 0.05) ( Figure S2c,d). For the other M3, M5, M6, and M7 subtypes, the survival comparison was unable to be performed due to no comparable patients in two subgroups.

Metabolism Pathway Activation with TRH Expression in t (8; 21) AML
To analyze the regulation of TRH expression, the RNA-seq data from RU RUNX1T1 AML patients (n = 62) [9] were analyzed and grouped into TRH-high ( and TRH-low (n = 31), based on the median level of TRH expression. Considering th of the fusion RUNX1-RUNX1T1 in t (8;21) AML, we analyzed and compared the RU RUNX1T1 fusion between these groups (Figure 2a). We found that there was no s cant correlation between the RNA expression of RUNX1-RUNX1T1 fusion transcri that of TRH (R = −0.16, p = 0.22). The Cistrome [26] database was used to predict the latory potential score of the top 20 transcription factors for TRH, which consis JARID2, FOXA1, SPI1, and ESR1 ( Figure 2b). Furthermore, using the RNA-seq d RUNX1-RUNX1T1 AML, we observed that there was a significant correlation be

Metabolism Pathway Activation with TRH Expression in t (8; 21) AML
To analyze the regulation of TRH expression, the RNA-seq data from RUNX1-RUNX1T1 AML patients (n = 62) [9] were analyzed and grouped into TRH-high (n = 31) and TRH-low (n = 31), based on the median level of TRH expression. Considering the role of the fusion RUNX1-RUNX1T1 in t (8;21) AML, we analyzed and compared the RUNX1-RUNX1T1 fusion between these groups ( Figure 2a). We found that there was no significant correlation between the RNA expression of RUNX1-RUNX1T1 fusion transcript and that of TRH (R = −0.16, p = 0.22). The Cistrome [26] database was used to predict the regulatory potential score of the top 20 transcription factors for TRH, which consisted of JARID2, FOXA1, SPI1, and ESR1 ( Figure 2b). Furthermore, using the RNA-seq data of RUNX1-RUNX1T1 AML, we observed that there was a significant correlation between TRH expression and JARID2 (R = 0.34, p = 0.0067), ESR1 (R = −0.26, p = 0.038) ( Figure 2c). However, there was no significant correlation between TRH expression and RUNX1 (R = −0.12, p = 0.37), nor RUNX1T1 (R = −0.2, p = 0.12) (Figure 2c).  Gene set enrichment analysis (GSEA) showed that the TRH-high group had significant activation of the oxidative phosphorylation pathway, citrate cycle, TCA cycle pathway, and the alanine aspartate and glutamate metabolism pathway, thus demonstrating the activation of metabolism-relevant pathways (Figures 2d and S3a). In addition, we observed an activation of spliceosome and RNA polymerase pathways ( Figure S3a). Further correlation analysis showed a significantly positive correlation between TRH and splicing factors SRSF2 and U2AF2 (Figure 2e). For the TRH-low group, GSEA exhibited an enrichment of the chemokine signaling pathway, cytokine-cytokine receptor interaction pathway, and the Fc gamma receptor-mediated phagocytosis pathway, indicating the activation of immune-related pathways ( Figure S3b).

Clinical Characteristics of TRH Expression in the TCGA AML Dataset
We further explored the clinical and genetic features of TRH expression in AML patients. With data from the TCGA AML project (n = 145), we found that AML patients older than 60 years had a lower expression of TRH (p = 0.05) (Figure 3a). There was no difference of TRH expression in either peripheral white blood count (WBC) or gender subgroups (Figures 3a and S4).

Validation of the Clinical Relevance of TRH Expression in Other AML Cohorts
To validate the clinical and genetic relevance of TRH expression in other AML cohorts we screened and obtained RNA-seq data of primary AML from Beat AML cohort (n = 223 [14]. We also observed a lower expression of TRH in AML patients older than 60 years old reaching a statistical significance (p = 0.027) (Figure 4a). AML patients with a higher pe We then analyzed the correlation of TRH expression with common driver mutations of AML. In the TCGA AML cohort, AML patients with KIT mutations and RUNX1 mutations had a significantly higher expression levels of TRH compared with the wild-type groups (p = 0.0023 and p = 0.047, respectively) (Figures 3b and S4). On the other hand, AML patients with the DNMT3A mutations and NPM1 mutations presented significantly lower expression of TRH (p = 0.0035 and p < 0.001, respectively) compared with the corresponding wild-type groups (Figure 3c). However, we did not observe a significant difference of TRH expression in terms of FLT3-ITD mutation, IDH1 mutation, IDH2 mutation, or TP53 mutation in TCGA AML cohort (Figures 3c and S4). We also analyzed the effect of methylation status on TRH expression (Figure 3d) and found that there was a significant negative correlation between methylation and the TRH expression (R = −0.40, p < 0.001).

Validation of the Clinical Relevance of TRH Expression in Other AML Cohorts
To validate the clinical and genetic relevance of TRH expression in other AML cohorts, we screened and obtained RNA-seq data of primary AML from Beat AML cohort (n = 223) [14]. We also observed a lower expression of TRH in AML patients older than 60 years old, reaching a statistical significance (p = 0.027) (Figure 4a). AML patients with a higher peripheral WBC had a significantly lower expression of TRH (p = 0.012) (Figure 4a). The genetic mutation pattern of the TRH expression in the Beat AML cohort was also analyzed. AML patients with KIT mutations had higher expression of TRH (p = 0.0004) (Figure 4b). In addition, patients with DNMT3A mutations and NPM1 mutations had lower expression levels of TRH (p < 0.0001) (Figure 4c). Moreover, patients with the KIT mutation had a higher expression of TRH (p = 0.0004) (Figure 4b). These results in the Beat AML cohort were consistent with the TCGA AML cohort. Nevertheless, we did not observe a correlation between TRH expression and RUNX1 mutation ( Figure S5).

Immune Status and Drug Resistance Correlation with TRH Expression in AML
Considering the activation of immune-related pathways in the TRH-low group, we then analyzed the immune infiltration in the bone marrow of RUNX1-RUNX1T1 AML patients using the CIBERSORT algorithm [22]. Correlation analysis showed that there was a significant negative correlation between TRH expression and monocytes (R = −0.28, p = 0.03) and eosinophils (R = −0.26, p = 0.04) (Figure 5a). In the TCGA AML cohort, the CIBERSORT results of bone marrow of AML bone marrow showed that the TRH expression had a significant correlation with a variety of immune cells, including B cells, monocytes, neutrophils, NK cells and CD4 T cells ( Figure S6). Consistent with the result in RUNX1-RUNX1T1 AML, we observed a significant negative correlation of TRH expression with monocytes (R = −0.18, p = 0.03) ( Figure S6).  In addition, we also obtained gene expression data from GSE6891 (n = 461) [1 and GSE17855 (n = 237) [17][18][19]. With the available information of targeted gene m tions, we found that AML patients with NPM1 mutations had significantly lower ex sion of TRH (p < 0.001) in both GSE6891 and GSE17855 (Figure 4d,e). In addition, patients with FLT3-ITD mutations had significantly lower expression of TRH in GSE (p < 0.001) (Figure 4d). AML patients with the KIT mutation had a higher expressi TRH in GSE17855 (p < 0.001) (Figure 4e). CIBERSORT results of bone marrow of AML bone marrow showed that the TRH expression had a significant correlation with a variety of immune cells, including B cells, monocytes, neutrophils, NK cells and CD4 T cells ( Figure S6). Consistent with the result in RUNX1-RUNX1T1 AML, we observed a significant negative correlation of TRH expression with monocytes (R = −0.18, p = 0.03) ( Figure S6). Drug resistance is one of the leading reasons contributing to the treatment failure of AML patients. To test the effect of TRH expression on drug resistance, data from the GDSC database [28] was analyzed. We found that the TRH-high group had a significantly lower IC50 for AML chemotherapeutics, including cytarabine (p = 0.014) and doxorubicin (p = 0.0018) (Figure 5b), which indicated greater sensitivity to chemotherapy. We also observed that the TRH-high group had a lower IC50 for ATRA (p = 0.015), which was used in acute promyelocytic leukemia patients (Figure 5b).

Improvement of Prognostic Stratification for the ELN Risk System
The European Leukemia Net (ELN) risk stratification system for AML [6] is commonly used. We further compared the TRH expression in different cytogenetic risk groups. In the TCGA AML cohort, we found that the favorable group of the ELN 2017 risk system had a significantly higher level of TRH expression than both the intermediate and adverse groups (Figure 6a). In the GSE6891 dataset, we also found that the good-risk group had a significantly higher expression of TRH (p < 0.001) (Figure 6b). Kaplan-Meier plots showed a significant survival difference between the high-and low-TRH groups in the inter- Drug resistance is one of the leading reasons contributing to the treatment failure of AML patients. To test the effect of TRH expression on drug resistance, data from the GDSC database [28] was analyzed. We found that the TRH-high group had a significantly lower IC50 for AML chemotherapeutics, including cytarabine (p = 0.014) and doxorubicin (p = 0.0018) (Figure 5b), which indicated greater sensitivity to chemotherapy. We also observed that the TRH-high group had a lower IC50 for ATRA (p = 0.015), which was used in acute promyelocytic leukemia patients (Figure 5b).

Discussion
In this work, based on what we discovered previously in RUNX1-RUNX1T1 AML, we further extended the research of TRH expression in AML. We revealed the genetic,

Discussion
In this work, based on what we discovered previously in RUNX1-RUNX1T1 AML, we further extended the research of TRH expression in AML. We revealed the genetic, immunological, and clinical features correlated with TRH expression in the TCGA AML cohort, and further validated these results in three other cohorts. Furthermore, we showed that TRH expression could complement and refine the accuracy of the ELN 2017 risk stratification system for AML.
TRH is known to be a component of the hypothalamic-pituitary-thyroid axis, and is involved in the regulation and release of thyroid-stimulating hormone. Recently, with single-cell RNA-seq, we revealed heterogeneous myeloblast populations and their relevant biomarkers in patients with RUNX1-RUNX1T1 AML [8,9]. TRH, as one of the relevant biomarkers, was proved to be an independent prognostic risk factor in RUNX1-RUNX1T1 AML. In this study, we found that the top activated pathways in the TRH-high group of RUNX1-RUNX1T1 AML were the oxidative phosphorylation pathway, the citrate cycle, the TCA cycle pathway, and the alanine aspartate and glutamate metabolism pathway. This revealed the aberrant activation of metabolism signaling pathways in the TRH-high group. In addition, the expression of TRH in t (8;21) AML was not regulated by the RUNX1-RUNX1T1 fusion. This might be due to genetic or epigenetic changes leading to the regulation of TRH expression. In fact, through the query of the RNA-seq data of pan cancer, TRH was barely expressed in a variety of cancers, which might be attributed to its methylation status [30,31].
With data from the MILE study [15,25], we analyzed and compared TRH expression in a variety of hematological cancers including ALL, CLL, CML, and MDS. Unlike AML patients, those with these hematological diseases presented relatively lower expression levels of TRH, suggesting that TRH was a characteristic biomarker for AML. We further explored the underlying correlation of TRH expression with genetic changes in AML. We observed lower expression levels of TRH when AML patients harbored DNMT3A or NPM1 mutations in four independent cohorts. However, when AML patients carried KIT mutations, they presented higher TRH expression levels. This may be due to the fact that the KIT mutation frequently occurred in patients with RUNX1-RUNX1T1 AML. In addition, TRH was reported to be highly expressed in RUNX1-RUNX1T1 AML patients by other studies [32,33].
Early identification of high-risk AML and administration of intensive chemotherapy and transplantation could greatly reduce the relapse rate. Thus, the European Leukemia Net has proposed a risk stratification system that incorporates cytogenetic mutational status and has been widely used clinically [6]. Recently, a novel AML prognostic score (APS) was proposed and demonstrated that RNA-seq might have additional advantages for clinical assessment [7]. In this study, we showed that different risk groups of the ELN 2017 risk system had differential expression of TRH, which could be used to further stratify AML patients. We also found that patients with higher expression of TRH in AML appeared to be more sensitive to chemotherapy. Future studies are still needed to explore the drugs that were more effective in the TRH-low group.
In conclusion, we explored the clinical and genetic features of a novel biomarker, TRH, in AML patients. Our preliminary analysis suggests that TRH could refine the ELN 2017 risk stratification system.

Supplementary Materials:
The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/biom12101359/s1. Figure S1. Expression of TRH in hematological cancer and normal tissues. (a) Expression of TRH in hematological cancer with data from MILE study via BloodSpot. (b) t-distributed stochastic neighbor embedding (tSNE) plots displaying the expression of TRH from samples of both fetal and adult tissue; Figure S2. Kaplan-Meier plots of TRH expression in ovarian serous cystadenocarcinoma (OV) and subgroup of AML in the TCGA project. Kaplan-Meier plots of disease-free survival (a) and overall survival (b) between high-and low-TRH groups in OV. Subgroup analysis (c-e) of high-and low-TRH groups in the TCGA AML project based on the FAB category; Figure S3. Gene set enrichment analysis (GSEA) of TRH expression in t (8;21) AML patients. (a) Top enriched pathways in the high-TRH (n = 31) group classified according to the median level of TRH expression from the 62 de novo t (8;21) AML patients. (b) Top enriched pathways in the low-TRH (n = 31) group classified according to the median level of TRH expression from the 62 de novo t (8;21) AML patients; Figure S4. Clinical and genetic correlation with TRH expression in TCGA AML datasets. Statistical significance was determined using the Wilcoxon test; Figure S5. Clinical and genetic correlation with TRH expression in the Beat AML datasets. Statistical significance was determined using the Wilcoxon test; Figure S6. Correlation of TRH expression with the proportion of immune cells in TCGA AML datasets. Spearman's correlation analysis and correlation coefficient (R) are shown.