To retrieve all RNA-seq differential gene expression studies on HD patients or models, we searched PubMed and data repositories in a systematic way. The following term was used to search PubMed with the Advanced Search function ((huntingt*[MeSH Terms] OR huntingt*[tiab]) AND (Sequence Analysis, RNA [MeSH Terms] OR High-Throughput RNA Sequencing [MeSH Terms] OR Gene Regulatory Networks/genetics* [MeSH Terms] OR Genomics/methods* [MeSH Terms] OR RNA/genetics* [MeSH Terms] OR Gene Expression Profiling [MeSH Terms] OR Transcriptome* [MeSH Terms] OR Transcriptome/genetics* [MeSH Terms] OR RNA-seq*[Title/Abstract] OR “transcriptome profiling” [tiab] OR “transcriptome analysis” [tiab] OR “transcriptional dysregulation” [tiab] OR “genome-wide expression profiling” [tiab] OR “genome-wide expression analysis” [tiab] OR “differential* express*” [tiab] OR “mRNA-seq” [tiab] OR “transcriptome sequencing” [tiab] OR “transcriptional signatures” [tiab] OR “transcriptional alterations” [tiab] OR “transcriptional changes” [tiab]) NOT (review [pt])) AND ((“2013/01/01” [Date - Publication]: “2020/08/15” [Date-Publication])). For data repositories, search terms were adapted to each database: 1) huntingt*[title] AND “transcriptome gene expression” [Filter] for BioProject (SRA and GEO); 2) huntingt* AND exptype:“RNA-seq of coding RNA” for ArrayExpress; 3) huntingt* AND RNA-seq, with filter “study” for European Nucleotide Archive; 4) huntingt* for European Genome-Phenome Archive. In total, 288 articles were identified through PubMed and 237 datasets were identified through data repositories. Then, the articles and datasets were screened with the following exclusion criteria in order of priority: 1) Non-HD study, 2) Erratum/Superseries, 3) Review, 4) Different analysis of same data, 5) Non-RNA study, 6) RNA isoform/noncoding RNA study, 7) HTT knockout/-down study, 8) Treatment study, 9) Microarray, 10) Non-RNA-seq study, 11) Non-peer reviewed and 12) DEG not available. The screening process is presented in a PRISMA flow diagram (Supplementary Figure S1).

For human HD datasets for which all filtered expressed genes were not provided (n = 12), lowly expressed genes were removed from the expression matrix based on the selection criteria defined in each paper. For each subcategory, background gene lists were created by taking the union of all expressed genes in the datasets within the corresponding category.

In total, 14 HD RNA-seq studies were performed on samples from either human cell culture models, blood samples or post-mortem brain samples. Cell culture models were derived from either embryonic stem cells or induced pluripotent stem cells, and samples were collected at the pluripotent, multipotent, or terminally differentiated stage. Two HD blood studies were included, one with monocytes from Grade 0–2 HD patients, one with platelets from all stages of HD. Post-mortem brain samples come from asymptomatic and symptomatic HD patients, specifically, caudate/putamen, motor cortex, prefrontal cortex and cingulate cortex. We analyzed 44 datasets after excluding 6 astrocyte datasets that didn’t analyze HD vs. WT differential expression. Then, based on sample origin, tissue type and disease stage, human HD datasets were grouped into 11 categories-cell culture datasets included pluripotent, progenitor, differentiated_neuronal, and differentiated_non_neuronal; post-mortem brain datasets included brain_early_striatal, brain_early_cortical, brain_late_striatal_neuronal_MSN, brain_late_striatal_neuronal_IN, brain_late_striatal_non_neuronal and brain_late_cortical; and the blood.

For non-human analyses, 14 HD RNA-seq studies from mouse, monkey and sheep models were included. HD mouse models included R6/1, R6/2, HdhQ150/175, BACHD, CHL2 and Swiss Webster mouse; mouse cell lines included E14 STHdhQ111 and BV2 microglia; other HD models included Rhesus-macaque pluripotent cell lines and OVT73 sheep. Then, we graded mouse HD datasets as early, intermediate and late stages based on prior literature. The datasets have been graded as follows. R6/1: 1 and 2 months-early, 7.5 months-late. R6/2: 1 month-early, 2–3.75 months-late. HdhQ150: 8 months-intermediate, 22 months-late. HdhQ175: 2 months-early, 6 months-intermediate, 10 months-late. CHL2: 12 months-intermediate. BACHD: 1 month-early.

