Natural genetic variation profoundly regulates gene expression in immune cells and dictates susceptibility to CNS autoimmunity.

Regulation of gene expression in immune cells is known to be under genetic control, and likely contributes to susceptibility to autoimmune diseases such as multiple sclerosis (MS). How this occurs in concert across multiple immune cell types is poorly understood. Using a mouse model that harnesses the genetic diversity of wild-derived mice, more accurately reflecting genetically diverse human populations, we provide an extensive characterization of the genetic regulation of gene expression in five different naive immune cell types relevant to MS. The immune cell transcriptome is shown to be under profound genetic control, exhibiting diverse patterns: global, cell-specific and sex-specific. Bioinformatic analysis of the genetically controlled transcript networks reveals reduced cell type specificity and inflammatory activity in wild-derived PWD/PhJ mice, compared with the conventional laboratory strain C57BL/6J. Additionally, candidate MS-GWAS (genome-wide association study candidate genes for MS susceptibility) genes were significantly enriched among transcripts overrepresented in C57BL/6J cells compared with PWD. These expression level differences correlate with robust differences in susceptibility to experimental autoimmune encephalomyelitis, the principal model of MS, and skewing of the encephalitogenic T-cell responses. Taken together, our results provide functional insights into the genetic regulation of the immune transcriptome, and shed light on how this in turn contributes to susceptibility to autoimmune disease.

Introduction

Over the last century, an increase in incidence and prevalence in many autoimmune diseases, such as multiple sclerosis (MS) [1], rheumatoid arthritis (RA) [2], type 1 diabetes [3], and systemic lupus erythematosus (SLE) [4], has been documented, and now these diseases impose a very significant public health burden [5]. The etiology of autoimmune disease is highly complex and multifactorial, owing both to increased genetic heterogeneity in human populations and diverse environmental influences. The contribution of the genetic component has been increasingly better defined, with early studies identifying the profound influence of the major histocompatibility complex (MHC) haplotypes, and linkage studies and the more recent genome-wide association studies (GWAS) identifying hundreds of additional disease-modifying loci [6, 7]. In concert with progress in human genetics, appropriate animal models are critical to a mechanistic understanding of complex autoimmune phenotypes. The reverse genetics revolution in the mouse has provided numerous critical insights into gene function, mostly through the use of knockout approaches. However, such “wreck-and-check” approaches yield only limited information applicable to the understanding of the impact of natural genetic variation at the population level. In this regard, classical quantitative trait locus (QTL) mapping studies in inbred mice are more useful, but these have been hampered by the limited genetic diversity of commonly employed laboratory mouse strains [8]. To overcome this limitation, so-called wild-derived inbred mouse strains have been established, such as PWD/PhJ (PWD), belonging to the Mus musculus musculus subspecies. These mice are genetically highly divergent from the standard laboratory strains, thereby more accurately modeling the greater evolutionarily-selected genetic diversity seen in human populations [9, 10]. Additionally, consomic strains of C57BL/6J (B6) mice carrying chromosomes from PWD (B6.ChrPWD) have been established [11], and have been useful in mapping QTL controlling various complex phenotypes [12, 13]. We have recently utilized this approach to begin to map QTL controlling susceptibility to experimental autoimmune encephalomyelitis (EAE), the principal autoimmune animal model of MS [14].

The consomic model carries predominantly the B6 genome, however, and thus it is also limited by the loss of many genome-wide epistatic interactions and trans-eQTL. Here, using the parental B6 and PWD strains of mice, we assessed the impact of natural genetic variation distinguishing these two strains on basal gene expression in five major immune cell types, as well as the outcomes in an autoimmune disease model and the associated immune responses. We found striking differences in basal immune cell gene expression that were genetically regulated and cell type-specific, and a smaller subset of genes whose expression was regulated in a sex-specific manner. Bioinformatic analyses identified several critical differentially regulated cellular pathways and processes, and predicted a dampened basal immune response in PWD compared to B6. Accordingly, we found that PWD mice were highly resistant to EAE induction, and exhibited altered encephalitogenic immune responses.

Results

Striking differences in immune cell gene expression exist between B6 and PWD mice

In order to understand how natural genetic variation and sex impact gene expression in the immune system, we performed extensive transcriptomic profiling of immune cells from male and female B6 and PWD mice. The following five major immune cell types were isolated from spleen and lymph nodes of naïve mice. B cells (BCs) were isolated by positive selection for the surface marker CD19. The remaining cell types were isolated by fluorescence-activated cell sorting (FACS) as follows: CD8 T cells (CD45+NK1.1−CD19−TCRβ+CD8+), CD4 T cells (CD45+NK1.1−CD19−TCRβ+CD4+CD25−), CD25+ T regulatory cells (Treg; CD45+NK1.1−CD19−TCRβ+CD4+CD25+); and myeloid antigen presenting cells (APCs; CD45+NK1.1−CD19−TCRβ−CD11b+CD11c+). RNA was isolated, and genome-wide gene expression was assessed using the Illumina Bead Array platform. Principal Component Analysis (PCA) was used to reduce the global expression differences between samples into a limited number of vectors, capturing a total of 64% of the variation in gene expression along three principal components (PCs) (Supplementary Fig. 1). Cell type-specific differences were captured primarily by PC #2 and to a lesser extent PC #3, which gave a clear separation between APC, BCs, and the three different T cell subsets, which clustered closer to each other, as expected. Strain-specific differences were captured primarily by PC #3, showing distinct separation between B6 and PWD samples. Sex-specific differences were much more subtle, and were captured partly by PC #1, although in most cases male and female samples clustered close together by strain and cell type. CD4 T cell and Treg samples showed the most sample-to-sample variability, which was also partly captured by PC #1.

