Human episomal HD and control iPSC lines were obtained from the NINDS Human Genetics Resource Center DNA and Cell Line Repository¹. For the analysis, we used three clonal HD lines with 71 CAG repeats (ND42228, ND42229, ND42230; derived from a 20-year-old patient), three juvenile HD clonal lines with 109 CAG repeats (ND42222, ND42223, ND42224; derived from a 9-year-old patient) control lines (two clonal lines with 17/18 (ND41654, ND41658) and one line with 21 (ND42245) CAG repeats. Human iPSCs were cultured in chemically defined conditions in Essential 8 medium (Life Technologies) and grown on recombinant human vitronectin-coated surfaces (VTN-N, Life Technologies). Cells were passaged using gentle dissociation with 0.5 mM EDTA in PBS.

After medium removal, iPS cells were washed once with PBS and subsequently covered with 1 mL of RNAzol RT RNA Isolation Reagent (GeneCopoeia, Inc.), scraped and frozen in -80°C. Upon thaw, total RNA isolation was performed according to the manufacturer’s protocol with 75% ethanol, isopropanol, and RNase-free water. Isolated total RNA was then treated with TURBO DNase (ambion) and purified with QIAquick^® Nucleotide Removal Kit (QIAGEN). Each reaction was performed in PCR probes and contained 10 μg of RNA sample, 10x TURBO Buffer, RNaseOUT (Invitrogen), DNase TURBO and RNase-free water. The reaction mix was incubated in 37°C for 30 min and transferred to 1.5 ml Eppendorf tubes. Purification was performed according to the manufacturer’s protocol. RNA content was measured on NanoDrop. RNA quality and integrity were validated using capillary electrophoresis (2100 Bioanalyzer, Agilent). Average RNA Integrity Number (RIN) for all samples was 9.78, where 10 is the highest score.

For RNA-seq, total RNA was extracted from 9 samples, including 3 from control iPSC lines, 3 from clonal iPSC lines with 71 CAG repeats and 3 from clonal iPSC lines with 109 CAG repeats. RNA was analyzed with Bioanalyzer 2100 and RNA 6000 Nano kit (Agilent Technologies, Santa Clara, CA, United States). Libraries were prepared from 5 μg of total RNA with RIN > 9.3, using KAPA Stranded mRNA-seq Kit (Kapa Biosystems, Wilmington, MA, United States) and NEBNext Multiplex Oligos for Illumina (New England Biolabs, Ipswich, MA, United States). 9 pM-indexed libraries were sequenced with the use of Illumina HiSeq2000, rapid run with paired-end 75 bp reads. On average, 65 mln reads were collected per library (Supplementary Figure S1).

NGSQCToolkit (v 2.2.3) (Patel and Jain, 2012) was used to generate quality metrics for assessment of FASTQ input files. Statistics generated by the above-mentioned tool was used to generate two charts which show the average quality score on the base positions for iPSC lines separately in R. Fastq files were then aligned to the GRCh38.p10 reference genome using STAR (v 2.5.3a) (Dobin et al., 2013) with default parameters suggested by QoRTs (v 1.3.0) (Hartley and Mullikin, 2015) used for genes and exons hits count. After alignment and hits counting, the differential gene expression was calculated using the Bioconductor package DESeq2 (v 1.14.1) (Love et al., 2014) and it consisted mainly of estimation of size factor and dispersion. The genes with significant differential expression were detected by DESeq2 pipeline and manually checked in order to find genes where outliers (automatic outliers detection with replacing their expression by mean across replicates is recommended for a larger number of replicates) influenced the mean expression of replicates. Dysregulated X and Y chromosome genes were excluded from DESeq2 analysis. After calculation of the differential expression by DESeq2, the shrinkage of effect size was performed (function lfcShrink of DESeq2 package). Then, for the results, the base mean expression against the log2 fold change (MA plots) were plotted with significant genes marked in red (plotMA from DESeq2 package). The gene counts were transformed to the log2 scale and normalized with respect to the library size (rlog function). Then, the Principal Component Analysis (PCA) was performed, and a chart was generated automatically by plotPCA DESeq2 function. For each normalized gene count for each sample the total mean gene expression was subtracted and based on the result, the heat map chart was prepared for manually selected list of genes (pheatmap function from pheatmap_1.0.10 package). All the RNA-seq data were uploaded to the GEO repository (accession number GSE124664).

Reverse transcription was performed using Maxima H Minus Reverse Transcriptase (Thermo Fisher) (200U per reaction) on 2 μg of RNA in 20 μl of total reaction according to the manufacturer’s protocol. For priming a mixture of random hexamers (25 pmol) and oligo(dT) 18 (25 pmol) was used. Additionally, RiboLock RNase inhibitor was added to the reaction mix (20U). Before adding the enzyme and the inhibitor, templates were denatured in 65°C for 5 min; after mixing all reaction reagents reaction was incubated for 10 min at 25°C followed by 15 min at 50°C. Resulting cDNA was further 10 times diluted with nuclease-free water and stored in -20°C. RT- controls were included. RT-PCR products from all HD lines were checked for reference gene (GAPDH) and several pluripotency markers, including SOX2, NANOG, OCT4 and LIN28A (Supplementary Figure S2).

qPCR was performed using HOT FIREPol^® EvaGreen^® qPCR Mix Plus (ROX) (SOLIS BIODYNE) on 1 μl of cDNA in 10 μl of total reaction volume. Reaction mix included 250/125 nM primers. Primers are listed in Supplementary Table S1. Thermocycling parameters were as follows: 12 min of initial denaturation at 95°C and 45 three-step cycles with 15 s denaturation at 95°C; 20 s annealing at 60–64°C and 20 s elongation at 72°C. The reaction was run on CFX96 instrument (Bio-Rad). Specificity of reaction for each primer pair was confirmed by agarose gel electrophoresis and EtBr staining. For each primer calibration curves were prepared within 10 1 to 10 -4 cDNA concentration; slope, y intercepts, PCR efficiency, r 2, linear dynamic range an CIs were calculated. Details for each primer set are included in Supplementary Table S1. Each qPCR plate included a calibrator, non-treated control (NTC) and samples were run in duplicates. For each cDNA template, a corresponding RT- sample was amplified; none of the analyzed samples had a positive signal of RT- reaction. Data was obtained and analyzed using CFX Manager 3.1 (Bio-Rad). Cq values were determined in software using the regression model and were exported to Excel for further analyses. GAPDH, PGK1, and C1orf43 were used as reference genes.

Log2 NRQ was calculated as follows for each gene. First, a mean Cq was calculated for each sample as an arithmetical mean of technical replicates Cq values. Then the mean Cq values were corrected for a gene’s amplification efficiency by its multiplication by Log2 of E. E values was calculated as a (%Efficiency ^∗ 0.01 + 1). Next, Relative Quantities (RQ) were calculated for each sample as the exponentiation of the efficiency corrected mean Cq of a sample subtracted from the arithmetic mean Cq of control samples mean Cqs, with the base of E, which is Eˆ[(arithmetic mean of control samples Cqs) – sample Cq]. Normalized Relative Quantities (NRQ) for each sample were calculated by dividing a gene of interest (GOI) RQ by Normalization Factor (NF). Normalization factor for each sample was calculated as a geometrical mean of Reference Genes (REF) RQs of that sample. Next, a base 2 logarithm was calculated from NRQ values and the resulting symmetrical Log2 NRQ value for each sample and GOI was visualized as a dot in a scatter plot together of Mean Log2 NRQ for each genotype. To calculate 95% CI in an RQ for genotype (geometrical mean of mean genotype Cqs) and genotype Standard deviation (from Mean genotype Cqs) were calculated. Based on these values NF and NRQ were calculated as before (genotype NF and genotype NRQ). Error propagation was calculated in subsequent steps: RQ SD = genotype SD ^∗ genotype RQ ^∗ ln(E), then SD NF = NF ^∗ sqrt(sum(SD RQ/X^∗ genotype RQ)ˆ2)REF1toX), where X is number of REF genes; SD NRQ = genotype NRQ ^∗ sqrt((SD NF/NF)ˆ2 + (SD RQ/RQ)ˆ2), then SE NRQ = SD NRQ/sqrt(genotype N), where N is number of samples per genotype, then SE Log2 NRQ = SE NRQ/(genotype NRQ ^∗ ln(2)). The final Combined-SE Log2 NRQ, which combines SE of a disease genotype with SE of a control genotype is sqrt ((SE log2 NRQ dis)ˆ2 + (SE log2 NRQ ctrl)ˆ2). Upper and lower 95% CI limits were as follows: CI = Log2 NRQ ± s^∗Combined-SE Log2 NRQ, where s is a scaling parameter value from a t-Student table for N-1 degrees of freedom and 0.05 confidence.

