Whole-transcriptome bioinformatics revealed HTRA3, KRT8, KRT17, and RHEX as novel targets in acute myeloid leukaemia.

Acute myeloid leukaemia (AML) is characterised by heterogeneous genomic signatures that vary among different patient groups. Hence, the current study aims to conduct a whole transcriptome analysis of a female patient with AML and a family history of the disease at the time of diagnosis. Genetic profiling has a useful impact on clinical management and treatment success of the disease as the complex genetic landscape of AML and differential responses to treatment might indicate inadequate therapeutic targeting. A 37 year old female patient with AML was admitted to the hospital complaining of general fatigue arthralgia and chest pain. AML diagnosis was confirmed by complete blood count and blood smears before being confirmed by cytogenetic analysis. Herein, we conducted whole-transcriptome sequencing analysis to assess differential gene expression profiles in patients and healthy control subjects. In addition, single nucleotide polymorphism/insertion-deletion analyses (SNP/INDEL) were performed to investigate gene variants in the present case. The results revealed a remarkable differential gene expression profile in AML compared to the corresponding control at the time of diagnosis, indicating that HTRA3, KRT8, KRT17, and RHEX are potential novel therapeutic targets. Additionally, a number of novel gene variants were also reported in the current study, as concluded from the SNP/INDEL analysis, which might be associated with disease risk assessment and probably affect prognosis. These genes and their new variants might be worth reporting to the scientific community for further exploration of AML.

Introduction

Acute myeloid leukaemia (AML) is a haematological malignancy characterised by sophisticated, heterogeneous genomic signatures, making its clinical management quite challenging (reviewed in). In this respect, genetic abnormalities are associated with disease prognosis and decisions on targeted therapies for patients. Hence, the exploration of whole–genome profiles may advocate the discovery of novel genomic variations that could be linked to disease progression and overall clinical management. In line with this, whole transcriptome ribonucleic acid (RNA) sequencing was recently employed to diagnose AML and detect molecular signatures of the disease in a group of patients with AML. Whole transcriptome sequencing results, in addition to providing information on the expression levels of featured genomic variants of AML, can also highlight novel variants on an individual basis. In concordance with this, it was reported that the expression levels of normal genes could influence disease prognosis and, consequently, clinical decisions. Likewise, another report showed that a single-cell mRNA sequencing platform helped categorise the phenotype of patients with AML according to the disease's genetic profile. Of note, patient karyotype and genetic alterations are known to be pivotal in patient stratification and risk assessment of AML. For example, patients with trisomy 8 are classified as those with intermediate risk and are susceptible to mutations in myeloid progenitor transcription factors, deoxyribonucleic acid (DNA) epigenetic modification genes and others, and RUNX1. In the context of speech about trisomy 8, it was reported that this karyotype abnormality is associated with a lower remission rate and more genetic changes in patients with AML, which leads to worse prognosis. Therefore, trisomy 8 can be considered as a driving force of genetic heterogeneity in patients with AML and an important contributor to leukaemogenesis, leading to overall poor prognosis. Also, among the well-known genetic anomalies in AML is the RUNX1/RUNX1T1 fusion, a variant that manipulates the expression of several genes that drive AML progression. Usually, RUNX1/RUNX1T1 coexists with mutations in other genes, including both KRAS and NRAS and ASXL1 and ASXL2, and its use to monitor minimal residual disease is of clinical importance. Bearing this information in mind, along with the poor prognosis of AML in some individuals, applying whole genome sequencing might help to understand the genetic landscape of AML on a personalised basis. Therefore, the aim of the current study was to analyse the transcriptomic profile of a patient with AML to explore differential gene expression and gene variants that might be correlated with the pathological phenotype of the disease, along with standard diagnostic markers.

Case presentation

In the current study, a 37 year old female patient with a family history of AML was admitted to King Fahd Hospital of the Imam Abdulrahman Bin Faisal University. The patient had a nephew who was diagnosed with AML at the age of 6 years. The patient complained of general fatigue, arthralgia, and chest pain.

