Single nucleotide polymorphisms in microRNA binding sites on the HOX genes regulate carcinogenesis: An in-silico approach.

Homeobox proteins, encoded by HOX genes, are transcriptional factors playing a crucial role in the master regulatory pathway in the cells. Any mutations in HOX genes will affect the expression of its allied proteins. Such mutations were correlated to the development of different cancer types. In this study, we found 15 HOX genes with a potential target to miRNA, which regulates the translation of the protein by binding to its mRNA through the 3'UTR region. Single nucleotide polymorphisms (SNPs) in this binding region could drastically affect the protein expression by affecting the number and the stability of miRNA-mRNA complexes. We found 77 miRNAs in 15 genes which were found to have altered binding efficiency because of 26 SNPs. After which, we tried to evaluate the impact of each of these SNPs on related HOX genes. Some SNPs such as SNP 15689 on the HOXB7 gene will decrease gene expression by creating or enhancing new binding sites for miRNA to mRNA, while other SNPs such as SNP 872760 on the HOXB5 gene will overexpress the gene by breaking or decreasing existing binding sites from miRNA to mRNA. Then we conducted an expression analysis to compare the mRNA expression profiles in normal and cancer tissue. Subsequently, we did an enrichment analysis followed by a network analysis to shed light on the metabolic function of the gene that could be affected by mutation and whether these mutations may affect other genes. For the first time, this study delivers information on the possible epigenetic regulation of HOX genes via the 77 miRNAs that have predicted target binding sites on HOX mRNAs, and SNPs may regulate those. Furthermore, we show that the HOX gene misregulation may influence other HOX and non-HOX genes, based on network analysis.

Introduction

Homeobox proteins (encoded by homeobox genes) are a crucial component of major regulatory pathways that have been studied intensively for their potential role in cancer development. These proteins are a transcriptional factor that mainly works as regulators of their downstream target genes by either activating or repressing the transcription of these genes [1]. Thirty-nine HOX genes are found on four different chromosomes and are classified into four genomic clusters (A, B, C, and D) [2,3]. Studies on normal and malignant tissues have revealed that HOX gene expression plays a role in the development of different cancers such as colorectal [4,5], lung [6], breast [7,8], and prostate cancer [9,10]. Interestingly, the expression of HOX genes could either be upregulated or be downregulated according to the specific type of cancer [2,11].

Previously it was thought that only the mutations in the intronic region of a gene might affect its expression, but recently, it is known that epigenetic modulations such as miRNA interactions also affect protein expression [12]. miRNA are the small nucleotide sequences of RNA that can bind to the 3′ UTR region of mRNA and regulate the protein expression through RNA Induced Silencing Complex (RISC), leading to mRNA silencing [[13], [14], [15]]. A single nucleotide polymorphism (SNP) mutation in the miRNA binding regions of the 3′ UTR of any mRNA can radically affect its expression in many ways, such as by creating a new binding site for miRNA or by enhancing the stability of miRNA-mRNA complexes to increase binding of RISC leading to the under-expression of that protein [16]. On the other hand, SNP mutations could have a counteracting effect by deleting existing binding sites or destabilizing miRNA-mRNA complexes which could abolish the binding affinity of RISC, leading to overexpression of that protein [17]. It is worth mentioning that the effect of SNPs in miRNA binding sites on mRNAs has been extensively studied on many different genes, such as the CEP gene family [18], OCA and OA gene [19], and colorectal cancer-associated genes [20], however, HOX genes and SNPs in their mRNAs have not been explored.

In this study, for the first time, we investigated the potential effect of SNPs mutation in the miRNA binding site of all the HOX genes and correlated its impact in different cancer types. Using in silico analysis, we filtered SNPs present in the 3′UTR region of HOX genes that can play a role in carcinogenesis, after which we characterized the mutation that could create/delete the binding site or enhance/decrease the stability of miRNA-mRNA complexes on HOX gene mRNAs. Then expression and enrichment analyses of HOX genes were performed to filter out potential tissues and the molecular and biological functions that the mutation could influence. Lastly, network analyses were used to get an idea about the effect of HOX gene SNPs on the other genes.

Materials and methods