The representative differentially expressed genes (DEGs) for all mouse and human categories were defined based on how often a gene was observed as significantly changed across the datasets in each category. Genes that are differentially expressed in at least two datasets were included in the category-representative DEGs. For subcategories with single dataset, all the genes were considered. As an exception, for the human blood category, union of DEGs from the two datasets within the blood category were considered as blood category-representative DEGs as the number of overlapping genes between the two blood datasets is 1.

Gene overlap analysis between pairwise comparisons were performed by using one-sided Fisher’s exact test. To test gene overlap among more than two categories, we used the SuperExactTest (Wang et al., 2015) package in R. Background gene list for each category was defined as the union of all expressed genes in each dataset. All analyses were performed in R.

Gene set enrichment analysis (GSEA) (Subramanian et al., 2005), for differentially expressed genes between HD and control samples, was performed for each HD human dataset (n = 44) with the GSEA method in the WebGestaltR package v.0.4.2 (Liao et al., 2019) for the following functional databases; KEGG pathway and Gene Ontology Biological Process. Differentially expressed genes were ranked based on their log fold changes with HD over wild type.

Common or unique gene sets were analyzed by online Gene Ontology enRIchment anaLysis and visuaLizAtion tool (GOrilla) (Eden et al., 2009). For pairwise comparisons, the intersection between the two datasets’ backgrounds was used as the background set. For multiple comparisons in the subcategories, the subcategory background was used as the background set. For unique genes in a subcategory or dataset, the subcategory or dataset background was used as the background set. Statistical significance was set at FDR <0.05.

To identify potential regulators for DEGs in each HD human dataset, we used hTFtarget database in which human TFs were curated from Chip-Seq experiments in different cell lines, tissues and cells (Zhang et al., 2020b). For each HD human dataset, DEGs were ranked based on absolute logFC and then top 500 DEGs were queried for potential regulators. For the datasets that have less than 500 DEGs, all DEGs were considered for the query. After obtaining potential regulators for each dataset, further filtering was performed, in which regulators that were not found in all expressed gene list of the corresponding dataset were removed. For the datasets that do not have all expressed gene lists, tissue specific gene lists were curated from the (Uhlen et al., 2015) and used for the filtering. An activity score was calculated for each regulator as a ratio of the number of DEGs a regulator controls to the total number of DEGs queried.

The overlap analysis between the TFs obtained from H3K27ac motif data from the HD iPSC Consortium 2017 dataset and the regulators obtained from hTFtarget database for the same study was performed by using Fisher’s exact test. All expressed genes from the same dataset were used as a background.

Human TFs (n = 2,765) and epigenetic modifiers (n = 167) were curated from The Human Transcription Factors website (Lambert et al., 2018) and dbEM database (Singh Nanda et al., 2016), respectively. Among 2,765 curated proteins classified as TF with a known motif, TF with an inferred motif, likely TF, ssDNA/RNA binding protein, or unlikely TF, we only considered the following classes; TF with a known motif, TF with an inferred motif, which revealed 1,639 known and human TFs. Then, differential expression status of the curated human TFs and epigenetic modifiers were checked in all 44 HD human datasets.

This review included 44 differentially expressed gene (DEG) datasets from 14 human HD studies (Agus et al., 2019; Al-Dalahmah et al., 2020; Consortium, 2017; Denis et al., 2019; Lee et al., 2020; Lim et al., 2017; Lin et al., 2016; Mehta et al., 2018; Miller et al., 2016; Ooi et al., 2019; Osipovitch et al., 2019; Ring et al., 2015; Smith-Geater et al., 2020; Switonska et al., 2018) (Table 1). The hierarchical clustering of all 44 datasets showed clustering of datasets irrespective of the differences in the origin of tissue/sample (Supplementary Figure S2) (Ooi et al., 2019; Lee et al., 2020). Thus, we based our analyses on grouping of the datasets according to tissue/cell type and disease stage.

