Transcriptome- and Epigenome-Wide Association Studies of Tic Spectrum Disorder in Discordant Monozygotic Twins

1. Introduction

Tourette syndrome (TS) is a childhood-onset multifactorial disorder characterized by motor and vocal tics [1,2,3]. Chronic tic disorder (CTD) differs from TS in that patients have either vocal or motor tics. As no clinical or genetic evidence exists to suggest that TS and CTD are distinct disorders, the term tic spectrum (TSD) disorder has been proposed to cover both [4].

Initial symptom presentation typically occurs between the ages of 3 and 8 years, starting with simple motor tics followed by phonic tics. Symptom severity commonly peaks at ages 8 to 12 years, and by adulthood most patients will experience symptom attenuation or even complete remission [1,3]. In a recent longitudinal follow-up study of 314 individuals with TSD, only 18% of patients experienced a complete remission after the age of 16 years, 60% still experienced some mild to moderate tics, and the remaining 23% had moderate to severe tics [5]. The reported prevalence of TSD varies from 0.1-6% with a male-to-female ratio of 4:1 [3,6]. Most individuals with TSD have comorbidities, of which the most frequent are obsessive-compulsive disorder (OCD) and attention-deficit hyperactivity disorder (ADHD), with as many as 70% of individuals affected by TSD having either disorder. Other common comorbidities include autism spectrum disorder (ASD), anxiety, depression, sleep disorders, migraine, and self-injurious behavior [1,3].

Large population-based and genome-wide association studies have suggested a heritability of TSD between 0.21-0.77 [7,8], with most of the genetic heritability explained by evolutionarily conserved SNPs [8]. It is largely accepted that TSD has a multifactorial etiology where genetic and environmental factors interact leading to neurobiological vulnerability [9]. Despite research advances, many gaps in knowledge about TSD etiology remain, including but not limited to understanding the environmental contribution to the development of the disorder. In this study, we investigated gene expression and DNA methylation data from eleven monozygotic twins either discordant or concordant for TSD to assess how environmental components contribute to altered gene regulation and the development of the disorder.

2. Materials and Methods

2.1. Monozygotic Twin TSD Cohort

Cross-linking of The Danish Twin Register and The Danish Psychiatric Central Register with subsequent genetic testing and interview-based diagnosis using DSM-5 led to the identification of fourteen monozygotic twin pairs where at least one twin had TSD, as previously described [10], resulting in a final cohort consisting of five twins discordant for TSD, four twins concordant for TSD, and two asymptomatic twins (Table S1 and Figure S1).

2.2. RNA Sequencing

Total RNA was isolated from PAXgene tubes. RNA sequencing was done using TruSeq Stranded Total RNA Library Prep Kit (Illumina, San Diego, CA, USA) and libraries were paired-end sequenced (2x125 bp) on a HiSeq 2500 sequencing system (Illumina).

Reads were mapped to GRCh38.p13 (GCA_000001405.28) using STAR [11]. The resulting count data served as input for downstream pre-processing and differential expression analysis using the edgeR [12] and limma [13] R packages. Raw counts were transformed into counts per million (CPM) and log-transformed (LCPM), genes with low expression were filtered using the filterByExpr function, and samples were normalized using the calcNormFactors function from the edgeR package with default settings. Heteroscedasticity was accounted for using voom from limma.

2.3. DNA Methylation

Genomic DNA was extracted from peripheral blood and subsequently subjected to bisulfite treatment using standard protocols. Bisulfite-converted DNA was analyzed on Infinium MethlationEPIC v1 arrays (Illumina) according to the manufacturer’s instructions. The data was analyzed using the R package minfi [14]. Exclusion criteria for probes were detection p-values > 0.01, SNPs at CpG sites, known cross-reactivity, and location on sex chromosomes. Samples were normalized using a stratified quantile normalization. Counts were converted to β values (β = M/M + U) and M-values (log2(M/U)) where M is the methylated intensity and U is the unmethylated intensity. β values were used for methylation quantification, and M-values were used for statistical analyses.

2.4. Identification of Differentially Expressed Genes and Differentially Methylated Probes

For both gene expression (RNA sequencing) and DNA methylation data, covariates were identified using principal component analysis (PCA), and quantile-quantile plots were used to assess differences between TSD and non-TSD individuals as well as to assess deviations from expected data trends and thus potential needs for model adjustment. To adjust for unwanted variation in the data such as cell type composition, batch effects, and library size, surrogate variables were identified and estimated using the sva R package [15].

