Bioinformatic analysis of microRNA expression in Huntington's disease.

Huntington's disease (HD) is an inherited, progressive neurodegenerative disease caused by a CAG expansion in the Huntingtin (HTT) gene and various dysfunctions of biological processes in HD have been proposed. Although monogenic, the exact pathogenesis of HD currently remains unclear. To identify the synergistic microRNA (miRNA) pattern in HD, the miRNA expression profile dataset GSE64977 and the gene expression profile dataset GSE64810 were downloaded. Programming software R was used to identify differentially expressed genes (DEGs) and differentially expressed miRNAs (DEMs). Target genes of DEMs were predicted using the TargetScan database. Gene ontology (GO) function of DEGs was generated using the FunRich and a miRNA‑mRNA interaction network was constructed using Cytoscape software. In total, 1,612 DEGs and 10 DEMs were identified. GO terms mainly included inflammatory response and immune response in DEGs. A total of 745 target genes were predicted from the DEMs and 33 overlaps were identified between these target genes and DEGs. The miRNA network demonstrated that hsa‑miR‑4488, hsa‑miR‑196a‑5p, and hsa‑miR‑549a had a high degree and may be involved with the pathogenesis and potential therapeutic targets of HD.

Introduction

Huntington's disease (HD) is a frequent and incurable hereditary neurodegenerative disorder that impairs motor and cognitive functions (1). With autosomal dominant inheritance, typical mid-life onset, and unrelenting, progressive motor, cognitive and psychiatric symptoms over 15–20 years, the impact of HD on patients and their families is devastating (2). Although caused by a dominantly inherited CAG trinucleotide repeat expansion in the HD gene on chromosome 4 (3), the pathogenesis of HD has not yet been fully elucidated. Previous studies showed the number of CAG repeats was associated with the age of onset; however, only 50 to 70% of the variation can be attributed to repeat size (4). In addition, significant variation in clinical phenotypes is not well explained (5). Consequently, all these variations indicate other pathogenic factors such as heredity affect the disease progression. Recently, studies have focused on miRNAs, the small non-coding RNAs that participate in transcriptional regulation and translational repression of target genes (6). In HD, a neurodegenerative disease caused by a trinucleotide repeat expansion, miRNAs can interact with RelA/NFkB, p53 (7), Mitofusin2 (8), TBP (9), REST, or RE1 (10,11). Dysregulation of miRNA may also impact CAG length and affect the progression or severity of HD (12). Conversely, a certain degree of genetic heterogeneity of HD that may exhibit different miRNA expression in some cases has also been reported (13). In addition to research on the pathological mechanism, the role of miRNA in the treatment has been a subject of interest. Currently, miRNA-based Huntingtin (HTT)-lowering therapy is one of the most advanced gene strategies (14), silencing the HD gene by injecting artificial miRNA into the striatum of Q140/Q140 mice and transgenic sheep models and achieving the expected effect (15,16). However, bioinformatics studies focusing on miRNA and mRNA expression profiles in HD patients and healthy controls have not been published to date. Therefore, identifying the miRNA-mRNA interactions, understanding their synergistic effects on the pathogenesis, and exploring possible therapeutic approaches for HD are important.

To elucidate the miRNAs and associated target genes and pathways involved in HD, we downloaded miRNA and mRNA expression profiles of HD patients and healthy controls from the Gene Expression Omnibus (GEO) database. The differentially expressed genes (DEGs) and differentially expressed miRNAs (DEMs) target genes were identified and a miRNA-mRNA regulatory network established.

Materials and methods

miRNA and mRNA expression profiles

