Systems analysis of primary Sjögren’s syndrome pathogenesis in salivary glands identifies shared pathways in human and a...

MISSING IMAGE

Material Information

Title:
Systems analysis of primary Sjögren’s syndrome pathogenesis in salivary glands identifies shared pathways in human and a mouse model
Physical Description:
Mixed Material
Language:
English
Creator:
Horvath, Steve
Zazmul-Hossain, Abu NM
Pollard, Rodney PE
Kroese, Frans GM
Vissink, Arjan
Publisher:
BioMed Central (Arthritis research and therapy)
Publication Date:

Notes

Abstract:
Introduction: Primary Sjögren’s syndrome (pSS) is a chronic autoimmune disease with complex etiopathogenesis. Despite extensive studies to understand the disease process utilizing human and mouse models, the intersection between these species remains elusive. To address this gap, we utilized a novel systems biology approach to identify disease-related gene modules and signaling pathways that overlap between humans and mice. Methods: Parotid gland tissues were harvested from 24 pSS and 16 non-pSS sicca patients and 25 controls. For mouse studies, salivary glands were harvested from C57BL/6.NOD-Aec1Aec2 mice at various times during development of pSS-like disease. RNA was analyzed with Affymetrix HG U133+2.0 arrays for human samples and with MOE430+2.0 arrays for mouse samples. The images were processed with Affymetrix software. Weighted-gene co-expression network analysis was used to identify disease-related and functional pathways. Results: Nineteen co-expression modules were identified in human parotid tissue, of which four were significantly upregulated and three were downregulated in pSS patients compared with non-pSS sicca patients and controls. Notably, one of the human disease-related modules was highly preserved in the mouse model, and was enriched with genes involved in immune and inflammatory responses. Further comparison between these two species led to the identification of genes associated with leukocyte recruitment and germinal center formation. Conclusion: Our systems biology analysis of genome-wide expression data from salivary gland tissue of pSS patients and from a pSS mouse model identified common dysregulated biological pathways and molecular targets underlying critical molecular alterations in pSS pathogenesis.
General Note:
Publication of this article was funded in part by the University of Florida Open-Access publishing Fund. In addition, requestors receiving funding through the UFOAP project are expected to submit a post-review, final draft of the article to UF's institutional repository at the University of Florida community, with research, news, outreach, and educational materials.
General Note:
Horvath et al. Arthritis Research & Therapy 2012, 14:R238 http://arthritis-research.com/content/14/6/R238; Pages 1-12
General Note:
doi:10.1186/ar4081 Cite this article as: Horvath et al.: Systems analysis of primary Sjögren’s syndrome pathogenesis in salivary glands identifies shared pathways in human and a mouse model. Arthritis Research & Therapy 2012 14:R238.

Record Information

Source Institution:
University of Florida
Holding Location:
University of Florida
Rights Management:
All rights reserved by the source institution.
System ID:
AA00013478:00001

Full Text


art htis
esearcn nerapy


-- 1


U)
00


o
xT 5-


ln
Module colors


A ,


-~ ~'~L


Systems analysis of primary Sj6gren's syndrome
pathogenesis in salivary glands identifies shared
pathways in human and a mouse model
Horvath et al.
O ioMed Central Horvath et al. Arthritis Research & Therapy 2012, 14:R238
http://arthritis-research.com/content/14/6/R238 (1 November 2012)






Horvath et al. Arthritis Research & Therapy 2012, 14:R238
http://arthritis-research.com/content/14/6/R238


Systems analysis of primary Sj6gren's syndrome

pathogenesis in salivary glands identifies shared

pathways in human and a mouse model

Steve Horvath t, Abu NM Nazmul-Hossainit, Rodney PE Pollard2, Frans GM Kroese2, Arjan Vissink2,
Cees GM Kallenberg2, Fred KL Spijkervet2 Hendrika Bootsma2, Sara A Michie3, Sven U Gorr4, Ammon B Peck5,
Chaochao Cail, Hui Zhoul and David TW Wong"


Abstract
Introduction: Primary Sjbgren's syndrome (pSS) is a chronic autoimmune disease with complex etiopathogenesis.
Despite extensive studies to understand the disease process utilizing human and mouse models, the intersection
between these species remains elusive. To address this gap, we utilized a novel systems biology approach to
identify disease-related gene modules and signaling pathways that overlap between humans and mice.
Methods: Parotid gland tissues were harvested from 24 pSS and 16 non-pSS sicca patients and 25 controls. For
mouse studies, salivary glands were harvested from C57BL/6.NOD AeclAec2 mice at various times during
development of pSS-like disease. RNA was analyzed with Affymetrix HG U133+2.0 arrays for human samples and
with MOE430+2.0 arrays for mouse samples. The images were processed with Affymetrix software. Weighted-gene
co-expression network analysis was used to identify disease-related and functional pathways.
Results: Nineteen co-expression modules were identified in human parotid tissue, of which four were significantly
upregulated and three were downregulated in pSS patients compared with non-pSS sicca patients and controls.
Notably, one of the human disease-related modules was highly preserved in the mouse model, and was enriched
with genes involved in immune and inflammatory responses. Further comparison between these two species led
to the identification of genes associated with leukocyte recruitment and germinal center formation.
Conclusion: Our systems biology analysis of genome-wide expression data from salivary gland tissue of pSS
patients and from a pSS mouse model identified common dysregulated biological pathways and molecular targets
underlying critical molecular alterations in pSS pathogenesis.


Introduction
Sj6gren's syndrome is a chronic, inflammatory autoim-
mune disease characterized by lymphocytic infiltration
of the exocrine glands, especially the salivary and lacri-
mal glands, leading to destruction of their functional
components. The disease may exist alone as primary
Sj6gren's syndrome (pSS) or in conjunction with
another autoimmune disorder as secondary Sj6gren's
syndrome [1]. The disease affects 0.5 to 1.0% of the

* Correspondence dtww@uclaedu
t Contributed equally
'School of Dentistry, Dental Research Institute, University of California at Los
Angeles, 10833 Le Conte Avenue, 73-017 CHS, Los Angeles, CA 90095-1668,
USA
Full list of author information is available at the end of the article


O BioMed Central


general population and shows a striking 9:1 female pre-
dominance [2,3]. Although dry mouth and dry eyes are
the hallmark symptoms of pSS, the disease can affect
almost any organ of the body and can cause substantial
morbidity [4]. Symptomatic therapy currently dominates
the overall management of affected individuals [5,6].
The diagnostic criteria, including the American-Eur-
opean Consensus Group (AECG) criteria involving
either serology or histopathology in conjunction with
oral and ocular dryness [7], are often inconsistently
applied, and are not optimal in differentiating pSS from
non-pSS sicca patients or addressing extraglandular
manifestations. As a result, diagnosis of Sjogren's syn-
drome often lags disease onset by 6 to 10 years. Recent


2012 Horvath et al.; licensee BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Common,
Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in
any medium, provided the original work is properly cited.


0 &






Horvath et al. Arthritis Research & Therapy 2012, 14:R238
http://arthritis-research.com/content/14/6/R238


attempts to address these issues the EULAR Sj6gren's
Syndrome Patient Reported Index, and a disease activity
index to evaluate systemic complications (EULAR Sj6g-
ren's Syndrome Disease Activity Index) [8,9] require
further large-scale validation. In addition to the difficulty
in disease diagnosis, the underlying pathophysiologic
mechanisms of Sj6gren's syndrome remain obscure.
A common model suggests interaction between
genetic susceptibility and environmental factors such as
viral infections for the immune cell activation and pro-
tracted inflammatory response resulting in glandular
dysfunction and systemic autoimmunity [6]. In addition,
recent studies in human and animal models have high-
lighted several components of the innate and adaptive
immune systems as well as nonimmunologic factors
[10-13]. However, there is a general lack of experimental
and research intersections between humans and mouse
models in terms of biological pathways and key molecu-
lar targets. Indeed, a wide variety of animal models of
Sj6gren's syndrome have been described; each capturing
certain aspects of the disease [14]. The most widely
used mouse model, NOD, exhibits CD4+ lymphocyte
infiltration, autoantibodies, xerostomia and a female-
dominant phenotype. However, the mice also develop
diabetes, which complicates the search for Sj6gren's syn-
drome specific changes of gene expression. Recently, the
modified NOD mouse model C57BL/6.NOD-AeclAec2
has been described. These mice develop the symptoms
of Sj6gren's syndrome but not diabetes [15], suggesting
that they will be a highly valuable model for understand-
ing the pathogenesis pSS.
In the present study, we applied weighted gene co-
expression network analysis (WGCNA) [16] to charac-
terize biological pathways and molecular targets asso-
ciated with pSS pathogenesis in both human parotid
tissue and the C57BL/6.NOD-AeclAec2 mouse model.
Our systems analysis of genome-wide expression data
from human (pSS) salivary gland tissue compared with
salivary gland tissue from a mouse model of pSS identi-
fied common biological pathways and molecular targets
that can pivotally contribute to critical molecular altera-
tions in pSS pathogenesis.

Materials and methods
Human studies
The study protocol was approved by the Institutional
Review Board of the University of California at Los
Angeles and the University Medical Center at Gronin-
gen. All patients were recruited from the Department of
Oral and Maxillofacial Surgery of the University Medical
Center at Groningen, were at least 21 years old, and
provided written informed consent to participate in the
study. The study included 24 patients who fulfilled the
2002 AECG criteria for pSS [17], 16 non-sicca patients


with subjective symptoms and objective signs of oral
and ocular dryness but not fulfilling the AECG criteria
for pSS, and 25 control patients with no subjective or
objective evidence of oral or ocular dryness. The patient
characteristics are presented in Table 1.
Parotid gland tissue specimens
Under local anesthesia, an incisional biopsy of one paro-
tid gland was performed as part of the diagnostic
workup of the pSS and non-pSS sicca patients following
AECG criteria, as described [7], and according to the
procedure described by Pijpe and colleagues [18]. Tissue
was harvested from the dorsal caudal lobe of the parotid
gland from the control patients during surgery for oral
or oropharyngeal squamous cell carcinoma.
After harvest, most of the specimen was snap-frozen
and delivered to the University of California at Los
Angeles on dry ice and then stored at -80C for gene
expression profiling. A minor part was fixed in 4% for-
malin, embedded in paraffin, sectioned, and stained with
H & E for histopathologic evaluation. The evaluation of
slides from pSS and non-pSS sicca patients was per-
formed independently by two oral pathologists, who
determined the focus score (>50 lymphocytes/4 mm2
glandular tissue) and other characteristics such as lym-
phoepithelial lesions, germinal centers, fibrosis, atrophy,
and so forth. The specimens from pSS patients showed
characteristic features (for example, focus score >1 and
presence of lymphoepithelial lesions) whereas those
from non-pSS sicca and control subjects did not show
any such characteristics. Histopathologic evaluation of
parotid tissue from controls revealed no carcinoma.
Gene expression profiling
The isolation of total RNA from snap-frozen parotid
gland tissues, in sample sizes ranging from 10 to 30 mg,
was performed using the RNeasy kit supplied with the
PARIS system (Ambion, Austin, TX, USA), according to
the manufacturer's protocol. Briefly, each frozen tissue
specimen was homogenized in 300 pl cold cell-disrup-
tion buffer and an equal volume of 2 x denaturing solu-
tion was added. The total lysate was centrifuged for 5
minutes at 10,000 x g to separate out the superficial
aqueous phase containing RNA and then 1.25 volumes
of 100% ethanol were added. The mixture was loaded
onto a spin column followed by washing twice with pro-
vided buffer. RNA was eluted at 100'C into 50 pl elution
buffer. The resultant RNA was subjected to RNase-free
DNase treatment followed by ethanol precipitation, dis-
solved into 15 pl DNase/RNase-free water, and quanti-
fied with a Nanodrop spectrophotometer (Nanodrop
Technology, Wilmington, DE, USA).
The RiboAmp RNA Amplification system (Molecular
Devices, Sunnyvale, CA, USA) was used to perform one
round of linear amplification of parotid gland tissue
mRNAs, taking 2.5 pg total RNA from each specimen as


Page 2 of 12







Horvath et al. Arthritis Research & Therapy 2012, 14:R238 Page 3 of 12
http://arthritis-research.com/content/14/6/R238



Table 1 Characteristics of primary Sjogren's syndrome, nonprimary Sjogren's syndrome sicca and control subjects
Characteristic pSS (n = 24) Non-pSS sicca (n = 16) Controls (n = 25)


Age (years)
Female:male
Ethnicity
Unstimulated whole saliva (ml/minute)
Stimulated whole saliva (ml/minute)
Schirmer's test (mm)


49.0 + 15.8
21:3
24 Caucasiar
0.10 + 0.14
0.31 + 0.35
9.2 + 76


52.8 + 15.4
13:3
16 Caucasia r
0.12 + 0.18
0.33 + 0.40
148 + 11.7


59.0 + 12.0
12:13
25 Caucasial


Focal score 3.5 + 1.2 0.5 + 1.0 0
syndrome A (positive:negative) 23:1 2:14 0:25
syndrome B (positive:negative) 19:5 0:16 0:25
Data presented as the mean + standard deviation. All primary Sjogren's syndrome (pSS) and non-pSS sicca patients were subjected to a complete American-
European Consensus Group diagnostic work-up. NA, not assessed (in controls being head and neck cancer patients without parotid involvement, no functional
tests were performed as all these patients were already scheduled for a neck dissection and functional tests were considered too much of a burden to the
patients at that time).


the template. The synthesized cDNA was transcribed to
cRNA and then biotinylated using the GeneChip
Expression 3'-Amplification reagents (Affymetrix, Santa
Clara, CA, USA) for in vitro labeling. Then 15 tg
labeled cRNA was fragmented into 50 to 200 bp frag-
ments and assessed for quality with a 2100 Bioanalyzer
(Agilent Technologies, Palo Alto, CA, USA).
The fragmented biotin-labeled cRNA from pSS, non-
pSS sicca and control subjects was hybridized overnight
to the Affymetrix HG U133 plus 2.0 arrays. After wash-
ing to remove the unbound transcripts, the hybridized
chips were stained prior to scanning. The acquired
images were processed with Affymetrix Microarray Suite
software and analyzed with the R software (which is
freely available at [19]).

Mouse studies
Mouse model of Sjogren's syndrome-like disease
C57BL/6.NOD-Aec1Aec2 mice were bred and main-
tained under specific pathogen-free conditions within
the mouse facility of the Department of Pathology at the
University of Florida, Gainesville. The breeding and use
of these animals for the present studies were approved
by the University of Florida Institutional Animal Care
and Use Committee. The animals were maintained on a
12-hour light-dark schedule and provided with food and
acidified water ad libitum.
Generation of salivary gland transcriptome data
Salivary glands were freshly excised from euthanized
mice (n = 5 per age group) at 4, 8, 12, 16, or 20 weeks
of age, snap-frozen in liquid nitrogen, and stored at -80
C until all glandular samples were obtained. Each sali-
vary gland was comprised of a submandibular, sublin-
gual, and parotid gland minus salivary lymph nodes.
Total RNA was isolated from salivary glands of each
mouse using the RNeasy Mini-Kit (Qiagen, Valencia,
CA, USA), in accordance with the manufacturer's proto-
col. Hybridizations were carried out with each of the 25


individual RNA samples using Affymetrix GeneChip
Mouse Genome 430 plus 2.0 Arrays, in accordance with
the manufacturer's instructions. Each GeneChip con-
tained 45,000 probe sets that analyzed the expression
level of over 39,000 transcripts and variants from over
34,000 well-characterized mouse genes.

Statistical procedures
Microarray data preprocessing
The Affymetrix U133 plus 2.0 microarray data were ana-
lyzed with functions and packages of the statistical soft-
ware R 2.12.0 and Bioconductor 2.7 (which can be
downloaded from [19]). Expression intensity values were
calculated for the 65 Human microarray CEL files (24
pSS, 16 non-pSS sicca and 25 controls) using the mas5
function of the Affy library. Potential array outliers were
identified by studying their interarray correlation. After
removing potential outliers in a completely unbiased
fashion, all 65 samples remained for downstream analy-
sis. After quantile normalization, the ComBat function
was used to correct for batch effects [20]. To summarize
multiple probe sets per gene, we used the default set-
tings of the collapseRows R function [21]. Each of the
20,718 genes on the array was thus represented by the
probe with the highest mean expression value.
Analogous preprocessing steps were applied to the 25
mouse microarray samples corresponding to five differ-
ent time points (weeks 4, 8, 12, 16 and 20, each time
point had five mice). We related mouse genes to human
data using a table of orthologous genes.
Weighted gene co-expression network analysis and
preservation analysis
The statistical analysis software (WGCNA R package)
and R tutorials for constructing a weighted gene co-
expression network can be found in the literature
[16,22]. The WGCNA package first calculates all pair-
wise Pearson's correlation coefficients across all samples.
In a signed weighted network, the resulting Pearson's






Horvath et al. Arthritis Research & Therapy 2012, 14:R238
http://arthritis-research.com/content/14/6/R238


correlation matrix is transformed into an adjacency
matrix [23]:

a(ij) = |0.5 + 0.5 cor(x(i),x(j))| )

The default value of the power f3 = 12 facilitates a
soft-thresholding approach that preserves the continu-
ous nature of the co-expression relationships [16]. As a
network dissimilarity measure we used 1 the topologi-
cal overlap measure as the input for average linkage
hierarchical clustering [24]. We used the dynamic
branch cutting method to define modules as branches of
the hierarchical clustering tree [25]. Unassigned back-
ground genes, outside each of the modules, were
denoted with the color grey. To group the 20,718 genes
into modules, we used the blockwiseModule R function
in the WGCNA R package.
Connectivity and module membership measures
Module membership (MM), also known as eigengene-
based connectivity, is a measure of intramodular con-
nectivity [26]. MM is defined as:

MM(i) = cor(x(i), ME)

where x(i) is the expression profile of ith gene and ME
is the eigengene (first principal component) of the given
module. We used the MM measure to select module
genes for a gene ontology (GO) enrichment analysis.
The MM measures of the human modules are reported
in Additional file 1.
Functional enrichment analysis
The Database for Annotation, Visualization, and Inte-
grated Discovery (available at [27]) and the Ingenuity
Pathways Analysis software (Ingenuity Systems available
at [28]) were used to determine whether sets of genes
(for example, preserved intramodular hub genes) were
significantly enriched with known GOs. The Ingenuity
software only reports uncorrected P values. The func-
tional enrichment results are reported in Additional file
2.
Module preservation analysis
To evaluate which human modules could also be found
in the mouse data, we used module preservation statis-
tics implemented in the modulePreservation R function
[17]. For each module in the reference data (for exam-
ple, human data), a permutation test leads to the Zsum-
mary statistic in the test data (here the mouse dataset).
Zsummary >10 indicates strong evidence of preservation
while Zsummary <2 indicates no preservation according
to the permutation test.
Standard differential expression analysis
To find differentially expressed genes between pairs of
groupings (involving pSS, sicca, controls), we used two
R functions in the WGCNA package (standardScree-
ningBinaryTrait and standardScreeningNumericTrait)


that report P values, false discovery rates (q values), fold
changes and other widely used statistics for selecting
differential expressed genes. Additional file 3 reports the
differential expression results in humans for studying
contrasting pSS versus controls, sicca patients versus
controls, and pSS versus sicca patients, and for relating
expression data to gender in control subjects. Further,
Additional file 3 reports the correlation of each gene
with an ordinal measure of diagnosis (0 = controls, 1 =
sicca, 2 = pSS), which allows the reader to identify
genes that positively or negatively increase with the dis-
ease progression in humans. We also correlated each
mouse gene with time (weeks) to find genes that are
increasing or decreasing as time progresses. These
results (P values, correlations) are reported in Additional
file 4. These results were also used to create Table 2.

