Integrative Meta-Analysis of Huntington’s Disease Transcriptome Landscape

Huntington’s disease (HD) is a neurodegenerative disorder with autosomal dominant inheritance caused by glutamine expansion in the Huntingtin gene (HTT). Striatal projection neurons (SPNs) in HD are more vulnerable to cell death. The executive striatal population is directly connected with the Brodmann Area (BA9), which is mainly involved in motor functions. Analyzing the disease samples from BA9 from the SRA database provides insights related to neuron degeneration, which helps to identify a promising therapeutic strategy. Most gene expression studies examine the changes in expression and associated biological functions. In this study, we elucidate the relationship between variants and their effect on gene/downstream transcript expression. We computed gene and transcript abundance and identified variants from RNA-seq data using various pipelines. We predicted the effect of genome-wide association studies (GWAS)/novel variants on regulatory functions. We found that many variants affect the histone acetylation pattern in HD, thereby perturbing the transcription factor networks. Interestingly, some variants affect miRNA binding as well as their downstream gene expression. Tissue-specific network analysis showed that mitochondrial, neuroinflammation, vasculature, and angiogenesis-related genes are disrupted in HD. From this integrative omics analysis, we propose that abnormal neuroinflammation acts as a two-edged sword that indirectly affects the vasculature and associated energy metabolism. Rehabilitation of blood-brain barrier functionality and energy metabolism may secure the neuron from cell death.


Introduction
Huntington's disease (HD) is a neurodegenerative disorder with a genetic inherent caused by glutamine expansion in the Huntingtin gene (HTT) on chromosome 4 [1][2][3][4]. Usually, the normal population has an average of 17-20 CAG (Glutamine) repeats in the HTT gene, but 40 or more CAG repeats are present in HD. HD patients show typical symptoms such as motor impairment, cognitive defects, chorea, and behavioral changes. Tetrabenazine and deep brain stimulation in pallidum reduce motor-related symptoms. However, in some HD patients, the cognitive deficit can affect memory association and leads to dementia [5]. Most of the current treatments are prescribed based on the symptoms. Selective vulnerability (in the striatum) is the predominant factor in the initial stages of the disease. The striatum is the nucleus of the basal ganglia and is involved in voluntary movement control, cognition, and decision-making. The striatum has major connectivity arbors to the frontal cortex, and in HD, prominent cell loss occurs in the cerebral cortex and striatum, affecting the cortico-striatal connection [6]. The cause of vulnerable neuronal death is unknown [7]. Expression studies can be performed on the striatum directly, but the highly degenerated tissue makes the interpretation difficult. Brodmann Area 9 (BA9) is the Dorsolateral Prefrontal Cortex (DLPFC), primarily connected with the executive subregion of the striatum and involved in motor planning and organization, regulation of intellectual function, and action. Analyzing the gene and transcript expression related to this area may aid in understanding the neuronal vulnerability in HD.
HD gene expression studies showed that heat shock proteins (HSP1A, HSP1B, and HSPB1) are upregulated in HD and havoc in the transcription factor networks affects the downstream gene expression [8]. Human iPSC cell studies demonstrated that downregulation in glutamate and GABA signaling, as well as destruction in calcium homeostasis, leads to apoptosis [9][10][11]. Lim et al. [12] reported that HD iPSC-derived brain microvascular endothelial cells have intrinsic abnormalities in angiogenesis and blood-brain barrier properties, as well as their signaling pathways. Expression studies from the putamen and cingulate cortex revealed that the dysregulation in the mitochondrial networks and astrocyte function collaterally induces inflammation [13,14]. Seredenina et al. [15] proposed that mutant HTT indirectly affects transcriptional regulation and miRNA binding that consecutively dysregulate the downstream gene expression. Epitranscriptome studies showed that repressing H3K4me3 promotors in human brain cells is beneficial in controlling the abnormal histone acetylation pattern [16,17]. Yildirium et al. [18] observed that aberrant transcriptional changes in H3K27 acetylation lead to changes in downstream expression patterns and an imbalance in histone acetylation and methylation patterns associated with HD. Dalamah et al. [13] performed single-nuclei experiments enriched for astrocyte dynamic states and identified four transcriptionally distinct states of astrocytes, which dysregulates calcium ion channel, synaptic transmission and oxidative phosphorylation.
Recent investigations discussed earlier focused on gene expression change and associated biological functional impact. These studies missed considering the transcript expression changes in HD conditions and how potential GWAS variants affect the expression by modulating regulatory elements and epigenetic motifs. In this work, we quantified the gene/transcript abundance and transcript proportions from BA9 control/HD samples to identify the interlink association between variants gene and transcript expression. In addition, we elucidated the variants from bulk RNA seq and compared them with genomewide association studies with various neurodegenerative disorders [19][20][21][22]. Further, we predicted the effect of variants on regulatory functions such as miRNA, transcription factor (TF) binding, and histone acetylation and methylation pattern.
We found that most variants affect the histone acetylation pattern and are upregulated in HD. This shows an imbalance in the regulatory network, which hampers the acetylation and methylation pattern. Predominant variants overlapped with transcripts and potentially changed their proportions. In addition, we identified variants affecting miRNA binding and associated target gene expression profiles. Differential gene expression studies revealed that mitochondrial, inflammatory-responsive genes are dysregulated at grade 3. Vasculature and angiogenesis-related functional genes are impacted at grade 4. From this study, we propose that neuroinflammation act as a bidirectional sword; it affects the endothelial cell and disrupts the BBB homeostasis. On the other hand, loss of vasculature affects energy metabolism and neurovascular coupling. Restoring vasculature and BBB homeostasis and reducing inflammatory signals may save the neuronal population from cell death.