The mRNA and miRNA expression profiles of HD were downloaded from the GEO database (www.ncbi.nlm.nih.gov/geo) and were termed GSE64810 (https://www.ncbi.nlm.nih.gov/gds/?term=GSE64810; accessed August 8, 2017) and GSE64977 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE64977; accessed August 9, 2017), respectively (6,17). In the GSE64810 dataset, 20 HD and 49 neurologically normal control samples from post-mortem human subjects were included. In the GSE64977 dataset, expression profiles were obtained from 28 HD and 36 neurologically normal control prefrontal cortex samples. Both mRNA and miRNA profiling was performed using the Illumina HiSeq 2000 (Homo sapiens) platform (Illumina, Inc., San Diego, CA, USA).

Differential expression analysis for mRNA and miRNA expression profiling

Raw microarray data were first preprocessed including background correction and normalization. The probes corresponding to multiple genes were abandoned and multiple probes corresponding to one gene or miRNA were used to calculate the average expression level. Limma package (www.bioconductor.org/packages/release/bioc/html/limma.html) was used to identify DEGs and DEMs between HD and control samples. A two-tailed Student's t-test was used for statistical analysis. The DEGs and DEMs were considered significantly differentially expressed if the P-value was <0.05 and the log fold change (FC) was >1. The miRNAs were further classified using hierarchical clustering analysis.

Prediction of miRNAs associated with DEGs and construction of miRNA gene regulatory network

Target genes regulated by DEMs were predicted using TargetScan (http://www.targetscan.org/vert_71/), which is used to predict biological targets of miRNAs by searching for the presence of conserved 8mer, 7mer, and 6mer sites that match the seed region of each miRNA (18). In mammals, predictions are ranked based on the predicted efficacy of targeting as calculated using cumulative weighted context ++ scores of the sites (19). In the present study, the miRNA target gene set with context ++ scores <-0.4 were then used for further analysis. A Venn diagram of DEGs overlapping with DEM target genes was constructed using FunRich (20). Furthermore, Cytoscape was utilized to visualize the miRNA gene regulatory network (21).

Functional enrichment analyses

Functional enrichment analyses were performed using the open software FunRich (http://funrich.org/faq), which is a stand-alone software tool used mainly for functional enrichment and interaction network analysis of genes and proteins (20). Statistical cut-off of enrichment analyses in FunRich software was set to default with a P-value <0.05 after Bonferroni correction.

Results

Differential expression analysis

Based on the gene expression analysis, 8 upregulated and 2 downregulated DEMs were identified in HD compared with the normal controls (Table I). Based on the mRNA expression profile analysis, 1,612 DEGs were identified, including 945 upregulated and 667 downregulated DEGs. The top 10 DEGs are presented in Table II. The hierarchical clustering heat map of miRNA shows the differences between HD and normal controls (Fig. 1).

Table I.

Differentially expressed miRNAs between HD patients and healthy controls.

miRNA	Log FC	P-value
hsa-miR-10b-5p	4.175	9.10×10–23
hsa-miR-196a-5p	2.430	5.85×10–22
hsa-miR-615-3p	1.678	2.12×10–17
hsa-miR-10b-3p	1.480	4.79×10–14
hsa-miR-196b-5p	1.439	6.62×10–11
hsa-miR-144-3p	1.045	2.77×10-6
hsa-miR-549a	1.111	2.39×10-5
hsa-miR-483-5p	1.256	2.49×10-4
hsa-miR-10a-5p	1.042	6.0×10-4
hsa-miR-4488	−1.193	2.54×10-3

P<0.05 was considered to indicate a statistically significant difference. miRNA, microRNA; FC, fold-change; HD, Huntington's disease.

Table II.

Top 10 most differentially expressed mRNAs between HD patients and healthy controls.

Gene	Log FC	P-value
PITX1	4.770	9.57×10–39
POU4F2	3.962	3.42×10–23
HAND1	3.703	1.46×10–17
HOXD9	3.657	1.22×10–18
SLC16A12	3.514	4.74×10–18
PITX2	3.404	1.66×10–12
BMP5	3.149	5.93×10–13
OGN	3.097	8.20×10–14
SLC22A2	3.089	6.93×10–11
IL1R2	3.038	3.35×10–12

P<0.05 was considered to indicate a statistically significant difference. FC, fold-change; HD, Huntington's disease.

Figure 1.

Hierarchical clustering heatmap of 10 DEMs in different datasets.

Functional enrichment analysis

Using FunRich, 48 enriched GO terms were obtained for DEGs. The 10 most significantly enriched GO terms are listed in Table III, including integral components of the plasma membrane, inflammatory response, plasma membrane, and immune response. Furthermore, the most significantly enriched GO terms including biological processes (BP), molecular function (MF), and cellular component (CC) were analyzed (Fig. 2).

Table III.

Top 10 GO functional annotation of differentially expressed genes.

GO ID	P-value	Term
GO:0005887	6.55×10–21	Integral component of plasma membrane
GO:0006954	1.54×10–19	Inflammatory response
GO:0005886	4.77×10–19	Plasma membrane
GO:0006955	3.86×10–19	Immune response
GO:0005576	4.51×10–18	Extracellular region
GO:0005615	1.43×10–17	Extracellular space
GO:0045087	1.65×10–15	Innate immune response
GO:0009952	1.61×10–14	Anterior/posterior pattern specification
GO:0043565	1.22×10–11	Sequence-specific DNA binding
GO:0042742	4.73×10–10	Defense response to bacterium

P<0.05 was considered to indicate a statistically significant difference. GO, Gene Ontology.

Figure 2.

The GO and biological pathway analysis of all DEGs. (A) Biological pathway analysis of DEGs, (B) BP of DEGs. (C) CC of DEGs. (D) MF of DEGs.

Synergistic miRNA network construction

Utilizing TargetScan, 745 target genes of DEMs and 33 overlaps were identified between the target genes and DEGs (Fig. 3). We identified 33 possible miRNA-mRNA target pairs. In the network, hsa-miR-4488, hsa-miR-196a-5p, and hsa-miR-549a had a high degree and may be involved in the pathogenesis and used as potential therapeutic targets for HD (Fig. 4).

Figure 3.

Venn diagram of DEGs overlapping with DEM target genes using FunRich.

Figure 4.

Interaction network of DEMs and target genes in HD. Interaction networks for (A) hsa-miR-4488, (B) hsa-miR-196a-5p, (C) hsa-miR-459A, The red and green colors represent upregulation and downregulation, respectively. Interaction network of DEMs and target genes in HD. (D) hsa-miR-144-3p, (E) hsa-miR-10b-3p, (F) hsa-miR-10a-5p and (G) hsa-miR-483-5p. The red and green colors represent upregulation and downregulation, respectively.

Discussion

HD is an inherited progressive neurodegenerative disorder that usually affects people in midlife (22). Previous studies have implicated that gene expression may be altered at more than one stage of RNA processing, translation, protein post-transcriptional trafficking or modifications (23). miRNAs are small non-coding RNAs that can influence diverse ranges of cellular processes (24). Recently, dysregulation of miRNAs has been reportedly associated with neurological and neurodegenerative disorders (25) and several studies have explored the functions of miRNAs in HD patients (10,26). Although altered expression of miRNAs has been reported in cellular models (10) and mouse models of HD (27) using quantified microarray technology, a bioinformatics study of miRNA and mRNA expression profiles in HD patients and healthy controls has not yet been performed. In the present study, we elucidated the synergistic effects of miRNAs on the pathogenesis of HD by constructing a miRNA-mRNA interaction network.

Using a bioinformatics approach to analyze differential expression profiles of mRNA and miRNA, HD patients showed significant differences in mRNA and miRNA expression when compared with neurologically normal controls. In total, 1,612 DEGs and 10 DEMs were identified. GO function (BP, CC, and MF) and biological pathway analysis of DEGs demonstrated the majority of genes were involved in some processes and pathways, such as signal transduction in BP (28), plasma membrane in CC (29), transcription factor activity in MF (30), and immune system in biological pathway (28). Changes in Ca2+ signaling and/or transduction systems and misfolded protein-plasma membrane interactions affected HD initiation and progression as reported by Fan and Shrivastava (29,31). Many previous studies showed that aberrant transcriptional regulations play important roles in the molecular pathogenesis of HD. The HTT gene interferes with some important transcription factors (32); among these, Sp1 and TAFII130 were shown disrupted in the early stage of HD progression (33). Furthermore, both innate and adaptive immune systems have been previously suggested to be activated during the progression of HD (34). Reportedly, HD is characterized by cellular and molecular features of inflammation (cytokine expression and microglia activation) and HTT mRNA expression in immune cells is on average higher than that observed in most organs (35). miRNAs may also influence the dysregulated production of both type 1 and type 2 cytokines observed at different stages of HD. This result is attracting attention as pathogenetic mechanism and as possible therapeutic approach with immunomodulation (36). Combined with the previous research conclusions, the present functional analysis may provide novel therapeutic targets or possible pathogenesis to be further studied.

In the context of the miRNA-mRNA interaction network, hsa-miR-4488, hsa-miR-196a-5p, and hsa-miR-549a had a high degree in this study (Table IV). Hsa-miR-4488 was downregulated and had 12 target genes in our study. Lee et al also found that miR-4488 was differentially downregulated in their HD patients (27). Among the target genes of hsa-miR-4488 in the present study, SPRED3 (37), IGF2 (38), HOXC6 (39), NANOS2 (40), SLC6A3 (41), TLX1 (42), and NT5C1A (43) have been reported to play various roles in the pathogenesis or development of HD. Although we could not determine the association of FAM167B, ECEL1, OTP, C10orf105, or MYLK2 with HD, further experimental studies are necessary to elucidate their possible relationships with HD. Previous studies also showed hsa-miR-196a plays an important role in the pathogenesis and progression of HD (44), Fu et al suggested that hsa-miR-196a altered the RIG-I-like receptor signaling pathway and the immune system, as well as changed the expression of several well-defined pathways of HD, such as apoptosis and cell adhesion (45). A previous study showed that HOX is indirectly involved in the neuroprotective response in HD and increased expression of HOX genes can enhance H3K27me3 or impair PcG repression (39). In our present study, hsa-miR-196a-5p was found to mainly interact with HOX genes, indicating hsa-miR-196a-5p may be involved in the neuroprotective response in HD. Reportedly, miR-196a dominantly altered the immune system or adaptive immune system (45). During the last decade, a hyperreactive immune system has been recognized as an important feature of HD pathogenesis. Macrophages in individuals with manifested HD and even pre-manifested HD have been reported to release more proinflammatory cytokines such as TNF-alpha (46) and prototypical anti-inflammatory cytokines such as IL-10 (47) and IL-13 (48). Similarly, Dobson et al showed that at baseline, monocytes from HD subjects released more cytokines than monocytes isolated from healthy volunteers, and this abnormality could be modulated by laquinimod, which exerts an immunomodulatory effect on isolated HD monocytes (49).

Table IV.

Some of the miRNAs associated with HD.

Author, year	miRNA	Change in expression	Role in HD	Target genes	(Refs.)
Burgunder J-M et al, 2014	miR-144-3p	Upregulated	Not found	GPR183, MAP3K8, MSX1	(4)
	miR-10b-3p	Upregulated	Influence CAG length in HD	HOXD1, ZG16B
	miR-10a-5p	Upregulated	Not found	HOXB3, TFAP2C
	miR-483-5p	Upregulated	Not found	MZB1
Traeger U et al, 2014; Shrivastava AN et al, 2017; van Hagen M et al, 2017; Fan MMY et al, 2007; Moumne L et al, 2013; Dunah AW et al, 2002; Squitieri F et al, 2006	miR-4488	Downregulated	Post-translational protein modification; post-transcriptional regulation by preventing translational initiation;	SPRED3, IGF2, HOXC6, NANOS2, SLC6A3, TLX1, NT5C1A, FAM167B, ECEL1, OTP, C10orf105, MYLK2	(28–34)
Butland SL et al, 2014; Salem L et al, 2016; Hoss AG et al, 2014; Rokavec M et al, 2014	miR-549a	Upregulated	Regulate transcription factors, affect autophagy and endoplasmic reticulum stress pathway	CCDC117, GMNN, LAPTM5, RGS18, TNFAIP8	(37–40)
Cheng P-H et al, 2013; Fu M-H et al, 2015; Zou J et al, 2015	miR-196a	Upregulated	Neuroprotective effect by interacting with HOX genes; affects the entire immune system leading to disorder of cytokine secretion	GPCPD1, HAND1, HOXA5, HOXA7, HOXA9, HOXB6, HOXB7, HOXC8	(44–46)

P<0.05 was considered to indicate a statistically significant difference. miR, microRNA; HD, Huntington's disease.

In addition, we found another miRNA, hsa-miR-549a, was upregulated in HD patients and could regulate five target genes. When reviewing the literature, we found that GMNN (50), TNFAIP18 (51), LAPTM5 (52), and RGS18 (53) played roles in the pathogenesis of neurodegenerative diseases through various mechanisms, including regulation of transcription factors and affecting autophagy and endoplasmic reticulum stress pathways and might be involved in the pathogenesis of HD. Other affected miRNAs found in our study including hsa-miR-144-3p, hsa-miR-10b-3p, hsa-miR-10a-5p, and hsa-miR-483-5p have been reported as differentially expressed in HD samples; hsa-miR-10b-3p showed a significant association with CAG length in HD (6). Regarding the other three miRNAs, more research is be needed to identify their exact relationships with HD.

Although extensive research has identified aberrantly expressed miRNAs in HD, the molecular mechanisms underlying the pathological implications remain largely unknown. Using expression profile datasets, we compared the genomic expression status of HD and revealed differentially expressed mRNAs and DEMs. We identified DEGs and constructed a miRNA-mRNA regulatory network. We found hsa-miR-4488, hsa-miR-196a-5p, and hsa-miR-549a had a high degree and may be involved in the pathogenesis of HD. Studies in this field could help improve the understanding of how miRNAs mediate the etiopathological mechanisms of HD. Since the neuroprotective effects of certain miRNAs have been demonstrated in animal studies, the therapeutic potential of miRNAs should be further investigated and followed by molecular validation.