Genetic modifiers of Huntington’s disease differentially influence motor and cognitive domains

Subjects and clinical phenotypes

HD subjects in our modifier GWAS and genome-wide genotyping were described previously.^{2; 3} These genetic studies have been approved by the Institutional Review Board of Mass General Brigham (formerly Partners Healthcare). The phenotypic data were from the components of the UHDRS as recorded by trained study examiners in the PREDICT-HD study and the Cohort study of the Huntington Study Group, the Registry study of the European Huntington’s Disease Network and the ENROLL-HD platform. For the vast majority of participants, the age-at-onset of motor signs was that recorded by trained raters. For a small number of participants, a rater estimate was not available, so a family estimate was used. TMS, which can range from 0 (normal) to 124 (severest impairment), is the sum of items from the motor-assessment scale. This comprises the examiners’ evaluation of, among others, oculomotor function, dysarthria, chorea, dystonia, bradykinesia, gait, and postural stability. The language independent Symbol Digit Modalities Test (SDMT), adapted from the Wechsler digit symbol subtest, aims to measure working memory, complex scanning, and processing speed (but can also be influenced by motor dysfunction). Participants match symbols with numbers in a timed process whereby appropriate numbers are placed below the matching symbols as quickly as possible, generating raw scores based upon the number of items correctly completed in 90 seconds.

Algorithmic prediction of age at HD landmarks

In our previous GWAS of HD onset, we relied on age at motor onset, which limits the sample size to those individuals who have passed a disease landmark whose assessment, by its nature, involves subjective judgment and, if judged retrospectively, potential recall bias. Consequently, we have explored the use of predictive modeling to generate GWAS phenotypes that integrate multiple individual standardized measures and reduce subjectivity to increase our power to detect associations. Computing a predicted phenotype for GWAS in HD is challenging. The number of individuals for which a phenotype is predicted should be maximized, while accounting for the vagaries of observational data. When the phenotype is a landmark event such as motor diagnosis, several difficulties present themselves. Participants can enter a study having experienced the event in the past, with the event timing being vague or unknown (left-censoring). Alternatively, individuals can enter the study and not experience the event in the time that they are observed (right-censoring). Visits tend to be annual and often irregular, so that the exact event timing can only be inferred to occur sometime between two visits (interval-censoring). An additional challenge is accounting for different ages at study entry when attempting to predict the age at which an event occurs. Age is the natural time metric for analysis with time 0 being birth. No participant has been observed at time 0 in an observational study, meaning that individuals have postponed or delayed entry (left-truncation), which requires proper identification of who is at risk of an event over what age interval. Finally, measurement of the event itself is tainted by error (e.g., natural waxing and waning of motor signs), and the amount of observed information (e.g., due to number of visits) among participants can vary widely. The above characteristics are very difficult to handle in a single statistical model. Therefore, for the primary predicted phenotypes of age-at-DCL4 and age-at-TFC6, we developed an algorithm for predicting age at each landmark event using both survival analysis and longitudinal analysis, described in Figure S1. We emphasize that the algorithm is intended to provide relatively accurate predictions for individuals in a manner that facilitates research and is not intended to make group-level inferences. From data on > 19,000 participants available through the Cohort, Predict-HD, Registry and Enroll-HD natural history studies, we predicted the age at each event and then utilized the subset of participants with available genome-wide genotyping data to perform association analyses.

In addition to age-at-DCL4 and age-at-TFC6, we defined additional landmarks based upon TMS or SDMT scores. We first determined the mean age-at-DCL4 and mean age-at-TFC6 for individuals at each CAG repeat length from 40 to 50 and compared these to the distribution of TMS and SDMT scores at each CAG length. For both TMS (scale: 0-107 in our sample) and SDMT (scale 0-110 in our sample), a score of 30 typically fell at an age between the median age at DCL4 and TFC6. Consequently, we selected this value for prediction of TMS and SDMT landmarks that occur within roughly the same age range of early manifest HD. Once the threshold value of 30 was decided, the estimated ages at which an individual progressed up to TMS = 30 (‘age-at-TMS30’) and down to SDMT = 30 (‘age-at-SDMT30’) was computed based on the third and fourth panels of Figure S1, except that linear mixed models for SDMT and TMS were used, respectively (rather than the model for cUHDRS from panel two).