Linear models were fitted using limma to assess differences between TSD and non-TSD individuals and to identify associations with Yale Global Tic Severity Scale (YGTSS) scores, incorporating monozygotic twin pair and surrogate variables as covariates. Correction for multiple testing was performed using the Benjamini-Hochberg method [16]. Two models were fitted for RNA-seq data; one compared TSD cases and controls, and the other tested the association of gene expression levels with YGTSS scores for all samples with YGTSS scores > 0, including concordant twins and asymptomatic controls. For the DNA methylation data, a single model was fitted to compare TSD cases and controls.

2.5. Expression Quantitative Trait Methylation (eQTM)

Expression quantitative trait methylation (eQTM) analysis was performed using the R package MatrixEQTL [17]. Gene expression input from RNA sequencing was represented as log2-transformed CPM values with an unadjusted p-value threshold of < 0.05. For DNA methylation, M-values for probes with unadjusted p-values < 0.005 were used. The analysis used a simple t-test and a linear (additive) model adjusted for monozygotic twin pair status. The p-value threshold for gene-CpG pairs was set to 0.05, and the analysis was run in cis with a maximum distance of 1x10⁶ base pairs.

2.6. Gene Set Enrichment Analysis (GSEA)

Gene set enrichment analysis (GSEA) was conducted on the results from the differential expression analysis using the clusterProfiler R package [18] to identify Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) terms associated with TSD. Log2-transformed fold changes was used as input. The geneSet size range was set to 3-800, all subontologies were included, p-value threshold was set to 0.05, correction for multiple testing was carried out using the Benjamini-Hochberg method, the seed parameter was set to TRUE, and the fast gene enrichment analysis was applied [19].

3. Results

3.1. Data Exploration and Statistical Model Generation

Principal component analyses (PCA) were carried out to explore technical variability in the RNA sequencing and DNA methylation data (Figures S2A and S3A), and quantile-quantile (Q-Q) plots were used to assess how the results of the linear models matched the expected null distribution and to identify biases in the results (Figures S2B and S3B). When adjusted for surrogate variables, a close approximation of linearity was achieved in the RNA sequencing data results, with p-values deviating from null only at the high end of the distribution, indicating true positive results (Figure S2B). The DNA methylation data showed a linear but deflated distribution of p-values which was partially corrected after surrogate variable adjustment (Figure S3B).

3.2. Differential Expression and Methylation Analyses

Seven genes were found to be differentially expressed in TSD compared to controls (FDR < 0.05, Figure 1A and Table S2). The fold changes of the seven genes ranged between a threefold to an eightfold change in expression (Table S2). When a less stringent FDR was applied (< 0.2), a total of 15 DEGs was detected (Table S2). Gene annotation for the 15 genes with FDR < 0.2 is shown in Table S3. None of the DNA methylation probes reached genome-wide significance (Figure 1A), however 29,868 probes had p-values < 0.05, and 1,981 had p-values < 0.005. The expression of one gene, RNY1, was found to be significantly associated with YGTSS score (FDR < 0.05, Table S4).

3.3. Expression Quantitative Trait Methylation (eQTM) Analysis

To identify methylation-regulated DEGs in TSD, an expression quantitative trait methylation (eQTM) analysis was performed based on LCPM and M-values. A total of 1,201 genes (unadjusted p < 0.05) and 1,981 probes (unadjusted p < 0.005) were included in the analysis. The eQTM analysis resulted in the identification of 236 potentially methylation-regulated genes with p-values < 0.05 (Table S6). The top 20 eQTMs have been characterized in Table S7. Gene expression and DNA methylation levels for the top 12 most significant eQTMs are visualized in Figure 1B. All eQTMs plotted had an unadjusted p-value < 2.5e-4. In all cases except one, the correlation between gene expression and DNA methylation was observed in all samples, and the average methylation changes ranged from 1.8% to 10.0% (Figure 1B).

3.4. Gene Set Enrichment Analysis (GSEA)

Gene set enrichment analysis (GSEA) was conducted to identify differentially regulated ontologies and pathways in TSD compared to non-TSD individuals. We identified 24 GO terms associated with TSD, of which 75% were pathways related to translation and mRNA processing while 25% were neurological GO terms (Table S8). As shown in the category net plot in Figure 2A, all the significantly altered pathways were downregulated, indicating that the effect on gene expression for TSD is mainly repressive. As for KEGG terms, four terms associate significantly with TSD (Figure 2B) – ribosome, nicotine addiction, nucleocytoplasmic transport, and signaling pathways regulating pluripotency of stem cells. If correction for multiple testing is omitted, 26 KEGG terms are identified (Table S9).