Two open source platforms were used for the bioinformatics analysis of significantly differentiated transcripts, including ConsensusPathDB (CPDB) (Herwig et al., 2016), and ClueGO (Cytoscape plug-in) (Bindea et al., 2009, 2013). A simple meta-analysis, considering data from several publications, was also performed. As considered to be significantly dysregulated, transcripts with padj < 0.05 were submitted to the analyses. As genes identifier type, HUGO Gene Nomenclature Comity symbols (HGNC) symbols were chosen.

Each cell line was collected for lysis in buffer containing 1 M triethylammonium bicarbonate (TEAB) and 0.1% SDS. The material was subjected to a threefold cycle of freezing and thawing followed by bath sonication for 3-min repeated three times while cooling the tube on ice in between the sonication. The material was centrifuged in 8,000 × g for 10 min and the supernatant was collected. Lysates were then ultra-filtrated using Amicon Ultra-2 Centrifugal Filter 3 kDa Devices (MERCK) to remove any remnants of culture media. The device was pre-rinsed with milli-Q water and 600 μL of lysate was subjected to a fivefold cycle of centrifugation with a fresh amount of TEAB buffer in 7,500 × g for 5 min. Proteins in clear concentrate were estimated using 2-D Quant Kit (GE Healthcare Life Sciences). Ten μg of total protein per sample were diluted with 15 μl of 50 mM NH4HCO3, reduced with 5.6 mM DTT for 5 min at 95°C followed by alkylation with 5 mM iodoacetamide for 20 min in the dark at RT. Subsequently, the proteins were digested with 0.2 μg of sequencing-grade trypsin (Promega) overnight at 37°C followed by label free quantitative proteomics.