Selection of subjects for analysis of TMS and SDMT phenotypic extremes

TMS and SDMT phenotypes were classified as extreme as previously described.^{25; 26} In brief, TMS and SDMT collected as part of the UHDRS were classified as extreme if they were outside of predefined cut-offs on a distribution of that phenotype. Statistically, extremes of the phenotype distribution were derived from quantile regression (a linear quantile mixed model approach).²⁷ For the longitudinal TMS and SDMT data, the quantitative phenotypes were classified as extreme at each visit if measures occurred outside of the cut-offs. The TMS data distribution was modelled as an exponential function of age stratified by CAG length. SDMT was examined on a linear function of age in ∼10-year age bins (18–29; 30–39; 40–49; 50–59; 60–69; 70–79; ≥ 80) controlling for CAG length and education. A single extreme measure may not mean the phenotype remains consistently extreme over time. Hence, participants were classified as being phenotypic extremes if they had an extreme phenotype, motor or cognitive, in the same direction at two or more consecutive visits. We defined the extremes as those below a 5% cut-off or above a 95% cut-off for TMS and below a 10% cut-off and above a 90% cut-off for SDMT in order to obtain roughly comparable numbers of genotyped subjects, since only a fraction of the study participants have genome-wide genotyping data.

Imputation and association analyses

Five batches of GWA data sets were produced and subsequently imputed using the HRC reference panel as described previously.³ Genotype imputation using the TOPMed reference panel was performed similarly. Briefly, the individual genotype data set (i.e., GWA1, 2, 3, 4, and 5) comprising subjects with genotype call rate greater than 90% and single nucleotide variants (SNVs) with call rate > 95% and MAF > 1% was subjected to quality control by the “HRC or 1000G Imputation preparation and checking” program. Additional QC was performed using the TOPMed Imputation Server, revealing genomic regions of 10 MB with at least one sample with call rate < 50%. Those low call rate samples at certain genomic regions were further excluded from genotype imputation. The final QC-passed typed data set was used for imputation by the TOPMed Imputation Server. Genomic coordinates of TOPMed imputation data were converted to GRCh37/hg19 to make this data set directly comparable to the HRC imputation data set. Finally, we removed SNVs with imputation r-square < 0.5 in any of the GWA data sets, HWE p-value < 1E-6 (excluding chr4:1-5000000 region containing HTT), or MAF < 0.1%. For a given phenotype, genome-wide association analysis was performed for individuals of European ancestry using a mixed effects model with a relationship matrix and a set of covariates such as sex, batch, and the first four principal component values from genetic ancestry analysis as covariates.³ For selected candidate regions with significant association signals, conditional analysis was performed using a fixed effect linear model that included one or more additional covariates representing the minor allele counts of conditioned SNVs. SNVs are expressed as chromosome_coordinate_reference-allele_alternative-allele based on the GRCh37/hg19 genome assembly. CAG repeat lengths for each subject were generated by a standard ABI fragment analysis assay.^{3; 28} For individuals uncovered in the TopMed imputation with HTT CAA-loss and CAACAG-duplication haplotypes, whose CAG lengths are systematically mis-estimated by the fragment analysis assay, the pure CAG repeat size was determined by MiSeq (Illumina) DNA sequencing analysis using the previously reported assay method.³

Gene-wide and Pathway Analysis

Gene-wide association analyses were carried out in MAGMA²⁹ on summary statistics from each GWAS using genotypes from 7,086 individuals from GWA3, 4 and 5 as a reference panel to estimate LD. As in our previous HRC-based analysis³, we used a window of 35kb upstream and 10kb downstream of gene positions (GRCh37/hg19) along with the ‘‘multi’’ analysis option, combining the mean SNV p-value with the top SNV p-value (corrected for number of SNVs and the LD between them). Pathway enrichment analyses were also performed in MAGMA, correcting for LD between genes, initially using ‘‘self-contained analysis’’ measuring overall association among genes in a pathway, and then using a more conservative ‘‘competitive’’ analysis, to compare association in genes within a pathway to those outside the pathway.