In the first set of analyses, expression data from males and females were pooled to assess the effect of strain. We found striking differences in gene expression between B6 and PWD. Using a conservative filter of false discovery rate (FDR) < 0.05 and a fold change (FC) > |2|, hundreds to thousands of genes were differentially expressed (DE), ranging from −267 to 170 fold difference in expression, with 865 DE genes in CD4 T cells, 1,044 in CD8 T cells, 1,357 in APCs, 1,097 in B cells, and 557 in Treg, representing a total of 2,764 unique DE genes in at least one cell type (Fig. 1 and Supplementary File 1). A pooled analysis of the five cell types revealed a core subset of 822 DE genes between B6 and PWD across all five cell types. Consistent with this, a survey of the top 20 upregulated (defined as having higher expression in PWD compared with B6) or top 20 downregulated (lower expression in PWD compared with B6) DE genes illustrates that some expression differences occurred uniformly across most cell types, e.g., Ifi27 or Actb, while others were cell type-specific (Fig. 1B). Examples of cell type-specific DE genes include CD4-specific Pdlim4; CD8-specific Klrd1; CD4- and CD8-specific Cd163 and Cldn10; Treg-specific Lcn10; APC-specific Fcer1g, Cd59a, and Chi3l3; and B cell-specific Blk and several Igk genes. Additionally, we noted the lower expression of several genes encoded in the mitochondrial genome in PWD compared to B6. Lastly, many hypothetical and/or uncharacterized genes were abundantly represented among the DE genes (e.g., Fig. 1B), suggesting that expression and likely function of these genes is genetically determined. Overall, these results demonstrate that the natural genetic variation distinguishing B6 and PWD mice results in a remarkable level of genetic control over the immune cell transcriptomes.

Figure 1

Genetic control of gene expression in immune cells

(A) Top significantly DE probes between B6 and PWD (for each strain biological replicates, male, n=3; female n=3) immune cells exhibiting |FC| > 4, at an FDR < 0.05 in any of the five cell types are shown. FC indicates expression in PWD relative to B6. Mean indicates the FC average across all five cell types, which was used to sort the order of the probes. The |FC| > 4 cutoff was chosen to highlight the significantly DE genes with the largest effect size, and to facilitate visualization of the data. (B) Top 40 DE genes in PWD vs. B6 cells, identified by top 20 maximum FC (upregulated) and top 20 minimum negative FC (down-regulated) in any cell type. (C) and (D) Distribution and overlap of DE genes between selected cell types is shown. The number of genes indicates DE genes passing the filter of |FC| > 2 and FDR < 0.05, darker shading indicates a higher number of genes.

The PWD immune cell transcriptome predicts a lower basal activation state, lower autoimmune susceptibility, and promiscuous cell type-specific gene expression

In order to predict the overall consequences of the altered transcriptome in PWD immune cells, we undertook a bioinformatic analysis of the DE transcripts across each of the five cell types. Pathway analysis revealed that DE transcripts were enriched in immune related pathways, as expected, such as dendritic cell maturation, lymphotoxin and JAK/STAT signaling, etc. Strikingly, using expression directionality (i.e. whether the transcripts were up- or downregulated in PWD relative to B6) to predict the directionality of change in impacted pathways, we found almost uniform predicted dampening of the activity of enriched canonical pathways in PWD cells, with the notable exceptions of p53, protein kinase A, and death receptor signaling, which were predicted to have increased activity (Fig. 2A). Similar results were obtained using upstream regulator analysis, which identified upstream regulation by central immune mediators such as interferon (IFN)γ, interleukin (IL)-2, nuclear factor κ B (NFκB), and IL-1β, most of which exhibited lower activity in PWD cells (Fig. 2B). The minority of upstream regulators that were predicted to have enhanced activity in PWD cells, e.g. KLF3, FOXP3, and IL-10RA, tended to be immune-regulatory or neutral. Altogether, these results predict a lower basal activation state in PWD immune cells.

Figure 2

Predicted activation state of PWD immune cells

(A) Canonical pathway analysis of differential gene expression between PWD and B6 immune cells. The top significantly enriched canonical pathways (p < 0.01) are shown. The heat map indicates the Z-score, indicative of predicted direction of change (orange – upregulated; blue, downregulated). (B) Upstream analysis of differential gene expression between PWD and B6 immune cells. Top predicted activators (genes and proteins only) (p < 0.01) are shown. Enrichment analysis of MS-GWAS (C) or MS-CIS (D) genes within the DE transcripts between B6 and PWD was performed as described in the Materials and Methods. Enrichment p values [displayed as negative log(p)] for each cell type are shown. (E) Gene set enrichment analysis of DE genes in CD4 T cells and APCs was performed as described in Materials and Methods. Top panels show genes upregulated in PWD, the bottom show downregulated genes. (F) Quantitation of cell lineage-specific gene expression (as a measure of enrichment) across different cell lineages was performed as described in the Materials and Methods. Significance of differences between genes upregulated in PWD relative to B6 vs. those downregulated is indicated as follows: *, p<0.05; $, p<0.01, #, p<0.001. Abbreviations, Stem – stem cells, GN – granulocytes, TC – T cells, Str – Stroma, ILC = innate-like lymphoid cells.