4. Discussion

In the present study, we included 11 twins (five twins discordant for TSD, four twins concordant for TSD, and two asymptomatic twins) and analyzed gene expression and DNA methylation data from five monozygotic twin pairs discordant for TSD. Studying discordant monozygotic twins allows us to investigate environmental or other non-genetic factors underlying disease susceptibility, as these individuals share nearly identical genetic backgrounds. We identified seven DEGs with FDR < 5% and 15 DEGs with FDR < 20% (Figure 1A, Table S2). Of the 15 DEGs, 11 were lncRNAs and pseudogenes, and most of them were at the time of writing poorly characterized (Table S3) [20,21,22]. To complement the TSD/non-TSD analysis, we investigated associations between gene expression and YGTSS score in 13 individuals affected by TSD. Of the top 20 DEGs associated with YGTSS score, 10 genes were involved in regulating immune responses [23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39], and four could be linked to neurological phenotypes [40,41,42,43,44]. Annotation for the YGTSS-associated DEGs can be found in Table S5 [45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67]. Immune abnormalities have been suggested in the pathogenesis of several movement and psychiatric disorders [68]. Association between immune-based mechanisms and tic exacerbation has been suggested [69]. Furthermore, an association between TSD and autoimmune disorders has previously been reported [3,70,71], and first-degree relatives to TSD patients have been shown to have an increased prevalence of autoimmune disorders such as systemic lupus erythematous (SLE) [70]. SLE patients are also 10-15 times more likely to have OCD, a common comorbidity of TSD [71]. The only gene that achieved genome-wide significance in the YGTSS analysis was RNY1, where an inverse relationship between RNY1 expression and YGTSS score was observed. RNY1 is a regulator of the protein Ro60, which is a common autoantigen in both SLE and Sjögren’s syndrome [23]. The earliest detectable autoantibodies in SLE are against Ro60, and autoantibodies precede symptoms, making it a clinically important biomarker [23]. Ro60 seems to have a protective effect against SLE, as Ro60 knockout mice models develop an SLE-like syndrome (Xue et al., 2003). As RNY1 assists Ro60 in its function, depletion of RNY1 might attenuate the protective effect of Ro60 and be an indication of autoimmunity. Further studies with larger cohorts are necessary to validate this hypothesis. Among the top twenty most significant DEGs associated with YGTSS we identified RNU4-2, BATF2, DUX4L26, and MIR4718, all which have been associated with neurodevelopmental or neurodegenerative phenotypes [73,74,75,76]. The DEG RP42 Homolog (RP42) Pseudogene (RP1-121G13.3) has been linked to ADHD, a common comorbidity for TSD patients [77] (Supplementary Table S7).

To identify DEGs regulated in cis by DNA methylation, we performed an eQTM analysis, from which we identified 12 cis eQTMs (Figure 1B). The most significant expression-methylation association was between chymotrypsin C (CTRC) and cg24014317. Chymotrypsin C is a protease with chymotrypsin-like specificity. The enzyme is also involved in calcium homeostasis and regulation of digestive enzymes such as trypsin and some carboxypeptidases [78,79]. Like its isozyme chymotrypsin, it is highly and almost exclusively expressed in the pancreas. Studies have shown that 60-65% of individuals with ASD have low levels of chymotrypsin, as well as low levels of circulating amino acids, suggesting that increased protease activity might be beneficial [80,81,82]. Pancreatic replacement therapy has recently been investigated as a treatment for ASD [83].

Finally, we identified dysregulated gene ontologies and pathways through gene set enrichment analysis (Figure 2). Most of the detected GO terms were related to translation and mRNA processing, with a select few being related to neurological processes, as one might expect from a neurological disorder (Supplementary Table S8). The identified KEGG terms were related to translation, general cell biology, stem cells, and addiction (Supplementary Table S9). The latter may reflect the ADHD comorbidity, as between 20-50% of individuals with ADHD have been reported to have a substance-use disorder [84,85]. Interestingly, all genes belonging to the significantly detected GO or KEGG terms were downregulated in TSD patients relative to controls, suggesting an overall gene repressive effect in the disorder.

A limitation of the present study was its small sample size, although the monozygotic twin setup was designed to increase the strength by limiting non-disease specific genetic variation and to allow for the identification of causative factors rooted in environmentally induced changes in gene expression and DNA methylation. Our study design thus ensured almost identical genetic background, such that significant changes in gene expression and/or DNA methylation can be attributed to non-genetic contributions to the pathology.