Results
Identification of gene co-expression modules in parotid
glands of pSS patients
A signed weighted gene co-expression network was con-
structed based on the 65 human parotid gland tissues
(24 pSS, 16 non-pSS sicca and 25 controls). The
WGCNA method clustered the 20,718 human genes
into 19 distinct gene co-expression modules. Since the
module detection is unbiased and does not make use of
GO information, each of the modules was initially
labeled with a unique color as an identifier (Figure 1).
To define a representative module expression profile
(referred to as the module eigengene), we summarized
the (standardized) gene expression profiles of the mod-
ule eigengene (= first principal component). The module
eigengene can be considered a weighted average of the
module gene expression profiles. To identify disease-
related modules, we correlate each module eigengene
with disease status.
Strikingly, seven out of 19 modules showed significant
differential expression between pSS and control samples,
which reflects the fact that thousands of genes are dif-
ferentially expressed between the two groups. In particu-
lar, we found a highly significant positive correlation
between disease and the Magenta (comprised of 576
genes) module eigengene (r = 0.67, P <2 x 10-7), the
Brown (2,502 genes) module eigengene (r = 0.6, P <6 x
10-6), the Light-Cyan (349 genes) module eigengene (r =
0.42, P <0.002), and the Grey60 (based on 216 genes)
module eigengene (r = 0.47, P <8 x 10-4). Since a posi-
tive correlation indicates that the corresponding module
genes are overexpressed in the diseased individuals, the
Magenta, Brown, Light-Cyan and Grey60 modules are
comprised of genes overexpressed in pSS samples com-
pared with both non-pSS sicca and controls (Figure 2).
On the contrary, we found a highly significant negative
correlation between disease status and the Turquoise


Page 4 of 12







Horvath et al. Arthritis Research & Therapy 2012, 14:R238
http://arthritis-research.com/content/14/6/R238


Table 2 Genes with concordant progression patterns in human and mouse disease
Progression GeneSymbol corHuman p.Human corMouse
+ ADA 0.57 5.5 x 107 0.62
+ AIF1 0.51 1.6 x 10 0.57
+ C1 QB 0.54 3.2 x 106 0.65
+ CASP3 0.59 1.8 x 107 0.52
+ CYBB 0.63 1.6 x 10' 0.56
+ ENTPD1 0.51 1.2 x 10 0.62
+ GZMA 0.63 1.7 x 10 0.70
+ GZMK 0.61 5.8 x 10 0.55
+ HLA-DQB1 0.53 6.0 x 10 0.65
+ HLA-DRB 0.53 5.9 x 106 0.74
+ IFIT1 0.59 2.6 x 107 0.59
+ ITGAX 0.51 1.2 x 10 0.63
+ LY86 0.58 4.4 x 107 0.53
+ PLEKHA2 0.63 1.5 x 10 0.52
+ STAT1 0.67 7.2 x 10 0 0.64
+ TLR7 0.63 1.6 x 10 0.63
ALDH3A2 -0.50 1.8 x 10 -0.63
ANGPTL4 -0.50 1.9 x 10 -0.58
ATP1B1 -0.56 1.0 x 106 -0.58
CPT1A -0.54 4.0 x 0-6 -0.76


p.Mouse


Page 5 of 12







MM.magenta


NEUI -0.52 9.6x 10" -0.53 6.9 x 10 -0.47
PALMD -0.50 2.2 x 10 -0.64 6.1 x 10-4 -0.65
PDHA1 -0.58 39 x 107 -0.63 8.0 x 104 -0.71
PPFIBP2 -0.56 1.2 x 06 -0.57 3.2 x 10' -0.48
corHuman denotes the correlation of the gene with human progression status (defined as 0 for controls, 1 for sicca, and 2 for primary Sjogren's syndrome).
corMouse denotes the correlation of the gene with mouse progression status (weeks). p.Human and p.Mouse report the correlation test P values in human and
mouse data, respectively. MM.magenta reports the module membership value of the gene.


O



So 0



















Figure 1 Hierarchical cluster trees for detecting modules based on weighted gene co-expression network analysis All array samples
(primary syndrome patients, sicca, and controls) were to find co-expression modules. Modules correspond to branches of the trees and
are assigned color categories as indicated by the color band underneath each tree (see Materials and methods). The two panels correspond to
two cluster trees resulting from the blockwiseModules WGCNA R function [21], which was used to circumvent computational challenges. Each
cluster tree thus corresponds to a block of genes.
clster tree thus corespondstoabock of genes.






Horvath et al. Arthritis Research & Therapy 2012, 14:R238
http://arthritis-research.com/content/14/6/R238


p= 9.1e-09


C
I
C)
wJ
aB


Page 6 of 12


p= 4.6e-06


diseaseProg


diseaseProg


p= 0.01


I-.


p= 0.014


diseaseProg


diseaseProg


Figure 2 Identification of primary Sj6gren's syndrome disease-related gene modules revealed by weighted gene co-expression
network analysis. Barplots showing the mean value of the module eigengene (y axis) versus human disease progression status (x axis) where 0
denotes healthy controls, I denotes nonprimary : syndrome (non-pSS) sicca patients, and 2 denotes pSS patients. Kruskal-Wallis test P
values above the plots show that the module eigengenes of all four modules Magenta (576 genes), Brown (2,502 genes), Grey60 (216 genes),
and Light-Cyan (349 genes) are I overexpressed in pSS patients compared with both non-pSS sicca patients and controls.


(based on 3,981 genes) module eigengene (r = -0.52,
P <10 -4), the Grey (based on 3,011 genes) module eigen-
gene (r = -0.41, P < 0.003), and the Salmon (based on
446 genes) module eigengene (r = -0.28, P <0.005). The
Turquoise, Grey and Salmon modules are thus comprised
of genes that are underexpressed in pSS samples com-
pared with both non-pSS sicca and controls (Figure 3).
All of these P values remain highly significant even after
carrying out the most stringent multiple comparison


adjustment (Bonferroni correction) for the number of
modules. This co-expression module-based analysis has a
major advantage over a standard differential gene expres-
sion analyses since it only relates 19 modules to pSS dis-
ease status, alleviating the multiple comparison problem
inherent in the raw data.
For each gene, Additional file 1 reports MM measures
(also known as eigengene-based connectivity) between
genes and modules.


q1T







Horvath et al. Arthritis Research & Therapy 2012, 14:R238
http://arthritis-research.com/content/14/6/R238


p= 0.00068


Page 7 of 12


p= le-04


diseaseProg


diseaseProg


p= 0.00029


0 1 2
diseaseProg
Figure 3 Primary Sj6gren's syndrome disease-related co-expression modules that are downregulated in primary Sj6gren's syndrome
patients. Barplots analogous to those described in Figure 2. The module eigengenes of two modules Turquoise (3,981 genes) and Grey (3,011
genes) are I underexpressed in primary syndrome (pSS) patients compared with both non-pSS sicca patients and controls.
Note that the Salmon module (446 genes) differs from the other two modules with respect to the non-pSS sicca patients, where the module
eigengene is underexpressed compared with controls.


Preservation of human Magenta module in a mouse
Sjogren's syndrome model
We applied the modulePreservation R function to assess
whether modules defined by the human data were pre-
served in C57BL/6.NOD-AeclAec2 mice, Sj6gren's syn-
drome-like disease model datasets. We selected the


C57BL/6.NOD-AeclAec2 mice for the comparative gene
expression analysis because it is a spontaneous disease
model, probably the best-defined mouse model to date
with respect to its disease profile, and has a comparative
control with the same genetic background that mimics
virtually all molecular and clinical aspects defined in






Horvath et al. Arthritis Research & Therapy 2012, 14:R238
http://arthritis-research.com/content/14/6/R238


humans. Importantly for this study, extensive transcrip-
tomic profiling data of salivary glands defining develop-
ment of Sj6gren's syndrome are available.
The module preservation statistics allow one to quantify
which aspects of within-module topology are preserved
between a reference network (human parotid glands) and
a test network (mouse salivary glands). We used two com-
posite preservation statistics (Zsummary and medianRank)
to assess overall module preservation. We evaluated the
preservation of the human modules in the entire mouse
dataset (that is, all weeks together). Strikingly, both statis-
tics revealed that the human Magenta module is the most
highly preserved module in the mouse data (for example,
Zsummary = 11, P <10-18). Further, Additional file 5
demonstrates that the Magenta module eigengene in the
mouse data shows increasing expression across times; that
is, the Magenta genes are significantly (P = 0.034) overex-
pressed in weeks 16 and 20. These results indicate that the
mouse model is appropriate when it comes to an aspect of
the human disease embodied in the Magenta module.

Gene ontology and pathway enrichment analysis
To study the ontology of genes in the seven human pSS
related co-expression modules, we used the Database for
Annotation, Visualization and Integrated Discovery and
the Ingenuity Pathway Analysis software. The enrichment
results are reported in Additional file 2. Here, we highlight
results for the two most significant modules: Magenta and
Turquoise.
For the Magenta module, the top-three over-repre-
sented subcategories within 'Gene Ontology Biological
Process' were immune response (115 genes out of 423
Magenta module genes in both human and mouse data-
sets, P = 1.58 x 10-52, Bonferroni-corrected P = 3.39 x
10-49), defense response (P = 2.69 x 10-27, Bonferroni-
corrected P = 5.76 x 10-24) and inflammatory response
(46 genes in both human and mouse dataset, P = 1.96 x
10-17, Bonferroni-corrected P = 4.20 x 10-14). Additional
enriched categories for the Magenta module include anti-
gen processing and presentation (P = 2.95 x 10 17), posi-
tive regulation of response to stimulus (P = 4.30 x 10 17),
cell adhesion molecules (P = 1.34 x 10-13), humoral
immune response (P = 2.38 x 10-13), response to wound-
ing (P = 5.68 x 10-13), lymphocyte activation (P = 2.93 x
10-11), and cell activation (P = 9.03 x 10-11).
Strikingly, four of the Magenta genes are implicated the
IL-4 signaling pathway (P = 0.025). Figure 4 shows an
Ingenuity Pathway Analysis network plot where genes
colored in grey are part of the Magenta module. Genes
identified include those that encode for Toll-like recep-
tors and their signal transduction molecule MyD88,
molecules that present antigen to natural killer cells (for
example, members of CD1), molecules involved in
immune cell recruitment and adherence (for example,


CCL21, CXCL10, CXCL12, CCR7 and SLAM7), and
molecules responsible for formation of germinal centers
in secondary lymphocyte organs,for example, lympho-
toxin alpha (LTA). CXCL10 (or IP10), CXCL12 (or stro-
mal cell-derived factor-1) and SLAM7 are of particularly
interest as they indicate recruitment of specific leukocyte
subsets in particular, macrophages, dendritic cells,
T and B lymphocytes, and natural killer cells. We have
independently validated these six genes using quantitative
PCR (Additional file 6).
For the Turquoise module where the genes were
underexpressed in pSS, the top-three over-represented
subcategories within 'Gene Ontology Biological Pro-
cess' were oxidation reduction (P = 7.89 x 10- Bonfer-
roni-corrected P = 3.74 x 10 -7), generation of precursor
metabolites and energy (P = 4.72 x 10-9, Bonferroni-cor-
rected p = 2.24 x 10-5) and cofactor metabolic process
(P = 1.45 x 10 6, Bonferroni-corrected P = 0.0069). The
category 'KEGG pathway' yielded the following top-
three categories: oxidative phosphorylation (P = 9.12 x
10-7, Bonferroni-corrected P = 1.77 x 10-4), Alzheimer's
disease (P = 5.66 x 10-6, Bonferroni-corrected P =
0.0011) and Huntington's disease (P = 3.69 x 10-5, Bon-
ferroni-corrected P = 0.0071).
Additional files 3 and 4 report the results from a stan-
dard differential expression analysis. These look-up
tables allow the user to identify genes strongly related to
disease status in humans and mice, respectively. For
example, Table 2 reports genes with concordant pro-
gression patterns in human and mouse disease. Posi-
tively related genes were selected by requiring a
correlation >0.5 with human progression status (defined
as 0 for controls, 1 for sicca, and 2 for pSS) and a corre-
lation >0.5 with mouse progression status (weeks). Simi-
larly, we selected the negatively correlated genes using a
correlation threshold of -0.5. The corresponding
P values are reported in Table 2. The MM column
(MM.magenta) shows that most of these genes have
concordant pattern with the Magenta module eigengene,
illustrating again that this module contains genes that are
involved in both human and mouse disease pathology.

Discussion
We have utilized a systems approach to identify the mole-
cular and cellular events underlying the pathogenic pro-
cess in parotid glands of pSS patients, and compared the
altered expressed genes and key signaling pathways with
those in the salivary glands of mice that mimic the features
of pSS. Using a robust statistical and bioinformatics tool,
WGCNA, coupled with GO and Ingenuity Pathways Ana-
lysis software, we have clustered hundreds of co-expressed
genes into 19 gene modules. Since the module definition
does not make use of GO information, the modules are
initially named by a color.


Page 8 of 12







Horvath et al. Arthritis Research & Therapy 2012, 14:R238
http://arthritis-research.com/content/14/6/R238


2000-2011 Ingenuity Systems, Inc Al rights reserved.
Figure 4 Ingenuity Pathway Analysis software network plot. Genes shaded grey were part of the magei
altered in the parotid gland of primary syndrome patients. Test statistics and P values of individual ger
file 3 (humans) and Additional file 4 (mouse).


Our analyses highlighted seven co-expressed modules
that were enriched with genes that could discriminate
pSS patients from unaffected individuals using parotid
gland tissue, a prime target in Sj6gren's syndrome
pathogenesis. Interestingly, four of the seven modules -
namely, Magenta, Brown, Light-Cyan and Grey60 were
positively correlated with the disease status, suggesting
that the genes contained in these modules are overex-
pressed in pSS. On the other hand, Turquoise, Gray and
Salmon modules were negatively correlated with the dis-
ease status, suggesting that the genes in these modules
are underexpressed in pSS.
We acknowledge that clinical samples are not gender-
matched between the Sj6gren's subjects and the non-
Sj6gren's subjects, which may be a potential source of


nta module and are
ies can be found in Additional


bias. We have found that gender has a negligible effect
on gene expression levels in the controls. None of the
genes were gender related at a false discovery rate
threshold of 0.05. Of the 5,461 genes that were differen-
tially expressed between pSS cases and controls at a
false discovery threshold of 0.05, only 140 genes showed
differential expression (P <0.05) between males and
females of the control samples. Further, we find that the
vast majority of module genes are on autosomes and
that there is no evidence that gender affects their
expression level in the control group.
When comparing human and mouse data, we found
that the Magenta module was the most highly conserved
module between these two species. Not surprisingly, this
module was enriched with genes involved in immunity


Page 9 of 12






Horvath et al. Arthritis Research & Therapy 2012, 14:R238
http://arthritis-research.com/content/14/6/R238


and inflammation the two cardinal events in pSS
pathogenesis. To our knowledge, we are the first to map
out important intersects between human and mouse
pertaining to key molecules and associated pathways in
pSS. The knowledge gained from this work will enhance
future target-based therapies for this devastating disor-
der. Our novel co-expression module-based comparison
of human and mouse models can be used to judge
whether a given mouse model mirrors the human dis-
ease at the transcriptional level. Future studies could
explore whether other mouse models allow for the pre-
servation of additional human disease-related modules.
Characterizing the molecular and cellular events dur-
ing the progression of pSS from a healthy or even a
non-pSS state remains an important challenge. Because
of the lack of specific molecular markers, it is difficult
to determine which healthy subject or sicca patient will
progress to pSS. The gain in knowledge of these events
can suggest which non-pSS individuals are at risk for
developing pSS. In addition, the knowledge base could
be used to intervene in the progression of disease by
target-based therapies. Surprisingly, such specific regi-
mens are not in place at present, still relying on a few
trials with B-cell-targeted and TNF-directed therapies
evolved from testing on other autoimmune diseases
[6,10,29-32]. Our systems-level analyses of high-
throughput gene expression data can demystify the criti-
cal molecular alterations in pSS pathogenesis.
WGCNA is a tool of systems biology analysis and has
proven to be instrumental in identifying biological path-
ways and key gene constituents in a number of diseases
[16,22,26,33,34]. Our module-based analysis not only
alleviated the multiple testing problems inherent in
microarray data analysis, but also identified biologically
plausible, pSS-related modules that are highly signifi-
cantly enriched with relevant GO categories.
The GO analyses revealed striking correlations, both
positive and negative, between gene modules and disease
status. The most positively correlated Magenta module
was significantly enriched with genes involved in antigen
processing and presentation, and immune and inflam-
matory responses. Antigen processing and presentation
are two important events in innate immunity in relation
to pSS, and our findings support these two phenomena
[31]. Interestingly, a few Magenta genes were found to
be part of the IL-4 signaling pathway, which intercon-
nected with the TNF pathway that has been implicated
in pSS disease progression. The negatively correlated
Turquoise module is enriched with genes underex-
pressed in pSS, including those involved in oxidation
reduction, generation of precursor metabolites and
energy, and cofactor metabolic process. pSS patients are
known to have reduced energy levels, while chronic fati-
gue syndrome is a common extraglandular manifestation


[35]. Our data seem in line with this clinical presenta-
tion of the disease.
A unique feature of the present study is the direct
comparison of the human gene expression data with
that of the Sj6gren's syndrome-susceptible C57BL/6.
NOD-AeclAec2 mouse model. Since only the human
disease-related Magenta module was highly significantly
preserved in the mouse data, we focused on the
Magenta module to identify disease-relevant pathways
and target genes common to both species. A majority of
the identified genes were part of the immune response,
whereas a minor fraction was part of the inflammatory
response. Further comparison between overexpressed
genes in human and increasingly expressed genes during
the disease time course in mouse led to the identifica-
tion of CD1, CCR7, CXCL10 (or IP10) CXCL12 (or
stromal cell-derived factor-1), SLAM7 and lymphotoxin
alpha (LTA). These genes are of particular interest as
they indicate recruitment of specific leukocyte subsets -
in particular, macrophages, dendritic cells, T and B lym-
phocytes, and natural killer cells or are involved in
germinal center formation [36]. These observations pro-
vide a compelling concept that our systems approach
can identify targets and pathways overlapped in human
and mouse, supporting the concept that these overlap-
ping genes and their associated pathways are critical for
pSS disease development and manifestations of clinical
disease. We have independently validated these six
genes using quantitative PCR (Additional file 6).

Conclusion
Weighted gene co-expression network analysis identified
a pSS-related co-expression module that relates to pSS
disease not only in human parotid gland but also in
mouse salivary gland gene expression data. The key
genes that are part of the human-mouse intersection
network are useful for elucidating critical pathways and
molecular alterations dysregulated in pSS pathogenesis,
for highlighting genes that could be a major focus of
rodent-based validation studies, and ultimately for devel-
oping therapeutic targets in this debilitating disease.
This systems biologic study also illustrates how com-
parative WGCNA analysis (based on module preserva-
tion statistics) can be used to assess the suitability of a
rodent model with respect to human disease-related
transcriptional changes.

Additional material

Additional file 1: a table presenting MM membership for human
data For each gene, the table reports the MM measure, which is also
known as eigengene-based connectivity Each module gives rise to its
own MM measure; for example, MM denotes the measure for the
magenta module Columns whose name starts MM report the Pearson
correlation coefficient between the gene expression value and the


Page 10 of 12








Horvath et al. Arthritis Research & Therapy 2012, 14:R238
http://arthritis-research.com/content/14/6/R238


respective module eigengene For each MM measure (a correlation
coefficient), one can also report a corresponding correlation test P value
based on the Student t test (see columns whose name starts pMM) For
example, p MM magenta reports a two-sided correlation test P value
based on the Student t distribution Column B reports the original
module assignment based on the hierarchical cluster tree but the
module membership measures were used to select genes for the
functional enrichment analysis
Additional file 2: a table presenting GO terms for genes in human
pSS-associated co-expression modules The table reports GO
categories, uncorrected test P values, and corresponding P values that
correct for multiple comparisons using the following methods'
Bonferroni, Benjamini, and the false discovery rate It also reports the
number of population hits and related count data used to calculate the
hypergeometric test P value Genes column reports the Affymetrix
probes that were hits for the corresponding GO term The magenta
module (highlighted in yellow) is of particular interest since it was
preserved in the mouse model
Additional file 3: a table presenting standard differential expression
analysis in human data. For each gene (transcript), the results of
several two-group comparison tests (carried out with the WGCNA
function standardScreeningBinaryTrait) are reported The first group
comparison test contrasts females versus males in control subjects only
(columns B through L) The second group comparison test compares
controls versus pSS (columns M through W) The third group comparison
test compares controls versus sicca (columns x through AH) The fourth
group comparison test compares pSS versus sicca (columns AI through
AS) Further, columns AT through AZ report the results of a correlation
test where each gene was correlated with an ordinal variable that
encodes disease status (0 for controls, 1 for sicca, and 2 for pSS) Any
column whose name is preceded by corPearson reports the Pearson
correlation coefficients where the binary grouping variable (for example,
female vs males) was coded as a binary numeric variable (1 for the first
group, 0 for the second group) The meaning of the first and second
groups can be learned from column tStudent, which reports the Student
ttest statistic For example, tStudent Fvs MGenderlnControls shows that
the first group corresponds to females and the second group to males
(among the control samples) Similarly, t Student Control vs pSS and t
Student Control vs Sicca show that the first group is comprised of
controls t Student pSSvs Sicca shows that the first group is comprised of
pSS patients Column meanFirstGroup reports the mean expression value
in the first group FoldChange column reports a signed fold-change
value defined by the ratio meanFirstGroup/meanSecondGroup if
meanFirstGroup >meanSecondGroup But if meanFirstGroup
meanSecondGroup/meanFirstGroup SE FirstGroup reports the standard
error in the first group AreaUnderROC reports the area under the
receiver operating characteristic curve pvalueStudent reports the Student
t test P value that corresponds to the Student t-test statistic q Student
reports the corresponding q value (local false discovery rate) calculated
with the value R package nPresentSamples reports the number of
nonmissing observations that were available
Additional file 4: a table presenting correlations of mouse
expression data with time This comma-delimited file reports the
results from a correlation test where each mouse gene (transcript) is
correlated with time (measured in weeks) Column corTime reports the
correlation coefficient, column ZCorTime reports the corresponding
Student t test statistic, and column pValueStudentCorTime reports the
corresponding two-sided Student t test P value Column
AreaUnderROCCorTime reports the area under the receiver operating
characteristic curve calculated with the function rcorrcens in the Hmisc R
package
Additional file 5: a figure showing Magenta module expression in
the mouse data For each week (x axis) the height of the bar shows the
mean of the magenta module eigengene value (+1 standard error) P
value calculated with the Kruskal-Wallis test, which is a nonparametric
group comparison test While the magenta module was defined based
on the human data, this plot shows how the corresponding module
eigengene relates to time course in the mouse data To define the