We classified 44 human DEG datasets into five categories based on sample source and disease stage (Figure 2). The five groups were further categorized into the following 11 subcategories: pluripotent, progenitor, differentiated neuronal and differentiated non-neuronal cells from cellular models; brain_early_striatal, brain_early_cortical, brain_late_cortical, brain_late_striatal_neuronal_MSN, brain_late_striatal_neuronal_interneuron and brain_late_striatal_non-neuronal from brain tissue and blood. The top upregulated and downregulated genes in each subcategory are shown in Table 2. There were no commonly differentially expressed genes (DEGs) across all the studies. We then identified frequently dysregulated genes in each subcategory by taking the genes that appeared in at least two of the datasets. Figure 3A shows a heatmap of the frequent DEGs across all the subcategories and shows that DEGs from early cortical tissues changed extensively in the late cortical tissues from symptomatic brains (Figure 3A, Supplementary Figure S3). The late striatal subcategories showed similar differential expression of many genes. Additionally, we observed some similarities among pluripotent, progenitor, differentiated_neuronal and differentiated_non_neuronal cellular subcategories.

Next, we performed Gene Set Enrichment Analysis (GSEA) (Subramanian et al., 2005) on 44 human DEG datasets using the ranked lists of DEGs based on their log fold changes, Gene ontology biological processes and KEGG pathways databases using WebGestalt (Liao et al., 2019).

GSEA analysis showed strong negative enrichment of genes associated with p53 signaling pathway and cancer in pluripotent cells in Switonska_hiPSC_PSC_109Q (Switonska et al., 2018). Contrastingly, Ring_hiPSC_PSC did not show any enrichment of p53 signaling pathways and cancer in pluripotent cells. It showed exclusively negative enrichment of gene sets related to metal ion homeostasis, immune response, regulation of transcription regulatory region DNA binding, JAK-STAT cascade and peptidyl-tyrosine phosphorylation (Ring et al., 2015). Regarding the progenitor subcategory, the neural stem cell dataset from Smith_hiPSC_MSN_NSC had the largest number of enriched gene gets, with positive enrichment in gene sets such as ribosome, spliceosome, DNA replication, base excision repair, RNA transport, mismatch repair, p53 signaling pathway, and cell cycle, and with negative enrichment in gene sets such as synaptic vesicle cycle, oxidative phosphorylation, Alzheimer disease, and Parkinson disease (Smith-Geater et al., 2020). Smith_hiPSC_MSN_adult showed strong positive enrichment of numerous gene sets such as metal ion homeostasis, lysosome, JAK-STAT cascade, regulation of transcription regulatory region, cAMP signaling pathway, nerve development, synaptic transmission, Interestingly, Mehta_hiPSC_CoNeuron datasets showed strong negative enrichment of cancer related genes on day-0 in differentiated neuronal cell culture that changed on day-80 and 100 showing positive enrichment and no enrichment on day-130 (Mehta et al., 2018) (Supplementary Dataset S7). Most differentiated_neuronal datasets showed positive enrichment of metal ion homeostasis and detoxification.

In brain, presymptomatic striatum and cortex showed positive regulation of gene sets associated with intrinsic apoptotic signaling pathway in response to endoplasmic reticulum stress, immune response, transcriptional misregulation in cancer, regulation of protein maturation, heat shock response. In symptomatic brain subcategory, the late cortex was positively enriched for various gene sets related to metal ion homeostasis, protein folding, p53 signaling, cancer and positive regulation of peptidyl-tyrosine phosphorylation. Late cortex and late_striatal_neuronal_MSNs showed common positive enrichment of glutathione metabolism, RNA polymerase, biosynthesis of unsaturated fatty acids, alanine, aspartate and glutamate metabolism, protein localization to chromosome/centromeric region, DNA replication, acidic amino acid transport, postsynapse assembly, purine metabolism, and synaptic vesicle cycle, while sharing negative enrichment in gene sets related to amyloid precursor protein metabolic process, cAMP signaling pathway, nerve development, alcoholism, Huntington disease, synaptic transmission, morphine addiction, long-term depression, Alzheimer disease and circadian entrainment (Figure 4, Supplementary Figure S5; Supplementary Dataset S7). In striatal interneurons, Sst/Npy-expressing (SN) GABAergic interneurons in Lee_CPu_IN_SN dataset and the Pvalb/Th-expressing (PT) GABAergic interneurons in Lee_CPu_IN_PT dataset revealed a contrasting enrichment profile of various gene sets (Figure 4, Supplementary Figure S5; Supplementary Dataset S7).