Pathway definition

The assignment of Gene Ontology (GO) terms to human genes was obtained from the “gene2go” file, downloaded from NCBI on March 11^th, 2020. “Parent” GO terms were assigned to genes using the ontology file downloaded from the Gene Ontology website on the same date. GO terms were assigned to genes based on experimental or curated evidence, excluding electronic or computational annotations. Reactome pathways were downloaded on April 26^th, 2020. Biocarta, KEGG and Pathway Interaction Database (PID) pathways were downloaded from v7.1 (March 2020) of the Molecular Signatures Database.³⁰ Analysis was restricted to GO terms containing between 10 and 2000 genes. No size restrictions were placed on the other gene sets, since these were fewer in number. This resulted in a total of 10,041 gene sets for analysis.

Statistics and software

TOPMed genotype imputation was based on the TOPMed Imputation Server. Mixed effects model GWAS was performed using the gemma program version 0.94 beta.³¹ Genomic inflation factors were calculated using the R package ‘GenABEL’ (version, 1.8-0).³² Estimated age-at-landmarks were developed using the R package flexsurv³³ for survival analysis and lme4 for longitudinal linear mixed models.³⁴ Quantiles for extreme values of the TMS and SDMT were determined using the R package lqmm.²⁷ Gene-wide and pathway analyses were performed with MAGMA v1.08. Gene Ontology annotations were downloaded from the NCBI and Gene Ontology websites. Reactome pathways were downloaded from the Reactome site. Biocarta, KEGG and Pathway Interaction Database (PID) pathways were downloaded from v7.1 of the Molecular Signatures Database.

TOPMed imputation-based age-at-onset GWAS

In our most recently reported HD age-at-onset GWAS, we analyzed 9,058 individuals of European ancestry whose genome-wide SNV data was imputed from array-typed genotypes using HRC data.³ This GWAS identified nine genome-wide significant loci tagged by SNVs of minor allele frequency (MAF) > 1%. Six of these loci harbor DNA maintenance genes (FAN1, MLH1, MSH3, PMS1, PMS2, LIG1) and, for three of these (FAN1, MSH3, LIG1), we further defined multiple haplotypes representing distinct modifier effects.^{3; 35} The other three loci harbor candidate genes (RRM2B, TCERG1, CCDC82) that might affect somatic instability indirectly or, alternatively, might influence subsequent toxicity/response mechanisms. We also reported SNVs tagging two infrequent (< 1%) HTT haplotypes that differ from canonical CAG expanded alleles. In canonical alleles, the CAG repeat is followed by a single CAACAG codon doublet but in the variant alleles the CAA interruption is either lost or the CAACAG doublet is duplicated.^{3; 36-38} These variants cause the length of the uninterrupted CAG repeat to be estimated incorrectly in a standard PCR-based fragment-sizing assay using canonical sequence standards. The initially mis-estimated CAG length resulted in incorrect assignment of expected age-at-onset, which then produced false association signals. After correction of the mis-estimated CAG sizes by direct sequencing in those subjects identified by SNVs tagging the false signal, the CAACAG duplication haplotype was no longer significant while the CAA loss haplotype continued to yield signal at the border of genome-wide significance (p-value 5E-08). Finally, we also reported several additional loci detected only by single SNVs at < 1% MAF but suggested that these might be spurious due to outlier effects and required further confirmation.

For the present study, we re-imputed genotypes genome-wide using TOPMed data²³ as a reference set and again performed GWAS of residual age-at-onset as a continuous phenotype for those 9,009 individuals whose data passed quality control. In preliminary analyses, a series of variants flanking the HTT CAG repeat location identified additional individuals carrying the rare CAA-loss and CAACAG-duplication haplotypes who had been missed by HRC imputation. We sequenced their CAG repeat region to correct the CAG repeat lengths and corresponding age-at-onset residuals prior to repeating a GWAS that yielded results largely comparable to HRC imputation for the 9,009 individuals (Figure S2). However, in the relationship between CAG repeat length and age-at-onset, the range of residual age-at-onset values is greater at lower than at higher repeat length (Figure S3). Consequently, to allow direct comparison with the Z score-based algorithmically-predicted landmarks (below), we calculated Z-scores for age-at-onset at each individual CAG repeat length for use as the continuous phenotype in GWAS.