Dataset and Pre-Processing
The transcriptomic data of the dorsolateral prefrontal cortex (Brodmann Area 9) from the tissue of HD patients and control samples were obtained from the SRA database (SRP051844) [23]. These are human post-mortem samples, and the age of death for HD patients varies from 40 to 75 and 36 to 106 for healthy controls. The HD samples were grouped according to their severity of neuropathological involvement, i.e., Vonsattel grading [24]. Vonsattel grading varies from 0 to 4 based on neuronal loss. Grade 0 has no observable abnormalities, and grade 4 has a 95% loss of neurons. The samples belong to grades 3 and grade 4. The RNA Seq data is in FASTQ format, and it is pre-processed to check the overall sequence quality. The transcriptomic data is prone to read errors, overrepresented sequences, or any experimental artifacts. It is essential to clean the raw data to avoid erroneous results in subsequent analyses, and processing of the raw data is performed using NGSQCToolkit [25], FastQ Groomer [26], Trim Galore [26], and FastQ Trimmer [27] tools. The data quality is assessed using the FASTQC tool, representing the base quality, GC content, duplicates, and adapter contaminations. Data is cleaned and prepared for spliced alignment. The detailed protocol is presented in Supplementary Figure S1.

Alignment: Reference Genome
The pre-processed reads are mapped with the human genome (hg38) using spliceaware aligners. RNA-Seq reads should be mapped carefully because there is a chance that reads spanning splice junctions are incorrectly mapped. Handling alternative splicing in neurodegenerative disorders is especially important. STAR (Spliced Transcripts Alignment to a Reference) [28] is a spliced aligner that uses a 2-way approach (searching and stitching) to optimize the mapping of reads to the reference genome. In this study, the % of uniquely mapped reads ranges from 89 to 94% in HD and 77 to 94% for controls (Supplementary Table S1).

SNP Calling
Mapped alignments were subjected to PCR duplicate removal using Rmdup [29] to avoid misinterpreting an error to be a variant. Variant calling is performed on these files using SAM tools-Mpileup [27] and VarScan [30], and the genomic locations of the variants are obtained using ANNOVAR [31]. We filtered the variants using the following criteria (i) Phred score > 30, (ii) Recalibrate the alignment near the variants, (iii) Ts/Tv > 1.5, and (iv) Minimum read depth is 10.
The variants uniquely found in HD samples were filtered out from the variants of control samples. The majority of these variants are found to be in the non-coding regions, which encompasses regulatory motifs for transcription factors and miRNA binding.

Effect of Variants on miRNA Binding
We compared the identified variants with miRsnpscore [42], a repository of variants of miRNA gene regulation, and polymiRTS [43], a database of variants in experimentally valid seed regions and target sites. The common variants are filtered for further analysis. We also checked whether the miRNAs are reported in previous studies of HD. The miRDB [44] is used in finding the miRNA and gene target interactions that are associated with miRNAs. Later, the gene targets are compared with differentially expressed genes.

Effect of Variants on Regulatory and Epigenetic Modifications
We compared the variants with various XQTL studies available at QTLbase [45] that are associated with neurodegenerative diseases. We performed this comparison to elucidate the role of identified variants on multiple aspects like splicing, expression, methylation and histone acetylation.

