Molecular Phenotyping of Immune Cells from Young NOD Mice Reveals Abnormal Metabolic Pathways in the Early Induction Phase of Autoimmune Diabetes

Islet leukocytic infiltration (insulitis) is first obvious at around 4 weeks of age in the NOD mouse – a model for human type 1 diabetes (T1D). The molecular events that lead to insulitis and initiate autoimmune diabetes are poorly understood. Since TID is caused by numerous genes, we hypothesized that multiple molecular pathways are altered and interact to initiate this disease. We evaluated the molecular phenotype (mRNA and protein expression) and molecular networks of ex vivo unfractionated spleen leukocytes from 2 and 4 week-old NOD mice in comparison to two control strains. Analysis of the global gene expression profiles and hierarchical clustering revealed that the majority (∼90%) of the differentially expressed genes in NOD mice were repressed. Furthermore, analysis using a modern suite of multiple bioinformatics approaches identified abnormal molecular pathways that can be divided broadly into 2 categories: metabolic pathways, which were predominant at 2 weeks, and immune response pathways, which were predominant at 4 weeks. Network analysis by Ingenuity pathway analysis identified key genes/molecules that may play a role in regulating these pathways. These included five that were common to both ages (TNF, HNF4A, IL15, Progesterone, and YWHAZ), and others that were unique to 2 weeks (e.g. MYC/MYCN, TGFB1, and IL2) and to 4 weeks (e.g. IFNG, beta-estradiol, p53, NFKB, AKT, PRKCA, IL12, and HLA-C). Based on the literature, genes that may play a role in regulating metabolic pathways at 2 weeks include Myc and HNF4A, and at 4 weeks, beta-estradiol, p53, Akt, HNF4A and AR. Our data suggest that abnormalities in regulation of metabolic pathways in the immune cells of young NOD mice lead to abnormalities in the immune response pathways and as such may play a role in the initiation of autoimmune diabetes. Thus, targeting metabolism may provide novel approaches to preventing and/or treating autoimmune diabetes.


Introduction
The earliest signs of autoimmune pathology in the pancreas of the NOD mouse (a model for human type 1 diabetes, T1D) occur at approximately 4 weeks of age [1]. At this stage accumulation of leukocytes is observed around the pancreatic islets. This accumulation (known as insulitis) progressively intensifies and becomes increasingly more invasive eventually leading to massive destruction of the insulin producing beta cells with clinical diabetes occurring at 12 weeks of age or later [2,3]. The molecular events leading to insulitis and that may influence initiation of the autoimmune pathology are poorly understood.
It is well established that genetic predisposition plays a major role in the etiology of TID. The strongest genetic determinant for T1D development is the MHC class II gene product, H2 g7 in NOD mice or HLA-DQ2 and -DQ8 in humans [2,4,5]. It is clear that MHC class II-restricted CD4 T-cells play an important role in inducing T1D [6][7][8][9]. In addition to the MHC locus, which is required for T1D development, over 20 non-MHC loci also contribute to the disease process [5,[10][11][12][13]. The molecular mechanisms by which most of these genes contribute to disease have not been elucidated. The focus in the field has been mainly to study candidate genes one by one. However, this approach has provided limited information on how the multiple genes interact to cause disease.
Since TID is caused by numerous genes, we hypothesized that multiple molecular pathways are altered and interact to initiate autoimmune pathology. To test this hypothesis, we conducted a detailed evaluation of the molecular phenotype (mRNA and protein expression) of untreated unfractionated spleen leukocytes from NOD mice in the pre-insulitis (early induction) phase of autoimmune diabetes. Applying a modern suite of multiple bioinformatics approaches, we identified abnormal molecular pathways that might play a role in the initiation of autoimmune pathology in NOD mice.

Ethics Statement
All the procedures involving animals were approved by the University of Tennessee (UT) Health Science Center and Veteran Affairs (VA) Medical Center Animal Care and Use Committee (Protocol Numbers: UT 1159/VA 00157). Breeder mice were purchased from the Jackson Laboratory. Animals were housed at the VA animal facility.

Sample collection
The procedures for collection of spleen leukocytes, and extraction of RNA and protein samples have been previously described [14]. Spleen leukocytes were collected from female NOD mice and females of two control strains, NON, and C57BL/ 6, at 2 weeks and 4 weeks of age (n = 5 for each strain and each age group, 30 samples in total). NON is the original diabetes-resistant control strain closely related genetically to the NOD strain, whereas C57BL/6 is a more distantly related strain with some immunogenetic similarities to NOD. The isolated cells were not subjected to any post-collection treatment or culture prior to RNA/protein extraction.

Transcriptome data collection
RNA expression was analyzed on Affymetrix expression arrays MOE430A and MOE430B as described previously [14]. Expression levels were evaluated using the Microarray Suite version 5.0 (MAS 5.0) using Affymetrix default analysis settings and global scaling as the normalization method (www.affymetrix.com). The trimmed mean target intensity of each array was arbitrarily set to 100. The two arrays were combined prior to statistical analysis to give a total of ,45,000 probe sets. Microarray data has been submitted to the GEO repository (Accession number: GSE37450). The microarray results were validated by quantitative Real-time PCR (qRT-PCR) using Roche's Universal library probe system (www.roche-applied-science.com). First strand cDNA was synthesized from 1 mg of RNA using Transcriptor First Strand cDNA synthesis kit (Roche) according to the manufacturer's instructions. The target genes, primers and probes are listed in Table S1. Expression values were normalized using glyceraldehyde-3phosphate dehydrogenase (Gapdh).

Proteome data collection
Protein expression was analyzed by 2D PAGE as decribed previously [14]. The protein gels were scanned and analyzed as previously described [15,16], except instead of using the PDQuest image analysis software to create the Match Set we used Nonlinear software (Progenesis). We identified a total of 382 protein spots as being expressed in at least one of the 30 gels.

Statistical analysis
We imported and further analyzed the RNA and protein expression datasets in GeneSpring software (version 7.3.1, Silicon Genetics, Redwood, CA) as previously described [14]. The datasets consisted of 15 samples for each age, 2 weeks and 4 weeks. The 45,000 probe sets present on the Affymetrix expression arrays were filtered to identify 24,959 that had a present/marginal expression flag in at least one of the samples. The lists of filtered probe sets and the total protein spots were subjected to statistical and hierarchical clustering analysis. In order to define lists of age-specific differentially expressed genes between strains, we performed one-way ANOVA on the expression datasets by age. For the transcriptome datasets, we applied two cut-off levels, p,0.005, p,0.05, both with Benjamini-Hochberg multiple test correction. The proteome datasets were analyzed at one level p,0.05, with no multiple test correction (MTC). The impact of potential false positives that may arise from not applying a MTC on the final conclusions drawn from these data would be minimized by downstream data mining analyses and comparison with the transcriptome data. We conducted hierarchical clustering on the resulting lists to identify genes or spots that were differentially expressed only in the NOD mice compared to both control strains. We identified proteins in differentially expressed spots by mass spectrometry according to methods previously described [15].

Data mining analyses
We subjected the lists of NOD differentially expressed mRNAs or proteins to data mining analyses using a suite of modern bioinformatics tools. The genes/proteins were grouped into functional classes to identify enriched gene ontology (GO) classes using WebGestalt Gene Set Analysis Toolkit (bioinfo.vanderbilt.edu/webgestalt) [17]. WebGestalt was also used for KEGG pathways analysis. Prior to uploading the lists into WebGestalt, the probe sets and proteins identifiers (ID) were converted to entrezgene ID using the gene ID conversion tool available at the DAVID (Database for Annotation, Visualization and Integrated Discovery) server (david.abcc.ncifcrf.gov) [18,19]. To identify molecular networks that are associated with the differentially expressed genes, we analyzed the gene lists using Ingenuity Pathway Analysis (IPA, www.ingenuity.com). The Affymetrix probe set IDs or GI numbers (for proteins) were uploaded into the server. Each identifier was mapped to its corresponding gene in the Ingenuity Pathways Knowledge Base (IPKB), a database compiled primarily from scientific literature. The mapped genes, called focus genes, were then algorithmically combined into molecular subnetworks together with other genes and/or endogenous chemicals from the IPKB (referred to as non-focus genes) based on information contained in the IPKB. The generated subnetworks are a graphical representation of the molecular relationships between genes/gene products and/or endogenous chemicals based on what is known about them from the literature. Herein, in the context of IPA, the term ''gene'' will be used for both genes/gene products as well as endogenous chemicals. Genes are represented as nodes of various shapes (shown in Figure S1) that represent the functional class of the gene product. The biological relationship between 2 nodes is represented as an edge (or line, also referred to in our study as a connection). Each connection is supported by at least one reference from the literature, a textbook, or from canonical information stored in the IPKB; and represents any known interaction between 2 genes (e.g. A induces B, C binds D, E inhibits F, and so on). Solid lines in the generated networks represent direct interactions (such as binding) while dashed lines represent indirect interactions (such as expression of gene A decreases expression of gene B). The Ingenuity package assesses the significance of each subnetwork with respect to random networks of the same size using the Fisher's exact test, and assigns a score based on the number of genes in the network and its size, the total number of genes analyzed, and the total number of genes in IPKB that could potentially be included in networks. We referred to the most highly scored preliminary networks as major networks. The biological processes and/or diseases most significantly associated to the genes in a particular network are identified using the IPA functional Analysis tool (Fisher's exact test). Because the Ingenuity server limits the number of genes per network to 35, bigger lists, like our transcriptome datasets, often create multiple networks. We think that this is an artificial division of the data which normally could be connected into one large network based on information available in the database. Therefore, we merged individual networks using the IPA ''merge networks'' tool. Due to the limitations of IPA (most important of which is the paucity or absence of information in the literature on particular genes) only a proportion of genes in the uploaded list may get incorporated into the generated network.