The TOPMed-based Z-score age-at-onset GWAS (Figure 2) detected all of the loci previously implicated at genome-wide significance by SNVs of > 1% MAF except the chromosome (chr) 11 CCDC82 locus, which in our previous GWAS was detected at genome-wide significance only by a dichotomous analysis of allele frequencies in the 30% extremes of residual age-at-onset.³ At MSH3, FAN1, and LIG1, where we have reported multiple modifier effects, SNVs tagging each of these effects again showed at least suggestive evidence of association (p-value < 1E-05) with age-at-onset (Table S1). With one exception, those peak SNVs in the current analysis that were not identical to the previously reported tag SNVs showed similar allele frequency, direction of effect and strong allele correlation (r-square > 0.97) with the prior tag.³ The peak SNV associated with onset-delay at the PMS2 locus showed a slightly lower correlation with the previously assigned 7AM1 tag SNV (r-square = 0.86). In addition, second modifier effects (here termed 2AM2 and 7AM2) were detected in the vicinity of PMS1 (p-value 3.67E-06 at 2_190639862_C_T) and PMS2 (p-value 2.44E-07 at 7_6056484_G_C) based upon suggestive significance of SNVs associated with onset-delaying and onset-hastening modification, respectively, opposite in direction to the 2AM1 and 7AM1 modifier effects. Notably, after the correction for uninterrupted CAG repeat length, the rare CAA-loss and CAACAG-duplication HTT-region haplotypes continued to support genome-wide and suggestive significant association, respectively (Figure S4, Table S1). This suggests either that these changes to the CAG repeat region sequence alter its inherent propensity for somatic instability, changing the rate of CAG expansion leading to onset, that they influence the later toxicity mechanism, or that there is a separate linked modifier variant in 4p16.3 on one or both of these extended haplotypes.^{3; 38} In addition to previously recognized loci, a novel genome-wide significant onset-delaying signal emerged on chr 7 (7BM1: p-value 7.29E-09 at 7_56241504_C_T; Figure 2; Table S1) near the testis-expressed NUPR2 (nuclear protein 2, transcriptional regulator), and upstream of both CHCHD2 (coiled-coil-helix-coiled-coil-helix domain containing 2), which encodes a negative regulator of apoptosis involved in response to mitochondrial stress, and PHKG1 (phosphorylase kinase catalytic subunit gamma 1), which regulates glycogenolysis. Finally, genome-wide analysis of SNVs with very low MAF (> 0.1% but < 1%) failed to reproduce any of the ‘hits’ represented by single rare SNVs reported previously from HRC imputation,³ reinforcing the need for caution in interpreting such loci.

• Download figure
• Open in new tab

Figure 2 GWAS of age-at-onset, age-at-DCL4 and age-at-TFC6

Stacked plots comparing GWAS are shown for the Z score of CAG-specific age-at-onset, the Z score from age-at-DCL4 prediction and the Z score from age-at-TFC6 prediction. Solid arrows indicate loci that yield genome-wide significant signal with SNV of > 1% MAF (with the exception of HTT), with candidate modifier genes all labeled in the top panel. One novel locus that awaits further confirmation is noted by a dashed arrow and its HD modifier effect designation (7BM1) based upon our standard nomenclature.³ Red labels indicate DNA maintenance genes.

Algorithmic prediction and GWAS of additional HD landmarks

To perform GWAS for additional landmarks during HD pathogenesis, we explored the use of predictive modeling, as illustrated in Figure S1, to produce estimates of the ages at which each individual had a 50% probability of reaching DCL4 and TFC6, respectively. Briefly, we generated a survival curve for the landmark phenotype (either DCL4 or TFC6) for each combination of CAG repeat length, sex, and education level (low = high school or less, or high = more than high school) to determine the median age at survival (i.e., the age at which 50% of the individuals have reached the landmark phenotype). We also modeled the mean trajectory for the cUHDRS to determine the cUHDRS value at the aforementioned median survival age. Individual cUHDRS trajectories extrapolated over the entire age range then predicted the individual’s age that corresponded which this predetermined cUHDRS value. The individual age-at-landmark prediction was then transformed to a Z-score by correcting for the median age and scaling by the standard deviation (SD) of all the ages for the individuals in the same CAG, sex and education cohort.