Differential Gene/Transcript/Transcript Proportion Expression Analysis
Differential Transcript Expression (DTE) is the transcript level expression in which changes in expression are observed in at least one transcript of the genes between control and disease samples. Differential transcript Usage (DTU) is used to examine the differential transcript proportion (transcript composition of a gene).
To identify the transcript abundance, we used a salmon pseudo lightweight transcriptomebased aligner [46]. Salmon is an alignment-free tool that does its job in a 2-step processindexing and quantification. It provides the count of overlapped reads with transcripts using HG38 transcriptome sequences. We used tximport to convert the transcript abundance to gene quantification. We used DESeq 2 [47] to identify differentially expressed genes. We filtered the DEG based on foldchange (|log 2 foldchange| > 1 and false discovery rate < 0.05). We categorized these genes based on the grade (grade 3/grade 4) of HD samples and compared them with already reported differentially expressed genes of HD from various literature sources and HumanMine [48]. We quantify transcript expression and proportions using DRIMSeq [49] and stageR [50] to understand transcript usage and expression in both HD and control samples.

Elucidating Interlinks between Variant Effect and Transcriptome Expression
The identified variants are majorly located in noncoding regions. To obtain further information about the effect of noncoding variants, we used an integrated tool called SNP nexus and compared variants with GWAS studies. The variants, which are not matching GWAS, are identified as novel variants. Further, we quantified gene/transcript expression from the BA9 data and compared the expression profiles with the overlapping variants. Interestingly, we found that most of the variant-associated genes are differentially expressed at the transcript level (Supplementary Figure S2).

Role of Variants in Epigenetic Modifications
Novel and GWAS-associated variants affect the regulatory functions and epigenetic changes. We compared the variants identified in the present study with various XQTL studies available at QTLbase that are associated with neurodegenerative diseases. We performed this comparison to elucidate the effect of variants in splicing, expression, methylation, and histone acetylation. In addition, we also found that these variant-associated genes are differentially expressed either at the transcript or gene level, as shown in Figure 1. Most of these variants may play an active role in histone modifications and associated expression change in grade 3 and grade 4, as reported in epitranscriptome studies, which also showed similar changes in the histone modification pattern in HD [16][17][18]. Table 1 depicts the effect of variants in epigenetic modification and associated expression patterns in gene and transcript levels in HD. We observed that these variants overlap with the transcript position (obtained with Ensembl [58]), and the transcripts are differentially expressed between control and HD subjects. We have compared our variants with already reported GWAS studies of neurodegenerative disorders, and the variants which are not reported in the GWAS studies are termed "novel." QTL studies showed that the IL18 variant modulates histone methylation, expression and splicing. This variant overlapped with ENST00000280357.11, ENST00000524595.5, ENST00000525547.5 and ENST00000534225.1. We also obtained a similar differential transcript expression pattern.  Table 1 depicts the effect of variants in epigenetic modification and associated expression patterns in gene and transcript levels in HD. We observed that these variants overlap with the transcript position (obtained with Ensembl [58]), and the transcripts are differentially expressed between control and HD subjects. We have compared our variants with already reported GWAS studies of neurodegenerative disorders, and the variants which are not reported in the GWAS studies are termed "novel." QTL studies showed that the IL18 variant modulates histone methylation, expression and splicing. This variant overlapped with ENST00000280357.11, ENST00000524595.5, ENST00000525547.5 and ENST00000534225.1. We also obtained a similar differential transcript expression pattern. IL-18 is involved in neuroinflammation and abundantly expressed by microglia and inflammasome. Recent transgenic mouse studies showed that inhibiting inflammasome reduces IL-18 activation and thus provides neuroprotection [59]. rs360720 (A/G) variant overlapped with the four transcripts and affected their expression pattern in HD BA9 samples. The SLC9A5 gene is involved in Na+/H+ exchange, interacts with the AMPA receptor, and controls ion equilibrium during neurotransmission. Novel variant-associated genes (rs61740454 (G/C)) affect the transcript expression and upregulated HD. MPRIP gene participates in regulating the actin fiber. Studies showed that increased gene expression in neuronal cells affects the actin-binding and stress fiber assembly. MPRIP variants overlap with the transcripts (ENST00000334209.9, ENST00000417669.6, ENST00000579361.1) and modulate the expression. The GAK1 gene plays a vital role in regulating Na+/K+ exchange and maintains the resting membrane potential [60]. The novel variant (rs17165089 (G/C)) affects the transcript expression. Epigenetic roadmap studies showed that it affects the histone modification pattern and is upregulated in HD. GOLGB1, located in the Golgi apparatus, is involved in protein metabolism and regulates the Golgi stress response. Epigenetic studies showed that a novel variant (rs77570895 (C/T)) influences methylation patterns and expression. In addition, this variant affected the transcript (ENST00000491690.1) expression and downregulated in HD. CYB5A is involved in electron transport and metabolism. PDassociated variant rs2032263 (G/C) disrupts the splicing, affecting the transcript expression (ENST00000299438.13, ENST00000397914.4) and also downregulated in HD conditions.
In summary, genes involved in ion homeostasis, stress response, and neuroinflammation are upregulated, and metabolism-related genes are downregulated in HD pathogenesis. In this study, we identified variants overlapping with the transcripts. Further, we found the interlink between the effect of variants in terms of transcript expression and epitranscriptomic modifications.