To test whether the observed genetic regulation of immune cell transcriptomes had any implications for human autoimmunity, we tested whether GWAS candidate genes for MS susceptibility (MS-GWAS) [15] were enriched within the DE gene sets for each cell type. Significant enrichment of MS-GWAS genes was found only for transcripts that were downregulated in PWD cells relative to B6, but not for those that were upregulated (Fig. 2C), suggesting that PWD cells express lower levels of autoimmune susceptibility genes. Interestingly, for the downregulated genes, the level of MS-GWAS enrichment varied across cell type, with APCs, CD8 T cells, and Tregs showing the highest level of enrichment. To verify the specificity of this observation, we assessed the enrichment of a GWAS candidate gene set from a related immune-mediated disease, inflammatory bowel disease (IBD) [16], which exhibits a significant amount of genetic overlap with MS [17], and a second set of candidate genes for a non-immune-mediated neurological disease, autism spectrum disorder (ASD) [18, 19]. We observed some cell type-specific significant enrichment of the IBD gene set, but unlike the case for the MS-GWAS gene set, the enrichment was less pronounced and exhibited no directionality (similar enrichment of genes up- or downregulated in PWD) (Supplementary Fig. 2A and B). With regard to the ASD candidate gene set, no significant enrichment was observed, supporting the specificity of our findings with MS and IBD candidate genes (Supplementary Fig. 2C). Additionally, we compared level of enrichment of the DE genes in PWD cells within the set of transcripts that were reported to be upregulated in CD4 T cells isolated from early onset MS patients (clinically isolated syndrome; MS-CIS) relative to healthy controls [20]. As was the case for the MS-GWAS gene set, genes upregulated in PWD compared to B6 showed no significant MS-CIS enrichment, while the downregulated genes exhibited robust enrichment (Fig. 2D). Here, the most significantly enriched subset was CD4 T cells, which was expected. As a specificity control, we also compared the enrichment of genes differentially expressed in ASD brains relative to healthy control brains [18, 21]. With the exception of the transcript set upregulated in PWD APCs, no significant enrichment was observed (Supplementary Fig. 2D and E). Taken together, these results predict a CNS autoimmunity-resistant phenotype for PWD immune cells, which is driven by differential expression of MS-GWAS and MS-CIS signature genes across different cell types.

To test how the altered gene expression pattern in PWD cells affects cell type-specific genes, we performed a gene set enrichment analysis using the ImmGen database, comparing the expression of DE genes in our dataset across multiple immune cell type-specific datasets. We found that DE genes that were upregulated in PWD CD4 T cells tended to have lower expression in T cells, and higher expression in non-T cells, e.g., myeloid lineage and stromal cells, whereas the downregulated genes tended to have a more T cell-like expression signature (Fig. 2E, top). The same was true for APCs, where upregulated transcripts in PWD tended to be expressed by non-myeloid/innate immune cells (e.g. T cells), and downregulated transcripts had a myeloid/innate immune-like signature (Fig. 2E, bottom). A global quantitative expression analysis supported these observations, revealing significant differences in lineage-specific gene expression between genes that were upregulated in PWD relative to B6 vs. those genes that were downregulated, typically in opposite directions across different cell lineages, e.g. higher expression of genes upregulated in PWD CD4 cells by innate immune cell lineages, and higher expression of downregulated in PWD CD4 cells by alpha-beta T cell lineage genes (Fig. 2F). This pattern also held true for other cell types, where downregulated transcripts in PWD cells typically belonged to the corresponding cell type, while the upregulated transcripts tended to be expressed by other cell types (data not shown). Additionally, this is supported by some significant enrichment of transcripts upregulated (but not downregulated) in PWD cells within the ASD brain transcript data set (Supplementary Fig. 2E). Altogether, these results demonstrate that PWD immune cells exhibit more promiscuous cell type-specific gene expression profiles, upregulating genes that are typically expressed by other cell types at the expense of cell type-specific genes.

Analysis of sex-specific gene expression and sex-by-strain interactions reveals minimal impact of sex on gene expression

The incidence and prevalence of many autoimmune diseases, such as MS, RA and SLE, exhibit a profound sexual dimorphism, with females being affected 3–10 times more often than males [22, 23]. The reasons for this are unclear, but it is thought that sex hormones and sex chromosomes influence gene expression in immune cells, which gives rise to sexual dimorphism in autoimmunity [22, 23]. To test this idea, and to see how it interacts with genetic control of gene expression, we compared the transcriptomes of immune cells isolated from male and female B6 and PWD mice. We first sought to identify genes that exhibited differential sex-specific expression as a function of strain (sex-by-strain interaction), i.e. those genes that exhibited a significantly different male:female (M:F) expression ratio in B6 vs. PWD (see Materials and Methods). Surprisingly, even using a relatively relaxed filter of FDR < 0.05, and |FC| > 1.5 (here the FC is in M:F ratio between PWD vs. B6), this analysis identified only two unique genes across all 5 cell types, Xist and Kdm5d (Fig. 3), encoded on the X and Y chromosomes, respectively, and well-known to exhibit sexually dimorphic expression (SDE) [24, 25]. While these two genes exhibited SDE in both strains (see below), Xist exhibited ~ 2.3 fold higher M:F ratio in PWD compared with B6, while Kdm5d exhibited ~2 fold lower M:F ratio (Fig. 3B and C). This pattern held true for these two genes across all 5 cell types, but did not reach the significance threshold in all cell types due to variability (data not shown). Lowering the stringency of our filter further (|FC| > 1.5, nominal p < 0.01) did identify a few more genes exhibiting strain-by-sex interaction (Supplementary Table 1), but given their low statistical significance, their impact is unclear. These results suggest that natural genetic variation exerts either a subtle or highly variable influence on sexually dimorphic gene expression in the immune system, which is in stark contrast to the profound genetic influence on the transcriptome independent of sex (e.g., Fig. 1).