Analysis of protein extracts was done on Dionex UltiMate 3000 RSLC nanoLC System coupled with QExactive Orbitrap mass spectrometer (Thermo Fisher Scientific). Peptides derived from in-solution digestion of proteins were separated on a reverse phase Acclaim PepMap RSLC nanoViper C18 column (75 μm × 25 cm, 2 μm granulation) using acetonitrile gradient (from 4 to 60%, in 0.1% formic acid) at 30°C and a flow rate of 300 nL/min (for 185 min). Mass spectra were acquired in a semi-targeted method using two analysis modes. First was a classic data-dependent mode with top 10 data-dependent MS/MS scans, and second, was scheduled MS/MS mode with inclusion list containing peptide sequences chosen based on protein targets selected from previous transcriptome analysis. The target value for the full scan MS spectra was set to 1e6 with a maximum injection time of 100 ms and a resolution of 70,000 at m/z 400. The 10 most intense ions charged two or more were selected with an isolation window of 2 Da and fragmented by a higher energy collisional dissociation with NCE 28. The ion target value for MS/MS was set to 5e4 with a maximum injection time of 100 ms and a resolution of 17,500 at m/z 400.

Protein identification was performed using UniProt human database (March 2017, 137404 entries) with a precision tolerance 10 ppm for peptide masses and 0.08 Da for fragment ion masses. For protein identification and quantification, all raw data obtained were analyzed using MaxQuant 1.5.3.30 (Max Planck Institute of Biochemistry, Munich) (Cox and Mann, 2008). Obtained normalized data were imported to Perseus 1.6.1.3 software (Max Planck Institute of Biochemistry, Munich) (Tyanova et al., 2016). All numeric values corresponding to protein intensity were transformed to a logarithmic scale, and all samples were grouped using categorical annotation. Missing values were then replaced by imputation and PCA analysis was performed. For protein differentiation, ANOVA test was performed with the p-value calculation used for protein truncation. The proteins were annotated for chromosome localization and all proteins translated from genes located on X and Y chromosomes were removed from the list. Then, for the clustering analysis, the data were normalized for each compound using the Z-score algorithm.