magenta module eigengene in the mouse data, human genes were
mapped to orthologous mouse genes
Additional file 6: a figure showing the quantitative real-time PCR
validation in the mouse data The results of the quantitative RT-PCR
validation analysis involving a select group of genes (CD1D1 ortholog of
CD1D, CCR7, CXCL10, CXCL12, SLAMF9, LTA) are reported Barplots in the
upper and lower panels correspond to the expression values measured
by microarrays and RT-PCR, respectively The bars are colored according
to time (weeks) To verify the selected gene expressions (n = 7), aliquots
of salivary gland RNA originally used for the microarray data were
pooled Each cDNA preparation was quantified by spectrophotometry
and PCR performed Quantifications were determined by ImageJ Relative
gene expression values yielded by the PCR arrays are compared directly
with data yielded by the Affymetrix 3' Expression Array GeneChip Mouse
Genome 430 20 arrays Pooling RNA from each time point prior to cDNA
preparation is thought to be the underlying reason for higher transcript
detection in a couple of RT-PCR reactions (for example, in CxcI12
samples)




Abbreviations
AECG American-European Consensus Group; bp base pair; EULAR European
League against Rheumatism; GO gene ontology; H & E hematoxylin and
eosin; IL interleukin; MM module membership; PCR polymerase chain
reaction; pSS primary Sjogren's syndrome; TNF tumor necrosis factor;
WGCNA weighted gene co-expression network analysis

Acknowledgements
The authors thank Prof JLN Roodenburg, Dr MAW Witjes and Dr KP
Schepman for providing parotid tissues from head and neck cancer patients
They also thank Nadine Spielmann for helping with data acquisition

Author details
'School of Dentistry, Dental Research Institute, University of California at Los
Angeles, 10833 Le Conte Avenue, 73-017 CHS, Los Angeles, CA 90095-1668,
USA Department of Oral and Maxillofacial Surgery, University Medical
Center Groningen, University of Groningen, PO Box 30001, Groningen, the
Netherlands School of Medicine, Stanford University, 300 Pasteur Drive,
R241 MC 5324, Stanford, CA 94305, USA 4School of Dentistry, University of
Minnesota, 7536 Moos HST, 515 Delaware Street SE, Minneapolis, MN 55455,
USA sDepartment of Pathology, Immunology and Laboratory Medicine,
School of Medicine, University of Florida College of Medicine, JHMHSC D6-
33D, 1600 SW Archer Road, Gainesville, FL 32610-0275, USA

Authors' contributions
SH, ANMN-H, AV, SAM, SUG, ABP and DTWW participated in the design of
the study ANMN-H, RPEP, FGMK, AV, CGMK, FKLS, HB and ABP were
involved with the data acquisition SH, CC, AV, ABP and HZ performed the
statistical analyses SH, ANMN-H, RPEP, FGMK, AV, CGMK, FKLS, SAM, SUG,
ABP and DTWW drafted the manuscript All authors read and approved the
final manuscript

Competing interests
DTWW is the co-founder of RNAmeTRIX Inc, a molecular diagnostic
company The remaining authors declare that they have no competing
interests

Received: 20 March 2012 Revised: 5 September 2012
Accepted: 7 September 2012 Published: 1 November 2012

References
1 Fox R Sjogren's syndrome. Lancet 2005, 366321-331
2 Jonsson R, Moen K, Vestrheim D, Szodoray P Current issues in Sjogren's
syndrome. Oral Dis 2002, 8 130-140
3 Delaleu N, Jonsson R, Koller MM Sjogren's syndrome. Eur J Oral Sci 2005,
113101-113
4 Meijer JM, Meiners PM, Huddleston Slater JJ, Spijkervet FK, Kallenberg CG,
Vissink A, Bootsma H Health-related quality of life, employment and


Page 11 of 12








Horvath et al. Arthritis Research & Therapy 2012, 14:R238
http://arthritis-research.com/content/14/6/R238