HOX gene dataset

Uniprot Knowledge Base (UniProtKB, http://www.uniprot.org/) was used to obtain the list of HOX genes [21].

Identification and analysis of cancer-causing SNPs in miRNA target sites

MirSNP database server (http://bioinfo.bjmu.edu.cn/mirsnp/search/) was used to acquire information about SNPs in the target site of miRNAs [22]. The miR-TS-SNPs attained from the MirSNP database after uploading the gene list obtained from UniPortKB as text format were classified into four groups based on their potential effect. (i) Create— in which the nucleotide base change will introduce a new miRNA binding site in the variant mRNA, (ii) enhance—where the effect of the mutation enhances the binding of the originally targeted miRNA to the variant mRNA (iii) break-in which the substitution will disrupt the miRNA binding site and (iv) decreased through which the derived allele reduces the binding efficacy of the originally targeted miRNA to the variant mRNA.

Obtaining HOX genes expression profile in different tissues and cancer types

Genotype-Tissue Expression (GTEx) online portal (https://www.gtexportal.org/home/) was used to assess HOX genes mRNA expression profile in normal tissues. For analysis of candidate gene expression in cancer tissue, the TCGA online portal (https://portal.gdc.cancer.gov/) was accessed [23]. This analysis helps us understand tissues that are more likely to have different gene expression patterns upon allele change. Moreover, the mRNA profile of HOX genes during cancer types was analyzed.

Functional annotation and enrichment analysis of HOX proteins

The WEB-based Gene Set Analysis Toolkit (WebGestalt, http://www.webgestalt.org/) was used to analyze the enrichment of functions by uploading the gene list in batch using gene symbol under “homo sapiens” using “hypergeometric” as a statistical method [24].

Network analysis of HOX genes with miR-TS-SNPs

For identifying the interaction network of HOX genes, GeneMANIA plugin (http://www.genemania.org/plugin/) was used to upload the HOX gene list Cytoscape [25,26]. This web-based tool uses the GeneMANIA algorithm, which incorporates different association networks into one composite network by a conjugate gradient optimization algorithm [27].

Results

SNPs in miRNA binding sites

Our findings revealed that out of 15 of 39 HOX genes have a potential target for miRNA. In total, 77 miRNAs in 15 genes were found to have altered binding efficiency dues to the effect of 26 SNPs. Some of these SNPs could enhance the binding between mRNA and miRNA complexes or create new binding sites; conversely, other SNPs could delete the binding site or decrease the binding between mRNA and miRNA complexes. This mutation could bring the mRNA of HOX genes to be under the regulation of new miRNA, thereby causing cancer, or they may decrease the number of natural miRNAs that regulate the gene expression leading to overexpression of HOX genes.

Table 1 showed a list of SNPs found to occur in the 3′UTR region of each gene. Also, the SNP id and miRNAs that target the genes, the type of nucleotide change in the allele, and the possible effect of that alteration on miRNA binding to mRNA were included. Many mutations are found in the database of MirSNP for HOX genes, but many of them rarely occurs in humans. Therefore, a population filter has been used in the database to include mutation with MAF more than 0.1 % in the population, and consequently, we obtained the 26 SNPs in 15 HOX genes listed in Table 1.

Table 1

List of SNPs associated with miRNA-mediated post-transcriptional alteration of HOX genes.

HOX genes	miRNA	rsID	Allele Change	Effect
HOXB2	hsa-miR-2277–5p	2,229,302	C = =>T	enhance
	hsa-miR-4706	2,229,302	C = =>T	create
	hsa-miR-4749–5p	2,229,302	C = =>T	create
	hsa-miR-4787–5p	2,229,302	C = =>T	break
	hsa-miR-4795–5p	1,042,815	C = =>T	create
	hsa-miR-542–5p	2,229,302	C = =>T	create
	hsa-miR-5681a	1,042,815	C = =>T	break
	hsa-miR-638	2,229,302	C = =>T	create
HOXB4	hsa-miR-3074–5p	6,504,397	C = =>T	create
	hsa-miR-4446–5p	6,504,397	C = =>T	break
	hsa-miR-4529–3p	6,504,397	C = =>T	decrease
HOXB5	hsa-miR-1304–3p	872,760	C = =>T	break
	hsa-miR-216a	872,760	C = =>T	decrease
	hsa-miR-7-5p	9299	A = =>G	decrease
HOXB7	hsa-miR-3935	15,689	C = =>T	enhance
	hsa-miR-4773	15,689	C = =>T	enhance
	hsa-miR-501–5p	15,689	C = =>T	create
	hsa-miR-5692 b	15,689	C = =>T	create
	hsa-miR-5692c	15,689	C = =>T	create
	hsa-miR-889	15,689	C = =>T	create
HOXB9	hsa-miR-1257	4,239,158	A = =>G	break
	hsa-miR-147a	4,239,158	A = =>G	create
	hsa-miR-1587	2,303,486	A = =>T	enhance
	hsa-miR-181a-2-3p	8844	A = =>G	break
	hsa-miR-3691–5p	4,239,158	A = =>G	create
	hsa-miR-3911	4,239,158	A = =>G	create
	hsa-miR-4252	8844	A = =>G	break
	hsa-miR-4693–5p	4,239,158	A = =>G	decrease
	hsa-miR-4752	4,239,158	A = =>G	create
	hsa-miR-4772–3p	2,303,486	A = =>T	enhance
	hsa-miR-888–3p	8844	A = =>G	break
HOXB13	hsa-miR-3131	3,107,299	A = =>G	create
	hsa-miR-373–3p	2,280,354	A = =>G	decrease
	hsa-miR-379–3p	3,107,299	A = =>G	decrease
	hsa-miR-411–3p	3,107,299	A = =>G	decrease
	hsa-miR-4305	2,280,354	A = =>G	create
	hsa-miR-4751	3,107,299	A = =>G	create
	hsa-miR-4777–5p	2,280,354	A = =>G	break
	hsa-miR-517–5p	2,280,354	A = =>G	decrease
	hsa-miR-520a-3p	2,280,354	A = =>G	enhance
	hsa-miR-611	3,107,299	A = =>G	create
HOXC6	hsa-miR-2467–3p	1,249,077	A = =>G	decrease
	hsa-miR-28–5p	1,249,077	A = =>G	break
	hsa-miR-3139	1,249,077	A = =>G	break
	hsa-miR-3184–5p	1,249,077	A = =>G	create
	hsa-miR-423–5p	1,249,077	A = =>G	create
	hsa-miR-4489	1,249,077	A = =>G	create
	hsa-miR-4633–3p	1,249,077	A = =>G	break
	hsa-miR-4731–5p	1,249,077	A = =>G	decrease
	hsa-miR-486–3p	1,249,077	A = =>G	decrease
	hsa-miR-708–5p	1,249,077	A = =>G	break
HOXC8	hsa-miR-4509	4,142,680	C = =>T	create
	hsa-miR-4639–5p	4,142,680	C = =>T	break
HOXC10	hsa-miR-1226–5p	4,237,810	C = =>T	break
	hsa-miR-3158–3p	4,237,810	C = =>T	break
	hsa-miR-4701–3p	4,237,810	C = =>T	enhance
	hsa-miR-4721	4,237,810	C = =>T	break
HOXC11	hsa-miR-3941	3,816,153	G = =>T	break
	hsa-miR-4731–3p	3,816,153	G = =>T	create
	hsa-miR-4801	3,816,153	G = =>T	create
	hsa-miR-521	3,816,153	G = =>T	break
	hsa-miR-5589–3p	3,816,153	G = =>T	decrease
HOXC13	hsa-miR-3145–5p	2,241,937	A = =>G	enhance
	hsa-miR-3656	4,759,058	A = =>C	create
	hsa-miR-3675–5p	4,759,058	A = =>C	decrease
	hsa-miR-3677–3p	1,710,538	A = =>G	break
	hsa-miR-4705	2,241,937	A = =>G	create
	hsa-miR-503	4,759,058	A = =>C	create
HOXD3	hsa-miR-1205	711,830	C = =>T	enhance
	hsa-miR-4783–5p	711,830	C = =>T	enhance
HOXD4	hsa-miR-1260a	1,063,657	C = =>T	create
	hsa-miR-1260 b	1,063,657	C = =>T	create
	hsa-miR-1303	1,063,657	C = =>T	decrease
HOXD9	hsa-miR-3185	711,822	A = =>G	enhance
HOXD11	hsa-miR-1292	6,745,764	A = =>G	enhance
	hsa-miR-196 b-3p	6,745,764	A = =>G	enhance
	hsa-miR-4475	863,678	A = =>C	break

Many factors would affect the type of impact of SNPs on miRNA binding. Some SNPs may create or delete (break) binding sites for miRNA, while others could enhance or decrease the stability of miRNA-mRNA complexes. For example, SNP 872760 would break or decrease the binding between 2 miRNAs (hsa-miR-1304–3p and hsa-miR-216a) and the mRNA of the HOXB5 gene resulting in overexpression by abolishing control of gene expression. In contrast, SNP 15689 will either create new or enhance existing binding sites for 6 miRNAs (hsa-miR-3935, hsa-miR-4773, hsa-miR-501–5p, hsa-miR-5692 b, hsa-miR-5692c and hsa-miR-889) on HOXB7 gene. This would increase the number and stability of miRNA binding to mRNA, which will consequently increase mRNA degradation, leading to the under-expression of the HOXB7 gene.

It is worth mentioning that one SNP could create new binding sites or enhance the stability for several miRNAs, and at the same time, it might have an opposite impact on other miRNA binding sites by either deleting the binding site or destabilizing miRNA-mRNA complexes. For example, SNP 2229302 on the HOXB2 gene will create a new binding site on four miRNAs (hsa-miR-4706, hsa-miR-4749–5p, hsa-miR-542–5p, and hsa-miR-638) also it will delete (break) existing binding site on another miRNA (hsa-miR-4787–5p). This implies that the numbers of miRNA targeting HOXB2 are more than the number of miRNAs that break the binding with mRNA, which will result in the under-expression of HOXB2.

mRNA expression profile analysis

GTEx analysis was performed to elaborate on potential tissues of the body that could be affected by HOX gene-related miRNA mutations. As seen in Fig. 1, the expression pattern of HOX genes (in normal and cancerous tissues is shown. It is expected that HOX genes would be expressed in almost all tissue of the body, but their expression profile was very low in the central nervous system, cardiovascular system, and whole blood. But a higher expression of most of the HOX genes was seen in the prostate, kidneys, and breast tissue (Fig. 1A). Interestingly, mutations HOXB13 has been reported to be associated with prostate cancer [28,29], so as we notice from Fig. 1 that the expression of HOXB13 is much less in prostate cancer compared to normal prostate, this decrease could be due to mutation in miRNA binding site leading to under-expression of HOXB13 either by creating a new binding site or stabilizing the miRNA-mRNA complexes.

Fig. 1

Expression profile of HOX genes which are predicted to be significantly influenced by miRNAs. (A) HOX mRNA expression in normal human tissues. (B–I) The expression of HOX genes is significantly influenced by miRNAs in various cancer tissues.

Furthermore, using TCGA RNA-Seq data, we present quantifying transcript abundances of the eight HOX genes (HOXB2, HOXB5, HOXB7, HOX13, HOXC6, HOXD3, HOXD4, and HOXD9), which were found significantly dysregulated in most of the cancers (Fig. 1B–I). SNPs present in the miRNA binding sites of HOXB2, HOXC6, and HOXD3 genes support either create or enhance the effect on the expression of these three genes, and these genes were significantly upregulated in pheochromocytomas/gangliomas (PCPG). Most of these HOX genes have several SNPs in their miRNA binding sites, making them genes of potential interest for further epigenetic studies.

Functional enrichment analysis

To shed light on the functioning of HOX genes, a WebGstalt server was used, which gives a comprehensive list of all molecular processes and biological functions that HOX genes play a role in. This analysis provides the reader with insights into the metabolic functions that possibly be affected by the mutation. Table 2 summarizes all functions of HOX genes.

Table 2

Functional annotation of HOX genes (queried in this study).

Molecular function/biological process (with GO category ID)	R	adjP
Anterior/posterior pattern specification ID:GO:0009952	65.65	2.36E-23
Regionalization ID:GO:0003002	44.19	3.47E-21
Pattern specification process ID:GO:0007389	31.61	2.73E-19
Skeletal system development ID:GO:0001501	32.76	1.15E-17
Embryonic skeletal system development ID:GO:0048,706	80.53	8.78E-15
Skeletal system morphogenesis ID:GO:0048,705	49.04	7.08E-13
Organ morphogenesis ID:GO:0009887	14.59	6.84E-12
Embryonic skeletal system morphogenesis ID:GO:0048,704	84.29	1.71E-11
Embryo development ID:GO:0009790	13.0	2.10E-11
Embryonic morphogenesis ID:GO:0048,598	19.70	8.35E-11
Sequence-specific DNA binding ID:GO:0043,565	22.48	4.57E-20
Sequence-specific DNA binding transcription factor activity ID:GO:0003700	15.15	6.13E-18
Nucleic acid binding transcription factor activity ID:GO:0001071	15.12	6.13E-18
DNA binding ID:GO:0003677	6.62	1.17E-12
Nucleic acid binding ID:GO:0003676	4.64	1.97E-10
Heterocyclic compound binding ID:GO:1,901,363	2.84	2.66E-07
Organic cyclic compound binding ID:GO:0097,159	2.81	2.66E-07
Binding ID:GO:0005488	1.28	3.14E-02
Protein binding ID:GO:0005515	0.42	1.0

GO gene ontology, R ratio of enrichment, adjP P value adjusted by multiple testing.

Gene network study

So far, we have shed light on mutations that could affect HOX gene expression and what potential consequences they may have on their function. We must find interaction among themselves and genes other than queried HOX genes that the mutations might influence. For this purpose, the GeneMania webserver was used to find the interaction between 15 HOX genes and other genes. Fig. 2 illustrates the network interaction of HOX genes with various genes; the interaction includes shared protein domain, co-expression, co-localization, and physical interaction. It is noticed that HOX genes interact with each other and with non-queried HOX genes such as HOXA2, HOXA3, HOXA4, HOXA5, HOXC4, HOXC5, and HOXD8, among others, also showed the weaker interaction with non HOX genes.

Fig. 2

The plot of interactions between HOX genes and other genes. This plot was generated using GeneMANIA plugin which is accessible through Cytoscape. The analysis of the network shows different interactions between HOX genes themselves and with other genes. Black spots represent the studies' HOX genes while gray spots represent genes with possible interaction.

Discussion

HOX gene alterations have been associated with the initiation and progression of different cancer types. They can interfere with cell proliferation, migration, invasion, and regulate stemness-related properties [30]. In this study, for the first time, we investigated the SNPs in the miRNA binding site of all the HOX genes and their mRNA expressions in different cancer types. Using in silico analysis, we searched for SNPs present in the miRNA binding sites of HOX genes that could create/delete the binding site or enhance/decrease the stability of miRNA-mRNA complexes. Then HOX gene enrichment analyses were performed to filter out potential biological and molecular functions. Lastly, we created a protein-protein interaction network analysis to understand the effect of HOX gene SNPs on the other genes. Brotto et al. have recently published enrichment analysis of HOX genes-induced positive or negative regulation of the expression of genes involved in cancer pathways [31]. However, we have performed a comprehensive analysis of all molecular processes and biological functions of the HOX genes. Our results provide the reader with insights into all the metabolic functions that could significantly affect mutations in HOX genes miRNA binding sites. Also, our results showed the interaction of HOX genes with other HOX genes or non-HOX genes, suggesting misregulation of one of the HOX genes due to SNPs in its miRNA binding sites could potentially affect other closely interactive genes.

Understanding the genetic regulation of early development was greatly enhanced by the discovery of the homeobox genes. HOX proteins exhibit their function as embryonic development regulators by acting through forming complexes with other transcription factors in addition to cofactor proteins [32]. These transcription factors' versatility results from their unique regulation through the epigenetic regulations (such as promoter methylation) and miRNA-mediated regulation of the HOX genes. However, for the first time, we report the presence of SNPs in the miRNA binding sites of HOX genes that could up or downregulate its expression. We found that SNP 872760 could decrease the binding of hsa-miR-1304–3p and hsa-miR-216a to the mRNA of the HOXB5 gene resulting in overexpression by abolishing control of gene expression. Further, SNP 15689 could enhance the existing binding sites for hsa-miR-3935, hsa-miR-4773, hsa-miR-501–5p, hsa-miR-5692 b, hsa-miR-5692c and hsa-miR-889 on HOXB7 gene. That would increase the degradation of HOXB7 mRNA leading to under-expression of the gene.

Genome-wide epigenetic status of genes and then gene expression could be regulated by miRNAs. Also, SNPs play significant roles in cancer development, treatment modalities, and disease prognosis. Hence, SNPs in the miRNA binding site add another layer of complexity in cancer pathogenesis [33]. Several researchers are now searching for the SNPs located in miRNA genes that are associated with cancer susceptibility. Such as, an SNP at the miR-1231 in the 3′-UTR of IFNRA1 was shown to influence the risk of developing hepatocellular carcinoma [34]. Similarly, we performed an extensive and comprehensive analysis for searching SNPs in miRNA targeting HOX genes. A lot of mutations were found, but many of them were unrelated to humans. Hence, we used the population filter with the mutation that has MAF<0.1 % and obtained the 26 SNPs in 15 different HOX genes. These SNPs could enhance or reduce the binding of miRNA and mRNA, ultimately affecting the HOX gene expression.

The HOX family of genes has a critical regulatory function in mammalian development, and its pathophysiological functions are now being studied intensely. HOX genes play a central role in cell proliferation, tissue maintenance, and renewal via stem cells and are highly expressed in cancer stem cells [35]. Previously it was shown that HOX genes are globally less frequently and less strongly expressed in cancer tissues. However, only 28 different cancer types were utilized in their results, excluding PCPG [36]. In contrast, we performed analysis across all the cancer types showing significant differential expression in various HOX genes. We presume that, probably due to the SNPs present in the miRNA binding sites of HOXB2, HOXC6, and HOXD3, these genes were upregulated in PCPG. Recent evidence showed that transcription factor HOXB9 has an essential role in cancer progression by promoting cellular motility, angiogenesis, and metastasis in different cancer types [37]. Therefore, targeting HOX genes could be a promising approach to modulate chemo-resistance. Also, the expression level of HOX genes could prove to be a potential biomarker for targeted and personalized therapies. However, the identification of HOX targets and their function in different types of cancer is incomplete. In this study, we reported the expression profile of all the HOX genes in normal tissues and pan-cancer analysis of mRNA expression levels for the most significant deregulated HOX genes.

In conclusion, for the first time, we report a possible role of miRNAs in cancer development at the protein level by studying the SNPs on miRNA target sites of HOX genes. Our study suggests that 15 genes were affected by miRNA due to SNPs in their 3′UTR regions. The impact of these SNPs may lead to either over or under expression of HOX genes that will consequently cause cancer. Also, we performed an enrichment analysis that delivers some information about the biological and molecular functions of HOX genes and their prominent role in post-transcriptional nucleic acid and protein binding, which emphasizes the importance of HOX genes in cellular regulation. Furthermore, our network analysis conceives that the mutational effect of any HOX gene would consequently affect other members of HOX genes in addition to other genes. The list of SNPs in predicted HOX genes miRNA sites provide us with the most potential HOX genes and SNPs to target as a direction for future studies.

Author contributions

WEH and JSM envisioned the concept, performed the in-silico analysis, reviewed the literature, and prepared the figures, tables and wrote the manuscript. MAE and JSM provided critical review and finalized the manuscript. All authors agree to the final version.

Funding information

JSM is supported by the seed grant (AJF2018036) from the , and the competitive research grant (1801090139) from the Research Institute of Medical and Health Sciences, .

Declaration of competing interest

The authors declare (1) no financial support for the submitted work from anyone other than their employer; (2) no financial relationships with commercial entities that might have an interest in the submitted work; (3) no spouses, partners, or children with relationships with commercial entities that might have an interest in the submitted work; and (4) no non-financial interests that may be relevant to the submitted work.