Figure 3

Genes exhibiting significant sex-by-strain interaction across five cell types

(A) Genes exhibiting sex-by-strain interactions were identified as outlined in the Materials and Methods. Genes passing the filter of |FC| > 1.5 and FDR < 0.05 are shown. FC represents the change in male:female ratio between PWD and B6, calculated as (FC PWDMale – PWDFemale) – (FC B6Male – B6Female), see Materials and Methods. Thus, a positive value is indicative of more male biased expression of a gene in PWD compared to B6. (C) and (D) Relative expression values of the indicated genes exhibiting significant sex-by-strain interaction in CD8 cells. Relative expression values represent log2-scaled normalized raw expression values. Error bars indicate standard error of the mean.

Next, we sought to identify transcripts whose expression was sexually dimorphic, independent of genetic background. Since genetic background exerted little interaction with SDE (see above), B6 and PWD data were pooled for this analysis. Using a filter of |FC| > 1.5 and FDR < 0.05, this analysis identified 4 genes exhibiting SDE: Y-encoded Kdm5d, Eif2s3y, Ddx3y, and X-encoded Xist (Table 1), all well-known to be expressed in a sexually dimorphic fashion. The SDE of these genes was similar across different cell types. Pooling the data from all 5 different cell types identified 2 additional genes exhibiting SDE: hemoglobin genes Hba-a1 and Hbb-b1. Further lowering the stringency of the filter to |FC| > 1.5 and nominal p <0.01 identified 16 additional genes: 2 on the X chromosome (Utx and Alas2), and the rest on autosomes (Supplementary Table 2). Interestingly, most of these genes exhibited SDE only in one cell type, and most tended to have higher expression in males.

Table 1

Genes exhibiting significant SDE

FC (M:F) is shown for those genes reaching the cutoff filter (|FC| > 1.5, FDR < 0.05) for a given cell type. “All” indicates a combined analysis of all 5 cell types. Direction and strength of FC is also indicated by color (blue, higher expression in males; red – higher expression in females). An absence of a FC value indicates a failure of a gene to reach the cutoff for a given filter. Chr, Chromosome.

gene	Chr	All	CD4	CD8	Treg	APC	BC
Eif2s3y	Y	16	14	18	14	16	21
Kdm5d	Y	1.7		1.8		1.5	2.2
Xist	X	−4.8		−6.5		−6.5	−7
Ddx3y	Y	3.1				2.9
Hbb-b1	7	1.9
Hba-a1	11	2.0

PWD mice display resistance to EAE and altered associated immune responses

Our gene expression results above suggested a dampened basal immune activation state in PWD mice compared to B6, as well as lower expression of MS susceptibility genes. This led us to hypothesize that this would result in decreased susceptibility to experimental autoimmune encephalomyelitis (EAE), the principal autoimmune model of MS. To test this hypothesis, B6 and PWD mice were immunized with mouse spinal cord homogenate (MSCH) in complete Freund’s adjuvant (CFA), together with pertussis toxin as an ancillary adjuvant. The primary EAE readout, cumulative disease score (CDS), differed significantly by strain independent of sex (two way ANOVA, sex, F = 1.2, p = 0.47; strain, F = 21.3, p = 0.004; sex*strain interaction, F = 0.29, p = 0.72), therefore EAE data for males and females were pooled by strain. Compared with B6, PWD mice were highly resistant to EAE, as illustrated by reduced disease incidence, CDS, and other EAE quantitative trait variables (Fig. 4A–F).

Figure 4

PWD mice are resistant to EAE

EAE was induced and evaluated B6 (female, n=6; male n=10) and PWD (female, n=5; male n=15) mice as described in the Materials and Methods. The following EAE quantitative traits were calculated: (A) incidence, (B) cumulative disease score (CDS), (C) severity index, (D) peak score, (E) days affected, and (F) day of onset. Significance of differences in (A) was determined by Fisher’s exact test. Significance of differences in (B–F) was determined by two-tailed Student’s t-test. Significance of differences between B6 and PWD is indicated using asterisks as follows: *, p < 0.05; **, p<0.01, ***, p<0.001. Error bars indicate standard error of the mean. The data are pooled from two independent experiments, both of which yielded similar results.

We next tested whether the relevant encephalitogenic T cell responses were affected in PWD mice. EAE and MS are thought to be initiated and driven by CNS autoantigen-reactive CD4 T cells of the Th1 or Th17 phenotype, identified by their signature cytokines, IFNγ and IL-17, respectively [26]. GM-CSF is another cytokine that can be produced by either Th1 or Th17 cells, and its expression correlates with their encephalitogenic potential. In contrast, FoxP3+ Treg cells are immune-regulatory in EAE. Therefore, we examined the expression of these 3 signature cytokines and the frequency of FoxP3+ Treg cells. B6 and PWD mice were immunized with MSCH as above, and T cell responses in the spleen and draining lymph nodes were assessed by flow cytometry and intracellular staining. We found that in the spleen, compared with B6, PWD mice had more GM-CSF+ CD4 T cells, and comparable numbers of IFNγ and IL-17 producers (Fig. 5A). Interestingly, PWD CD8 T cells in the spleen produced significantly lower amounts of IFNγ, but higher amounts of IL-17 (Fig. 5B). In contrast, in the draining lymph nodes, PWD CD4 and CD8 T cells produced much lower amounts of all three cytokines compared with B6 (Fig. 5C and D). Treg frequency largely followed the magnitude of the effector T cell responses, with PWD mice having more FoxP3+ Treg in the spleen, but fewer in the draining lymph nodes (Fig. 5E). These findings suggest that PWD mice are capable of mounting a potent T cell response, but it is weaker in the lymph nodes compared to spleen, where B6 mice exhibit a much more robust T cell response. The reduced effector T cell responses in PWD do not appear to be due to an enhanced Treg expansion, since the FoxP3+ Treg frequency is proportionally to the effector T cell responses in both strains.

