If you don't remember your password, you can reset it by entering your email address and clicking the Reset Password button. You will then receive an email that contains a secure link for resetting your password
If the address matches a valid account an email will be sent to __email__ with instructions for resetting your password
Child Study Center, Yale University School of Medicine, New Haven, ConnecticutProgram in Neurodevelopment and Regeneration, Yale University School of Medicine, New Haven, Connecticut
Child Study Center, Yale University School of Medicine, New Haven, ConnecticutProgram in Neurodevelopment and Regeneration, Yale University School of Medicine, New Haven, Connecticut
Child Study Center, Yale University School of Medicine, New Haven, ConnecticutProgram in Neurodevelopment and Regeneration, Yale University School of Medicine, New Haven, Connecticut
Child Study Center, Yale University School of Medicine, New Haven, ConnecticutDepartment of Psychiatry, Yale University School of Medicine, New Haven, Connecticut
Child Study Center, Yale University School of Medicine, New Haven, ConnecticutProgram in Neurodevelopment and Regeneration, Yale University School of Medicine, New Haven, Connecticut
Program in Neurodevelopment and Regeneration, Yale University School of Medicine, New Haven, ConnecticutDepartment of Pathology, Yale University School of Medicine, New Haven, Connecticut
Yale Kavli Institute for Neuroscience, Yale University School of Medicine, New Haven, ConnecticutDepartment of Neurobiology, Yale University School of Medicine, New Haven, Connecticut
Address correspondence to Flora M. Vaccarino, M.D., Yale University School of Medicine, Child Study Center, 230 South Frontage Road, New Haven, CT 06520
Program in Neurodevelopment and Regeneration, Yale University School of Medicine, New Haven, ConnecticutYale Kavli Institute for Neuroscience, Yale University School of Medicine, New Haven, ConnecticutDepartment of Neurobiology, Yale University School of Medicine, New Haven, Connecticut
Genome-wide association studies have not revealed any risk-conferring common genetic variants in Tourette syndrome (TS), requiring the adoption of alternative approaches to investigate the pathophysiology of this disorder.
Methods
We obtained the basal ganglia transcriptome by RNA sequencing in the caudate and putamen of nine TS and nine matched normal control subjects.
Results
We found 309 downregulated and 822 upregulated genes in the caudate and putamen (striatum) of TS individuals. Using data-driven gene network analysis, we identified 17 gene coexpression modules associated with TS. The top-scoring downregulated module in TS was enriched in striatal interneuron transcripts, which was confirmed by decreased numbers of cholinergic and gamma-aminobutyric acidergic interneurons by immunohistochemistry in the same regions. The top-scoring upregulated module was enriched in immune-related genes, consistent with activation of microglia in patients’ striatum. Genes implicated by copy number variants in TS were enriched in the interneuron module, as well as in a protocadherin module. Module clustering revealed that the interneuron module was correlated with a neuronal metabolism module.
Conclusions
Convergence of differential expression, network analyses, and module clustering, together with copy number variants implicated in TS, strongly implicates disrupted interneuron signaling in the pathophysiology of severe TS and suggests that metabolic alterations may be linked to their death or dysfunction.
Tourette syndrome (TS) is a childhood-onset neuropsychiatric disorder characterized by motor and vocal tics with a population prevalence estimated to be 1% to 3% and frequent comorbidity with obsessive-compulsive disorder and attention-deficit/hyperactivity disorder (
The Tourette Syndrome Association International Consortium for Genetics A complete genome screen in sib pairs affected by Gilles de la Tourette syndrome. The Tourette Syndrome Association International Consortium for Genetics.
). Rare mutations have been identified in two genes (SLITRK1 and HDC), but it remains unclear whether and to what degree these mutations provide insights applicable to the majority of TS cases (
) are more prevalent in TS patients. In a small proportion of cases, this disorder is associated with an infectious event and accompanied by immunoreactive changes in the peripheral blood (
). Postmortem analyses of brain tissue from patients with lifelong TS demonstrated reductions in immunohistochemical markers for cholinergic and parvalbumin+-gamma-aminobutyric acid (GABA)ergic interneurons in the caudate and putamen (
). However, whether interneurons are lost or whether the cells simply downregulate protein marker expression due to some other unrelated mechanism is unknown.
To analyze the neurobiology of TS in an unbiased fashion, we characterized the transcriptome of the caudate and putamen of nine TS individuals with documented unremitted symptoms compared with nine matched normal control (NC) subjects. Differential gene expression and unbiased, data-driven gene coexpression network analyses allowed us to identify genes and gene networks that distinguish brain tissue of individuals with chronic TS from that of unaffected individuals. Our findings provide the first direct evidence for two separate neurobiological factors playing a role in the disorder, decreased signaling of multiple classes of striatal interneurons and dysregulated activation of the immune system.
Methods and Materials
Subjects
TS brain samples of cases with refractory symptoms in adulthood, collected by the Tourette Syndrome Association, were obtained from the Harvard Brain Tissue Resource Center. Samples of normal control subjects were obtained from the Harvard Brain Tissue Resource Center and from the Department of Pathology, Yale University (Table S1). Cases with known or suspected neurological diseases and prolonged hypoxia were excluded. Four of the nine pairs were female subjects. There were no statistical differences in age, postmortem interval, or RNA quality between the TS and control groups (Table S1).
RNA Isolation, Library Construction, and Sequencing
Dissected regions were stored at −80°C until use, when sections were cut at 50 µm. During cryosectioning, prescoring the surface of the block every 150 µm permitted separate collection of the caudate and putamen. Total RNA was isolated using RNeasy Plus Mini kit (Qiagen, Germantown, Maryland) according to the manufacturer’s instruction, using 20 mg of tissue per genomic DNA eliminator column. RNA was assessed using Nanodrop spectrometer (Thermo Scientific, Suwanee, Georgia) and Bioanalyzer (Agilent Technologies, Wilmington, Delaware) to check quantity and quality. Paired-end complementary DNA libraries were prepared introducing a six-base sequence index into the adapter at the polymerase chain reaction (PCR) stage of sample preparation for Illumina multiplexed sequencing, according to the manufacturer’s specification. Samples were sequenced using an Illumina Hi Seq2000 sequencing machine (Illumina, Inc., San Diego, California) to an average depth of 30 million reads per sample.
). We then used coverageBed in BEDtools to process the BED files and MRFquantifier in RSEQtools to process the mapped read format files. This resulted in estimated expression levels for about 36,000 genes.
) to assess differential gene expression (DGE) at a false discovery rate (FDR) <.05. We accounted for matched pair design by using the generalized linear model capability of edgeR. Blocking was applied to individuals in the case of caudate- or putamen-specific differential expression analysis or to both individuals and regions in the case of the analysis of the full dataset. See Supplementary Methods in Supplement 1 for more details.
Co-expression Network Analysis
We used the Weighted Gene Co-expression Network Analysis (WGCNA) package (
). The analysis produced a network of 20 modules, corresponding to about 12,000 genes. Permutation analysis was used to verify the modules’ genes to be coexpressed beyond chance. The WGCNA function userListEnrichment was used to test the modules for overrepresentation in BrainLists, DGE, and copy number variants (CNVs) (see Supplementary Methods in Supplement 1 for details).
Gene Set Enrichment Analysis
We used the Gene Set Enrichment Analysis (GSEA) package (
) to identify the WGCNA modules that were enriched in differentially expressed genes (DEGs). The gene expression matrix was used as input to the package, normalizing the RPKM data by log2 (i.e., log2 [RPKM + .5]) transformation. The network modules were used as test gene sets, according to the format required by the package. The cutoff for significant enrichment was FDR q value <.05.
Functional Enrichment Analysis
Ingenuity Pathway Analysis (IPA) (Ingenuity Systems, Redwood City, California; www.ingenuity.com) and ConsensusPathDB (
) were used to test differentially expressed genes and WGCNA modules for overrepresentation in Gene Ontologies, Canonical Pathways, Upstream Regulators, and Networks.
Sample Clustering
The function heatmap.2 was used to make the heat map in Figure S1 in Supplement 1. We used the Canberra method within the distance function. The DESeq package (
RNA was reverse transcribed using SuperScript III (Life Technologies, Grand Island, New York), and expression levels were evaluated with Applied Biosystems OneStep (Life Technologies) quantitative real-time PCR (qPCR) using SYBR Power Green master mix and a standard protocol that included melting curve analysis. The raw data were normalized to the geometric mean of two reference genes, GAPDH and IPO8. Differential gene expression was determined using the ddCt method (
Stability of endogenous reference genes in postmortem human brains for normalization of quantitative real-time PCR data: Comprehensive evaluation using geNorm, NormFinder, and BestKeeper.
). Genes chosen for validation included randomly selected genes, as well as cholinergic-related and inflammatory/immune response-related genes. Primers are listed in Table S2.
Immunohistochemistry
Formalin-fixed tissue samples from the basal ganglia were sectioned at 50 µm and processed for nitric oxide synthase 1 (NOS1), CD45RO (protein tyrosine phosphatase, receptor type, C [PTPRC]), and ionized calcium-binding adapter molecule 1 (IBA) immunostaining as previously described (
Differentially Expressed Genes in the Basal Ganglia in TS
Whole genome transcripts in the caudate and putamen of nine TS subjects and nine matched normal control subjects (Table S1) were obtained by RNA-Seq and analyzed for differences between patients and control groups using DGE and co-expression network analysis; for an overview, see Figure 1. DGE analysis of the combined caudate and putamen dataset (considering n = 9 matched pairs, blocking for individual pairs and for region within each pair; see Methods and Materials) revealed 1131 genes differentially expressed between TS and NC, of which 309 were downregulated and 822 were upregulated (FDR < .05; Table S3a). Unsupervised hierarchical clustering of the top 100 differentially expressed genes revealed co-regulated patterns of gene expression across samples with segregation of TS and NC, independent of postmortem interval, RNA quality, medication use, and other confounding variables (Figure S1 in Supplement 1). Only a single TS subject (TS4) more closely resembled the NCs in the unsupervised hierarchical clustering analysis. This subject had no first-degree relatives with TS but did have a positive history for a number of perinatal risk factors. This observation again provides clear evidence for heterogeneity in the pathobiology and etiology of TS. Validation by qPCR showed the same direction fold change for 13 of 13 genes examined with an overall correlation of .88 between qPCR and DGE; p value: 6.4 × 10−5 (Figure S2 in Supplement 1).
Figure 1General scheme of our analyses. Analysis flow diagram for RNA-Seq data analysis. Boxes: preanalysis and postanalysis data; arrows: methods applied to preanalysis data to get postanalysis data/results. Red text: results. RNA-Seq reads are converted into read counts and the reads per kilobase of transcript per million mapped reads, or RPKM, for each transcript in the Gencode annotation. EdgeR is applied to counts to infer differential gene expression (DGE). Weighted Gene Co-expression Network Analysis, or WGCNA, in combination with Gene Set Enrichment Analysis, or GSEA, is applied to RPKM to infer modules of coexpressed genes enriched in differentially expressed genes. FDR, false discovery rate; immuno, immunohistochemistry; qPCR, quantitative real-time polymerase chain reaction.
Downregulation of Transcripts Involved in Neuronal Signaling, Specifically Striatal Interneurons
Functional annotation of the downregulated genes by IPA revealed a signature of neuronal signaling dysfunction in the striatum of TS individuals (Tables S4 and S5). In the striatum, projection neurons of the direct and indirect pathways comprise at least 70% to 80% of all neurons (
). The striatum also has four types of interneurons: the tonically active cholinergic neurons, and three classes of GABA interneurons, expressing NOS1/neuropeptide Y (NPY)/somatostatin (SST), parvalbumin (PVALB), or CALRETININ (CALB2 or CR) (
). To find out whether transcript downregulation was specific for particular neuron types, we performed manual classification of the DGE according to neurotransmitter class (Table 1). Genes encoding for abundant neuropeptides expressed by striatal projection neurons, such as prodynorphin and tachykinin precursor 1, were not changed, whereas there was significant downregulation for genes related to cholinergic interneurons, including acetylcholine neurotransmitter synthesis and transport (CHAT, SLC5A7, SLC18A3), cholinergic receptors (CHRM2), and the homeodomain genes LHX8 and GBX2, critical for the development of forebrain cholinergic neurons (
) (Table 1). Several genes related to striatal GABAergic interneurons also showed reduced expression, including NOS1, NPY, and SST but not PVALB, CALB2/CR (Table 1).
Table 1Hand-Curated Subset of Differentially Expressed Genes
FC
FDR
Cholinergic
CHAT
Choline O-acetyltransferase
−2.31
1.39E-07
CHRM2
Cholinergic receptor, muscarinic 2
−1.99
4.32E-04
SLC18A3
Solute carrier family 18 (vesicular acetylcholine), member 3
−2.17
4.04E-05
SLC5A7
CHT1; solute carrier family 5 (choline transporter), member 7
−2.55
6.65E-08
LHX8
LIM homeobox 8
−1.90
1.39E-03
NTRK1
Neurotrophic tyrosine kinase, receptor, type 1
−1.91
1.17E-03
GBX2
Gastrulation brain homeobox 2
−1.65
4.23E-02
GABAergic
GAD1
Glutamate decarboxylase 1 (brain, 67 kDa)
−1.64
1.21E-02
GABRA1
Gamma-aminobutyric acid A receptor, alpha 1
−1.47
2.31E-02
GABRA3
Gamma-aminobutyric acid A receptor, alpha 3
−1.96
5.85E-04
GABRG2
Gamma-aminobutyric acid A receptor, gamma 2
−1.50
2.61E-02
GABRQ
Gamma-aminobutyric acid A receptor, theta
−1.82
1.11E-02
NOS1
Nitric oxide synthase 1 (neuronal)
−1.58
6.40E-03
NPY
Neuropeptide Y
−1.66
8.15E-03
NPY2R
Neuropeptide Y receptor Y2
−2.81
2.33E-04
SST
Somatostatin
−1.59
4.42E-03
Dopaminergic
DRD1
Dopamine receptor D1
−1.50
3.20E-02
DRD5
Dopamine receptor D5
−2.10
2.47E-03
Others of Interest
SLITRK5
SLIT and NTRK-like family, member 5
−1.54
1.08E-02
OPRK1
Opioid receptor, kappa 1
−1.60
3.48E-02
OPRM1
Opioid receptor, mu 1
−1.43
4.03E-02
HTR2C
5-hydroxytryptamine (serotonin) receptor 2C, G protein-coupled
−1.62
7.44E-03
Genes are listed under each category they represent, along with fold change and FDR corrected p values.
Previously, we reported decreased numbers of cells immunostained for CHAT and PVALB in the striatum of TS; however, NOS/NPY/SST interneurons were not assessed (
). Hence, we quantified the density of immunohistochemically identified NOS1+ cells in series of postmortem tissue sections from seven TS and seven control cases, three of which were new cases and four that had been also evaluated using RNA-Seq. We found a significant decrease in NOS1+ interneurons both in the caudate (−37.6%, p = 1.25E-2) and in the putamen (−42.2%, p = 1.24E-2) of TS versus control subjects (Figure 2A–C), confirming the transcriptome analyses. No significant change in PVALB transcripts was detected in the transcriptome. However, we detected a 75% decrease in messenger RNA for potassium voltage-gated channel subfamily C member 1 (KV3.1) (KCNC1) approaching statistical significance (FDR = 6.8E-2). KCNC1 is a potassium channel selectively expressed by PVALB interneurons that confers fast re-polarization and high firing properties typical of these cells (
Developmental expression and functional characterization of the potassium-channel subunit Kv3.1b in parvalbumin-containing interneurons of the rat hippocampus.
). Additional GABAergic signaling genes downregulated in TS included the GABA synthetic enzyme glutamate decarboxylase 1 (GAD1) and several GABAergic receptor subunits (Table 1).
Figure 2Differential nitric oxide synthase 1 (NOS+) cell density, elevated protein tyrosine phosphatase, receptor type, C (CD45/PTPRC) cell density, and microglial activation in Tourette syndrome (TS) patients. (A–C) Immunostaining for NOS+ cells in the caudate of representative normal control (NC) (A) and TS cases (B). Scale bars = 100 µm. (C) Stereological estimates of NOS1+ cell density (cells per mm3) in the caudate and putamen of normal control (n = 7) and TS (n = 7) cases. (D–F) Immunostaining for CD45/PTPRC+ cells in the caudate of representative normal control (D) and Tourette syndrome cases (E). Scale bars = 100 µm. (F) Stereological estimates of PTPRC+ cell density (cells per mm3) in the caudate of control (n = 10) and TS (n = 8) cases. (G–I) Immunostaining for ionized calcium-binding adapter molecule 1 (IBA+) cells in caudate of representative normal control (G) and Tourette syndrome cases (H). Scale bars = 100 µm. (I) IBA+ cell body size (µm2) and dendrite length (µm) in the caudate in control (n = 7) and TS (n = 5) cases. *p < .05, **p < .01 one-tailed Student t test. Cd, caudate; Pt, putamen.
Upregulation of Transcripts and Proteins Related to the Monocyte/Microglia Lineage
Functional annotation of the genes upregulated in TS by IPA revealed inflammatory response and immune cell development and trafficking signatures (Tables S4 and S5). Upregulated transcripts included tumor necrosis factor alpha, interleukin (IL)-6, and IL-12, previously found altered in TS in peripheral studies (
); toll-like receptors that modulate regulatory T cell function; several major histocompatibility complex class II molecules; immunoglobulin chains; and several chemokine (C-C motif) receptor tumor necrosis factor-, and IL-family cytokines that have been implicated in autoimmune disorders. Pathway analysis conducted with MetaCore (Thomson Reuters, New York, New York) yielded results similar to IPA for both upregulated and downregulated transcripts (Table S6).
We considered the possibility that an infiltrate of immune cells from the periphery might explain the upregulation of immune transcripts in brain. However, the transcriptome data did not show a significant upregulation of B cell and T cell specific transcripts, whereas the messenger RNA for phosphotyrosine phosphatase C (CD45/PTPRC), a transcript expressed by all differentiated hematopoietic cells including lymphocytes, monocytes, and macrophages (
), was increased by approximately twofold in TS brains (Table S3c). We examined PTPRC+ cells in the caudate of TS and control brains by immunohistochemistry. We found a significant increase in density of PTPRC+ cells in TS (+153%, p = 7.30E-3) (Figure 2D–F). As several microglia-related transcripts (IBA, CD68) were upregulated in TS subjects, we considered that increased number of cells of the macrophage/microglia lineage might be underlying the inflammatory changes. Evaluation of IBA protein expression in the caudate revealed no change in the density of IBA+ cells (data not shown). However, we found a significant decrease in IBA+ cell body size (p = 1.27E-2) and process length (p = 2.53E-2), indicative of microglia activation (Figure 2G–I). Together, the results suggest that the primary source of inflammation may be activated brain microglia, although we cannot exclude the involvement of other cells, i.e., peripheral monocytes.
We next intersected our 1131 DEGs with existing microarray transcriptome studies, three in blood (
GABA- and acetylcholine-related gene expression in blood correlate with tic severity and microarray evidence for alternative splicing in Tourette syndrome: A pilot study.
). We found that not much overlap exists between our DGE discovered by RNA-Seq and previous microarray studies in TS tissue. Of the 10 genes found to be nominally differentially expressed in the putamen of three TS brains in the Hong et al. study (
) (see their Table 4), only PTP4A3 was significantly changed in our dataset, and of the 13 genes found to be nominally differentially expressed in the TS blood by exon array (
) (see their Supplemental Table 1 and Supplemental Table 2), only one gene, MYO7A, intersected with our dataset. None of the 51 genes nominally differentially expressed in the second blood study (
GABA- and acetylcholine-related gene expression in blood correlate with tic severity and microarray evidence for alternative splicing in Tourette syndrome: A pilot study.
) (see their Supplemental Table 1) intersected with our dataset. Finally, of the 14 genes significantly differentially expressed in the third TS blood study (
) (see their Table 1), only one gene, LILRA3, intersected with our 1131 DEGs.
Differentially Expressed Genes in the Caudate and Putamen in TS
In addition to performing a combined analysis, we also analyzed separately the caudate and putamen to identify possible region-specific contributions to the overall signature. We found 434 DEGs in the caudate and 368 in the putamen (FDR < .05; Table S3b and c). Intersection of the DEGs between the caudate and putamen (see Venn Diagram in Figure S3 in Supplement 1) revealed 198 DEGs specific to the caudate and 132 specific to the putamen (Table S7a), hinting at a possible differential regional involvement of these regions in the disease. Functional analysis of the region-specific lists (listed in Table S7a) showed strong involvement of the caudate-specific DEGs in immunological disease, inflammatory response, cell-to-cell signaling and interaction (Table S7b) and a moderate involvement of the putamen-specific DEGs in neurological disorders, cell-to-cell signaling/interaction, and tissue morphology (Table S7c). However, intersection of the DEGs from the caudate, putamen, and combined analysis (see Venn Diagram in Figure S3 in Supplement 1) revealed that analyzing the caudate and putamen in isolation, as compared with the overall combined analysis, identifies, respectively, only 18 and 27 DEGs (Figure S3 in Supplement 1), suggesting that the combined analysis captures most of the individual region DEGs and, as expected, produces about 600 to 700 more DEGs as a consequence of the increased power, thus allowing more discovery. In conclusion, while the difference in DEGs between caudate and putamen may be a true biological signature, we cannot rule out the possibility of a statistical artifact due to lower number of samples. We therefore focus on the combined dataset throughout the rest of this article.
Co-expression Network Analysis
Preliminary analyses suggested that the level of gene expression for CHAT was strongly correlated with that of other cholinergic genes (Pearson’s mean r = .71), as well as with GABAergic genes (Pearson’s mean r = .74), including GAD1, several GABA receptor subunits, NOS1, NPY, and SST, across all TS subject and control subjects. In contrast, the correlation between CHAT and immune system-related transcripts displayed a poor correlation coefficient (Pearson’s mean r = .16) (Figure S4 in Supplement 1).
We next used network inference and analysis to explore in an unbiased fashion correlations among transcripts that may provide clues to disease pathogenesis. To this end, we applied WGCNA (
) to the entire transcriptome dataset (see Figure 1 for an overview of the analyses). This analysis yielded 20 modules where transcripts in each module were significantly correlated in overall pattern of expression across all samples. After permutation analysis, 17 modules of co-expressed genes remained significant (FDR < .05; Table S8).
To test whether the WGCNA modules were differentially expressed between TS and NC, we performed GSEA (
), where the modules were used as test gene sets and gene expression levels were used as input. Ten modules were significantly enriched in upregulated or downregulated genes (Table 2; Table S9). Their member genes, together with their fold change in expression in TS versus control subjects from DGE analysis, are listed in Table S10. Two modules, magenta and turquoise, were the top scoring upregulated and downregulated modules by GSEA normalized enrichment score (Table S9). As expected, the magenta and turquoise modules showed the strongest statistical overlap with the upregulated and downregulated genes resulting from DGE analysis: 297 of the 822 upregulated genes were in the magenta module (p = 2.5E-318), while 116 of the 309 downregulated genes were in the turquoise module (p = 3.4E-43) (Table S11).
Table 2Biological Classification for the Modules Significantly Enriched in Differentially Expressed Genes
Module Name/Downregulated/Upregulated
Biological Classification
Database
Corrected p Values
Turquoise (Down)
Neuron_probable__Cahoy
BrainLists
3.33E-267
Synaptic transmission
GO:Biofunctions
2.96E-81
Glutamatergic synapse – Homo sapiens
CP:KEGG
2.87E-14
Blue (Down)
blue_M16_Neuron__CTX
BrainLists
6.71E-79
Cellular metabolic process
GO:Biofunctions
6.62E-21
Ubiquitin mediated proteolysis – Homo sapiens
CP:KEGG
8.23E-05
Pink (Down)
DownWithAlzheimers_Liang__ADvsCT_inCA1
BrainLists
1.42E-12
Nucleic acid metabolic process
GO:Biofunctions
5.53E-26
Gene Expression
CP:Reactome
1.35E-17
Green (Down)
turquoise_M14_Nucleus__HumanMeta
BrainLists
8.90E-30
Cellular macromolecule metabolic process
GO:Biofunctions
2.72E-10
Spliceosome – Homo sapiens
CP:KEGG
2.83E-06
Brown (Down)
brown_M3_Astrocytes__HumanMeta
BrainLists
2.79E-03
Nucleic acid metabolic process
GO:Biofunctions
4.63E-27
Gene Expression
CP:Reactome
1.45E-08
Lightcyan (Down)
Astrocyte_probable__Cahoy
BrainLists
1.43E-02
Homophilic cell adhesion
GO:Biofunctions
2.80E-20
Signaling by NOTCH
CP:Reactome
9.37E-04
Tan (Down)
blue_M2_Oligodendrocytes__HumanMeta
BrainLists
9.82E-127
Ensheathment of neurons
GO:Biofunctions
3.86E-08
Metabolism of lipids and lipoproteins
CP:Reactome
2.76E-02
Red (down)
NA
BrainLists
NA
Nucleobase-containing compound biosynthetic process
GO:Biofunctions
8.12E-15
Generic Transcription Pathway
CP:Reactome
8.51E-11
Magenta (Up)
pink_M10_Microglia(Type1)__HumanMeta
BrainLists
1.29E-107
Immune response
GO:Biofunctions
5.40E-79
Staphylococcus aureus infection - Homo sapiens
CP:KEGG
2.02E-29
Purple (Up)
Astrocyte_probable__Cahoy
BrainLists
2.90E-49
Cell adhesion
GO:Biofunctions
4.43E-06
Propanoate metabolism – Homo sapiens
CP:KEGG
1.11E-02
Modules are shown in decreasing order of GSEA normalized enrichment score (see Table S8) for each of the downregulated or upregulated groups. The Biological Classification column lists the most significant lists from overrepresentation analysis using list categories shown in the Database column. Their sources are WGCNA for BrainLists and ConsensusPathDB (
We next examined the biological role of the 10 differentially expressed modules by enrichment analysis for brain-related cell types and diseases using the BrainLists database from WGCNA, as well as for Gene Ontology and Canonical Pathways categories using GO:Biofunctions, CP:KEGG, and CP:Reactome from ConsensusPathDB (
). The modules’ biological classifications are summarized in Table 2 and described in detail in Tables S12 and S13. There were two neuronal modules, turquoise and blue, which were the two top scoring downregulated modules in TS. The turquoise module was enriched in PVALB cell type categories and various neurotransmission-related biofunctions (i.e., glutamatergic, dopaminergic, GABAergic, and cholinergic synapses, calcium and potassium channels); indeed, the majority of the striatal interneuron-related transcripts in Table 1 are members of the turquoise module. The blue module was enriched in cellular catabolic pathways and protein ubiquitination. The magenta module, the top scoring for upregulated genes, was enriched in microglia cell type categories and immune-related pathways. There were two astrocyte modules, the purple and the lightcyan; the first was enriched in cell adhesion and astrocyte-related metabolic pathways, whereas the second was enriched in protocadherin-gamma molecules (Table 2 and Tables S12 and S13).
To understand whether the DGE genes impact specific subgroups or more connected genes within modules, we estimated the gene scaled intramodular connectivity (
) and on this basis determined the 10 genes that were most connected to all the other genes within each module (hub genes) (shown in the center of each module in Figure 3). We found that the microglia/immune module (magenta) is perturbed dramatically at gene expression level, with many highly interconnected genes, including the hub genes, being overexpressed (Figure 3, red color). Conversely, the gene expression changes in the turquoise, blue, lightcyan, and purple modules involved nonhub, less interconnected genes (Figure 3, red and blue colors for overexpressed and underexpressed genes, respectively).
Figure 3Differentially coexpressed gene networks in Tourette syndrome. The center of the figure shows the global network where the nodes (green circles) are Weighted Gene Co-expression Network Analysis modules and the edges are significant correlations between the modules eigengenes. The four inserts show the detailed network within four relevant modules. Top left: blue (neuronal metabolism) module. Top right: magenta (immune/microglia) module. Bottom right: lightcyan (cell adhesion/protocadherin) module. Bottom left: turquoise module (synaptic transmission/interneuron). Circles: network genes. Large circles: top five genes with highest intramodular connectivity (hubs). Squares: network genes overlapping with copy number variant genes. Red: upregulated gene. Blue: downregulated gene. Concentric circle arrangement in the case of the turquoise (bottom left) and magenta (top right) modules reflects intramodular connectivity: genes closer to the center have higher intramodular connectivity than genes further away from the center.
Lack of Correlation Between Interneuron-Related and Immune-Related Gene Expression
To reveal how modules interact within the network, we first computed the eigengene, which is the first principal component of a module’s member genes expression levels. We then applied hierarchical clustering to the 17 modules’ eigengenes (Figure S5 in Supplement 1) and found that the two astrocytes modules, lightcyan and purple, formed a distinct group. Furthermore, the two neuronal modules, turquoise and blue, enriched in interneuron and neuronal catabolism biofunctions, respectively, clustered together, whereas the microglial module (magenta) significantly correlated with the nucleus (green), transcription (red), and oligodendrocyte (tan) modules but was quite independent from the neuronal and astrocyte modules. This confirms the previous suggestion that the immune response may be decoupled from the changes in interneurons. The overall network between the modules, where the nodes are modules and edges are the significant correlations, is represented in the center of Figure 3.
Significant Overlap of Genes Revealed by Co-expressed Network Analysis With Previous TS Genomic Study
To gain further insights into the present findings on the transcriptome, we overlapped our DGE and WGCNA findings with the list of genes intersecting rare CNVs from a family-based TS study (
). We first intersected the CNV genes from the TS group with the ones from the NC group, identified the common genes, and generated a list of TS-specific, NC-specific, and common CNV genes, which we then used for subsequent analysis. We found that only the TS-specific CNV genes, and none of the common or NC-specific genes, overlapped with our datasets. Remarkably, only the turquoise (interneurons) and the lightcyan (protocadherin-gamma/cell adhesion) modules significantly overlapped with TS-specific CNV genes (p < 3.1E-3 and p < 3.1E-8, respectively), whereas no enrichment was found for any other module, including the magenta, which was the top-scoring module for upregulated genes. The TS CNV genes in the lightcyan module intersected the gamma-protocadherin hub genes (Figure 3, squares). In contrast, CNVs in the turquoise module involved more peripheral, less interconnected genes (Figure 3, squares); however, differently than the lightcyan, these CNVs genes included many genes that were differentially expressed (blue and red squares), suggesting significant convergence of DGE, gene network, and CNV analyses in the turquoise module.
Discussion
Here, we performed the first unbiased and comprehensive characterization of gene expression changes in the basal ganglia of patients with TS. We identified 1131 differentially expressed genes in the striatum of individuals with lifelong TS. In comparison, fewer differentially expressed genes were found when caudate and putamen were analyzed separately, and even fewer were specific to each region, suggesting that biological changes are mostly shared across the striatum. We were unable to confirm that transcripts previously found to be nominally differentially expressed in the blood in a subset of TS subjects were differentially expressed in the basal ganglia of the nine subjects included in this study. The observed lack of overlap with previous studies in blood further stresses the need to analyze brain samples, as opposed to blood samples (
GABA- and acetylcholine-related gene expression in blood correlate with tic severity and microarray evidence for alternative splicing in Tourette syndrome: A pilot study.
). There was only one gene overlapping with a previous microarray study of the TS putamen, where 10 genes were nominally differentially expressed in three TS subjects (
), possibly reflecting the low sample size of their study. Interestingly, we obtained a much higher yield in differential gene expression despite a still limited number of brain samples (n = 9 matched pairs, two regions each). This reflects, in part, the careful sample selection and experimental design and, in part, the increased power of RNA-Seq over microarray analyses. Co-expression network analysis, a well-known and established data reduction technique, also addresses power issues, by reducing the dimensionality of the problem from tens of thousands of genes to a few modules of coexpressed genes. The fact that our results are consistent across different, yet complementary, data analytical approaches (i.e., DGE and WGCNA) supports the idea that the number of samples in our study is sufficient to infer differences between the Tourette and NC. Hierarchical clustering revealed segregation of gene expression differences according to diagnosis but not drug treatment, comorbidities, postmortem interval, and RNA quality. Nevertheless, the small sample size limited our ability to control for these variables and future studies will need to consider these potentially confounding factors. Notably, our correlation between qPCR and RNA-Seq was close to .9, in good agreement to previous comparative reports of the two methods, wherein correlation analysis and sometimes additional validation subjects were employed due to the intrinsic higher variability of qPCR assays (
WGCNA identified eight downregulated modules and two upregulated modules in TS. The two top-scoring downregulated, turquoise and blue, were characterized, respectively, by interneuron signaling and neuronal catabolic pathways. The two upregulated modules, magenta and purple, were characterized, respectively, by microglia signaling and astrocyte metabolic pathways (Table 2). The fact that the turquoise module, enriched for genes critical for neurotransmitter production, signaling, and function of cholinergic and GABAergic interneurons, was the most significantly enriched in downregulated genes is consistent with immunohistochemical findings suggesting that striatal interneuron proteins are decreased in TS. While two of these classes of interneurons were previously reported to be decreased in TS based on antibody labeling of one or two identifying proteins (
), the current data show that multiple genes involved in the comprehensive function of these cells are decreased, providing the first confirmation that indeed cholinergic and GABAergic interneurons may be functionally immature or missing in TS. Although GABA interneurons represent less than 4% of the total neuron population in the striatum (
Immunocytochemical localization of choline acetyltransferase within the rat neostriatum: A correlated light and electron microscopic study of cholinergic neurons and synapses.
), they exert a prominent role in feed-forward inhibition driven by corticostriatal input, which predicts that their reduction can result in severe disruptions to cortically driven internal striatal activity and ultimate circuit output (
Unexpectedly, a module greatly enriched in inflammatory response transcripts (magenta) was the most significantly enriched in upregulated genes. This finding is intriguing, given previous reports describing peripheral immunological changes in TS patients (
Antibodies against neural, nuclear, cytoskeletal, and streptococcal epitopes in children and adults with Tourette’s syndrome, Sydenham’s chorea, and autoimmune disorders.
). It is important to note that the cases used in this study did not meet criteria for pediatric autoimmune neuropsychiatric disorders associated with streptococcal infection (
Identification of children with pediatric autoimmune neuropsychiatric disorders associated with streptococcal infections by a marker associated with rheumatic fever.
). We did not find any evidence for a white cell infiltrate in brain tissue, but rather, the changes were consistent with a macrophage/monocyte mediated inflammatory response, supported by an increase in cells positive for the generic marker PTPRC and local reactive microglia within the brain.
A possible hypothesis linking the above findings is that an intrinsically upregulated immune response in striatal tissue is the cause of the neuronal dysfunction, either because interneurons may be directly killed by the immune cells or because of secondary inflammatory alterations. Another possibility is that immune activation would be the consequence of neuronal dysfunction, i.e., secondary to interneuron death. Against both hypotheses is the fact that the transcripts in the magenta microglial module were not correlated with those in the turquoise interneuron module, suggesting that the immune and neuronal events were not causally related to each other. Instead, the magenta module was most correlated with two modules with “nucleus” and “transcription” functional annotations, whereas the turquoise neuronal module was most correlated with the other neuronal module (blue), functionally annotated by catabolic pathways and protein ubiquitination suggesting that neuronal demise may be secondary to metabolic alterations.
Interestingly, only the turquoise interneuron module and lightcyan protocadherin module (but not the magenta immune module) were significantly enriched in genes implicated in the etiology of TS by the largest CNV study published in TS patients to date (
). These two modules, however, are not correlated with each other: the turquoise is most correlated with the blue neuronal catabolism module, as mentioned above, while the lightcyan module is most correlated with the other astrocyte module (purple), which is enriched in organic acid metabolism as well as cell adhesion (Table S13), suggesting a unique role for astrocytes in TS pathogenesis. Furthermore, differently than in the lightcyan module, the CNV genes in the turquoise module include many differentially expressed neurotransmitter-related transcripts (i.e., CHAT, several GABA receptors, and the high affinity choline transporter). This convergence again emphasizes that genes in the turquoise module may provide fundamental clues in TS etiology/pathogenesis. The microglia/immune module (magenta) had widespread perturbation in gene expression levels, including many hub gens, yet lacked overlap with CNV genes. This may reflect a downstream effect of the disease, although it cannot be excluded that any undetected genetic defect may be present at the level of the regulators of transcripts within this module. Overall, the lack of correlation between the magenta and the other differentially expressed modules that are enriched in genes overlapping with CNVs suggests that immune and neuronal events are not pathophysiologically dependent on each other; however, we cannot exclude that a previous relationship might have existed developmentally. Specifically, it is worth noting that microglia are intimately involved in circuit formation during central nervous system development and circuit refinement (
In summary, a convergence of results from DGE and network analyses of the present transcriptome data, together with CNVs previously implicated in TS, strongly implicates disrupted basal ganglia neuronal signaling, particularly interneurons (i.e., cholinergic and NOS1+-GABAergic populations) in the pathophysiology of severe, persistent TS and suggests that metabolic alterations may be linked to their death or dysfunction. Concomitantly, but decoupled from the neuronal changes, we found a significant increase in immune and inflammatory transcripts that appear to be endogenous to the central nervous system. Developmentally relevant modeling may help us to untangle the cause-effect mechanisms underlying these phenomena.
Acknowledgments and Disclosures
This work was supported by grants from the National Tourette Syndrome Association, the National Institutes of Health (Grant Nos. NS054994 to FMV and MH081896 to NS), and from the Brain and Behavior Research Foundation (National Alliance for Research on Schizophrenia and Depression) to FMV. YK-S is currently affiliated with the Department of Neural Regenerative Medicine, Sapporo Medical University School of Medicine, Sapporo, Japan. DP is currently affiliated with the Institute of Mathematics and Informatics, Bulgarian Academy of Sciences, Sofia, Bulgaria.
We thank Yuka Imamura-Kawasawa for her technical advice on RNA isolation; the National Tourette Syndrome Association, George Tejada, Louis Fernandes, and the Harvard Brain Tissue Resource Center for providing Tourette and control tissue; and Nancy Thompson for help with the clinical phenotyping. The Harvard Brain Tissue Resource Center is supported in part by National Institutes of Health Grant R24 MH068855.
Dr. Leckman has received support from the following: National Institutes of Health (salary and research funding), Tourette Syndrome, Grifols (research funding), Klingenstein Third Generation Foundation (medical student fellowship program), John Wiley and Sons (book royalties), McGraw Hill (book royalties), and Oxford University Press (book royalties). All other authors report no biomedical financial interests or potential conflicts of interest.
The Tourette Syndrome Association International Consortium for Genetics
A complete genome screen in sib pairs affected by Gilles de la Tourette syndrome. The Tourette Syndrome Association International Consortium for Genetics.
Stability of endogenous reference genes in postmortem human brains for normalization of quantitative real-time PCR data: Comprehensive evaluation using geNorm, NormFinder, and BestKeeper.
Developmental expression and functional characterization of the potassium-channel subunit Kv3.1b in parvalbumin-containing interneurons of the rat hippocampus.
GABA- and acetylcholine-related gene expression in blood correlate with tic severity and microarray evidence for alternative splicing in Tourette syndrome: A pilot study.
Immunocytochemical localization of choline acetyltransferase within the rat neostriatum: A correlated light and electron microscopic study of cholinergic neurons and synapses.
Antibodies against neural, nuclear, cytoskeletal, and streptococcal epitopes in children and adults with Tourette’s syndrome, Sydenham’s chorea, and autoimmune disorders.
Identification of children with pediatric autoimmune neuropsychiatric disorders associated with streptococcal infections by a marker associated with rheumatic fever.
00551-4&id=gr1.jpg){kind=link}
00551-4&id=gr2.jpg){kind=link}
00551-4&id=gr3.jpg){kind=link}