This study utilized publicly accessible sequencing data obtained from HD gene-positive and neurologically normal individuals in the Sequence Read Archive (SRA) public repository database. Consequently, formal ethical approval was deemed unnecessary for this investigation.

The sequencing data of the dorsolateral pre-frontal cortex (BA9) from both HD gene-positive and control samples was obtained from the SRA database (BioProject PRJNA271929, available on https://www.ncbi.nlm.nih.gov/bioproject/PRJNA271929). This BioProject is comprised of 69 postmortem samples, featuring 20 from HD gene-positive individuals aged between 40 and 75 years (with a mean of 58.2 ± 10.4 years) and 49 from neurologically normal control human individuals aged between 36 and 106 years (with a mean of 68.3 ± 15.8 years). The selection of this BioProject for our study was based on the following criteria: (i) the presence of a satisfactory number of case-control samples, (ii) the high quality of sequencing data, and (iii) the prior analysis of these samples by Labadorf et al. (Labadorf et al., 2015), followed by subsequent reanalysis, including additional BA9 and/or CAU samples (Labadorf et al., 2018; Agus et al., 2019). Demographic data from the samples are described in Supplementary Excel S1.

To assess our hypothesis that the inappropriate use of older control samples can affect DEGs identification, we employed propensity score matching (PSM) at a 1:1 ratio (HD gene-positive subject: neurologically normal subject). PSM is a quasi-experimental method, initially introduced by Rosenbaum and Rubin (1983), that aims to align affected/treated and control groups based on a targeted feature to enhance comparability. For this purpose, we utilized the MatchIt package in R (Ho et al., 2007). To validate the outcomes derived from the package, the age distribution was visually examined before and after PSM through both boxplot and density plot analyses. Additionally, the statistical confirmation of age distribution equality was carried out using a t-Student test, all executed within the R environment.

The RNA-Seq data, formatted in FASTQ, underwent pre-processing through FastQC (Andrews, 2010) and MultiQC (Ewels et al., 2016) tools to ensure sequencing quality. Subsequently, the pre-processed reads were mapped to the latest human genome (hg38) using STAR (Spliced Transcript Alignment to a Reference), a splice-aware aligner designed for accurately aligning reads to the reference genome (Dobin et al., 2013). Transcript abundance was then estimated utilizing the feature count read summarization program (Liao et al., 2014). The abundance estimates from all samples were consolidated into a unified expression matrix and normalized using the DESeq2 package v1.10.1 (Love et al., 2014), following the methodology outlined by Labadorf et al. (2015).

To gain deeper insights into how the age of control samples may impact the identification of differentially expressed genes, the normalized read counts of HD gene-positive individuals (n = 20) were individually compared with both the entire control sample set (n = 49) and the PSM-selected control samples (n = 20, hereafter referred to as Age-matched control samples). The differential expression analysis was performed using the DESeq2 method (Love et al., 2014), generating a comprehensive list of genes along with their respective levels of differentiation and statistical significance.

Mapped genes were categorized based on their abundance, p-value, and fold change (|log2FC|) into six distinct categories, as outlined in Table 2; Figure 1. Genes classified as upregulated (URG) or downregulated (DRG) were considered differentially expressed genes (DEGs), while other categories were designated as non-differentially expressed genes (NDEGs). Notably that both zero count genes (ZCG; genes with row counts = 0) and low count genes (LCG; genes with normalized counts < 10 in less than 70% of the smaller category) are excluded from the differential gene expression analysis using DESeq2.

For a comprehensive comparative assessment of the entire transcriptome between HD gene-positive individuals and all neurologically normal subjects, we employed the Uniform Manifold Approximation and Projection (UMAP) dimensionality reduction technique (McInnes et al., 2018), which proved itself to produce insightful reduced dimensions to represent genomic data (Dorrity et al., 2020). Utilizing three distinct components from UMAP, we subsequently applied the Density-based Spatial Clustering of Applications with Noise (DBSCAN) algorithm to unveil primary grouping patterns (Kriegel et al., 2011); the advantage of this approach is the fact that no prior knowledge on the ideal number of clusters is required; also, DBSCAN generates clusters based on density of data in a particular region and thus automatically recognizes anomalous or outlying observations. The UMAP analysis was performed in CRAN-R environment through UMAP package, while the subsequent application of DBSCAN was applied using the scikit-learn-based DBSCAN implementation for Python 2.6.

Recognizing that the biological effects on the HD phenotype hinge on the intricate interplay among DEGs, we subjected the ranked list of genes—comprising both URGs and DRGs—derived from the analysis using all control samples (n = 49) and age-matched controls (n = 20) to distinct functional enrichment analyses. These analyses were conducted using the Gene Set Enrichment Analysis (GSEA) to explore biological pathways (BP) and molecular functions (MF).

To assess the influence of aging on genes that have become DEGs or NDEGs, we employed a linear regression model on the median of the reads by age across all control samples. This approach allowed us to examine trends and determine their statistical significance. These analyses were conducted within the R environment, utilizing the stats package.

Given that the raw p-values obtained through the Wald test (employed for DEG identification by DESeq2) are subjected to correction for multiple testing using the Benjamin and Hochberg (BH) method to control the false-discovery rate (FDR) (Benjamin and Hochberg, 1995) and considering that age matching primarily affects the genes identified as DEGs, we further hypothesized that age matching could potentially influence DEG discovery by altering the p-adjusted values. To investigate this hypothesis, we compared the raw p-values and the p-adjusted values of genes that ceased to be identified as DEGs upon excluding control samples from older neurologically normal individuals.

Quality control analysis stands as an indispensable facet of RNA-Seq analysis (Consortium, 2011; Conesa et al., 2016; Chung et al., 2021; Faustinoni-Neto et al., 2023). Therefore, we meticulously evaluated the quality control of the 69 BioSamples by utilizing FastQC and MultiQC tools. The outcomes unequivocally affirmed that all BioSamples within the BioProject PRJNA271929 exhibited a median per-sequence Phred score of 36, signifying a sequencing accuracy surpassing 99.9%. Furthermore, the samples displayed a %GC content of 49%, indicating the absence of potential contaminations. These findings reassert the high quality of all 69 BioSamples, aligning with prior validations conducted by Labadorf et al. (2015).

As anticipated, the PSM method effectively selected control samples from neurologically normal individuals with a comparable age distribution to HD gene-positive individuals (Figure 2). Consequently, two distinct control groups were established: (i) the “All Controls” group, encompassing all neurologically normal individuals (n = 49, with an average age of 68.3 ± 15.8 years), replicating the control cohort employed in the study by Labadorf et al. (2015), and (ii) the “Age-Matched” group, formed by samples from neurologically normal individuals exhibiting a similar age distribution (n = 20, with an average age of 58.2 ± 10.4 years) to that of the HD gene-positive individuals (n = 20, with an average age of 58.8 ± 10.5 years), representing the “virtual” age-matched cohort.

To assess the potential impact of utilizing control samples from older individuals on differential expression analysis, we initially compared the entire transcriptome of HD gene-positive individuals with that of all neurologically normal individuals using UMAP.

The three UMAP-components effectively grouped most samples into seven distinct density-based clusters, revealing three overall patterns (Figure 3A). Notably, clusters C3 and C4 on the left side predominantly comprised young control individuals. In contrast, the transcriptome of HD gene-positive individuals (clusters C6 and C7) occupied an intermediary position on the Cartesian plane between neurologically normal young individuals and more distant older control individuals (aged over 70 years, Figures 3C,D). This outcome strongly suggests that age significantly influences the BA9 transcriptomic profile, emphasizing that using older individuals as controls may not be appropriate for identifying DEGs in HD.

To assess the impact of including older control individuals on the identification of DEGs, we conducted a comparison by analyzing the number of genes identified across various classes (ZCG, LCG, EEG, NSLFCG, URG, and DRG) in HD gene-positive individuals relative to (i) the entire control group (n = 49) and (ii) age-matched controls (n = 20) using a cross table. The results revealed that excluding older individuals from the control group led to the reclassification of 1,915 genes of interest. Among these, 1,523 (79.5%) were no longer considered DEGs (putative false-positive genes), while 392 (20.5%) transitioned to being classified as DEGs (putative false-negative genes) (Figure 4, highlighted in color in Table 3).

Recognizing that different genes can concurrently regulate multiple molecular functions in various biological processes, the DEGs identified in HD gene-positive individuals concerning all controls (n = 49) or age-matched controls (n = 20) were independently subjected to functional enrichment analysis using GSEA in terms of Geneontology (GO). The results revealed that protein-coding DEGs from all controls enriched for 162 biological pathways, while those obtained exclusively from age-matched controls enriched for 137 biological pathways (Figure 5; Supplementary Excel S2). Upon comparison, it was observed that 77 pathways (Supplementary Table S1) were no longer identified as enriched when samples from older neurologically normal individuals were excluded from the control group (Figure 5).

Among the pathways excluded is the one associated with cellular response to heat (Supplementary Table S1), encompassing heat shock genes previously reported as upregulated in HD^10–12. Conversely, the omission of samples from older neurologically normal controls enriched 52 new pathways (Figure 5; Supplementary Table S2). A total of 85 pathways remained consistent between the analyses using all and age-matched controls (Figure 5; Supplementary Table S3). Interesting, we also observed that the genes that were no longer classified as DEG (highlighted in red in Table 3) are directly involved in inflammatory process (Supplementary Table S4), reinforcing that the aging exacerbates (neuro)inflammatory process which can lead to false-positive results. These findings underscore that the inclusion of samples from older control individuals not only leads to the misidentification of DEGs but also exerts a detrimental impact on functional enrichment.

Considering the involvement of DEGs in HD pathophysiology, mediated by protein interactions within each biological pathway, we identified the protein-coding genes associated with pathways no longer enriched (235 genes, 77 excluded pathways) and those enriched in new pathways (222 genes, 52 pathways) identified in the analysis using age-matched controls. To assess the impact of aging on the expression of these genes, we categorized them based on their normalized counts into four groups reflecting their abundance levels (Table 4). This categorization facilitated the visualization of gene expression levels across age groups. We compared normalized counts of these genes (per category) between HD gene-positive and neurologically normal individuals aged less than 60, 60–75, and older than 75 years (the maximum age observed in the HD group, Figure 2). The rate estimates of expression changes across age in the control groups for the genes that became DEGs showed statistical significance.

These results highlight that aging has a more pronounced effect on genes that became DEGs, reducing the expression levels of DRGs and increasing the expression levels of URGs in HD gene-positive individuals compared to neurologically normal controls aged under 60 years (Figure 6). Notably, the level of gene expression abundance correlates with the magnitude of changes observed across age groups.

These findings collectively affirm that aging significantly impacts gene expression in BA9. This underscores the importance of utilizing control samples with a comparable age distribution to HD gene-positive individuals, emphasizing that such an approach is preferable and imperative for accurately identifying DEGs with potential utility as pharmaceutical or prognostic targets. However, it’s worth noting that aging exhibits negligible effects on genes that are no longer URGs and have subtle effects on genes that are no longer DRGs (Figure 6).

The impact of age matching on p-value adjustment is evident in the reclassification of numerous DEGs to NDEGs (Supplementary Figure S1). To illustrate these findings, we examined the normalized counts of two genes (previously identified as upregulated by Labadorf et al. (2015)), HILPDA and SERPINH1, which were no longer classified as DEGs in the analysis with age-matched controls. The raw p-values obtained using both all controls and age-matched controls were statistically significant (p-value < 0.05). However, the p-value adjustment applied in DESeq2, which controls the false-discovery rate (FDR) and considers age-matched controls, resulted in a non-significant adjusted p-value (p-adjust > 0.05, Figure 7). This effect is observed in 3.206 genes (2.686 protein-coding genes), from which 478 genes had log2FC > 0.58 and 310 had log2FC < −0.58 (Supplementary Excel S3). These results underscore that the absence of age matching leads to the identification of false-positive DEGs due to p-value adjustment, reinforcing the findings in Figure 6.

Furthermore, we demonstrated that age matching also influences the log2 fold change (log2FC), altering the gene abundance levels of control samples, as depicted in Figure 6. As anticipated, this effect was more pronounced for genes that became DEGs and those that are no longer DEGs (Supplementary Figure S1) than those that remain classified as DEGs (Supplementary Figure S1). Collectively, these findings strongly support the contention that the lack of age matching leads to the misidentification of DEGs.

In our final assessment, we examined the expression levels, based on normalized counts, of genes previously identified as URGs in HD gene-positive individuals by Labadorf et al. (2015), using both all controls and age-matched controls. Interestingly, we noted that the heat shock protein-coding genes identified as URGs in HD by Labadorf et al. (2015) exhibit relatively low expression, with a median expression in HD gene-positive individuals falling below the commonly used threshold for defining gene expression (count > 10, Figures 8A,B). This observation suggests that, in addition to age matching in the control group, normalized counts should be further scrutinized with adjusted p-values and log2FC for accurate DEG identification.

Consequently, through the integration of adjusted p-values, log2FC, and normalized counts, we identified ADRA2B as upregulated in HD gene-positive individuals (Figure 8C), thereby emerging as a potential candidate for HD treatment. This underscores the importance of a comprehensive analytical approach for robust and accurate DEG identification. The genes identified as putative DEGs in HD (in relation to age-matched controls) can be visualized in Figure 9; Supplementary Excel S3.