The landmark predictions were made using longitudinal clinical data from 19,173 participants in the Predict-HD, COHORT, Registry and Enroll-HD studies. Genome-wide genotyping data are currently available for ∼7,400 of these participants, ∼5,200 of which could be analyzed in the age-at-onset GWAS. The predicted age-at-DCL4 and age-at-TFC6 phenotype Z-scores were highly correlated with each other and both were less correlated with the age-at-onset Z-score (Table 1). We tested the capacity of these phenotypes to capture modifier loci by age-at-landmark GWAS directly comparable to the age-at-onset Z-score analysis. As shown in Figure 2 and in Table S1, despite the smaller sample size for both age-at-DCL4 (N = 6,697) and age-at-TFC6 (N = 6,900) landmarks compared to age-at-onset (N = 9,009), the predicted phenotypes revealed genome-wide significant signals at all but one of the previously reported loci.

Two intriguing aspects of the age-at-landmark GWAS were the striking differences compared to the age-at-onset GWAS in the level of significance at some loci and the changes in their relative significance compared to other loci in the same analysis. The modifier signals for some loci, such as MSH3, PMS2, CCDC82 and FAN1 were increased substantially, while others, such as PMS1, TCERG1, RRM2B and LIG1 were not. For most loci, the source of the signal could be equated confidently to the same modifier effect captured by age-at-onset because the peak signal was attributed to the same SNV or to one very highly correlated with it (r-square > 0.97), although exceptions at MSH3, PMS2 and LIG1 hinted at a potentially greater complexity underlying these modifier effects (Table S1).

The change in relative signal was most dramatic at MSH3. In the age-at-onset GWAS, the common 15AM2 delaying modifier at FAN1 was far more significant (peak p-value 7.23E-28 at 15_31241346_G_A) than the common 5AM1 hastening modifier at MSH3 (peak p-value 1.05E-11 at 5_79913275_G_A). By contrast, while both 5AM1 and 15AM2 showed greater significance in the age-at-DCL4 and age-at-TFC6 analyses than the maximum achieved in the larger age-at-onset sample, the signal for 5AM1 was disproportionately higher (peak p-values 2.69E-40 and 1.42E-39, respectively, at 5_79950781_A_G), becoming comparable to 15AM2 (peak p-values 9.03E-38 and 3.45E-41, respectively, at 15_31241346_G_A). This change in absolute and relative signals was not due to differences in the sample size or in the individuals tested, since analysis of only those 5,009 individuals where all three phenotypes were available generated similar results (5AM1: peak p-values 9.18E-07 for age-at-onset, 6.26E-32 for age-at-DCL4 and 1.06E-33 for age-at-TFC6, at 5_79950781_A_G; 15AM2: peak p-values 8.55E-14 for age-at-onset, 5.33E-33 for age-at-DCL4, 8.48E-34 for age-at-TFC6, at 15_31241346_G_A).

Relative p-value improvement in GWAS of the cUHDRS-predicted landmarks was also particularly notable at PMS2 (Figure 2), again indicating that the age-at-onset and predicted landmark phenotypes had a differential capacity to reveal effects of particular modifiers. In this instance, the top SNV associated with landmark delay (7_6041836_T_A) had a much higher MAF than the 7AM1 age-at-onset modifier tag (MAF ∼41% vs. ∼15%), pointing to an additional modifier effect and thereby prompting the naming of a third modifier effect (7AM3) at the locus (Table S1).

