article: Inappropriate activation or inadequate regulation of CD4+ and CD8+ T cells may contribute to the initiation and progression of multiple autoimmune and inflammatory diseases. Studies on disease-associated genetic polymorphisms have highlighted the importance of biological context for many regulatory variants, which is particularly relevant in understanding the genetic regulation of the immune system and its cellular phenotypes. Here we show cell type-specific regulation of transcript levels of genes associated with several autoimmune diseases in CD4+ and CD8+ T cells including a trans-acting regulatory locus at chr12q13. 2 containing the rs1131017 SNP in the RPS26 gene. Most remarkably, we identify a common missense variant in IL27, associated with type 1 diabetes that results in decreased functional activity of the protein and reduced expression levels of downstream IRF1 and STAT1 in CD4+ T cells only. Altogether, our results indicate that eQTL mapping in purified T cells provides novel functional insights into polymorphisms and pathways associated with autoimmune diseases.
  T cells are essential elements of the adaptive immune response [1]. CD4+ T cells, together with an appropriate cytokine environment, are required for the activation and differentiation of CD8+ T cells that mediate defense and pathogen clearance during various infections [2]. The role of CD4+ T cells is also necessary for B cells and macrophages to execute their protective functions. T cells are activated and start to differentiate in response to complex stimuli involving cytokines, membrane receptors and transcription factors. Signaling through T cell receptor and costimulatory molecules induces their rapid proliferation and clonal expansion, and differentiation into effector and memory T cell subtypes [3]. Faulty activation or inadequate regulation of CD4+ and CD8+ T cells may contribute to the initiation and progression of multiple autoimmune diseases, including type 1 diabetes (T1D), rheumatoid arthritis (RA), autoimmune thyroiditis, systemic lupus erythematosus (SLE), multiple sclerosis, psoriasis, inflammatory bowel disease, as well as allergy and asthma [4,5]. In addition, a distinct lineage of CD4+ T cells, regulatory T cells, have a central role in the control of the pathogenesis of autoimmune and inflammatory diseases [6]. Genome-wide association studies (GWAS) have identified thousands of single nucleotide polymorphisms (SNPs) associated with various immune related diseases [7–12]. Studies of the consequences of the risk alleles in biologically relevant contexts for different diseases have stressed the importance of the tissue and cell type-specificity of the regulatory variants [13–15]. We have recently shown that by mapping eQTLs in a sufficiently large sample set, it is possible to identify cell type-specific effects in whole blood, but the challenge to distinguish the cells responsible for the associations remains [16,17]. The power of studying expression-associated genetic variants in purified cells types has now been illustrated for B cells, monocytes, neutrophils, T regulatory cells and CD4+ T cells [18–21], which allowed the identification of functional roles for several polymorphisms at autoimmune and even neurodegenerative disease loci. Meanwhile, some particularly interesting cis- and trans-eQTLs identified in whole blood could not be attributed to these cell types, and require an expanded survey of cells involved in immune response. To this end, we purified CD4+ and CD8+ T cells from the peripheral blood of 313 healthy individuals for genome-wide mapping of genetic variation affecting the expression of genes involved in immune response. Our analysis characterizes both the extent of genetic control of gene expression in T cells and its cellular specificity, as well as variation in isoform levels of genes and their association with epigenetic changes. We show that integrating knowledge from GWAS with eQTLs enables us to clarify the functional consequences of disease-associated variants and assess the enrichment of autoimmune response. Finally, we highlight a common T1D-associated missense variant in IL27 affecting the STAT1 and IRF1 pathway in CD4+ T cells. Our analysis provides insights into the basic processes of the regulation of gene expression in T cells and advances our understanding about the pathways involved with disease susceptibility in the adaptive immune system.
  We purified CD4+ and CD8+ T cells from peripheral blood mononuclear cells (PBMCs) of 313 healthy European individuals from the Estonian Biobank [22]. The purified cells were subjected to genome-wide gene expression analysis, genotyping and imputation using the 1000 Genomes reference panel. After stringent quality control and filtering, close to 6 million SNPs, and expression data from 38,839 probes within 23,704 genes were included in the analysis. DNA methylation of 450,000 CpG sites was determined in the same purified cells of 100 of the individuals [23] to add explanations to observed differential expression between cell types. To characterize the extent of genetic control of gene expression in T cells and its cellular specificity, we first tested the association between SNPs and gene expression within 1Mb intervals, referred to as cis-eQTLs. In total, we identified cis regulatory SNPs for 2,605 genes in CD4+ T cells and 2,056 genes in CD8+ T cells at probe-level false discovery rate (FDR) < 0. 05 (Fig 1A, S1 Table), with an overlap of 1,637 genes and estimated replication rate π1 = 0. 99. The high replication rate reflects the similarity of the two cell types compared, and also highlights the limitations of arbitrary cut-off levels for significance. Moreover, about half of the significant eQTLs detected in a large meta-analysis of peripheral blood from 5000 individuals [16] could be replicated in the CD4+ and CD8+ T cells (π1 = 0. 51 and 0. 45, respectively, S1 Fig), indicating a high level of T cell-specific eQTL effects [24,25]. As examples of cis-eQTLs in CD4+ and CD8+ T cells, we observed evidence of genetic regulation of long coding RNA RP11−534L20. 5 located 7 kb downstream of the IKBKE gene which has an essential role in regulating inflammatory responses to viral infections; STAT6, a transcriptional activator involved in T cell differentiation, and STX2 (syntaxin 2), involved in intracellular transport of vesicles (Fig 1B). We next used a Bayesian model averaging and hierarchical modelling [26] that combines information across genes to jointly call eQTLs in one framework to assess the specificity and proportion of shared eQTLs in CD4+ and CD8+ T cells. Of 3871 genes associated with eQTLs, all showed very strong posterior probability to be shared between CD4+ and CD8+ T cells (S2 Table). One of the strengths of the multi-tissue joint cis-eQTL analysis is its capability to account for incomplete power [26,27]. As our sample size for both of the cell types was similar, we noted the advantage of Bayesian methods in particular for genes with eQTLs that have modest effects. Moreover, for eQTLs where the absolute effect size was higher in one cell type, there was also a tendency for higher expression levels in that gene in the given cell type (chi-square test P-value < 2x10-16, S2 Fig). For a standard tissue-by-tissue analysis, higher effect sizes and small standard deviations lead to more statistical evidence against the null, thus resulting in lower P-values. The effect of the magnitude of the expression levels is exemplified by the effect sizes and P-values of the probes in the inhibitory immune checkpoint gene, CTLA4, known to be associated with T1D and other endocrine autoimmune diseases [28,29]. CTLA4 is covered by three probes all showing consistent eQTL effects shared by both cell types in the multi-tissue analysis (Fig 2A). In the tissue-by-tissue analysis, expression levels of CTLA4 were associated with 122 SNPs (S3A Fig), including T1D and RA-associated variants rs231775 (+49A/G) [28], rs3087243 (CT60A/G) [28,30], rs231806 (MH30C/G) [28], and rs231735 [8] (pair-wise r2 between the SNPs ranging from 0. 51 to 0. 90, and strong linkage disequilibrium (LD, r2 > 0. 8) between MH30 and CT60, S3C Fig). The CT60 polymorphism is also in strong LD with length variation of the (AT) n repeat in the 3’UTR of CTLA4 with long variants of the repeat associating with lower CTLA4 mRNA expression in autoimmune T cell clones [31]. Our results extend this finding to the polyclonal population of CD4+ and CD8+ T cells, where the disease risk alleles were associated with lower expression levels of CTLA4 (S3B Fig). Genetic control of the expression of the alternatively spliced third exon, that is omitted from the secretory variant of CTLA4 [28], was significantly detectable in CD4+ T cells only, whereas genetic association with the expression of the fourth exon was significant in both cell types (Fig 2B and 2E). Further investigation of the CTLA4 gene locus revealed differential methylation between CD4+ and CD8+ T cells (Fig 2C). Specifically, methylation of a CpG site near the second exon of CTLA4 (cg26091609) was inversely correlated with the expression levels of the three different probes in the gene (Fig 2D), possibly explaining the expression differences in CD4+ and CD8+ T cells, and thus failing to detect significant eQTL effects in the tissue-by-tissue analysis. Interestingly, divergent DNA methylation between CD4+ and CD8+ T cells was similarly detectable for an additional 73 out of 187 genes covered by more than one probe, and with seemingly isoform-specific eQTL effects (S3 Table). Among the 74 genes, multi-tissue analysis resulted in the identification of eQTL effects for two or more probes in 66 genes which were shared between the cell types. These observations suggest that the differential splicing effects, plausibly regulated by epigenetic DNA methylation, influence the gene expression levels, and thus the power to significantly detect eQTLs in the tissue-by-tissue analysis. For the analysis of SNPs affecting the expression of distal genes (>5 Mb apart), referred to as trans-eQTLs, we selected all 4,638 genome-wide significant (P-value < 5x10-8) SNPs from the GWAS catalog [32] (accessed 24/03/2015). After correcting gene expression levels for cis-eQTL effects, we identified 36 and 40 GWAS SNPs associated with the expression levels of 209 and 378 distant genes in CD4+ and CD8+ T cells, respectively (overlap of 21 SNPs and 133 genes; Fig 3A, S4 Table). The functions of the genes associated with the trans-acting GWAS SNPs highlighted their role in T1D (Ingenuity pathway analysis, P = 4. 39x10-5, Fig 3A) and mTOR signaling (P = 3. 84x10-3, Fig 3A) in CD4+ and CD8+ T cells, respectively. Strikingly, we observed a trans-acting regulatory locus at chr12q13. 2 in CD4+ and CD8+ T cells with broad-range impact on hundreds of genes. Five SNPs in that region were previously implicated in B cell-specific trans-associations by Fairfax et al. [18]. In T cells, we identified the lead SNP in this region as rs1131017 located in the 5´UTR of the ribosomal small subunit protein RPS26 gene (S5 Table). The rs1131017 SNP is in LD with eight GWAS SNPs implicated in T1D [34,36,40], vitiligo [33,35] and other autoimmune and inflammatory diseases [30,37–39]. The pair-wise r2 values of the nine SNPs are shown in S4B Fig and their eQTL effects in S4C Fig, illustrating the strong LD (r2 > 0. 8) between the top trans-eQTL SNP rs1131017 and several of the GWAS SNPs. In addition to RPS26, the trans-acting region contains the CDK2, RAB5B, SUOX, IKZF4 and ERBB3 genes (S4A Fig) and associates with the expression levels of 187 and 351 genes in CD4+ and CD8+ T cells, respectively, with an overlap of 124 genes (Fig 3B). Many of these genes are highly expressed and have specific roles in T cells such as CTLA4, GZMA, GZMB, GZMH, GNLY and CD8A (Fig 3B). We also identified a trans-eQTL (GWAS SNP rs4788084[T]) on chr16p11. 2 close to the IL27 gene associated with lower expression of IRF1 (P = 1. 84x10-9) and STAT1 (P = 2. 91x10-8) in CD4+ T cells only. The rs4788084[T] is associated with lower expression of STAT1 in peripheral blood [16] and reduced risk for T1D [7,11]. To further explore the effects of genetic variants in the IL27 region, we mapped trans-eQTLs for all SNPs in the chromosome 16: 28,2–29,1 Mb region (829 SNPs in total, S6 Table, Fig 4A). This analysis revealed an even stronger trans-eQTL signal for a missense SNP rs181206[G] within the IL27 gene (Fig 4B) and decreased expression levels of IRF1, STAT1 and REC8 (Fig 4C). The rs181206[G] allele, with a minor allele frequency of 29% in the 1000 Genomes phase 3 European population, causes the amino acid Leu119Pro change in the alpha-helical domain of IL-27 (S5 Fig). Strong LD (r2 ranges from 0. 62 in the Estonian population to 0. 84 in the other European populations, 1000 Genomes phase 3) between rs181206 and the GWAS SNP rs4788084 supports the protective role of the rs181206[G] allele for T1D among individuals of European ancestry. Of note, the effect of the missense SNP on IRF1 and STAT1 remains after regressing out the effect of the GWAS SNP, but not the other way around. Moreover, Bayesian test for colocalisation [41] adds strength to the hypothesis that T1D susceptibility and changes of expression in IRF1 and STAT1 are associated with the region, and share a single causal variant with the most posterior support for rs181206[G] (S7 Table). Importantly, despite the relatively small sample size, our focus on T cell subtypes enabled us to identify this effect as specific to CD4+ T cells, as we did not detect significant trans-eQTL effects in this region in CD8+ T cells (S6 Fig). We also confirmed the absence of the signal in B cells and monocytes by re-analyzing the data from Fairfax et al. [18] (S6 Fig). Notably, the trans-eQTL locus had no effects on the expression level of IL27 in cis (S7 Fig). In agreement with the cell type-specific effect, we found higher expression levels of the IL27RA and IL6ST (gp130) genes, which together act as a receptor for the IL-27 cytokine, in CD4+ cells in comparison to CD8+ T cells (Fig 4D). Given the strong positive correlation between the expression levels of IRF1 and STAT1, we used structural equation modeling to determine whether the IL27 SNP rs181206 affects both genes independently or via each other. Overall, the best-fitting scenario suggested IRF1 to mediate the SNP and STAT1 relationship (model 1, S8 Table). The finding was supported by a simulation experiment (S8 Fig) and suggested a mechanism for the effect of IL-27 on IRF1 and STAT1 expression (S9 Fig). In order to confirm the functional effect of the IL27 rs181206[G] allele we investigated its effect on IRF1 and STAT1 expression in human PBMCs by additional experiments. IL-27 is produced by innate immune cells, and after forming a heterodimer with EBI3, it interacts with its receptor IL27RA and activates the STAT1/STAT3 pathway in T cells [42]. After binding to interferon stimulated response elements (ISRE), the STAT1/STAT3 pathway induces transcription of several interferon-induced genes, including IRF1 and STAT1 itself. We cloned cDNA variants of the IL-27 wild-type (Leu119) and missense (Pro119), as well as EBI3. After transfection into HEK293 cells, we combined the cell supernatants containing either IL-27 Leu119 or Pro119 protein with an equal amount of EBI3 protein and studied their effect on IRF1 and STAT1 expression by real-time PCR in human PBMCs from four healthy individuals. As shown in Fig 4E, the missense SNP resulting in Pro119 in IL-27 induced lower STAT1 and IRF1 transcript levels compared to IL-27 Leu119. The comparison of the fixed effects of the rs181206[G] allele resulted in highly significant estimates: β^mutant = -1. 11 (P = 5. 52x10-13) for IRF1 and β^mutant = -1. 28 (P = 4. 86x10-10) for STAT1. This result supports our trans-eQTL analysis, suggesting that the Pro119 in the IL27 gene is the causal SNP of these associations.
  We here report eQTL mapping in purified CD4+ and CD8+ T cells and reveal multiple effects on regulation of genes associated with autoimmune diseases. The eQTL studies in cell subtypes have both demonstrated a high level of specificity as their effects vary across cell types, as well as remarkable sensitivity despite their several fold smaller sample sizes [24,25]. Indeed, we also observed many cis-eQTLs identified in CD4+ and CD8+ T cells that were shared according to the joint analysis, but were seemingly cell type-specific by the tissue-by-tissue analysis. The genes with such notable patterns include CTLA4, known for its inhibitory role in T cell mediated immune responses and association with many organ-specific autoimmune diseases [29]. The autoimmunity-susceptible SNPs in this region (rs231775, rs3087243, and rs231806) have been associated with various functional effects including CTLA4 transcript levels, splicing, production of the soluble form of CTLA-4 and posttranslational modifications [28,43–45]. In addition, our results show differential DNA methylation in the CTLA4 region to be associated with the potential isoform-specific differences in CD4+ and CD8+ T cells. This is likely due to negative correlation between DNA methylation and gene expression, which eventually leads to lower expression levels and moderate effect sizes. To confirm whether the observed effects are due to true differences in isoforms instead of the quality of the designed probes at capturing variation in gene expression levels, a follow-up study with RNA-seq data is needed. The observed trans-acting regulatory region on chr12q13. 2 in CD4+ and CD8+ T cells indicates the important role of genetic variants in that region affecting the expression levels of over a hundred genes across the genome. Interestingly, the region includes the RPS26 gene, which may constitute a mechanism for the detected trans-eQTL effects. Our lead SNP rs1131017 is located in an oligopyrimidine tract of the 5´UTR of the RPS26 mRNA and its T1D risk allele rs1131017[C] correlates positively with RPS26 expression levels [46,47]. The oligopyrimidine tract controls the translation of many mammalian ribosomal protein genes [48], and the effect of the SNP on RPS26 ribosomal distribution has been reported [49]. RPS26 is a main component of the ribosomal region involved in the recruitment of cellular mRNA during translational initiation and in maintenance of the path of mRNA molecules to the ribosomal exit site [50]. Hence, it is conceivable that changed RPS26 protein levels may affect the stability or translational efficiency of a large number of cytosolic mRNAs. Nevertheless, the exact functional role of rs1131017[C] in T1D remains to be identified, as well as the effect of other candidate genes in this region such as IKZF4, ERBB3 [51] and CDK2 [52]. Strikingly, we identified a common missense variant in cytokine IL27 as a significant trans-eQTL for IRF1 and STAT1 in CD4+ T cells. The effect of IL-27 in T cells is regarded as anti-inflammatory but it has also been shown as a growth and survival factor for T cells [42]. As a heterodimer with EBI3, IL-27 activates the STAT1/STAT3 pathway that induces the transcription of several interferon-induced genes, including IRF1 and STAT1 itself. Our model suggests that IRF1 mediates the effect of the IL27 SNP on expression levels of STAT1. Moreover, our functional studies with the mutated form of IL-27 confirmed its decreased capacity to activate the STAT1 pathway, and we showed that a potential causal variant (the missense variant rs181206) for T1D susceptibility and changes in IRF1 and STAT1 expression is shared. Furthermore, our findings are supported by studies of a T1D mouse model with high levels of IL-27 and delayed T1D onset after treatment with an IL-27 blocking antibody [53]. Altogether these results suggest that the rs181206[G] variant of the IL27 gene confers protection against T1D through the inhibited expression of IRF1 and STAT1 in CD4+ T cells. Our results also suggest that IL-27 may promote autoimmunity toward pancreatic islets via the upregulation of the STAT1 pathway. In conclusion, the analysis of genetic modulators of gene expression profiles as intermediate phenotypes between human traits and underlying genetic variation offers new instruments to refine our understanding of disease susceptibility. Moreover, eQTL studies in purified cell types instead of whole blood enable us to establish the specific mechanisms and pathways involved in diseases progression and create the basis for future explorations and drug interventions.
  The aims of this study were to perform cis- and trans-eQTL mapping in purified CD4+ and CD8+ T cells and evaluate their pathogenic implications for autoimmune mechanisms. Based on similar studies [18], we considered the sample size of 300 individuals to be sufficient to ensure statistical power to detect eQTL effects in purified cells. This study participants were healthy donors of the Estonian Genome Center of the University of Tartu [22]. In total, 313 subjects were selected for the study, with median age 54 (standard deviation 17. 8), 154 females and 159 males. The study was approved by the Ethics Review Committee of Human Research of the University of Tartu, Estonia (permission no 206/T-4, date of issue 25. 08. 2011) and it was carried out in compliance with the Helsinki Declaration. A written informed consent to participate in the study was obtained from each individual prior to recruitment. All methods were carried out in accordance with approved guidelines. DNA from the samples were genotyped using HumanOmniExpress BeadChips (Illumina), according to the manufacturer’s instructions. We imputed both datasets using the 1000 Genomes project reference by using IMPUTE v2 [54], resulting in 5,879,386 autosomal SNPs with a minor allele frequency (MAF) of > 0. 05 for downstream analyses. CD4+ and CD8+ T cells were extracted from the peripheral blood mononuclear cells (PBMC) by consecutive positive separation using microbeads (CD4+ #130-045-101; CD8+ #130-045-201) and AutoMACS technology (Miltenyi Biotec) according to the manufacturer' s protocol. Gene expression data was generated using HumanHT-12v4 BeadChips (Illumina), according to the standard protocol. Preprocessing and quality control of the data was done using R [55] and the Bioconductor packages lumi [56] and arrayQualityMetrics [57]. The number of samples retained for further analysis was 293 for CD4+ and 283 for CD8+ T cells, 303 unique individuals. The effect of the mutant and wild-type IL27 was experimentally tested in HEK 293 cells. The procedures are described in details in the S1 File. The effects of SNPs on local (cis-eQTL) and distant (trans-eQTL) genes were determined via eQTL mapping as described in eQTL mapping analysis cookbook developed by the University Medical Center Groningen at the Genetics Department and the Genomics Coordination Center (https: //github. com/molgenis/systemsgenetics/wiki/eQTL-mapping-analysis-cookbook) and previously in [16]. For cis-analysis the distance between the probe midpoint and SNP genomic location was up to 1 Mb and for trans-analysis the distance was more than 5 Mb or the probe and SNP were on different chromosomes. Only SNPs with a MAF > 0. 05, call rate > 0. 95 and a Hardy-Weinberg equilibrium P-value > 0. 001 were included in the analyses. We used the Spearman correlation coefficient (Spearman’s rho) to detect associations between the coding allele of the SNP (directly genotyped or imputed allele dosages) and the variations in gene expression levels (residual gene expression levels obtained after corrections as described in “Gene expression quality control and normalization” in Supplementary Methods). To control for multiple testing, we applied more conservative probe-level false discovery rate (FDR) procedure. The eQTL mapping procedures are described in details in the S1 File. To estimate the sharing of cis-eQTLs between the two cell types and peripheral blood we used the π1 statistic (Storey and Tibshirani q-value approach [58] implemented in the R qvalue package, default settings used). Based on the list of P-values, the overall proportion of true null hypotheses (π0) is estimated, i. e. those following the Uniform (0,1) distribution. An estimate of true alternative tests is π1 = 1 – π0, i. e. π1 is the proportion of true positives. The reported replication rate (proportion of sharing) between two groups (CD4+ or CD8+ T cells or pheripheral blood) is the π1 statistic estimated by taking the significant SNP-probe list from one group and using the corresponding P-value distribution in the other group. We used Bayesian framework for multi-tissue joint eQTL analysis by Flutre et al. [26] in CD4+ and CD8+ T cells implemented in the eQtlBma program. The cis candidate region was defined as +/- 1 Mb from probe midpoint. As an input to the program, residual gene expression levels obtained after corrections as described above and the corresponding genotypes were used. Firstly, the raw Bayes factors for each configuration (eQTL active in CD4+, eQTL active in CD8+, eQTL active in both cell types) were computed assessing the support in the data for each probe-SNP pair being an eQTL. Secondly, information across all genes were combined by the hierarchical model with an EM algorithm to get maximum-likelihood estimates of the configuration probabilities. Thirdly, by Bayesian model averaging using the raw Bayes factors weighted by the estimated configuration probabilities, the posterior probabilities for the eQTL to be active in a given configuration were obtained. As a final result, only the best SNP per probe is picked, based on the posterior probability for a SNP to be “the” eQTL for a probe. To obtain the posterior probabilities, an estimate of the probability for a probe to have no eQTL in any cell type (π0) is needed. This was estimated using Storey and Tibshirani q-value approach [58] via probe-level P-values obtained by 10,000 permutations. The Ingenuity Pathway Analysis (IPA) tool was used to establish associations with processes and diseases. The standard IPA enrichment analysis was performed with the set of all human genes on the Illumina array as background. The R package coloc [41] was used to assess whether two association signals are consistent with a shared causal variant. Assuming a single causal SNP for both of the traits in a region, posterior support for following hypothesis is estimated: For trait1, we chose T1D susceptibility and obtained necessary summary statistics from Onengut-Gumuscu et al. [59] from www. t1dbase. org (accessed 13/09/2016) webpage. For trait2, we chose the association with IRF1 and STAT1 expression in CD4+ T cells. Following Giambartolomei et al. [41], we set prior probability (p1) that a variant is associated with trait1 to 10−4, prior probability (p2) for trans-eQTL effect to 10−5, and prior probability (p12) that a variant is associated with both traits to 10−6. To show the consistency of the results, we varied p2 and p12 values from 10−4 to 10−5 and 10−5 to 10−6, respectively. We compared the plausibility of potential models for the effects of the IL27 SNP on gene expression levels of STAT1 and IRF1 using structural equation modeling (SEM). First, to test the model’s concordance to the true population we used the chi-square test. Next, to assess the goodness of fit of the resulting models, we used the following measures: Comparative Fit Index (CFI), Tucker-Lewis Index (TLI) and Root Mean Square Error of Approximation (RMSEA). The CFI and TLI goodness of fit measures indicate good models for values higher than 0. 9. RMSEA values less or equal to 0. 05 indicate reasonable fit between the model and the data. To compare the models, we used the Akaike information criterion (AIC). Lower values of that theoretic measure indicate better models. We applied the Sobel test to test whether a mediator carries the influence of the causal variable to the outcome, i. e. testing whether the indirect effect of the causal variable on the outcome via the mediator is significantly different from zero. The test statistic was calculated using the Sobel formula [60]: ab/b2sa2+a2sb2, where a is the effect of the causal variable on the mediator and b is the effect of the mediator on the outcome and sa and sb are the corresponding standard errors of a and b. We conducted a simulation study to investigate the relationship between the coding allele of the SNP and expression levels of two genes, IRF1 and STAT1. Therefore, we performed 1000 simulations using a sample size of n = 1000. We generated SNP genotypes with the minor allele frequency being the sum of two random binary traits from a binomial distribution B (1,0. 37), IRF1 expression levels depending only on the SNP (irf1 = 0. 903 − 1. 237 × snp + N (0,1) ) and STAT1 expression levels depending only on IRF1 expression levels (stat1 = 0. 031 − 0. 534 × irf1 + N (0,1) ). Then using simulated data, we estimated the significance of the SNP allelic effect and IRF1 expression levels on STAT1 expression levels and the SNP allelic effect and STAT1 expression levels on IRF1 expression levels by linear models. Graphs were generated using R [55] base packages and ggplot2 [61]. The LD plots were generated using LDheatmap [62], the regional association plots were generated with Locus Zoom [63] using hg19/1000 Genomes Nov 2014 EUR for showing LD between the variants, and the circos plots were drawn using the Circos visualization tool [64]. We used linear mixed effects models to estimate the mean of the differences in methylation and expression levels between CD4+ and CD8+ T cells assessing random intercepts for each of the individuals adjusted for sex, age and batch effects (chip and position on chip). The effect of the mutant SNP on IRF1 and STAT1 transcript levels was also evaluated by linear mixed effects models.
summary: Variation in regulatory regions as well as coding regions of the genome can affect the expression of genes. Many of these variants have been associated with different diseases and other traits, but the underlying biological pathways are often left unknown. Analysing the effect of single nucleotide polymorphisms (SNPs) on gene expression levels, referred to as expression quantitative trait loci (eQTL), in specific cell types can be used to gain insight into specific mechanisms of disease. By analyzing eQTLs in CD4+ and CD8+ T cells, essential elements of adaptive immune response, we identified both cis- and trans-acting SNPs in genes associated with several autoimmune diseases.
keywords: blood cells
  genome-wide association studies
  medicine and health sciences
  immune cells
  gene regulation
  immunology
  genome analysis
  cytotoxic t cells
  epigenetics
  dna
  molecular biology techniques
  dna methylation
  chromatin
  research and analysis methods
  white blood cells
  genomics
  gene mapping
  animal cells
  chromosome biology
  gene expression
  t cells
  chromatin modification
  dna modification
  molecular biology
  genetic loci
  biochemistry
  cell biology
  nucleic acids
  genetics
  biology and life sciences
  cellular types
  computational biology
  human genetics
title: Pathogenic implications for autoimmune mechanisms derived by comparative eQTL analysis of CD4+ versus CD8+ T cells