Results
Identification of differentially expressed mRNAs in spleen leukocytes from NOD mice at the preinsulitis stage of autoimmune diabetes To identify differentially expressed mRNAs and proteins in spleen leukocytes from NOD mice at the preinsulitis stage, we analyzed transcriptome and proteome expression in 2 and 4 weekold NOD mice. Each dataset was analyzed by age. We compared gene expression in NOD mice to that in two control strains, NON and C57. We were interested in genes that were differentially expressed in NOD relative to both controls (hereinafter referred to as NOD differentially expressed genes or mRNAs or proteins). We judged that this strategy would increase the likelihood that the identified genes are linked to initiation of the autoimmune pathology seen in NOD mice.
To identify differentially expressed mRNAs, we used expression microarrays and performed a one-way ANOVA on ,25,000 filtered probe sets. We identified 175 and 189 probe sets (genes) with highly significant strain differences (p,0.005, Benjamini-Hochberg multiple test correction) at 2 and 4 weeks, respectively. We then subjected these two gene lists to hierarchical clustering as previously described [20]. The hierarchical clusters of the 2 week ( Figure 1) and 4 week ( Figure 2) datasets each revealed six distinct patterns of expression: NOD low, NOD high, NON low, NON high, C57 low, and C57 high. At 2 weeks, a total of 51 genes (marked by black bars in Figure 1, and listed in Table 1) were differentially expressed in NOD mice compared to both control strains: 45 (88%) were of lower expression while 6 (12%) were of higher expression in NOD. At 4 weeks, a total of 76 genes (marked by black bars in Figure 2 and listed in Table 2) were differentially expressed in NOD mice compared to both control strains: 67 (88%) were of lower expression while 9 (12%) were of higher expression in NOD mice. Twenty-five genes (88% of which was of Figure 1. Hierarchical clustering of 175 genes with highly significant expression differences between strains at 2 weeks. The list was generated by one-way ANOVA of 24,959 probe sets (p,0.005, Benjamini-Hochberg multiple test correction). A total of 51 genes were differentially expressed in NOD mice relative to both control strains (NON and C57BL/6, C57): 45 were of lower expression (L) and 6 were of higher expression (H) in NOD mice. The color intensity of the rectangles representing each gene for each sample indicates the degree of increase (red) or decrease (blue) of the gene expression signal relative to the mean signal intensity (yellow). doi:10.1371/journal.pone.0046941.g001 lower expression in NOD mice) were common to both 2 and 4 week old mice. Microarray results were validated by qRT-PCR. Results confirming the microarray data for six genes (H60a, Trim5, Trim12, Ly6c1, Snx6 and Sp110) are shown in Figure S2. The major functional categories (as determined by us based on information obtained from literature searches) included immune response, apoptosis, ubiquitination, transporters, transcription and protein synthesis, and enzymes (including those involved in phosphorylation) at both ages. Several genes did not fall in one major group, thus were listed under ''Other''. Probe sets of unknown identification/function are listed under ''Unknown''. Of note, a number of these genes have been reported to be linked to T1D risk. For example, Ly6c1 and Trim 5 belong to a previously identified IRF7-driven expression network that is linked to TID risk [21]. Ly6c1, Trim 12 and Trim 30 are linked to a defect in T cell negative selection and have been previously identified as T1D candidate genes [22].
To identify genes likely associated directly with initiation of insulitis, we compared the trascriptomes of 2 week old NOD mice with those of 4 week old mice, with reference to corresponding comparative analyses of the control strains. Statistical analysis of the mRNA expression datasets by unpaired t-test (p,0.001, no MTC) followed by Venn diagram analysis identified 66 probe sets (representing 65 different genes, Table 3), 122 probesets (representing 120 genes, Table S2) and 105 probesets (representing 104 genes; Table S3) that were uniquely differentially expressed in spleen leukocytes of NOD, NON and C57BL/6 mice, respectively. Expression of ,62% [40] of the NOD-specific age differences was upregulated at 4 weeks compared to expression at 2 weeks of age (Table 3), in contrast to 45% and ,26% in NON and C56BL/6 mice, respectively. Unlike the gene lists of the control strains, the list of NOD mice consisted of several molecules involved in immune activation (e.g. Plcd4, Mapk4, Ifit3, Pdcd1lg2, Ifih1, IL17rb).
Identification of differentially expressed proteins in spleen leukocytes from NOD mice at the preinsulitis stage of autoimmune diabetes To identify differentially expressed proteins, we analyzed protein extracts that were simultaneously obtained from the same samples that generated the RNA as described previously [14]. Supplemental figure S3 shows an example of the 2D-gels that were  produced by the protein extracts. We conducted a one-way ANOVA on a total of 382 spots to identify proteins with statistically significant (p,0.05) expression differences between strains. We found that 19 and 26 spots had strain differences at 2 and 4 weeks, respectively, of which 11 and 8, respectively, were differentially expressed in NOD compared to both control strains.
Of the proteins differentially expressed in NOD only, we identified 9 at two weeks (6 lower, 3 higher in NOD) and 7 at 4 weeks (4 lower, 3 higher) ( Table 4). In some cases, we identified two proteins from a single spot (423, 441, 523), likely due to proximity in their respective molecular weights. In other cases, a single protein was identified from two different spots (Lmnb1 and Lcp1 from 423 and 532, and Hspa8 from 471 and 601), likely due to post-translational modifications of the protein. Although, no proteins were common to the proteome datasets of 2 and 4 week old mice, our functional categorization showed that apoptosis (the major functional category represented) was common to both age groups.