disability in patients with Sjogren's syndrome. (Oxford)
2009, 48'1077-1082
5 Pijpe J, Meijer JM, Bootsma H, van der Wal JE, Spijkervet FK, Kallenberg CG,
Vissink A, Ihrler S Clinical and histologic evidence of salivary gland
restoration supports the efficacy of rituximab treatment in Sjogren's
syndrome. Arthntis Rheum 2009, 603251-3256
6 Ramos-Casals M, Tzioufas AG, Stone JH, Siso A, Bosch X Treatment of
primary Sjogren syndrome: a systematic review. JAMA 2010, 304452-460
7 Vitali C, Bombardieri S, Jonsson R, Moutsopoulos HM, Alexander EL,
Carsons SE, Daniels TE, Fox PC, Fox Rl, Kassan SS, Pillemer SR, Talal N,
Weisman MH Classification criteria for Sjogren's syndrome: a revised
version of the European criteria proposed by the American-European
Consensus Group. Ann Rheum Dis 2002, 61 554-558
8 Seror R, Ravaud P, Mariette X, Bootsma H, Theander E, Hansen A, Ramos-
Casals M, Dorner T, Bombardieri S, Hachulla E, Brun JG, Kruize AA,
Praprotnik S, Tomsic M, Gottenberg JE, Devauchelle V, Devita S,
Vollenweider C, Mandl T, Tzioufas A, Carsons S, Saraux A, Sutcliffe N, Vitali C,
Bowman SJ EULAR Sjogren's Syndrome Patient Reported Index (ESSPRI):
development of a consensus patient index for primary Sjogren's
syndrome. Ann Rheum Dis 2011, 70'968-972
9 Seror R, Ravaud P, Bowman SJ, Baron G, Tzioufas A, Theander E,
Gottenberg JE, Bootsma H, Mariette X, Vitali C EULAR Sjogren's syndrome
disease activity index: development of a consensus systemic disease
activity index for primary Sjogren's syndrome. Ann Rheum Dis 2010,
69'1103-1109
10 Mariette X, Gottenberg JE Pathogenesis of Sjogren's syndrome and
therapeutic consequences. Curr Opin Rheumatol 2010, 22471-477
11 Chiorini JA, Cihakova D, Ouellette CE, Caturegli P Sjogren syndrome:
advances in the pathogenesis from animal models. J Autoimmun 2009,
33'190-196
12 Delaleu N, Nguyen CQ, Peck AB, Jonsson R Sjogren's syndrome: studying
the disease in mice. Arthritis Res Ther 2011, 13'217
13 Hu S, Zhou M, Jiang J, Wang J, Elashoff D, Gorr S, Michie SA, Spijkervet FK,
Bootsma H, Kallenberg CG, Vissink A, Horvath S, Wong DT Systems biology
analysis of Sjogren's syndrome and mucosa-associated lymphoid tissue
lymphoma in parotid glands. Arthritis Rheum 2009, 6081-92
14 Soyfoo MS, Steinfeld S, Delporte C Usefulness of mouse models to study
the pathogenesis of Sjogren's syndrome. Oral Dis 2007, 13'366-375
15 Cha S, Nagashima H, Brown VB, Peck AB, Humphreys-Beher MG Two NOD
Idd-associated intervals contribute synergistically to the development of
autoimmune exocrinopathy (Sjogren's syndrome) on a healthy murine
background. Arthritis Rheum 2002, 46'1390-1398
16 Zhang B, Horvath S A general framework for weighted gene co-
expression network analysis. Stat Appl Genet Mol Biol 2005, 4Article 7
17 Davidoff AM, Pence JC, Shorter NA, Iglehart JD, Marks JR Expression of p53
gene in human neuroblastoma and neuroepithelioma-derived cell lines.
Oncogene 1992, 7'127-133
18 Pijpe J, Kalk WW, van der Wal JE, Vissink A, Kluin PM, Roodenburg JL,
Bootsma H, Kallenberg CG, Spijkervet FK Parotid gland biopsy compared
with labial biopsy in the diagnosis of patients with primary Sjogren's
syndrome. (Oxford) 2007, 46335-341
19 R project. i I ... i /]
20 Johnson WE, Li C, Rabinovic A Adjusting batch effects in microarray
expression data using empirical Bayes methods. Biostatistics 2007,
8118-127
21 Miller JA, Cai C, Langfelder P, Geschwind DH, Kurian SM, Salomon DR,
Horvath S Strategies for aggregating gene expression data: the
collapseRows R function. BMC bioinformatics 2011, 12322
22 Langfelder P, Horvath S WGCNA: an R package for weighted correlation
network analysis. BMC Bioinformatics 2008, 9559
23 Mason MJ, Fan G, Plath K, Zhou Q, Horvath S Signed weighted gene co-
expression network analysis of transcriptional regulation in murine
embryonic stem cells. BMC Genomics 2009, 10327
24 Yip AM, Horvath S Gene network interconnectedness and the
generalized topological overlap measure. BMC Bioinformatics 2007, 8'22
25 Langfelder P, Zhang B, Horvath S Defining clusters from a hierarchical
cluster tree: the Dynamic Tree Cut package for R. Bioinformatics 2008,
24719-720
26 Horvath S, Dong J Geometric interpretation of gene co-expression
network analysis. PLoS Comput Biol 2008, 4'e000117
27 DAVID Bioinformatics Resources 6.7. [http//david abcc ncifcrfgov/]


28 Ingenuity Systems. [http//www ingenuity com]
29 Meiners PM, Vissink A, Kallenberg CG, Kroese FG, Bootsma H Treatment of
primary Sjogren's syndrome with anti-CD20 therapy (rituximab). A
feasible approach or just a starting point? Expert Opin Biol Ther 2011,
11 1381-1394
30 Kallenberg CG, Vissink A, Kroese FG, Abdulahad WH, Bootsma H What have
we learned from clinical trials in primary Sjogren's syndrome about
pathogenesis? Arthntis Res Ther 2011, 13'205
31 Hansen A, Lipsky PE, Dorner T Immunopathogenesis of primary Sjogren's
syndrome: implications for disease management and therapy. Curr Opin
Rheumatol 2005, 17'558-565
32 Rotondi M, Chiovato L, Romagnani S, Serio M, Romagnani P Role of
chemokines in endocrine autoimmune diseases. Endocrne Rev 2007,
28492-520
33 Horvath S, Zhang B, Carlson M, Lu KV, Zhu S, Felciano RM, Laurance MF,
Zhao W, Qi S, Chen Z, Lee Y, Scheck AC, Liau LM, Wu H, Geschwind DH,
Febbo PG, Kornblum HI, Cloughesy TF, Nelson SF, Mischel PS Analysis of
oncogenic signaling networks in glioblastoma identifies ASPM as a
molecular target. Proc Nati Acad Sci USA 2006, 103'17402-17407
34 Langfelder P, Luo R, Oldham MC, Horvath S Is my network module
preserved and reproducible? PLoS Comput Biol 2011, 7e1001057
35 Segal B, Thomas W, Rogers T, Leon JM, Hughes P, Patel D, Patel K,
Novitzke J, Rohrer M, Gopalakrishnan R, Myers S, Nazmul-Hossain A,
Emamian E, Huang A, Rhodus N, Moser K Prevalence, severity, and
predictors of fatigue in subjects with primary Sjogren's syndrome.
Arthritis Rheum 2008, 59'1780-1787
36 Peck AB, Saylor BT, Nguyen L, Sharma A, She JX, Nguyen CQ, Mclndoe RA
Gene expression profiling of early-phase Sjogren's syndrome in C57BL/6.
NOD-Aec1Aec2 mice identifies focal adhesion maturation associated
with infiltrating leukocytes. Invest Ophthalmol Vis So 2011, 525647-5655

doi:10.1186/ar4081
Cite this article as: Horvath et al Systems analysis of primary Sjogren's
syndrome pathogenesis in salivary glands identifies shared pathways in
human and a mouse model. Arthritis Research & Therapy 2012 14 R238


O BloMed Central


Page 12 of 12


Submit your next manuscript to BioMed Central
and take full advantage of:

* Convenient online submission
* Thorough peer review
* No space constraints or color figure charges
* Immediate publication on acceptance
* Inclusion in PubMed, CAS, Scopus and Google Scholar
* Research which is freely available for redistribution


Submit your manuscript at
www.biomedcentral.com/submit




Full Text

PAGE 1

Additional File 6 : Verification of selected gene expression Previous verifications of microarray gene expression data using randomly selected genes have been carried out using semi quantitative polymerase chain reaction (RT PCR) analysis [ 34 35 ] quantitative real time PCR (rt PCR) analysis [ 35 36 ] and comparative microarray analyses using SABioscie nce arrays [ 37 ] These analyses have provided excellent correlations. To verify a set a selected gene expressions (n=7), a liquots of salivary gland RNA originally used for the microarray data were pooled, then prepared as described elsewhere [ 36 ] Each cDNA preparation was quantified by spect rophotometry and PCR performed. Quantifications were determined by ImageJ. Relative gene expression values yielded by the PCR arrays are compared directly with data yielde Mouse Genome 430 2.0 arrays. Pooling RNA from each time point prior to cDNA preparation is thought to be the underlying reason for higher transcript detection in a couple of RT PCR rea ctions, e.g., in Cxcl1 2 samples.



PAGE 1

4 8 12 16 20 Magenta eigengene expression in mice p= 0.034 Week mean ME.magenta -0.10 0.00 0.10


xml version 1.0 encoding UTF-8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchema-instance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd
INGEST IEID EBZWQCJWH_39P2F0 INGEST_TIME 2013-03-05T20:11:31Z PACKAGE AA00013478_00001
AGREEMENT_INFO ACCOUNT UF PROJECT UFDC
FILES



PAGE 1

Systems analysis of primary Sjgren's syndrome pathogenesis in salivary glands identifies shared pathways in human and a mouse modelHorvath et al Horvath etal ArthritisResearch&Therapy 2012, 14 :R238 http://arthritis-research.com/content/14/6/R238(1November2012)

PAGE 2

RESEARCHARTICLE OpenAccessSystemsanalysisofprimarySjgren ssyndrome pathogenesisinsalivaryglandsidentifiesshared pathwaysinhumanandamousemodelSteveHorvath1,AbuNMNazmul-Hossain1,RodneyPEPollard2,FransGMKroese2,ArjanVissink2, CeesGMKallenberg2,FredKLSpijkervet2,HendrikaBootsma2,SaraAMichie3,SvenUGorr4,AmmonBPeck5, ChaochaoCai1,HuiZhou1andDavidTWWong1*AbstractIntroduction: PrimarySjgren ssyndrome(pSS)isachronicautoimmunediseasewithcomplexetiopathogenesis. Despiteextensivestudiestounderstandthediseaseprocessutilizinghumanandmousemodels,theintersection betweenthesespeciesremainselusive.Toaddressthisgap,weutilizedanovelsystemsbiologyapproachto identifydisease-relatedgenemodulesandsignalingpathwaysthatoverlapbetweenhumansandmice. Methods: Parotidglandtissueswereharvestedfrom24pSSand16non-pSSsiccapatientsand25controls.For mousestudies,salivaryglandswereharvestedfromC57BL/6.NODAec1Aec2 miceatvarioustimesduring developmentofpSS-likedisease.RNAwasanalyzedwithAffymetrixHGU133+2.0arraysforhumansamplesand withMOE430+2.0arraysformousesamples.TheimageswereprocessedwithAffymetrixsoftware.Weighted-gene co-expressionnetworkanalysiswasusedtoidentifydisease-relatedandfunctionalpathways. Results: Nineteenco-expressionmoduleswereidentifiedinhumanparotidtissue,ofwhichfourweresignificantly upregulatedandthreeweredownregulatedinpSSpatientscomparedwithnon-pSSsiccapatientsandcontrols. Notably,oneofthehumandisease-relatedmoduleswashighlypreservedinthemousemodel,andwasenriched withgenesinvolvedinimmuneandinflammatoryresponses.Furthercomparisonbetweenthesetwospeciesled totheidentificationofgenesassociatedwithleukocyterecruitmentandgerminalcenterformation. Conclusion: Oursystemsbiologyanalysisofgenome-wideexpressiondatafromsalivaryglandtissueofpSS patientsandfromapSSmousemodelidentifiedcommondysregulatedbiologicalpathwaysandmoleculartargets underlyingcriticalmolecularalterationsinpSSpathogenesis.IntroductionSjgren ssyndromeisachronic,inflammatoryautoimmunediseasecharacterizedbylymphocyticinfiltration oftheexocrineglands,especiallythesalivaryandlacrimalglands,leadingtodestructionoftheirfunctional components.Thediseasemayexistaloneasprimary Sjgren ssyndrome(pSS)orinconjunctionwith anotherautoimmunedisor derassecondarySjgren s syndrome[1].Thedisease affects0.5to1.0%ofthe generalpopulationandshowsastriking9:1femalepredominance[2,3].Althoughdrymouthanddryeyesare thehallmarksymptomsofpSS,thediseasecanaffect almostanyorganofthebodyandcancausesubstantial morbidity[4].Symptomatictherapycurrentlydominates theoverallmanagementofaffectedindividuals[5,6]. Thediagnosticcriteria,includingtheAmerican-EuropeanConsensusGroup(AECG)criteriainvolving eitherserologyorhistopat hologyinconjunctionwith oralandoculardryness[7],areofteninconsistently applied,andarenotoptimalindifferentiatingpSSfrom non-pSSsiccapatientsoraddressingextraglandular manifestations.Asaresult,diagnosisofSjgren ssyndromeoftenlagsdiseaseonsetby6to10years.Recent *Correspondence:dtww@ucla.edu Contributedequally1SchoolofDentistry,DentalResearchInstitute,UniversityofCaliforniaatLos Angeles,10833LeConteAvenue,73-017CHS,LosAngeles,CA90095-1668, USA FulllistofauthorinformationisavailableattheendofthearticleHorvath etal ArthritisResearch&Therapy 2012, 14 :R238 http://arthritis-research.com/content/14/6/R238 2012Horvathetal.;licenseeBioMedCentralLtd.ThisisanopenaccessarticledistributedunderthetermsoftheCreativeCommons AttributionLicense(http://creativecommons.org/licenses/by/2.0),whichpermitsunrestricteduse,distribution,andreproductionin anymedium,providedtheoriginalworkisproperlycited.

PAGE 3

attemptstoaddresstheseissues-theEULARSjgren s SyndromePatientReportedIndex,andadiseaseactivity indextoevaluatesystemiccomplications(EULARSjgren sSyndromeDiseaseActivityIndex)[8,9]-require furtherlarge-scalevalidation.Inadditiontothedifficulty indiseasediagnosis,theunderlyingpathophysiologic mechanismsofSjgren ssyndromeremainobscure. Acommonmodelsuggestsinteractionbetween geneticsusceptibilityandenvironmentalfactorssuchas viralinfectionsfortheimmunecellactivationandprotractedinflammatoryresponseresultinginglandular dysfunctionandsystemicautoimmunity[6].Inaddition, recentstudiesinhumanandanimalmodelshavehighlightedseveralcomponentsoftheinnateandadaptive immunesystemsaswellasnonimmunologicfactors [10-13].However,thereisagenerallackofexperimental andresearchintersectionsbetweenhumansandmouse modelsintermsofbiologicalpathwaysandkeymoleculartargets.Indeed,awidevarietyofanimalmodelsof Sjgren ssyndromehavebeendescribed;eachcapturing certainaspectsofthedisease[14].Themostwidely usedmousemodel,NOD,exhibitsCD4+lymphocyte infiltration,autoantibodie s,xerostomiaandafemaledominantphenotype.However,themicealsodevelop diabetes,whichcomplicatesthesearchforSjgren ssyndromespecificchangesofgeneexpression.Recently,the modifiedNODmousemodelC57BL/6.NODAec1Aec2 hasbeendescribed.Thesemicedevelopthesymptoms ofSjgren ssyndromebutnotdiabetes[15],suggesting thattheywillbeahighlyvaluablemodelforunderstandingthepathogenesispSS. Inthepresentstudy,weappliedweightedgenecoexpressionnetworkanalysis(WGCNA)[16]tocharacterizebiologicalpathwaysandmoleculartargetsassociatedwithpSSpathogenesisinbothhumanparotid tissueandtheC57BL/6.NODAec1Aec2 mousemodel. Oursystemsanalysisofgenome-wideexpressiondata fromhuman(pSS)salivaryglandtissuecomparedwith salivaryglandtissuefromamousemodelofpSSidentifiedcommonbiologicalpathwaysandmoleculartargets thatcanpivotallycontributetocriticalmolecularalterationsinpSSpathogenesis.MaterialsandmethodsHumanstudiesThestudyprotocolwasapprovedbytheInstitutional ReviewBoardoftheUniversityofCaliforniaatLos AngelesandtheUniversityM edicalCenteratGroningen.AllpatientswererecruitedfromtheDepartmentof OralandMaxillofacialSurgeryoftheUniversityMedical CenteratGroningen,wereatleast21yearsold,and providedwritteninformedconsenttoparticipateinthe study.Thestudyincluded24patientswhofulfilledthe 2002AECGcriteriaforpSS[17],16non-siccapatients withsubjectivesymptomsandobjectivesignsoforal andoculardrynessbutnotfulfillingtheAECGcriteria forpSS,and25controlpatientswithnosubjectiveor objectiveevidenceoforaloroculardryness.Thepatient characteristicsarepresentedinTable1.ParotidglandtissuespecimensUnderlocalanesthesia,anincisionalbiopsyofoneparotidglandwasperformedaspartofthediagnostic workupofthepSSandnon-pSSsiccapatientsfollowing AECGcriteria,asdescribed[7],andaccordingtothe proceduredescribedbyPijpeandcolleagues[18].Tissue washarvestedfromthedorsalcaudallobeoftheparotid glandfromthecontrolpatientsduringsurgeryfororal ororopharyngealsquamouscellcarcinoma. Afterharvest,mostofthespecimenwassnap-frozen anddeliveredtotheUniversityofCaliforniaatLos Angelesondryiceandthenstoredat-80Cforgene expressionprofiling.Aminorpartwasfixedin4%formalin,embeddedinparaffin,sectioned,andstainedwith H&Eforhistopathologicevaluation.Theevaluationof slidesfrompSSandnon-pSSsiccapatientswasperformedindependentlybytwooralpathologists,who determinedthefocusscore( 50lymphocytes/4mm2glandulartissue)andothercharacteristicssuchaslymphoepitheliallesions,germinalcenters,fibrosis,atrophy, andsoforth.ThespecimensfrompSSpatientsshowed characteristicfeatures(forexample,focusscore 1and presenceoflymphoepitheliallesions)whereasthose fromnon-pSSsiccaandcontrolsubjectsdidnotshow anysuchcharacteristics.Histopathologicevaluationof parotidtissuefromcontrolsrevealednocarcinoma.GeneexpressionprofilingTheisolationoftotalRNAfromsnap-frozenparotid glandtissues,insamplesizesrangingfrom10to30mg, wasperformedusingtheRNeasykitsuppliedwiththe PARISsystem(Ambion,Austin,TX,USA),accordingto themanufacturer sprotocol.Briefly, eachfrozentissue specimenwashomogenizedin300 lcoldcell-disruptionbufferandanequalvolumeof2denaturingsolutionwasadded.Thetotallysatewascentrifugedfor5 minutesat10,000 g toseparateoutthesuperficial aqueousphasecontainingRNAandthen1.25volumes of100%ethanolwereadded.Themixturewasloaded ontoaspincolumnfollowedbywashingtwicewithprovidedbuffer.RNAwaselutedat100Cinto50 lelution buffer.TheresultantRNAwassubjectedtoRNase-free DNasetreatmentfollowedbyethanolprecipitation,dissolvedinto15 lDNase/RNase-freewater,andquantifiedwithaNanodropspectrophotometer(Nanodrop Technology,Wilmington,DE,USA). TheRiboAmpRNAAmplificationsystem(Molecular Devices,Sunnyvale,CA,USA)wasusedtoperformone roundoflinearamplificationofparotidglandtissue mRNAs,taking2.5 gtotalRNAfromeachspecimenasHorvath etal ArthritisResearch&Therapy 2012, 14 :R238 http://arthritis-research.com/content/14/6/R238 Page2of12

PAGE 4

thetemplate.ThesynthesizedcDNAwastranscribedto cRNAandthenbiotinylatedusingtheGeneChip Expression3 -Amplificationreagents(Affymetrix,Santa Clara,CA,USA)for invitro labeling.Then15 g labeledcRNAwasfragmentedinto50to200bpfragmentsandassessedforqualit ywitha2100Bioanalyzer (AgilentTechnologies,PaloAlto,CA,USA). Thefragmentedbiotin-labeledcRNAfrompSS,nonpSSsiccaandcontrolsubjectswashybridizedovernight totheAffymetrixHGU133plus2.0arrays.Afterwashingtoremovetheunboundtranscripts,thehybridized chipswerestainedpriortoscanning.Theacquired imageswereprocessedwithAffymetrixMicroarraySuite softwareandanalyzedwiththeRsoftware(whichis freelyavailableat[19]).Mousestudies MousemodelofSjgren ssyndrome-likediseaseC57BL/6.NODAec1Aec2 micewerebredandmaintainedunderspecificpathogen-freeconditionswithin themousefacilityoftheDepartmentofPathologyatthe UniversityofFlorida,Gainesville.Thebreedinganduse oftheseanimalsforthepresentstudieswereapproved bytheUniversityofFloridaInstitutionalAnimalCare andUseCommittee.Theanimalsweremaintainedona 12-hourlight-darkscheduleandprovidedwithfoodand acidifiedwater adlibitum .GenerationofsalivaryglandtranscriptomedataSalivaryglandswerefreshlyexcisedfromeuthanized mice( n =5peragegroup)at4,8,12,16,or20weeks ofage,snap-frozeninliquidnitrogen,andstoredat-80 Cuntilallglandularsampleswereobtained.Eachsalivaryglandwascomprisedofasubmandibular,sublingual,andparotidglandminussalivarylymphnodes. TotalRNAwasisolatedfro msalivaryglandsofeach mouseusingtheRNeasyMini-Kit(Qiagen,Valencia, CA,USA),inaccordancewiththemanufacturer sprotocol.Hybridizationswerecarriedoutwitheachofthe25 individualRNAsamplesusingAffymetrixGeneChip MouseGenome430plus2.0Arrays,inaccordancewith themanufacturer sinstructions.EachGeneChipcontained45,000probesetsthatanalyzedtheexpression levelofover39,000transcriptsandvariantsfromover 34,000well-characterizedmousegenes.Statisticalprocedures MicroarraydatapreprocessingTheAffymetrixU133plus2.0microarraydatawereanalyzedwithfunctionsandpackagesofthestatisticalsoftwareR2.12.0andBioconductor2.7(whichcanbe downloadedfrom[19]).Expressionintensityvalueswere calculatedforthe65HumanmicroarrayCELfiles(24 pSS,16non-pSSsiccaand2 5controls)usingthemas5 functionoftheAffylibrary.Potentialarrayoutlierswere identifiedbystudyingtheiri nterarraycorrelation.After removingpotentialoutliersinacompletelyunbiased fashion,all65samplesremainedfordownstreamanalysis.Afterquantilenormalization,theComBatfunction wasusedtocorrectforbatcheffects[20].Tosummarize multipleprobesetspergene,weusedthedefaultsettingsofthecollapseRowsRfunction[21].Eachofthe 20,718genesonthearraywasthusrepresentedbythe probewiththehighestmeanexpressionvalue. Analogouspreprocessingstepswereappliedtothe25 mousemicroarraysamplescorrespondingtofivedifferenttimepoints(weeks4,8,12,16and20,eachtime pointhadfivemice).Werelatedmousegenestohuman datausingatableoforthologousgenes.Weightedgeneco-expressionnetworkanalysisand preservationanalysisThestatisticalanalysissoftware(WGCNARpackage) andRtutorialsforconstructingaweightedgenecoexpressionnetworkcanbefoundintheliterature [16,22].TheWGCNApackagefirstcalculatesallpairwisePearson scorrelationcoefficientsacrossallsamples. Inasignedweightednetwork,theresultingPearson s Table1CharacteristicsofprimarySjgren ssyndrome,nonprimarySjgren ssyndromesiccaandcontrolsubjectsCharacteristic pSS( n =24) Non-pSSsicca( n =16) Controls( n =25) Age(years) 49.015.8 52.815.4 59.012.0 Female:male 21:3 13:3 12:13 Ethnicity 24Caucasian 16Caucasian 25Caucasian Unstimulatedwholesaliva(ml/minute) 0.100.14 0.120.18 NA Stimulatedwholesaliva(ml/minute) 0.310.35 0.330.40 NA Schirmer stest(mm) 9.27.6 14.811.7 NA Focalscore 3.51.2 0.51.0 0 Sjgren ssyndrome-A(positive:negative) 23:1 2:14 0:25 Sjgren ssyndrome-B(positive:negative) 19:5 0:16 0:25Datapresentedasthemeanstandarddeviation.AllprimarySjgren ssyndrome(pSS)andnon-pSSsiccapatientsweresubjectedtoacompleteAmericanEuropeanConsensusGroupdiagnosticwork-up.NA,notassessed(incontrolsbeingheadandneckcancerpatientswithoutparotidinvolvement,nofunctional testswereperformedasallthesepatientswerealreadyscheduledforaneckdissectionandfunctionaltestswereconsideredtoomuchofaburdentothe patientsatthattime).Horvath etal ArthritisResearch&Therapy 2012, 14 :R238 http://arthritis-research.com/content/14/6/R238 Page3of12

PAGE 5

correlationmatrixistransformedintoanadjacency matrix[23]: a ( ij ) = | 0.5+0.5 cor ( x ( i ) x ( j )) | ) Thedefaultvalueofthepower b =12facilitatesa soft-thresholdingapproachthatpreservesthecontinuousnatureoftheco-expressionrelationships[16].Asa networkdissimilaritymeasureweused1-thetopologicaloverlapmeasureastheinputforaveragelinkage hierarchicalclustering [24].Weusedthedynamic branchcuttingmethodtodefinemodulesasbranchesof thehierarchicalclusteringtree[25].Unassignedbackgroundgenes,outsideeachofthemodules,were denotedwiththecolorgrey.Togroupthe20,718genes intomodules,weusedtheblockwiseModuleRfunction intheWGCNARpackage.ConnectivityandmodulemembershipmeasuresModulemembership(MM),a lsoknownaseigengenebasedconnectivity,isameasureofintramodularconnectivity[26].MMisdefinedas: MM ( i ) =cor ( x ( i ) ,ME ) where x(i )istheexpressionprofileof i thgeneandME istheeigengene(firstprincipalcomponent)ofthegiven module.WeusedtheMMmeasuretoselectmodule genesforageneontology(GO)enrichmentanalysis. TheMMmeasuresofthehumanmodulesarereported inAdditionalfile1.FunctionalenrichmentanalysisTheDatabaseforAnnotation,Visualization,andIntegratedDiscovery(availableat[27])andtheIngenuity PathwaysAnalysissoftware(IngenuitySystemsavailable at[28])wereusedtodeterminewhethersetsofgenes (forexample,preservedintramodularhubgenes)were significantlyenrichedwithknownGOs.TheIngenuity softwareonlyreportsuncorrected P values.ThefunctionalenrichmentresultsarereportedinAdditionalfile 2.ModulepreservationanalysisToevaluatewhichhumanmodulescouldalsobefound inthemousedata,weusedmodulepreservationstatisticsimplementedinthemodulePreservationRfunction [17].Foreachmoduleinthereferencedata(forexample,humandata),apermutationtestleadstotheZsummarystatisticinthetestdata(herethemousedataset). Zsummary>10indicatesstrongevidenceofpreservation whileZsummary<2indicatesnopreservationaccording tothepermutationtest.StandarddifferentialexpressionanalysisTofinddifferentiallyexpressedgenesbetweenpairsof groupings(involvingpSS,sicca,controls),weusedtwo RfunctionsintheWGCNApackage(standardScreeningBinaryTraitandstandardScreeningNumericTrait) thatreport P values,falsediscoveryrates( q values),fold changesandotherwidelyusedstatisticsforselecting differentialexpressedgenes.Additionalfile3reportsthe differentialexpressionresultsinhumansforstudying contrastingpSSversuscontrols,siccapatientsversus controls,andpSSversussiccapatients,andforrelating expressiondatatogenderincontrolsubjects.Further, Additionalfile3reportsthecorrelationofeachgene withanordinalmeasureofdiagnosis(0=controls,1= sicca,2=pSS),whichallowsthereadertoidentify genesthatpositivelyornegativelyincreasewiththediseaseprogressioninhumans.Wealsocorrelatedeach mousegenewithtime(weeks)tofindgenesthatare increasingordecreasingastimeprogresses.These results( P values,correlations)arereportedinAdditional file4.TheseresultswerealsousedtocreateTable2.ResultsIdentificationofgeneco-expressionmodulesinparotid glandsofpSSpatientsAsignedweightedgeneco-expressionnetworkwasconstructedbasedonthe65humanparotidglandtissues (24pSS,16non-pSSsiccaand25controls).The WGCNAmethodclusteredthe20,718humangenes into19distinctgeneco-exp ressionmodules.Sincethe moduledetectionisunbiasedanddoesnotmakeuseof GOinformation,eachofthemoduleswasinitially labeledwithauniquecolorasanidentifier(Figure1). Todefinearepresentativemoduleexpressionprofile (referredtoasthemoduleeigengene),wesummarized the(standardized)geneexpressionprofilesofthemoduleeigengene(=firstprincipalcomponent).Themodule eigengenecanbeconsideredaweightedaverageofthe modulegeneexpressionprofiles.Toidentifydiseaserelatedmodules,wecorrelateeachmoduleeigengene withdiseasestatus. Strikingly,sevenoutof19modulesshowedsignificant differentialexpressionbetweenpSSandcontrolsamples, whichreflectsthefactthatthousandsofgenesaredifferentiallyexpressedbetweenthetwogroups.Inparticular,wefoundahighlysignificantpositivecorrelation betweendiseaseandtheMagenta(comprisedof576 genes)moduleeigengene( r =0.67, P <210-7),the Brown(2,502genes)moduleeigengene( r =0.6, P <6 10-6),theLight-Cyan(349genes)moduleeigengene( r = 0.42, P <0.002),andtheGrey60(basedon216genes) moduleeigengene( r =0.47, P <810-4).Sinceapositivecorrelationindicatesthatthecorrespondingmodule genesareoverexpressedinthediseasedindividuals,the Magenta,Brown,Light-CyanandGrey60modulesare comprisedofgenesoverexpressedinpSSsamplescomparedwithbothnon-pSSsiccaandcontrols(Figure2). Onthecontrary,wefoundahighlysignificantnegative correlationbetweendiseasestatusandtheTurquoiseHorvath etal ArthritisResearch&Therapy 2012, 14 :R238 http://arthritis-research.com/content/14/6/R238 Page4of12

PAGE 6

Table2GeneswithconcordantprogressionpatternsinhumanandmousediseaseProgressionGeneSymbolcorHumanp.HumancorMousep.MouseMM.magenta +ADA0.575.510-70.62 9.210-40.70 + AIF1 0.51 1.610-50.57 3.110-30.67 + C1QB 0.54 3.210-60.65 4.610-40.71 + CASP3 0.59 1.810-70.52 8.210-30.59 + CYBB 0.63 1.610-80.56 3.410-30.82 + ENTPD1 0.51 1.210-50.62 8.510-40.66 + GZMA 0.63 1.710-80.70 1.010-40.72 + GZMK 0.61 5.810-80.55 4.110-30.63 +HLA-DQB10.536.010-60.65 5.010-40.48 + HLA-DRB1 0.53 5.910-60.74 2.610-50.89 + IFIT1 0.59 2.610-70.59 1.910-30.43 + ITGAX 0.51 1.210-50.63 6.810-40.68 + LY86 0.58 4.410-70.53 6.410-30.78 + PLEKHA2 0.63 1.510-80.52 7.110-30.87 + STAT1 0.67 7.210-100.64 5.210-40.60 + TLR7 0.63 1.610-80.63 6.610-40.61 -ALDH3A2-0.501.810-5-0.63 8.010-4-0.64 ANGPTL4 -0.50 1.910-5-0.58 2.610-3-0.40 ATP1B1 -0.56 1.010-6-0.58 2.510-3-0.62 CPT1A -0.54 4.010-6-0.76 8.710-6-0.55 FAF1 -0.54 4.110-6-0.55 4.510-3-0.58 NEO1 -0.52 9.610-6-0.53 6.910-3-0.47 PALMD -0.50 2.210-5-0.64 6.110-4-0.65 PDHA1 -0.58 3.910-7-0.63 8.010-4-0.71 -PPFIBP2-0.561.210-6-0.57 3.210-3-0.48corHumandenotesthecorrelationofthegenewithhumanprogressionstatus(definedas0forcontrols,1forsicca,and2forprimarySjgren ssyndrome). corMousedenotesthecorrelationofthegenewithmouseprogressionstatus(weeks).p.Humanandp.Mousereportthecorrelationtest P valuesinhumanand mousedata,respectively.MM.magentareportsthemodulemembershipvalueofthegene. Figure1 Hierarchicalclustertreesfordetectingmodulesbasedonweightedgeneco-expressionnetworkanalysis .Allarraysamples (primarySjgren ssyndromepatients,sicca,andcontrols)weretofindco-expressionmodules.Modulescorrespondtobranchesofthetreesand areassignedcolorcategoriesasindicatedbythecolorbandunderneatheachtree(seeMaterialsandmethods).Thetwopanelscorrespondto twoclustertreesresultingfromtheblockwiseModulesWGCNARfunction[21],whichwasusedtocircumventcomputationalchallenges.Each clustertreethuscorrespondstoablockofgenes. Horvath etal ArthritisResearch&Therapy 2012, 14 :R238 http://arthritis-research.com/content/14/6/R238 Page5of12

PAGE 7

(basedon3,981genes)moduleeigengene( r =-0.52, P <10-4),theGrey(basedon3,011genes)moduleeigengene( r =-0.41, P <0.003),andtheSalmon(basedon 446genes)moduleeigengene( r =-0.28, P <0.005).The Turquoise,GreyandSalmonmodulesarethuscomprised ofgenesthatareunderexpressedinpSSsamplescomparedwithbothnon-pSSsiccaandcontrols(Figure3). Allofthese P valuesremainhighlysignificantevenafter carryingoutthemoststringentmultiplecomparison adjustment(Bonferronicorrection)forthenumberof modules.Thisco-expressionmodule-basedanalysishasa majoradvantageoverastandarddifferentialgeneexpressionanalysessinceitonlyrelates19modulestopSSdiseasestatus,alleviatingthemultiplecomparisonproblem inherentintherawdata. Foreachgene,Additionalfile1reportsMMmeasures (alsoknownaseigengene-basedconnectivity)between genesandmodules. Figure2 IdentificationofprimarySjgren ssyndromedisease-relatedgenemodulesrevealedbyweightedgeneco-expression networkanalysis .Barplotsshowingthemeanvalueofthemoduleeigengene( y axis)versushumandiseaseprogressionstatus( x axis)where0 denoteshealthycontrols,1denotesnonprimarySjgren ssyndrome(non-pSS)siccapatients,and2denotespSSpatients.Kruskal-Wallistest P valuesabovetheplotsshowthatthemoduleeigengenesofallfourmodules-Magenta(576genes),Brown(2,502genes),Grey60(216genes), andLight-Cyan(349genes)-aresignificantlyoverexpressedinpSSpatientscomparedwithbothnon-pSSsiccapatientsandcontrols. Horvath etal ArthritisResearch&Therapy 2012, 14 :R238 http://arthritis-research.com/content/14/6/R238 Page6of12

PAGE 8

PreservationofhumanMagentamoduleinamouse Sjgren ssyndromemodelWeappliedthemodulePreservationRfunctiontoassess whethermodulesdefinedbythehumandatawerepreservedinC57BL/6.NODAec1Aec2 mice,Sjgren ssyndrome-likediseasemodeldatasets.Weselectedthe C57BL/6.NODAec1Aec2 miceforthecomparativegene expressionanalysisbecause itisaspontaneousdisease model,probablythebest-definedmousemodeltodate withrespecttoitsdiseaseprofile,andhasacomparative controlwiththesamegeneticbackgroundthatmimics virtuallyallmolecularand clinicalaspectsdefinedin Figure3 PrimarySjgren ssyndromedisease-relatedco-expressionmodulesthataredownregulatedinprimarySjgren ssyndrome patients .BarplotsanalogoustothosedescribedinFigure2.Themoduleeigengenesoftwomodules-Turquoise(3,981genes)andGrey(3,011 genes)-aresignificantlyunderexpressedinprimarySjgren ssyndrome(pSS)patientscomparedwithbothnon-pSSsiccapatientsandcontrols. NotethattheSalmonmodule(446genes)differsfromtheothertwomoduleswithrespecttothenon-pSSsiccapatients,wherethemodule eigengeneisunderexpressedcomparedwithcontrols. Horvath etal ArthritisResearch&Therapy 2012, 14 :R238 http://arthritis-research.com/content/14/6/R238 Page7of12

PAGE 9

humans.Importantlyforthisstudy,extensivetranscriptomicprofilingdataofsalivaryglandsdefiningdevelopmentofSjgren ssyndromeareavailable. Themodulepreservationstatisticsallowonetoquantify whichaspectsofwithin-moduletopologyarepreserved betweenareferencenetwork(humanparotidglands)and atestnetwork(mousesalivaryglands).Weusedtwocompositepreservationstatistics(ZsummaryandmedianRank) toassessoverallmodulepreservation.Weevaluatedthe preservationofthehumanmodulesintheentiremouse dataset(thatis,allweekstogether).Strikingly,bothstatisticsrevealedthatthehumanMagentamoduleisthemost highlypreservedmoduleinthemousedata(forexample, Zsummary=11, P <10-18).Further,Additionalfile5 demonstratesthattheMagentamoduleeigengeneinthe mousedatashowsincreasingexpressionacrosstimes;that is,theMagentagenesaresignificantly( P =0.034)overexpressedinweeks16and20.Theseresultsindicatethatthe mousemodelisappropriatewhenitcomestoanaspectof thehumandiseaseembodiedintheMagentamodule.GeneontologyandpathwayenrichmentanalysisTostudytheontologyofgenesinthesevenhumanpSS relatedco-expressionmodules,weusedtheDatabasefor Annotation,Visualizationan dIntegratedDiscoveryand theIngenuityPathwayAnalysissoftware.Theenrichment resultsarereportedinAdditionalfile2.Here,wehighlight resultsforthetwomostsignificantmodules:Magentaand Turquoise. FortheMagentamodule,thetop-threeover-representedsubcategorieswithin GeneOntology-Biological Process wereimmuneresponse(115genesoutof423 Magentamodulegenesinbothhumanandmousedatasets, P =1.580-52,Bonferroni-corrected P =3.39 10-49),defenseresponse( P =2.6910-27,Bonferronicorrected P =5.7610-24)andinflammatoryresponse (46genesinbothhumanandmousedataset, P =1.96 10-17,Bonferroni-corrected P =4.2010-14).Additional enrichedcategoriesfortheMagentamoduleincludeantigenprocessingandpresentation( P =2.950-17),positiveregulationofresponsetostimulus( P =4.300-17), celladhesionmolecules( P =1.340-13),humoral immuneresponse( P =2.3810-13),responsetowounding( P =5.680-13),lymphocyteactivation( P =2.93 10-11),andcellactivation( P =9.0310-11). Strikingly,fouroftheMagentagenesareimplicatedthe IL-4signalingpathway( P =0.025).Figure4showsan IngenuityPathwayAnalysisnetworkplotwheregenes coloredingreyarepartoftheMagentamodule.Genes identifiedincludethosethatencodeforToll-likereceptorsandtheirsignaltran sductionmoleculeMyD88, moleculesthatpresentantigentonaturalkillercells(for example,membersofCD1),moleculesinvolvedin immunecellrecruitmentan dadherence(forexample, CCL21,CXCL10,CXCL12,CCR7andSLAM7),and moleculesresponsibleforfor mationofgerminalcenters insecondarylymphocyteorg ans,forexample,lymphotoxinalpha(LTA).CXCL10(orIP10),CXCL12(orstromalcell-derivedfactor-1)andSLAM7areofparticularly interestastheyindicaterecruitmentofspecificleukocyte subsets-inparticular,macr ophages,dendriticcells, TandBlymphocytes,andnaturalkillercells.Wehave independentlyvalidatedthesesixgenesusingquantitative PCR(Additionalfile6). FortheTurquoisemodulewherethegeneswere underexpressedinpSS,thetop-threeover-represented subcategorieswithin GeneOntology-BiologicalProcess wereoxidationreduction( P =7.8910-11,Bonferroni-corrected P =3.740-7),generationofprecursor metabolitesandenergy( P =4.7210-9,Bonferroni-correctedp=2.2410-5)andcofactormetabolicprocess ( P =1.450-6,Bonferroni-corrected P =0.0069).The category KEGGpathway yieldedthefollowingtopthreecategories:oxidativephosphorylation( P =9.12 10-7,Bonferroni-corrected P =1.770-4),Alzheimer s disease( P =5.660-6,Bonferroni-corrected P = 0.0011)andHuntington sdisease( P =3.6910-5,Bonferroni-corrected P =0.0071). Additionalfiles3and4reporttheresultsfromastandarddifferentialexpressionanalysis.Theselook-up tablesallowtheusertoidentifygenesstronglyrelatedto diseasestatusinhumansandmice,respectively.For example,Table2reportsgeneswithconcordantprogressionpatternsinhumanandmousedisease.Positivelyrelatedgeneswere selectedbyrequiringa correlation>0.5withhumanprogressionstatus(defined as0forcontrols,1forsicca,and2forpSS)andacorrelation>0.5withmouseprogressionstatus(weeks).Similarly,weselectedthenegativelycorrelatedgenesusinga correlationthresholdof-0.5.Thecorresponding P valuesarereportedinTable2.TheMMcolumn (MM.magenta)showsthatmostofthesegeneshave concordantpatternwiththeMagentamoduleeigengene, illustratingagainthatthismodulecontainsgenesthatare involvedinbothhumanandmousediseasepathology.DiscussionWehaveutilizedasystemsapproachtoidentifythemolecularandcellulareventsunderlyingthepathogenicprocessinparotidglandsofpSSpatients,andcomparedthe alteredexpressedgenesandkeysignalingpathwayswith thoseinthesalivaryglandsofmicethatmimicthefeatures ofpSS.Usingarobuststatisticalandbioinformaticstool, WGCNA,coupledwithGOandIngenuityPathwaysAnalysissoftware,wehaveclusteredhundredsofco-expressed genesinto19genemodules.Sincethemoduledefinition doesnotmakeuseofGOinformation,themodulesare initiallynamedbyacolor.Horvath etal ArthritisResearch&Therapy 2012, 14 :R238 http://arthritis-research.com/content/14/6/R238 Page8of12

PAGE 10

Ouranalyseshighlightedsevenco-expressedmodules thatwereenrichedwithgenesthatcoulddiscriminate pSSpatientsfromunaffectedindividualsusingparotid glandtissue,aprimetargetinSjgren ssyndrome pathogenesis.Interestingly,fourofthesevenmodulesnamely,Magenta,Brown,Light-CyanandGrey60-were positivelycorrelatedwiththediseasestatus,suggesting thatthegenescontainedinthesemodulesareoverexpressedinpSS.Ontheotherhand,Turquoise,Grayand Salmonmoduleswerenegativelycorrelatedwiththediseasestatus,suggestingthatthegenesinthesemodules areunderexpressedinpSS. WeacknowledgethatclinicalsamplesarenotgendermatchedbetweentheSjgren ssubjectsandthenonSjgren ssubjects,whichmaybeapotentialsourceof bias.Wehavefoundthatgenderhasanegligibleeffect ongeneexpressionlevelsinthecontrols.Noneofthe genesweregenderrelatedatafalsediscoveryrate thresholdof0.05.Ofthe5,461genesthatweredifferentiallyexpressedbetweenpSScasesandcontrolsata falsediscoverythresholdof0.05,only140genesshowed differentialexpression( P <0.05)betweenmalesand femalesofthecontrolsamples.Further,wefindthatthe vastmajorityofmodulegenesareonautosomesand thatthereisnoevidencethatgenderaffectstheir expressionlevelinthecontrolgroup. Whencomparinghumanandmousedata,wefound thattheMagentamodulewasthemosthighlyconserved modulebetweenthesetwospecies.Notsurprisingly,this modulewasenrichedwithgenesinvolvedinimmunity Figure4 IngenuityPathwayAnalysissoftwarenetworkplot .Genesshadedgreywerepartofthemag entamoduleandaresignificantly alteredintheparotidglandofprimarySjgren ssyndromepatients.Teststatisticsand P valuesofindividualgenescanbefoundinAdditional file3(humans)andAdditionalfile4(mouse). Horvath etal ArthritisResearch&Therapy 2012, 14 :R238 http://arthritis-research.com/content/14/6/R238 Page9of12

PAGE 11

andinflammation-thetwocardinaleventsinpSS pathogenesis.Toourknowledge,wearethefirsttomap outimportantintersectsbetweenhumanandmouse pertainingtokeymoleculesandassociatedpathwaysin pSS.Theknowledgegainedfromthisworkwillenhance futuretarget-basedtherapiesforthisdevastatingdisorder.Ournovelco-expressionmodule-basedcomparison ofhumanandmousemodelscanbeusedtojudge whetheragivenmousemodelmirrorsthehumandiseaseatthetranscriptionallevel.Futurestudiescould explorewhetherothermousemodelsallowforthepreservationofadditionalhumandisease-relatedmodules. CharacterizingthemolecularandcellulareventsduringtheprogressionofpSSfromahealthyorevena non-pSSstateremainsanimportantchallenge.Because ofthelackofspecificmolecularmarkers,itisdifficult todeterminewhichhealthysubjectorsiccapatientwill progresstopSS.Thegaininknowledgeoftheseevents cansuggestwhichnon-pSSindividualsareatriskfor developingpSS.Inaddition,theknowledgebasecould beusedtointerveneinthe progressionofdiseaseby target-basedtherapies.Surprisingly,suchspecificregimensarenotinplaceatprese nt,stillrelyingonafew trialswithB-cell-targeted andTNF-directedtherapies evolvedfromtestingonotherautoimmunediseases [6,10,29-32].Oursystems-levelanalysesofhighthroughputgeneexpressiondatacandemystifythecriticalmolecularalterationsinpSSpathogenesis. WGCNAisatoolofsystemsbiologyanalysisandhas proventobeinstrumentalinidentifyingbiologicalpathwaysandkeygeneconstituentsinanumberofdiseases [16,22,26,33,34].Ourmodule-basedanalysisnotonly alleviatedthemultipletestingproblemsinherentin microarraydataanalysis,but alsoidentifiedbiologically plausible,pSS-relatedmodulesthatarehighlysignificantlyenrichedwithrelevantGOcategories. TheGOanalysesrevealedstrikingcorrelations,both positiveandnegative,betweengenemodulesanddisease status.ThemostpositivelycorrelatedMagentamodule wassignificantlyenrichedwithgenesinvolvedinantigen processingandpresentation,andimmuneandinflammatoryresponses.Antigenprocessingandpresentation aretwoimportanteventsininnateimmunityinrelation topSS,andourfindingssupportthesetwophenomena [31].Interestingly,afewMagentageneswerefoundto bepartoftheIL-4signalingp athway,whichinterconnectedwiththeTNFpathwaythathasbeenimplicated inpSSdiseaseprogression.Thenegativelycorrelated TurquoisemoduleisenrichedwithgenesunderexpressedinpSS,includingthoseinvolvedinoxidation reduction,generationofprecursormetabolitesand energy,andcofactormetabolicprocess.pSSpatientsare knowntohavereducedenergylevels,whilechronicfatiguesyndromeisacommonextraglandularmanifestation [35].Ourdataseeminlinewiththisclinicalpresentationofthedisease. Auniquefeatureofthepresentstudyisthedirect comparisonofthehumangeneexpressiondatawith thatoftheSjgren ssyndrome-susceptibleC57BL/6. NODAec1Aec2 mousemodel.Sinceonlythehuman disease-relatedMagentamodulewashighlysignificantly preservedinthemousedata,wefocusedonthe Magentamoduletoidentifydisease-relevantpathways andtargetgenescommontobothspecies.Amajorityof theidentifiedgeneswerepartoftheimmuneresponse, whereasaminorfractionwaspartoftheinflammatory response.Furthercomparisonbetweenoverexpressed genesinhumanandincreasinglyexpressedgenesduring thediseasetimecourseinmouseledtotheidentificationofCD1,CCR7,CXCL10(orIP10)CXCL12(or stromalcell-derivedfactor-1),SLAM7andlymphotoxin alpha(LTA).Thesegenesareofparticularinterestas theyindicaterecruitmentofspecificleukocytesubsetsinparticular,macrophages,dendriticcells,TandBlymphocytes,andnaturalkillercells-orareinvolvedin germinalcenterformation[36].Theseobservationsprovideacompellingconceptthatoursystemsapproach canidentifytargetsandpathwaysoverlappedinhuman andmouse,supportingtheconceptthattheseoverlappinggenesandtheirassociatedpathwaysarecriticalfor pSSdiseasedevelopmentandmanifestationsofclinical disease.Wehaveindependentlyvalidatedthesesix genesusingquantitativePCR(Additionalfile6).ConclusionWeightedgeneco-expressionnetworkanalysisidentified apSS-relatedco-expressionmodulethatrelatestopSS diseasenotonlyinhumanparotidglandbutalsoin mousesalivaryglandgeneexpressiondata.Thekey genesthatarepartofthehuman-mouseintersection networkareusefulforelucidatingcriticalpathwaysand molecularalterationsdysreg ulatedinpSSpathogenesis, forhighlightinggenesthatcouldbeamajorfocusof rodent-basedvalidationstudies,andultimatelyfordevelopingtherapeutictargetsinthisdebilitatingdisease. ThissystemsbiologicstudyalsoillustrateshowcomparativeWGCNAanalysis(basedonmodulepreservationstatistics)canbeusedto assessthesuitabilityofa rodentmodelwithrespecttohumandisease-related transcriptionalchanges.AdditionalmaterialAdditionalfile1:atablepresentingMMmembershipforhuman data .Foreachgene,thetablereportstheMMmeasure,whichisalso knownaseigengene-basedconnectivity.Eachmodulegivesrisetoits ownMMmeasure;forexample,MMdenotesthemeasureforthe magentamodule.ColumnswhosenamestartsMMreportthePearson correlationcoefficientbetweenthegeneexpressionvalueandtheHorvath etal ArthritisResearch&Therapy 2012, 14 :R238 http://arthritis-research.com/content/14/6/R238 Page10of12

PAGE 12

respectivemoduleeigengene.ForeachMMmeasure(acorrelation coefficient),onecanalsoreportacorrespondingcorrelationtest P value basedontheStudent t test(seecolumnswhosenamestartsp.MM).For example,p.MM.magentareportsatwo-sidedcorrelationtest Pvalue basedontheStudent t distribution.ColumnBreportstheoriginal moduleassignmentbasedonthehierarchicalclustertreebutthe modulemembershipmeasureswereusedtoselectgenesforthe functionalenrichmentanalysis. Additionalfile2:atablepresentingGOtermsforgenesinhuman pSS-associatedco-expressionmodules.ThetablereportsGO categories,uncorrectedtest P values,andcorresponding P valuesthat correctformultiplecomparisonsusingthefollowingmethods: Bonferroni,Benjamini,andthefalsediscoveryrate.Italsoreportsthe numberofpopulationhitsandrelatedcountdatausedtocalculatethe hypergeometrictest P value.GenescolumnreportstheAffymetrix probesthatwerehitsforthecorrespondingGOterm.Themagenta module(highlightedinyellow)isofparticularinterestsinceitwas preservedinthemousemodel. Additionalfile3:atablepresentingstandarddifferentialexpression analysisinhumandata.Foreachgene(transcript),theresultsof severaltwo-groupcomparisontests(carriedoutwiththeWGCNA functionstandardScreeningBinaryTrait)arereported .Thefirstgroup comparisontestcontrastsfemalesversusmalesincontrolsubjectsonly (columnsBthroughL).Thesecondgroupcomparisontestcompares controlsversuspSS(columnsMthroughW).Thethirdgroupcomparison testcomparescontrolsversussicca(columnsthroughAH).Thefourth groupcomparisontestcomparespSSversussicca(columnsAIthrough AS).Further,columnsATthroughAZreporttheresultsofacorrelation testwhereeachgenewascorrelatedwithanordinalvariablethat encodesdiseasestatus(0forcontrols,1forsicca,and2forpSS).Any columnwhosenameisprecededbycorPearsonreportsthePearson correlationcoefficientswherethebinarygroupingvariable(forexample, femalevs.males)wascodedasabinarynumericvariable(1forthefirst group,0forthesecondgroup).Themeaningofthefirstandsecond groupscanbelearnedfromcolumnt.Student,whichreportstheStudent t -teststatistic.Forexample,t.Student.F.vs.MGenderInControlsshowsthat thefirstgroupcorrespondstofemalesandthesecondgrouptomales (amongthecontrolsamples).Similarly,t.Student.Control.vs.pSSandt. Student.Control.vs.Siccashowthatthefirstgroupiscomprisedof controls.t.Student.pSS.vs.Siccashowsthatthefirstgroupiscomprisedof pSSpatients.ColumnmeanFirstGroupreportsthemeanexpressionvalue inthefirstgroup.FoldChangecolumnreportsasignedfold-change valuedefinedbytheratiomeanFirstGroup/meanSecondGroupif meanFirstGroup>meanSecondGroup.ButifmeanFirstGroup
PAGE 13

disabilityinpatientswithSjogren ssyndrome. Rheumatology(Oxford) 2009, 48 :1077-1082. 5.PijpeJ,MeijerJM,BootsmaH,vanderWalJE,SpijkervetFK,KallenbergCG, VissinkA,IhrlerS: Clinicalandhistologicevidenceofsalivarygland restorationsupportstheefficacyofrituximabtreatmentinSjogren s syndrome. ArthritisRheum 2009, 60 :3251-3256. 6.Ramos-CasalsM,TzioufasAG,StoneJH,SisoA,BoschX: Treatmentof primarySjogrensyndrome:asystematicreview. JAMA 2010, 304 :452-460. 7.VitaliC,BombardieriS,JonssonR,MoutsopoulosHM,AlexanderEL, CarsonsSE,DanielsTE,FoxPC,FoxRI,KassanSS,PillemerSR,TalalN, WeismanMH: ClassificationcriteriaforSjogren ssyndrome:arevised versionoftheEuropeancriteriaproposedbytheAmerican-European ConsensusGroup. AnnRheumDis 2002, 61 :554-558. 8.SerorR,RavaudP,MarietteX,BootsmaH,TheanderE,HansenA,RamosCasalsM,DornerT,BombardieriS,HachullaE,BrunJG,KruizeAA, PraprotnikS,TomsicM,GottenbergJE,DevauchelleV,DevitaS, VollenweiderC,MandlT,TzioufasA,CarsonsS,SarauxA,SutcliffeN,VitaliC, BowmanSJ: EULARSjogren sSyndromePatientReportedIndex(ESSPRI): developmentofaconsensuspatientindexforprimarySjogren s syndrome. AnnRheumDis 2011, 70 :968-972. 9.SerorR,RavaudP,BowmanSJ,BaronG,TzioufasA,TheanderE, GottenbergJE,BootsmaH,MarietteX,VitaliC: EULARSjogren ssyndrome diseaseactivityindex:developmentofaconsensussystemicdisease activityindexforprimarySjogren ssyndrome. AnnRheumDis 2010, 69 :1103-1109. 10.MarietteX,GottenbergJE: PathogenesisofSjogren ssyndromeand therapeuticconsequences. CurrOpinRheumatol 2010, 22 :471-477. 11.ChioriniJA,CihakovaD,OuelletteCE,CaturegliP: Sjogrensyndrome: advancesinthepathogenesisfromanimalmodels. JAutoimmun 2009, 33 :190-196. 12.DelaleuN,NguyenCQ,PeckAB,JonssonR: Sjogren ssyndrome:studying thediseaseinmice. ArthritisResTher 2011, 13 :217. 13.HuS,ZhouM,JiangJ,WangJ,ElashoffD,GorrS,MichieSA,SpijkervetFK, BootsmaH,KallenbergCG,VissinkA,HorvathS,WongDT: Systemsbiology analysisofSjogren ssyndromeandmucosa-associatedlymphoidtissue lymphomainparotidglands. ArthritisRheum 2009, 60 :81-92. 14.SoyfooMS,SteinfeldS,DelporteC: Usefulnessofmousemodelstostudy thepathogenesisofSjogren ssyndrome. OralDis 2007, 13 :366-375. 15.ChaS,NagashimaH,BrownVB,PeckAB,Humphreys-BeherMG: TwoNOD Idd-associatedintervalscontributesynergisticallytothedevelopmentof autoimmuneexocrinopathy(Sjogren ssyndrome)onahealthymurine background. ArthritisRheum 2002, 46 :1390-1398. 16.ZhangB,HorvathS: Ageneralframeworkforweightedgenecoexpressionnetworkanalysis. StatApplGenetMolBiol 2005, 4 :Article17. 17.DavidoffAM,PenceJC,ShorterNA,IglehartJD,MarksJR: Expressionofp53 geneinhumanneuroblastomaandneuroepithelioma-derivedcelllines. Oncogene 1992, 7 :127-133. 18.PijpeJ,KalkWW,vanderWalJE,VissinkA,KluinPM,RoodenburgJL, BootsmaH,KallenbergCG,SpijkervetFK: Parotidglandbiopsycompared withlabialbiopsyinthediagnosisofpatientswithprimarySjogren s syndrome. Rheumatology(Oxford) 2007, 46 :335-341. 19. Rproject. [http://cran.r-project.org/]. 20.JohnsonWE,LiC,RabinovicA: Adjustingbatcheffectsinmicroarray expressiondatausingempiricalBayesmethods. Biostatistics 2007, 8 :118-127. 21.MillerJA,CaiC,LangfelderP,GeschwindDH,KurianSM,SalomonDR, HorvathS: Strategiesforaggregatinggeneexpressiondata:the collapseRowsRfunction. BMCbioinformatics 2011, 12 :322. 22.LangfelderP,HorvathS: WGCNA:anRpackageforweightedcorrelation networkanalysis. BMCBioinformatics 2008, 9 :559. 23.MasonMJ,FanG,PlathK,ZhouQ,HorvathS: Signedweightedgenecoexpressionnetworkanalysisoftranscriptionalregulationinmurine embryonicstemcells. BMCGenomics 2009, 10 :327. 24.YipAM,HorvathS: Genenetworkinterconnectednessandthe generalizedtopologicaloverlapmeasure. BMCBioinformatics 2007, 8 :22. 25.LangfelderP,ZhangB,HorvathS: Definingclustersfromahierarchical clustertree:theDynamicTreeCutpackageforR. Bioinformatics 2008, 24 :719-720. 26.HorvathS,DongJ: Geometricinterpretationofgeneco-expression networkanalysis. PLoSComputBiol 2008, 4 :e1000117. 27. DAVIDBioinformaticsResources6.7. [http://david.abcc.ncifcrf.gov/]. 28. IngenuitySystems. [http://www.ingenuity.com]. 29.MeinersPM,VissinkA,KallenbergCG,KroeseFG,BootsmaH: Treatmentof primarySjogren ssyndromewithanti-CD20therapy(rituximab).A feasibleapproachorjustastartingpoint? ExpertOpinBiolTher 2011, 11 :1381-1394. 30.KallenbergCG,VissinkA,KroeseFG,AbdulahadWH,BootsmaH: Whathave welearnedfromclinicaltrialsinprimarySjogren ssyndromeabout pathogenesis? ArthritisResTher 2011, 13 :205. 31.HansenA,LipskyPE,DornerT: ImmunopathogenesisofprimarySjogren s syndrome:implicationsfordiseasemanagementandtherapy. CurrOpin Rheumatol 2005, 17 :558-565. 32.RotondiM,ChiovatoL,RomagnaniS,SerioM,RomagnaniP: Roleof chemokinesinendocrineautoimmunediseases. EndocrineRev 2007, 28 :492-520. 33.HorvathS,ZhangB,CarlsonM,LuKV,ZhuS,FelcianoRM,LauranceMF, ZhaoW,QiS,ChenZ,LeeY,ScheckAC,LiauLM,WuH,GeschwindDH, FebboPG,KornblumHI,CloughesyTF,NelsonSF,MischelPS: Analysisof oncogenicsignalingnetworksinglioblastomaidentifiesASPMasa moleculartarget. ProcNatlAcadSciUSA 2006, 103 :17402-17407. 34.LangfelderP,LuoR,OldhamMC,HorvathS: Ismynetworkmodule preservedandreproducible? PLoSComputBiol 2011, 7 :e1001057. 35.SegalB,ThomasW,RogersT,LeonJM,HughesP,PatelD,PatelK, NovitzkeJ,RohrerM,GopalakrishnanR,MyersS,Nazmul-HossainA, EmamianE,HuangA,RhodusN,MoserK: Prevalence,severity,and predictorsoffatigueinsubjectswithprimarySjogren ssyndrome. ArthritisRheum 2008, 59 :1780-1787. 36.PeckAB,SaylorBT,NguyenL,SharmaA,SheJX,NguyenCQ,McIndoeRA: Geneexpressionprofilingofearly-phaseSjogren ssyndromeinC57BL/6. NOD-Aec1Aec2miceidentifiesfocaladhesionmaturationassociated withinfiltratingleukocytes. InvestOphthalmolVisSci 2011, 52 :5647-5655.doi:10.1186/ar4081 Citethisarticleas: Horvath etal .: SystemsanalysisofprimarySjgren s syndromepathogenesisinsalivaryglandsidentifiessharedpathwaysin humanandamousemodel. ArthritisResearch&Therapy 2012 14 :R238. Submit your next manuscript to BioMed Central and take full advantage of: Convenient online submission Thorough peer review No space constraints or color gure charges Immediate publication on acceptance Inclusion in PubMed, CAS, Scopus and Google Scholar Research which is freely available for redistribution Submit your manuscript at www.biomedcentral.com/submit Horvath etal ArthritisResearch&Therapy 2012, 14 :R238 http://arthritis-research.com/content/14/6/R238 Page12of12


!DOCTYPE art SYSTEM 'http:www.biomedcentral.comxmlarticle.dtd'
ui ar4081
ji 1478-6354
fm
dochead Research article
bibl
title p Systems analysis of primary Sjögren's syndrome pathogenesis in salivary glands identifies shared pathways in human and a mouse model
aug
au id A1 ce yes snm Horvathfnm Steveinsr iid I1 email shorvath@mednet.ucla.edu
A2 Nazmul-Hossainmi NMAbuabu0001@ucla.edu
A3 PollardPERodneyI2 r.p.e.pollard@umcg.nl
A4 KroeseGMFransf.g.m.koese@med.umcg.nl
A5 VissinkArjana.vissink@umcg.nl
A6 KallenbergGMCeesc.g.m.kallenberg@umcg.nl
A7 SpijkervetKLFredf.k.l.spijkervet@umcg.nl
A8 BootsmaHendrikah.bootsma@umcg.nl
A9 MichieASaraI3 smichie@stanford.edu
A10 GorrUSvenI4 sugorr@umn.edu
A11 PeckBAmmonI5 peck@pathology.ufl.edu
A12 CaiChaochaocaichaochao@gmail.com
A13 ZhouHuihzhou@dentistry.ucla.edu
A14 ca WongTWDaviddtww@ucla.edu
insg
ins School of Dentistry, Dental Research Institute, University of California at Los Angeles, 10833 Le Conte Avenue, 73-017 CHS, Los Angeles, CA 90095-1668, USA
Department of Oral and Maxillofacial Surgery, University Medical Center Groningen, University of Groningen, P.O. Box 30.001, Groningen, the Netherlands
School of Medicine, Stanford University, 300 Pasteur Drive, R241 MC 5324, Stanford, CA 94305, USA
School of Dentistry, University of Minnesota, 7536 Moos HST, 515 Delaware Street SE, Minneapolis, MN 55455, USA
Department of Pathology, Immunology and Laboratory Medicine, School of Medicine, University of Florida College of Medicine, JHMHSC D6-33D, 1600 SW Archer Road, Gainesville, FL 32610-0275, USA
source Arthritis Research & Therapy
issn 1478-6354
pubdate 2012
volume 14
issue 6
fpage R238
url http://arthritis-research.com/content/14/6/R238
xrefbib pubidlist pubid idtype doi 10.1186/ar4081pmpid 23116360
history rec date day 20month 3year 2012revrec 592012acc 792012pub 1112012
cpyrt 2012collab Horvath et al.; licensee BioMed Central Ltd.note This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
abs
sec st Abstract
Introduction
Primary Sjögren's syndrome (pSS) is a chronic autoimmune disease with complex etiopathogenesis. Despite extensive studies to understand the disease process utilizing human and mouse models, the intersection between these species remains elusive. To address this gap, we utilized a novel systems biology approach to identify disease-related gene modules and signaling pathways that overlap between humans and mice.
Methods
Parotid gland tissues were harvested from 24 pSS and 16 non-pSS sicca patients and 25 controls. For mouse studies, salivary glands were harvested from C57BL/6.NOD-it Aec1Aec2 mice at various times during development of pSS-like disease. RNA was analyzed with Affymetrix HG U133+2.0 arrays for human samples and with MOE430+2.0 arrays for mouse samples. The images were processed with Affymetrix software. Weighted-gene co-expression network analysis was used to identify disease-related and functional pathways.
Results
Nineteen co-expression modules were identified in human parotid tissue, of which four were significantly upregulated and three were downregulated in pSS patients compared with non-pSS sicca patients and controls. Notably, one of the human disease-related modules was highly preserved in the mouse model, and was enriched with genes involved in immune and inflammatory responses. Further comparison between these two species led to the identification of genes associated with leukocyte recruitment and germinal center formation.
Conclusion
Our systems biology analysis of genome-wide expression data from salivary gland tissue of pSS patients and from a pSS mouse model identified common dysregulated biological pathways and molecular targets underlying critical molecular alterations in pSS pathogenesis.
bdy
Introduction
Sjögren's syndrome is a chronic, inflammatory autoimmune disease characterized by lymphocytic infiltration of the exocrine glands, especially the salivary and lacrimal glands, leading to destruction of their functional components. The disease may exist alone as primary Sjögren's syndrome (pSS) or in conjunction with another autoimmune disorder as secondary Sjögren's syndrome abbrgrp abbr bid B1 1. The disease affects 0.5 to 1.0% of the general population and shows a striking 9:1 female predominance B2 2B3 3. Although dry mouth and dry eyes are the hallmark symptoms of pSS, the disease can affect almost any organ of the body and can cause substantial morbidity B4 4. Symptomatic therapy currently dominates the overall management of affected individuals B5 5B6 6.
The diagnostic criteria, including the American-European Consensus Group (AECG) criteria involving either serology or histopathology in conjunction with oral and ocular dryness B7 7, are often inconsistently applied, and are not optimal in differentiating pSS from non-pSS sicca patients or addressing extraglandular manifestations. As a result, diagnosis of Sjögren's syndrome often lags disease onset by 6 to 10 years. Recent attempts to address these issues the EULAR Sjögren's Syndrome Patient Reported Index, and a disease activity index to evaluate systemic complications (EULAR Sjögren's Syndrome Disease Activity Index) B8 8B9 9 require further large-scale validation. In addition to the difficulty in disease diagnosis, the underlying pathophysiologic mechanisms of Sjögren's syndrome remain obscure.
A common model suggests interaction between genetic susceptibility and environmental factors such as viral infections for the immune cell activation and protracted inflammatory response resulting in glandular dysfunction and systemic autoimmunity 6. In addition, recent studies in human and animal models have highlighted several components of the innate and adaptive immune systems as well as nonimmunologic factors B10 10B11 11B12 12B13 13. However, there is a general lack of experimental and research intersections between humans and mouse models in terms of biological pathways and key molecular targets. Indeed, a wide variety of animal models of Sjögren's syndrome have been described; each capturing certain aspects of the disease B14 14. The most widely used mouse model, NOD, exhibits CD4sup + lymphocyte infiltration, autoantibodies, xerostomia and a female-dominant phenotype. However, the mice also develop diabetes, which complicates the search for Sjögren's syndrome specific changes of gene expression. Recently, the modified NOD mouse model C57BL/6.NOD-Aec1Aec2 has been described. These mice develop the symptoms of Sjögren's syndrome but not diabetes B15 15, suggesting that they will be a highly valuable model for understanding the pathogenesis pSS.
In the present study, we applied weighted gene co-expression network analysis (WGCNA) B16 16 to characterize biological pathways and molecular targets associated with pSS pathogenesis in both human parotid tissue and the C57BL/6.NOD-Aec1Aec2 mouse model. Our systems analysis of genome-wide expression data from human (pSS) salivary gland tissue compared with salivary gland tissue from a mouse model of pSS identified common biological pathways and molecular targets that can pivotally contribute to critical molecular alterations in pSS pathogenesis.
Materials and methods
Human studies
The study protocol was approved by the Institutional Review Board of the University of California at Los Angeles and the University Medical Center at Groningen. All patients were recruited from the Department of Oral and Maxillofacial Surgery of the University Medical Center at Groningen, were at least 21 years old, and provided written informed consent to participate in the study. The study included 24 patients who fulfilled the 2002 AECG criteria for pSS B17 17, 16 non-sicca patients with subjective symptoms and objective signs of oral and ocular dryness but not fulfilling the AECG criteria for pSS, and 25 control patients with no subjective or objective evidence of oral or ocular dryness. The patient characteristics are presented in Table tblr tid T1 1.
tbl hint_layout double Table 1caption Characteristics of primary Sjögren's syndrome, nonprimary Sjögren's syndrome sicca and control subjectstblbdy cols 4
r
c left
b Characteristic
pSS (n = 24)
Non-pSS sicca (n = 16)
Controls (n = 25)
cspan
hr
Age (years)
49.0 ± 15.8
52.8 ± 15.4
59.0 ± 12.0
Female:male
21:3
13:3
12:13
Ethnicity
24 Caucasian
16 Caucasian
25 Caucasian
Unstimulated whole saliva (ml/minute)
0.10 ± 0.14
0.12 ± 0.18
NA
Stimulated whole saliva (ml/minute)
0.31 ± 0.35
0.33 ± 0.40
NA
Schirmer's test (mm)
9.2 ± 7.6
14.8 ± 11.7
NA
Focal score
3.5 ± 1.2
0.5 ± 1.0
0
Sjögren's syndrome A (positive:negative)
23:1
2:14
0:25
Sjögren's syndrome B (positive:negative)
19:5
0:16
0:25
tblfn
Data presented as the mean ± standard deviation. All primary Sjögren's syndrome (pSS) and non-pSS sicca patients were subjected to a complete American-European Consensus Group diagnostic work-up. NA, not assessed (in controls being head and neck cancer patients without parotid involvement, no functional tests were performed as all these patients were already scheduled for a neck dissection and functional tests were considered too much of a burden to the patients at that time).
Parotid gland tissue specimens
Under local anesthesia, an incisional biopsy of one parotid gland was performed as part of the diagnostic workup of the pSS and non-pSS sicca patients following AECG criteria, as described 7, and according to the procedure described by Pijpe and colleagues B18 18. Tissue was harvested from the dorsal caudal lobe of the parotid gland from the control patients during surgery for oral or oropharyngeal squamous cell carcinoma.
After harvest, most of the specimen was snap-frozen and delivered to the University of California at Los Angeles on dry ice and then stored at -80°C for gene expression profiling. A minor part was fixed in 4% formalin, embedded in paraffin, sectioned, and stained with H & E for histopathologic evaluation. The evaluation of slides from pSS and non-pSS sicca patients was performed independently by two oral pathologists, who determined the focus score (≥50 lymphocytes/4 mm2 glandular tissue) and other characteristics such as lymphoepithelial lesions, germinal centers, fibrosis, atrophy, and so forth. The specimens from pSS patients showed characteristic features (for example, focus score ≥1 and presence of lymphoepithelial lesions) whereas those from non-pSS sicca and control subjects did not show any such characteristics. Histopathologic evaluation of parotid tissue from controls revealed no carcinoma.
Gene expression profiling
The isolation of total RNA from snap-frozen parotid gland tissues, in sample sizes ranging from 10 to 30 mg, was performed using the RNeasy kit supplied with the PARIS system (Ambion, Austin, TX, USA), according to the manufacturer's protocol. Briefly, each frozen tissue specimen was homogenized in 300 μl cold cell-disruption buffer and an equal volume of 2 × denaturing solution was added. The total lysate was centrifuged for 5 minutes at 10,000 × g to separate out the superficial aqueous phase containing RNA and then 1.25 volumes of 100% ethanol were added. The mixture was loaded onto a spin column followed by washing twice with provided buffer. RNA was eluted at 100°C into 50 μl elution buffer. The resultant RNA was subjected to RNase-free DNase treatment followed by ethanol precipitation, dissolved into 15 μl DNase/RNase-free water, and quantified with a Nanodrop spectrophotometer (Nanodrop Technology, Wilmington, DE, USA).
The RiboAmp RNA Amplification system (Molecular Devices, Sunnyvale, CA, USA) was used to perform one round of linear amplification of parotid gland tissue mRNAs, taking 2.5 μg total RNA from each specimen as the template. The synthesized cDNA was transcribed to cRNA and then biotinylated using the GeneChip Expression 3'-Amplification reagents (Affymetrix, Santa Clara, CA, USA) for in vitro labeling. Then 15 μg labeled cRNA was fragmented into 50 to 200 bp fragments and assessed for quality with a 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA).
The fragmented biotin-labeled cRNA from pSS, non-pSS sicca and control subjects was hybridized overnight to the Affymetrix HG U133 plus 2.0 arrays. After washing to remove the unbound transcripts, the hybridized chips were stained prior to scanning. The acquired images were processed with Affymetrix Microarray Suite software and analyzed with the R software (which is freely available at B19 19).
Mouse studies
Mouse model of Sjögren's syndrome-like disease
C57BL/6.NOD-Aec1Aec2 mice were bred and maintained under specific pathogen-free conditions within the mouse facility of the Department of Pathology at the University of Florida, Gainesville. The breeding and use of these animals for the present studies were approved by the University of Florida Institutional Animal Care and Use Committee. The animals were maintained on a 12-hour light-dark schedule and provided with food and acidified water ad libitum.
Generation of salivary gland transcriptome data
Salivary glands were freshly excised from euthanized mice (n = 5 per age group) at 4, 8, 12, 16, or 20 weeks of age, snap-frozen in liquid nitrogen, and stored at -80°C until all glandular samples were obtained. Each salivary gland was comprised of a submandibular, sublingual, and parotid gland minus salivary lymph nodes. Total RNA was isolated from salivary glands of each mouse using the RNeasy Mini-Kit (Qiagen, Valencia, CA, USA), in accordance with the manufacturer's protocol. Hybridizations were carried out with each of the 25 individual RNA samples using Affymetrix GeneChip Mouse Genome 430 plus 2.0 Arrays, in accordance with the manufacturer's instructions. Each GeneChip contained 45,000 probe sets that analyzed the expression level of over 39,000 transcripts and variants from over 34,000 well-characterized mouse genes.
Statistical procedures
Microarray data preprocessing
The Affymetrix U133 plus 2.0 microarray data were analyzed with functions and packages of the statistical software R 2.12.0 and Bioconductor 2.7 (which can be downloaded from 19). Expression intensity values were calculated for the 65 Human microarray CEL files (24 pSS, 16 non-pSS sicca and 25 controls) using the mas5 function of the Affy library. Potential array outliers were identified by studying their interarray correlation. After removing potential outliers in a completely unbiased fashion, all 65 samples remained for downstream analysis. After quantile normalization, the ComBat function was used to correct for batch effects B20 20. To summarize multiple probe sets per gene, we used the default settings of the collapseRows R function B21 21. Each of the 20,718 genes on the array was thus represented by the probe with the highest mean expression value.
Analogous preprocessing steps were applied to the 25 mouse microarray samples corresponding to five different time points (weeks 4, 8, 12, 16 and 20, each time point had five mice). We related mouse genes to human data using a table of orthologous genes.
Weighted gene co-expression network analysis and preservation analysis
The statistical analysis software (WGCNA R package) and R tutorials for constructing a weighted gene co-expression network can be found in the literature 16B22 22. The WGCNA package first calculates all pair-wise Pearson's correlation coefficients across all samples. In a signed weighted network, the resulting Pearson's correlation matrix is transformed into an adjacency matrix B23 23:
display-formula m:math name ar4081-i1 xmlns:m http:www.w3.org1998MathMathML m:mrow
m:mi a
m:mo stretchy false (
i
j
)
=
|
m:mn 0.5
+
0.5
*
m:mtext cor
(
x
(
i
)
,
x
(
j
)
)
m:msup
|
mathvariant normal β
)
The default value of the power β = 12 facilitates a soft-thresholding approach that preserves the continuous nature of the co-expression relationships 16. As a network dissimilarity measure we used 1 the topological overlap measure as the input for average linkage hierarchical clustering B24 24. We used the dynamic branch cutting method to define modules as branches of the hierarchical clustering tree B25 25. Unassigned background genes, outside each of the modules, were denoted with the color grey. To group the 20,718 genes into modules, we used the blockwiseModule R function in the WGCNA R package.
Connectivity and module membership measures
Module membership (MM), also known as eigengene-based connectivity, is a measure of intramodular connectivity B26 26. MM is defined as:
ar4081-i2
MM
(
i
)
=
cor
(
x
(
i
)
,
ME
)
where x(i) is the expression profile of ith gene and ME is the eigengene (first principal component) of the given module. We used the MM measure to select module genes for a gene ontology (GO) enrichment analysis. The MM measures of the human modules are reported in Additional file supplr sid S1 1.
suppl
Additional file 1
text a table presenting MM membership for human data. For each gene, the table reports the MM measure, which is also known as eigengene-based connectivity. Each module gives rise to its own MM measure; for example, MM denotes the measure for the magenta module. Columns whose name starts MM report the Pearson correlation coefficient between the gene expression value and the respective module eigengene. For each MM measure (a correlation coefficient), one can also report a corresponding correlation test P value based on the Student t test (see columns whose name starts p.MM). For example, p.MM.magenta reports a two-sided correlation test P value based on the Student t distribution. Column B reports the original module assignment based on the hierarchical cluster tree but the module membership measures were used to select genes for the functional enrichment analysis.
file ar4081-S1.XLS
Click here for file
Functional enrichment analysis
The Database for Annotation, Visualization, and Integrated Discovery (available at B27 27) and the Ingenuity Pathways Analysis software (Ingenuity Systems available at B28 28) were used to determine whether sets of genes (for example, preserved intramodular hub genes) were significantly enriched with known GOs. The Ingenuity software only reports uncorrected P values. The functional enrichment results are reported in Additional file S2 2.
Additional file 2
a table presenting GO terms for genes in human pSS-associated co-expression modules. The table reports GO categories, uncorrected test P values, and corresponding P values that correct for multiple comparisons using the following methods: Bonferroni, Benjamini, and the false discovery rate. It also reports the number of population hits and related count data used to calculate the hypergeometric test P value. Genes column reports the Affymetrix probes that were hits for the corresponding GO term. The magenta module (highlighted in yellow) is of particular interest since it was preserved in the mouse model.
ar4081-S2.XLS
Click here for file
Module preservation analysis
To evaluate which human modules could also be found in the mouse data, we used module preservation statistics implemented in the modulePreservation R function 17. For each module in the reference data (for example, human data), a permutation test leads to the Zsummary statistic in the test data (here the mouse dataset). Zsummary >10 indicates strong evidence of preservation while Zsummary <2 indicates no preservation according to the permutation test.
Standard differential expression analysis
To find differentially expressed genes between pairs of groupings (involving pSS, sicca, controls), we used two R functions in the WGCNA package (standardScreeningBinaryTrait and standardScreeningNumericTrait) that report P values, false discovery rates (q values), fold changes and other widely used statistics for selecting differential expressed genes. Additional file S3 3 reports the differential expression results in humans for studying contrasting pSS versus controls, sicca patients versus controls, and pSS versus sicca patients, and for relating expression data to gender in control subjects. Further, Additional file 3 reports the correlation of each gene with an ordinal measure of diagnosis (0 = controls, 1 = sicca, 2 = pSS), which allows the reader to identify genes that positively or negatively increase with the disease progression in humans. We also correlated each mouse gene with time (weeks) to find genes that are increasing or decreasing as time progresses. These results (P values, correlations) are reported in Additional file S4 4. These results were also used to create Table T2 2.
Additional file 3
a table presenting standard differential expression analysis in human data. For each gene (transcript), the results of several two-group comparison tests (carried out with the WGCNA function standardScreeningBinaryTrait) are reported. The first group comparison test contrasts females versus males in control subjects only (columns B through L). The second group comparison test compares controls versus pSS (columns M through W). The third group comparison test compares controls versus sicca (columns × through AH). The fourth group comparison test compares pSS versus sicca (columns AI through AS). Further, columns AT through AZ report the results of a correlation test where each gene was correlated with an ordinal variable that encodes disease status (0 for controls, 1 for sicca, and 2 for pSS). Any column whose name is preceded by corPearson reports the Pearson correlation coefficients where the binary grouping variable (for example, female vs. males) was coded as a binary numeric variable (1 for the first group, 0 for the second group). The meaning of the first and second groups can be learned from column t.Student, which reports the Student t-test statistic. For example, t.Student.F.vs.MGenderInControls shows that the first group corresponds to females and the second group to males (among the control samples). Similarly, t.Student.Control.vs.pSS and t.Student.Control.vs.Sicca show that the first group is comprised of controls. t.Student.pSS.vs.Sicca shows that the first group is comprised of pSS patients. Column meanFirstGroup reports the mean expression value in the first group. FoldChange column reports a signed fold-change value defined by the ratio meanFirstGroup/meanSecondGroup if meanFirstGroup >meanSecondGroup. But if meanFirstGroup ar4081-S3.CSV
Click here for file
Additional file 4
a table presenting correlations of mouse expression data with time. This comma-delimited file reports the results from a correlation test where each mouse gene (transcript) is correlated with time (measured in weeks). Column corTime reports the correlation coefficient, column ZCorTime reports the corresponding Student t-test statistic, and column pValueStudentCorTime reports the corresponding two-sided Student t test P value. Column AreaUnderROCCorTime reports the area under the receiver operating characteristic curve calculated with the function rcorr.cens in the Hmisc R package.
ar4081-S4.CSV
Click here for file
Table 2Genes with concordant progression patterns in human and mouse disease7
Progression
GeneSymbol
corHuman
p.Human
corMouse
p.Mouse
MM.magenta
+
ADA
0.57
5.5 × 10-7
0.62
9.2 × 10-4
0.70
+
AIF1
0.51
1.6 × 10-5
0.57
3.1 × 10-3
0.67
+
C1QB
0.54
3.2 × 10-6
0.65
4.6 × 10-4
0.71
+
CASP3
0.59
1.8 × 10-7
0.52
8.2 × 10-3
0.59
+
CYBB
0.63
1.6 × 10-8
0.56
3.4 × 10-3
0.82
+
ENTPD1
0.51
1.2 × 10-5
0.62
8.5 × 10-4
0.66
+
GZMA
0.63
1.7 × 10-8
0.70
1.0 × 10-4
0.72
+
GZMK
0.61
5.8 × 10-8
0.55
4.1 × 10-3
0.63
+
HLA-DQB1
0.53
6.0 × 10-6
0.65
5.0 × 10-4
0.48
+
HLA-DRB1
0.53
5.9 × 10-6
0.74
2.6 × 10-5
0.89
+
IFIT1
0.59
2.6 × 10-7
0.59
1.9 × 10-3
0.43
+
ITGAX
0.51
1.2 × 10-5
0.63
6.8 × 10-4
0.68
+
LY86
0.58
4.4 × 10-7
0.53
6.4 × 10-3
0.78
+
PLEKHA2
0.63
1.5 × 10-8
0.52
7.1 × 10-3
0.87
+
STAT1
0.67
7.2 × 10-10
0.64
5.2 × 10-4
0.60
+
TLR7
0.63
1.6 × 10-8
0.63
6.6 × 10-4
0.61
-
ALDH3A2
-0.50
1.8 × 10-5
-0.63
8.0 × 10-4
-0.64
-
ANGPTL4
-0.50
1.9 × 10-5
-0.58
2.6 × 10-3
-0.40
-
ATP1B1
-0.56
1.0 × 10-6
-0.58
2.5 × 10-3
-0.62
-
CPT1A
-0.54
4.0 × 10-6
-0.76
8.7 × 10-6
-0.55
-
FAF1
-0.54
4.1 × 10-6
-0.55
4.5 × 10-3
-0.58
-
NEO1
-0.52
9.6 × 10-6
-0.53
6.9 × 10-3
-0.47
-
PALMD
-0.50
2.2 × 10-5
-0.64
6.1 × 10-4
-0.65
-
PDHA1
-0.58
3.9 × 10-7
-0.63
8.0 × 10-4
-0.71
-
PPFIBP2
-0.56
1.2 × 10-6
-0.57
3.2 × 10-3
-0.48
corHuman denotes the correlation of the gene with human progression status (defined as 0 for controls, 1 for sicca, and 2 for primary Sjögren's syndrome). corMouse denotes the correlation of the gene with mouse progression status (weeks). p.Human and p.Mouse report the correlation test P values in human and mouse data, respectively. MM.magenta reports the module membership value of the gene.
Results
Identification of gene co-expression modules in parotid glands of pSS patients
A signed weighted gene co-expression network was constructed based on the 65 human parotid gland tissues (24 pSS, 16 non-pSS sicca and 25 controls). The WGCNA method clustered the 20,718 human genes into 19 distinct gene co-expression modules. Since the module detection is unbiased and does not make use of GO information, each of the modules was initially labeled with a unique color as an identifier (Figure figr fid F1 1). To define a representative module expression profile (referred to as the module eigengene), we summarized the (standardized) gene expression profiles of the module eigengene (= first principal component). The module eigengene can be considered a weighted average of the module gene expression profiles. To identify disease-related modules, we correlate each module eigengene with disease status.
fig Figure 1Hierarchical cluster trees for detecting modules based on weighted gene co-expression network analysis
Hierarchical cluster trees for detecting modules based on weighted gene co-expression network analysis. All array samples (primary Sjögren's syndrome patients, sicca, and controls) were to find co-expression modules. Modules correspond to branches of the trees and are assigned color categories as indicated by the color band underneath each tree (see Materials and methods). The two panels correspond to two cluster trees resulting from the blockwiseModules WGCNA R function 21, which was used to circumvent computational challenges. Each cluster tree thus corresponds to a block of genes.
graphic ar4081-1
Strikingly, seven out of 19 modules showed significant differential expression between pSS and control samples, which reflects the fact that thousands of genes are differentially expressed between the two groups. In particular, we found a highly significant positive correlation between disease and the Magenta (comprised of 576 genes) module eigengene (r = 0.67, P <2 × 10-7), the Brown (2,502 genes) module eigengene (r = 0.6, P <6 × 10-6), the Light-Cyan (349 genes) module eigengene (r = 0.42, P <0.002), and the Grey60 (based on 216 genes) module eigengene (r = 0.47, P <8 × 10-4). Since a positive correlation indicates that the corresponding module genes are overexpressed in the diseased individuals, the Magenta, Brown, Light-Cyan and Grey60 modules are comprised of genes overexpressed in pSS samples compared with both non-pSS sicca and controls (Figure F2 2). On the contrary, we found a highly significant negative correlation between disease status and the Turquoise (based on 3,981 genes) module eigengene (r = -0.52, P <10-4), the Grey (based on 3,011 genes) module eigengene (r = -0.41, P < 0.003), and the Salmon (based on 446 genes) module eigengene (r = -0.28, P <0.005). The Turquoise, Grey and Salmon modules are thus comprised of genes that are underexpressed in pSS samples compared with both non-pSS sicca and controls (Figure F3 3). All of these P values remain highly significant even after carrying out the most stringent multiple comparison adjustment (Bonferroni correction) for the number of modules. This co-expression module-based analysis has a major advantage over a standard differential gene expression analyses since it only relates 19 modules to pSS disease status, alleviating the multiple comparison problem inherent in the raw data.
Figure 2Identification of primary Sjögren's syndrome disease-related gene modules revealed by weighted gene co-expression network analysis
Identification of primary Sjögren's syndrome disease-related gene modules revealed by weighted gene co-expression network analysis. Barplots showing the mean value of the module eigengene (y axis) versus human disease progression status (x axis) where 0 denotes healthy controls, 1 denotes nonprimary Sjögren's syndrome (non-pSS) sicca patients, and 2 denotes pSS patients. Kruskal-Wallis test P values above the plots show that the module eigengenes of all four modules Magenta (576 genes), Brown (2,502 genes), Grey60 (216 genes), and Light-Cyan (349 genes) are significantly overexpressed in pSS patients compared with both non-pSS sicca patients and controls.
ar4081-2
Figure 3Primary Sjögren's syndrome disease-related co-expression modules that are downregulated in primary Sjögren's syndrome patients
Primary Sjögren's syndrome disease-related co-expression modules that are downregulated in primary Sjögren's syndrome patients. Barplots analogous to those described in Figure 2. The module eigengenes of two modules Turquoise (3,981 genes) and Grey (3,011 genes) are significantly underexpressed in primary Sjögren's syndrome (pSS) patients compared with both non-pSS sicca patients and controls. Note that the Salmon module (446 genes) differs from the other two modules with respect to the non-pSS sicca patients, where the module eigengene is underexpressed compared with controls.
ar4081-3
For each gene, Additional file 1 reports MM measures (also known as eigengene-based connectivity) between genes and modules.
Preservation of human Magenta module in a mouse Sjögren's syndrome model
We applied the modulePreservation R function to assess whether modules defined by the human data were preserved in C57BL/6.NOD-Aec1Aec2 mice, Sjögren's syndrome-like disease model datasets. We selected the C57BL/6.NOD-Aec1Aec2 mice for the comparative gene expression analysis because it is a spontaneous disease model, probably the best-defined mouse model to date with respect to its disease profile, and has a comparative control with the same genetic background that mimics virtually all molecular and clinical aspects defined in humans. Importantly for this study, extensive transcriptomic profiling data of salivary glands defining development of Sjögren's syndrome are available.
The module preservation statistics allow one to quantify which aspects of within-module topology are preserved between a reference network (human parotid glands) and a test network (mouse salivary glands). We used two composite preservation statistics (Zsummary and medianRank) to assess overall module preservation. We evaluated the preservation of the human modules in the entire mouse dataset (that is, all weeks together). Strikingly, both statistics revealed that the human Magenta module is the most highly preserved module in the mouse data (for example, Zsummary = 11, P <10-18). Further, Additional file S5 5 demonstrates that the Magenta module eigengene in the mouse data shows increasing expression across times; that is, the Magenta genes are significantly (P = 0.034) overexpressed in weeks 16 and 20. These results indicate that the mouse model is appropriate when it comes to an aspect of the human disease embodied in the Magenta module.
Additional file 5
a figure showing Magenta module expression in the mouse data. For each week (x axis) the height of the bar shows the mean of the magenta module eigengene value (±1 standard error). P value calculated with the Kruskal-Wallis test, which is a nonparametric group comparison test. While the magenta module was defined based on the human data, this plot shows how the corresponding module eigengene relates to time course in the mouse data. To define the magenta module eigengene in the mouse data, human genes were mapped to orthologous mouse genes.
ar4081-S5.PDF
Click here for file
Gene ontology and pathway enrichment analysis
To study the ontology of genes in the seven human pSS related co-expression modules, we used the Database for Annotation, Visualization and Integrated Discovery and the Ingenuity Pathway Analysis software. The enrichment results are reported in Additional file 2. Here, we highlight results for the two most significant modules: Magenta and Turquoise.
For the Magenta module, the top-three over-represented subcategories within 'Gene Ontology Biological Process' were immune response (115 genes out of 423 Magenta module genes in both human and mouse datasets, P = 1.58 × 10-52, Bonferroni-corrected P = 3.39 × 10-49), defense response (P = 2.69 × 10-27, Bonferroni-corrected P = 5.76 × 10-24) and inflammatory response (46 genes in both human and mouse dataset, P = 1.96 × 10-17, Bonferroni-corrected P = 4.20 × 10-14). Additional enriched categories for the Magenta module include antigen processing and presentation (P = 2.95 × 10-17), positive regulation of response to stimulus (P = 4.30 × 10-17), cell adhesion molecules (P = 1.34 × 10-13), humoral immune response (P = 2.38 × 10-13), response to wounding (P = 5.68 × 10-13), lymphocyte activation (P = 2.93 × 10-11), and cell activation (P = 9.03 × 10-11).
Strikingly, four of the Magenta genes are implicated the IL-4 signaling pathway (P = 0.025). Figure F4 4 shows an Ingenuity Pathway Analysis network plot where genes colored in grey are part of the Magenta module. Genes identified include those that encode for Toll-like receptors and their signal transduction molecule MyD88, molecules that present antigen to natural killer cells (for example, members of CD1), molecules involved in immune cell recruitment and adherence (for example, CCL21, CXCL10, CXCL12, CCR7 and SLAM7), and molecules responsible for formation of germinal centers in secondary lymphocyte organs,for example, lymphotoxin alpha (LTA). CXCL10 (or IP10), CXCL12 (or stromal cell-derived factor-1) and SLAM7 are of particularly interest as they indicate recruitment of specific leukocyte subsets in particular, macrophages, dendritic cells, T and B lymphocytes, and natural killer cells. We have independently validated these six genes using quantitative PCR (Additional file S6 6).
Figure 4Ingenuity Pathway Analysis software network plot
Ingenuity Pathway Analysis software network plot. Genes shaded grey were part of the magenta module and are significantly altered in the parotid gland of primary Sjögren's syndrome patients. Test statistics and P values of individual genes can be found in Additional file 3 (humans) and Additional file 4 (mouse).
ar4081-4
Additional file 6
a figure showing the quantitative real-time PCR validation in the mouse data. The results of the quantitative RT-PCR validation analysis involving a select group of genes (CD1D1 ortholog of CD1D, CCR7, CXCL10, CXCL12, SLAMF9, LTA) are reported. Barplots in the upper and lower panels correspond to the expression values measured by microarrays and RT-PCR, respectively. The bars are colored according to time (weeks). To verify the selected gene expressions (n = 7), aliquots of salivary gland RNA originally used for the microarray data were pooled. Each cDNA preparation was quantified by spectrophotometry and PCR performed. Quantifications were determined by ImageJ. Relative gene expression values yielded by the PCR arrays are compared directly with data yielded by the Affymetrix 3' Expression Array GeneChip Mouse Genome 430 2.0 arrays. Pooling RNA from each time point prior to cDNA preparation is thought to be the underlying reason for higher transcript detection in a couple of RT-PCR reactions (for example, in Cxcl12 samples).
ar4081-S6.DOCX
Click here for file
For the Turquoise module where the genes were underexpressed in pSS, the top-three over-represented subcategories within 'Gene Ontology Biological Process' were oxidation reduction (P = 7.89 × 10-11, Bonferroni-corrected P = 3.74 × 10-7), generation of precursor metabolites and energy (P = 4.72 × 10-9, Bonferroni-corrected p = 2.24 × 10-5) and cofactor metabolic process (P = 1.45 × 10-6, Bonferroni-corrected P = 0.0069). The category 'KEGG pathway' yielded the following top-three categories: oxidative phosphorylation (P = 9.12 × 10-7, Bonferroni-corrected P = 1.77 × 10-4), Alzheimer's disease (P = 5.66 × 10-6, Bonferroni-corrected P = 0.0011) and Huntington's disease (P = 3.69 × 10-5, Bonferroni-corrected P = 0.0071).
Additional files 3 and 4 report the results from a standard differential expression analysis. These look-up tables allow the user to identify genes strongly related to disease status in humans and mice, respectively. For example, Table 2 reports genes with concordant progression patterns in human and mouse disease. Positively related genes were selected by requiring a correlation >0.5 with human progression status (defined as 0 for controls, 1 for sicca, and 2 for pSS) and a correlation >0.5 with mouse progression status (weeks). Similarly, we selected the negatively correlated genes using a correlation threshold of -0.5. The corresponding P values are reported in Table 2. The MM column (MM.magenta) shows that most of these genes have concordant pattern with the Magenta module eigengene, illustrating again that this module contains genes that are involved in both human and mouse disease pathology.
Discussion
We have utilized a systems approach to identify the molecular and cellular events underlying the pathogenic process in parotid glands of pSS patients, and compared the altered expressed genes and key signaling pathways with those in the salivary glands of mice that mimic the features of pSS. Using a robust statistical and bioinformatics tool, WGCNA, coupled with GO and Ingenuity Pathways Analysis software, we have clustered hundreds of co-expressed genes into 19 gene modules. Since the module definition does not make use of GO information, the modules are initially named by a color.
Our analyses highlighted seven co-expressed modules that were enriched with genes that could discriminate pSS patients from unaffected individuals using parotid gland tissue, a prime target in Sjögren's syndrome pathogenesis. Interestingly, four of the seven modules namely, Magenta, Brown, Light-Cyan and Grey60 were positively correlated with the disease status, suggesting that the genes contained in these modules are overexpressed in pSS. On the other hand, Turquoise, Gray and Salmon modules were negatively correlated with the disease status, suggesting that the genes in these modules are underexpressed in pSS.
We acknowledge that clinical samples are not gender-matched between the Sjögren's subjects and the non-Sjögren's subjects, which may be a potential source of bias. We have found that gender has a negligible effect on gene expression levels in the controls. None of the genes were gender related at a false discovery rate threshold of 0.05. Of the 5,461 genes that were differentially expressed between pSS cases and controls at a false discovery threshold of 0.05, only 140 genes showed differential expression (P <0.05) between males and females of the control samples. Further, we find that the vast majority of module genes are on autosomes and that there is no evidence that gender affects their expression level in the control group.
When comparing human and mouse data, we found that the Magenta module was the most highly conserved module between these two species. Not surprisingly, this module was enriched with genes involved in immunity and inflammation the two cardinal events in pSS pathogenesis. To our knowledge, we are the first to map out important intersects between human and mouse pertaining to key molecules and associated pathways in pSS. The knowledge gained from this work will enhance future target-based therapies for this devastating disorder. Our novel co-expression module-based comparison of human and mouse models can be used to judge whether a given mouse model mirrors the human disease at the transcriptional level. Future studies could explore whether other mouse models allow for the preservation of additional human disease-related modules.
Characterizing the molecular and cellular events during the progression of pSS from a healthy or even a non-pSS state remains an important challenge. Because of the lack of specific molecular markers, it is difficult to determine which healthy subject or sicca patient will progress to pSS. The gain in knowledge of these events can suggest which non-pSS individuals are at risk for developing pSS. In addition, the knowledge base could be used to intervene in the progression of disease by target-based therapies. Surprisingly, such specific regimens are not in place at present, still relying on a few trials with B-cell-targeted and TNF-directed therapies evolved from testing on other autoimmune diseases 610B29 29B30 30B31 31B32 32. Our systems-level analyses of high-throughput gene expression data can demystify the critical molecular alterations in pSS pathogenesis.
WGCNA is a tool of systems biology analysis and has proven to be instrumental in identifying biological pathways and key gene constituents in a number of diseases 162226B33 33B34 34. Our module-based analysis not only alleviated the multiple testing problems inherent in microarray data analysis, but also identified biologically plausible, pSS-related modules that are highly significantly enriched with relevant GO categories.
The GO analyses revealed striking correlations, both positive and negative, between gene modules and disease status. The most positively correlated Magenta module was significantly enriched with genes involved in antigen processing and presentation, and immune and inflammatory responses. Antigen processing and presentation are two important events in innate immunity in relation to pSS, and our findings support these two phenomena 31. Interestingly, a few Magenta genes were found to be part of the IL-4 signaling pathway, which interconnected with the TNF pathway that has been implicated in pSS disease progression. The negatively correlated Turquoise module is enriched with genes underexpressed in pSS, including those involved in oxidation reduction, generation of precursor metabolites and energy, and cofactor metabolic process. pSS patients are known to have reduced energy levels, while chronic fatigue syndrome is a common extraglandular manifestation B35 35. Our data seem in line with this clinical presentation of the disease.
A unique feature of the present study is the direct comparison of the human gene expression data with that of the Sjögren's syndrome-susceptible C57BL/6.NOD-Aec1Aec2 mouse model. Since only the human disease-related Magenta module was highly significantly preserved in the mouse data, we focused on the Magenta module to identify disease-relevant pathways and target genes common to both species. A majority of the identified genes were part of the immune response, whereas a minor fraction was part of the inflammatory response. Further comparison between overexpressed genes in human and increasingly expressed genes during the disease time course in mouse led to the identification of CD1, CCR7, CXCL10 (or IP10) CXCL12 (or stromal cell-derived factor-1), SLAM7 and lymphotoxin alpha (LTA). These genes are of particular interest as they indicate recruitment of specific leukocyte subsets in particular, macrophages, dendritic cells, T and B lymphocytes, and natural killer cells or are involved in germinal center formation B36 36. These observations provide a compelling concept that our systems approach can identify targets and pathways overlapped in human and mouse, supporting the concept that these overlapping genes and their associated pathways are critical for pSS disease development and manifestations of clinical disease. We have independently validated these six genes using quantitative PCR (Additional file 6).
Conclusion
Weighted gene co-expression network analysis identified a pSS-related co-expression module that relates to pSS disease not only in human parotid gland but also in mouse salivary gland gene expression data. The key genes that are part of the human-mouse intersection network are useful for elucidating critical pathways and molecular alterations dysregulated in pSS pathogenesis, for highlighting genes that could be a major focus of rodent-based validation studies, and ultimately for developing therapeutic targets in this debilitating disease. This systems biologic study also illustrates how comparative WGCNA analysis (based on module preservation statistics) can be used to assess the suitability of a rodent model with respect to human disease-related transcriptional changes.
Abbreviations
AECG: American-European Consensus Group; bp: base pair; EULAR: European League against Rheumatism; GO: gene ontology; H & E: hematoxylin and eosin; IL: interleukin; MM: module membership; PCR: polymerase chain reaction; pSS: primary Sjögren's syndrome; TNF: tumor necrosis factor; WGCNA: weighted gene co-expression network analysis.
Competing interests
DTWW is the co-founder of RNAmeTRIX Inc., a molecular diagnostic company. The remaining authors declare that they have no competing interests.
Authors' contributions
SH, ANMN-H, AV, SAM, SUG, ABP and DTWW participated in the design of the study. ANMN-H, RPEP, FGMK, AV, CGMK, FKLS, HB and ABP were involved with the data acquisition. SH, CC, AV, ABP and HZ performed the statistical analyses. SH, ANMN-H, RPEP, FGMK, AV, CGMK, FKLS, SAM, SUG, ABP and DTWW drafted the manuscript. All authors read and approved the final manuscript.
bm
ack
Acknowledgements
The authors thank Prof. JLN Roodenburg, Dr MAW Witjes and Dr KP Schepman for providing parotid tissues from head and neck cancer patients. They also thank Nadine Spielmann for helping with data acquisition.
refgrp Sjogren's syndromeFoxRLancet2005366321lpage 33110.1016/S0140-6736(05)66990-5link fulltext 16039337Current issues in Sjogren's syndromeJonssonRMoenKVestrheimDSzodorayPOral Dis2002813014010.1034/j.1601-0825.2002.02846.x12108757Sjogren's syndromeDelaleuNJonssonRKollerMMEur J Oral Sci200511310111310.1111/j.1600-0722.2004.00183.x15819815Health-related quality of life, employment and disability in patients with Sjogren's syndromeMeijerJMMeinersPMHuddleston SlaterJJSpijkervetFKKallenbergCGVissinkABootsmaHRheumatology (Oxford)2009481077108210.1093/rheumatology/kep141Clinical and histologic evidence of salivary gland restoration supports the efficacy of rituximab treatment in Sjogren's syndromePijpeJMeijerJMBootsmaHvan der WalJESpijkervetFKKallenbergCGVissinkAIhrlerSArthritis Rheum2009603251325610.1002/art.2490319877054Treatment of primary Sjogren syndrome: a systematic reviewRamos-CasalsMTzioufasAGStoneJHSisoABoschXJAMA201030445246010.1001/jama.2010.101420664046Classification criteria for Sjogren's syndrome: a revised version of the European criteria proposed by the American-European Consensus GroupVitaliCBombardieriSJonssonRMoutsopoulosHMAlexanderELCarsonsSEDanielsTEFoxPCFoxRIKassanSSPillemerSRTalalNWeismanMHAnn Rheum Dis20026155455810.1136/ard.61.6.554pmcid 175413712006334EULAR Sjogren's Syndrome Patient Reported Index (ESSPRI): development of a consensus patient index for primary Sjogren's syndromeSerorRRavaudPMarietteXBootsmaHTheanderEHansenARamos-CasalsMDornerTBombardieriSHachullaEBrunJGKruizeAAPraprotnikSTomsicMGottenbergJEDevauchelleVDevitaSVollenweiderCMandlTTzioufasACarsonsSSarauxASutcliffeNVitaliCBowmanSJAnn Rheum Dis20117096897210.1136/ard.2010.14374321345815EULAR Sjogren's syndrome disease activity index: development of a consensus systemic disease activity index for primary Sjogren's syndromeSerorRRavaudPBowmanSJBaronGTzioufasATheanderEGottenbergJEBootsmaHMarietteXVitaliCAnn Rheum Dis2010691103110910.1136/ard.2009.110619293702219561361Pathogenesis of Sjogren's syndrome and therapeutic consequencesMarietteXGottenbergJECurr Opin Rheumatol20102247147710.1097/BOR.0b013e32833c36c520671520Sjogren syndrome: advances in the pathogenesis from animal modelsChioriniJACihakovaDOuelletteCECaturegliPJ Autoimmun20093319019610.1016/j.jaut.2009.09.009343915419800762Sjogren's syndrome: studying the disease in miceDelaleuNNguyenCQPeckABJonssonRArthritis Res Ther20111321710.1186/ar3313321887121672284Systems biology analysis of Sjogren's syndrome and mucosa-associated lymphoid tissue lymphoma in parotid glandsHuSZhouMJiangJWangJElashoffDGorrSMichieSASpijkervetFKBootsmaHKallenbergCGVissinkAHorvathSWongDTArthritis Rheum200960819210.1002/art.24150271869019116902Usefulness of mouse models to study the pathogenesis of Sjogren's syndromeSoyfooMSSteinfeldSDelporteCOral Dis20071336637510.1111/j.1601-0825.2007.01376.x17577322Two NOD Idd-associated intervals contribute synergistically to the development of autoimmune exocrinopathy (Sjogren's syndrome) on a healthy murine backgroundChaSNagashimaHBrownVBPeckABHumphreys-BeherMGArthritis Rheum2002461390139810.1002/art.1025812115247A general framework for weighted gene co-expression network analysisZhangBHorvathSStat Appl Genet Mol Biol20054Article1716646834Expression of p53 gene in human neuroblastoma and neuroepithelioma-derived cell linesDavidoffAMPenceJCShorterNAIglehartJDMarksJROncogene199271271331741160Parotid gland biopsy compared with labial biopsy in the diagnosis of patients with primary Sjogren's syndromePijpeJKalkWWvan der WalJEVissinkAKluinPMRoodenburgJLBootsmaHKallenbergCGSpijkervetFKRheumatology (Oxford)200746335341R projecthttp://cran.r-project.org/Adjusting batch effects in microarray expression data using empirical Bayes methodsJohnsonWELiCRabinovicABiostatistics2007811812710.1093/biostatistics/kxj03716632515Strategies for aggregating gene expression data: the collapseRows R functionMillerJACaiCLangfelderPGeschwindDHKurianSMSalomonDRHorvathSBMC bioinformatics20111232210.1186/1471-2105-12-322316694221816037WGCNA: an R package for weighted correlation network analysisLangfelderPHorvathSBMC Bioinformatics2008955910.1186/1471-2105-9-559263148819114008Signed weighted gene co-expression network analysis of transcriptional regulation in murine embryonic stem cellsMasonMJFanGPlathKZhouQHorvathSBMC Genomics20091032710.1186/1471-2164-10-327272753919619308Gene network interconnectedness and the generalized topological overlap measureYipAMHorvathSBMC Bioinformatics200782210.1186/1471-2105-8-22179705517250769Defining clusters from a hierarchical cluster tree: the Dynamic Tree Cut package for RLangfelderPZhangBHorvathSBioinformatics20082471972010.1093/bioinformatics/btm56318024473Geometric interpretation of gene co-expression network analysisHorvathSDongJPLoS Comput Biol20084e100011710.1371/journal.pcbi.1000117244643818704157DAVID Bioinformatics Resources 6.7http://david.abcc.ncifcrf.gov/Ingenuity Systemshttp://www.ingenuity.comTreatment of primary Sjogren's syndrome with anti-CD20 therapy (rituximab). A feasible approach or just a starting point?MeinersPMVissinkAKallenbergCGKroeseFGBootsmaHExpert Opin Biol Ther2011111381139410.1517/14712598.2011.60535221819314What have we learned from clinical trials in primary Sjogren's syndrome about pathogenesis?KallenbergCGVissinkAKroeseFGAbdulahadWHBootsmaHArthritis Res Ther20111320510.1186/ar3234315764021371351Immunopathogenesis of primary Sjogren's syndrome: implications for disease management and therapyHansenALipskyPEDornerTCurr Opin Rheumatol20051755856510.1097/01.bor.0000172801.56744.c316093833Role of chemokines in endocrine autoimmune diseasesRotondiMChiovatoLRomagnaniSSerioMRomagnaniPEndocrine Rev20072849252010.1210/er.2006-0044Analysis of oncogenic signaling networks in glioblastoma identifies ASPM as a molecular targetHorvathSZhangBCarlsonMLuKVZhuSFelcianoRMLauranceMFZhaoWQiSChenZLeeYScheckACLiauLMWuHGeschwindDHFebboPGKornblumHICloughesyTFNelsonSFMischelPSProc Natl Acad Sci USA2006103174021740710.1073/pnas.0608396103163502417090670Is my network module preserved and reproducible?LangfelderPLuoROldhamMCHorvathSPLoS Comput Biol20117e100105710.1371/journal.pcbi.1001057302425521283776Prevalence, severity, and predictors of fatigue in subjects with primary Sjogren's syndromeSegalBThomasWRogersTLeonJMHughesPPatelDPatelKNovitzkeJRohrerMGopalakrishnanRMyersSNazmul-HossainAEmamianEHuangARhodusNMoserKArthritis Rheum2008591780178710.1002/art.24311310697819035421Gene expression profiling of early-phase Sjogren's syndrome in C57BL/6.NOD-Aec1Aec2 mice identifies focal adhesion maturation associated with infiltrating leukocytesPeckABSaylorBTNguyenLSharmaASheJXNguyenCQMcIndoeRAInvest Ophthalmol Vis Sci2011525647565510.1167/iovs.11-7652317605421666236


xml version 1.0 encoding utf-8 standalone no
mets ID sort-mets_mets OBJID sword-mets LABEL DSpace SWORD Item PROFILE METS SIP Profile xmlns http:www.loc.govMETS
xmlns:xlink http:www.w3.org1999xlink xmlns:xsi http:www.w3.org2001XMLSchema-instance
xsi:schemaLocation http:www.loc.govstandardsmetsmets.xsd
metsHdr CREATEDATE 2012-12-17T20:08:14
agent ROLE CUSTODIAN TYPE ORGANIZATION
name BioMed Central
dmdSec sword-mets-dmd-1 GROUPID sword-mets-dmd-1_group-1
mdWrap SWAP Metadata MDTYPE OTHER OTHERMDTYPE EPDCX MIMETYPE textxml
xmlData
epdcx:descriptionSet xmlns:epdcx http:purl.orgeprintepdcx2006-11-16 xmlns:MIOJAVI
http:purl.orgeprintepdcxxsd2006-11-16epdcx.xsd
epdcx:description epdcx:resourceId sword-mets-epdcx-1
epdcx:statement epdcx:propertyURI http:purl.orgdcelements1.1type epdcx:valueURI http:purl.orgeprintentityTypeScholarlyWork
http:purl.orgdcelements1.1title
epdcx:valueString Systems analysis of primary Sjogren's syndrome pathogenesis in salivary glands identifies shared pathways in human and a mouse model
http:purl.orgdctermsabstract
Abstract
Introduction
Primary Sjögren's syndrome (pSS) is a chronic autoimmune disease with complex etiopathogenesis. Despite extensive studies to understand the disease process utilizing human and mouse models, the intersection between these species remains elusive. To address this gap, we utilized a novel systems biology approach to identify disease-related gene modules and signaling pathways that overlap between humans and mice.
Methods
Parotid gland tissues were harvested from 24 pSS and 16 non-pSS sicca patients and 25 controls. For mouse studies, salivary glands were harvested from C57BL/6.NOD-Aec1Aec2 mice at various times during development of pSS-like disease. RNA was analyzed with Affymetrix HG U133+2.0 arrays for human samples and with MOE430+2.0 arrays for mouse samples. The images were processed with Affymetrix software. Weighted-gene co-expression network analysis was used to identify disease-related and functional pathways.
Results
Nineteen co-expression modules were identified in human parotid tissue, of which four were significantly upregulated and three were downregulated in pSS patients compared with non-pSS sicca patients and controls. Notably, one of the human disease-related modules was highly preserved in the mouse model, and was enriched with genes involved in immune and inflammatory responses. Further comparison between these two species led to the identification of genes associated with leukocyte recruitment and germinal center formation.
Conclusion
Our systems biology analysis of genome-wide expression data from salivary gland tissue of pSS patients and from a pSS mouse model identified common dysregulated biological pathways and molecular targets underlying critical molecular alterations in pSS pathogenesis.
http:purl.orgdcelements1.1creator
Horvath, Steve
Nazmul-Hossain, Abu NM
Pollard, Rodney PE
Kroese, Frans GM
Vissink, Arjan
Kallenberg, Cees GM
Spijkervet, Fred KL
Bootsma, Hendrika
Michie, Sara A
Gorr, Sven U
Peck, Ammon B
Cai, Chaochao
Zhou, Hui
Wong, David TW
http:purl.orgeprinttermsisExpressedAs epdcx:valueRef sword-mets-expr-1
http:purl.orgeprintentityTypeExpression
http:purl.orgdcelements1.1language epdcx:vesURI http:purl.orgdctermsRFC3066
en
http:purl.orgeprinttermsType
http:purl.orgeprinttypeJournalArticle
http:purl.orgdctermsavailable
epdcx:sesURI http:purl.orgdctermsW3CDTF 2012-11-01
http:purl.orgdcelements1.1publisher
BioMed Central Ltd
http:purl.orgeprinttermsstatus http:purl.orgeprinttermsStatus
http:purl.orgeprintstatusPeerReviewed
http:purl.orgeprinttermscopyrightHolder
Steve Horvath et al.; licensee BioMed Central Ltd.
http:purl.orgdctermslicense
http://creativecommons.org/licenses/by/2.0
http:purl.orgdctermsaccessRights http:purl.orgeprinttermsAccessRights
http:purl.orgeprintaccessRightsOpenAccess
http:purl.orgeprinttermsbibliographicCitation
Arthritis Research & Therapy. 2012 Nov 01;14(6):R238
http:purl.orgdcelements1.1identifier
http:purl.orgdctermsURI http://dx.doi.org/10.1186/ar4081
fileSec
fileGrp sword-mets-fgrp-1 USE CONTENT
file sword-mets-fgid-0 sword-mets-file-1
FLocat LOCTYPE URL xlink:href ar4081.xml
sword-mets-fgid-1 sword-mets-file-2 applicationpdf
ar4081.pdf
sword-mets-fgid-3 sword-mets-file-3 applicationvnd.ms-excel
AR4081-S2.XLS
sword-mets-fgid-4 sword-mets-file-4 textcsv
AR4081-S4.CSV
sword-mets-fgid-5 sword-mets-file-5
AR4081-S3.CSV
sword-mets-fgid-6 sword-mets-file-6
AR4081-S5.PDF
sword-mets-fgid-7 sword-mets-file-7
AR4081-S1.XLS
sword-mets-fgid-8 sword-mets-file-8 applicationvnd.openxmlformats-officedocument.wordprocessingml.document
AR4081-S6.DOCX
structMap sword-mets-struct-1 structure LOGICAL
div sword-mets-div-1 DMDID Object
sword-mets-div-2 File
fptr FILEID
sword-mets-div-3
sword-mets-div-4
sword-mets-div-5
sword-mets-div-6
sword-mets-div-7
sword-mets-div-8
sword-mets-div-9