Clinical examination showed a well-developed, well-nourished, and alert female. Lymphadenopathy or organomegaly was not observed. She had scattered petechiae, but no purpura or ecchymosis. A complete blood count (CBC) workout with differential revealed a state of thrombocytopenia, which manifested as petechiae along with three consecutive blast counts in blood smears of 40%, 30%, and 41%. Routine haematology investigations and cytogenetic reports have been revisited in our previous publication. In brief, cytogenetics revealed 47, XX, +8[18]/46, XX[1] karyotype and 88% trisomy 8 with RUNX1 partner transcriptional corepressor 1 (RUNX1T1) abnormality, both of which are characteristics of AML. A post-treatment bone marrow aspirate on 25/9/2018 showed a marked suppression of trilineage haematopoiesis with a relative increase in lymphocytes and macrophages. Differential counts showed erythroid precursors (7%), blast cells (2%), segments (2%), lymphocytes (86%), and plasma cells (2%). CBC performed on the same day showed marked leukopenia, neutropenia, and thrombocytopenia. The patient was then recommended allogeneic stem cell bone marrow transplantation, after which bone marrow aspiration and biopsy were again investigated on 19/2/2019. The analysis indicated intergranular haemorrhage with crushing and tissue loss. An average cellularity of 10% was also reported in trilineage haematopoiesis, and no blast cells were recorded. The patient died in 4/2019. Patient consent was documented prior to participation.

Materials and Methods