Integration of the transcriptome and proteome data
Inspection of the lists of proteins or transcripts that were differentially expressed in NOD relative to control strains revealed no overlap between them. While expression of a gene at the mRNA and protein level does not always correlate, other factors may be at play in our study, including: a) protein spots detected on 2D gels do not represent the total expression of the gene at the protein level but often just one of many post-translationally modified form of the protein; b) the protein analysis only included proteins within the pI range of 4-7; c) the percentage of the total proteome covered was much smaller than that the percentage of the transcriptome; and d) different statistical thresholds were applied to the two datasets (Benjamini-Hochberg, p,0.005 vs. no MTC, p,0.05). To this end, we analyzed gene expression at the RNA level of the 16 proteins identified (at 2 and 4 weeks combined) using the statistical cut-off level that was applied to the proteome data set. We found that 4 (Hsp90ab1, Hspa9, S100a5 and Serpinb1a) were also differentially expressed at the transcript level (Table 4). However, the majority were significantly differen-tially expressed only at the protein level. Thus, investigation at the proteome level allowed us to capture additional genes whose expression may be altered in NOD mice spleen leukocytes. Overall, the data indicate that the vast majority of genes (88%) that are abnormally expressed in NOD spleen leukocytes in comparison to control strains are repressed and that regulation of abnormal expression may occur at either the RNA or protein level.

Gene ontology (GO) enrichment analysis
To further understand the biological processes associated with abnormalities in the NOD mouse, we performed Gene ontology (GO) analysis of the NOD-specific strain differences using WebGestalt. Because GO analysis is better suited for larger gene lists, we used a slightly less stringent analysis for the transcriptome datasets to obtain expanded lists. Such bigger lists minimize the number of false negatives and increase the efficiency of GO classification. Thus, one-way ANOVA at cut-off p,0.05 with Benjamini-Hochberg followed by hierarchical clustering identified 187 (87% were of lower expression in NOD) and 327 (91% were of lower expression in NOD) NOD differentially expressed genes at 2 and 4 weeks, respectively. All the genes, including FC and adjusted p-values, are listed in Tables S4 (2 weeks) and S5 (4 weeks). Functional clustering of these genes revealed that the top ranked major biological process was metabolism ( Table 5). The major unique biological process terms at 2 weeks were biosynthesis, amine metabolism and RNA processing, while those at 4 weeks were regulation of biological processes, immune system process (further categorized into defense response, and antigen processing and presentation), and apoptosis. Under the molecular function GO (Table 6), the common (i.e. enriched in both the 2 and 4 week gene lists) top ranked term was transferase activity. Catalytic activity was uniquely enriched at 2 weeks, while binding (accounting for more than half of the genes) was uniquely enriched at 4 weeks. Kinase activity was also a significant category but only at 4 weeks. The most salient features under the cellular component category (Table 6) were significant enrichment for mitochondrion at 2 weeks, and nucleus and external side of plasma  membrane (including MHC protein complex) at 4 weeks. Gene Ontology analysis of the NOD differentially expressed proteins revealed involvement of the cytoskeleton, protein binding and response to stress. In summary, the GO analyses suggest that spleen leukocytes from NOD mice at the preinsulitis stage have impaired metabolism, which was evident at both 2 and 4 weeks, and impaired immune response and apoptosis, evident only at 4 weeks. The analysis also suggests abnormalities involving the mitochondria in 2 week old NOD mice.