Figure 5

PWD mice display skewed peripheral immune responses

B6 (n=6) and PWD (n=9) mice were immunized as in Fig. 3. At day 10 post-immunization, cells were isolated from the spleen (A, B, E) and draining lymph nodes (C, D, E), restimulated with PMA/Ionomycin (except in E), stained for surface markers, followed by fixation and intracellular staining for the indicated cytokines (A–D) or FoxP3 (E), and flow cytometric analysis. Percentages of cytokine positive cells among the live CD19-TCRβ+CD4+ (A, C, E), or CD19-TCRβ+CD8+ (B, D) populations are shown. Significance of differences was determined by two-way ANOVA with Bonferroni’s multiple comparison test. The data represent one independent experiment.

We next examined the immune response in the relevant target organ for EAE, the CNS. Mice were immunized with MSCH as above, and immune cells were isolated from the CNS on D30 post-EAE induction, and analyzed by flow cytometry. Compared with B6, PWD mice had a profound reduction in the infiltration of immune cells into the CNS, as demonstrated by reduced numbers of CD45+ and TCRβ+ cells (Fig. 6A and B). In addition to the reduced CNS T cell numbers in PWD mice, a lower proportion of CD4 T cells in the CNS produced IFNγ and IL-17 (Fig. 6C). Collectively these results suggest that the EAE resistance of PWD mice is associated with altered and/or reduced encephalitogenic T cell responses in peripheral lymphoid organs, and the inability of these T cells to enter the CNS efficiently.

Figure 6

PWD mice display reduced immune responses in the CNS during EAE

B6 (n=4) and PWD (n=6) mice were immunized as in Fig. 3. At day 30 post-immunization, mononuclear cells were isolated from the CNS by Percoll gradient, counted, and enumerated by flow cytometry. Numbers of cells positive for the indicated markers were calculated by multiplying the total number of isolated mononuclear cells by the percentage of CD45+ cells (A) or by the percentage of CD45+TCRβ+ cells (B). In (C), mononuclear cells were restimulated with PMA/Ionomycin and analyzed by intracellular staining and flow cytometry, as in Fig. 4. Percentages of CD45+TCRβ+CD4+ cells positive for the indicated cytokines are shown. Significance of differences was determined using the Student’s t-test. The data represent one independent experiment.

Discussion

While standard inbred laboratory strains of mice exhibit a low level of genetic diversity, the wild-derived PWD strain is highly divergent compared with the standard B6 strain [9, 10], similar perhaps to the genetic differences between ethnically distinct human populations such as, for example, Europeans, Africans, or Asians. Importantly, our mouse model eliminates environmental factors that have profound influences on gene expression in human populations [27], and allows for the study of genetic control only. Our results demonstrate that the level of this genetic control over gene expression in immune cells is profound, with thousands of genes differentially expressed between B6 and PWD strains at baseline, some cell-specific and others conserved across different cell types. Of note, some of the genes highly upregulated in PWD cells compared to B6 included anti-viral or interferon-induced genes, such as Mx1 and Mx2, several genes encoding IFITM (interferon-induced transmembrane) and IFIT (interferon-induced proteins with tetratricopeptide repeats) proteins, (with the notable exception of Ifi27, which was highly downregulated in PWD), which may reflect a loss of evolutionary pressure exerted by viral infection in the laboratory B6 strain (Fig. 1B and Supplementary file 1). Moreover, laboratory strains of mice, unlike wild-derived mice, carry a non-functional and often poorly expressed alleles of Mx1 and Mx2, which results in their heightened susceptibility to influenza virus [28-30]. This is consistent with the low expression of these two genes in B6 compared with PWD cells in our data set.

Another group of genes which were highly upregulated in PWD cells were Actb and Actn2, encoding beta actin and the actin-regulatory protein actinin alpha 2, respectively (Fig. 1B). This may reflect the requirement for a higher cytoskeleton-dependent mobility for leukocytes in PWD mice, again potentially selected by higher evolutionary pressure exerted by pathogens. Additionally, the cytoskeleton and its associated proteins are an important regulator of intracellular signaling and immune effector function (beyond motility) in lymphocytes and other immune cells [31-33].

Several DE genes are encoded by the mitochondrial genome, with lower expression in PWD cells. This may reflect different metabolic states or requirements in the two different strains, suggesting a lower mitochondria-dependent metabolic demand in the naïve state of PWD immune cells, which may be important in the conservation of the available resources that may be limited due to scarce food sources in the wild. Notably, the dynamics of metabolic state(s) in immune cells has recently emerged as a critical regulator of immune cell effector function [34, 35]. The DE of mitochondrial genes is also consistent with the idea that the non-recombining nature of the mitochondrial genome is likely to result in highly divergent genome sequence and expression profiles between these two distantly related strains.