The relative p-value improvement with the algorithmically predicted landmarks was also evident beyond the DNA maintenance genes. Both age-at-DCL4 and age-at-TFC6 robustly detected the CCDC82 locus at far greater significance than age-at-onset. The HTT CAA-loss and CAACAG-duplication haplotypes also showed increased significance in the predicted-landmark GWAS. Surprisingly, the predicted landmarks yielded only suggestive significant signal at RRM2B, a locus that emerged as genome-wide significant in our first HD GWAS² and was confirmed in our most recent larger GWAS³ (as well as in the age-at-onset analysis in Figure 2). Overall, the age-at-DCL4 and age-at-TFC6 analyses appeared to track a subset of established modifier loci and did not uncover evidence for a frequent modifier that acts after the DCL4 landmark by a new mechanism. However, along with the age-at-onset GWAS, the age-at-DCL4 and age-at-TFC6 GWAS did nominate several new loci represented by alleles at MAF < 1% (Table S2) that, like the new chr7 locus seen only in the age-at-onset analysis, will require future confirmation as bona fide modifiers since each could result from a small number of phenotypic outliers.

Dichotomous analysis of TMS and SDMT extremes

The selectively greater capacity of the predicted landmarks to detect particular loci, such as MSH3 and PMS2, relative to other known age-at-onset modifiers pointed to a difference in the nature of these phenotypes compared to the observed age-at-onset phenotype. Age-at-onset was typically determined by expert raters based upon motor signs characteristic of HD whereas the algorithmic prediction of age-at-DCL4 and age-at-TFC6 made use of motor, cognitive and functional measures in the cUHDRS in refining the age-at-landmark predictions. Furthermore, the age-at-onset analysis was limited to individuals who had the event of onset, whereas the prediction-based analysis had no such constraint. Consequently, we reasoned that, as opposed to landmarks that reflect a combination of motor and cognitive impairment, a direct comparison of two measures with no constraining event, the TMS and SDMT, aimed primarily (though not exclusively) at motor and cognitive domains, respectively, might reveal a differential impact of the modifiers on these domains.

In an initial comparison, we examined the longitudinal records of the Registry study of the European Huntington’s Disease Network and of the Enroll-HD study and identified those individuals who, for a given CAG repeat length, consistently fell within the extremes of the TMS or SDMT score distributions. We then performed a dichotomous GWAS comparing SNV allele frequencies between those extreme high and low scoring individuals for whom genome-wide genotyping data were available (for TMS (N = 920) and for SDMT (N = 1,200)) (Figure 3). The FAN1 locus showed genome-wide significance for both phenotypes with peak-scoring SNVs being associated with lower TMS and higher SDMT score (top SNV p-value TMS: 4.54E-17 at 15_31241346_G_A; SDMT: 1.23E-12 at 15_31279883_CTTTT_C) indicative of less disruption of both motor and cognitive domains. By contrast, the MSH3 locus failed to show even suggestive significance for the TMS phenotype but yielded two distinct peak signals in the SDMT extreme analysis for SNV’s whose minor alleles were associated with better and worse performance, respectively (top SNV p-values 6.64E-16 at 5_79915020_G_A and 8.54E-15 at 5_79947431_A_G). No other locus showed genome-wide significance in either of these small, highly selected TMS-and SDMT-extreme samples.

• Download figure
• Open in new tab

Figure 3 Dichotomous GWAS of Total Motor Score and Symbol Digit Modalities Test extremes

Stacked plots of dichotomous GWAS comparing allele frequencies among HD subjects in the Registry study who represented consistent phenotypic extremes with respect to TMS (top) or SDMT (bottom). Arrows mark the two loci that yielded genome-wide significant signal in at least one of the analyses with modifier gene labels denoted in the top panel.

Analysis of TMS and SDMT score landmarks

To address the differential action of the modifier loci on TMS and SDMT in a manner directly comparable to our age-at-landmark predictions of DCL4 and TFC6, we then used algorithmic prediction similar to Figure S1 panels 3 and 4 (using linear mixed models for SDMT and TMS rather than cUHDRS) to compute estimated ages in early manifest HD at which each individual progressed up to TMS = 30 (N = 6,897) and down to SDMT = 30 (N = 6,727). The resulting age-at-TMS30 and age-at-SDMT30 landmark Z-scores were more highly correlated with the age-at-DCL4 and age-at-TFC6 Z scores than with each other or with the Z-score for age-at-onset (Table 1).