KEGG pathway enrichment analysis
To identify well characterized molecular pathways that were significantly over-represented in the gene lists of NOD-specific strain differences, we performed KEGG pathway analysis. KEGG pathways are manually drawn maps representing well known molecular interaction and reaction networks. We conducted this analysis also using the expanded transcriptome gene lists described above under GO analysis. The results are shown in Table 7. Consistent with the GO analysis, the predominant category was metabolic pathways, which was highly significantly enriched at 2 weeks (p = 4.98e-8) and slightly less so at 4 weeks (p = 2.65e-5). The lists were also more significantly enriched for pathways involved in RNA splicing (spliceosome) at 2 weeks (p = 2.76e-6) than at 4 weeks (p = 0.0032). Both lists were also enriched for immune function pathways, including natural killer cell mediated cytotoxicity, antigen processing and presentation and T cell receptor signaling. Conversely to the metabolic pathways, the immune function pathways were more significantly enriched at 4 weeks than at 2 weeks. Additionally, each of the lists also had unique significantly enriched pathways, including: at 2 weeks, further metabolic pathways categorizations (such as purine/ pyrimidine metabolism, amino sugar and nucleotide sugar metabolism, and amino acid metabolism); and at 4 weeks, immune cell signaling pathways (such as graft-versus-host disease, B cell receptor signaling, cell adhesion molecules, etc) and autoimmune diseases (Type 1 diabetes mellitus and autoimmune thyroid disease). KEGG pathway analysis of the combined list of differentially expressed proteins revealed enrichment for 4 pathways (shown in Table 8) also pointing to the same themes as the transcriptome data.
Overall, the KEGG pathway analyses suggest that two broad molecular pathways -metabolic and immune response pathwaysare altered in the various immune cell types (natural killer cells, antigen presenting cells (B cells, macrophages, dendritic cells) and T cells (CD4 and CD8 T-cells). Alterations in the metabolic pathways appear to be predominant at 2 weeks, while those in the immune response and signaling pathways appear to be predominant at 4 weeks. Thus, the data suggest that impairment of the metabolic pathways may be the primary defect contributing to initiation of autoimmune pathology in NOD mice.

Ingenuity pathway analysis
To identify the genes and/or endogenous chemicals that may play a role in regulating the abnormal molecular pathways in NOD mice, we performed Ingenuity pathway analysis (IPA). IPA generates de novo networks from the uploaded gene lists based on information from the literature contained in the IPA knowledge  Table 3. Transcripts significantly differentially expressed in spleen leukocytes between 2 weeks and 4 weeks of age in NOD mice but not in two control strains.  Tables 1  and 2) generated several networks each. Each respective set of networks were merged to create one molecular network for each age (Figure 3 and Figure 4, for 2 and 4 weeks, respectively). The transcriptome network for the 2 week old mice dataset clustered around several central genes including MYC/MYCN, HNF4A, YWHAZ, TNF, IL2, IL15, TGFB1, progesterone, and Ca2+. That for the 4 week old mice dataset also clustered around several genes including TP53, YWHAZ, HNF4A, TNF, IFNG, IL15, beta-estradiol, progesterone, Akt, PRKCA, IL12, and HLA-C. IPA analyses of the lists of NOD differentially expressed proteins at 2 weeks (9 proteins) and at 4 weeks (7 proteins) generated one network each. The proteome network for the 2 week old mice centered around TP53, MYC, APP, IL1B, Ca2+ and HSPA8 ( Figure 5). That for the 4 week-old NOD mice clustered around TNF, CASP3, HSPD1, D-glucose and Jnk ( Figure 6). The major biological functions significantly represented by the transcriptome networks included cancer (indicative of apoptosis/cellular proliferation) and immunological disease/response ( Table 9). The top IPA functional categories represented by the proteome networks were cell proliferation, cell death and oxidative stress (Table 9). To facilitate an objective comparison between networks and to enable us gain further insights into the genes that may play a central role in each network, we ranked the central genes according to the total number of connections linked to each gene (Table 10). This analysis revealed that five central genes (TNF, HNF4A, IL15, Progesterone, and YWHAZ) were common to the transcriptome networks. A greater number of signaling and/or immune response genes (such as IFNG, EGFR, NFKB, Akt, PKRCA, IL12, and HLA-C) were present in the transcriptome network for 4 week-old mice than in that of 2 week-old mice. Inspection of the central genes between the proteome networks indicated that none were shared between the two networks (Table 10). However, a few molecules (TP53, MYC and Ca2+) were shared between the transcriptome and proteome networks of the 2 week old mice. Moreover, TNF was central to three networks, both the transcriptome networks and the proteome network for the 4 week-old mice. Taken together, IPA analyses of networks generated from NOD-specific strain differences identified several key factors that may play a role in regulating the molecular pathways that are abnormal in spleen leukocytes from NOD mice at the preinsulitis stage of autoimmune diabetes.
Because the gene list of NOD-specific age differences ( Table 3) was suggestive of immune activation, we conducted IPA on the three lists of strain-specific age differences (Tables 3, S2 and S3). The functional categories associated with the three topmost major networks generated from each of these lists are shown in Table 11. Only the NOD networks are associated with immunological disease and inflammatory response. The three networks for each strain were merged to create one molecular network for each strain ( Figures S4, S5 and S6 for NOD, NON and C57BL/6, respectively). The network for NOD age differences centered around several cytokines/cytokine receptors (IFNG, IL12, interferon alpha, IL17RB) and acute phase signaling molecules (PI3K, Akt, ERK, p38MAPK, NFKB). Interestingly, the network of NON-specific age differences also had IFNG as a central gene (but lacked the other cytokines observed in the NOD network) and only NFKB of the acute phase signaling molecules seen in the NOD network. However, the NON network also had TGFB, STAT1 and Ca2+ as central molecules. Strikingly, the network generated from the C57BL/6-specific age differences conspicuously lacked the cytokines seen in the NOD network (except interferon alpha) but had all the acute phase signaling molecules seen in the NOD network. In summary, in NOD mice, ,60% of the differentially expressed genes are upregulated at 4 weeks as compared to at 2 weeks, and a vast majority of these genes appear to be involved in the acute phase signaling cascade (PI3K, Akt, ERK, p38MAPK, NFKB, Figure S4), possibly driven by IFNG and other cytokines. In contrast, in C57BL/6 mice, ,75% of the differentially Table 4. Proteins with significant expression differences in spleen leukocytes of 2 week and 4 week old NOD mice compared to two control strains and their expression differences at RNA level. Statistical analysis of protein expression data by 1-way-ANOVA (p,0.05) followed by hierarchical clustering identified 11 and 8 spots that were differentially expressed in spleen leukocytes of 2 week and 4 week old NOD mice, respectively, compared to two control strains (NON and C57BL/6); H, indicates spots that were of higher expression in NOD, all other spots were of lower expression in NOD mice. We identified proteins in 9 and 7 of the spots at 2 and 4 weeks, respectively. We also analyzed expression differences of the identified proteins at the RNA level using the microarray (mRNA) dataset; NS, indicates no significance in gene expression between NOD and control strains. Fold change (FC) was calculated by ratio of means of expression in NOD mice versus controls; FC of spots/genes of lower expression in NOD are indicated by negative values. doi:10.1371/journal.pone.0046941.t004 expressed genes are downregulated at 4 weeks compared to at 2 weeks. Intriguingly, a large majority of these genes are also involved in the acute phase signaling cascade ( Figure S6). These results would suggest that the acute phase signaling pathway is activated in the NOD mice (inasmuch as the majority of the genes are upregulated) whereas it is suppressed in the C57BL/6 mice (because the majority of the genes are downregulated at 4 weeks). Thus, the data suggest immune activation in the NOD mice between 2 weeks and 4 weeks possibly driven by the inflammatory cytokine IFNG.

