Alteration of transcriptomic networks in adoptive-transfer experimental autoimmune encephalomyelitis.

Adoptive transfer experimental autoimmune encephalomyelitis (AT-EAE) is an inflammatory demyelination that recapitulates in mouse spinal cord (SC) the human multiple sclerosis disease. We now analyze previously reported cDNA array data from age-matched young female adult control and passively myelin antigen-sensitized EAE mice with regard to organizational principles of the SC transcriptome in autoimmune demyelination. Although AT-EAE had a large impact on immune response genes, broader functional and chromosomal gene cohorts were neither significantly regulated nor showed significant changes in expression coordination. However, overall transcriptional control was increased in AT-EAE and the proportions of transcript abundances were perturbed within each cohort. Striking likenesses and oppositions were identified in the coordination profiles of genes related to myelination, calcium signaling, and inflammatory response in controls that were substantially altered in AT-EAE. We propose that up- or down-regulation of genes linked to those targeted by the disease could potentially compensate for the pathological transcriptomic changes.

Introduction

Experimental autoimmune encephalomyelitis (EAE), a well-established animal model that recapitulates many clinical and pathophysiological aspects of multiple sclerosis (MS, Baranzini et al., 2005) is a T-cell-mediated inflammatory demyelinating disease of the central nervous system (CNS). Like MS, EAE is physiologically characterized by disturbed axonal conduction leading to motor and sensory impairment. Pathological findings in both MS and EAE include infiltration of activated peripheral inflammatory cells into the CNS, loss of myelin and oligodendrocytes, edema, and axonal damage. Astrocytic hypertrophy and astrogliotic scarring are also prominent features of both MS and EAE, and these changes in astrocytes may contribute to conditions that lead to lasting axonal damage, oligodendrocytic loss, and insufficient or absent remyelination. Adoptive transfer (AT-) EAE is induced by injecting myelin antigen-sensitized T-cells obtained from previously immunized syngeneic donor animals. Although more laborious, this form of EAE has a highly synchronous onset and disease course (Mokhtarian et al., 1984; Pitt et al., 2000). In AT-EAE, myelin basic protein reactive T-cells initiate an immune attack against CNS myelin, characterized by cytokine production and cell death.