For the blood monocyte dataset- Miller_blood_mono, there was positive enrichment in gene sets such as cytokine-cytokine receptor interaction, chemokine signaling pathway, IL-17 signaling pathway, JAK-STAT signaling pathway and pathways in cancer (Figure 4, Supplementary Figure S5). Concordant enrichment in gene sets, such as, chemokine signaling pathway and IL-17 signaling pathway was found in early cortical dataset from Agus et al., 2019, as well as in most of the late cortical datasets and late striatal non-neuronal datasets from Lee et al., 2020. Positive enrichment in the gene set JAK-STAT signaling pathway was also found in the neural stem cell and adult-onset HD datasets from Smith-Geater et al., 2020, and most of the late cortical datasets and late striatal non-neuronal datasets (Smith-Geater et al., 2020). For the blood platelet dataset, the only enriched gene set was related to detoxification in the positive direction, and the same positive enrichment was found in neural stem cell and adult-onset HD datasets from Smith-Geater et al., 2020, the Day 60, Day 80 and Day 100 iPSC-derived cortical neuron datasets from Mehta et al., 2018, all the late cortical datasets and both of the interneuron datasets. Meanwhile, negative enrichment in the gene set detoxification was observed in the hiPSC pluripotent stem cell dataset from Ring et al., 2015, the neural dataset from HD iPSC Consortium 2017 as well as in the astrocyte, the microglia, the OPC, the oligodendrocyte, the endothelial cell and the ependymal cell datasets from Lee et al., 2020 (Ring et al., 2015; Consortium, 2017; Mehta et al., 2018; Lee et al., 2020; Smith-Geater et al., 2020) (Supplementary Dataset S7).

To investigate potential regulators responsible for differential expression of genes in HD, we searched for the potential regulators of the DEGs in each dataset using human transcription factors target (hTFtarget) database (Zhang et al., 2020a). An activity score based on the number of DEGs co-regulated by a regulator in a dataset were assigned for each regulator in the respective dataset.

In the circular heatmap depicting the regulators and their activity scores (Figure 5, Supplementary Dataset S8), one third of the heatmap (indicated as group A) depicts the most common set of regulators across the categories and the datasets. In this set of most common regulators, we observed regulators such as, SP4, NRF1, BRD4, RYBP, ZFP64, HDAC1, FOXO1, CREB1, NFKB1, BRD2, SP1, CTCF and JUN some of which were previously implicated in HD pathology (De Souza et al., 2016; Dunah et al., 2002; Steffan et al., 2000). Gene ontology analysis of these transcription factors (TFs) showed significant enrichment of genes related to transcriptional misregulation in cancer, pathways in cancer, Huntington’s disease, longevity regulating pathway and cell cycle (Supplementary Dataset S5). Regulators grouped in B show strong scores for all except for differentiated_non_neuronal, Mehta_hiPSC_CoNeurons, pluripotent cell and presymptomatic brain. Some of the less common regulators (Figure 5, group C), such as PCGF1, KLF4, ELF3 and EGR2 showed strong scores in all categories except the datasets in the late striatal category from Lee et al., 2020. Interestingly, the DEGs from the blood platelet dataset from Denis et al. were enriched for a group of high-scoring TFs (Figure 5, group G), namely POLR2A, TAF3, BHLHE40, REST, TEAD4, RBBP5, RING1, GTF3C5 and CBFB. Of note, to test the reliability of our method, we checked the overlap between transcription factors (TFs) from the available H3K27ac motif data from a previous publication (Consortium, 2017) and the potential regulators predicted by hTFtarget for the same study. We found a significant overlap of 39 regulators (p-value = 3.79e-10; Fisher’s exact test; Supplementary Dataset S8).