Strand-specific RNA-seq of the whole transcriptome was performed using clonal lines from 2 HD patients with 71 or 109 CAG repeats in exon 1 of the HTT gene, and 2 healthy individuals, with 17/18 and 21 CAG repeats, to comprehensively identify mRNAs related to HD. A DESeq2 pipeline for transcripts of HD was developed to identify significantly dysregulated mRNAs. During the RNA-seq data analysis, we compared 6 HD iPSC lines derived from both patients (71Q and 109Q) with control lines and also compared separately each set of three HD lines from one patient with three control lines. As a result of such an approach, we generated statistical values for three comparison groups, HD vs. WT, HD71Q vs. WT, and HD109Q vs. WT (Supplementary Table S3). Heat map diagram, PCA graph and MA plots were generated for differential gene expression analysis (Figures 1A–D and Supplementary Figures S3, S4).

We set the padj cutoff value to create lists of significantly dysregulated mRNAs in each comparison group. Transcripts with padj < 0.05 were considered to be differentially expressed. In the first group, comparing HD iPSC lines (71Q and 109Q) with control lines, 107 significantly dysregulated mRNAs were identified (31 downregulated and 76 upregulated). In HD71Q iPSC lines 198 differentially expressed mRNAs were identified (64 downregulated and 134 upregulated). In the last group, in which HD109Q lines were compared with control lines, 217 significantly dysregulated mRNAs were identified (111 downregulated and 106 upregulated). The number of transcripts differentially expressed in both HD lines and in each line separately is shown on the Venn diagram (Figures 1E,F). We have focused our experimental validation and further bioinformatics analyses on these mRNAs (Supplementary Table S3). Lists of 30 the most dysregulated transcripts for each group were shown in Tables 2–4.

A total of 17 differentially expressed mRNAs were selected to verify RNA-seq data by quantitative real-time PCR (qPCR) (Figure 2). We have selected 10 mRNAs which were differentially expressed in both HD iPSC lines (OTOGL, TRIM69, CNTNAP3B, MEIOB, C3, PARP12, XDH, CDKN1A, ZFP30, WDR72), 4 mRNAs which were differentially expressed in HD71Q iPSC lines (PIWIL2, HIST1H3C, FAM65B, PDGFB) and 3 mRNAs differentially expressed in HD109Q iPSC line (TP53, PHLDA3, TRIM22). Validation of transcripts dysregulated in both HD lines confirmed results obtained from the RNA-seq analysis for 7 out of 10 analyzed mRNAs. Validation of chosen transcripts dysregulated only in HD71Q lines was completely consistent with up- and downregulation in RNA-seq data. Last experimental analysis, considering mRNAs dysregulated only in HD109Q lines also showed consistency with RNA-seq (Figure 2).

We performed a mass spectrometry analysis to validate dysregulated transcripts at the protein level (Figure 3 and Supplementary Table S4). As a result, we identified 65 differentiating proteins, however, the proteins from chromosomes X and Y were excluded. Among the proteins which demonstrated statistically significant change, there were TP53 and ZFP30 which were also found in the pool of dysregulated transcripts and which showed the same direction of the level change. These proteins are connected with intrinsic apoptotic signaling pathway in response to DNA damage (TP53) and DNA binding transcription factor activity (ZFP30). Moreover, in support of our transcriptional changes, we found that many of the non-statistically significant proteins demonstrated dysregulation pattern similar to significantly dysregulated transcripts.