There is a rich literature regarding gene expression changes both during normal development and differentiation of myelinating cells (e.g., D'Antonio et al., 2006; Garbay et al., 1998; Gokhan et al., 2005; Jiang et al., 2005) and also the changes that occur acutely (within a day), delayed (within a week), and chronically (several weeks and longer thereafter) following spinal cord (SC) injury (e.g., Carmel et al., 2001; Di Giovanni et al., 2003, see also comprehensive overview in Bareyre and Schwab, 2003). Genomic effects of SC injury have also been compared with other conditions related to CNS demyelinating disease such as MS (e.g., Lock et al., 2002; Lock and Heller, 2003) and mouse-EAE (e.g., Matejuk et al., 2003; Xu et al., 2003). Moreover, several studies have shown altered expression of genes responsible for myelination in a variety of CNS pathologies (e.g., Haroutunian et al., 2007; Kumar et al., 2006).

This report extends the analysis of our previously published cDNA microarray study (Brand-Schieber et al., 2005) in which we compared the transcriptomes of SC from female syngeneic, age-matched control and AT-EAE SJL/J mice at the peak of clinical disability (clinical index 4 = hind- and front-limb paralysis). The disease was induced by injecting MBP-primed immune cells into syngeneic recipients. We reported that, compared to healthy controls, SC of AT-EAE mice at the peak of disability displayed axonal dystrophy, extensive infiltration of the lumbar ventral white matter with CD11β-immunoreactive monocytes, and a threefold decrease in mRNA encoding the gap junction protein connexin43 (Gja1, connexin43, Cx43).

A subset of 3776 distinct genes with known protein products whose expression levels were adequately quantified in all arrays and averaged for all quantifiable spots probing the same gene, was selected for further analysis reported in this paper. In order to test the hypothesis that coordinated expression with Cx43 might account for some of the observed altered patterns of gene expression, we have compared the coordination profile of Gja1 to those of the immune response and myelination genes.

Materials and Methods

Data set

Briefly, 60 μg total RNA, extracted in trizol from each set of two SC of control (C) and AT-EAE (E) mice was reverse transcribed into cDNA using fluorescent dUTPs [Cy3-dUTP (green, g) or Cy5-dUTP (red, r)]. The labeled cDNAs were hybridized overnight at 50°C with four 27 k cDNA microarrays produced by the Microarray Facility of the Albert Einstein College of Medicine (http://129.98.70.229) in the combinations: C1(r)C2(g), C3(r)C4(g), E1(r)E2(g), E3(r)E4(g) (“multiple yellow” strategy, see Iacobas et al., 2006a). The images were acquired and primarily analyzed with GenePix™ Pro 4.1 software (www.axon.com) then normalized according to our in-house developed algorithm. The input data for subsequent analyses were the ratios between the normalized net red/green fluorescence of that spot and the average total net red/green fluorescence of all valid spots in the slide.

The microarray study was performed according to the standards of the Microarray Gene Expression Data Society (MGED) and data complying with the “Minimum Information About Microarray Experiments” (MIAME, Brazma et al., 2001) have been deposited in the National Center for Biotechnology Information Gene Expression Omnibus database (http://www.ncbi.nlm.nih.gov/geo) as series GSE2446.

Detection of significant regulation

We considered a gene as significantly regulated in AT-EAE SC compared to control if the absolute fold change >1.5× and the p-value of the heteroscedastic t-test (two-sample unequal variance) applied to the means of the background subtracted normalized fluorescence values in the four biological replicas of the compared transcriptomes was <0.05. Following previously published procedures (Iacobas et al., 2007a; Iacobas et al., 2007b), we have also tested whether large groups of genes sharing same chromosomal location or encoding functionally similar proteins (hereafter termed cohorts) exhibit specific interactions that might regulate and/or perturb the cohort expression profile. Thus, a gene cohort was considered as significantly regulated if the average expression level of the composing genes was changed by factor larger than 1.5 and the p-value with the Bonferroni adjustment applied to multiple comparisons (Draghici, 2003; Iacobas et al., 2005a; Stekel, 2003) was less than 0.05. The expression of a gene cohort was considered as significantly perturbed if the standard deviation of the fold changes within that group exceeded 1.5. High standard deviation is interpreted as indicating that the proportions of transcript abundances within the cohort were considerably altered, thereby potentially introducing “bottlenecks” in pathway dynamics by perturbed transcriptomic “stoichiometry” (Iacobas et al., 2007a; Iacobas et al., 2007b) which may also change the probability distribution of their outcomes. As in previous papers (Iacobas et al., 2005a; Iacobas et al., 2007a; Iacobas et al., 2007b), the genes were again classified in the following functional cohorts: CSD, Cell cycle, shape, differentiation, death; CYT, cytoskeleton; ENE, energy metabolism; JAE, cell junction, adhesion, extracellular matrix; RNA, RNA processing; SIG, cell signaling; TIC, transport of small molecules and ions into or out of the cells; TRA, transcription; TWC, transport of ions/molecules within the cells; UNK, function not yet assigned.

Variability of transcript abundance and TRA control

The relative estimated (TRA) variability (REV) and the gene expression stability (GES) of both control and experimental specimens were computed as previously described (Iacobas et al., 2003). Since basic transcriptional mechanisms are expected to be similar in a homogeneous set of mice, whereas the local conditions may differ from mouse to mouse, lower REV values are interpreted as indicating lower sensitivity to the local conditions, most probably resulting from increased TRA control, while higher REV values indicate reduced TRA control. Therefore, GES scores reveal the priorities of control over the transcript abundance, with GES = 100 indicating the most stably expressed gene (highest priority in transcript abundance control) and GES = (100/number of quantified genes) the least controlled gene. We organized a database containing the REV and GES values in both conditions and identified the genes for which the system preserved or significantly changed the priorities of TRA control in AT-EAE mice. Both analyses of expression variability and TRA control were extended to the cohorts of genes to test whether the average variability of the cohort was significantly changed in AT-EAE, and whether patterns of hierarchy in expression control among cohorts (significantly different average GES values) were altered in AT-EAE.

TRA coordination

In order to gain insight into whether expression of individual genes is linked to each other so that the protein amounts would respect the “stoichiometry” of the biochemical reactions, we performed analysis of TRA coordination using the pair-wise Pearson correlation coefficient between sets of expression levels in biological replicas (procedure in Iacobas et al., 2005a) of both individual genes and gene cohorts to detect specific interactions that may be responsible for perturbation of functional pathways induced by the AT-EAE. At 5% statistical significance, two genes were considered as synergistically expressed if the Pearson correlation coefficient ρ > 0.9 , antagonistically expressed if ρ > −0.9, and independently expressed if |ρ| > 0.05. For each gene, we determined the synergome, antagome, and exclusome (i.e., sets of synergistically, antagonistically, or independently expressed partners of a given gene), as well as the coordination profile (i.e., the set of the correlation coefficients between expression levels of that gene and each other gene in the biological replicas). Together, the synergome and the antagome of a gene forms its expressome. From the point of view that genes whose protein products work together in functional pathways should be coordinately expressed (Iacobas et al., 2007b), the expressome of a gene indicates the extent of the transcriptomic network related to that gene, while its exclusome indicates the delimitations of this network with respect to other networks. Finally, we compared the coordination profiles of selected genes by computing their overlap (OVL), an indicator ranging from −100 to 100% (Iacobas et al., 2007b), with high positive values indicating likeness, high negative values opposition, and values close to zero indicating neutrality of the coordination profiles. We termed the gene pairs with OVL > 80 or OVL < −80 as “see-saw” partners (Iacobas et al., 2007a).

Results

Expression regulation

As previously reported (Brand-Schieber et al., 2005), we found 1433 regulated genes, encompassing all functional classes and located on all chromosomes, thereby indicating a high degree of complexity of the transcriptomic alterations induced by AT-EAE. Affected pathways were identified using GenMapp and MappFinder database searches (www.genemapp.org; Dahlquist et al., 2002; Doniger et al., 2003). High up-regulated genes were associated with immune and defense response, moderate up-regulated genes were associated with the proinflammatory response, antigen presentation and processing, immune cell migration, and endosome transport. Low (but still significant) up-regulation included genes related to cholesterol metabolism and cytokine biosynthesis. As antigen presentation and cytokine biosynthesis are components of the immune response, our analysis reveals that these related pathways are up-regulated to a major extent in EAE. Most prominently down-regulated GO categories included heterochromatin, acid phosphatase, cytoplasm organization and biogenesis, mitochondrial inner membrane presequence translocase, regulation of coagulation and chemokine, and cytokine-mediated signaling (zinc finger CCCH-type containing 15, Zc3h15 and suppressor of cytokine signaling 5, Socs5).

We found that no single cohort of genes (open bars in Figures 1A and 1B) with regard to chromosomal location or function of the encoded protein was entirely regulated unidirectionally, up- and down-regulations of individual genes being roughly balanced within each cohort. An example is shown in Figure 1C for part of the subcategory JAE. However, this balance did not extend to the subcohorts: for example, the expression levels of 54 of 56 significantly regulated genes (96%) related to immune response (part of the JAE cohort, see Table 1) showed significant increases, with histocompatibility 2 class II antigen A alpha (H2-Aa) exhibiting the highest fold change (26.4×).

Figure 1

Expression regulation of gene cohorts in AT-EAE mouse spinal cord. A–B. Log2ratios (negative for down-regulation) and standard deviation of expression ratios of gene cohorts sharing the same molecular function () or chromosomal location (). Note that while the average fold change was less than 1.5 for all cohorts, the standard deviation of the fold change exceeded 1.5. Highest standard deviations were obtained for JAE genes and for those located on chromosome 17. () Heatmap of the regulated JAE-J genes. C1–C4 = logratios in the four controls, E1–E4 = logratios in the four AT-EAE mice. All logratios were computed with respect to the average expression level in control mice. Red/green/black color indicates up/down/no regulation in the respective set. Note both the reproducibility (green, red, or black color for all four sets of the same condition) and the variability (non-uniform nuance of the color) of the gene expression patterns among the animals of the same condition.

Table 1

Regulation of immune response genes in AT-EAE spinal cord.

Name	Symbol	X	p
Histocompatibility 2, class II antigen A, alpha	H2-Aa	26.37	0.000
Complement component 1, q subcomponent, alpha polypeptide	C1qa	10.03	0.000
Guanylate nucleotide-binding protein 2	Gbp2	9.10	0.000
Complement component 1, q subcomponent, beta polypeptide	C1qb	8.45	0.000
Proteosome (prosome, macropain) subunit, beta type 9 (large multifunctional protease 2)	Psmb9	8.33	0.000
Histocompatibility 2, class II antigen A, beta 1	H2-Ab1	6.67	0.000
Histocompatibility 2, L region	H2-L	6.50	0.000
Histocompatibility 2, K region	H2-K	6.33	0.000
Macrophage expressed gene 1	Mpeg1	6.08	0.000
Fc receptor, IgG, low affinity III	Fcgr3	5.52	0.001
Histocompatibility 2, D region locus 1	H2-D1	5.30	0.001
Chemokine (C–C motif) ligand 9	Ccl9	5.26	0.001
Interferon-induced transmembrane protein 3	Ifitm3	4.46	0.001
B lymphoma Mo-MLV insertion region 1	Bmi1	4.37	0.001
Complement component 3	C3	4.33	0.005
Interferon-induced transmembrane protein 1	Ifitm1	4.23	0.000
Histocompatibility 2, class II, locus DMa	H2-DMa	4.22	0.001
Interleukin 2 receptor, gamma chain	Il2rg	3.55	0.007
Tumor necrosis factor, alpha-induced protein 2	Tnfaip2	3.53	0.000
T-cell receptor alpha, variable 22.1	Tcra-V22.1	3.51	0.000
CD44 antigen	Cd44	3.24	0.000
Proteosome (prosome, macropain) subunit, beta type 8 (large multifunctional protease 7)	Psmb8	3.24	0.001
LPS-responsive beige-like anchor	Lrba	3.22	0.002
SAM domain and HD domain, 1	Samhd1	3.04	0.000
Interferon regulatory factor 1	Irf1	3.04	0.000
Interleukin 7 receptor	Il7r	3.03	0.001
Interferon consensus sequence-binding protein 1	Icsbp1	3.01	0.000
B-cell leukemia/lymphoma 2-related protein A1d	Bcl2a1d	2.99	0.000
Complement component factor h	Cfh	2.91	0.002
Interleukin 1 family, member 6	Il1f6	2.83	0.003
Interferon activated gene 203	Ifi203	2.79	0.001
2′-5′-Oligoadenylate synthetase-like 2	Oasl2	2.71	0.006
SLAM family member 8	Slamf8	2.66	0.000
Fc receptor, IgE, high affinity I, gamma polypeptide	Fcer1g	2.63	0.005
Immunoglobulin superfamily, member 7	Igsf7	2.61	0.004
Interferon-induced transmembrane protein 3-like	Ifitm3l	2.56	0.000
Immediate early response 3	Ier3	2.55	0.027
Chemokine (C–C motif) ligand 22	Ccl22	2.46	0.012
Histocompatibility 2, Q region locus 7	H2-Q7	2.44	0.004
Histocompatibility 2, complement component factor B	H2-Bf	2.44	0.000
Small chemokine (C–C motif) ligand 11	Ccl11	2.41	0.000
IL2-inducible T-cell kinase	Itk	2.39	0.000
Histocompatibility 13	H13	2.38	0.024
Interferon regulatory factor 6	Irf6	2.30	0.001
Proteasome (prosome, macropain) 28 subunit, alpha	Psme1	2.27	0.003
CD47 antigen (Rh-related antigen, integrin-associated signal transducer)	Cd47	2.22	0.000
Chemokine-like factor super family 3	Cklfsf3	2.17	0.002
Interleukin 11 receptor, alpha chain 1	Il11ra1	2.17	0.000
ICOS ligand	Icosl	2.07	0.000
T-complex-associated testis expressed 3	Tcte3	2.07	0.000
CD1d1 antigen	Cd1d1	2.01	0.001
Chemokine (C motif) ligand 1	Xcl1	2.00	0.006
Histocompatibility 2, class II antigen E beta	H2-Eb1	1.78	0.042
Ia-associated invariant chain	Ii	1.76	0.023
Immunoglobulin superfamily, member 8	Igsf8	−2.20	0.000
Integrin alpha FG-GAP repeat containing 1	Itfg1	−2.28	0.000

X, expression ratio (negative for down-regulation) in AT-EAE with respect to control; p, p-value.

Similarity of the regulomes of AT-EAE SC and Cx43 null brain

Since Cx43 was down-regulated in AT-EAE SC by about threefold, we checked whether this alteration of Cx43 expression had similar effects on other genes as observed in Cx43 null brain (Iacobas et al., 2005a), where significantly regulated genes were also located on all chromosomes and encoded proteins of all major functional categories, extending beyond those that might be expected to depend on junctional communication. Indeed, we found substantial OVL between the regulomes of AT-EAE SC and that of Cx43 null brain with respect to their controls (part of which are listed in Table 2), with 84% of the 585 significantly altered genes in both samples exhibiting the same type of regulation.

Table 2

Representative examples of significantly regulated JAE genes in EAE spinal cord (EAE) and in Cx43 null brain (BR).

GB ACC	Name	Symbol	EAE	BR
AA266651	Afamin	Afm	3.29	1.64
AW552701	Amyloid beta (A4) precursor protein	App	−3.00	−5.47
AA544881	Attractin	Atrn	1.65	1.53
AU044759	Cadherin 22	Cdh22	2.51	1.68
AW537209	Calsyntenin 1	Clstn1	−2.26	−2.24
AI573427	Catenin beta	Catnb	−2.26	−2.79
AA154812	Claudin 10	Cldn10	1.96	1.59
AI425965	Contactin 1	Cntn1	−1.87	−1.77
AU040950	Ependymin 2	Epdm2-pending	−4.13	−1.69
W58845	Fc receptor, IgE, high affinity I, gamma polypeptide	Fcer1g	2.63	1.52
AA277329	Histocompatibility 2, class II antigen E beta	H2-Eb1	1.78	1.50
AA061908	Integrin alpha M	Itgam	3.52	1.59
AW545236	Mitochondrial ribosomal protein L28	Mrpl28	−1.72	−1.92
C79334	Musculin	Msc	2.01	1.64
AW537800	N-ethylmaleimide sensitive fusion protein attachment protein alpha	Napa	−4.17	−1.82
AI323974	Neuropilin	Nrp	1.77	1.64
C85793	Putative neuronal cell adhesion molecule	Punc	1.90	1.60
AI414315	Rhomboid, veinlet-like 4 (Drosophila)	Rhbdl4	1.51	1.59
AI327207	Tenascin C	Tnc	2.52	1.78
AA249976	Xeroderma pigmentosum, complementation group C	Xpc	1.64	1.58

Values in right-most columns indicate fold change (negative for down-regulation). Note that although fold changes differ, all listed genes were regulated in the same direction in both experimental groups.

Regulation of expression variability

Figure 1C illustrates the variability of gene expression among animals of the same condition (note the non-uniform color nuances in the heatmap representation). We found that the overall TRA variability increased from 42.5% in the control mice to 49.7% in the AT-EAE mice, indicating a significant (p < 0.0001) loss of the overall TRA control in the diseased mice. This observation of overall increased variability in EAE mice compared to controls was robust for all functional cohorts (Figure 2A) and all chromosomal locations (Figure 2B) with a slight bias toward CYT genes, perhaps reflecting the ongoing alterations in the tissue of these animals. However, both control and EAE specimens had uniform control stringencies among functional categories (Figure 2C) and chromosomal locations (not shown) as indicated by the roughly uniform distributions of the GES scores.

Figure 2

Regulation of the mean relative expression variability (REV) of gene cohorts sharing the same molecular function (A) or chromosomal location (B) in spinal cords of AT-EAE mice compared to control mice. () Average gene expression stability (GES) of functional gene cohorts. Note the uniformity of the two REV distributions both among functional classes: REV(C) = (42.9 ± 2.3)%, REV(E) = (50.2 ± 2.3)% and among chromosomes: REV(C) = (42.5 ± 2.1)%, REV(E) = (49.7 ± 2.3)% and that all functional and chromosomal cohorts became less stably expressed (higher REV) in AT-EAE mice, with the highest increments in the transport into cell (23.8%) and chromosome 13 (34.1%). Cytoskeleton genes were the most unstably expressed in both conditions. Observe the quasi-uniform distribution of GES values among functional categories for both control and AT-EAE specimens.

Table 3 presents examples of very stably and very unstably expressed genes in the control mice that significantly preserved or changed their stability classes in AT-EAE mice.

Table 3

Examples of genes that conserved (SS and UU in STAB column) or significantly changed (SU and US in STAB column) their stability or instability in AT-EAE spinal cord with respect to controls.

STAB	Name	Symbol	CHR	FUNC	GES-C	GES-E	REG
SS	Phosphofructokinase, platelet	Pfkp	13	ENE	98.04	98.17	↓
SS	Burkitt lymphoma receptor 1	Blr1	9	SIG	98.62	98.78	↑
SS	Ethanol induced 2	Etohi2	15	UNK	96.13	95.74	↓
SS	Cadherin 2	Cdh2	18	JAE	99.21	98.57	↓
SS	Ephrin A5	Efna5	17	JAE	97.38	98.04	↓
SS	START domain containing 7	Stard7	2	UNK	95.23	94.17	↓
SS	Bmi1 upstream gene	Bup	2	UNK	96.64	98.07
SS	Gene trap locus 3	Gtl3	8	TRA	98.44	99.92	↓
SS	Muscleblind-like 3 (Drosophila)	Mbnl3	X	CSD	96.72	98.25	↓
SS	Coatomer protein complex, subunit epsilon	Cope	8	TWC	95.47	97.25
SS	Translation factor sui1 homolog	Gc20-pending	9	RNA	99.31	97.35	↓
SS	DiGeorge syndrome critical region gene 2	Dgcr2	16	TWC	96.61	98.81
SS	Ribosomal protein L13	Rpl13	8	RNA	95.66	99.13	↓
SS	High mobility group box 1	Hmgb1	5	TRA	97.19	93.59	↓
SS	ADP-ribosylation factor-like 4	Arl4	12	SIG	98.23	94.44	↑
SU	Ornithine aminotransferase	Oat	7	ENE	96.40	8.85
SU	UDP-N-acetyl-alpha-D-galactosamine:polypeptide N-acetylgalactosaminyltransferase 4	Galnt4	10	ENE	88.29	0.45
SU	Histone cell cycle regulation defective interacting protein 5	Hirip5	6	UNK	93.19	5.11	↓
SU	Chloride intracellular channel 1	Clic1	17	TIC	99.87	11.10
SU	Interleukin-1 receptor-associated kinase 1	Irak1	X	SIG	88.90	0.05
SU	Tubulin, alpha 2	Tuba2	15	CYT	97.83	8.92
SU	RE1-silencing transcription factor	Rest	5	TRA	90.84	1.72
SU	Acetyl-coenzyme A acyltransferase 1	Acaa1	9	ENE	99.97	10.70
SU	Src homology 2 domain-containing transforming protein D	Shd	17	SIG	95.60	6.04
SU	T-complex-associated testis expressed 1	Tcte1	17	TWC	90.86	1.22
SU	Coatomer protein complex, subunit zeta 1	Copz1	15	TWC	95.42	5.75
SU	Sorting nexin 4	Snx4	16	TWC	94.65	4.69
SU	Cleavage and polyadenylation specificity factor 1	Cpsf1	15	RNA	95.21	4.10
SU	Cystathionase (cystathionine gamma-lyase)	Cth	3	ENE	92.24	0.95
SU	Chronic myelogenous leukemia tumor antigen 66	Cml66-pending	15	JAE	94.33	2.97
US	RPEL repeat containing 1	Rpel1	13	UNK	0.69	99.76
US	F-box and WD-40 domain protein 5	Fbxw5	2	UNK	1.69	99.47
US	Ankyrin repeat domain 10	Ankrd10	8	TRA	4.42	98.28	↑
US	Chromobox homolog 3 (Drosophila HP1 gamma)	Cbx3	6	TRA	3.68	96.03
US	DEAD (Asp-Glu-Ala-Asp) box polypeptide 24	Ddx24	12	TRA	4.40	96.00
US	Protein phosphatase 2 (formerly 2A), regulatory subunit B (PR 52), alpha isoform	Ppp2r2a	14	SIG	2.49	95.37
US	Myotubularin-related protein 1	Mtmr1	X	CYT	1.03	95.26
US	GNAS (guanine nucleotide-binding protein, alpha stimulating) complex locus	Gnas	2	SIG	3.20	94.84
US	Coatomer protein complex, subunit gamma 2, antisense 2	Copg2as2	6	TWC	1.06	92.43
US	Four jointed box 1 (Drosophila)	Fjx1	2	UNK	2.15	92.03
US	Glycogen synthase kinase 3 beta	Gsk3b	16	SIG	4.05	91.68
US	BH3 interacting domain death agonist	Bid	6	CSD	1.38	91.42
US	Solute carrier family 14 (urea transprorter), member 1	Slc14a1	18	TIC	0.56	91.39
US	Eukaryotic translation initiation factor 2 alpha kinase 4	Eif2ak4	2	RNA	3.26	90.84
US	Cytochrome b-561	Cyb561	11	ENE	2.99	90.23
UU	Nectin-lke 1	Necl1-pending	1	UNK	0.21	6.94
UU	Neural precursor cell expressed, developmentally down-regulated gene 4	Nedd4	9	ENE	4.98	6.73
UU	Brain protein 17	Brp17	1	UNK	4.82	6.62
UU	MAD homolog 1 (Drosophila)	Madh1	8	SIG	2.89	5.61
UU	Myeloid/lymphoid or mixed-lineage leukemia 5	Mll5	5	TRA	0.13	4.93
UU	Estrogen receptor 1 (alpha)	Esr1	10	SIG	1.54	4.56
UU	Mpv17 transgene, kidney disease mutant	Mpv17	5	ENE	1.80	4.50
UU	Beta-1,3-glucuronyltransferase 3 (glucuronosyltransferase I)	B3Gat3	19	ENE	4.29	2.60
UU	Ring finger protein 103	Rnf103	6	TRA	3.10	2.38
UU	Metastasis suppressor 1	Mtss1	15	CYT	3.15	1.67
UU	Sparc/osteonectin, cwcv and kazal-like domains proteoglycan 1	Spock1	13	JAE	3.73	1.01
UU	A disintegrin and metalloproteinase domain 19 (meltrin beta)	Adam19	11	JAE	2.07	0.72
UU	Guanine nucleotide-binding protein, alpha 11	Gna11	10	SIG	1.19	0.69
UU	Transducin-like enhancer of split 4, E(spl) homolog (Drosophila)	Tle4	19	TRA	0.90	0.64
UU	Seri drox m l transferas mitochon	Shmt2	10	ENE	0.29	0.37

CHR, chromosome location; FUNC, functional class; GES-C and GES-E are the gene expression stabilities in control C and AT-EAE (E) spinal cord while arrows in the REG column indicate whether the gene was significantly up- or down-regulated in E specimens.

The average gene in control SC was found to be synergistically expressed with 467 genes (i.e., 12.4% of 3776), antagonistically expressed with 401 (10.6 %), and independently expressed with 231 (6.1%). As illustrated in Figure 3A, these numbers were not significantly altered in the AT-EAE mice: average synergome size = 436 (11.5%), antagome = 397 (10.6%), and exclusome = 241 (6.4%). However, the networks of genes that were connected to one another were markedly altered, as indicated by the lack of correlation between the control and AT-EAE expressomes and exclusomes. In both conditions, we found remarkably bimodal distributions of coordination frequencies (Figure 3B).

Figure 3

Transcription coordination. () The 95% confidence intervals of the percentages of synergistically (Syn), antagonistically (Ant), and independently (Ind) expressed gene pairs in spinal cord of control (C) and AT-EAE (E) mice. Note that the intervals remained practically unchanged in AT-EAE. () Histograms of the expression coordinations (SYN + ANT) in control and AT-EAE spinal cord. Note the two modal distributions in both conditions, with two distinct groups of genes: one group in which most genes are coordinately expressed with 4–8% of the other genes and the second group in which most genes are coordinately expressed with 40–44% of the other genes. () Examples in which percentages of statistically significant coordination partners of myelination genes differed between control and AT-EAE spinal cords. Chk, choline kinase; Fyn, Fyn proto-oncogene; Hmgcr, 3-hydroxy-3-methylglutaryl-coenzyme A reductase; Mal, myelin and lymphocyte protein T-cell differentiation protein; Mpz, myelin protein zero; Nsmaf, neutral sphingomyelinase activation-associated factor; Pdgfra, platelet derived growth factor receptor alpha polypeptide; Pdgfrb, platelet derived growth factor receptor beta polypeptide; Pdgfrl, platelet-derived growth factor receptor-like; Plp, proteolipid protein; Pmp22, peripheral myelin protein; Qk, quaking; Scd1, stearoyl-coenzyme A desaturase 1; Scd2, stearoyl-coenzyme A desaturase 2; Smpd1, sphingomyelin phosphodiesterase 1 acid lysosomal; Smpdl3a, sphingomyelin phosphodiesterase acid-like 3A.

Some genes exhibited very low numbers of coordinated partners as compared to the average gene; for example, T-complex-associated testis expressed 3 (Tcte3) and P2rx3 have expressomes covering only 2.4 and 2.5% of the sampled transcriptome. In contrast, other genes had large numbers of coordinated partners, such as Spastin (Spast) with an expressome of 45.3% of the transcriptome, and ATP synthase H+ transporting mitochondrial F1 complex gamma polypeptide 1 (Atp5c1), whose expressome encompassed 45.2% of the transcriptome. A large diversity of coordination degrees was found among the genes involved in myelination. Thus, Pmp22 was coordinately expressed with over 43% of the selection, while Nsmaf was coordinated with only 5% of the selection (Figure 3C).

Although coordination of the average gene was not significantly changed in AT-EAE mice, the coordinations within gene cohorts were markedly affected. Thus, the coordination of 13 out of 16 quantified genes related to myelination decreased significantly, while the coordination of stearoyl-coenzyme A desaturases 1 and 2 (Scd1, Scd2) and sphingomyelin phosphodiesterase 1, acid lysosomal (Smpd1) increased significantly.

When the average coordination degree was computed for gene cohorts, we found rather uniform distributions with regard to chromosomal location and functional categories in control as well as in AT-EAE samples (Figures 4A–4D), similar to our findings in brain of the nenatal mouse (Iacobas et al., 2007a). If both coordination and variability of gene expression reflect control mechanisms to promote transcriptome stability, we might expect that these parameters would be mathematically related. Indeed, we found that the coordination decreased exponentially with the expression control (GES), a robust observation for both conditions (Figures 4E and 4F).

Figure 4

Average percentages of synergistically (Syn), antagonistically (Ant), and independently (Ind) expressed gene pairs in functional categories (A and B), chromosomes (C and D) and GES percentiles for control (A, C, E) and AT-EAE (B, D, F) spinal cords. Note the quasi-uniform distributions of the percentages in functional categories and chromosomal locations in both conditions, the not-significant alteration of the percentages in AT-EAE mice and the inverse exponential relationship between synergistic and antagonistic coordination and expression stability (GES percentile).

Table 4 presents the genes with the largest and the smallest synergomes, antagomes, and exclusomes in control SC and the corresponding values in AT-EAE mice, while Table 5 presents the sizes of the synergomes, antagomes, and exclusomes of the quantified immune response genes in controls.

Table 4

Examples of genes with high (H) and low (L) percentages of synergistic (SYN), antagonistic (ANT), and independent (IND) expression partners within the control (C) spinal cord and the corresponding values for AT-EAE (E) mice.

Type	Name	Symbol	CHR	FUNC	SYN-C	ANT-C	IND-C	SYN-E	ANT-E	IND-E
H-SYN	A disintegrin and metalloprotease domain 10	Adam10	9	JAE	28.02	16.76	3.18	16.79	16.55	4.53
H-SYN	Hermansky–Pudlak syndrome 1 homolog (human)	Hps1	19	TWC	27.97	16.87	3.31	18.99	18.64	4.10
H-SYN	Myotubularin-related protein 1	Mtmr1	X	CYT	27.97	17.19	3.39	5.46	1.48	8.10
H-SYN	Potassium channel tetramerisation domain containing 2	Kctd2	11	TIC	27.97	17.08	3.50	8.05	4.18	5.88
H-SYN	Bromodomain containing 3	Brd3	2	TRA	27.94	16.66	3.39	4.77	2.01	14.14
H-SYN	Carbamoyl-phosphate synthetase 2, aspartate transcarbamylase, and dihydroorotase	Cad-pendin	5	ENE	27.94	16.98	3.39	17.51	18.30	3.84
H-SYN	CBF1 interacting corepressor	Cir-pending	2	TRA	27.94	16.90	3.44	14.30	10.99	4.69
H-SYN	Enoyl-coenzyme A hydratase, short chain, 1, mitochondrial	Echs1	7	ENE	27.94	16.74	3.20	23.28	15.04	4.08
H-SYN	Nuclear receptor subfamily 3, group C, member 1	Nr3c1	18	CSD	27.94	16.60	3.26	23.54	16.53	4.10
H-SYN	Nuclear receptor coactivator 1	Ncoa1	12	CSD	27.91	17.16	3.42	8.08	4.90	4.74
L-SYN	Glucose phosphate isomerase 1	Gpi1	7	ENE	0.85	2.38	21.82	4.40	1.43	16.37
L-SYN	Exostoses (multiple)-like 3	Extl3	14	CSD	0.87	2.01	16.82	14.88	23.17	4.05
L-SYN	Phosphoglycerate dehydrogenase like 1	Phgdhl1	14	ENE	0.87	2.75	15.89	5.99	10.51	5.16
L-SYN	Cytotoxic T lymphocyte-associated protein 2 beta	Ctla2b	12	JAE	0.90	1.67	25.85	19.12	14.22	3.55
L-SYN	HLA-B-associated transcript 8	Bat8	17	TRA	0.90	2.38	21.64	17.51	13.06	3.60
L-SYN	Nuclear factor of kappa light chain gene enhancer in B-cells inhibitor, alpha	Nfkbia	12	TRA	0.90	1.69	25.50	12.39	15.39	4.37
L-SYN	Phospholipase D3	Pld3	7	SIG	0.93	2.91	16.18	15.84	16.45	4.74
L-SYN	Acetyl-coenzyme A acetyltransferase 1	Acat1	9	ENE	0.95	2.75	12.24	20.10	17.45	3.84
L-SYN	Fucosidase, alpha-L-1, tissue	Fuca	4	ENE	0.95	1.88	13.45	18.96	10.30	4.45
L-SYN	Glyoxalase 1	Glo1	17	ENE	0.95	2.78	11.81	12.98	21.66	4.16
H-ANT	Spastic paraplegia 7 homolog (human)	Spg7	8	UNK	17.27	28.05	3.50	14.80	23.09	4.00
H-ANT	Protein tyrosine phosphatase, non-receptor type 9	Ptpn9	9	SIG	17.13	28.02	3.42	1.14	3.87	16.58
H-ANT	Dehydrogenase/reductase (SDR family) member 1	Dhrs1	14	ENE	16.95	27.99	3.58	3.95	8.58	6.09
H-ANT	Kinesin family member 23	Kif23	9	TWC	16.84	27.97	3.31	16.63	17.32	4.18
H-ANT	RNA and export factor-binding protein 2	Refbp2	1	TRA	16.92	27.97	3.34	16.74	22.70	3.65
H-ANT	ADP-ribosylation factor 1	Arf1	11	TWC	16.31	27.94	3.23	17.06	23.97	3.73
H-ANT	Fibroblast growth factor (acidic) intracellular-binding protein	Fibp	19	CSD	17.16	27.91	3.50	2.78	2.07	6.41
H-ANT	HECT domain containing 1	Hectd1	12	UNK	17.06	27.91	3.47	2.91	4.48	7.42
H-ANT	ATP synthase, H+ transporting, mitochondrial F1 complex, gamma polypeptide 1	Atp5c1	2	TWC	17.32	27.89	3.44	16.95	18.03	4.58
H-ANT	Nuclear RNA export factor 1 homolog (S. cerevisiae)	Nxf1	19	TWC	17.27	27.89	3.55	5.61	3.63	7.63
L-ANT	B-cell receptor-associated protein 29	Bcap29	12	TWC	2.28	0.85	19.94	22.83	13.98	4.32
L-ANT	Neighbor of Punc E11	Nope	9	TIC	2.22	0.85	20.15	7.42	5.80	4.95
L-ANT	v-rel reticuloendotheliosis viral oncogene homolog A (avian)	Rela	19	TRA	1.83	0.85	17.88	23.73	16.84	4.03
L-ANT	Fat-specific gene 27	Fsp27	6	CSD	2.41	0.87	14.43	19.92	10.81	4.21
L-ANT	Ligand of numb-protein X 2	Lnx2	5	SIG	2.75	0.87	15.86	4.29	6.59	3.34
L-ANT	ras homolog gene family, member H	Arhh	5	SIG	2.20	0.87	17.74	2.54	0.95	13.59
L-ANT	CpG-binding protein	Cgbp-pend	18	TRA	1.85	0.90	23.68	18.51	19.39	4.05
L-ANT	Hydroxyacid oxidase 1, liver	Hao1	2	ENE	2.30	0.90	18.96	10.17	5.72	5.72
L-ANT	miRNA containing gene	Mirg	12	UNK	2.04	0.90	24.36	9.24	5.88	4.74
L-ANT	Sphingosine-1-phosphate phosphatase 1	Sgpp1	12	ENE	2.25	0.90	19.33	12.45	10.54	4.69
H-IND	Anterior gradient 2 (Xenopus laevis)	Agr2	12	UNK	1.03	1.51	26.11	6.99	3.20	7.73
H-IND	Argininosuccinate lyase	Asl	5	ENE	1.67	0.93	25.85	7.36	4.98	6.44
H-IND	Cytotoxic T-lymphocyte-associated protein 2 beta	Ctla2b	12	JAE	0.90	1.67	25.85	19.12	14.22	3.55
H-IND	ATP-binding cassette, sub-family G (WHITE), member 4	Abcg4	9	TIC	1.14	1.38	25.74	6.44	10.20	5.03
H-IND	Sorting nexin 10	Snx10	6	TWC	1.01	1.67	25.72	16.98	23.83	4.10
H-IND	myc-induced nuclear antigen	Mina	16	CSD	1.75	1.03	25.48	20.34	10.91	4.16
H-IND	Acetyl-coenzyme A acyltransferase 2 (mitochondrial 3-oxoacyl-coenzyme A thiolase	Acaa2	18	ENE	1.35	1.30	25.32	17.64	23.57	3.76
H-IND	Ena-vasodilator stimulated phosphoprotein	Evl	12	SIG	1.01	1.77	25.32	8.63	4.45	4.79
H-IND	Solute carrier family 25 (mitochondrial carrier; ornithine transporter), member 15	Slc25a15	8	TWC	1.59	1.11	25.19	11.10	6.25	5.40
H-IND	AE-binding protein 2	Aebp2	6	TRA	1.80	0.95	24.97	1.85	4.82	16.45
L-IND	Acidic (leucine-rich) nuclear phosphoprotein 32 family, member E	Anp32e	3	UNK	26.30	13.22	2.83	4.79	1.48	16.31
L-IND	Fatty acid amide hydrolase	Faah	4	ENE	27.20	14.38	2.83	18.46	9.90	4.56
L-IND	Meiotic recombination 11 homolog A (S. cerevisiae)	Mre11a	9	ENE	26.24	13.11	2.83	24.02	16.98	5.48
L-IND	BCL2-associated athanogene 4	Bag4	8	CSD	25.98	12.95	2.89	15.76	8.18	5.08
L-IND	KDEL containing protein 1	Kdel1	1	TRA	26.22	13.19	2.89	24.02	16.98	3.34
L-IND	Nuclear fragile X mental retardation protein interacting protein	Nufip1	14	TWC	24.55	11.57	2.89	13.16	6.81	5.22
L-IND	T-box 2	Tbx2	11	TRA	26.01	12.71	2.89	19.23	20.37	3.42
L-IND	Disrupted in bipolar disorder 1 homolog (human)	Dibd1	9	TRA	27.09	14.27	2.91	9.19	5.27	5.61
L-IND	Phosphatidylglycerophosphate synthase 1	Pgs1-pend	11	UNK	14.19	27.12	2.91	1.11	2.01	7.71
L-IND	Solute carrier family 4, sodium bicarbonate cotransporter, member 7	Slc4a7	14	TIC	26.22	13.29	2.91	5.01	2.38	10.67

CHR, chromosome location; FUNC, primary function performed by the encoded protein.

Table 5

The extent of the synergomes (SYN), antagomes (ANT), expresomes (EXP = SYN + ANT), and exclusomes (EXC) of the immune response genes in the control spinal cord.

Name	Symbol	SYN	ANT	EXP	EXC
B-cell receptor-associated protein 31	Bcap31	17.46	27.23	44.69	1.91
Interferon-induced transmembrane protein 3	Ifitm3	16.58	27.79	44.37	1.64
Histocompatibility 2, class II antigen E beta	H2-Eb1	27.50	15.84	43.34	1.56
Chemokine-like factor super family 3	Cklfsf3	27.34	15.52	42.86	1.48
B-lymphoma Mo-MLV insertion region 1	Bmi1	16.45	25.72	42.17	1.70
Histocompatibility 13	H13	24.64	17.35	41.99	1.93
CD1d1 antigen	Cd1d1	16.13	24.74	40.87	1.64
CD44 antigen	Cd44	16.87	21.91	38.78	2.04
Proteosome (prosome, macropain) subunit, beta type 8 (large multifunctional protease 7)	Psmb8	23.34	15.18	38.52	1.85
Cd300D antigen	Cd300d	16.50	20.13	36.64	2.04
Chemokine (C–C motif) ligand 22	Ccl22	23.52	10.49	34.01	1.70
Histocompatibility 2, L region	H2-L	18.07	14.52	32.58	2.12
B-cell leukemia/lymphoma 2-related protein A1d	Bcl2a1d	12.13	20.29	32.42	2.17
Small chemokine (C–C motif) ligand 11	Ccl11	11.87	19.36	31.23	1.59
Interferon regulatory factor 6	Irf6	18.57	9.99	28.56	1.91
Interleukin 2 receptor, gamma chain	Il2rg	19.05	8.77	27.81	1.75
Histocompatibility 2, class II, locus DMa	H2-DMa	16.72	8.34	25.06	2.01
Proteosome (prosome, macropain) subunit, beta type 9 (large multifunctional protease 2)	Psmb9	16.72	8.08	24.79	2.09
Histocompatibility 2, class II antigen A, alpha	H2-Aa	12.37	10.62	22.99	1.85
SAM domain and HD domain, 1	Samhd1	12.13	9.64	21.77	1.83
LPS-responsive beige-like anchor	Lrba	14.17	7.02	21.19	2.09
Proteasome (prosome, macropain) 28 subunit, alpha	Psme1	9.56	11.15	20.72	1.80
Fc receptor, IgE, high affinity I, gamma polypeptide	Fcer1g	9.99	9.30	19.28	2.86
Histocompatibility 2, class II antigen A, beta 1	H2-Ab1	12.40	5.85	18.25	2.12
Complement component 1, q subcomponent, alpha polypeptide	C1qa	5.06	12.56	17.62	2.04
Interferon-induced transmembrane protein 3-like	Ifitm3l	9.83	6.86	16.69	2.52
Interferon-induced transmembrane protein 1	Ifitm1	4.03	11.36	15.39	1.83
2′-5′-oligoadenylate synthetase-like 2	Oasl2	10.70	4.03	14.73	2.52
Chemokine (C–C motif) ligand 9	Ccl9	6.83	5.96	12.79	2.75
ICOS ligand	Icosl	6.44	6.15	12.58	2.23
Interferon activated gene 203	Ifi203	8.53	3.39	11.92	2.25
Histocompatibility 2, Q region locus 7	H2-Q7	8.03	3.79	11.81	2.28
Tumor necrosis factor, alpha-induced protein 2	Tnfaip2	7.60	3.55	11.15	1.96
Complement component 1, q subcomponent, beta polypeptide	C1qb	7.50	2.70	10.20	2.46
CD47 antigen (Rh-related antigen, integrin-associated signal transducer)	Cd47	5.48	4.56	10.04	2.86
Histocompatibility 2, K1, K region	H2-K1	5.93	3.74	9.67	2.12
Chemokine (C motif) ligand 1	Xcl1	5.43	2.15	7.58	3.02
T-cell receptor alpha, variable 22.1	Tcra-V22.1	4.53	2.78	7.31	2.97
SLAM family member 8	Slamf8	4.90	2.04	6.94	3.15
Interleukin 7 receptor	Il7r	3.50	1.93	5.43	7.26
Histocompatibility 2, complement component factor B	H2-Bf	1.77	3.31	5.09	7.84
IL2-inducible T-cell kinase	Itk	3.34	1.72	5.06	4.64
Interferon consensus sequence-binding protein 1	Icsbp1	2.70	2.15	4.85	6.78
Complement component factor h	Cfh	2.83	1.93	4.77	6.91
Complement component 3	C3	2.60	2.09	4.69	3.18
Histocompatibility 2, D region locus 1	H2-D1	2.20	1.72	3.92	4.79
Interleukin 1 family, member 6	Il1f6	2.75	1.01	3.76	10.68
Interferon regulatory factor 1	Irf1	2.68	1.09	3.76	8.00
Interleukin 11 receptor, alpha chain 1	Il11ra1	1.14	2.60	3.74	9.17
Guanylate nucleotide-binding protein 2	Gbp2	1.88	1.48	3.36	5.32
Macrophage expressed gene 1	Mpeg1	1.14	1.59	2.73	9.99
Immediate early response 3	Ier3	1.06	1.59	2.65	11.15
Fc receptor, IgG, high affinity I	Fcgr3	1.17	1.40	2.57	7.87
T-complex-associated testis expressed 3	Tcte3	1.01	1.40	2.41	12.29

Table 6 presents examples of genes that preserved or significantly changed the expressome size in AT-EAE SC. Thus, the B-cell receptor-associated protein 31 (Bcap31), that is preferentially associated with the membrane antigen receptor IgD (Adachi et al., 1996) maintained its high coordination degree in both conditions, whereas other genes showed substantial alterations in coordinated expression.

Table 6

Examples of genes that preserved (HH and LL) or significantly changed (HL and LH) the coordination (SYN + ANT) degree in AT-EAE (E) spinal cord as compared to controls (C).

Type	Name	Symbol	CHR	FUNC	SYN-C	ANT-C	IND-C	SYN-E	ANT-E	IND-E
HH	B-cell receptor-associated protein 31	Bcap31	X	TWC	17.45	27.22	3.42	19.25	22.17	3.39
HH	E4F transcription factor 1	E4f1	17	TRA	27.04	17.66	3.47	22.46	18.94	3.65
HH	Ferrochelatase	Fech	18	ENE	27.78	17.27	3.50	21.98	19.23	3.44
HH	Fibroblast growth factor receptor 1	Fgfr1	8	CSD	27.28	17.53	3.50	23.73	17.43	3.84
HH	Gamma-glutamyl carboxylase	Ggcx	6	ENE	27.89	16.79	3.31	22.80	17.61	3.18
HH	Golgi apparatus protein 1	Glg1	8	TWC	27.36	17.24	3.47	23.17	18.19	3.26
HH	Immediate early response 5	Ier5	1	CSD	17.24	27.44	3.36	18.22	23.33	3.39
HH	Mpv17 transgene, kidney disease mutant	Mpv17	5	ENE	27.75	16.84	3.31	24.02	17.32	3.55
HH	Neural proliferation, differentiation and control gene 1	Npdc1	2	UNK	27.86	17.13	3.47	21.61	18.83	3.55
HH	Nuclear cap-binding protein subunit 2	Ncbp2	16	UNK	27.73	17.27	3.55	23.70	18.06	3.36
HH	Placental-specific protein 1	Plac1	X	UNK	17.40	27.20	3.50	23.33	18.59	3.20
HH	Ring finger protein 103	Rnf103	6	TRA	27.81	17.11	3.58	24.05	17.88	3.44
HH	SEC61, gamma subunit	Sec61g	11	TWC	17.40	27.33	3.50	18.06	23.41	3.68
HH	Transducin-like enhancer of split 4, E(spl) homolog (Drosophila)	Tle4	19	TRA	27.54	17.29	3.58	24.02	17.45	3.10
HH	Ubiquitination factor E4B, UFD2 homolog (S. cerevisiae)	Ube4b	4	ENE	27.30	17.29	3.44	21.85	19.07	3.23
HL	Bromodomain containing 3	Brd3	2	TRA	27.94	16.66	3.39	4.77	2.01	14.14
HL	Calcium-binding protein, intestinal	Cai	6	SIG	27.75	17.16	3.28	4.32	1.43	11.94
HL	Capicua homolog (Drosophila)	Cic	7	TWC	16.98	27.73	3.50	0.64	1.93	25.56
HL	Drebrin-like	Dbnl	11	CYT	17.11	27.81	3.65	2.36	2.36	5.61
HL	Ferredoxin reductase	Fdxr	11	ENE	27.09	17.61	3.58	1.56	2.30	7.47
HL	Fibroblast growth factor (acidic) intracellular-binding protein	Fibp	19	CSD	17.16	27.91	3.50	2.78	2.07	6.41
HL	FLN29 gene product	Fln29-pend	5	UNK	27.91	17.43	3.50	0.87	0.82	14.35
HL	HECT domain containing 1	Hectd1	12	UNK	17.06	27.91	3.47	2.91	4.48	7.42
HL	Inhibitor of kappaB kinase epsilon	Ikbke	1	SIG	27.86	17.13	3.52	1.40	0.79	13.51
HL	Myotubularin-related protein 1	Mtmr1	X	CYT	27.97	17.19	3.39	5.46	1.48	8.10
HL	Phosphatidylinositol-4-phosphate 5-kinase, type II, alpha	Pip5k2a	2	ENE	16.90	27.75	3.58	1.69	1.72	17.61
HL	Procollagen C-endopeptidase enhancer 2	Pcolce2	9	JAE	17.16	27.65	3.55	3.60	2.38	15.12
HL	Protein tyrosine phosphatase, non-receptor type 9	Ptpn9	9	SIG	17.13	28.02	3.42	1.14	3.87	16.58
HL	Pyruvate dehydrogenase kinase, isoenzyme 4	Pdk4	6	ENE	27.81	17.21	3.44	2.89	0.85	9.06
HL	von Hippel–Lindau-binding protein 1	Vbp1	X	TWC	16.84	27.83	3.20	1.62	2.67	8.61
LH	Abhydrolase domain containing 3	Abhd3	18	ENE	0.95	1.93	18.88	17.21	23.86	3.20
LH	Acetyl-coenzyme A acyltransferase 2 (mitochondrial 3-oxoacyl-coenzyme A thiolase)	Acaa2	18	ENE	1.35	1.30	25.32	17.64	23.57	3.76
LH	ATP synthase, H+ transporting, mitochondrial F0 complex, subunit c (subunit 9), isoform 2	Atp5g2	15	TWC	1.22	1.62	19.92	17.29	24.02	3.47
LH	Hematopoietically expressed homeobox	Hhex	19	TRA	1.24	1.22	19.15	23.65	17.16	3.84
LH	Immediate early response 3	Ier3	17	CSD	1.06	1.59	22.75	16.23	23.41	4.21
LH	Interleukin-1 receptor-associated kinase 1	Irak1	X	SIG	1.38	1.11	19.54	19.31	22.48	3.36
LH	Polymerase (DNA directed), beta	Polb	8	TRA	1.77	1.11	14.96	23.54	16.26	3.97
LH	RAD9 homolog (S. pombe)	Rad9	19	CSD	1.88	0.95	21.35	23.15	18.86	3.05
LH	Resistin like beta	Retnlb	16	UNK	1.69	1.19	16.39	21.74	19.20	3.50
LH	Retinol-binding protein 4, plasma	Rbp4	19	TWC	1.93	0.98	17.45	23.62	18.30	3.73
LH	RME8 protein	Rme8-pend	9	UNK	1.59	1.11	22.80	16.92	23.09	3.23
LH	Sorting nexin 10	Snx10	6	TWC	1.01	1.67	25.72	16.98	23.83	4.10
LH	Tyrosine 3-monooxygenase/tryptophan 5-monooxygenase activation protein, beta polypeptide	Ywhab	2	TWC	1.48	1.22	12.00	19.17	22.25	3.36
LH	UDP-N-acetyl-alpha-D-galactosamine:polypeptide N-acetylgalactosaminyltransferase 4	Galnt4	10	ENE	1.91	1.01	22.67	23.25	18.64	3.39
LH	v-rel reticuloendotheliosis viral oncogene homolog A (avian)	Rela	19	TRA	1.83	0.85	17.88	23.73	16.84	4.03
LL	AE-binding protein 2	Aebp2	6	TRA	1.80	0.95	24.97	1.85	4.82	16.45
LL	CDC26 subunit of anaphase promoting complex	Cdc26-pen	4	CSD	1.80	1.03	12.66	3.81	2.15	13.72
LL	COP9 (constitutive photomorphogenic) homolog, subunit 4 (Arabidopsis thaliana)	Cops4	5	TRA	1.09	1.59	22.62	3.34	2.94	10.62
LL	Fc receptor, IgG, low affinity III	Fcgr3	1	JAE	1.17	1.40	15.63	2.17	2.44	5.88
LL	Glypican 1	Gpc1	1	JAE	1.22	1.32	18.91	2.67	2.54	5.56
LL	Hairy/enhancer-of-split related with YRPW motif-like	Heyl	4	TRA	1.93	0.93	24.21	0.82	0.74	16.34
LL	Leucine-zipper-like transcriptional regulator, 1	Lztr1	16	TRA	1.77	1.14	20.31	2.17	2.44	4.13
LL	Msx2 interacting nuclear target protein	Mint-pendin	4	UNK	1.62	1.17	21.80	4.63	2.01	14.30
LL	NADH dehydrogenase (ubiquinone) Fe-S protein 1	Ndufs1	1	ENE	1.30	1.24	15.52	2.15	1.40	12.50
LL	Piwi like homolog 2 (Drosophila)	Piwil2	14	CSD	1.30	1.24	24.76	2.73	0.90	6.41
LL	Pre-B-cell leukemia transcription factor 3	Pbx3	2	TRA	1.40	1.51	12.05	2.04	1.88	19.33
LL	Solute carrier family 5 (sodium-dependent vitamin transporter), member 6	Slc5a6	5	TIC	1.11	1.40	16.18	0.50	0.87	18.19
LL	TBP-interacting protein	Tp120a-pe	10	TRA	1.69	1.24	9.38	2.04	2.41	18.03
LL	von Hippel–Lindau syndrome homolog	Vhlh	6	TWC	1.17	1.43	16.10	4.58	1.56	15.20
LL	WD repeat domain 9	Wdr9	16	TRA	1.24	1.59	10.54	1.59	4.13	9.38

H, high coordination; L, low coordination, first symbol in Type column indicating the coordination degree in the C and second one in the E extract; CHR, chromosome location; FUNC, primary function performed by the encoded protein; SYN, synergistic expression; ANT, antagonistic expression; IND, independent expression.

Transcriptomic “see-saws”

We compared the coordination profiles of selected genes involved in immune response, myelination, and calcium signaling by computing their OVL. The rationale for selecting these genes was that EAE is a rodent immune-cell-mediated inflammatory demyelinating disease and that calcium signaling (Iacobas et al., 2006b) is directly related to myelination (Butt, 2006; Fields, 2006). Although as expected, the coordination see-saws were exceptions, as most gene pairs exhibited neutral coordination profiles; in all functional categories, we identified genes with striking likeness or opposition, as illustrated for Pmp22 in Figure 5 in controls. Figure 5 illustrates also how the “see-saw” partnership (here of Pmp22) in control SC was altered by disease.

Figure 5

See-saw partners of peripheral myelin protein (Pmp22) in control spinal cord and alteration by AT-EAE. ( and ) Immune response partners. ( and ) Platelet derived growth factors. ( and ) Calcium signaling genes. The correlation coefficients of the indicated genes with each other gene were plotted against those of Pmp22 with each other gene. In controls, red color indicates likeness, green neutrality, and blue opposition. The table presents the overlap (OVL) of the coordination profiles in both conditions. Note that the AT-EAE turned the significant similarity and opposition of the coordination profiles into neutrality. Cklfsf3, chemokine-like factor super family 3; Ifitm3, interferon-induced transmembrane protein 3; Ier3, immediate early response 3; Pdgfra, platelet derived growth factor receptor, alpha polypeptide; Pdfrl, platelet-derived growth factor, receptor-like; Pdgfc, platelet-derived growth factor, C polypeptide; Itpr1, inositol 1,4,5-triphosphate receptor 1; Grinl1a, glutamate receptor, ionotropic, N-methyl D-aspartate-like 1A; P2rx3, purinergic receptor P2X, ligand-gated ion channel, 3. The overlap scores (OVL) of the Pmp22 in the two conditions were:

Table 7 lists the most strikingly similar or opposite see-saw partners within genes responsible for myelination, calcium signaling, and immune response. The pair H2-DMa: Psmb9 (histocompatibility 2, class II, locus DMa: proteosome (prosome, macropain) subunit, beta type 9 (large multifunctional protease 2)) has the highest likeness, while the pair Scd2:Plcl2 (stearoyl-coenzyme A desaturase 2: phospholipase C-like 2) has the highest opposition as coordination profiles. Pmp22 shared with Pdgfra 956 of its 1039 synergistic partners and 561 of its 593 antagonistic partners, and with Itpr1 929 synergistic and 522 antagonistic partners. In addition, 667 synergistic partners of Pmp22 are antagonistic for Pdgfrl and 979 are antagonistic for Grinl1a, while 442 antagonistic partners of Pmp22 are synergistic for Pdgfrl and 557 are synergistic for Grinl1a. The values of the OVL in legend of Figure 5 confirmed the higher similarity and opposition of the coordination profiles of the illustrated genes.

Table 7

The most like (OVL > 90) and opposite (OVL < −90) see-saw partners among genes involved in myelination (MYE), calcium signaling (CAS), and immune response (IRS) in spinal cord of adult control female mice.

Pair	OVL	Pair	OVL	Symbol	Name
H2-DMa:Psmb9	99.45	Scd2:Plcl2	−99.84	Bcap31	B-cell receptor-associated protein 31
Pdgfra:Pdia4	99.17	H2-Eb1:Ifitm3	−97.48	Bcl2a1b	B-cell leukemia/lymphoma 2-related protein A1b
Pmp22:Cmtm3	98.52	Pmp22:Grinl1a	−95.59	Bmi1	B lymphoma Mo-MLV insertion region 1
Plp:Fyn	97.60	Cd1d1:H13	−94.84	C1qb	Complement component 1, q subcomponent, beta polypeptide
Bmi1:Cd1d1	97.34	Pdgfra:Grinl1a	−94.70	Cab39	Calcium-binding protein 39
Cd44:Cd300d	97.14	Ifitm1:Oasl2	−94.44	Capns1	Calpain, small subunit 1
H2-Ab1:Lrba	96.84	Mal:Tuba8	−94.25	Ccl11	Small chemokine (C–C motif) ligand 11
Pdgfrl:Mal	96.22	Cmtm3:Ifitm3	−94.20	Cd1d1	CD1d1 antigen
Pdgfrb:Cab39	96.19	Bmi1:H13	−93.79	Cd300d	Cd300D antigen
Pdgfrb:H2-Aa	95.44	Pmp22:Ifitm3	−93.57	Cd44	CD44 antigen
Cmtm3:H2:Eb1	95.10	C1qb:H2-Ab1	−93.47	Cherp	Calcium homeostasis endoplasmic reticulum protein
Bcap31:Ifitm3	94.93	C1qb:Lrba	−93.42	Cmtm3	CKLF-like MARVEL transmembrane domain containing 3
Il2rg:Psmb9	94.78	Bcap31:H2-Eb1	−92.57	Fyn	Fyn proto-oncogene
Pmp22:S100a6	94.67	Pdgfra:Plcd1	−92.04	Grinl1a	Glutamate receptor, ionotropic, N-methyl D-aspartate-like 1A
H2-Dma:Il2rg	94.44	Cd44:H13	−92.01	H13	Histocompatibility 13
Pmp22:H2:Eb1	94.01	Bcl2a1b:Irf6	−91.93	H2-Aa	Histocompatibility 2, class II antigen A, alpha
Pmp22:Pdia4	93.90	Pdgfrl:S100g	−91.91	H2-Ab1	Histocompatibility 2, class II antigen A, beta 1
Scd2:Snx9	93.85	Ccl11:Psmb8	−91.89	H2-DMa	Histocompatibility 2, class II, locus DMa
Mpz:Hmgcr	93.82	Pdgfrl:Tuba8	−91.79	H2-Eb1	Histocompatibility 2, class II antigen E beta
Pdgfrb:H2-L	93.70	Scd2:Tuba8	−91.78	H2-K1	Histocompatibility 2, K1, K region
Mal:Qk	93.49	Pmp22:Hspb1	−91.34	H2-L	Histocompatibility 2, L region
Pmp22:Capn3	93.37	Pdgfra:Hspb1	−90.96	H2-Q7	Histocompatibility 2, Q region locus 7
Irf6:Lrba	93.34	Bcap31:Psmb8	−90.51	Hmgcr	3-Hydroxy-3-methylglutaryl-coenzyme A reductase
Pdgfrl:Qk	93.28	Bcap31:Cmtm3	−90.23	Hspb1	Heat shock protein 1
Pmp22:Pdgfra	93.23	Bcap31:Pmp22	−90.16	Ifitm1	Interferon-induced transmembrane protein 1
Pdgfra:S100a6	93.19	H13:Cd300d	−90.02	Ifitm3	Interferon-induced transmembrane protein 3
Pdgfra:Itpr1	92.66			Ifitm3l	Interferon-induced transmembrane protein 3-like
Pdgfrb:Camk2g	92.49			Il2rg	Interleukin 2 receptor, gamma chain
Mal:Snx9	92.39			Irf6	Interferon regulatory factor 6
Pmp22:Itpr1	92.30			Itpr1	Inositol 1,4,5-triphosphate receptor 1
Pdgfrb:Fyn	92.06			Lrba	LPS-responsive beige-like anchor
Pmp22:Capns1	91.58			Mal	Myelin and lymphocyte protein, T-cell differentiation protein
Pdgfra:Capns1	91.54			Mpz	Myelin protein zero
Pmp22:Pde4a	91.54			Oasl2	2′-5′-Oligoadenylate synthetase-like 2
Pdgfra:Plcg2	91.43			Pdgfra	Platelet-derived growth factor receptor, alpha polypeptide
H2:Ab1:Irf6	91.32			Pdgfrb	Platelet-derived growth factor receptor, beta polypeptide
Pdgfra:Capn3	91.31			Pdgfrl	Platelet-derived growth factor receptor-like
Plp:Cab39	90.92			Pdia4	Protein disulfide isomerase-associated 4
Pdgfrl:Snx9	90.90			Plcd1	Phospholipase C, delta 1
Scd2:Cherp	90.66			Plcl2	Phospholipase C-like 2
H2-Aa:H2-L	90.47			Plp	Proteolipid protein (myelin)
Pdgfrl:Cd44	90.45			Pmp22	Peripheral myelin protein
Pdgfrl:S100a11	90.35			Psmb8	Proteosome (prosome, macropain) subunit, beta type 8 (large multifunctional protease 7)
Bcl2a1b:Cd1d1	90.33			Psmb9	Proteosome (prosome, macropain) subunit, beta type 9 (large multifunctional protease 2)
Bcl2a1b:Bmi1	90.32			Qk	Quaking
Plp-H2-Aa	90.30			S100a11	S100 calcium-binding protein A11 (calizzarin)
H2:K1:H2:Q7	90.16			S100a6	S100 calcium-binding protein A6 (calcyclin)
H2-Ab1-Tnfaip2	90.06			S100g	S100 calcium-binding protein G
Pdgfrb:Plp	90.05			Scd2	Stearoyl-coenzyme A desaturase 2
				Snx9	Sorting nexin 9
				Tnfaip2	Tumor necrosis factor, alpha-induced protein 2
				Tuba8	Tubulin, alpha 8

Discussion

Confirmation of array data

Our microarray data are in surprisingly high qualitative agreement with those previously obtained by another group on SC of spontaneous EAE mice through the use of Affymetrix arrays (where some regulations were confirmed by qRT-PCR; Matejuk et al., 2003). As illustrated in Table 8, when common hits were compared between data sets, only 2 of 40 genes (granulin and mevalonate diphospho decarboxylase) showed opposite regulation. Even though magnitudes of expression changes varied between these studies, this qualitative OVL is remarkable, given reports of poor correlations between results provided by different microarray platforms even for the same extracts (e.g., Knudtson et al., 2002).

Table 8

Common hits for EAE spinal cord in this study and that by Matejuk et al., (13) obtained through Affymetrix (Affy) and qRT-PCR (QPCR) techniques.

Name	Symbol	cDNA	p-Value	Affy	p-Value	QPCR
Beta-2 microglobulin	B2m	5.59	0.001	10.8	0.010
Cathepsin C	Ctsc	4.14	0.001	21.1	0.000	12.0
Cathepsin S	Ctss	14.23	0.000	8.3	0.000
Cathepsin Z	Ctsz	5.93	0.000	8.2	0.000
CCAAT/enhancer-binding protein (C/EBP), beta	Cebpb	1.32	0.119	4.6	0.010
CD44 antigen	Cd44	3.24	0.000	6.2	0.020
CD53 antigen	Cd53	1.54	0.404	10.9	0.020
Ceruloplasmin	Cp	4.33	0.012	4.6	0.010
Complement component 1, q subcomponent, alpha polypeptide	C1qa	10.03	0.000	12.2	0.000
Complement component 1, q subcomponent, beta polypeptide	C1qb	8.45	0.000	9.9	0.000
Complement component 3	C3	4.33	0.005	11.4	0.000	8.0
Cysteine-rich protein 1 (intestinal)	Crip1	2.39	0.001	8.8	0.030
EGF-like module containing, mucin-like, hormone receptor-like sequence	Emr1	2.85	0.001	6.3	0.010
Fc receptor, IgG, high affinity I	Fcgr1	1.20	0.388	5.6	0.010	5.0
Granulin	Grn	−2.47	0.003	4.1	0.000
GTP-binding protein 4	Gtpbp4	2.01	0.004	135.8	0.000
Histocompatibility 2, class II antigen A, alpha	H2-Aa	26.37	0.000	424.2	0.000	124.6
Histocompatibility 2, class II antigen A, beta 1	H2-Ab1	6.67	0.000	76.0	0.000
Histocompatibility 2, class II antigen E beta	H2-Eb1	1.78	0.042	66.1	0.000
Histocompatibility 2, D region locus 1	H2-D1	5.30	0.001	15.8	0.000
Histocompatibility 2, K region	H2-K	6.33	0.000	14.4	0.000
Histocompatibility 2, Q region locus 7	H2-Q7	2.44	0.004	19.7	0.000
Ia-associated invariant chain	Ii	1.76	0.023	81.1	0.000
Inositol polyphosphate-5-phosphatase B	Inpp5b	2.74	0.000	15.0	0.000
Interferon activated gene 203	Ifi203	2.79	0.001	4.6	0.020
Interferon consensus sequence-binding protein 1	Icsbp1	3.01	0.000	6.8	0.000
Interferon regulatory factor 1	Irf1	3.04	0.000	10.3	0.000
Lymphocyte specific 1	Lsp1	1.86	0.022	9.2	0.000
Lysosomal-associated protein transmembrane 5	Laptm5	5.49	0.001	5.7	0.000
Mevalonate (diphospho) decarboxylase	Mvd	1.76	0.001	−5.4	0.000
mutS homolog 3 (E. coli)	Msh3	2.11	0.001	5.1	0.010
Proteasome (prosome, macropain) 28 subunit, alpha	Psme1	2.27	0.003	6.3	0.000
SAM domain and HD domain, 1	Samhd1	3.04	0.000	6.8	0.030
Signal transducer and activator of transcription 1	Stat1	6.09	0.000	15.2	0.010	8.0
Signal transducer and activator of transcription 6	Stat6	4.15	0.000	15.8	0.000
Suppressor of cytokine signaling 3	Socs3	1.43	0.017	26.3	0.030
T-cell receptor alpha, variable 22.1	Tcra-V22.1	3.51	0.000	6.0	0.030
Tumor necrosis factor, alpha-induced protein 2	Tnfaip2	3.53	0.000	5.0	0.020
TYRO protein tyrosine kinase-binding protein	Tyrobp	2.80	0.000	7.9	0.000
Vesicle-associated membrane protein 8	Vamp8	3.21	0.001	5.8	0.000

Note that all 40 but 2 (enhanced) genes were found to be regulated in the same sense by both cDNA microarray and Affymetrix studies; exceptions are given in bold. The fold changes obtained by the two platforms was very similar in the case of Laptm5 (98%), Cp (97%), C1qb (92%) and C1qa (90%). Five common hits were also confirmed by Matejuk et al., .

Connexin43 is a key gene of the SC transcriptome

The remarkable OVL between the regulomes of AT-EAE SC and Cx43 null brain indicates that in pathological conditions in which Cx43 is regulated a subset of the altered transcriptome may be attributable to the alteration in Cx43. This finding provides additional support for the hypothesis that Cx43 is a node of gene expression regulation, where its expression is closely tied to that of many other genes (Iacobas et al., 2007a; Iacobas et al., 2007b; Spray and Iacobas, 2007). Although mechanisms responsible for the regulation of other genes by Cx43 expression remain to be completely understood, they likely include regulation of signal molecule exchange between coupled cells and binding to Cx43 cytoplasmic domains of molecules with TRA factor activity (see Kardami et al., 2007).

Gene cohorts were perturbed but not regulated

Cohort analysis of quantified genes has the advantage of providing both manageable and statistically significant ontological information. The expression analysis of gene cohorts is significantly more accurate than that of individual genes due to averaging of technical noise in expression levels of individual genes. In addition, it provides a measure of tendency toward particular patterns of expression. Although we found that no gene cohort was significantly regulated in AT-EAE SC with respect to control (due to a rough balance between up- and down-regulation of individual genes), all cohorts were significantly perturbed, meaning that the disease significantly altered the proportions of transcript abundances within each cohort. Presumably, these “stoichiometric” perturbations contribute to phenotypic alterations that lead to neurological impairment in AT-EAE, including the inflammation, demyelination, axonal damage, and cell death that are associated with infiltration of myelin-specific inflammatory Th1 CD4+ T-cells and, subsequently, activated monocytes into the CNS (Ercolini and Miller, 2006). Interestingly, the highest perturbation was that of JAE genes and chromosome 17 due to the up-regulation of all nine histocompatibility 2 genes involved in immune defense response that were quantified (H2-Aa, H2-Ab1, H2-Eb1, H2-Dma, H2-Bf, H2-D1, H2-K, H2-L, H2-Q7), consistent with the inflammatory nature of the disease. Remarkably, mouse chromosome 17 has the largest homology with human chromosome 6p21, the locus containing the major histocompatibility complex known to influence susceptibility to MS (Yeo et al., 2007).

AT-EAE SC has a loose control of gene expression

Since AT-EAE represents a massive immunological response resulting in gross tissue changes, we examined whether this condition altered the overall inter-animal gene expression variability as compared to healthy mice, expecting that AT-EAE mice would have greater variability due to differences in disease status or in immune attack. Indeed, although for 44.3% of the studied genes the TRA variability decreased, the overall REV value increased from 42.5% in the control mice to 49.7% in the AT-EAE mice, indicating a significant (p < 0.0001) loss of the overall TRA control in the diseased mice. In addition, the sets of REV values in the two conditions were independent (ρ < 10−23) suggesting that multiple mechanisms may contribute to alter the TRA control stringency of individual genes.

Among the most stably expressed genes in control SC whose expression was not affected by AT-EAE was cadherin 2 (Cdh2), which is required for regulating presynaptic function at glutamatergic synapses (Jungling et al., 2006). In contrast, interleukin-1 receptor-associated kinase 1 (Irak1), critical for the induction of EAE as well as for the activation and expansion of autoreactive T-cells (Deng et al., 2003), was very stably expressed in control but became very unstably expressed in AT-EAE, in line with the dynamic nature of the inflammatory process in the disease.

AT-EAE remodels the expression coordination network in SC

The coordination analysis revealed complex interlinkages of gene expression in both control and AT-EAE SC and that the pathology not only alters expression of individual genes but also perturbs functional pathways and rearranges transcriptomic interlinkages However, AT-EAE did not change significantly the average extents of the synergomes, antagomes, and exclusomes, in contradiction with what we have observed in Cx43 null brain were the expressome sizes diminished by more than 20% with respect to the wildtype brain (Iacobas et al., 2007a).

Adam10, the disintegrin and metalloprotease that modulates the cell adhesion role of Pcdh gamma (Reiss et al., 2006) displayed the largest synergome (28.0%), while spastic paraplegia 7 homolog (Spg7), encoding a protein involved in anterograde axon cargo transport (Ferreirinha et al., 2004), had the largest antagome (28%). Remarkably, peripheral myelin protein (Pmp22) was found to be the most coordinately expressed gene involved in myelination (43.2%), while neutral sphingomyelinase activation-associated factor (Nsmaf) was found to be the least coordinated (5.0%). The very large expressome of Pmp22 is consistent with its role in mediating the interaction of Schwann cells with the extracellular environment (Amici et al., 2006; Amici et al., 2007; Berger et al., 2006) and the consequences of its misexpression in generating a family of hereditary peripheral neuropathies (Amici et al., 2007; Sereda and Nave, 2006).

Among the genes responsible for the immune response, interferon-induced transmembrane protein 3 (Ifitm3), involved in the negative regulation of cell proliferation (Ropolo et al., 2004), had the largest expressome (44.4%) and Fc receptor IgG low affinity III (Fcgr3) the smallest (2.5%). Ifitm3 also displayed the largest antagome (27.8%), while histocompatibility 2 class II antigen E beta (H2-Eb1), involved in antigen processing via MHC class II (Alfonso et al., 2001), had the largest synergome (27.5%). Immediate early response 3 (Ier3) appeared to have the largest exclusome (11.2%). Complex gene–gene interactions within the inflammatory pathway have also been reported by other groups. For example, Motsinger et al., 2007 explored 51 single nucleotide polymorphisms (SNPs) in 36 candidate genes within the inflammatory pathway, finding that multi-locus models successfully predicted MS disease risk with high accuracy.

The bimodal distributions of the expressome sizes of SC in both control and disease (Figure 3B), was more prominent than that found in the brain (Iacobas et al., 2005a), and suggest the existence of two categories of genes in terms of coordinations with other genes.

The perspective of “see-saw” partnership

The redundancy provided by similar or opposite coordination profiles may offer the possibility to compensate for functional effects of alteration in gene expression through regulation of interlinked partners. The striking similarity or opposition with regard to coordination profile of Pmp22 with genes such as Pdgfra and Pdgfrl (Figure 5B) may explain why some of the Pdgf genes have been reported to act individually and/or cooperatively in spontaneous remyelination (Murtie et al., 2005). A particularly interesting result of our study was the finding of strikingly similar or opposite partners of Pmp22 in the immune response gene cohort (Cklfsf3 and Ifitm3, Figure 5A), intercellular calcium signaling (inositol 1,4,5-triphosphate receptor 1 (Itpr1) and glutamate receptor ionotropic N-methyl D-aspartate-like 1A (Grinl1a)) (Figure 5C). Itpr1 is responsible for calcium mobilization from the endoplasmic reticulum calcium stores (Iacobas et al., 2006b), while Grinl1a allows calcium ions from the extracellular space to diffuse into the cell when activated by glutamate. The existence of these “see-saw” partners of myelination genes suggests the possibility of correcting myelin defects through stimulating or inhibiting genes responsible for calcium signaling and immune response due to their similar or opposite interlinkage with thousands of genes. These relations between the coordination profiles can be substantially altered in pathological conditions as presented in Figures 5D–5F. These findings of network alterations add a novel concept to expression analysis, in which gene expression profiles can be considered not only by whether expression of genes in a similar functional pathway are affected, but also with regard to how their linkage to one another is altered in a disease state.

We have previously reported (Brand-Schieber et al., 2005) that the gene encoding the gap junction protein Cx43, the most abundant connexin expressed in astrocytes, was among the down-regulated genes in AT-EAE, a result recently confirmed in a guinea pig model of EAE (Roscoe et al., 2007). We concluded that in addition to damage of myelinating glia, altered astrocyte connectivity is a prominent feature of inflammatory demyelination. In order to test the hypothesis that coordinated expression with Cx43 might account for some of the observed altered patterns of gene expression, we have compared the coordination profile of Gja1 to those of the immune response and myelination genes quantified in this experiment. Thus, the coordination profile of Gja1 had a remarkable likeness with those of: histocompatibility 2 Q region locus 7 (H2-Q7, OVL = 87.4%), Psmb9 (OVL = 84.8%), H2-DMa (OVL = 84.5%), interleukin 2 receptor gamma chain (Il2rg, OVL = 83.6%), myelin protein zero (Mpz, OVL = 82.4%), complement component 1 q subcomponent beta polypeptide (C1qb, OVL = 81.6%), and proteolipid protein (myelin) (Plp, OVL = 75.5%). These associations may suggest that the consequences of altered Cx43 expression in damaged SC, evoked by endogenous mechanisms after traumatic SCI (Theriault et al., 1997) might be compensated by overexpression of these other genes. In this regard, we particularly highlight the myelination genes Mpz and Plp because of the prevalence of white matter disturbance in oculodentodigital dysplasia syndrome (Loddenkemper et al., 2002), which is caused by dysfunctional Cx43 mutations (Shibayama et al., 2005). As a corollary of this hypothesis, we might predict that pro-myelinating treatments would result in increased Cx43 expression in SC. Such a result has been recently reported as a consequence of treatment with two drugs in a guinea pig model of EAE (Roscoe et al., 2007).

Conclusion

We found that AT-EAE had a strong impact on the transcriptomic organization of the SC, perturbing most of the functional pathways, relaxing the TRA control, and altering the expression coordination of numerous genes, including those involved in myelination, immune response, and calcium signaling. In addition, identification of what we termed “see-saw” partners may provide alternative therapeutic targets for specific gene-related diseases.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.