The GWAS for these two new landmark phenotypes (Figure 4 and Table S1) reinforced the capacity of distinct disease phenotypes to detect individual modifier loci differentially. For example, the MSH3 locus showed relatively greater significance than FAN1 for age-at-SDMT30. The relative differences in modifier peaks for age-at-TMS30 and age-at-SDMT30 were not due to the slightly different samples tested, since p-values calculated based upon only the 6,716 participants who were scored for both age-at-TMS30 and age-at-SDMT30 did not substantially alter the results (e.g., joint set vs. full sample p-values for 5_79950781_A_G (5AM1): age-at-TMS30 [3.62E-19 vs 1.35E-23] and age-at-SDMT30 [4.96E-37 vs 4.70E-37]; for 15_31241346_G_A (15AM2): age-at-TMS30 [2.51E-41 vs 7.41E-45] and age-at-SDMT30 [7.19E-22 vs 1.34E-21]. As expected, the different information content of these landmark phenotypes was also evident in conditional analyses, where the age-at-TFC6 signal at MSH3 was preferentially reduced by conditioning on age-at-SDMT30 whereas at FAN1, it was preferentially reduced by conditioning on age-at-TMS30 (Figure 5). Taken together with the SDMT extreme analysis, these results confirmed that the integration of SDMT into the cUHDRS, which was used in the prediction of age-at-DCL4 and age-at-TFC6, at least partially explains the relatively increased MSH3 signal in those GWAS compared to that for age-at-onset.

• Download figure
• Open in new tab

Figure 4 GWAS of Total Motor Score and Symbol Digit Modalities Test landmarks

Stacked plots comparing GWAS analysis of age-at-TMS30 and age-at-SDMT30. Arrows indicate loci that yielded genome-wide significant signal for variants of >1% MAF (with the exception of HTT), with the candidate genes labeled in the top panel. Red labels indicate DNA maintenance genes.

• Download figure
• Open in new tab

Figure 5 Age-at-TFC6 signal at MSH3 and FAN1 conditioned on other phenotypes

Regional plots displaying age-at-TFC6 association significance (-log10(p-value)) at MSH3 (left, chr 5) and FAN1 (right, chr 15) loci for the 4,879 participants for whom age-at-onset, age-at-DCL4, age-at-TFC6, age-at-TMS30 and age-at-SDMT30 were all available are shown above equivalent plots where the analysis was conditioned on each of the other age-at-landmark phenotypes. Each circle represents a different SNV with colored symbols indicating tag SNVs for particular modifier effects derived from the analysis of all 6,900 participants with age-at-TFC6 data: 5AM1, 5_79950781_A_G (red), 5AM2, 5_80086504_A_G (green); 5AM3, 5_79955360_G_GT (purple, peak SNV after conditioning on 5AM1 and 5AM2); 15AM1, 15_31202961_G_A (red), 15AM2, 15_31241346_G_A (green); 15AM3, 15_31197995_C_T (purple); 15AM5, 15_31282611_T_C (cyan).

The age-at-TMS30 and age-at-SDMT30 analyses also supported additional complexity of modifier effects at some loci where the top SNVs were not highly correlated across phenotypes. At MSH3, the common landmark-hastening 5AM1 modifier effect was captured by highly correlated peak SNVs across all four predicted landmarks (r-square > 0.97) while the rare 5AM2 delaying effect was revealed by the same tag SNV for all phenotypes except age-at-SDMT30, where it did not achieve suggestive significance (Table S1). However, the top common SNVs associated with landmark delay differed for each of the four predicted landmarks, with pairwise r-square values from 0.44 to 0.99, pointing to potential differences in the sources of the modifier effects. Consequently, we first performed association analyses conditioning simultaneously on the 5AM1 and 5AM2 peak SNVs to identify the most significant remaining modification signal for each phenotype. Each landmark yielded a different top SNV tagging a common landmark-delaying effect (Table 2) but only those for age-at-DCL4 and age-at-TFC6 were very highly correlated. We compared separate association analyses across all four predicted landmark phenotypes conditioning each time simultaneously on the 5AM1 and 5AM2 peak SNVs plus one of the top remaining SNVs defined in the previous step. Each of these conditional analyses showed greatly reduced signal for all phenotypes, indicating that the common causative landmark-delaying variants are not independent (Table 2). However, weak to suggestive significant SNVs remained, suggesting the possibility that relative strength of the effect on motor and cognitive circuits differs based on multiple variants in the local haplotype.

Of the DNA maintenance genes other than MSH3, age-at-SDMT30 detected only PMS2 and FAN1 at genome-wide significance. However, at PMS2, age-at-SDMT30 also revealed a second landmark-hastening modifier effect (7AM4), tagged at suggestive significance by a missense SNV (7_6026530_C_T; p.M622I) associated with altered PMS2-MLH1 binding.³⁹ Conditioning on either of the PMS2 landmark-hastening SNVs resulted in increased significance for the other, confirming that they represent different modifier effects (Table S3).

The detection of MLH1, TCERG1, LIG1 and RRM2B at genome-wide significance by age-at-TMS30 but not by age-at-SDMT30 suggests that the respective modifier alleles act preferentially on the motor domain, explaining their capture in the age-at-onset GWAS. However, these findings also suggest that TMS as a measure of motor dysfunction does not provide information identical to the expert rater’s subjective assessment of motor function when determining age-at-onset. In particular, the CCDC82 locus was robustly detected in continuous analysis by the age-at-TMS30 landmark, yet in the continuous analysis of age-at-onset it achieved only suggestive significance (Table S1).

Notably, based on pure, uninterrupted CAG length, the HTT CAA-loss haplotype showed genome-wide significance for both age-at-TMS30 and age-at-SDMT30, consistent with that variation at the site of the CAG repeat influencing both motor and cognitive domains. In contrast, the HTT CAACAG-duplication haplotype failed to achieve suggestive significance for either phenotype. The age-at-TMS30 and, especially, the age-at-SDMT30 GWAS also revealed several loci with apparent genome-wide significant association detected only by infrequent alleles (MAF < 1%) (Table S2). However, only one of these alleles achieved even suggestive significance in the age-at-onset, age-at-DCL4 or age-at-TFC6 GWAS suggesting that most or all are likely to result from the coincidence of rare alleles in a small number of phenotypic outliers.

Gene-wide association and pathway analyses

We also performed gene-wide association analysis across the landmark phenotypes (Table S4). As expected from the single SNV analyses, the genes at several modifier loci, particularly the MSH3 region, showed increased signal for some or all predicted landmarks compared to age-at-onset. Notably, the age-at-SDMT30 results showed the greatest divergence from age-at-onset in the gene-wide analysis, being far more significant at the MSH3 locus but much less significant at MLH1, LIG1, RRM2B, and TCERG1. Pathway analyses based upon these data again implicated mismatch repair, as was the case in our HRC-imputation GWAS (Table S5). Mismatch repair pathways topped the “self-contained” association list for all measures, with significance being greatest for age-at-TFC6. Mismatch repair terms were also most significant (p-value < E-5) in “competitive” analysis of age-at-TFC6 results and nominally significant with the other measures. Also, nucleotide excision repair topped the list in the competitive age-at-SDMT30 analysis, driven largely by FAN1, LIG1, and POLD1, suggesting that genome-wide significant SNVs may emerge at the latter locus in future expanded GWAS analyses. Interestingly, terms related to regulation autophagy (GO:0016239 “positive regulation of macroautophagy” and GO:0010508 “positive regulation of autophagy”) were also significant in the age-at-TFC6 competitive analysis. Although HTT was the most significant individual gene in those gene sets, its exclusion did not dramatically reduce the enrichment p-values (e.g., age-at-TFC6 competitive / self-contained enrichments p-values with HTT: GO:0016239 1.40E-07 / 7.64E-07, GO:0010508 4.00E-06 / 3.50E-08; excluding HTT: GO:0016239 3.76E-07 / 9.05E-06, GO:0010508 8.44E-06 / 2.8E-07) arguing that the pathway enrichment is due to a number of modestly significant genes rather than a few highly significant genes.