Using the web interface of CPDB we conducted the interaction- and pathway-centric analysis of list of differentially expressed mRNAs in each comparison group. Induced network modules analysis generates network in which genetic, biochemical and protein interactions between given genes are shown. It also includes genes that are not in the submitted list but connect two or more seed genes with each other and have many connections within the induced network module. Such bioinformatics approach helped us with visualizing biological associations between dysregulated transcripts and their connections with other genes, which may have an impact on HD pathogenesis.

Our first analysis considered mRNAs differentially expressed in both iPSC lines (Figure 4A). Proteins for 25 out of 107 dysregulated transcripts were assigned to the generated network. Main observation from the visualized protein complex is that most of seed genes are associated with each other through intermediate nodes which were not present in the input list. Although no obvious center node can be distinguished, PARK2 seems to be crucial for the whole network as it unites all of the other protein complexes. PARK2 gene is associated with synaptic vesicle exocytosis and central nervous system development. Among other interesting proteins highlighted in the network, there is PIK3R1 whose gene is strongly upregulated in both HD iPSC lines and necessary for the insulin-stimulated increase in glucose uptake, but it is also connected with axon guidance and negative regulation of the apoptotic process. Different proteins in the generated network are also involved in signal transduction and nervous system development, like ARHGAP8, DPYSL4, and FLRT2. Other biological processes in which visualized proteins are involved are RNA splicing and cell motility.

Second analysis considered mRNAs differentially expressed only in HD71Q iPSC lines (Figure 4B). Proteins for 37 out of 198 seed genes were visualized during the analysis. The analysis shows the presence of histones (HIST1H3C, HIST1H2BB, HIST1H1A) and other transcription regulators (ZNF257, BCAS3) in the network. All of the downregulated genes for which proteins are present in the main network are also connected with transcription regulation, like MCM5 and RNF20. Almost all of the upregulated nodes mentioned above are associated with positive regulation of transcription.

The third analysis, considering mRNAs differentially expressed in HD109Q iPSC lines (Figure 4C), shows the crucial role of downregulated transcripts and gene regulatory interactions. Proteins for 41 out of 217 genes from the input list were visualized in the network. The most important protein present in the generated network, with the highest number of edges, is TP53. Many of its interactors are, like TP53 itself, associated with apoptotic signaling, such as tumor necrosis factor receptor superfamily members, TNRFSR10B, TNFRSF10C, and TNFRSF10D but also BAX, CDKN1A, PLK2, HSPA1A, and others.

Over-representation and enrichment analyses were performed with significantly dysregulated mRNAs for each comparison group. We used two bioinformatics tools to identify overrepresented gene ontology and pathway-based terms – CPDB and ClueGO.

In each analysis, we have focused on identifying pathway and gene ontology-based gene sets considering molecular functions, biological processes and cellular components (Figure 5). First analysis included the list of differentially expressed mRNAs in both HD iPSC lines. The most overrepresented GO terms, with the biggest number of genes associated with the term were metal ion binding, regulation of cytokine production and GTPase activator activity (Figure 5A). As for the pathway-based analysis, phospholipase c signaling pathway was distinguished as the one with the lowest p-value (Supplementary Table S5). Phospholipase-C is known to be key signaling proteins in the cellular action of many hormones, neurotransmitters, growth factors, and other extracellular stimuli. The input overlap members in this pathway are PLCB1 and PIK3R1. These two genes are also members of pathways associated with Joubert syndrome, a brain development disorder characterized by the agenesis or underdevelopment of the cerebellar vermis and also by the malformed brain stem. Other identified pathways include selective serotonin reuptake inhibitor pathway, reelin signaling pathway, Beta2 integrin cell surface interactions, downstream signaling of activated FGFR1 and apoptosis.

The most overrepresented GO terms obtained due to the analysis of significantly dysregulated transcripts in HD71Q iPSC line included positive regulation of cellular process, metal ion binding and regulation of cellular component organization (Figure 5B). Cellular components with the lowest p-value and the higher number of submitted transcripts at the same time referred to ruffle, germ plasm and chromatoid body. As for the pathway-based analysis, transcriptional regulation of white adipocyte differentiation is the one with the lowest p-value (Supplementary Table S5). The input overlap members in this pathway are LPL, CEBPA, and ANGPTL4. Other identified pathways are p73 transcription factor network, inflammatory mediator regulation of TRP channels, DNA damage response and lncRNA mediated mechanisms of therapeutic resistance.