Discussion
Insulitis, the earliest sign of autoimmune pathology in the NOD mouse model for autoimmune (Type 1) diabetes, is first obvious at around 4-5 weeks of age [1]. The purpose of our study was to gain new insights into the molecular events that lead to this process and that may influence initiation of autoimmune diabetes. Since autoimmune diabetes is a multigenic disease, the present study tested the hypothesis that alterations in multiple molecular pathways lead to initiation of autoimmune pathology. We conducted a detailed analysis of genome-wide gene expression abnormalities in the spleen leukocytes (i.e. the effector tissue) of young NOD mice in the preinsulitis stage of autoimmune diabetes at ages 2 and 4 weeks. This study differs from our previous efforts [14] in that gene expression in the NOD mouse was compared to two controls (NON and C57BL/6) instead of one control (C57BL/ 6). This allowed us to eliminate C57BL/6-specific differences (i.e. C57 low and C57 high, Figure 1 and 2) and to focus on only those genes where NOD expression differed from that of C57BL/6 and an additional control (NON); thus, increasing the stringency of our attempt at identifying diabetes-related abnormalities. Additionally, we conducted a 1-w ANOVA in the current study instead of a 2-w ANOVA that we performed in the previous study [14]. This allowed us to put a stronger focus on differences unique to each specific age, a significant number of which could have been Table 5. Enriched biological process Gene ontology (GO) categories for the lists of differentially expressed transcripts in spleen leukocytes from 2 and 4 week old NOD mice compared to two control strains. One-way ANOVA (p,0.05, Benjamini-Hochberg) of mRNA expression data followed by hierarchical clustering identified 187 and 327 probe sets whose expression was altered in spleen leukocytes from NOD mice compared to two control strains (NON and C57BL/6) at 2 weeks and 4 weeks, respectively. excluded by a 2-w ANOVA. Furthermore, we applied a modern suite of bioinformatics approaches to gain the fullest possible insights from these datasets. The current study confirms the results of our previous study, which showed a defect in apoptosis/cell proliferation. Furthermore, it provides new insights into the molecular events that occur prior to insulitis in a time course fashion. Thus, we identified abnormalities in two broad pathways: metabolic, which were predominant at 2 weeks, and immune response pathways, which were predominant at 4 weeks. Surprisingly, analysis of the global gene expression profiles and hierarchical clustering revealed that the vast majority (,90%) of the abnormally expressed genes in NOD spleen leukocytes compared to two control strains were repressed. A previous study [23] also investigating gene expression changes in various tissues (pancreatic lymph nodes, spleen and peripheral blood cells) and at various ages (10 days, 4 weeks, 12 weeks, 16 weeks, or 20 weeks) in NOD mice compared to NOD.B10 mice (in which a nonpermissive MHC haplotype is imposed onto NOD to silence autoimmune disease) also reported that the great majority of identified NOD differentially expressed genes were down regulated (except only with peripheral blood cells at 4 weeks). It is important to note that, unlike our study, this study did not identify any differentially expressed genes in the spleen tissue at ages 2 weeks and 4 weeks. This could be due to their use of the whole spleen tissue rather than the leukocyte fraction as we did in our study and/or their narrowing changes down to only the diabetes susceptibility genes in the MHC region (excluding differences influenced by all the other chromosomal regions conferring diabetes susceptibility). Indeed, the authors noted that accurate identification and characterization of genes involved in complex autoimmune diseases like T1D may also require analysis of expression in specific cells of the tissues.
Gene ontology (GO) analyses of the NOD differentially expressed genes identified in our study revealed alterations in several processes that were common to both ages (Tables 5 and 6), including abnormalities in metabolism and enzymatic activity. Abnormalities in genes involved in metabolism (including the mitochondrion, carbohydrate and amine metabolism) and biosynthesis were predominant at 2 weeks of age while abnormalities in the immune system process (including response to external stimulus, antigen processing and presentation) were predominant at 4 weeks of age. KEGG pathway analysis (Tables 7 and 8), likewise, identified that abnormalities in metabolic pathways were predominantly enriched at 2 weeks, including spliceosome, One-way ANOVA (p,0.05, Benjamini-Hochberg) of mRNA expression data followed by hierarchical clustering identified 187 and 327 probe sets whose expression was altered in spleen leukocytes from NOD mice compared to two control strains (NON and C57BL/6) at 2 weeks and 4 weeks, respectively.  One-way ANOVA (p,0.05, Benjamini-Hochberg) of mRNA expression data followed by hierarchical clustering identified 187 and 327 genes whose expression was altered in spleen leukocytes from NOD mice compared to two control strains (NON and C57BL/6) at 2 weeks and 4 weeks, respectively. The lists of differentially expressed genes were analyzed in WebGestalt (http://bioinfo.vanderbilt.edu/webgestalt) for enriched KEGG pathways using the hypergeometric test (p,0.01, Benjamini-Hochberg); dashes (-) indicate that the corresponding categories were not significantly enriched at the respective ages. Pathways represented by $2 genes are shown. Bold text indicates categories that were common to both age groups; ECM, extracellular membrane. doi:10.1371/journal.pone.0046941.t007 Table 8. Enriched KEGG Pathways for the lists of differentially expressed proteins in spleen leukocytes from 2 week and 4 week old NOD mice compared to two control strains. One-way ANOVA (p,0.05) of protein expression data followed by hierarchical clustering and protein identification revealed 9 and 7 proteins whose expression was altered in spleen leukocytes from NOD mice compared to two control strains (NON and C57BL/6) at 2 weeks and 4 weeks of age, respectively. The two lists of differentially expressed proteins were combined and protein identifiers (IDs) were converted to gene IDs. purine/pyrimidine metabolism, amino acid and sugar metabolism. Alterations in the immune system function and signaling pathways including natural killer cell mediated cytotoxicity, antigen processing and presentation, B and T cell receptor signaling, graft-versus-host disease, Fc gamma R-mediated phagocytosis, type 1 diabetes mellitus, etc. were predominantly enriched at 4 weeks. Spleen leukocytes are naturally in a resting state with low level of protein synthesis until they encounter an activation signal, which leads to rapid growth followed by proliferation. This immune activation results in a dramatic demand for energy and biosynthetic precursors that is met by upregulation of metabolism. Because metabolism was a significantly enriched process in the NOD differentially expressed genes, yet the vast majority of these genes were repressed, we hypothesize that metabolism is repressed in NOD mice as compared to the two control strains. Preautoimmune metabolic changes have been previously reported in female NOD mice that later progressed to autoimmune diabetes [24]. Dysregulation of lipid and amino acid metabolism was observed in these mice. Interestingly, dysregulation of lipid and amino acid metabolism has also been observed in the serum of children who later progressed to human type 1 diabetes in the period prior to seroconversion to islet autoantibody positivity [25].
Thus, NOD autoimmune diabetes and human T1D may have some similarities in this respect. Furthermore, it should be noted that the elimination of autoreactive lymphocytes is normally done by apoptosis, an active process that requires ATP provided by metabolism.
Significantly enriched immune response pathways that we identified are consistent with the idea of an activation of immune system activity. However, we think that this process is repressed in NOD mice as compared to the two control strains because the vast majority of NOD differentially expressed genes were repressed. As stated above, a previous study also investigating gene expression changes in NOD mice [23] found the same result. However, they did not pursue this observation any further in their report. It is interesting to note, however, that many of the genes identified in this study were also found within the previously identified diabetes susceptibility regions, the Idd chromosomal loci. We think that the idea of an apparent paradoxical immune suppression in NOD mice may help explain previous observations that have been made in these mice: i) that immunosuppression exacerbates autoimmune diabetes in NOD mice [26]; and conversely, ii) that immunostimulation circumvents diabetes in NOD mice [27]. Notwithstanding that NOD mice have important differences in their immune system compared to humans, this idea of paradoxical immune week-old mice. The merged network was generated from the list of transcripts differentially expressed in 2 week-old NOD mice compared to both control strains, NON and C57BL/6. The list was selected from a hierarchical cluster of 175 genes that had highly significant expression differences between strains at 2 weeks. The genes derived from our uploaded gene list (focus genes) are represented by gray icons. White icons represent genes (or endogenous chemicals) derived solely from the IPA knowledge base and that could be algorithmically connected to the focus genes. doi:10.1371/journal.pone.0046941.g003 suppression appears to have support from human type 1 diabetes studies [28,29]. Elo et al. [28] found a high-level of repression of immune response genes (including MHC Class I and Class II molecules, T-cell receptor signaling, actin filament system and NF-kB signaling) in children that had developed T1D-associated autoantibodies (prediabetes) who eventually progressed to clinical T1D. Similar to our study, they also found evidence of impairment of transcription, transport, metabolism and cell proliferation. Additionally, another study found a repression of all genes that were differentially expressed in whole untreated CD4 T-cells from new onset T1D patients compared to matched controls [29]. This repression affected key immune functions, such as immune signaling pathways, and cell cycle regulation. In functional studies, these authors found that CD4 T-cell surface expression levels of two key adhesion molecules, LFA-1 and P-selectin, were significantly lower in T1D patients compared to controls. Furthermore, this study demonstrated that a lower percentage of CD4 T-cells from patients entered into mitosis as compared to those from controls suggesting inhibition of CD4 T-cell proliferation in T1D patients. They concluded that CD4 T-cells in new onset T1D have a hyporesponsive immune cell phenotype.
Whereas a generalized repression of the immune system in autoimmune diabetes may seem counterintuitive in a disease caused by activation of the immune system against target self antigens, it should be remembered that tolerance to self antigens requires the activation of self-reactive lymphocytes and their elimination by apoptosis as a result of this activation. Furthermore, elimination of self-reactivity by lymphocytes in the periphery also requires activation of regulatory mechanisms (such as regulatory T cells). Our data are consistent with the hypothesis that a defect in the mitochondria (Table 6) of NOD spleen leukocytes coupled with impaired metabolism leads to a defect in antigen processing/ presentation and T cell stimulation (Tables 5,6,7,8), ultimately resulting into generalized immune suppression. Islet-specific autoimmunity likely occurs due to defective tolerance mechanisms in the NOD mice either failing to eliminate islet-reactive T cells by apoptosis or neutralize them by regulatory T-cells activity. To this end, it is worth noting that some of the identified NOD differentially expressed genes (Ly6c1, Trim 12, and Trim30) are known to be linked to a defect in T cell negative selection. Furthermore, the MHC class I-related glycoprotein H60 (H60a), which we identified as being highly expressed in NOD, but virtually unexpressed in control strains (Tables 1 and 2; Figure S1), is known to play a role in immune regulatory mechanisms. H60a can mediate strong suppressive effects on T cell proliferation [30]. It has also been shown to play a role in NK cell tolerance in experimental autoimmune encephalomyelitis [31]. Additionally, a recent study, using a mouse model of spontaneous autoimmune arthritis [32], provided evidence suggesting that efficient suppression of autoimmune diseases may require polyclonal specificities rather than single antigen-specific regulatory T cell activity. week-old mice. The merged network was generated from the list of transcripts differentially expressed in 4 week-old NOD mice compared to both control strains, NON and C57BL/6. The list was selected from a hierarchical cluster of 189 genes that had highly significant expression differences between strains at 4 weeks. The genes derived from our uploaded gene list (focus genes) are represented by gray icons. White icons represent genes (or endogenous chemicals) derived solely from the IPA knowledge base and that could be algorithmically connected to the focus genes. doi:10.1371/journal.pone.0046941.g004 Network analysis by IPA identified several genes/molecules (central genes, Table 10) that we think may play a key role in regulating the identified abnormal processes. Thus, five central genes common to both transcriptome networks were identified, including the proinflammatory cytokine, TNF, a transcription factor (HNF4A), a cytokine (IL15) and two cellular activation/ signaling molecules (YWHAZ and progesterone). In addition to the common central genes, unique central genes were identified in the respective networks, e.g. cytokines IL2 and TGFB1, and transcription factors MYC/MYCN in 2 week transcriptome network; and cytokines IFNG and IL12, transcription factors/ regulators beta-estradiol, TP53, NFKB, and AR, and cellular activation/signaling EGFR, Akt, PRKCA in 4 week transcriptome network. Although, several of these genes have been previously known to be associated with T1D development, and with immune function, many are novel.
TNF-a, a proinflammatory, multifunctional cytokine was present in three of the four networks underscoring its importance in the initiation of autoimmune diabetes. TNF-a in synergy with IFN-c (which was present in the 4 week transcriptome network) is known to be important in inducing autoimmunity in the developing immune system [33][34][35]. Also, another report indicated that progression to insulitis in the NOD mouse was connected to changes in T cell activation and signaling associated with TNFa and another proinflammatory cytokine, IL-6 [36]. Our identification of IL2 as a central gene (at 2 weeks) further supports involvement of IL2 signaling in T1D pathogenesis. The defect in IL2 may contribute to the impairment of regulatory T-cells [37]. The NFKB pathways are essential for the function of antigen presenting cells and related defects have previously been reported in T1D patients [38]. The pregnancy hormones, beta-estradiol and progesterone are known to lead to improvement in T1D in some patients as evidenced by a decline in insulin requirement during pregnancy [39,40]. Estrogens are known to play a role in modulating immune function of various immune cells and may also play a role in regulating glucose homeostasis and plasma lipids [41]. Progesterone is known to modulate macrophage activation [42].
Myc, an important factor in regulating apoptosis and cellular growth, was central in both the transcriptome and proteome based networks of the 2 week old mice. It has been hypothesized that defects in the process of apoptosis might explain why autoreactive leukocytes persist in the NOD mouse [43][44][45]. Indeed, upregulation of Myc in lymphocytes from NOD mice (in contrast to Myc downregulation in C57BL/6) was associated with their resistance to glucocorticoid-induced apoptosis [46]. Additionally, NOD mice treated with the streptococcal wall component OK432, which restores dexamethasone-induced apoptosis and is associated with downregulation of Myc, do not develop diabetes. Recently, Wang et al. [47] revealed a novel function of Myc as an essential regulator of activation-induced T cell metabolic reprogramming. T cell activation (and other immune cells) leads to a metabolic Figure 5. Proteome network created by Ingenuity Pathway Analysis from the dataset of 2 week-old mice. A single network was generated from the list of proteins that were differentially expressed in 2 week-old NOD mice compared to both control strains, NON and C57BL/6. The protein spots from which these proteins were identified were selected from the hierarchical cluster of 19 protein spots that had significant expression differences between strains at 2 weeks. The gene names derived from our uploaded lists (focus genes) are represented by gray icons. White icons represent genes (or endogenous chemicals) derived solely from the IPA knowledge base and that could be algorithmically connected to the focus genes. doi:10.1371/journal.pone.0046941.g005 switch to increased glucose and glutamine metabolism. These authors showed that Myc plays a key role in upregulation of glucose, glutamine, and polyamine metabolism and that its deficiency prevented T-cell proliferation. In addition to Myc as a metabolic regulator, the Akt-mTOR (mammalian target of rapamycin) pathway also promotes posttranslational regulation of glycolytic metabolism [48]. Akt was central in the transcriptome network at 4 weeks. Table 9. Functional categories associated with the major Ingenuity Pathway Networks of differentially expressed transcripts/ proteins in spleen leukocytes from 2 week and 4 week old NOD mice.  week-old mice. A single network was generated from the list of proteins that were differentially expressed in 4 week-old NOD mice compared to both control strains, NON and C57BL/6. The protein spots from which these proteins were identified were selected from the hierarchical cluster of 26 protein spots that had significant expression differences between strains at 4 weeks. The gene names derived from our uploaded lists (focus genes) are represented by gray icons. White icons represent genes (or endogenous chemicals) derived solely from the IPA knowledge base and that could be algorithmically connected to the focus genes. doi:10.1371/journal.pone.0046941.g006 In addition to the above components, our analysis identified components previously unrecognized to play a role in immunological function, including, p53, HNF4A, and AR (androgen receptor). Tumor suppressor p53 is widely known to act in response to excessive stresses and abnormalities in cell physiology. Recent studies, however, suggest that p53 has additional important roles in normal cellular physiology [49]. To this end, the p53 pathway is tightly involved in the homeostatic regulation of metabolic processes. In these novel functions, p53 helps to align metabolic processes with the proliferation and energy status of the cells, to maintain optimal mode of glucose metabolism and to boost the energy efficient mitochondrial respiration in response to ATP deficiency. Additionally, under normal physiological conditions, p53 mobilizes defenses against ''physiological'' levels of reactive oxygen species (ROS) to maintain homeostasis. Thus, a possible hypothesis would be that the novel functions of p53 play a role in the physiological immunological process of immune cell activation in response to stimuli, including antigenic stimulation. Our data suggest that these p53 functions are abnormal in NOD mice. HNF4A, a nuclear receptor transcription factor, has been linked to developmental and metabolic functions, and to several diseases, including maturity-onset diabetes of the young and type 2 diabetes [50]. Although HNF4A is known to be expressed in the spleen (among several other tissues), its role in the immunological process has yet to be fully established [51]. Understanding this role may be important as HNF4A could be a good candidate for a drug target. Indeed, this idea has been further reinforced by the recent identification of linoleic acid as the long-sought for ligand for HNF4A [52]. HNF4A controls the expression of numerous genes, including well known target genes involved in regulating cellular metabolism and homeostasis, as well as many new direct target genes involved in regulation of the cell cycle and immune function [53,54]. AR is a nuclear receptor transcription factor involved in the regulation of many different physiological processes, including developmental, immunological, metabolic and homeostatic processes [55][56][57][58][59]. Classical AR signaling occurs via binding of androgens to the AR, which then translocates to the nucleus and activates transcription leading to growth promotion and differentiation. However, AR can also be activated in an androgenindependent fashion via growth factors or cytokines. AR dysfunction is associated with a wide variety of pathologies including prostate cancer, type 1 diabetes and metabolic disorders, such as type 2 diabetes.
In summary, our study identified abnormal molecular pathways in the spleen leukocytes from young NOD mice in the preinsulitis (i.e. early inductive) stage of autoimmune diabetes. The pathways associated with the NOD-specific strain differences can be divided broadly into 2 categories, metabolic pathways (predominant at 2 weeks) and immune response pathways (predominant at 4 weeks). The identified putative key regulators of the metabolic pathways include Myc and HNF4A at 2 weeks, and beta-estradiol, p53, Akt, HNF4A and AR at 4 weeks. The molecular pathways associated with NOD-specific expression differences between 2 week old and 4 week old mice suggest an immune activation phenotype in these mice. We hypothesize that the NOD-specific strain differences we identified represent a defect in immune tolerance predisposing NOD mice to autoimmune diabetes development while the age differences represent the additional changes that accompany initiation of insulitis. It is worth noting that since spleen leukocytes constitute a heterogeneous cell population the molecular networks identified by IPA are likely representative of the predominant cell types. Thus, future experiments to dissect altered molecular networks within individual cell types are warranted. Indeed such studies may reveal age-specific gene expression abnormalities that could have been masked in the current study. Notwithstanding, our data on untreated unfractionated spleen leukocytes suggest that that abnormalities in regulation of metabolic pathways in the immune cells of young NOD mice lead to abnormalities in the immune response pathways and as such may play a role in the initiation of autoimmune diabetes. Indeed our study suggests that targeting metabolism may be a novel strategy to modulate autoimmunity and to prevent and/or treat autoimmune diabetes. Table 11. Functional categories associated with the 3 topmost major IPA Networks of transcripts uniquely differentially expressed in spleen leukocytes between 2 weeks and 4 weeks of age in NOD, NON and C57BL/6 mice.   Figure S4 Transcriptome network created by Ingenuity Pathway Analysis from genes uniquely differentially expressed in NOD mice between 2 weeks and 4 weeks of age. The merged network was generated from the list of transcripts uniquely differentially expressed in spleen leukocytes of NOD mice between 2 weeks and 4 weeks of age in comparison to two control strains, NON and C57BL/6. It represents the three topmost major networks. The genes derived from our uploaded gene list (focus genes) are represented by gray icons. White icons represent genes (or endogenous chemicals) derived solely from the IPA knowledge base and that could be algorithmically connected to the focus genes.
(TIFF) Figure S5 Transcriptome network created by Ingenuity Pathway Analysis from genes uniquely differentially expressed in NON mice between 2 weeks and 4 weeks of age. The merged network was generated from the list of transcripts uniquely differentially expressed in spleen leukocytes of NON mice between 2 weeks and 4 weeks of age in comparison to NOD and C57BL/6. It represents the three topmost major networks. The genes derived from our uploaded gene list (focus genes) are represented by gray icons. White icons represent genes (or endogenous chemicals) derived solely from the IPA knowledge base and that could be algorithmically connected to the focus genes.
(TIFF) Figure S6 Transcriptome network created by Ingenuity Pathway Analysis from genes uniquely differentially expressed in C57BL/6 mice between 2 weeks and 4 weeks of age. The merged network was generated from the list of transcripts uniquely differentially expressed in spleen leukocytes of C57BL/6 mice between 2 weeks and 4 weeks of age in comparison to NOD and NON mice. It represents the three topmost major networks. The genes derived from our uploaded gene list (focus genes) are represented by gray icons. White icons represent genes (or endogenous chemicals) derived solely from the IPA knowledge base and that could be algorithmically connected to the focus genes. (TIFF) Table S1 Genes and primers/probes used for quantitative Real-time PCR. (DOC)