Next, we investigated the differential expression of TFs and epigenetic modifiers in HD using the human TFs (Lambert et al., 2018) and dbEM (Singh Nanda et al., 2016) databases, respectively, in each of the 44 datasets. Of 1,639 human TFs and 167 human epigenetic modifiers, we found differential expression of 1,388 TFs and 155 epigenetic modifiers in HD in at least one of the 44 datasets (Supplementary Dataset S9). Figure 6 shows the top common up-regulated and down-regulated TFs and epigenetic modifiers in HD in all 44 human datasets. The top common up-regulated TFs were TCF4, ZBTB20, ZHX2 and BCL6, while the most common downregulated TFs were SON, ZNF302, ZKSCAN1 and ZMAT1. The most common up-regulated epigenetic modifiers were IN O 80D, CREBBP, HDAC4 and KDM4B, while the most common downregulated epigenetic modifiers were PRMT2, PRMT8, JMJD1C and TAF1. Interestingly, in MSNs, three of the top upregulated epigenetic modifiers, KDM4B, KMT2C and SMARCD1 were uniquely downregulated, while the top downregulated epigenetic modifier genes SETD5 and PRMT8 were upregulated (Figure 6).

Finally, we investigated whether the predicted regulators of DEGs were themselves differentially expressed in HD in the 44 human DEG datasets. In total, 204 out of the predicted 557 regulators were differentially expressed in at least one of the datasets (Supplementary Dataset S9). Among the differentially expressed regulators, TCF4, FOSL2, BCL6, ZBTB16, FOXP1, KLF15, RXRA, CUX1, CREBBP and NFIA were the top upregulated regulators that were also predicted to be the regulators of DEGs in HDare, while the top downregulated ones were ZKSCAN1, MAX, E2F3, BCL11A, EGR1, FOXG1, TP53, HMG20A, JMJD1C and STAT4 (Supplementary Dataset S9).

In order to evaluate the suitability of non-human models to study transcriptional changes in HD patients, we compared the non-human data to human mRNA changes. For this, we retrieved 14 studies and 79 datasets (Table 5; Supplementary Figure S1). Among 14 studies, 12 were related to mouse models (Ng et al., 2013; Vashishtha et al., 2013; Crotti et al., 2014; Mielcarek et al., 2014; Achour et al., 2015; Langfelder et al., 2016; Miller et al., 2016; Handley et al., 2017; Jacobsen et al., 2017; Pan et al., 2018; Goodnight et al., 2019; Radulescu et al., 2019; Yildirim et al., 2019; Lee et al., 2020; Wertz et al., 2020), one study was related to sheep (Handley et al., 2017) and one to monkey (Goodnight et al., 2019).

Similar to human data, in the mouse data, we checked the expression changes of 45 previously published key HD genes (Seredenina and Luthi-Carter, 2012; Vashishtha et al., 2013). Among the 45 key HD genes, the commonly dysregulated genes in mice were Ryr1, Scn4b, Gpr6, Hrh3, Foxp1, Gpr88, Rgs9, Pde10a, Arpp21, Adora 2a, Adcy5, Drd2, Pcp4 (Supplementary Figure S6A). Moreover, the available monkey dataset showed genes such as CAMKV, KCNIP2, OPRD1, RGS9, HRH3, ADORA2A, KLF16, CNR1, OPRK1, FOXP1, NR4A1, TESC, GSN and RARB that were dysregulated in monkey model of HD (Supplementary Figure S6B). These findings provide evidence that the HD animal models faithfully recapitulate human HD regarding its molecular pathology.