Last analysis, which focused on differentially expressed mRNAs in HD109Q iPSC line, revealed the regulation of biological and cellular process and a great number of apoptosis-related terms (Figure 5C). As for the pathway-based analysis, direct p53 effectors and p53 signaling p53 signaling pathway are the most significant among identified pathways (Supplementary Table S5). Other include terms like validated transcriptional targets of TAp63 isoforms, DNA damage response, signaling pathways in glioblastoma, apoptosis modulation and signaling, Wnt signaling pathway and pluripotency and a few viral infection pathways.

Protein complex-based sets of mRNAs, whose protein products are members of the same annotated protein complex, are shown in Supplementary Table S6.

ClueGO overrepresentation and enrichment analysis of transcripts significantly dysregulated in both HD iPSC lines revealed positive regulation of humoral immune response and ion transport by P-type ATPases as GO terms with the lowest p value (Figure 6 and Supplementary Table S7) C3 gene, which is highlighted as associated with the first biological process, is in the top 10 dysregulated transcripts in both HD iPSC lines. Next to C3, there are few other genes that are present in many generated clusters, such as ZP3, KLK5, ATP10B, ATP8B3, and ATP12A.

Analysis of transcripts differentially expressed in HD71Q iPSC lines showed the presence of processes connected with DNA damage response among which CDKN1A, GADD45B, DDB2, and RAD51 play major roles (Figure 6 and Supplementary Table S8). Other annotations included regulation of establishment of cell polarity, PPAR signaling pathway, negative regulation of transforming growth factor beta receptor signaling pathway, Prader-Willi and Angelman Syndrome, removal of licensing factors from origins, mononuclear cell migration and inflammatory mediator regulation of TRP channels.

The last ClueGO analysis, considered transcripts significantly dysregulated in HD109Q iPSC lines, generated a large cluster connected with p53 signaling pathway (Figure 6 and Supplementary Table S9). Six genes associated with these processes are also in the list of the top 30 dysregulated transcripts in HD109Q iPSC lines. Different processes highlighted during the analysis, which were specific to HD109Q iPSC lines, included Pancreatic cancer, interleukin-3, 5 and GM-CSF signaling, Leishmaniasis, Systemic lupus erythematosus and a large cluster associated to TNF-related factors activation by TP53.

The performed meta-analysis showed that among genes significantly dysregulated in HD iPSC lines several are shared with iPS-derived or ES-derived NSCs and/or neurons (Figures 7A–D and Supplementary Table S10). Five out of 107 genes found to be altered in both HD lines overlap with genes from neurons. Two of them, LHFP and FLRT2 are strongly downregulated in both HD lines, with log2FoldChange = -4. Among 16 genes which overlap between HD109Q lines and neurons 9 are strongly downregulated in our work, with TP53 on top of the list. Names of genes overlapping between all three comparison groups were listed in Figure 7D.

The meta-analysis was also performed to establish whether genes altered in our study are also shared with iPSCs obtained in other studies (Figures 7E–G). The analysis revealed that 8, 15, and 14 genes dysregulated in both HD, 71Q and 109Q lines are shared with iPSCs, respectively. Among eight genes which overlap between both HD lines and iPSCs from other studies, ZNF208 and ZNF257 were the most dysregulated in our analyses with log2FoldChange 9,24 and 7,15, respectively. Lists of overlapping genes created during the analysis included also those associated with HD-altered neurodevelopmental pathways, like TGFβ (TGFBI), as well as p53 pathway (CDKN1A, GADD45B), cell adhesion (ANK1) and calcium signaling (CALCRL, ANXA2). Lists of overlapping genes are included in Supplementary Table S10. Only four genes overlap between all HD lines and ESCs, including ACTA1, DPYSL4, HSPA5, and PCK2.