The ImmGen Consortium recently published a large microarray-based study profiling gene expression in two immune cell types, CD4 T cells and granulocytes, across a panel of 39 inbred strains of (male) mice, including several wild-derived strains [36]. The results of this analysis are in line with ours, demonstrating a high level of genetic control over gene expression, and widespread variation in gene expression across different strains, with the largest differences seen between wild-derived and conventional laboratory strains. Interestingly, many of the most profoundly DE genes in wild-derived mice in the ImmGen study match ours, e.g., Ifitm1, Ifitm2, Ifi27, Cd163, Klrd1, Anxa3, Cd59a, Chi3l3 (Fig. 1B), as well as Tlr1 and Tlr7 (Supplementary File 1).

Our published work utilizing B6.ChrPWD consomic strains to map EAE QTL revealed striking sex differences in the genetic control of EAE [14]. Based on these differences, we expected to find large sex differences in immune cell transcriptomes. However, in stark contrast to the dramatic effect of genetic background on gene expression, we found very few genes exhibiting significant SDE or sex-by-strain interactions, with the majority of these localized on the sex chromosomes. Interestingly, Kdm5d, an ancestral single-copy gene which resides on the short arm of mouse chromosome Y [37], was the only gene found to exhibit significant strain-by-sex interaction. This highlights Kdm5d as a potential candidate gene responsible for the EAE phenotype in B6.ChrYPWD consomic mice, which display augmented EAE susceptibility [14, 38]. Collectively, these results are consistent with previous reports of SDE across non-immune tissues, where relatively subtle differences in SDE of autosomal genes have been detected [24, 25, 39]. Taken together, our results suggest that genetic background has a profound influence on gene expression, but its influence on SDE is relatively subtle. This is also in line with the findings from human GWAS, e.g., in MS, where no sex-specific autosomal candidate genes have been reported to date, although notably, loci on sex chromosomes have not been included for technical reasons [6, 40–42]. Given the dramatic phenotypic differences between the sexes in immunity [23], these findings are surprising, yet they suggest that sex may exert a relatively minor influence on transcript expression levels in the naïve state of immune cells. It is possible that more profound differences in gene expression are observed after immune activation. It is also possible that there is a higher level of sex-specific control at post-transcriptional levels, which is supported by the robust SDE of two Y-linked translation regulator genes, Eif2s3y and Ddx3y, observed in our study (Table 1) and many others [37]. Both of these genes, like the other handful of Y-linked ancestral single-copy ubiquitously expressed genes (2 to 4 other genes in the mouse, including Kdm5d), have been proposed to function as dosage-sensitive regulators of gene expression, translation, and protein stability, and as such likely play essential roles in male viability, development, and sexual dimorphism in health and disease far beyond male gamete function and sexual differentiation [37, 43]. Moreover, the notion that bigger sex differences can be seen at the protein level is also supported by a recent systems proteomics approach, which revealed that sex is a major factor in determining protein levels of autosomal genes [44].

The recent explosion in GWAS has identified hundreds of genetic variants associated with complex polygenic diseases, including autoimmune diseases. Integrating these data with expression QTL (eQTL) studies (see below) has suggested that many GWAS candidates modify autoimmune susceptibility by driving differential expression in autoimmune disease [45-47]. Using MS as a prototypical autoimmune disease, we show enrichment of MS-GWAS genes in our DE gene sets, with striking directionality: only those genes that are exhibit lower expression in PWD are significantly enriched with MS-GWAS genes. This suggests that the B6-PWD genetic model appropriately models natural genetic variation that is relevant to human autoimmune disease. It also predicts reduced susceptibility of PWD mice to CNS autoimmunity, a hypothesis that is supported by our functional data using the EAE model. It is important to note that most of the candidate MS-GWAS genes are typically identified by imputation analysis using the nearest SNP marker with a significant effect, therefore it is likely that not all current candidates represent the true MS genes, and improved fine-mapping and candidate gene identification continues to be a work-in-progress [6]. Nonetheless, the results from our model support the functional importance of at least a majority of the current GWAS candidates included in our analysis. Future studies can include similar analyses of our data sets using emerging results from follow-up GWAS and fine-mapping studies, which should prove informative.

Several recent seminal studies in humans have examined the effect of natural genetic variability on gene expression in adaptive and innate immune cells in large cohorts of genetically diverse individuals [45-48]. Thousands of cell type-specific and non-cell type-specific expression eQTL were identified, some of which were restricted to specific ethnic groups, and others that were shared across ethnic groups. While it is difficult to compare our results to these studies directly, it is clear that natural genetic variation exerts a strong influence on gene expression in immune cells in both humans and mice. In the PWD:B6 mouse comparison, this influence is very robust, since we are able to eliminate variability introduced by environmental influences and heterogeneous genetic backgrounds in the human studies [27]. This also highlights the utility of the mouse model in studying gene-by-environment interactions in a setting where the genetic and environmental factors can be tightly controlled and manipulated, to support or refute cause-effect relationships, which are more challenging to assess in human studies. Such future studies will complement human studies, providing a better mechanistic understanding of the genetic basis of complex diseases and their environmental modulators.

Materials and Methods

RNA isolation and microarray analysis of RNA expression

For the microarray analysis on basal expression differences in immune cell subsets, three biological replicates for each strain and sex combination were created by pooling cells from three different individual naïve 8–10 week old mice into each biological replicate. Cells were isolated from Liberase/DNase I-digested spleens and combined with total cells from lymph nodes (axillary, brachial, and inguinal) for each mouse. B cells were isolated using the EasySep B cell positive selection kit and EasySep magnet (STEMCELL Technologies, Inc., Canada). The remaining live cells from the flow through were purified by fluorescently activated cell sorting using fluorophore conjugated antibodies against cell surface markers as follows: CD4 T cells (CD45+NK1.1−CD19−TCRβ+CD4+CD25−); CD8 T cells (CD45+NK1.1−CD19−TCRβ+CD8+); Treg cells (CD45+NK1.1−CD19−TCRβ+CD4+CD25+); APCs (CD45+NK1.1−CD19−TCRβ−CD11b+CD11c+). Antibodies were purchased from BioLegend, San Diego, CA, USA; catalog were numbers as follows: CD45, NK1.1, CD19, CD4, CD25, CD8, CD11b, CD11c, TCRβ; 103112, 108707, 115534, 100531, 102016, 101206, 117319, 109222, respectively.. High quality RNA was isolated using the Qiagen RNeasy Plus Mini Kit and the transcriptomes analyzed using Illumina BeadArray technology, using 45,281 unique probes.