While different HD mouse models show disease progression at different ages, we sorted all the datasets along a disease progression timeline based on the literature regarding HD mouse models (see methods). For this, three studies were excluded as they were cell culture models and could not be confidently assigned to a disease stage. Then, we categorized the mouse datasets as early, intermediate and late stages (Figure 7A), in which the latter two were comparable to the symptomatic stage in humans. Out of 1,037, 13788 and 12609 total DEGs (Supplementary Dataset S10), 58, 5201 and 6985 genes were frequently dysregulated in early, intermediate and late disease stage of HD mice, respectively, of which only 47 genes were common across all the three disease stages (Figure 7B). Some of the common genes were Scn4b, Penk, Arpp21, Kcnab1, Ngef, Ppp4r4, Gpx6, Slc4a11, Rgs4, Rgs14, Adora2a, Gpr153, Plk5 and Wnt4 (Supplementary Dataset S10). In early stage, genes such as Clspn, Gpr153, Gpx6, Krt9, Penk and Plk5 were frequently dysregulated, while in intermediate stage, Penk, Atp2b1, Scn4b, Pde10a and Rbfox1 were dysregulated to name a few. In late stage, Pcp4, Scn4b, Pde10a, Penk Atp2b1, Hspa8 and Nrxn3 were dysregulated among many others (Supplementary Dataset S10). The top upregulated and downregulated genes in each dataset are shown in Table 6. Additionally, we compared single cell transcriptomic analysis of medium spiny neurons and interneurons with bulk striatal data. We observed that of 4,936 differentially expressed genes in MSNs, 2,522 genes were unique in medium spiny neurons and out of 173 in interneurons, 77 genes were unique to interneurons (Figures 7C,D).

To evaluate how well the mouse models recapitulate the transcriptional changes in human HD, we compared early disease stages in mice with presymptomatic human brain, and intermediate and late stages in mice with symptomatic human brain. In early/presymptomatic stages, we observed no overlapping DEGs in the striatum and only one overlapping DEG, Serpina1 in the cortex (Figures 7E,F, Supplementary Dataset S10). The overlap between the mouse and human data was comparably large with 1,616, 85 and 192 overlapping DEGs between intermediate-late mouse MSNs, interneurons and cortex with corresponding symptomatic HD human brain subcategories (Figures 7G–I). Some of the genes common between human and mouse are Mtf2, Mtus2, Slc25A, Smarcal1, Snap25, Smarca2, Mt-co3, Camk4 (Supplementary Dataset S10). Moreover, the expression changes of the three common heat shock genes from the human subcategories, Dnajb1, Hspa1b and Hspb1, were variable in the mouse datasets. While the expression of Hspa1b and Hspb1 weren’t altered in majority of the striatal mouse datasets, Dnajb1 showed upregulation in some datasets and downregulation in some others in HD mice Table 4).

We also compared the mRNA data included in our analysis with two key microarray studies in the field, Hodges et al. (2006) and Kuhn et al. (2007). As reported in Hodges et al. (2006), we observed the majority of dysregulated genes in caudate or striatum, which they showed in their study to be independent of neuronal loss in this area. Out of the listed top 30 upregulated and down regulated genes, 15 and 20 genes, respectively were common in our data (Supplementary Dataset S6). Dysregulation of genes associated mainly with heat shock proteins, chaperonins, protein folding and response to unfolded protein were in consistency with the microarray study in human HD caudate. In accordance with Hodges et al. (2006), we also report that the dysregulated genes in MSNs are mostly associated with neurotransmitter receptor activity, as well as with those conveying signals from excitatory amino acids like glutamate receptor signalling pathway. Furthermore, Kuhn et al. (2007) demonstrated that the different genetic models of mouse could replicate the transcriptomic changes of human HD caudate, however, regional specificity in most mouse models was considered to be very low. The study showed stronger similarity in downregulated genes between mouse models and human HD. We observed in our study 43 genes out of 60 top mouse genes that were concordant with human caudate in the microarray analysis of Kuhn et al., 2007. Some of the common genes between our study and Kuhn et al. (2007) are RASGRP2, RGS14, HOMER1, PENK, PPP3CA, PLCB1, PPP1R1A, CAMK2A and RGS4 (Supplementary Dataset S6). Overall, our finding of genes dysregulated in human and non-human HD models using RNA-seq were consistent with the findings from previous key microarray studies.