At the time of diagnosis, RNA was extracted from the blood of the patient and a healthy female control for whole-transcriptome RNA sequencing., Briefly, samples were withdrawn into PAXgene® blood RNA tubes (#762165, preanalytix, Hombrechtikon, Switzerland), and RNA was extracted according to the protocol described in the PAXgene Blood RNA Kit (#762164, preanalytix, Hombrechtikon, Switzerland). Sample quality was first checked after thawing on ice, centrifugation, and mixing for analysis by Agilent 2100 bioanalyzer to assess concentration, RNA integrity number, and estimation of 28s/18s ribosomal RNA ratio. Total RNA fragments were analysed using an Agilent RNA 6000 nano kit (5067–1511, Agilent Technologies, Waldbronn, Germany). Following quality checks, a library of cDNA was constructed for the sequencing platform (BGISEQ-500 platform) as follows: 1) mRNA enrichment: DNA oligo dT nucleotides were used as they are complementary to sequences in rRNA, then these hybrids were degraded by RNaseH to deplete rRNA; 2) RNA fragmentation and reverse transcription: This step encompasses fragmentation of RNA and reverse transcription to double-strand cDNA by N6 random primers; 3) End repair, addition of poly A tailing and adaptor ligation: The synthesized cDNA was subjected to end-repair and then 3′ adenylated. Adaptors were ligated to the 3’ adenylated ends of cDNA; 4) polymerase chain reaction (PCR) amplification: The ligation products were amplified by PCR to enrich the purified cDNA; 5) Denaturation and cyclisation: PCR products were denatured by heat and the single strand DNA was cyclised by splint oligo and DNA ligase; and 6) sequencing on the BGISEQ-500 platform (supplementary Figure S1).

An average of 23.28 million reads per sample were generated with 95.25% mapping of the reference genome. Only clean reads (without unknown bases, without sequencing adapters, and bases with quality scores <15) were included in the analysis. Bioinformatics was performed to understand the differential gene expression profiles by performing gene expression quantification. Essentially, clean reads were mapped to the reference genome using hierarchical indexing for spliced alignment of transcripts (HISAT). Gene expression levels were calculated using RNA-Seq by Expectation-Maximization (RSEM), which estimates the relative abundance of the transcripts. The expression of a transcript is expressed as fragments per kilobase of transcript per million mapped reads (FPKM), which is proportional to the number of cDNA fragments that originate from it, helping to identify the up- and down-regulated genes in differentially expressed genes (DEGs). NOIseq algorithms, were then employed to screen for differentially expressed genes, where genes with log2fold change values ≥ 1 and ≥80% probability were considered upregulated, while those with log2fold change ≤ −1 and a probability ≥80% were considered downregulated in comparison to the corresponding control. This was followed by single nucleotide polymorphism/insertion deletion (SNP/INDEL) analysis, which was done on RNA sequencing data using the genome analysis toolkit. HtrA serine peptidase 3 (HTRA3; forward primer: CGTGGTGGAGAAGATCGCAC; reverse primer: CATGATGAAGCCAGAACCGC), Keratin 8 (KRT8; forward primer: TCCTCGAGCTTCTCCCGAGT; reverse primer: CTGGTTGACCGTAACTGCGG), Keratin 17 (KRT17; forward primer: CAGGAGATTGCCACCTACCG; reverse primer: CCACAATGGTACGCACCTGA), and regulator of haemoglobinisation and erythroid cell expansion (RHEX; forward primer: AAAATGACTCCCCGCTGGAC; reverse primer: GCTGGCTCAGACAGAGCAGG) expression was confirmed by quantitative real time PCR (qRT-PCR). Primers were designed and prepared using BGI Genomics (Shenzhen, Guangdong, China). RNA was reverse transcribed, and qRT-PCR was performed in a 16 μL final reaction volume containing 1 μL template cDNA, 8 μL 2x SYBR Green PCR master mix, 0.2 μL of forward and reverse primers (50 pM/μL) and 6.6 μl RNase/DNase–free water. The ABI ViiA 7 real-time PCR cycler was used with the following program: (95 °C for 2 min) and 40 cycles of (94 °C for 10s, 59 °C for 10s, 72 °C for 40s). Ct values were collected using a corrected threshold value. The 2−ΔΔCt method was used for relative quantification. β-actin was used as the reference gene.

Results

Bioinformatics showed overexpression of many genes at the time of diagnosis compared to the corresponding control profile, suggesting their potential use as candidate biomarkers for the disease (the total number of differentially expressed genes (DEGs) was 2539, as shown in Figure 1).

Figure 1

Library construction and stepwise procedure of RNA sequencing.

Genes with a probability of expression ≥80% and a log2fold change ≥5 were considered for screening. Of those, HTRA3, KRT8, KRT17, and RHEX were among the top 15% of the upregulated genes and met the criteria of DEG screening for probability and log2fold change values, which were further validated by qRT-PCR (Table 1 and Figure 2) and are discussed here because of their reported carcinogenic roles and lack of knowledge on their contribution to leukaemogenesis and AML. Additionally, gene ontology and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed and are summarised in the supplementary file with the presentation of gene ontology terms of the candidate genes reported herein (Figure S1).

Table 1

RNA sequencing and qRT-PCR fold change values of candidate genes in the patient with AML.

Gene name	Log2 fold change	Fold change	Probability	qRT-PCR fold change
HTRA3	5.66	50.5	83.6%	1099.2
KRT8	5.33	40.22	96.4%	18.5
KRT17	5.93	60.9	89.8%	29.67
RHEX	5.04	32.9	98%	54.55

Figure 2

The number of differentially expressed genes between the control and patient with AML. The number of up- and down-regulated genes is depicted as a pheatmap (A) showing the log10 transformed gene expression level; the darker the colour, the higher the expression. The MA plot (B) shows log2 transformed mean expression level on the x-axis and log2 transformed fold change on the y-axis; the red dots represent the up-regulated DEGs, the blue dots represent the downregulated DEGs, and the grey dots represent the non-DEGs.

In addition, RNA sequencing analysis confirmed the diagnosis by ascertaining the overexpression of AML genomic signatures in comparison to the corresponding control (Table 2), including core-binding factor subunit beta (CBFB), CCAAT/enhancer-binding protein alpha (CEBPA), fms-like tyrosine kinase 3 (FLT3), nucleophosmin 1 (NPM1), and runt-related transcription factor 1 (RUNX1). This finding suggests the potential use of RNA sequencing as a diagnostic technique.

Table 2

Differential gene expression of well-known genetic signatures in AML.

Gene	Log2 fold change	Actual fold change	Probability
RUNX1	3.48	11.15	96%
CBFB	1.85	3.6	85.17%
FLT3	5.59	48.16	98%
NPM1	3.14	8.81	95.88%
CEBPA	3.65	12.55	0.97%

Furthermore, SNP/INDEL analyses that are summarised in Table 3 helped identify novel variants of these genes that could be further validated as risk/benefit assessment factors in AML. In summary, the position of the variant and gene sequence (exon/intron) were presented, and genomic variants were verified using the SNP database (dbSNP). Novel variants were indicated by (N) and were particularly reported in RHEX as intron and exon INDEL and single-nucleotide exon variants, in addition to one novel single-nucleotide exon variant in HTRA3.

Table 3

SNP/INDEL analysis of candidate target genes.

INDEL analysis
Gene name	Reference sequence	AML case sequence	Location	Insertion/deletion	position	Chromosome
RHEX	C	CA	Intron (N)	Insertion	206241994	1
	GAC	G	Intron(N)	Deletion	206247131
	G	GAAA	Intron(R)	Insertion	206251589
	C	CA	Intron(N)	Insertion	206265713
	CT	C	Intron(N)	Deletion	206269767
	C	CT	Intron(N)	Insertion	206276043
	A	AT	Exon(N)	Insertion	206285883
	T	TAA	Intron (R)	Insertion	206289128
	A	ATTT	Intron(N)	Insertion	206295869
SNP analysis
RHEX	A	T	Exon (N)	206283483		1
	C	G		206283527
	A	C		206283551
	C	T		206283776
	C	G		206284047
	C	T		206287944
	G	A		206288251
	G	A		206288452
KRT 8	T	C	Exon (R)	53294381		12
	A	G		53298769
HTRA3	A	G	Exon (R)	8297770		4
	C	T	Exon (N)	8297883

N = novel and R = report.

Discussion

Acute myeloid leukaemia is characterised by a complex genomic signature, which contributes to its complicated stratification and the pressing need to categorise patients at the time of diagnosis for successful clinical management. This requires costly diagnostic platforms to identify the molecular profiles of every patient. Therefore, there is a continuous need to search for unique biomarkers and novel target genes in AML that can be verified as potential therapeutic candidates. In this context, employing whole-transcriptome sequencing platforms at different disease stages in different patient groups could yield significant results. In the current report, whole-transcriptome sequencing confirmed the overexpression of well-known genetic signatures of AML, suggesting the potential of this technology to be used as a single diagnostic tool for the disease. In addition, bioinformatics analysis indicated overexpression of other genes at the time of diagnosis in comparison to a normal corresponding control; of these, HTRA3, KRT8, KRT17, and RHEX were validated by qRT-PCR (Table 1), which confirmed their overexpression in the patient as compared to the corresponding control. Furthermore, SNP and INDEL analyses revealed many novel genetic variants in the HTRA3 and RHEX genes that could be verified as risk/benefit factors in patients with AML (Table 3).

HTRA3, which is upregulated in patients with AML at the time of diagnosis, was recently reported as a potential biomarker for gastric cancer, and its expression was found to be correlated with cancer stage and histopathological grade. In contrast, HTRA3 has been recognised previously for its anti-tumour potential as a mediator of apoptosis in lung cancer as well as pancreatic cancer; this indicates the possible multifaceted nature of this gene in different cancers. For this reason, it is important to highlight its overexpression in our study on AML, which might attract researchers’ attention to explore its role in this disease.

With respect to KRT8, its expression was linked to gastric cancer progression and poor prognosis and promoted the invasion capacity of tumours. In a similar context, RNA-sequencing-based data analysis confirmed the overexpression of KRT8 in lung adenocarcinoma and lung squamous cell carcinoma patients and found that its overexpression was associated with poor overall survival and recurrence-free survival. On the other hand, an anti-metastatic role for KRT8 was recently reported in lung cancer cells, wherein cisplatin treatment upregulated KRT8 in cancer-associated fibroblasts, an action which inhibited cell invasion and attenuated metastatic potential. To the best of our knowledge, having indicated the involvement of KRT8 in human cancers, its role and overexpression in AML has not been reported. Tracing KRT8 in patients with AML might then be introduced as a novel therapeutic target in this blood cancer. KRT17 was also upregulated in patients with AML, which is in line with a similar report, implicating this gene in chemoresistance and metastasis in cervical cancer. Likewise, KRT17 gene transcripts and protein products were correlated with clinicopathological characteristics of colorectal cancer patients, as it was recently reported as a prognostic marker to predict disease recurrence in stage II patients. In concordance with this, silencing KRT17 with siRNA produced anticancer effects in gastric cancer cell types. Referring to AML, a statistical study indicated KRT17 downregulation as a marker of the disease, which contradicts our findings. To this end, we have to note the scarcity of literature on the relationship between KRT8 and KRT17 and AML.

The regulator of erythroid progenitor proliferation, RHEX, was recently identified as a marker of lymph node metastasis in early-stage endometrial cancer. Lack of reports on the potential role of RHEX in human cancers, including AML and other types of blood cancers, makes the findings of the current study interesting and indicates that the role of RHEX overexpression in AML is worth further exploration in different disease models and a larger population of patients with AML.

Conclusions

In summary, transcriptomic analysis could represent an important platform for stratifying patients on an individual basis by helping to identify novel genomic signatures that might be targeted during the course of treatment and result in an overall improvement in clinical outcomes. Herein, the technique highlighted HTRA3, KRT8, KRT17, and RHEX as potential novel targets in AML for further validation to assess how monitoring or targeting these gene products could serve the clinical management of AML.

Having cast a beam of light on the potential role of the four candidate genes in AML, the exact role of these genes should be further explored in a larger sample and validated using in-vitro and in-vivo models. Similarly, the novel exonic and intronic variants presented here should be validated to determine their significance in AML.

Source of funding

This work was supported by the Deanship of Scientific Research at the (project number: 2017-CAMS-099).

Conflict of interest

The authors have no conflict of interest to declare.

Ethical approval

Informed consent was obtained from all participants involved in the study, and the study was reviewed and approved by the Institutional Review Board of Imam Abdulrahman Bin Faisal University (IRB# 2017-03-147).

Authors contributions

O S. E.: Study design and conceptualisation, data collection, and drafting of the manuscript; A. A., F. A., and S. A.: Writing and approval of the final draft; K. A: Bioinformatics analysis and approval of the final draft. All authors have critically reviewed and approved the final draft and are responsible for the content and similarity index of the manuscript.