Illumina BeadArray

100 ng of RNA was amplified and converted to cRNA using Illumina TotalPrep-96 RNA Amplification kit (Ambion, Carlsbad, CA, USA). 500 ng of cRNA was used for hybridization onto the Illumina Whole-Genome Gene Expression Direct Hybridization microarray (Illumina MouseWG-6 v2.0 R2 Expression BeadChip; Illumina, San Diego, CA, USA). 45,281 probes from the microarray were included for analysis. All microarray data were uploaded to Gene Expression Omnibus, under accession number GSE85418.

Statistical analyses of microarray data

Probe-level intensities were calculated using the lumi and limma packages in R specifically for Illumina arrays (http://www.basic.northwestern.edu/publications/lumi/lumi.pdf), including background-correction and quantile normalization for each probe set and sample. Summarized intensity data were imported into Partek Genomics Suite®, version 6.6 (Copyright © 2009, Partek Inc., St. Louis, MO, USA) for multivariate and univariate analyses. Principal Component Analysis (PCA), using the covariance matrix, was performed to 1) look for outlier samples that would potentially introduce latent variation into the analysis of differential expression across sample groups, and to 2) assess sample-based differential expression within and between sample groups. One outlier sample (female B6 Treg) was identified by PCA and excluded from all analyses. Univariate linear modeling of sample groups was performed using ANOVA as implemented in Partek Genomics Suite. The magnitude of the response (fold change (FC) calculated using the least square mean) and the p-value associated with each probe set and binary comparison are calculated, as well a “step-up,” adjusted p-value for the purpose of controlling the false discovery rate (FDR) [49]. In all analyses, the FDR was at least five times larger than the nominal/uncorrected p-value. For downstream analysis, e.g., identification of number of DE genes, pathway analysis, etc. multiple probes probing for the same gene were averaged.

For identification of sex-by-strain interactions, to identify genes where male:female expression ratio significantly different between B6 and PWD cells, the following comparison was made for each cell type, to calculate FC, p-values, and FDR: (PWDMale minus PWDFemale) minus (B6Male minus B6Female), which is algebraically equivalent to (PWDMale and B6Female) minus (B6Male and PWDFemale). FC here represents the fold change in male:female ratio between PWD and B6.

Bioinformatic analyses

Pathway analysis was performed using Ingenuity Pathway Analysis (IPA) software. The expression dataset for all 5 cell types was uploaded into IPA and filtered by |FC| > 2 and FDR < 0.05, then subsequently analyzed using the Core Analysis function in IPA, followed by the Comparison Analysis function to compare across the 5 cell types, as follows. The Canonical Pathway function was used to identify the top canonical pathways (p < 0.01, Z score > 0.5) affected by the DE genes between B6 and PWD. The Upstream Analysis function was similarly used to identify top upstream predicted regulators (p < 0.01, Z score > 2). Top twenty pathways (ranked by Z-score) were shown.

Enrichment of DE genes in the MS-GWAS candidate gene list was performed in IPA software as follows. The current published best list of MS-GWAS candidate genes [15] was imported into IPA. The Core Analysis function was used to determine the significance of enrichment of DE genes (up- or downregulated separately) within the MS-GWAS list. The same procedure was carried out on the following data sets/gene lists: 1) a list of transcripts reported to be upregulated in CD4 T cells from MS-CIS subjects vs. controls [20], 2) a GWAS candidate gene set for IBD [16], 3) a set of candidate genes for ASD [18, 19], and 4) a set of genes differentially expressed in ASD brains relative to healthy control brain [18, 21].