Effect of Variants on Transcription Factor Binding and Histone Methylation/Acetylation Modification
The variants located in the noncoding region affect the regulatory elements responsible for transcription factor binding (TF), which indirectly affects the downstream gene expression. In Table 2, we listed some crucial variants and affected TFs that are differentially expressed and showed changes in the histone acetylation and methylation profile.
The gene LRRN4CL is responsible for the positive regulation of apoptosis and participates in the tyrosine kinase receptor signaling pathway. This gene is upregulated in HD conditions, and its expression pattern is not reported in the literature. Comparing this variant with QTL, studies showed this variant potentially affects expression and splicing patterns. This GWAS variant rs2512561 (G/A) disrupts the TF binding of NFKB1 and plays a bidirectional role in neuronal survival/death by providing neurotrophic factors and neuroinflammatory response. NFKB1 is downregulated in HD BA9 samples. The B3GNT8 gene present in the exosome membrane helps to transport proteins. PD-associated variant (rs284662 (C/T) affects the histone acetylation and methylation process, affecting the expression pattern. The B3GNT8 gene is upregulated in HD, and its expression pattern is not reported before in HD. The CPZ gene participates in WNT signaling, which plays a crucial role in cell survival/apoptosis. The CPZ novel variant rs13121547 (G/T) affects the histone acetylation pattern, which in turn upregulates the gene expression in HD, and this gene variant/expression is not reported so far in HD. In addition, B3GNT8 and CPZ variants affect the transcriptional repressor CTCF binding, which regulates the histone acetylation pattern. This gene is downregulated in HD, which may lead to the overactivation of the histone acetylation pattern. From our variant analysis, we showed most of the variants (77% and 68% variants in grade 3 and grade 4, respectively), are affecting the histone acetylation process (Figure 1). SLC13A4 is involved in transporter activity, and the GWAS variant (rs3112355 (C/T)) affects the histone acetylation pattern, and this gene is upregulated in HD. This gene variant affects the TF binding of Prrx2 and is involved in vesicle trafficking; this gene is downregulated in BA9 HD, and PRRX2 expression is not explored in HD. PTGDR and PDGFD genes participate in vasodilation and cell migration. GWAS variants (rs4898758, rs11226059) affect the expression pattern, and the MAFF and RUNX1 TF binding is involved in cell proliferation and NOTCH signaling. These genes are upregulated in HD BA9 samples. The DSP gene is involved in the developmental process and the regulation of bundle of His cells to Purkinje myocyte communication. It disrupts the binding of the SOX10 transcription factor, which affects angiogenesis regulation. In this study, we elucidate the relationship between the variant-associated genes and their effect on epigenetic and transcription factor binding, which affects the downstream genes.
Interestingly, about 9.8% of the variants disrupting transcription factor binding affect the binding of novel TF Prrx2, whose binding is affected by 185 and 76 variants in grades 3 and 4, respectively. Prrx2 [61] is a novel gene that is upregulated in HD. It is a gene of paired-related homeobox 2 involved in ERK cascade, PI3K/AKT signaling, and WNT/ ß-catenin signaling pathway and is crucial for cell survival and migration. ERK activation [62] is observed in HD patients, affecting BDNF signaling, apoptosis, glutamate signaling, and EGF signaling. Aberration of PRRX2 may affect neuronal cell survival. In essence, the identified variants affect the transcription factor binding and their respective expression pattern through epigenetic modification.

Effect of Variants on miRNA Binding and Its Expression
Noncoding variants potentially affect the binding motifs of miRNA. We compared our variants with miRsnpscore, a repository of variants affecting miRNAs gene regulation, and polymiRTS, a database of variants in experimentally valid seed regions and target sites. The common variants are filtered for further analysis and examined whether the miRNAs are reported in previous studies of HD. miRDB [44] is used to find the miRNA and gene target interactions. Later, the gene targets are compared with differentially expressed genes. The variants affecting the miRNA binding and their expression are shown in Table 3.
DRAM1 is downregulated in HD. The FARP1 gene participates in synaptic retrograde transport, and the GWAS variant (rs2282048) affects the miR-100 binding and its expression. CSNK1A1 variant modulates the binding of miR-374a, and this miRNA is upregulated in HD, and target gene expression is downregulated. This gene is involved in WNT signaling and is a potential biomarker for Alzheimer's disease. TMEM43 is involved in cardiovascular disease and plays a major role in metabolic pathways; this gene variant changes the binding of miR-30a and is upregulated in HD, and TMEM43 is downregulated in HD BA9 samples. Furthermore, we found other potential target genes for the abovediscussed upregulated miRNA. We noticed that all downstream target genes for these miRNAs were also downregulated in HD.

Differential Gene Expression/Transcript Expression
We performed transcriptome-based quantification using salmon. In grade 4, we found 33 differentially expressed genes reported in HD and 38 novel genes in our study. In grade 3, we found eight differentially expressed genes already reported in HD and 27 novel genes.  We found that the PITX2 gene is important for GABA interneuron development and is upregulated in HD, as reported previously. The HOXA10 gene participates in histone deacetylation and is found to be upregulated in grade 3. CXCR2 and B2M genes play a crucial role in immune response, and this gene is upregulated in grade-3 HD. We have listed all the differentially expressed genes of grade 3 in Supplementary Table S2. In summary, mitochondrial-related genes are downregulated 3, and immune response genes are upregulated in grade 3. This may lead to energy dysfunction and associated neuroinflammation at the early stages of the disease. The normalized TPM count plots for some of the above-discussed genes are represented in Figure 2.

Grade-4 Specific Gene Expression
In grade 4, we found that the TBX15, WNT6 and SFRP2 genes are related to cell fate commitment, and WNT signaling is upregulated and participates in apoptosis. A novel gene, SFRP2, is reported to be negatively regulating the WNT signaling pathway [63], which is crucial for blood vessel formation and essential for vasculature maintenance [64]. OSR1 gene is involved in smooth muscle and pericyte regulation, which are important for blood pressure and vasculature morphogenesis.
CCR2 and GPR183 are involved in astrocyte migration, and these genes are upregulated in HD. The CD177, H19, IGHG3, IGKV3-20 and IGLC2 genes are involved in the inflammatory response, and these genes are upregulated in HD. USP17L2 gene is involved in MAPK and cell surviving signaling, which is downregulated in HD. NPHS1 plays a crucial role in actin filament and mediated transport, and this gene is downregulated in HD. MIR149 is involved in angiogenesis, and vasculature downregulates in HD. The above-discussed gene expression is not reported earlier in the literature. We have listed all the differentially expressed genes of grade 4 in Supplementary Table S3.
In both grades, the MMP9, HOXD10, BARX1, HOXA13 and CPZ genes are upregulated and participate in various functions such as regulation of cytochrome C regulation/apoptosis, motor neuron cell fate regulation and endothelial cell migration, and WNT signaling, respectively. The normalized TPM count plots for some of the above-discussed genes are shown in Figure 3.

Grade-4 Specific Gene Expression
In grade 4, we found that the TBX15, WNT6 and SFRP2 genes are related to cell fate commitment, and WNT signaling is upregulated and participates in apoptosis. A novel gene, SFRP2, is reported to be negatively regulating the WNT signaling pathway [63], which is crucial for blood vessel formation and essential for vasculature maintenance [64]. OSR1 gene is involved in smooth muscle and pericyte regulation, which are important for blood pressure and vasculature morphogenesis.
CCR2 and GPR183 are involved in astrocyte migration, and these genes are upregulated in HD. The CD177, H19, IGHG3, IGKV3-20 and IGLC2 genes are involved in the inflammatory response, and these genes are upregulated in HD. USP17L2 gene is involved in MAPK and cell surviving signaling, which is downregulated in HD. NPHS1 plays a crucial role in actin filament and mediated transport, and this gene is downregulated in HD. MIR149 is involved in angiogenesis, and vasculature downregulates in HD. The above-discussed gene expression is not reported earlier in the literature. We have listed all the differentially expressed genes of grade 4 in Supplementary Table S3.
In both grades, the MMP9, HOXD10, BARX1, HOXA13 and CPZ genes are upregulated and participate in various functions such as regulation of cytochrome C regulation/apoptosis, motor neuron cell fate regulation and endothelial cell migration, and WNT signaling, respectively. The normalized TPM count plots for some of the above-discussed genes are shown in Figure 3. Genes 2022, 13, x FOR PEER REVIEW 13 of 23

Comparison of Gene Expression (BA9/HD) with Literature Bulk RNA-Seq
To validate the expression pattern, we also compared the identified gene expression profiles of BA9 HD samples with other published bulk RNA seq studies in HD and the results are shown in Figure 4.

Comparison of Gene Expression (BA9/HD) with Literature Bulk RNA-Seq
To validate the expression pattern, we also compared the identified gene expression profiles of BA9 HD samples with other published bulk RNA seq studies in HD and the results are shown in Figure 4.
It is crucial to know that the genes in HD samples are exhibiting differential expression as in other reported bulk RNA Seq studies. This helps to confirm that the up/downregulated genes of our study are also expressed in other studies. Most of the gene expression profiles are matched with the reported studies.

Comparison of Gene Expression (BA9/HD) with GTEX (BA9) and Other Tissue
We compared the novel gene expression pattern in HD with huge control samples from BA9 and various tissues (GTEX consortium data) to examine whether the disease sample gene expression differs from the control tissues. The comparison of novel gene expression patterns with GTEX and other tissue control samples is shown in Figure 5.
It is observed that the gene expression pattern of HD differs from that of the controls and also normal samples from other brain tissues. We found that the control BA9 TPM values match the GTEX BA9 and exhibit similar expression patterns. The novel gene expression in HD differed from that of the other control expression, which shows that these genes are potentially dysregulated in HD BA9.  It is crucial to know that the genes in HD samples are exhibiting differential expression as in other reported bulk RNA Seq studies. This helps to confirm that the up/downregulated genes of our study are also expressed in other studies. Most of the gene expression profiles are matched with the reported studies.

Gene Co-Expression Network and Function Module Network for Differentially Expressed Genes
We have constructed a co-expression network ( Figure 6) for the differentially expressed genes using COEXPRESSdb [53]. Genes that are co-expressed together tend to have similar functions. Here we discuss some of the DEGs based on the network properties (degree and betweenness centrality). MRC1 is a novel HD gene co-expressed with CD163, VSIG4, and MS4A4A, and these genes are involved in vasculature regulation. ANXA2 is an upregulated gene co-expressed with S100A4 and S100A6 genes, and these genes are active participants in interleukin signaling that promotes neuroinflammation in HD [65]. GZMK is also co-expressed with PTGDR, which has a major role in MAPK cascade and immune response. GZMK and PTGDR are upregulated, and their involvement in MAPK signaling can be considered as active participation in cell survival mechanisms [66]. OSR1 gene is involved in pericyte migration and maintains the vasculature, and this gene co-expressed with the LRRN4CL, FOXF2 and OSR2 genes. Novel TF FOXD2 is co-expressed with FOXD1 and WNT6. We compared the novel gene expression pattern in HD with huge control samples from BA9 and various tissues (GTEX consortium data) to examine whether the disease sample gene expression differs from the control tissues. The comparison of novel gene expression patterns with GTEX and other tissue control samples is shown in Figure 5. It is observed that the gene expression pattern of HD differs from that of the controls and also normal samples from other brain tissues. We found that the control BA9 TPM values match the GTEX BA9 and exhibit similar expression patterns. The novel gene expression in HD differed from that of the other control expression, which shows that these genes are potentially dysregulated in HD BA9.

Gene Co-Expression Network and Function Module Network for Differentially Expressed Genes
We have constructed a co-expression network ( Figure 6) for the differentially expressed genes using COEXPRESSdb [53]. Genes that are co-expressed together tend to have similar functions. Here we discuss some of the DEGs based on the network properties (degree and betweenness centrality). MRC1 is a novel HD gene co-expressed with

Tissue-Specific Function-Interaction Network
We constructed a functional interaction network using novel differentially expressed genes and transcription factors (Figure 7) and associated downstream gene expression. We have used ReactomefiVIZ [51] to construct a functional interaction network. This tool uses interactions obtained from pathways in the KEGG and Reactome knowledgebase and merges them with predicted interactions using a machine learning approach to annotate function interactions. Figure 7 shows that TCF12 plays a significant role in cell fate determination. This TF binding was affected by EP300 and FOXD1 novel variants. Novel DEGs such as HOXA11, HOXD8, HOXD9, NKX2-5, HOXD10, and FOXD1 genes are activated by TCF12. All these downstream genes are dysregulated in HD. We predict that novel differentially expressed genes KRT19 and TBX15 may be activated by ESR1 and EP300, respectively, based on the information obtained from function interactions using the ReactomeFI tool [51]. EP300 also regulates the expression of TCF12, an upregulated transcription factor in HD. EP300 is a variant-associated gene known for its function in histone acetylation [67] and plays a significant role in epigenetic modifications and gene expression [68].
ANXA2 is an upregulated gene co-expressed with S100A4 and S100A6 genes, and these genes are active participants in interleukin signaling that promotes neuroinflammation in HD [65]. GZMK is also co-expressed with PTGDR, which has a major role in MAPK cascade and immune response. GZMK and PTGDR are upregulated, and their involvement in MAPK signaling can be considered as active participation in cell survival mechanisms [66]. OSR1 gene is involved in pericyte migration and maintains the vasculature, and this gene co-expressed with the LRRN4CL, FOXF2 and OSR2 genes. Novel TF FOXD2 is coexpressed with FOXD1 and WNT6.

Tissue-Specific Function-Interaction Network
We constructed a functional interaction network using novel differentially expressed genes and transcription factors ( Figure 7) and associated downstream gene expression. We have used ReactomefiVIZ [51] to construct a functional interaction network. This tool uses interactions obtained from pathways in the KEGG and Reactome knowledgebase and merges them with predicted interactions using a machine learning approach to annotate function interactions.  Figure 7 shows that TCF12 plays a significant role in cell fate determination. This TF binding was affected by EP300 and FOXD1 novel variants. Novel DEGs such as HOXA11, HOXD8, HOXD9, NKX2-5, HOXD10, and FOXD1 genes are activated by TCF12. All these downstream genes are dysregulated in HD. We predict that novel differentially expressed genes KRT19 and TBX15 may be activated by ESR1 and EP300, respectively, based on the information obtained from function interactions using the ReactomeFI tool [51]. EP300 also regulates the expression of TCF12, an upregulated transcription factor in HD. EP300 is a variant-associated gene known for its function in histone acetylation [67] and plays a significant role in epigenetic modifications and gene expression [68].
These interaction networks revealed the functional relationship between DEGs and TFs. Variant genes are also indirectly affecting these functional links as they alter the binding of transcription factors that are majorly involved in the histone acetylation process. Modulation in the TF network affects the downstream target genes. These interaction networks revealed the functional relationship between DEGs and TFs. Variant genes are also indirectly affecting these functional links as they alter the binding of transcription factors that are majorly involved in the histone acetylation process. Modulation in the TF network affects the downstream target genes.

Tissue-Specific Functional Module Network
Differentially expressed genes were subjected to tissue-specific functional module enrichment analysis. Modules represent different functions in which the genes are involved. The network of tissue-specific function modules for BA9 is shown in Figure 8. These mod-ules consist of genes responsible for angiogenesis, vasculature development, epithelial cell differentiation, immune response, and G-protein coupled receptor signaling pathway. Most of the genes involved in the vasculature and blood vessel-related functions are dysregulated in HD pathogenesis. This study showed that most of the genes related to immune response are differentially expressed, which may lead to aberrant neuronal inflammation. We also found that genes involved in vasculature maintenance are dysregulated, which may affect the blood-brain integrity and add stress burden leading to inflammation.

Proposed Mechanism for Neurodegeneration
Neuroinflammation plays a bifacial role in cell survival and death. Massive amounts of neuroinflammatory cytokines released from hypoxic/injured neurons activate the inflammatory microglia and reactive astrocytes (Figure 9). By analyzing the microenvironment, microglia and astrocytes provide trophic or cytotoxic signals. This massive cytokine rush affects the endothelial structure and distorts tight junctions and pericytes in blood vessels. Vascular abnormalities affect smooth muscle control that eventually strikes vasodilation and constriction. Loss of coordination between neurons, blood vessels and astrocytes hampers the energy metabolism process. From the transcriptome and functional enrichment analysis, we propose that aberration in the vascular region and neuroinflammation affects the structural integrity, which hampers the metabolic pathways, making the associated neurons vulnerable in HD.

Proposed Mechanism for Neurodegeneration
Neuroinflammation plays a bifacial role in cell survival and death. Massive amounts of neuroinflammatory cytokines released from hypoxic/injured neurons activate the inflammatory microglia and reactive astrocytes (Figure 9). By analyzing the microenvironment, microglia and astrocytes provide trophic or cytotoxic signals. This massive cytokine rush affects the endothelial structure and distorts tight junctions and pericytes in blood vessels. Vascular abnormalities affect smooth muscle control that eventually strikes vasodilation and constriction. Loss of coordination between neurons, blood vessels and astrocytes hampers the energy metabolism process. From the transcriptome and functional enrichment analysis, we propose that aberration in the vascular region and neuroinflammation affects the structural integrity, which hampers the metabolic pathways, making the associated neurons vulnerable in HD.

Limitations
This limitation of the study is the small size of data (population of HD samples) obtained from BA9. Although we have compared the quantified HD expression with GTEX and other bulk RNA seq data, we still need a large population to understand the variant effects in expression patterns.

Conclusions
In this study, we explored the effect of variants on transcript expression and splicing patterns. Most of the identified variants are located in the noncoding region and modify the regulatory elements. Novel variants associated with the following genes, GAK1, SLC9A5, MPRIP, GOLGB1 and GAK1, affect the histone acetylation pattern. However, these variants are overlapped with their transcript and are upregulated in HD at both gene and transcript levels. Seven novel variants impact the transcription factor binding of HOXA5, FOXD1 and PRRX2 and play a crucial role in homeostasis and neuromuscular signaling pathways. These TFs are dysregulated in HD and have not been reported so far. Interestingly, we found a novel variant in the CPZ gene that affects the histone acetylation pattern and is upregulated in HD. The CPZ variant modulates the CTCF binding and is a well-known histone acetylation repressor. The CTCF gene is downregulated in HD. We also found the variant genes affecting miRNAs and their downstream target genes. Identified variants affect the miRNA binding. These miRNAs are upregulated in HD, and their downstream target genes are downregulated in HD.
Differential gene expression analysis showed that mitochondrial function-related genes and inflammatory-responsive genes are dysregulated at grade 3. Blood vessel morphogenesis, cell survival signaling, and angiogenesis genes are downregulated in grade 4; Inflammatory related genes are tremendously upregulated in grade 4. 33 genes of grade 4 differentially expressed genes are upregulated in order to activate Aryl hydrocarbon

Limitations
This limitation of the study is the small size of data (population of HD samples) obtained from BA9. Although we have compared the quantified HD expression with GTEX and other bulk RNA seq data, we still need a large population to understand the variant effects in expression patterns.

Conclusions
In this study, we explored the effect of variants on transcript expression and splicing patterns. Most of the identified variants are located in the noncoding region and modify the regulatory elements. Novel variants associated with the following genes, GAK1, SLC9A5, MPRIP, GOLGB1 and GAK1, affect the histone acetylation pattern. However, these variants are overlapped with their transcript and are upregulated in HD at both gene and transcript levels. Seven novel variants impact the transcription factor binding of HOXA5, FOXD1 and PRRX2 and play a crucial role in homeostasis and neuromuscular signaling pathways. These TFs are dysregulated in HD and have not been reported so far. Interestingly, we found a novel variant in the CPZ gene that affects the histone acetylation pattern and is upregulated in HD. The CPZ variant modulates the CTCF binding and is a well-known histone acetylation repressor. The CTCF gene is downregulated in HD. We also found the variant genes affecting miRNAs and their downstream target genes. Identified variants affect the miRNA binding. These miRNAs are upregulated in HD, and their downstream target genes are downregulated in HD.
Differential gene expression analysis showed that mitochondrial function-related genes and inflammatory-responsive genes are dysregulated at grade 3. Blood vessel morphogenesis, cell survival signaling, and angiogenesis genes are downregulated in grade 4; Inflammatory related genes are tremendously upregulated in grade 4. 33 genes of grade 4 differentially expressed genes are upregulated in order to activate Aryl hydrocarbon receptor (Ahr) transcription factor that is responsible for the degeneration of motor neurons and behavioral activity related neurons. The removal of Ahr is beneficial and gives a neuroprotective effect to the HD samples, as previously reported [69].
We explored the BA9-specific function interaction networks to understand the relationship between novel differentially expressed genes and transcription factors. Tissue-specific functional enrichment showed that predominant genes related to vascular architecture are dysregulated in HD. From this detailed analysis, we propose that abnormal neuroinflammation affects endothelial cell integrity, affecting BBB homeostasis.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/genes13122385/s1, Figure S1: Protocol followed for variant calling, differential gene expression and network analysis of RNA-Seq Huntington's data; Table  S1: Dataset information along with uniquely mapped reads rate; Figure S2: Variants matching with differential gene expression and transcript expression; Figure S3: Transcript expression and proportion plots of grade 3 and grade 4 variant associated genes; Table S2: Differentially expressed genes of grade 3; Table S3: Differentially expressed genes of grade 4.  Data Availability Statement: Data will be made available upon request.