Abstract

   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 [33][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 [34][2], [35][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 [36][2], [37][4], [38][5]. It is clear that
   MHC class II-restricted CD4 T-cells play an important role in inducing
   T1D [39][6]–[40][9]. In addition to the MHC locus, which is required
   for T1D development, over 20 non-MHC loci also contribute to the
   disease process [41][5], [42][10]–[43][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.

Materials and Methods

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 [44][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 [45][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 ([46]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: [47]GSE37450). The microarray results were validated
   by quantitative Real-time PCR (qRT-PCR) using Roche's Universal library
   probe system ([48]www.roche-applied-science.com). First strand cDNA was
   synthesized from 1 µg 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 [49]Table S1. Expression
   values were normalized using glyceraldehyde-3-phosphate dehydrogenase
   (Gapdh).

Proteome data collection

   Protein expression was analyzed by 2D PAGE as decribed previously
   [50][14]. The protein gels were scanned and analyzed as previously
   described [51][15], [52][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 [53][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 [54][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) [55][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 entrez-gene ID using the gene ID conversion tool
   available at the DAVID (Database for Annotation, Visualization and
   Integrated Discovery) server (david.abcc.ncifcrf.gov) [56][18],
   [57][19]. To identify molecular networks that are associated with the
   differentially expressed genes, we analyzed the gene lists using
   Ingenuity Pathway Analysis (IPA, [58]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 [59]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 week-old 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
   [60][20]. The hierarchical clusters of the 2 week ([61]Figure 1) and 4
   week ([62]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
   [63]Figure 1, and listed in [64]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 [65]Figure 2 and listed in
   [66]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 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 [67]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 [68][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 [69][22].

Figure 1. Hierarchical clustering of 175 genes with highly significant
expression differences between strains at 2 weeks.

   [70]Figure 1
   [71]Open in a new tab

   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).

Figure 2. Hierarchical clustering of 189 genes with highly significant
expression differences between strains at 4 weeks.

   [72]Figure 2
   [73]Open in a new tab

   The list was generated by one-way ANOVA of 24,959 probe sets (p<0.005,
   Benjamini-Hochberg multiple test correction). A total of 76 genes were
   differentially expressed in NOD mice relative to both control strains
   (NON and C57BL/6): 67 were of lower expression (L) and 9 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).

Table 1. Transcripts with highly significant expression differences in spleen
leukocytes of 2 week old NOD mice compared to two control strains.

   Probe Set ID Gene Symbol Fold Change Adjusted p-value Description
   Immune response
   1457769_at (H) H60A 29.75 2.41e-5 histocompatibility 60a
   1443858_at Trim5 −20.14 0.000656 tripartite motif protein 5
   1442338_at Tnfrsf14 −8.68 0.00367 tumor necrosis factor receptor
   superfamily, member 14
   1450355_a_at Capg −5.9 8.87e-5 capping protein (actin filament),
   gelsolin-like
   1421571_a_at Ly6c1 −5.15 0.00029 lymphocyte antigen 6 complex, locus C1
   1437937_at (H) Ccbp2 2.77 0.000834 Chemokine binding protein 2
   1449991_at (H) Cd244 2.74 0.00099 CD244 molecule, natural killer cell
   receptor 2B4
   Apoptosis
   1423747_a_at Pdk1 −5.14 6.47e-5 pyruvate dehydrogenase kinase,
   isoenzyme 1
   1438454_at Pten −4.22 0.000953 Phosphatase and tensin homolog
   Ubiquitination
   1437432_a_at Trim12 −26.06 3.39e-5 tripartite motif protein 12
   1437007_x_at Usp39 −7.80 6.47e-5 ubiquitin specific protease 39
   1446608_at Cbl −5.98 0.00487 Casitas B-lineage lymphoma
   1428592_s_at Usp38 −2.38 0.000396 ubiquitin specific protease 38
   Transporters
   1451602_at Snx6 −17.31 0.000834 sorting nexin 6
   1417600_at Slc15a2 −6.88 0.00179 solute carrier family 15 (H+/peptide
   transporter), member 2
   1416629_at Slc1a5 −4.53 0.00428 solute carrier family 1, member 5
   1418473_at Cutc −1.79 0.00391 cutC copper transporter homolog (E.coli)
   Transcription
   1439780_at Rpl7l1 −6.12 0.000895 ribosomal protein L7-like 1
   1449648_s_at Rpo1-1 −3.86 0.00109 RNA polymerase 1-1 (Ingenuity symbol:
   POLR1C)
   1424170_at Phf5a −2.16 0.000608 PHD finger protein 5A
   Protein biosynthesis
   1454714_x_at Phgdh −10.28 3.11e-5 phosphoglycerate dehydrogenase
   1448472_at Vars −4.14 0.000382 valyl-tRNA synthetase 2
   1417762_a_at Rpl8 −2.13 0.00348 ribosomal protein L8
   1455887_at Alg8 −1.17 0.000555 Asparagine-linked gylcosylation 8
   homolog
   Phosphorylation
   1446130_at Pctk2 −7.98 0.00453 PCTAIRE-motif protein kinase 2
   1416770_at Stk25 −2.98 0.000365 serine/threonine kinase 25 (yeast)
   1448450_at Ak2 −1.60 0.00206 adenylate kinase 2
   Enzymes
   1421802_at Ear1 −103.15 0.00138 eosinophil-associated, ribonuclease A
   family, member 1
   1416494_at Ndufs5 −6.89 0.000102 NADH dehydrogenase (ubiquinone) Fe-S
   protein 5
   1427441_a_at Suclg2 −2.64 0.000354 succinate-Coenzyme A ligase,
   GDP-forming, beta subunit
   1451277_at Zadh2 −1.50 0.00236 zinc binding alcohol dehydrogenase,
   domain containing 2
   Other
   1417461_at (H) Cap1 25.17 3.96e-5 CAP, adenylate cyclase-associated
   protein 1 (yeast)
   1456635_at Sp110 −9.89 3.07e-5 Sp110 nuclear body protein
   1416479_a_at Tmem14c −9.30 2.41e-5 transmembrane protein 14C
   1435553_at Pdzd2 −9.25 0.000938 PDZ domain containing 2
   1416503_at Lxn −3.32 0.00231 Latexin
   1434441_at C9orf21 −3.23 7.15e-5 chromosome 9 open reading frame
   21(peroxiredoxin (PRX)-like 2 family)
   1436061_at Chaf1a −3.03 0.00139 chromatin assembly factor 1, subunit A
   (p150)
   1435745_at 5031439G07Rik −2.94 0.00164 chromosome 22 open reading frame
   9 (C22ORF9)
   1435660_at (H) LOC434484 2.79 0.000271 similar to SP140 nuclear body
   protein; nuclear body protein Sp140
   1418072_at Hist1h2bi −2.77 0.000993 histone cluster 1, H2bi
   1459827_x_at Hps1 −2.40 0.00083 Hermansky-Pudlak syndrome 1 homolog
   (human)
   1436836_x_at Cnn3 −2.00 7.15e-5 calponin 3, acidic
   Unknown
   1426936_at LOC215866 −15.28 0.000834 Mus musculus, cDNA clone
   IMAGE:3493956, mRNA
   1417822_at D17H6S56E-5 −7.76 6.43e-7 DNA segment, Chr 17, human D6S56E
   5
   1440214_at A630001G21Rik −4.04 0.00083 RIKEN cDNA A630001G21 gene
   [74]Open in a new tab

   Statistical analysis of mRNA expression data by 1-way-ANOVA
   (Benjamini-Hochberg, p<0.005) followed by hierarchical clustering
   identified 51 probe sets (representing 47 different genes) that were
   differentially expressed in spleen leukocytes of 2 week old NOD mice
   compared to two control strains (NON and C57BL/6); H, indicates 5 genes
   (∼11%) with significantly higher expression in NOD relative to
   controls; the rest, 42 (∼89%) had significantly lower expression in
   NOD. Fold change (FC) was calculated by ratio of means of expression in
   NOD mice versus controls; FC of genes of lower expression in NOD are
   indicated by negative values. We assigned the genes to selected
   functional categories based on information obtained by searches of
   public databases. Genes highlighted with bold font were differentially
   expressed at both 2 and 4 weeks.

Table 2. Transcripts with highly significant expression differences in spleen
leukocytes of 4 week old NOD mice compared to two control strains.

   Probe Set ID Gene Symbol Fold Change Adjusted p-value Description
   Immune response
   1457769_at (H) H60A 26.55 7.49e-5 histocompatibility 60a
   1443858_at Trim5 −17.29 0.00179 tripartite motif protein 5
   1442338_at Tnfrsf14 −6.78 0.000529 tumor necrosis factor receptor
   superfamily, member 14
   1421571_a_at Ly6c1 −4.61 1.76e-5 lymphocyte antigen 6 complex, locus C1
   1427799_x_at (H) Igk 4.46 0.000109 immunoglobulin kappa chain variable
   21 (V21)-4
   1450355_a_at Capg −3.95 3.64e-6 capping protein (actin filament),
   gelsolin-like
   1451965_at Igk −3.65 0.000255 Mouse mRNA for immunoglobulin kappa
   variable region
   1424948_x_at (H) H2-D1 2.89 0.00039 major histocompatibility complex,
   class I, C
   1422188_s_at Tcrg-V4 −2.81 0.00156 T-cell receptor gamma, variable 4
   1426112_a_at Cd72 −2.55 0.00206 CD72 antigen
   1425477_x_at (H) H2-Ab1 2.40 0.00402 major histocompatibility complex,
   class II, DQ beta 2
   1425436_x_at Klra7 −2.30 0.000858 killer cell lectin-like receptor,
   subfamily A, member 7
   1448706_at Ttrap −2.08 0.00344 Traf and Tnf receptor associated protein
   1442219_at Ms4a4b −2.07 0.00206 membrane-spanning 4-domains, subfamily
   A, member 4B
   1456426_at Clec2i −1.93 0.0029 C-type lectin domain family 2, member i
   1460204_at Tec −1.89 0.00292 cytoplasmic tyrosine kinase, Dscr28C
   related (Drosophila)
   1419481_at Sell −1.88 0.00402 selectin, lymphocyte
   1460431_at Gcnt1 −1.80 0.00234 glucosaminyl (N-acetyl) transferase 1,
   core 2
   1448862_at Icam2 −1.66 0.0031 intercellular adhesion molecule 2
   Apoptosis
   1431708_a_at Tia1 −12.96 2.25e-5 TIA1 cytotoxic granule-associated RNA
   binding protein 1
   1423747_a_at Pdk1 −6.65 1.99e-5 pyruvate dehydrogenase kinase,
   isoenzyme 1
   1455287_at Cdk6 −2.74 0.00276 Cyclin-dependent kinase 6
   1437393_at Prkca −2.28 0.00282 protein kinase C, alpha
   1434312_at Arf6 −1.98 0.00342 ADP-ribosylation factor 6
   1419706_a_at Akap12 −1.88 0.00319 A kinase (PRKA) anchor protein
   (gravin) 12
   Ubiquitination
   1437432_a_at Trim12 −51.45 6.17e-6 tripartite motif protein 12
   1437007_x_at Usp39 −5.28 0.00123 ubiquitin specific protease 39
   1424857_a_at Trim34 −2.64 6.22e-5 tripartite motif protein 34
   Transporters
   1451602_at Snx6 −8.74 0.00377 sorting nexin 6
   1416629_at Slc1a5 −2.97 0.000477 solute carrier family 1, member 5
   1426775_s_at (H) Scamp1 2.33 0.00269 secretory carrier membrane protein
   1
   1433898_at Slc25a30 −1.93 0.00261 solute carrier family 25, member 30
   1460224_at Snx2 −1.63 0.00427 sorting nexin 2
   Transcription
   1440845_at Dmtf1 −12.13 3.64e-6 cyclin D binding myb-like transcription
   factor 1
   1441823_at Ria17 −7.04 0.000744 retinoic acid induced 17/Zinc finger,
   MIZ-type containing 1
   1449648_s_at Rpo1-1 −3.19 0.00241 RNA polymerase 1-1
   1417961_a_at Trim30 −2.99 0.00123 tripartite motif protein 30
   1428662_a_at Hopx −2.24 0.00354 homeobox only domain, MGI:1916782
   1433599_at Baz1a −1.96 0.00343 bromodomain adjacent to zinc finger
   domain 1A
   Protein biosynthesis
   1454714_x_at Phgdh −5.33 3.14e-5 phosphoglycerate dehydrogenase
   1417762_a_at Rpl8 −2.01 0.00015 ribosomal protein L8
   1420994_at B3gnt5 −1.89 0.00292 UDP-GlcNAc:betaGal
   beta-1,3-N-acetylglucosaminyltransferase 5
   Phosphorylation
   1451253_at Pxk −2.99 0.00142 PX domain containing serine/threonine
   kinase
   1416770_at Stk25 −2.35 0.00214 serine/threonine kinase 25 (yeast)
   Enzymes
   1416494_at Ndufs5 −7.24 1.33e-5 NADH dehydrogenase (ubiquinone) Fe-S
   protein 5
   1448942_at (H) Gng11 3.64 0.00247 guanine nucleotide binding protein (G
   protein), gamma 11
   1447458_at St3gal4 −2.92 0.00214 ST3 beta-galactoside
   alpha-2,3-sialyltransferase 4
   1453109_at Arsk −2.87 0.00113 Arylsulfatase K
   1455734_at Crbn −2.41 0.00214 cereblon
   1436494_x_at Trmt1 −2.23 0.00478 Trm1 tRNA methyltransferase 1homolog
   (yeast)
   1417316_at Them2 −1.81 0.00214 thioesterase superfamily member 2
   1424043_at Ppil4 −1.58 0.00223 peptidylprolyl isomerase
   (cyclophilin)-like 4
   Other – partial list; other genes are listed below in the legend.
   1417461_at (H) Cap1 29.41 0.00017 CAP, adenylate cyclase-associated
   protein 1 (yeast)
   1456635_at Sp110 −13.82 2.10e-5 Sp110 nuclear body protein
   1416479_a_at Tmem14c −7.51 6.15e-5 Transmembrane protein 14C
   1436061_at Chaf1a −7.05 0.00185 chromatin assembly factor 1, subunit A
   (p150)
   [75]Open in a new tab

   Statistical analysis of mRNA expression data by 1-way-ANOVA
   (Benjamini-Hochberg, p<0.005) followed by hierarchical clustering
   identified 76 probe sets (representing 72 different genes) that were
   differentially expressed in spleen leukocytes of 4 week old NOD mice
   compared to two control strains (NON and C57BL/6); H, indicates 8 genes
   (∼11%) with significantly higher expression in NOD relative to
   controls; the rest, 68 (∼89%) had significantly lower expression in
   NOD. Fold change (FC) was calculated by ratio of means of expression in
   NOD mice versus controls; FC of genes of lower expression in NOD are
   indicated by negative values. We assigned the genes to selected
   functional categories based on information obtained by searches of
   public databases. Other genes under “Other” include: Fcrla, Otu1,
   Ankrd16, Treml2, Plekha2, Fam3c, 0610039P13Rik, Fam65b, Pyhin1,
   5031439G07Rik, C9orf21, LOC434484 (H); Unknown genes include:
   D17H6S56E-5, A630001G21Rik, A530030E21Rik. Genes highlighted with bold
   font were differentially expressed at both 2 and 4 weeks.

   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, [76]Table 3), 122 probesets (representing 120 genes,
   [77]Table S2) and 105 probesets (representing 104 genes; [78]Table S3)
   that were uniquely differentially expressed in spleen leukocytes of
   NOD, NON and C57BL/6 mice, respectively. Expression of ∼62% [79][40] of
   the NOD-specific age differences was upregulated at 4 weeks compared to
   expression at 2 weeks of age ([80]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).

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.

   Probe Set ID Gene Symbol Fold Change Regulation at 4 wks p-value
   Description
   1459452_at Cml2 27.32 up 0.000959 camello-like 2
   1418086_at Ppp1r14a 15.68 up 0.001 protein phosphatase 1, regulatory
   (inhibitor) subunit 14A
   1424631_a_at Ighg 11.77 up 7.21e-5 Immunoglobulin heavy chain (gamma
   polypeptide)
   1435987_x_at 1110059G02Rik 8.3 up 0.000645 RIKEN cDNA 1110059G02 gene
   1440771_at Zkscan1 −7.59 down 0 zinc finger with KRAB and SCAN domains
   1
   1431640_at 4933431J24Rik 7.37 up 0.00056 RIKEN cDNA 4933431J24 gene
   1432928_at Plcd4 6.93 up 5.25e-6 phospholipase C, delta 4
   1425385_a_at Ighm 6.9 up 0.000779 immunoglobulin heavy constant mu
   1437986_x_at Fuk 6.63 up 9.22e-5 fucokinase
   1444140_at Pum1 −6.43 down 0.000796 pumilio 1 (PUF family RNA-binding
   protein)
   1452358_at Rai2 5.67 up 7.14e-5 retinoic acid induced 2
   1426174_s_at Ighg 5.45 up 0.000444 Immunoglobulin heavy chain (gamma
   polypeptide)
   1423243_at Mpp1 −4.79 down 0.000209 membrane protein, palmitoylated 1,
   55kDa
   1432782_at 4933435C09Rik 4.63 up 4.78e-5 RIKEN cDNA 4933435C09 gene
   1441524_at Tec 4.53 up 0.001085 cytoplasmic tyrosine kinase, Dscr28C
   related (Drosophila)
   1439855_at Tmtc1 −4.48 down 0.0003 transmembrane and tetratricopeptide
   repeat containing 1
   1426547_at Gc 4.21 up 1.41e-5 group-specific component (vitamin D
   binding protein)
   1435367_at Mapk4 −4.08 down 0.000806 mitogen-activated protein kinase 4
   1453957_a_at Igf2bp3 −4.03 down 0.00029 insulin-like growth factor 2
   mRNA binding protein 3
   1441710_at 4930479D17Rik −3.57 down 0.001082 RIKEN cDNA 4930479D17 gene
   1447934_at C12orf5 −3.55 down 0.000111 chromosome 12 open reading frame
   5
   1451492_at Sla2 3.53 up 7.55e-5 Src-like-adaptor 2
   1441955_s_at Paip1 3.30 up 0.000908 poly(A) binding protein interacting
   protein 1
   1444416_at CenpA −3.23 down 0.000564 centromere protein A
   1432162_s_at Ptar1 3.10 up 0.000509 protein prenyltransferase alpha
   subunit repeat containing 1
   1453760_at Mier1 −3.01 down 0.000385 mesoderm induction early response
   1 homolog (Xenopus laevis)
   1432530_s_at Boll 2.97 up 0.000906 bol, boule-like (Drosophila)
   1430717_at Fam83e 2.84 up 0.000592 family with sequence similarity 83,
   member E
   1449306_at Hsf2 2.67 up 0.000948 heat shock transcription factor 2
   1449306_at Rps4y2 −2.62 down 0.000292 ribosomal protein S4, Y-linked 2
   1434405_at Fnip1 −2.46 down 0.000305 folliculin interacting protein 1
   1459227_at 953001807Rik 2.40 up 0.000204 RIKEN cDNA 9530018I07 gene
   1429754_a_at 1700110I01Rik 2.40 up 0.000381 predicted gene 5957
   1426041_a_at Fgd4 −2.34 down 0.000524 FYVE, RhoGEF and PH domain
   containing 4
   1442315_at Foxd2 2.28 up 0.000931 forkhead box D2
   1446698_at Lpp −2.27 down 0.000523 LIM domain containing preferred
   translocation partner in lipoma
   1449153_at Mmp12 2.22 up 0.000549 matrix metallopeptidase 12
   (macrophage elastase)
   1437230_at Kcna1 −2.22 down 0.001029 potassium voltage-gated channel,
   shaker-related subfamily, member 1 (episodic ataxia with myokymia)
   1426218_at Glcci1 −2.19 down 0.000131 glucocorticoid induced transcript
   1
   1424698_s_at Gca −2.13 down 0.001086 grancalcin, EF-hand calcium
   binding protein
   1420731_a_at Csrp2 −2.05 down 9.09e-5 cysteine and glycine-rich protein
   2
   1427471_at Fbxl3 −2.02 down 0.000372 F-box and leucine-rich repeat
   protein 3
   1427718_a_at Mdm2 −1.97 down 0.001089 Mdm2, p53 E3 ubiquitin protein
   ligase homolog (mouse)
   1420330_at Clec4e 1.94 up 0.000866 C-type lectin domain family 4,
   member E
   1435856_x_at Smarcb1 −1.89 down 0.000589 SWI/SNF related, matrix
   associated, actin dependent regulator of chromatin, subfamily b, member
   1
   1416107_at Hmp19 1.83 up 0.001088 HMP19 protein, neuron specific gene
   family member 2
   1425294_at Slamf8 1.83 up 0.000142 SLAM family member 8
   1443884_at Thada 1.80 up 0.00066 thyroid adenoma associated
   1449025_at Ifit3 1.73 up 0.000955 interferon-induced protein with
   tetratricopeptide repeats 3
   1419634_a_at Ghrh 1.73 up 0.000456 growth hormone releasing hormone
   1450290_at Pdcd1lg2 1.67 up 0.000798 programmed cell death 1 ligand 2
   1426276_at Ifih1 1.64 up 0.000163 interferon induced with helicase C
   domain 1
   1454799_at AgpaT9 −1.64 down 0.000332 1-acylglycerol-3-phosphate
   O-acyltransferase 9
   1460651_at Lat 1.61 up 0.000971 linker for activation of T cells
   1416087_at Ap1s1 −1.50 down 0.000297 adaptor-related protein complex 1,
   sigma 1 subunit
   Genes with fold change of <1.5 and their regulation at 4 weeks are
   listed in the legend.
   [81]Open in a new tab

   Statistical analysis of mRNA expression data by unpaired t-Test
   (p<0.001) followed by Venn diagram analysis identified 66 probe sets
   (representing 65 different genes) that were uniquely differentially
   expressed in spleen leukocytes of NOD mice between 2 and 4 weeks of age
   in comparison to two control strains, NON and C57BL/6. Fold change (FC)
   was calculated by ratio of means of expression at 2 weeks versus 4
   weeks; FC of genes of lower expression at 4 weeks (compared to 2 weeks)
   are indicated by negative values. Genes with FC<1.5 (and their FC
   values) included: Syt11 (1.47), Snap47 (1.46), Cklf (−1.45), Cd274
   (1.43), Heatr7a (1.42), Fam76b (−1.41), Npsr1 (1.39), Nrn1 (1.38),
   Eef2k (1.24), Tomm34 (1.21), and IL17rb (1.2). Expression of
   approximately 62% (40) of the genes was upregulated at 4 weeks compared
   to expression at 2 weeks of age.

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 [82][14]. Supplemental
   [83]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) ([84]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.

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.

   Spot ID Fold Change (FC) Proteins identified p-value at RNA level (FC)
   Description
   2 week
   149 (H) 3.08 Clip1 NS Cap-Gly domain containing linker protein 1
   190 −2.39 Hsp90ab1 0.02 (−1.26) Heat shock protein1, beta
   532 −1.82 Lmnb1 and Lcp1 NS Lamin B1 and Lymphocyte cytosolic protein
   1, respectively
   441 −1.80 Tuba1 and Vim NS Tubulin alpha 1 and Vimentin, respectively
   24 (H) 1.78 Hspa9 0.03 (−1.22) Heat shock protein 9A
   423 −1.74 Lmnb1 and Lcp1 NS Lamin B1 and Lymphocyte cytosolic protein
   1, respectively
   471 −1.65 Hspa8 NS Heat shock protein 8
   601 −1.50 Hspa8 NS Heat shock protein 8
   618 (H) 1.33 Gopc NS Golgi associated PDZ and coiled-coil motif
   containing
   45 (H) 1.74 unidentified unknown not available
   702 (H) 1.36 unidentified unknown not available
   4 week
   74 (H) 6.61 Arhgdib NS Rho, GDP dissociation inhibitor (GDI) beta
   72 (H) 6.44 S100a5 0.03 (1.78) S100 calcium binding protein A5
   132 −3.45 Ctnna1 NS Catenin alpha 1
   295 (H) 2.11 Hspd1 NS Heat shock protein 1(chaperonin)
   430 −2.00 Ppp5c NS Protein phosphatase 5, catalytic subunit
   390 −1.91 Serpinb1a 0.006 (−1.37) Serine (or cysteine) proteinase
   inhibitor, clade B, member 1a
   225 −1.88 Wdr1 NS WD repeat domain 1
   485 −1.83 unidentified NS not available
   [85]Open in a new tab

   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.

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 ([86]Table 4). However, the majority were
   significantly differentially 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 [87]Tables S4 (2
   weeks) and [88]S5 (4 weeks). Functional clustering of these genes
   revealed that the top ranked major biological process was metabolism
   ([89]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
   ([90]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 ([91]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.

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.

   2 wk (187 genes) 4 wk (327 genes)
   GO category # of genes in category Adjusted p-value # of genes in
   category Adjusted p-value
   metabolism 84 3.78e-6 122 1.26e-3
   macromolecule metabolism 58 9.39e-7 76 1.10e-3
   cellular macromolecule metabolism 39 4.83e-4 57 1.50e-3
   protein metabolism 42 1.82e-4 58 3.62e-3
   biopolymer metabolism 40 1.01e-5 52 1.03e-3
   cellular protein metabolism 39 3.57e-4 56 1.81e-3
   protein modification 21 4.79e-3 34 7.59e-4
   biopolymer modification 23 1.67e-3 34 1.59e-3
   cellular catabolism 10 4.08e-3 - -
   micromolecule biosynthesis 14 1.83e-3 - -
   cellular carbohydrate metabolism 7 9.99e-3 - -
   biopolymer catabolism 7 4.57e-3 - -
   cellular macromolecule catabolism 8 1.65e-3 - -
   biosynthesis 19 4.66e-3 - -
   amine metabolism 11 1.10e-4 - -
   aminoacid and derivative metabolism 7 9.30e-3 - -
   macromolecule catabolism 9 1.10e-3 - -
   mRNA metabolism 7 2.26e-3 - -
   RNA processing 8 9.03e-3 - -
   chromatin modification 5 8.91e-3 - -
   regulation of biological process - - 67 5.83e-3
   immune system process - - 20 1.50e-3
   antigen processing and presentation - - 8 1.12e-5
   membrane lipid metabolism - - 7 4.32e-3
   chemotaxis - - 7 6.73e-4
   actin polymerization and/or depolymerization - - 5 4.90e-4
   apoptosis - - 15 9.15e-3
   [92]Open in a new tab

   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. The lists of differentially expressed genes were
   analyzed in WebGestalt ([93]http://bioinfo.vanderbilt.edu/webgestalt)
   for enriched GO categories using hypergeometric test (p<0.01,
   Benjamini-Hochberg); dashes (-) indicate that the corresponding
   categories were not significantly enriched at the respective ages.
   Categories represented by ≥5 genes are shown. Bold text indicates
   categories that were common to both age groups.

Table 6. Enriched molecular function/cellular component Gene Ontology (GO)
categories for the lists of differentially expressed transcripts in spleen
leukocytes from 2 week and 4 week old NOD mice compared to two control
strains.

   2 wk (187 genes) 4 wk (327 genes)
   GO category # of genes in category Adjusted p-value # of genes in
   category Adjusted p-value
   Molecular function
   transferase activity 29 2.54e-3 43 5.94e-3
   catalytic activity 68 3.65e-5 - -
   transcription factor binding 6 8.06e-3 - -
   nuclease activity 5 9.24e-3 - -
   binding - - 167 8.48e-3
   protein binding - - 84 4.69e-3
   zinc ion binding - - 42 2.79e-3
   sugar binding - - 9 3.64e-3
   transferase activity, transferring phosphorous containing groups - - 29
   2.68e-4
   kinase activity - - 23 4.53e-3
   endopeptidase inhibitor activity - - 7 9.64e-3
   Cellular component
   intracellular 81 7.13e-6 131 8.96e-4
   intracellular membrane-bound organelle 64 2.46e-5 - -
   cytoplasm 41 1.43e-3 - -
   intracellular organelle part 28 6.45e-3 - -
   mitochondrion 15 1.36e-3 - -
   nucleus - - 68 8.16e-3
   ribonucleoprotein complex - - 13 6.21e-3
   cell surface - - 10 5.97e-5
   MHC protein complex - - 6 1.76e-5
   [94]Open in a new tab

   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. The lists of differentially expressed genes were
   analyzed in WebGestalt ([95]http://bioinfo.vanderbilt.edu/webgestalt)
   for enriched GO categories using hypergeometric test (p<0.01,
   Benjamini-Hochberg); dashes (-) indicate that the corresponding
   categories were not significantly enriched at the respective ages.
   Categories represented by ≥5 genes are shown. Bold text indicates
   categories that were common to both age groups.

   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 [96]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
   [97]Table 8) also pointing to the same themes as the transcriptome
   data.

Table 7. Enriched KEGG Pathways for the lists of differentially expressed
transcripts in spleen leukocytes from 2 week and 4 week old NOD mice compared
to two control strains.

   2 wk 4 wk
   KEGG Pathway # of genes Adjusted p-value # of genes Adjusted p-value
   Metabolic pathways 19 4.98e-8 21 2.65e-5
   Spliceosome 7 2.79e-6 5 0.0032
   Natural killer cell mediated cytotoxicity 4 0.0027 11 5.56e-9
   T cell receptor signaling 5 0.0004 7 5.80e-5
   Antigen processing and presentation 3 0.0087 7 2.35e-5
   Purine metabolism 5 0.0007 - -
   Amino sugar and nucleotide sugar metabolism 3 0.0013 - -
   RNA degradation 3 0.0042 - -
   ECM-receptor interaction 3 0.0059 - -
   Chronic myeloid leukemia 3 0.0059 - -
   Pyrimidine metabolism 3 0.0087 - -
   Glycine, serine and threonine metabolism 2 0.0093 - -
   Graft-versus-host disease - - 8 9.04e-8
   B cell receptor signaling - - 6 5.80e-5
   Fc gamma R-mediated phagocytosis - - 6 0.0002
   Cell adhesion molecules (CAMs) - - 7 0.0002
   Type I diabetes mellitus - - 5 0.0003
   Autoimmune thyroid disease - - 5 0.0004
   Glycosphingolipid biosynthesis – lacto and neolacto series - - 3 0.001
   Ribosome - - 5 0.0023
   Complement and coagulation cascades - - 4 0.0025
   Phosphatidlylinositol signaling system - - 4 0.0028
   Endocytosis - - 6 0.0033
   Notch signaling pathway - - 3 0.0091
   MAPK signaling pathway - - 6 0.0099
   [98]Open in a new tab

   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 ([99]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.

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.

   2 wk and 4 wk combined (16 proteins)
   KEGG Pathway # of genes Adjusted p-value
   RNA degradation 2 0.0004
   Antigen processing and presentation 2 0.0006
   MAPK signaling pathway 2 0.0032
   Pathways in cancer 2 0.0036
   [100]Open in a new tab

   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. Gene IDs were uploaded and analyzed in WebGestalt
   ([101]http://bioinfo.vanderbilt.edu/webgestalt) for enriched KEGG
   pathways using the hypergeometric test (p<0.01, Benjamini-Hochberg).
   Pathways represented by ≥2 proteins are shown.

   Overall, the KEGG pathway analyses suggest that two broad molecular
   pathways – metabolic and immune response pathways – are 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 base (see “[102]Materials and Methods”
   for details). IPA analyses of the lists of NOD-specific strain
   differences derived at the most stringent threshold
   (Benjamini-Hochberg, p<0.005; [103]Tables 1 and [104]2) generated
   several networks each. Each respective set of networks were merged to
   create one molecular network for each age ([105]Figure 3 and
   [106]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 ([107]Figure 5). That for the 4
   week-old NOD mice clustered around TNF, CASP3, HSPD1, D-glucose and Jnk
   ([108]Figure 6). The major biological functions significantly
   represented by the transcriptome networks included cancer (indicative
   of apoptosis/cellular proliferation) and immunological disease/response
   ([109]Table 9). The top IPA functional categories represented by the
   proteome networks were cell proliferation, cell death and oxidative
   stress ([110]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
   ([111]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
   ([112]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.

Figure 3. Transcriptome network created by Ingenuity Pathway Analysis from
the dataset of 2 week-old mice.

   [113]Figure 3
   [114]Open in a new tab

   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.

Figure 4. Transcriptome network created by Ingenuity Pathway Analysis from
the dataset of 4 week-old mice.

   [115]Figure 4
   [116]Open in a new tab

   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.

Figure 5. Proteome network created by Ingenuity Pathway Analysis from the
dataset of 2 week-old mice.

   [117]Figure 5
   [118]Open in a new tab

   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.

Figure 6. Proteome network created by Ingenuity Pathway Analysis from the
dataset of 4 week-old mice.

   [119]Figure 6
   [120]Open in a new tab

   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.

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.

   Network # of focus genes[121]^a Top IPA functional categories
   Transcriptome – 2 week
   1 13 cancer, cell morphology, immunological disease
   2 13 cell morphology, cellular development, hematological system
   development and function
   Transcriptome – 4 week
   1 17 developmental disorder, cancer, cell death
   2 13 cancer, cell cycle, gastrointestinal disease
   3 12 cellular development, drug metabolism, lipid metabolism
   Proteome – 2 week
   1 9 free radical scavenging, molecular transport, cancer
   Proteome – 4 week
   1 6 cell death, cancer, cell morphology
   [122]Open in a new tab

   One-way ANOVA of expression data sets followed by hierarchical
   clustering identified transcripts/proteins differentially expressed in
   spleen leukocytes from NOD mice compared to two control strains (NON
   and C57BL/6). The lists of differentially expressed genes were analyzed
   by Ingenuity Pathway Analysis ([123]www.ingenuity.com).
   ^a

   Genes that were present in the uploaded lists.

Table 10. Central Genes in the Ingenuity Pathway Networks of differentially
expressed transcripts/proteins in spleen leukocytes of 2 and 4 week old NOD
mice.

   Central Gene Total Connections (# of Focus genes[124]^a) Focus
   Genes[125]^a
   Transcriptome – 2 week
   MYC 20 (5) PDK1,AK2, LXN, PHF5A, VARS
   TNF 13 (4) H60A, PHGDH, AK2, PDZD2
   MYCN 13 (2) RPL8, PHGDH
   TGFB1 11 (2) CAP1, CNN3
   HNF4A 10 (4) CAP1, POLR1C, AK2, CUTC,
   IL2 10 (3) LY6C1, POLR1C, CD244
   HRAS 10 (2) CAPG, LXN,
   IL15 7 (2) TRIM5, CD244
   Progesterone 7 (2) SLC1A5, LXN
   Ca2+ 7 (2) STK25, CNN3,
   PDGFRB 5 (1) SNX6
   YWHAZ 4 (2) C22ORF9, VARS
   Transcriptome – 4 week
   IFNG 22 (9) H60A, LY6C1, PHGDH, CDK6, CD72, HLA-DQB2, HLA-C, PRKCA,
   SELL
   Beta-estradiol 19 (5) BAZ1A, CDK6, HLA-C, PRKCA, TREML2
   TNF 18 (5) H60A, PHGDH, AKAP12, HLA-C, PRKCA
   TP53 18 (6) AKAP12, CDK6, FAM3C, PRKCA, SCAMP1, TEC
   EGFR 17 (4) SNX6, PRKCA, SCAMP1, SNX2,
   NFKB 16 (3) B3GNT5, HLA-C, PRKCA
   Akt 15 (4) PDK1, HLA-C, ICAM2, PRKCA
   PRKCA[126]^a 14 (4) AKAP12, CAPG, PRKCA, SELL
   IL15 12 (5) TRIM5, CDK6, ICAM2, GCNT1, SELL
   HNF4A 11 (7) CAP1, POLR1C, BAZ1A, HLA-C, THEM2, TIA1, TTRAP
   FSH 11 (3) STK25, AKAP12, CDK6
   Progesterone 10 (3) SLC1A5, HLA-DQB2, HOPX
   AR 10 (2) CDK6, ZMIZ1
   IL12 8 (3) LY6C1, SELL, TRIM30
   HLA-C[127]^a 7 (2) HLA-C, KLRA7
   YWHAZ 6 (2) C22ORF9, PRKCA
   CDK6 (FG) 6 (1) CDK6
   Amino acids 4 (4) B3GNT5, PRKCA, PXK, TEC
   Proteome – 2 week
   TP53 9 (4) HSPA8, HSPA9, HSP90AB1, VIM
   MYC 7 (2) LMNB1, TUBA1A
   APP 7 (2) CLIP1, HSPA8
   Ca2+ 6 (2) HSPA8, TUBA1
   IL1B 6 (1) LCP1
   HSPA8[128]^a 5 (2) HSPA8, HSP90AB1
   Proteome – 4 week
   TNF 8 (2) CTNNA, HSPD1
   CAPS3 7 (2) ARHGDI, HSPD1
   D-glucose 7 (2) HSPD1, WDR1
   Jnk 7 (2) HSPD1, PPP5C
   HSPD1[129]^a 5 (2) CTNNA1, HSPD1
   [130]Open in a new tab
   ^a

   Genes that came from the list of differentially expressed genes in NOD
   mice compared to control strains (NON and C57BL/6). The connections to
   each gene (i.e. edges representing gene-gene relationships) in the
   molecular networks were manually counted. The central genes were then
   ranked from highest to lowest based on the number of total connections;
   cut-off, 4 total connections. Genes indicated in Bold text are central
   or focus genes that were common to both ages.

   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 ([131]Table 3)
   was suggestive of immune activation, we conducted IPA on the three
   lists of strain-specific age differences ([132]Tables 3, [133]S2 and
   [134]S3). The functional categories associated with the three topmost
   major networks generated from each of these lists are shown in
   [135]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
   ([136]Figures S4, [137]S5 and [138]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, [139]Figure S4), possibly driven by IFNG and other
   cytokines. In contrast, in C57BL/6 mice, ∼75% of the differentially
   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 ([140]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.

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.

   Network # of focus genes[141]^a Top IPA functional categories
   NOD
   1 15 cellular function and maintenance, immunological disease
   2 13 cell death, DNA replication, recombination, repair, gene
   expression
   3 10 inflammatory response, renal/urological disease, infectious
   disease
   NON
   1 22 hematological system development and function, tissue morphology,
   cell-mediated immune response
   2 20 cell cycle, cancer, skeletal and muscular disorders
   3 17 cellular assembly and organization, cell morphology and cellular
   development
   C57BL/6
   1 19 cellular morphology, cellular assembly and organization, cellular
   function and maintenance
   2 17 cellular movement, immune cell trafficking, cell-to-cell signaling
   and interaction
   3 16 post-translational modification, cell cycle, cellular assembly and
   organization
   [142]Open in a new tab

   Statistical analysis of mRNA expression data by unpaired t-Test
   (p<0.001) followed by Venn diagram analysis identified transcripts
   uniquely differentially expressed in spleen leukocytes between 2 weeks
   and 4 weeks of age in NOD mice and two control strains, NON and
   C57BL/6. The lists of differentially expressed genes were analyzed by
   IPA (Ingenuity Pathway Analysis; [143]www.ingenuity.com).
   ^a

   Genes that were present in the uploaded lists.

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 [144][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 [145][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, [146]Figure 1 and [147]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 [148][14]. This allowed us to put a
   stronger focus on differences unique to each specific age, a
   significant number of which could have been 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 [149][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 non-permissive 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 ([150]Tables 5 and [151]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 ([152]Tables
   7 and [153]8), likewise, identified that abnormalities in metabolic
   pathways were predominantly enriched at 2 weeks, including spliceosome,
   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. Pre-autoimmune metabolic changes have been
   previously reported in female NOD mice that later progressed to
   autoimmune diabetes [154][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 [155][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
   [156][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
   [157][26]; and conversely, ii) that immunostimulation circumvents
   diabetes in NOD mice [158][27]. Notwithstanding that NOD mice have
   important differences in their immune system compared to humans, this
   idea of paradoxical immune suppression appears to have support from
   human type 1 diabetes studies [159][28], [160][29]. Elo et al.
   [161][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
   [162][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 ([163]Table 6) of NOD
   spleen leukocytes coupled with impaired metabolism leads to a defect in
   antigen processing/presentation and T cell stimulation ([164]Tables 5,
   [165]6, [166]7, [167]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 ([168]Tables 1 and [169]2;
   [170]Figure S1), is known to play a role in immune regulatory
   mechanisms. H60a can mediate strong suppressive effects on T cell
   proliferation [171][30]. It has also been shown to play a role in NK
   cell tolerance in experimental autoimmune encephalomyelitis [172][31].
   Additionally, a recent study, using a mouse model of spontaneous
   autoimmune arthritis [173][32], provided evidence suggesting that
   efficient suppression of autoimmune diseases may require polyclonal
   specificities rather than single antigen-specific regulatory T cell
   activity.

   Network analysis by IPA identified several genes/molecules (central
   genes, [174]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 proinflammatory, multifunctional cytokine was present in three
   of the four networks underscoring its importance in the initiation of
   autoimmune diabetes. TNF-α in synergy with IFN-γ (which was present in
   the 4 week transcriptome network) is known to be important in inducing
   autoimmunity in the developing immune system [175][33]–[176][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 TNF-α and another proinflammatory cytokine, IL-6 [177][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 [178][37]. The NFKB
   pathways are essential for the function of antigen presenting cells and
   related defects have previously been reported in T1D patients
   [179][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 [180][39], [181][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 [182][41]. Progesterone is known to
   modulate macrophage activation [183][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 [184][43]–[185][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 [186][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. [187][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 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 [188][48]. Akt was central in the
   transcriptome network at 4 weeks.

   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
   [189][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 [190][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 [191][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
   [192][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 [193][53], [194][54].
   AR is a nuclear receptor transcription factor involved in the
   regulation of many different physiological processes, including
   developmental, immunological, metabolic and homeostatic processes
   [195][55]–[196][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
   androgen-independent 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.

Supporting Information

   Figure S1

   Symbols in Ingenuity Pathway Analysis networks and their meaning.

   (PPT)
   [197]Click here for additional data file.^ (344KB, ppt)
   Figure S2

   Validation of microarray data by quantitative real-time PCR.
   Quantitative real-time PCR was performed for H60a, Trim5, Trim12,
   Ly6c1, Snx6 and Sp110. Fold change values represent expression levels
   of genes in control strains (NON and C57BL/6 mice) relative to NOD
   mice. The results are consistent with the microarray data results.

   (PPT)
   [198]Click here for additional data file.^ (129.5KB, ppt)
   Figure S3

   Sample 2D gel images of protein extracts from spleen leukocytes of one
   control strain (C57BL/6, A) and NOD mice (B). The approximate pI and
   molecular mass in kDa is given on the x and y axes, respectively.

   (PPT)
   [199]Click here for additional data file.^ (304.5KB, ppt)
   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)
   [200]Click here for additional data file.^ (1.1MB, tif)
   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)
   [201]Click here for additional data file.^ (1.1MB, tif)
   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)
   [202]Click here for additional data file.^ (1MB, tif)
   Table S1

   Genes and primers/probes used for quantitative Real-time PCR.

   (DOC)
   [203]Click here for additional data file.^ (29KB, doc)
   Table S2

   Transcripts significantly differentially expressed in spleen leukocytes
   between 2 weeks and 4 weeks of age only in NON mice compared to NOD and
   C57BL/6 mice, and their p-values and fold changes.

   (XLS)
   [204]Click here for additional data file.^ (34KB, xls)
   Table S3

   Transcripts significantly differentially expressed in spleen leukocytes
   between 2 weeks and 4 weeks of age only in C57BL/6 mice compared to NOD
   and NON mice, and their p-values and fold changes.

   (XLS)
   [205]Click here for additional data file.^ (31KB, xls)
   Table S4

   Transcripts significantly differentially expressed in spleen leukocytes
   of 2 week old NOD mice with their ANOVA p-values (Benjamini-Hochberg,
   p<0.05) and fold-changes relative to two control strains, NON and
   C57BL/6. A. 162 probesets of significantly lower expression in NOD
   mice; and B. 25 probesets of significantly higher expression in NOD
   mice.

   (XLS)
   [206]Click here for additional data file.^ (45KB, xls)
   Table S5

   Transcripts significantly differentially expressed in spleen leukocytes
   of 4 week old NOD mice with their ANOVA p-values (Benjamini-Hochberg,
   p<0.05) and fold-changes relative to two control strains, NON and
   C57BL/6. A. 298 probesets of significantly lower expression in NOD
   mice; and B. 29 probesets of significantly higher expression in NOD
   mice.

   (XLS)
   [207]Click here for additional data file.^ (61KB, xls)

Funding Statement

   This work was supported by National Institute of Diabetes and Digestive
   and Kidney Diseases (NIDDK), National Institutes of Health (NIH) grant
   DK62103, and American Diabetes Association (ADA) grant ADA 7-11-BS-49
   to ICG, and NIH instrument grants RR10522 and RR14583 to DMD. The
   funders had no role in study design, data collection and analysis,
   decision to publish, or preparation of the manuscript.

References