Cell type specific gene set enrichment analysis was performed using the ImmGen database, using MyGeneSet function (http://rstats.immgen.org/MyGeneSet/). The top 200 (ranked by FC; FDR < 0.05) upregulated genes in PWD (relative to B6) were used in the W Plot function, then the same procedure was repeated for the top 200 downregulated genes. To generate quantitative comparisons of enrichment, the expression of the top 200 up- and downregulated genes for each of the five cell subsets was analyzed across the ImmGen cell type-specific data set (version 1). A global average was obtained by averaging the expression of all 200 genes for a given ImmGen subtype, then by determining the average of all of these cell subtypes within a specific lineage/category, e.g. monocytes. This average expression thus serves as a quantitative measure of enrichment of gene expression within a particular ImmGen population, and this measure was compared between 200 top genes upregulated in PWD relative to B6 vs. 200 top genes downregulated in PWD, for each of the 5 cell types analyzed in our study. Significance of differences was determined in GraphPad Prism, using the two-stage linear step-up procedure of Benjamini, Krieger and Yekutieli, with Q = 5%.

Animals and Induction and evaluation of EAE

C57BL/6J and PWD/PhJ mice were purchased from Jackson Laboratories (Bar Harbor, Maine, USA) and bred and housed in the vivarium at the University of Vermont. The experimental procedures used in this study were approved by the Animal Care and Use Committee of the University of Vermont.

EAE was induced in male and female B6 and PWD mice as follows. Mice were injected subcutaneously with 0.1 ml of emulsion containing 2.5 mg of MSCH in PBS and 50% CFA (Sigma, USA) on day 0 and day 7. CFA was supplemented with 4 mg/ml Mycobacterium tuberculosis H37Ra (Difco, USA). On day 0 and day 2 mice also received an i.p. injection of 200 ng pertussis toxin (List Laboratories, USA) as an ancillary adjuvant. Starting on day 10, mice were scored visually as previously described [50]. Briefly, the clinical scores were as follows: 0.5 - partial loss of tail tone, 1 – full loss of tail tone, 2 - loss of tail tone and weakened hind limbs, 3 - hind limb paralysis, 4 - hind limb paralysis and incontinence, 5 - quadriplegia or death. EAE scoring was not performed in a blinded fashion, since B6 and PWD mice are visually distinct. EAE quantitative traits were calculated as previously described [51], as follows. The incidence of EAE was recorded as positive for any mouse with clinical signs of EAE for 1 or more consecutive days. Cumulative disease score (CDS) was calculated as the sum of all daily scores over the course of 30 days. Days affected was calculated as the number of days an animal displayed a clinical score > 0, and day of onset was the day a clinical score > 0 was first observed (not calculated for animals without clinical signs). Severity index (assessed in affected animals only) was generated by averaging the clinical scores for each animal over the number of days that it exhibited clinical symptoms. Peak score represents the maximum daily score.

Flow cytometry

For intracellular cytokine staining ex vivo, mice were immunized for EAE induction as above. Spleen and draining (for the immunization site) lymph nodes (axillary, brachial, and inguinal) were harvested on day 10 post-immunization, and cells were stimulated with 5 ng/ml of PMA, 250 ng/ml of ionomycin (Sigma-Aldrich, USA) and brefeldin A (Golgi Plug reagent; BD Biosciences) for 4 hours. Cells were then stained with the UV-Blue Live/Dead fixable stain (Invitrogen, USA) and then surface stained for the following markers: CD4, CD8, and TCRβ. Cells were then fixed with 1% paraformaldehyde (Sigma-Aldrich, USA), permeabilized with buffer containing 0.2% saponin and stained with anti-IL-17A, anti-IFNγ, and anti-GM-CSF (Biolegend, USA).

For surface marker analysis and Foxp3 staining, unstimulated isolated cells were stained directly ex vivo with the UV-Blue Live/Dead fixable stain and then surface labeled for different combinations of following markers: CD25, CD19, CD4, CD8, and TCRβ (Biolegend, USA) and fixed using the Foxp3 fixation/permeabilization buffer (eBioscience, USA), followed by intracellular staining for Foxp3. Antibodies used for flow cytometry were directly conjugated to fluorophores and obtained commercially (Biolegend, USA, catalog numbers were as follows: CD19, CD4, CD25, CD8, CD11b, TCRβ, IL-17A, IFNγ, GM-CSF; 115534, 100531, 102016, 101206, 109222, 506904, 5050813, 505404, respectively. Anti-FoxP3 antibody was purchased from eBiosciences (USA), catalog number 12-5773-82.

Labeled cells were analyzed using an LSR II cytometer (BD Biosciences). Compensation was calculated using appropriate single color controls. Data were analyzed using FlowJo software (Tree Star Inc, Ashland, OR).

CNS-infiltrating mononuclear cell isolation

Animals were perfused with PBS and brains and spinal cords were removed. A single cell suspension was obtained and passed through a 70 μm strainer. Mononuclear cells were obtained by Percoll gradient (37%/70%) centrifugation and collected from the interphase. For intracellular cytokine analysis, cells were washed and stimulated with 5 ng/ml of PMA, 250 ng/ml of ionomycin the presence of brefeldin A (Golgi Plug reagent, BD Bioscience) for 4 hours. Cells were labeled with the UV-Blue Live/Dead fixable stain (Invitrogen) followed by surface staining (CD45, CD11b, CD4, CD8, TCRγδ, and TCRβ). Afterwards, cells were fixed, permeabilized and stained for intracellular IL-17A, IFNγ, and GM-CSF as described above. Alternatively, unstimulated CNS cells were surface labeled for CD45, CD11b, TCRβ, CD4, CD8, then fixed and stained for FoxP3, as above.

General Statistical Analyses

Statistical analyses not pertaining to microarray data were carried out using GraphPad Prizm software, version 6. Details of the analyses are provided in the Figure Legends. All statistical tests were two-sided, and adjustments for multiple comparisons were made as indicated. All center values represent the mean, and error bars represent the standard error of the mean. P-values below 0.05 were considered significant. Sample sizes for animal experiments were chosen based on previous experience with similar analyses. No randomization was used to assign animals to different treatment groups since no differential treatment was performed between the two different strains or the two sexes.

Supplementary Table 1. Additional genes exhibiting significant sex-by-strain interaction across five cell types. Genes exhibiting sex-by-strain interactions were identified as outlined in the Materials and Methods. Genes passing the filter of |FC| > 1.5 and nominal p < 0.01 are shown.

Supplementary Table 2. Additional genes exhibiting significant SDE. FC (M:F) is shown for those genes reaching the cutoff filter (|FC| > 1.5, nominal p < 0.01) for a given cell type. “All” indicates a combined analysis of all five cell types. Direction and strength of FC is also indicated by color (blue, higher expression in males; red – higher expression in females). An absence of a FC value indicates a failure of a gene to reach the cutoff for a given filter. Chr, Chromosome.