|Home | About | Journals | Submit | Contact Us | Français|
Neural cultures derived from Huntington’s disease (HD) patient-derived induced pluripotent stem cells were used for ‘omics’ analyses to identify mechanisms underlying neurodegeneration. RNA-seq analysis identified genes in glutamate and GABA signaling, axonal guidance and calcium influx whose expression was decreased in HD cultures. One-third of gene changes were in pathways regulating neuronal development and maturation. When mapped to stages of mouse striatal development, the profiles aligned with earlier embryonic stages of neuronal differentiation. We observed a strong correlation between HD-related histone marks, gene expression and unique peak profiles associated with dysregulated genes, suggesting a coordinated epigenetic program. Treatment with isoxazole-9, which targets key dysregulated pathways, led to amelioration of expanded polyglutamine repeat-associated phenotypes in neural cells and of cognitive impairment and synaptic pathology in HD model R6/2 mice. These data suggest that mutant huntingtin impairs neurodevelopmental pathways that could disrupt synaptic homeostasis and increase vulnerability to the pathologic consequence of expanded polyglutamine repeats over time.
HD is characterized by motor abnormalities, psychiatric symptoms and cognitive deficits. It is caused by a CAG repeat expansion in the huntingtin gene (HTT) encoding a polyglutamine tract expansion in the huntingtin protein (HTT)1. CAG repeats of 40 or more cause symptoms and those above ~60 cause juvenile-onset HD. Mutant HTT (mHTT) is implicated in a spectrum of cellular aberrations2. Neuropathology includes cortical atrophy and loss of striatal medium spiny neurons1. Age of clinical onset is primarily predicted by the extent of CAG repeat expansion and considered to result from cumulative toxic insults, together with environmental and other genetic factors3.
There may be deficiencies of neurodevelopment and differentiation of neural stem and progenitor cells in the adult striatum throughout life4,5. Adult neurogenesis appears to be impaired in the striata of HD patients, with increased cell proliferation6 and an absence of adult-born neurons7, suggesting that neurogenesis may be initiated but maturation is impaired. Neuroimaging scans of premanifest HD-affected brains detect changes in striatal, cortical and whole-brain volume before symptoms8–10. Head circumference is less in premanifest HD subjects than unaffected subjects11.
Patient-derived induced pluripotent stem cells (iPSCs), somatic cells reprogrammed to a pluripotent state, provide an important resource for deciphering mechanisms underlying neurological disease (for example, ref. 12) and aid in the design of new therapeutics. iPSCs can be differentiated into multipotent neural stem progenitor cells that produce mature neural subpopulations (for example, ref. 13). We previously reported that differentiated HD-derived iPSCs display expanded CAG-associated phenotypes14.
We used neural cells differentiated from HD patient-derived iPSC lines with juvenile-onset CAG repeat expansions (60 and 109 repeats) to explore the effects of pathological HTT alleles by using parallel cultures for unbiased omics analyses. This discovery-based approach revealed consistent deficits related to neurodevelopment and adult neurogenesis, suggesting that specific gene networks represent potential therapeutic interventions. We tested a small molecule, isoxazole-9 (Isx-9), that targets several of the dysregulated gene networks. Isx-9 normalized CAG repeat-associated phenotypes in both juvenile-repeat and adult-onset repeat HD iPSC-derived neural cultures, as well as cognition and synaptic pathology in R6/2 mice, demonstrating that these deficits might be reversed and synaptic homeostasis improved.
iPSC lines from non-disease and HD fibroblasts with CAG repeats from 21 to 33 (non-disease) and 60 and 109 (juvenile-onset HD range) were reprogrammed using non-integrating episomal factors (Supplementary Table 1)15 and pluripotency and normal karyotypes confirmed. iPSC lines were differentiated for 56 d into mixed neural cultures containing neurons, glia and neural progenitors as described14. Mixed cultures broadly represent cell types and progenitors in human brain, as opposed to pure cultures of a single cell type. Lines similarly expressed germ-lineage and cell type–specific markers (Supplementary Figs. 1 and 2). Cell quantification in HD and non-disease lines was performed (Supplementary Fig. 1a–e), with total counted cells expressing glial fibrillary acidic protein (GFAP, 16.1%), neuronal markers TUJ1 (26.5%) and MAP2ab (14.8%), and striatal marker DARPP-32 (14.3%) (Supplementary Fig. 1b–e), as an average among lines. Oligodendrocyte (Supplementary Fig. 1f), endoderm (SOX17, FOXA2), mesoderm (MYO1), and microglia (IBA1) (Supplementary Fig. 2a, b) markers were absent. iPSC lines had a similar overall composition from staining data, including nestin (Supplementary Fig. 1a). At earlier times, nestin-expressing neural progenitor cells may persist longer in the expanded 109Q repeat line, and nestin-positive cells were more susceptible to cell stress after acute brain-derived neurotrophic factor (BDNF) withdrawal15, suggesting early and subtle effects exerted by expanded CAG repeats.
We used unbiased whole-genome and multi-platform approaches in parallel iPSC cultures. Whole-transcriptome analysis (RNA-seq) was performed, and principal component analysis showed separation of HD 109Q and 60Q lines (two clones each) from non-disease 33Q, 28Q and 21Q lines, indicating minimal variability within groups and a maximal variance explained by disease and non-disease differences (Fig. 1a). Statistical analysis using the Bioconductor package DESeq (RNA-seq differential expression) revealed 1,869 differentially expressed genes (DEGs) between HD and non-disease lines (Supplementary Table 2). Hierarchical clustering (Fig. 1b) of log2-transformed gene expression values showed that samples clustered by repeat length, with distinct separation and expression patterning among the samples. Independent clonal lines derived from individual subjects clustered tightly together, with low variability within groups (Fig. 1b).
Ingenuity pathway analysis (IPA; Supplementary Fig. 2c) was used to investigate biological changes by examining genes from DESeq analysis. Of the 1,869 DEGs, 543 (29%) (Supplementary Table 2) centered on development. The top three were cellular development, nervous system development and function and tissue development (Fig. 2a). Over half of the DEGs (59%) were associated with nervous system development and function.
Some of these neural developmental genes and their physical and regulatory interactions are depicted in Figure 2b (Supplementary Table 3), with NEUROD1 as a highly enriched hub. HTT interacts physically with the products of two genes that are among the most downregulated in the RNA-seq data set: NEUROD1 and GAD1 (refs. 16,17). Expression of the proneuronal basic helix-loop-helix gene NEUROD1 was decreased across all HD lines (Supplementary Table 2). NEUROD1 regulates neurodevelopment18 and adult neurogenesis19. Of genes in this network (Fig. 2b and Supplementary Table 3), several encode proteins that regulate NEUROD1 expression (POU4F, NEUROG2, ASCL and REST) and are dysregulated or predicted to have aberrant activities at the protein level. RNA expression levels for several, including NEUROD1, were validated by quantitative PCR (qPCR) (Supplementary Fig. 3a–f).
IPA also predicted the activation or inhibition of upstream regulators (Supplementary Table 4), including several implicated in HD, such as increased REST activation and decreased BDNF signaling20. REST-mediated epigenetic remodeling regulates the developmental switch of synaptic NMDA receptors from GluN2B to GluN2A subunits, thus decreasing GRIN2B levels21; GRIN2B was downregulated in the RNA-seq data set. Several genes essential for neurogenesis, including NEUROD1, NEUROG2 and ASCL1 and the microRNA miR-124, are predicted to have inhibited activity. Genes encoding proteins involved in the transforming growth factor (TGF)-β pathway (TGFB2, TGFB3, TGFB3R) were upregulated in the differentiated HD iPSCs, with predicted activation of upstream regulators TGFβ1 and TGFβ1R. TGFβ signaling is implicated as a ‘switching factor’ whose levels and temporal activity require tight regulation during development to determine neuronal cell fate22. This signaling network was disrupted in neural stem cells from HD patient iPSCs23. Highlighted in this network were genes overlapping the IPA-predicted regulators and genes that are major hubs with direct connections to other genes in the larger network, suggesting an association with HD.
IPA also identified dysregulated pathways relevant to neuronal development and maturation, including axonal guidance, Wnt signaling, Ca2+ signaling, neuronal CREB signaling, and glutamate and GABA receptor signaling (Fig. 3a), that are altered in HD models or human HD2. Axonal guidance pathways regulate connectivity, integrating actions of guidance cues and receptors. Four families of molecules and receptors provide axon guidance cues24, including netrins, slits, ephrins and semaphorins, and genes within each family are dysregulated in the RNA-seq data set (Supplementary Table 2), with most downregulated.
Intracellular calcium signaling promotes axonal and dendritic outgrowth, connecting alterations in axonal guidance and calcium-signaling pathways, and is critical for neuronal synaptic activity–regulated transcription and maturation25. Widespread dysregulation of Ca2+ signaling pathways in the HD samples (Fig. 3b) included genes for members of the NMDA- and AMPA-type glutamate receptors, nicotinic acetylcholine receptor subunits, several subunits of the voltage-gated Ca2+ channel CACNA1, and the plasma membrane Ca2+-ATPase, in addition to the calcium sensors and downstream effectors CAMKII, CALM and CREB (Fig. 3b). Among downregulated DEGs encoding proteins in the glutamate and GABA signaling pathways are glutamate decarboxylases and GAD1 and GAD2, which encode GABAergic neuron markers. Several other DEGs in the HD iPSC-neural cultures are involved with GABA synthesis, release, reuptake or degradation. As glutamate is transported into neurons via SLC1A3 and SLC1A6, downregulating the corresponding genes would be predicted to yield reduced glutamate substrate for GABA synthesis. This is significant: GABAergic neurons are the most vulnerable to mHTT toxicity1. The dysregulated genes display integrated features and potential effects on nervous system development and function, including clustering of DEGs encoding proteins involved in CREB, Wnt, axonal guidance signaling and synaptic function (Supplementary Fig. 3g, h).
We next compared the 1,869 DEGs from the RNA-seq analysis to orthologs (1,647 mouse genes) involved in mouse striatal development derived from microarray mouse developmental data using identifiers from Xspecies, a cross-species ortholog identifier and retrieval tool. This comparison yielded 679 overlapping genes. Hierarchical clustering showed that the HD60Q and HD109Q RNA-seq samples clustered with the proliferative germinal zone samples whereas non-disease samples clustered with the postmitotic mantle zone samples (Fig. 4a). Thus, differentiated HD lines may mature slower than non-disease lines or improperly shut off gene expression in earlier stages of development.
Enriched GO terms for the set of 679 genes (Fig. 4b) included synaptic transmission (GO:0007268), nervous system development (GO:0007399), CNS development (GO:0007417), and a series of more specific processes, such as axonal guidance (GO:0007411), cell adhesion and locomotor behavior (GO:0007626); these were similar to the categories enriched in the RNA-seq analysis. Using the hierarchical cluster algorithm with average linkage for gene expression values, we obtained 13 clusters from the 679 genes (Supplementary Fig. 4). Activation z-scores predicted that most of the processes have decreased or inhibited function in the HD lines (Supplementary Fig. 5a). Analysis with Metacore to visualize functional connections between genes in the sets yielded 296 nodes and also pointed to developmentally described transcription factors (DLX2, c-MYC, DACH1, MEF2A, hASH1) and relevant pathway regulators (TGFβ receptor) (Fig. 4c). Gene cluster I (Supplementary Fig. 5b) contains NEUROD1 and related genes that are expressed more highly in non-disease repeat samples.
Finally, we compared HD DEGs to genes that change in expression during human striatal maturation26. Genes implicated in development of human striatum are regulators of genes differentially expressed in our HD neural cells (Supplementary Fig. 5c), indicating a core network of genes that both contribute to the maturation of the human striatum and are altered in HD iPSC-derived neural cells.
Widespread transcriptional and epigenetic changes suggest that epigenetic dysregulation represents a primary pathognomonic molecular feature of HD27, and chromatin remodeling has been implicated in neural development and synaptic plasticity28. We performed chromatin immunoprecipitation and sequencing (ChIP-seq) for three histone marks: H3K4me3 (promoter-specific activation), H3K27ac (enhancer) and H3K36me3 (in actively transcribed gene bodies). In HD and non-disease lines, most H3K4me3 peaks were in promoter and intronic regions and most H3K27ac and H3K36me3 peaks were in intergenic regions. The number of peaks did not seem to be related to mHTT repeat length.
We reported an epigenetic pattern in wild-type mice that marks genes with altered expression in R6/2 mice, which express an expanded repeat containing HTT fragment29. Applying a similar analysis to the iPSC data, we found that, in agreement with results in mouse tissue, genes with altered expression in high-repeat-number HD lines were more likely to be members of a particular cluster (class 1) for H3K4me3 (P = 3.7 × 10−6 for genes with higher expression in HD, P = 1.1 × 10−36 for lower expression). This profile was characterized by a broad peak just downstream of the TSS in non-disease cell lines (Fig. 5a). We observed a similar pattern for H3K27ac (P = 6.3 × 10−25 for higher expression, P = 3.5 × 10−45 for lower expression) (Fig. 5b). In both, GO enrichment analysis showed that class 1 genes were enriched for neuronal terms, such as neuron differentiation (FDR = 2.2 × 10−24 H3K4me3, FDR = 8.7 × 10−31 H3K27ac) (Supplementary Table 5a, b), supporting the idea that this profile marks a coherent group of genes relevant to the disease and dysregulated in HD lines. For H3K36me3, genes upregulated in high-repeat cell lines were likely to be members of a cluster (class 5) that in non-disease cell lines had a peak closest to the TSS of the gene (P = 5.7 × 10−7 for higher expression, P = 1.1 × 10−4 for lower expression) (Fig. 5c). GO analysis of H3K36me3 class 5 genes showed enrichment for neurological terms, such as synapse (FDR = 7.1 × 10−8) (Supplementary Table 5c). These profiles were observed in histone mark data from non-disease cell lines, marking genes that are affected in the HD lines. Therefore, the findings provide an important mechanistic clue by revealing the potential vulnerability of promoters with a specific signature that leads to transcriptional dysregulation in HD. These classes of genes and the enzymes that regulate their histone marks could therefore be selective targets for therapies attempting to restore transcriptional regulation to the normal state.
Changes in histone marks associated with polyglutamine expansion revealed widespread epigenetic changes near relevant genes. We compared pooled data from all expanded repeat HD cell lines and non-disease samples. Genes near H3K4me3 peaks with higher levels in HD cell lines were associated with cell-adhesion terms, including cadherin (FDR = 1.4 × 10−21) and calcium-dependent cell-cell adhesion (FDR = 6.7 × 10−5) (Supplementary Table 5d). Cell adhesion machineries mediate several neuronal functions, including neurite outgrowth and synaptogenesis30. Genes within 10 kb of H3K27ac peaks that were stronger in HD lines were associated with structural GO terms, such as actin cytoskeleton (FDR = 7.7 × 10−5) (Supplementary Table 5e), and genes associated with higher levels of H3K27ac in the non-disease lines were associated with neurological GO terms, such as synaptic transmission (FDR = 1.4 × 10−11) and neuron differentiation (FDR = 4.7 × 10−9) (Supplementary Table 5e). These results are similar to the RNA-seq results and suggest that disease-dependent modification of epigenetic components subtly affects cell lineage determination and perturbs neuronal fate specification. IPA analysis identified pathways and functional categories in common with RNA-seq results, including calcium signaling, GABA receptor signaling, axonal guidance and nervous system development and function.
We observed a highly significant correlation between the ChIP-seq signal and RNA-seq (Pearson P < 1 × 10−10 in all cases), suggesting a tight link between these two processes (Supplementary Fig. 6a–g). When each epigenetic class is individually considered, correlations between levels of ChIP-seq signal and RNA-seq are striking (Supplementary Fig. 6a–f). For example, H3K27Ac levels and transcription changes for class 1 genes in the 21Q cell line had a correlation coefficient of 0.7 (P < 1 × 10−10). In every cell line, class 1 genes in H3K27ac and class 5 genes in H3K36me3 showed the highest correlation between ChIP-seq counts and mRNA expression (Supplementary Fig. 6). Tracks for the gene ELAV3 provide an example of this correlation (Supplementary Fig. 6g).
We next looked for sequence motifs that were significantly enriched near those H3K27Ac peaks that were stronger in the HD lines or, separately, in non-disease lines. The top motifs are shown (Supplementary Fig. 6h).
GO enrichment analysis of transcription factor motifs with a P value <0.05 highlighted several differences. The primary DNA-binding regions used by these transcription factors appeared to differ, with motifs near non-disease-biased H3K27ac peaks enriched for zinc-finger (FDR = 7.6 × 10−24) and helix-loop-helix (FDR = 2.6 × 10−4) DNA-binding terms, but those under HD-biased peaks were enriched for the GO terms homeobox (FDR = 7.0 × 10−59), forkhead (FDR = 5.3 × 10−17) and leucine zipper (FDR = 3.2 × 10−11) (Supplementary Fig. 6i). Within HD and non-HD-biased H3K27ac peaks, categories of transcription factors were enriched for the term embryonic morphogenesis (HD FDR = 2.4 × 10−12, non-disease FDR = 2.2 × 10−5) and similar embryonic development terms. Interestingly, transcription factors whose motifs were found near HD-biased peaks were enriched for terms related to neuronal development, including forebrain (FDR = 4.0 × 10−5), diencephalon (FDR = 0.043) and regulation of nervous system development (FDR = 0.026), suggesting altered neurodevelopmental pathways in HD lines. Enrichment of certain motifs highlights pathways that may affect neuronal development. For instance, NEUROD is a member of a cluster of motifs enriched near peaks that are stronger in non-disease cell lines (P = 3.1 × 10−22) (Supplementary Fig. 6i), consistent with downregulation of the NEUROD1 transcript in HD lines. REST similarly shows enrichment in non-disease cell lines (P = 2.0 × 10−7), indicating that REST-regulated genes have lower H3K27ac in HD cell lines, consistent with REST being a predicted upstream regulator of differentially expressed genes (Supplementary Table 4). Thus, master regulators of nervous system development and function may be altered in HD lines.
To determine whether NEUROD1 repression depends on HTT in HD neural cultures, we evaluated whether an HTT-directed antisense oligonucleotide (ASO)15 would normalize expression in an expanded repeat line (109Q). HTT silencing in HD neural cultures increased NEUROD1 expression (Fig. 6a). We asked whether transcriptional alterations in the HD neural cells are modulated by enhanced expression of NEUROD1 in the same line. Lentiviral transduction of NEUROD1 cDNA into differentiated 109Q iPSCs increased expression of selected genes based on altered gene expression (Fig. 6b and Supplementary Fig. 7a). Genes upstream of NEUROD1, such as POU4F2, were not upregulated. In general, increases were not as significant in an unexpanded line (Fig. 6b). The extent to which an increase could be detected was limited by transduction efficiency, which was 30–40% depending on the line. Therefore, we evaluated whether modulation by a small molecule would de-repress transcriptionally repressed genes identified by RNA-seq.
Isx-9 upregulates genes, including NEUROD1, potentially through Ca2+ influx and other mechanisms, and it induces neuronal differentiation of cortical and subventricular zone (SVZ) progenitors31. We evaluated two sets of differentiated non-disease (21Q, 33Q) and HD (60Q, 109Q) iPSC lines to determine whether a neurogenic profile would be reestablished after treatment with 20 μM Isx-9, a dose that induces neurogenesis in a variety of cells, adult mouse whole brain and SVZ progenitors31. Expression of NEUROD1 (Fig. 6c and Supplementary Fig. 7b, c) and other selected genes in categories including neurodevelopment and neuronal function (POU4F2), Ca2+ signaling (CALB1, CAMK4, ATP2B2), GABA synthesis (GAD1) and potassium and calcium channels involved in neuronal excitability (KCNQ3 and CACNA1A, respectively) were evaluated (Supplementary Fig. 7b, c). qPCR showed that Isx-9 treatment activated transcription of these genes. Isx-9 treatment of differentiated HD and non-disease cells increased protein levels of NEUROD1 (Fig. 6d–e and Supplementary Fig. 8a), CALB1 and CAMK4 (Supplementary Figs. 7d and 8c, d). This was not due simply to differences in cell composition, as early differentiation cell markers tested (doublecortin and nestin) did not show differences in treated versus untreated (Supplementary Fig. 9a)—as expected, since treatment was begun at day 37, when cells are past the neuroblast stage.
To establish whether Isx-9 modulates mHTT-related phenotypes, we used a series of functional assays. BDNF withdrawal causes reduced cell viability of differentiated HD iPSC cultures in a CAG repeat-dependent manner14. Differentiated iPSCs were transferred to BDNF-free media, with or without 20 μM Isx-9, and compared to cultures with 20 ng/ml BDNF. BDNF withdrawal produced cell death in the HD (109Q) line (Fig. 7a), but not a non-disease (21Q) line (Fig. 7a and Supplementary Fig. 9b). Cell death rescue was nearly complete when differentiated HD iPSC cultures were treated with 20 μM Isx-9 during BDNF withdrawal (Fig. 7a). Notably, rescue depended on the presence of NEUROD1, as knockdown with a NEUROD1 short hairpin RNA prevented the Isx-9 effect (Fig. 7b).
We also wanted to test patient-derived HD neural cells with CAG repeats representing adult-onset HD. Using a similar differentiation method that yielded numerous MAP-2-positive cells with ~5–15% DARPP-32 and low numbers of proliferating cells, we examined survival of control and HD neural cells by robotic microscopy (Fig. 7c and Supplementary Figs. 9c–f and 10a–h)14. Robotic microscopy uses a custom-built automated platform to collect epifluorescence or confocal images from individual cells over their lifetime with high sensitivity to detect phenotypic differences in cells in the absence of stress. We subjected HD iPSC-derived neural cells with 46–48 CAG repeats (46Qn1 and 46Qn10) and 53–58 CAG repeats (53Qn3 and 53Qn5) to robotic microscopy. Survival analysis showed greater cumulative risk of death in HD cells than in 18Qn2 and 18Qn6 lines (Fig. 7c and Supplementary Fig. 9c–f). Due to a low number of samples (n = 2), the 18Qn6 clone was not significantly different from the Q46 clones by itself (Supplementary Fig. 9d); however, combining it with its sister clone, we detected survival differences (Fig. 7c). Adding 20 μM Isx-9 increased survival of all HD clonal lines (Fig. 7c and Supplementary Fig. 9c–f). Isx-9 also increased survival of the control line, but this effect appeared less effective and robust (Supplementary Fig. 9).
We also examined whether HD neural cells and controls differed morphologically. In human HD, mHTT elicits a progressive degenerative morphology in dendrites and branches or axons32. We examined the neurite-like processes in control and HD neural cells. Processes were measured on day ~39 of differentiation. For all HD lines, the neurite-like processes of 46Q, 53Q and 109Q lines were longer than for control 18Q and 28Q lines (Fig. 7d and Supplementary Fig. 10i–p). These results suggest that a component of mHTT pathology may be an alteration in the developing circuits, possibly due to loss of axonal guidance cues and signaling events that instruct connectivity, such as those identified in the RNA-seq analysis, which lead to an overproliferative phenotype. Remarkably, adding Isx-9 restored the neurite-like process length in the HD induced neurons to closer to that of controls, but did not affect the controls (Fig. 7d), suggesting Isx-9 corrects aberrant signaling and guidance cues specific to HD iPSC-derived neural cells.
RNA-seq analysis and functional annotation suggest impairment to genes involved in synaptic signaling and learning and memory (Supplementary Fig. 5a), and previous work suggests Isx-9 improves memory29 in mice33. Thus, we tested whether Isx-9 treatment could provide neuroprotection in R6/2 transgenic mice34. We first determined whether motor impairment in the mice was improved with Isx-9 treatment, but observed no rescue of rotarod, pole test or grip-strength deficits (Supplementary Fig. 11). However, when we evaluated recognition memory, Isx-9 treated R6/2 mice had better cognitive performance (exploratory preference score) than vehicle-treated counterparts (69.3% versus 54.62%) (Fig. 8a). Non-transgenic Isx-9 treated mice showed no improvement in this task, suggesting the treatment is selective for the mutant transgene context in R6/2 mice.
We then investigated whether Isx-9 reduces synaptic pathology in treated R6/2 mice, as synaptic loss has been reported in other HD mice (for example, ref. 35). Mice were sacrificed 5 weeks after initiating treatment, and the tissue was dissected for synaptic staining. Cortico-striatal synapses in the dorsal striatum were identified through colocalization of the presynaptic marker VGLUT1 and postsynaptic density protein HOMER1. R6/2 mice treated with vehicle alone showed fewer cortico-striatal synapses than nontransgenic littermate controls (Fig. 8b, c). Remarkably, R6/2 mice treated with Isx-9 had more synapses than those treated with vehicle alone, but numbers of synapses in non-transgenic mice were not affected, (Fig. 8b, c). These results demonstrate rescue of synaptic pathology after treatment with Isx-9.
Here we analyzed differentiated iPSCs from HD and unaffected subjects using unbiased multi-omics and bioinformatics. We identified alterations in genes and gene pathways that are associated with a pathological CAG repeat length. These iPSCs represent early stages of neuronal development, allowing elucidation of pathogenic mechanisms that affect disease onset or progression. Comparing DEGs and genes associated with changes in epigenetic motifs in differentiated HD iPSCs with highly expanded, juvenile-onset range CAG repeats suggests that alterations in an integrated epigenetic network affect neuronal development and could impair adult neurogenesis. Notably, Isx-9 provides neuroprotection to both juvenile and adult-onset range repeat lines, restores the structure of neuronal processes in HD cells and improves cognition and synaptic pathology in R6/2 mice, suggesting these pathogenic pathways may be targeted therapeutically.
Changes in gene expression determined by RNA-seq and comparison to mouse developmental transcription implicate altered neurodevelopment and maturation in HD pathophysiology. Neuronal and neurodevelopmental pathways are altered in differentiated HD neural cells, including those regulated by Wnt and NEUROD1, and affect expression of related genes, such as ASCL1 (Mash1), GAD1 and POU4F2. Levels of PSD2 and GAD1, which are among the most downregulated genes in differentiated HD lines, are essential for development of human striatum26.
Of particular relevance is the prediction that REST is the most activated upstream regulator of the differential gene expression in the HD cells. REST is implicated in neurotoxicity of mHTT20 and in dysregulation of genes, such as BDNF. It is also a master regulator of neurogenesis and neurodevelopment, potentially as a hub for recruiting epigenetic chromatin-modifying enzymes, and may promote neural stem or progenitor cell self-renewal but restrict generation and maturation of neurons (for example, ref. 36), consistent with reduced maturation in the HD neural iPSCs. Additional targets of REST include GABAA receptor genes. GABA receptors are among the earliest neurotransmitter systems to emerge during development, and this epigenetic regulation by REST may contribute to unbalanced generation of GABAergic and glutamatergic neurons, thus altering neurodevelopment. REST also regulates transition from neuronal progenitors to mature neurons, in part by regulating miR-124 (ref. 37), which is predicted to be inhibited in differentiated HD lines.
Our data support epigenetic ‘signatures’ that contribute to differential regulation of transcription in HD. Genes with a particular ‘shape’ or profile of H3K4me3 and H3K27ac peaks in the differentiated non-disease lines were enriched for DEGs in the RNA-seq analysis of HD lines. This tight association was reported in an HD mouse model and for selected genes in human autopsy samples29,38. We suggest an emerging theme that groups of genes marked by distinct epigenetic patterns associate with or drive specific properties, such as tissue-specific activities and developmental processes. A genome-wide analysis of H3K4me3 in neuronal nuclei from HD versus non-disease prefrontal human tissue resulted in enrichment for genes in neuronal development and neurodegeneration, notably HES4, which targets ASCL1 (Mash1)39. Thus, epigenetic patterns may represent critical mechanisms underlying pathogenesis and may identify enzymes and pathways for development of therapeutic approaches.
We observed decreased NEUROD1 expression in differentiated HD lines, which was alleviated by total HTT knockdown (HTT ASO) and deficits in other regulatory pathways (for example, TGFβ, ASCL1) whose central roles would likely affect embryonic and/or adult neurogenesis in HD. Impaired neurogenesis in HD and a role for NEUROD1 are supported by previous work describing aberrant hippocampal neurogenesis in R6/2 mice associated with reduced NeuroD1 (ref. 33). In human HD-affected SVZ overlying the caudate nucleus, a variety of proliferating progenitor stem cells display a heterogeneity of GABAA receptor subunits, supporting the presence of immature neural cells40. A series of reports (for example, ref. 6) highlighted an increase in SVZ progenitor cell populations in HD embryonic stem cells, mouse models and human tissue, suggesting upregulation of cues that simulate progenitor cell proliferation and neurogenesis. In the early stages of HD, postmortem brain sections show evidence for proliferative changes, such as increases in numbers and sizes of dendritic spines32 and in proliferation and size of the SVZ and number of newly born neurons6,41. Further, in developing murine ESCs, ESCs with expanded repeats have more elaborate and longer processes than ESCs with normal CAG repeat lengths42. This is consistent with our observation that neurites are in fact longer in adult onset HD neural cells. However, while adult neurogenesis appears to initiate within the striatum, these newborn neurons and interneurons are absent in the striatum of advanced-stage HD patients7, suggesting that maturation of neuronal precursors may not proceed appropriately.
We previously showed a selective early persistence of nestin-positive immature neural progenitor cells in differentiated juvenile onset HD iPSCs that were vulnerable to BDNF withdrawal–induced toxicity15. Rodent studies support the notion of deficits in striatal development in HD (for example, refs. 43,44). Even for adult-onset repeat ranges, autopsied brain tissue shows impairments that likely occurred during development and may arise from epigenetic factors39,45 eliciting dysregulation of developmental genes in human HD brain, such as homeobox genes. Impairments arising from an expanded polyglutamine expressed only during development46 may have consequences later in life, leading to phenotypes comparable to those in models with continuous mHTT expression. Finally, coexpression network analysis (WGCNA) in HD knock-in mice revealed that the module with the strongest association with CAG length and most extensive gene expression dysregulation contained striatal medium spiny neuron identity genes47.
Clinical symptoms take years to develop in adult-onset HD, possibly reflecting cognitive reserves and synaptic plasticity in the human brain. However, neurodevelopmental effects may disrupt homeostasis of the system, establishing a vulnerability to later effects of mutant HTT over time. While effects may not show up until adulthood, subclinical transcriptional alterations could create a different starting point in HD relative to non-HD brain. These pathways may contribute to aberrations in structural connectivity that influence clinical phenotypes48. Of note, children carrying the APOE ε4 risk allele for AD display altered brain development that may establish a neuroanatomical vulnerability to AD49.
Despite the complexity of the pathways identified using omics-based methods, the results identified an existing small molecule that could target the pathways affected. Isx-9 restored expression of selected dysregulated gene and protein expression patterns, and provided neuroprotection after BDNF withdrawal, the latter dependent on NEUROD1 gene expression. Isx-9 promotes Ca2+ influx in neural stem or progenitor cells in part through NMDA receptor signaling31. In the developing brain and during adult neurogenesis, electrical activity is essential for the differentiation, maturation and integration of neuronal circuitry (for example, ref. 50) This activity-dependent neurogenesis (termed excitation–neurogenesis coupling) is mediated by GABA-induced depolarization, Ca2+ influx and induction of neurogenic gene expression50. In this study, we observed HD mutation-associated changes in expression of genes encoding ion channels, GABA synthesis and transport, Ca2+ influx and basic helix-loop-helix proneural genes. The rescue of phenotypic deficits in cells for juvenile- and adult-onset repeat lengths, normalization of neurite length in adult CAG-repeat neural cells, and improved cognition and rescue of synaptic pathology with Isx-9 in R6/2 mice suggest potential strategies to overcome aberrations in developmental networks and abrogate neuronal vulnerability and neurodegeneration in HD. The beneficial effects of Isx-9 in adult mice could indicate a role for developmental pathways that may be linked in adulthood to overall plasticity or response to injury. It will be important to evaluate potential efficacy of Isx-9 in full-length mHTT mouse models. Finally, future iPSC-based studies will clarify the extent of the epigenetic alterations in HD and define the mechanisms for rescue.
HD and non-disease repeat iPSCs were generated and characterized as described15. Unaffected and apparently healthy non-disease (GM05400, 03814 and 02183) and HD (GM09197) human fibroblast cell lines were obtained from the Coriell Institute for Medical Research, under their consent and privacy guidelines as described on their website (http://catalog.coriell.org/). The parental fibroblast for the 109CAG repeat HD line (Coriell, ND39258) was obtained from Johns Hopkins University under Russell Margolis’s IRB protocol #NA00018358. All procedures were performed in accordance with the institutional review board’s guidelines at the Cedars-Sinai Medical Center under the auspices of IRB-SCRO protocols Pro00028429 (Transplantation of iPSC-derived human neural progenitors), Pro00021505 and Pro00032834. Upon iPSC generation at Cedars Sinai, they were renamed using established nomenclature15 as CS00iCTR-21nXX, CS14iCTR-nXX, CS25iCTR-18nXX, CS83iCTR-33nXX, CS13iHD-43nXX, CS03iHD53nXX, and CS109iHD-109nXX to reflect (1) the last two digits of parental lines identifier, (2) non-disease or HD line, (3) CAG repeat number and (4) clone number, represented here by XX15,51. Fibroblasts previously reprogrammed using integrating viral reprogramming methods14 were reprogrammed into nonintegrating and virus-free iPSC lines using the Amaxa Human Dermal Fibroblast Nucleofector Kit to express episomal plasmids with six factors: OCT4, SOX2, KLF4, L-MYC, LIN28, and p53 shRNA (Addgene)52. This method has an advantage over viral transduction because exogenously introduced genes do not integrate and are instead expressed episomally in a transient fashion. Specifically, fibroblasts (0.8 × 106 cells per nucleofection) were harvested, centrifuged at 200g for 5 min and resuspended carefully in Nucleofector Solution (VPD-1001, Lonza), and the U-023 program was applied. All cultures were maintained under standard oxygen conditions (5% O2) during reprogramming, which further enhances the efficiency of iPSC generation. The medium was kept in place for 48 h and gradually changed to chemically defined mTeSR1 medium containing small molecules to enhance reprogramming efficiency. The small molecules used were (1) sodium butyrate (0.5 mM), (2) glycogen synthase kinase 3β inhibitor of the Wnt/β-catenin signaling pathway (CHIR99021, 3 μM), (3) MEK pathway inhibitor (PD 0325901, 0.5 μM) and (4) selective inhibitor of TGFβ type I receptor ALK5 kinase, type I activin/nodal receptor ALK4 and type I nodal receptor ALK7 (A 83-01, 0.5 μM). Individual iPSC colonies with ES/iPSC-like morphology appeared at days 25–32, and those with best morphology were mechanically isolated, transferred to 12-well plates with fresh Matrigel matrix, and maintained in mTeSR1 medium. The iPSC clones were further expanded and scaled up for further analysis, expanded and cryopreserved.
Human iPSCs were rigorously characterized at the Cedars-Sinai iPSC core using several assays. G-band karyotyping ensured a normal karyotype, and genomic DNA PCR confirmed the absence of episomal plasmid genes, as described15,52,53. Pluripotency was assessed by immunostaining with surface and nuclear pluripotency markers for subsequent flow cytometry quantification (>80% SSEA4 and Oct3/4 double positivity) by quantitative RT-PCR of endogenous pluripotency genes and by gene-chip and bioinformatics-based PluriTest assays54. Spontaneous embryoid body differentiation confirmed the capacity to form all germ layers.
Spheres were incubated in colcemid (100 ng/mL; Life Technologies) for 30 min at 37 °C and then dissociated using TrypLE for 10 min. They were then washed in phosphate-buffered saline (PBS) and incubated at 37 °C in 5 mL of hypotonic solution (1 g KCl, 1 g sodium citrate in 400 mL of water) for 30 min. The cells were centrifuged for 2.5 min at 234g and resuspended in fixative (methanol: acetic acid 3:1) at room temperature for 5 min. This was repeated twice, and finally cells were resuspended in 500 μL of fixative solution and submitted to the Cedars-Sinai Clinical Cytogenetics Core for G-band karyotyping.
iPSC colonies grown on Matrigel in TeSR medium (feeder-free) were scraped into EGF- and FGF2-containing (100 ng/ml each) medium (70:30 DMEM:F12 plus 2% B27) and grown as floating neural progenitor spheres for at least nine passages55. Cells for immunocytochemistry or BDNF withdrawal were then be plated on PLO/laminin coated coverslips. Cells for omics studies and western blotting were plated on Matrigel-coated six-well plates. The cells were then differentiated in DMEM:F12 with 1% N2 for 7 d. For the next 21 d, cells were differentiated in 20 ng/ml BDNF (20 ng/ml; Peprotech 450-02), rhShh (200 ng/ml; R&D 1845-SH) and Dkk1 (100 ng/ml; R&D 1096-DK-010) to promote a rostral forebrain fate. Cells were then matured in 20 ng/ml BDNF, dibutyryl cyclic AMP (0.5 mM; Sigma D0260) and valproic acid (0.5 mM; Sigma P4546) for the rest of the differentiation. Medium was half-changed three times per week or as needed.
Cells were fixed in 4% paraformaldehyde (PFA) at room temperature, rinsed with PBS, and permeabilized with 5% normal goat and/or donkey serum containing 0.2% Triton X-100 for 30 min at room temperature. Cells were then labeled with primary antibodies Tuj 1:1,000 (Sigma T8660)56, DARPP32 1:400 (Cell Signaling 19A3)55, Ki-67 1:400 (EMD MAB4190MI)57, GFAP 1:1,000 (DAKO Z0334)58,59, Map2ab 1:500 (Sigma M1406)60, myosin 1:20 (Developmental Studies Hybridoma Bank D7F2)61, Iba1 1:500 (Wako 019-19741)62, Sox17 1:500 (R&D AF1924)63, FoxA2 1:500 (Abnova H00003170-M01)64, O4 1:50 (StemCell 01416)65, nestin (EMD Millipore ABD69, 1:5,000)66, DCX (Aves Labs DCX, 1:1,000)67 or PDGFR 1:500 (Santa Cruz SC-338)68 for 60 min at room temperature or overnight at 4 °C, and then with the appropriate fluorescently tagged secondary antibodies (1:500; Life Technologies A-21202, A-31571, A-11005, A-21206, A-21207, A-31573) for 60 min at room temperature. Hoechst dye was used to label nuclei. Slides were blinded and cells were counted using unbiased stereological software (Stereoinvestigator). Counting criteria included estimating the population of all cells on the coverslip by colocalization of a given marker with the nuclear marker DAPI. Sampling was performed using 20× magnification. Number of experimental replicates is indicated in the figure legend. To estimate the total number of cells on each coverslip of a given population, the fractionator was used to project a given number of cells for the entire population. To ensure that the sampling size was correct it was ensured that each count had a second estimated CE (Schmitz-Hof) of less than 0.1. Significance was measured via one-way ANOVA with a Bonferroni post-analysis comparison.
HD and control i-neurons were differentiated using a similar protocol as described above. iPSCs were made into EZ spheres55 and switched in to DMEM + 0.5× N2 and 0.5× B-27 for 5 d. On day 5, the medium was changed to DMEM + 0.5× N2 and 0.5× B-27 plus 25 ng/mL BDNF. On day 7, medium was changed to DMEM + 0.5× N2 and 0.5× B-27 plus 25 ng/mL BDNF, 10 nM purmorphamine and 100 ng/ml Dickkopf 1 for 3 weeks. On day ~28, cells were transfected with mApple driven off of the MAP-2 promoter12 using Amaxa Nucleofection and placed in to DMEM + 0.5× N2 and 0.5× B-27 supplemented with 25 ng/mL BDNF, 0.5 mM dibutyryl cyclic AMP and 0.5 mM valproic acid on growth factor–reduced Matrigel. On day ~32, 20 μM Isx-9 was added to the medium and the culture was subjected to robotic microscopy and imaged daily for 9 d as described14. Medium was changed every 2–3 d as needed. Images were used to follow individual cells over time and cell death was recorded to assess the effect of Isx-9 on HD and control cells. We used Cox proportional hazards analysis to reveal the cumulative risk of death of the HD lines and the controls. We found evidence of nonproportionality for a subset of the lines, and so we used log-rank tests to assess the differences of survival between the HD and controls lines and the effects of Isx-9. All P values are reported from the log-rank test; however, we report the HRs as an estimate of the hazard from the Cox proportional hazards model. Data shown in Figure 7 are all experiments combined into one curve. The individual survival curves are shown in Supplementary Figure 10. For the 46Q survival curves: 18Qn2/n6 + DMSO, n = 124 cells,18Qn2/n6 + Isx-9, n = 143 cells (4 experiments); 46Qn1/n10 + DMSO, n = 702 cells, 46Qn1/n10 + Isx-9, n = 783 cells (6 experiments). For the 53Q survival curves: Q18n2/6 + DMSO, n = 490 cells, Q18n2/6 + Isx-9, n = 452 cells (7 experiments); 53Qn3/5 + DMSO, 1,047 cells; 53Qn3/5 + Isx-9, n = 1,107 cells (7 experiments).
On day 39 of differentiation, lengths of processes from neurons were measured using Neuron J. To control for intra-experiment variability, each individual process length was normalized by dividing by the average length of the corresponding control in each experiment. In the cases where two controls were used in an experiment, the average of the two was used to normalize. Actual process lengths averaged between ~100 and 400 μm per cell, with an average of 1–3 processes per cell. The data were tested for normality by the Shapiro-Wilk normality test and the variances were not normal; P < 0.0001. Therefore, a one-way ANOVA on ranks (a Kruskal–Wallis test) was performed by GraphPad and P values were reported by a Dunn’s multiple comparisons test. Q18n2 (254 cells, 7 experiments), Q18n2 + Isx-9 (187 cells, 6 experiments), Q18n6 (208 cells, 4 experiments), Q18n6 + Isx-9 (290 cells, 4 experiments), Q28n6 (152 cells, 4 experiments), Q28n6 + Isx-9 (149 cells, 4 experiments), Q46n1 (540 cells, 7 experiments), Q46n1 + Isx-9 (708 cells, 7 experiments), Q46n10 80 cells, 1 experiment), Q46n10 + Isx-9 (168 cells, 1 experiment), Q53n3 (574 cells, 7 experiments), Q53n3 + Isx-9 (640 cells, 7 experiments), Q53n5 (285 cells, 5 experiments), Q53n5 + Isx-9 (387 cells, 5 experiments), Q109n4 (125 cells, 2 experiments), Q109n4 + Isx-9 (119 cells, 2 experiments) and Q109n5 (244 cells, 8 experiments), Q109n5-9 (201 cells, 7 experiments).
Neuronal cells were stained as described above using antibodies against MAP-2 (1:1,000, Abcam, ab5392 against all three isoforms), DARPP32 (1:200, Cell Signaling 19A3) and Ki-67 (1:200, Fisher Scientific MAB419MI). Secondary antibodies were goat anti-chicken for MAP-2 (Life Technologies catalog no. A-21437) and donkey anti-rabbit for DARPP-32 and Ki-67 (Life Technologies catalog no. A-21206).
RNA was isolated from cell pellets with a Qiagen RNeasy Kit with QIAshredder homogenization. RNA-seq libraries were made with 1 μg of RNA using the Illumina Tru-Seq v2 protocol. Final PCR products were run on 2% agarose E-Gels (Invitrogen) and libraries were size-selected at 300 bases. Libraries were sequenced on one lane of the Hi-Seq 2000 at 50 bases. Reads were trimmed to 40 bases and mapped to the hg18 genome with Bowtie (0.12.7) with settings --best --m 1 --v 2. Counts per gene were calculated by counting reads mapping to constitutive exons; counts for each gene were then analyzed with the R package DESeq (1.0.6) to identify differentially expressed genes with a single comparison of HD versus non-disease. A 10% false discovery rate cutoff and log2 difference of 0.5 between non-disease and HD conditions was used for significance. Outliers were further excluded by restricting the residual variance quotients to less than 10. Counts were also used to calculate reads per kilobase of exon per million mapped reads (RPKM) values for each gene. Data are shown in RPKMs. Data were analyzed through the use of QIAGEN’s Ingenuity Pathway Analysis.
Total RNA was isolated using TRI Reagent (Sigma-Aldrich T9424) following the manufacturers’ protocol. Ten microliters of total RNA at a concentration of 200 ng/μL (2 μg in total) for each sample was reverse-transcribed with random primers using the High-Capacity RNA-to-cDNA Kit (Life Technologies 4387406). Ten microliters of retrotranscription cocktail (2 μL of 10× reverse transcription buffer, 2 μL of random primers, 1 μL of dNTP mix; 1 μL of MultiScribe reverse transcriptase) were added to each sample (20 μL total volume). After gentle mixing, the samples were incubated for 10 min at room temperature followed by 2 h at 37°, 10 min on ice and 10 min at 75°.
Customized Openarrays (Life Technologies) containing 112 TaqMan probes that specifically detect all isoforms of each of 106 genes, selected from the literature, and six housekeeping genes (18S, B2M, HPRT1, HSP90AB1, RPL13A, UBC) as reference genes, were produced and validated. cDNAs were loaded onto the custom Openarrays and run as recommended by the manufacturer on the QuantStudio 12K Flex Real-Time PCR system by Servei Veterinari de Genètica Molecular at Faculty of Veterinary in Universitat Autònoma de Barcelona (Spain).
Relative gene expression values were calculated with 2−Ct (ref. 69) using Expression Suite Software 1.0.3 (Life Technologies). RQ minimum and maximum values (error bars) were calculated with a confidence level of 95%, using Benjamini-Hochberg false discovery rate to adjust P values. Maximum allowed Ct included in calculations is 34 and Cq confidence > 0.8. Multivariate Student’s t-test was applied and values of P < 0.05 were considered statistically significant. Error bars are presented in all graphs as s.d. Gene expression profile data are represented in graphics as relative quantity of 2−Ct of 21Q normalized to others.
Total RNA was passed through RNeasy spin columns (Qiagen) for further cleanup. Oligo(dT)-primed cDNA synthesis was performed on 200–300 ng of total RNA using the Superscript III RT Kit (Life Technologies). Primers for qPCR were designed by Origene and ordered through Eurofins MWG Operon. qPCR was performed using SYBR Green (Bio Rad) in a 384-well format on a ViiA 7 real-time PCR system (Applied Biosystems 4309155). Gene expression for each sample was normalized to RPLPO, and expression under treatment conditions further normalized to untreated samples. Statistical analysis of differences in gene expression between vehicle and Isx-9 treatment was performed using an unpaired two-tailed t-test in GraphPad Prism. Primers: NEUROD1 forward: GGTGCCTTGCTATTCTAAGACGC, reverse: GCAAAGCGTCTGAACGAAGGAG, CALB1 forward: TTTCCTGCTGCTCTTCCGATGC, reverse: GCTCCTCAGTTTCTATGAAGCCA, CAMK4 forward: GTTCTTCTTCGCCTCTCACATCC, reverse: CTGTGACGAGTTCTAGGACCAG, POU4F2 forward: TATGCGGAGAGCCTGTCTTCCA, reverse: TCTTGCTCTGGGAGACGATGTC, ATP2B2 forward: CATCCTCAACGAACTCACCTGC, reverse: GATATTGTCGCCAGTGACCATGC, CACNA1A forward: CTGGTAGCCTTTGCCTTCACTG, reverse: ACACAGCCTTGAGCTTTGGCAG, CDK5R1 forward: TCATCTCCGTGCTGCCTTGGAA, reverse: CTCATTGTTGAGGTGCGTGATGT, CHRNA3 forward: TGGAGACCAACCTGTGGCTCAA, reverse: CAGCACAATGTCTGGCTTCCAG, CHRNA4 forward: CGTCCAGTACATTGCAGACCAC, reverse: TCCAGAGGAAGATGCGGTCGAT, GAD1 forward: TGTCCAGGAAGCACCGCCATAA, reverse: TCCTTGACGAGAATGGCAGAGC, GRIA1 forward: GGATGCTCTTTCAGGACCTGGA, reverse: GTAGTGGTAGCCGATGCCATTC, KCNQ3 forward: CGTCTGATTGCCGCCACCTTTT, reverse: TTCTGACGGTGTTGCTCCTGCA.
Mouse development samples spanning embryonic day 12.5 to 18.5 and covering the lateral ganglionic eminence, germinal zone and mantle zone were obtained by laser microdissection and hybridized to Affymetrix chip Mouse_ST1 with ~27,800 geneprobe sets. Using R software with the Bioconductor, affy, and limma packages, microarray data were subjected to quality metrics assessment (arrayQualitymetrics) and robust multiarray normalization. A linear fit algorithm with Bayesian correction was applied grouping biological replicates of samples and designing a contrast matrix that compared zones and regions, giving as a result 3,633 differentially expressed genes (DEGs) that complied with the criteria fold change > 2, adjusted P < 0.01 (Benjamini-Hochberg procedure).
A list of 1,869 human genes was compared with our microarray data once the gene symbol and Uniprot identifiers for mouse genes were added via the Xspecies identifier obtained from original human Refseq on BRM software from Pacific Northwest National Laboratory70, yielding 1,674 matching human–mouse IDs. Matches were further cross-checked with our list of 3,633 DEGs, which resulted in a list consisting of 679 genes. We separately normalized the values given for human samples and the log2 of intensity of probe for mouse data using Genesis Software71. Further heat map representation ranging from −3 to 3 s.d. was thus possible, and hierarchical clustering with average linkage for samples and genes was performed using Euclidian distance.
GeneCodis web-based software was used to obtain GO enrichment of gene lists by modules. The Thomson Reuters platform Metacore was used to obtain curated literature-described relationships between gene products and to perform external network analysis.
ChIP libraries were prepared from the cell lines and the resulting sequences were analyzed as described29. We used antibodies that specifically recognize H3K4me3 (Millipore, Billerica, MA, cat #17-614, 3 μl per ChIP, http://www.merckmillipore.com/DE/en/product/ChIPAb%2B-Trimethyl-Histone-H3-%28Lys4%29---ChIP-Validated-Antibody-and-Primer-Set%2C-rabbit-monoclonal,MM_NF-17-614#anchor, ref. 29), H3K27ac (Abcam, Cambridge, UK, cat #ab4729, 10 μl per ChIP, http://www.abcam.com/histone-h3-acetyl-k27-antibody-chip-grade-ab4729-references.html), and H3K36me3 (Abcam, cat #ab9050, Cambridge, UK, 5 μl (5 μg) per ChIP, http://www.abcam.com/histone-h3-tri-methyl-k36-antibody-chip-grade-ab9050.html)72. Sequences were aligned to the hg19 genome from UCSC using Bowtie with an m parameter of 1 (ref. 29). To call peaks, MACS73 was run on each cell line sample individually, with an mfold parameter of 10, 30 and P-value threshold of 10−5. As MACS results for individual lines were very similar, sequences for all HD cell lines were combined, as were all non-disease cell lines for each histone mark. MACS was run on these samples against a uniform background, with the same parameters. Differential peaks between the non-disease sample and the high-repeat sample were called using these MACS peaks as the input to MAnorm with the --overlap-independent parameter, m value of 1, and P-value threshold of 0.01 (ref. 74). Peaks were mapped to all gene TSSs within 10 kb, according to the RefSeq database, and then lists of those nearby genes were run through DAVID Functional Annotation Clustering tool.
Transcription start sites were taken from the RefSeq database and filtered so that only those sites with no other sites within 7 kb were considered. If more than one TSS was reported for one gene, the TSS that had the most H3K4me3 reads in its immediate vicinity in our cell lines was chosen as the primary TSS. Raw data from each non-disease cell line were first normalized using Fixseq75, and then reads in the region around each TSS were counted and stored as a vector. These vectors were binned, normalized and then clustered by the k-means algorithm using Euclidean distance and complete linkage. Genes were assigned to the cluster to which the vector of reads around its TSS was assigned. Genes in each cluster were analyzed with DAVID Functional Annotation Clustering tool. The average vector of those assigned to each cluster was calculated and shown in Figure 3a. The lists of genes in each cluster was compared to a list of genes that change expression in HD cell lines, and the hypergeometric test was used to assess the significance of the overlap.
Reads from H3K27ac ChIP-seq in the non-disease lines and the four HD lines were combined and MACS and MAnorm were run using the same parameters as above. Resulting differential regions were split into CpG-high and CpG-poor regions, where CpG-poor is defined as having an observed-to-expected CpG ratio, calculated as (number of CpGs)/[(number of C’s × number of G’s)/length of sequence], below 40%, and THEME76 was used to find enriched motifs in the CpG-poor regions of HD-biased peaks and non-disease-biased peaks. Our analysis was restricted to areas of CpG-poor DNA sequence, as this often results in more tissue-specific motif discoveries77, and we disregarded transcription factors not expressed across all six cell lines.
The set of motif hypotheses is derived from all vertebrate-specific scoring matrices (PSSMs) from TRANSFAC78, filtered for sufficient information content (>9 total bits). P-values were calculated using THEME by using the other set of peaks as background sequences (i.e. HD-enriched peaks as foreground and non-disease-enriched peaks as background, and vice versa). THEME calculated the over-representation of the motif in the foreground compared to the background sequences. Motif hypotheses were then clustered according to similar binding site motifs, and the top P value for all motifs in a cluster is reported.
Neural cultures were generated from iPSC monolayer grown on Matrigel-coated six-well plates, using a 37-d protocol79 HTT knockdown was achieved by supplementation of the culture medium with HTT-specific (TCTCTATTGCACATTCCA) or control (CCTTCCCTGAAGGTTCCTCC) ASO15 such that neurons received four doses of 1 μM ASO over 7 d immediately before harvest at day 37. Total RNA was isolated from two wells each of untreated, control-ASO-treated and HTT-ASO-treated neural differentiation with TRIzol and 400 ng of RNA used for RT-PCR, with three technical replicates performed for each gene. For qPCR analysis, a statistical difference in gene expression between control ASO and HTT ASO treatment was determined in GraphPad Prism using an unpaired two-tailed t-test.
Differentiated cultures were transduced at day 39 of culture with 1 × 106 infectious units human NEUROD1 overexpression virus (pLenti-GIII-EF1α vector) (Applied Biological Materials) or 100 ng GFP control (pFUGW-eGFP) (Addgene; Cambridge, MA). All cultures received a complete medium change daily and were harvested at day 53. Gene expression was determined by RT-qPCR as described, beginning with 500 ng total RNA.
Cell cultures received daily medium changes containing 20 μM Isx-9 (Tocris, Minneapolis, MN), beginning at day 37, until harvested at 56 days in vitro (DIV). DMSO was used as the vehicle control. Two wells of each of untreated, vehicle-treated and Isx-9-treated neural differentiation were pooled and total RNA isolated using TRIzol (Life Technologies). RNA (300 ng) was used for RT-PCR, and three technical replicates were performed for each gene. For qPCR analysis, a statistical difference in gene expression between vehicle and Isx-9 treatment was determined in GraphPad Prism using an unpaired two-tailed t-test.
Cells were lysed with RIPA buffer (Sigma, cat# R0278) supplied with 1% Triton X100 (Sigma, cat# 93443) and 0.5% Protease Inhibitor Cocktail (Set III, Calbiochem, now a part of MilliporeSigma, cat# 539134) and samples were sonicated. Protein concentrations were determined using Pierce BCA Protein Assay Kit (ThermoFisher, cat# 23225), according to manufacturer’s directions. Approximately 15 μg of protein was loaded into NuPAGE 4–12% Bis-Tris protein gels (ThermoFisher, cat# NP0321BOX), separated with 200 V for approximately 1 h, and electrotransferred to PVDF membrane in the XCell II Blot Module for wet (tank) transfer (ThermoFisher) for 90 min at 30 V. The membrane was blocked in 5% non-fat dry milk in TBS buffer, pH 7.4 (Quality Biological, cat# 351-086-101) plus 0.1% Tween 20 (Sigma-Aldrich) for 1 h at room temperature and then exposed to primary antibodies. Antibodies for detection of NEUROD1, CALB1, and CAMKIV were purchased from Cell Signaling: NeuroD (D90G12) rabbit mAb: cat# 7019, 1:1,000 dilution (for WB validation, see https://www.cellsignal.com/products/primary-antibodies/neurod-d90g12-rabbit-mab/7019); Calbindin (D1I4Q) XP rabbit mAb: cat# 13176, 1:1,500 dilution (for WB validation and references, see https://www.cellsignal.com/products/primary-antibodies/calbindin-d1i4q-xp-rabbit-mab/13176?_=1487290902164&Ntt=13176&tahead=true)80–82; CaMK IV: rabbit, cat# 4032, 1:1,000 dilution (for WB validation and references, see https://www.cellsignal.com/products/primary-antibodies/camkiv-antibody/4032; Product Citations)83,84. Equal loading was verified by western blotting of actin (α-actin: 1:10,000, mouse, Sigma-Aldrich, clone AC-15, cat# A5441; for WB validation see http://www.sigmaaldrich.com/catalog/product/sigma/a5441?lang=en®ion=US). Densitometry was done using ImageJ software and statistical analysis was performed using one-way ANOVA for four independent experiments.
Nuclear condensation assay of iPSC-neurons was performed as described14. Cultures were supplemented with 20 μM Isx-9 or 20 ng/ml BDNF. Briefly, cells were automatically analyzed under an inverted fluorescence microscope (Axiovert 200, Zeiss) and images were digitized from 49 independent fields per well. Quantification of cell survival was done using the Volocity software (PerkinElmer) by automated measurement of the average intensity of Hoechst 33342–stained nuclei of transfected cells visualized by eGFP expression. Cells were considered viable when their intensity was lower than 200% of the control intensity. Statistical analysis was performed using one-way ANOVA for at least three independent experiments.
The HD iPSC derived EZ-spheres were differentiated in 24-well plates coated with Matrigel (BD), transduced with pGFP-sh lentiviral particles carrying vectors encoding either scrambled control or NEUROD1 shRNA (KD) (Origene) for 24 h, and then transferred to neural induction medium (NIM) without additives or supplied with 20 μM Isx-9 for 48 h. Lentivirus was used at concentration 100 ng/ml. The nuclear condensation assay performed as described14.
All animal procedures were performed in accordance with National Institutes of Health and University of California guidelines. R6/2 transgenic (~120 ± 5 CAG repeat) mice and their non-transgenic littermates obtained from the Jackson Laboratory. Isx-9 was prepared as 2 mg/ml of 30% vehicle (2-hydroxypropyl-β-cyclodextrin in sterile milliQ-purified H2O). R6/2 and age-matched non-transgenic mice (n = 10) were obtained at 5 weeks of age and given a previously established dose of Isx-9 (20 mg/kg) by daily intraperitoneal injection after a brief acclimation period. Mice were tested in behavior tasks each week thereafter and then sacrificed at 10 weeks of age. All mice were housed on a 12 h light/dark schedule with ad libitum access to food and water.
Mice were assigned to groups in a semirandomized manner. Behavioral tests were performed at 6, 7, 8, 9 or 10 weeks of age depending on the task. Researchers were blind to which mice had been injected during experiment testing and data collection. To minimize experimenter variability a single investigator conducted each behavioral test. Rotarod, pole test and grip strength were performed as described85.
Each mouse was exposed to two objects and allowed to explore for a period of time, after which one of the objects was replaced by another. If memory is functioning normally, the mouse will spend more time exploring this novel object than it does exploring the familiar object. If exploration of all objects is the same, this can be interpreted as a memory deficit. A white plastic tub was used for testing and a mounted digital camera recorded movement. Two different objects were tested to ensure that the animals do not have an innate bias toward one or another: a 100 mL glass beaker turned upside down versus a two-hump tall Mega Lego block. Lighting was adjusted to 25 lux using black trash bags taped over the overhead lights. Mice were acclimated to the testing room then placed in the tubs without any objects to freely explore the area for 5 min (habituation to box). The next day after acclimation mice were presented the same object (either two beakers or two blocks) for 5 min. On the third day mice were presented with a new object replacing one old. Using the recorded data the amount of time that the animal interacts with each object was scored in seconds. An exploratory preference score (time spent exploring the novel object divided by the total time spent exploring box objects × 100) was calculated. An exploratory preference score of 50% indicates chance performance, and a higher preference score indicates the presence of memory.
Mice were transcardially perfused with 4% PFA and, following dissection, tissue was postfixed for a further 2 h before being placed in a 30% sucrose solution for 48 h. The tissue was subsequently embedded in a 2:1 mixture of 30% sucrose:OCT and frozen. Sections (15 μm) were cut with a cryostat before being mounted onto glass slides. Sections were then incubated with blocking solution (10% normal goat serum and 0.4% Triton) followed by incubation with antibodies to Homer-1 (1:200, Synaptic Systems 106003) and VGLUT1 (1:1000, Millipore AB5905) in the blocking solution (both antibodies have been used for this type of analysis86). Subsequently sections were washed and incubated with appropriate secondary antibodies (Invitrogen #A11034 and A11076), before being washed again and coverslips mounted onto the slides with mounting medium (Vector Shield). To visualize synaptic staining sections were imaged on a Zeiss LSM 700 confocal microscope at 63× magnification. Digital pictographs were captured and Image J software was used to quantify colocalized pre and postsynaptic puncta. All data analysis was carried out blind to the conditions of the experiment.
Images were taken every 24 h for 9 d as described87,88 on an inverted microscope (Nikon Ti Eclipse) that contains the PerfectFocus system, a 20× Plan Fluor S 0.45 NA ELWD objective and a 16-bit electron multiplying charge-coupled device (EMCCD) cooled to −70 °C (Andor iXon 888). Illumination was delivered from a Sutter Lambda XL Xenon lamp. Semrock BrightLine full-multiband filter sets for DAPI, FITC and TRITC were used for excitation and detection. The illumination, filter wheels, focusing, stage movements and image acquisitions were fully automated and coordinated with publicly available (ImageJ and μManager) software and a custom-written image acquisition plug-in. During longitudinal experiments, plates were kept in a robotically controlled incubator (Liconic STX44) and were delivered to the microscope stage that was housed in an OKO-Labs environmental chamber by an automated robotic arm (Peak Robotics). Scheduling and automation was controlled with Green Buttton Go (Biosero).
Images from differentiated i-neurons were montaged using custom scripts written in Pipeline Pilot PipelinePilot (Accelrys, San Diego, CA). The images were then tracked in MATLAB using custom-based algorithms developed in our lab. These programs allowed us to keep track of the amount of time each cell lived and died during the course of 9 d. Cells were tracked only if they had neuronal morphology as reported14. Cell death was ascribed to cells that showed morphological features of cell death such as cellular blebbing, cell rupture or loss of fluorescent signal87. Survival time for each cell was denoted as the last time the cell was seen alive. Scripts written in R’s survival package were used to generate cumulative risk of death curves and to perform Cox proportional hazards analysis88 or log-rank test89 to assess the relative risk of death between the HD i- neurons and controls and the effect of Isx-9 on survival. The cumulative risk of death curves were assessed for violations of proportionality by evaluating both the cumulative risk of death curves and by using the using cox.zph function in the R survival package.
There were no violations of proportionality except when analyzing the Q18n2 + DMSO and Q18n2 + Isx-9 and 53Qn3 + DMSO survival probabilities in the Q53 data sets. We found strong evidence for nonproportional hazards between the two groups. In particular, we observed crossing Kaplan–Meier survival curves and a cox.zph with P values of 2.97 × 10−5 for Q18n2 + DMSO, 1.38 × 10−4 for Q18n2 + Isx-9 and 2.00 × 10−2 53Qn3 + DMSO. Further inspection of the data revealed that the cause for this was that one of the repeated experiments for group did not match the typical trends for the other runs of the experiment. While for each individual run of the experiment the proportional hazard assumption appeared appropriate, when combining all runs into a single experiment group the outlier experiment was the main cause of the nonproportional hazards. Rather than drop the outlier group, we used a log-rank test that does not require the assumptions of the proportional hazards to evaluate the differences in survival between the Q53 and control data set. To estimate the approximate HR between groups, we used the HR values from the Cox proportional hazards model for this data set.
Multiple statistical models were used during these analyses and are detailed further in the in the figure legends, Methods above, and below. A Supplementary Methods Checklist is available with specific details on statistical tests and methodology used for each assay. Samples sizes were chosen based on previous expertise and knowledge from past experiments using similar model systems, mice and iPSCs, and current accepted standards based on literature review (for example, refs. 14,29). A power analysis was conducted for animal studies and is detailed below. For most data, data distribution was assumed to be normal but this was not formally tested. Principal component analysis was used to assess variance in global gene expression between HD and control iPSC-derived neural cells. Where appropriate, individual data points are displayed to allow visualization of distribution. Where appropriate, samples and animals were randomized for experimentation, and the investigator was blinded to groups and treatments. These details are described below.
To ensure that the sampling size was correct, it was ensured that each count had a second estimated CE (Schmitz-Hof) of less than 0.1. Significance was measured via one-way ANOVA with a Bonferroni post-analysis comparison.
Western blot results (densitometry) were normalized to actin. Statistical significance between Isx-9 and untreated cells (NIM) was determined by one-way ANOVA.
Counts per gene were calculated by counting reads mapping to constitutive exons; counts for each gene were then analyzed with the R package DESeq (1.0.6) to identify differentially expressed genes with a single comparison of HD versus non-disease. A 10% false discovery rate cutoff and log2 difference of 0.5 between non-disease and HD conditions was used for significance. Outliers were further excluded by restricting the residual variance quotients to less than 10.
Relative gene expression values were calculated with 2−Ct (ref. 69) using Expression Suite Software 1.0.3 (Life Technologies). RQ minimum and maximum values (error bars) were calculated with a confidence level of 95%, using Benjamini-Hochberg false discovery rate to adjust P values. Maximum allowed Ct included in calculations is 34 and Cq confidence > 0.8. Multivariate Student’s t-test was applied and values of P < 0.05 were considered statistically significant. Error bars are presented in all graphs as s.d. Gene expression profile data are represented in graphics as relative quantity of 2−Ct of 21Q normalized to others.
Statistical analysis of differences in gene expression between vehicle and Isx-9 treatment was performed using an unpaired two-tailed t-test. A statistical difference in gene expression between control ASO and HTT ASO treatment was determined in GraphPad Prism using an unpaired two-tailed t-test.
Using R software with the Bioconductor, affy and limma packages, microarray data were subjected to quality metrics assessment (arrayQualitymetrics) and robust multiarray normalization. A linear fit algorithm with Bayesian correction was applied grouping biological replicates of samples and designing a contrast matrix that compared zones and regions giving as a result 3,633 differentially expressed genes (DEGs) that complied with criteria fold change > 2, adjusted P < 0.01 (Benjamini-Hochberg procedure).
Peaks were called and differential regions were determined using MACS and MAnorm, respectively, with statistical thresholds as described above. Functional annotation terms for nearby genes were assigned FDRs by the DAVID Functional Annotation Clustering tool. Motifs near differential H3K27ac peaks were assigned P values by the THEME software. Genes were then assigned to epigenetic profile classes, and statistically significant overlap between these classes and transcriptionally dysregulated genes was determined using the hypergeometric test, where the background was all genes detected by RNA-seq, and P values were calculated by a two-tailed Fisher’s exact test.
Our experiments were done blindly using a coding system for the different treatments. The placement of the treatments on the plates was randomized with each experiment. All statistics are shown for independent experiments. No power analysis was done a priori for these experiments since we had performed power analysis for similar conditions previously14.
In mouse behavior tasks statistical tests used one-way ANOVA followed by Tukey HSD test with Scheffé, Bonferroni and Holm multiple comparison calculation also performed post hoc. Data met the assumptions of the statistical tests used and P values less than 0.05 were considered significant. In figures for behavior tasks individual data points are shown and data distribution was assumed to be normal, but this was not formally tested. To determine our group and trial sizes in the mouse experiments we used G Power (http://www.psycho.uni-duesseldorf.de/abteilungen/aap/gpower3/). Assessment of differences in outcome were based upon previous experience85 and published results for R6/2 (ref. 90). All mice were randomly assigned to treatment groups and behavior tasks performed in a random manner, with the individual performing the experiment blinded to group genotypes and treatment status. R6/2 and age-matched NT male mice (n = 10/group) were obtained at 5 weeks of age and tested in behavior tasks each week thereafter until sacrifice at the end of 10 weeks of age. In the novel object recognition task, results from 2 mice each from the NT groups and 4 mice each from the R6/2 groups had to be excluded because an inappropriate novel object was used in the assays. The novel object was different from that used with the other mice.
Statistical differences between groups were analyzed with an unpaired two-tailed t-test. No statistical methods were used to predetermine sample sizes, but our sample sizes are similar to those reported90. Data distribution was assumed to be normal, but no formal test was carried out. No animals or data points were excluded from the analysis.
Scripts written in R’s survival package were used to generate cumulative risk of death curves and to perform Cox proportional hazards analysis88 or log-rank test89 to assess the relative risk of death between the HD i- neurons and controls and the effect of Isx-9 on survival. The cumulative-risk-of-death curves were assessed for violations of proportionality by evaluating both the cumulative risk of death curves and using the cox.zph function in the R survival package. While for each individual run of the experiment the proportional hazard assumption appeared appropriate, when combining all runs into a single experiment group the outlier experiment was the main cause of the non-proportional hazards. Rather than drop the outlier group, we used a log-rank test that does not require the assumption of proportional hazards to evaluate the differences in survival between the Q53 and control data set. To estimate the approximate HR between groups, we used the HR values from the Cox proportional hazards model for this data set.
Scripts for the k-means clustering of epigenetic profiles and overlap of these profiles with transcriptional data are available on GitHub at https://github.com/fraenkel-lab/HD-iPSC-Consortium-NN. Any additional code used in this study is available from the corresponding author upon reasonable request.
ChIP-seq and RNA-seq data is available in the Gene Expression Omnibus database under accession code GSE95344. Lists of genes in each epigenetic class used in Figure 5 are available on GitHub at https://github.com/fraenkel-lab/HD-iPSC-Consortium-NN. Any additional data that support the findings of this study are available from the corresponding author upon reasonable request.
We thank the patients and their families for their essential contributions to this research. We also thank E. Cattaneo and J. Arjomond for discussions of the data, D. Merry for critique of the manuscript and data, S. Svendsen for editorial assistance, J. Dunn for technical culture assistance, G. Vatine (Cedars-Sinai Medical Center, Los Angeles) for the iPSC-derived oligodendrocyte precursors, M. Godoy and G. Gowing (Cedars-Sinai Medical Center, Los Angeles) for the rat muscle and R. Barrett (Cedars-Sinai Medical Center, Los Angeles) for the definitive endoderm positive control. We also thank F. Bennet and Ionis Pharmaceuticals for providing the HTT ASO. Primary support for this work was from NIH NS078370 (L.M.T., C.N.S., J.F.G., M.E.M., C.A.R. and S.F.) and from the CHDI Foundation (J.M.C, P.J.K. and N.D.A.). Additional support was provided by NIH: U54 NS091046 NeuroLINCS center (L.M.T., C.N.S., E.F.); NIH NS089076 (L.M.T., D.E.H., E.F.); P50NS16367 (HD Center Without Walls); NIH R01GM089903 (E.F.); NIH NS101996-01 (S.F.); NIH R01NS084298 (B.S.); American Heart Association, CIRM and NRSA fellowships (R.G.L.); the Hereditary Disease Foundation (V.B.M.); the Taube-Koret Center and the Hellman Family Foundation (S.F.); the UCI Institute for Clinical and Translational Science (L.M.T.); Huntington’s Disease Society of America (L.L.S); HD CARE (L.M.T.) and the NIH Biotechnology Training Program Fellowship (T32GM008334, A.J.K.). Additional support was provided by grants from the Ministerio de Economia y Competitividad (SAF 2014-57160-R to JA; SAF2015-66505-R to J.M.C.), from the ISCIII-Subdirección General de Evaluación and European Regional Development Fund (ERDF) (RETICS to JMC (RD12/0019/0002; Red de Terapia Celular); ADVANCE(CAT) with the support of ACCIÓ (Catalonia Trade & Investment; Generalitat de Catalunya) and the European Community under the Catalonian ERDF operational program 2014-2020), Spain, from the European Union FP7 (P.J.K. and N.D.A.) and from the Ser Cymru Life Sciences & Health Network in Drug Discovery Programme (M.W.S.). This work was made possible, in part, through access to the Genomic High Throughput Facility Shared Resource of the Cancer Center Support Grant (CA-62203) at the University of California, Irvine. Support also included computing resources from National Science Foundation grant DB1-0821391 and sequencing support from National Institutes of Health grant P30-ES002109.
AUTHOR CONTRIBUTIONSDesigned the experiments: R.G.L., L.L.S., D.K.W., J.C.R., S.T.W., L.M.T., J.A., J.M.C., N.D.A., P.J.K., J.A.K., S.F., F.Y., D.E.H., E.F., J.F.G., M.E.M., S.S.A., N.A., C.A.R., V.B.M. and C.N.S. Generated iPSC lines in study: L.O., A.S., L.L., B.M. and D. Sareen. iPSC culture and neuronal differentiation: V.B.M., L.L.S., A.R.K., J.T.S., C.M.T., S.S.A., J.A.K., H.M. and M.D. Carried out experiments: R.G.L. and T.A.G., RNA-seq; L.L.S., Isx-9 qPCR and NEUROD1 overexpression; A.L.L., M.P., A.G.G.D.-B., M.S. and P.S., mouse neurodevelopment studies; A.G.G.D.-B., comparison between mouse and human data; V.B.M. and C.M.T., cell counts; V.B.M., immunocytochemistry; F.Y., R.S.A. and M.B., ChIP; S.S.A., Cell Titer-Glo cell survival assay; S.S.A., N.A. and L.L.S., NEUROD1 knockdown; N.A. and J.S., cell culture and transfection of mouse primary neurons and nuclear condensation assay; S.S.A., Western analysis; E.M., J.A.K., M.D. and H.M., Isx-9 neuron assays; D.K.W., Isx-9 synaptic assays; J.C.R. and D. Sharifabad, mouse Isx-9 studies. Analyzed the data: R.G.L., L.L.S., D.K.W., B.S., J.C.R., M.S.C., S.T.W., L.M.T., J.A.K., M.D., H.M., L.J., D.K.W., C.A.-B., S.F., A.J.K., T.A.G., F.Y., C.W.N., P.M., D.E.H., E.F., J.F.G., M.E.M., S.S.A., N.A., C.A.R., V.B.M. and C.N.S. Wrote the manuscript: R.G.L., L.L.S., D.K.W., J.C.R., S.T.W., L.M.T., J.M.C., N.D.A., P.J.K., J.A.K., K.H., S.F., A.J.K., T.A.G., P.M., D.E.H., E.F., M.E.M., J.F.G., S.S.A., C.A.R., V.B.M. and C.N.S. A list of authors by individual consortium group appears in the Supplementary Note.
COMPETING FINANCIAL INTERESTS
The authors declare no competing financial interests.
Note: Any Supplementary Information and Source Data files are available in the online version of the paper.