Evidence-Based Complementary and Alternative MedicineVolume 2012 (2012), Article ID 203043, 17 pagesdoi:10.1155/2012/203043

Research Article

Network-Based Gene Expression Biomarkers for Cold and Heat Patterns of Rheumatoid Arthritis in Traditional Chinese Medicine

Cheng Lu,¹ Xuyan Niu,¹ Cheng Xiao,² Gao Chen,1,3 Qinglin Zha,⁴ Hongtao Guo,⁵ Miao Jiang,¹ and Aiping Lu1,5

¹Institute of Basic Research in Clinical Medicine, China Academy of Chinese Medical Science, No. 16, Nanxiaojie, Dongzhimennei, Beijing 100700, China²Institute of Clinical Medicine, China-Japan Hospital, Beijing 100029, China³School of Life Science, Hubei University, Wuhan, Hubei 430062, China⁴National Research Center of TCM, Jiangxi University of Chinese Medicine, Nanchang, Jiangxi 330006, China⁵E-Institute of Shanhai Municipal Education Commission, Shanghai TCM University, Shanghai 201203, China

Received 27 August 2011; Revised 2 November 2011; Accepted 9 December 2011

Academic Editor: Lixing Lao

Copyright © 2012 Cheng Lu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

In Traditional Chinese Medicine (TCM), patients with Rheumatoid Arthritis (RA) can be classified into two main patterns: cold-pattern and heat-pattern. This paper identified the network-based gene expression biomarkers for both cold- and heat-patterns of RA. Gene expression profilings of CD4+ T cells from cold-pattern RA patients, heat-pattern RA patients, and healthy volunteers were obtained using microarray. The differentially expressed genes and related networks were explored using DAVID, GeneSpring software, and the protein-protein interactions (PPI) method. EIF4A2, CCNT1, and IL7R, which were related to the up-regulation of cell proliferation and the Jak-STAT cascade, were significant gene biomarkers of the TCM cold pattern of RA. PRKAA1, HSPA8, and LSM6, which were related to fatty acid metabolism and the I-κB kinase/NF-κB cascade, were significant biomarkers of the TCM heat-pattern of RA. The network-based gene expression biomarkers for the TCM cold- and heat-patterns may be helpful for the further stratification of RA patients when deciding on interventions or clinical trials.

1. Introduction

Rheumatoid arthritis (RA) is a systemic autoimmune disease characterized by chronic inflammation throughout the body, particularly in diarthrodial joints [1, 2], and inflammatory cell infiltration, including activated CD4+ T cells in peripheral blood [3]. There is substantial evidence that CD4+ T cells play a key role in inflammatory processes in the pathogenesis of RA [4, 5] through their ability to stimulate the secretion of proinflammatory cytokines, which induces immunoglobulin production and matrix metalloproteinase secretion and causes osteoclastogenesis [6]. A hallmark pathological feature of RA is the infiltration and accumulation of T cells in the synovium of joints [7]. Because the shared epitope in human leukocyte antigen-DR genes is found in approximately 80% of RA patients, dysregulated CD4+ T-cell activation and function have been investigated based on available genetic predisposition evidence [8, 9].

Pattern classification, a key aspect of traditional Chinese medicine (TCM) therapeutic theory, is based on symptoms, tongue appearance, and pulse feelings [10]. The efficacy of TCM is based entirely on pattern differentiation because the pattern guides the prescription [11–13]. An increasing number of medical researchers are recognizing that the combination of disease diagnosis in biomedicine and pattern classification in TCM is essential for clinical practice. This model is commonly practiced in China because it produces positive clinical effects [14]. In many clinical studies that employed TCM pattern classification, desirable outcomes were achieved [15–23]. Moreover, the correlations between TCM patterns and biomedical parameters are receiving an increasing amount of attention during the mechanism exploration step of TCM pattern classification. Such correlations include the linkage between sex hormones and Kidney deficiency syndrome in chronic nephritis [24], the linkage between the C-reactive protein (CRP) and cold and hot syndromes in rheumatoid arthritis (RA) [25], the association between the Qi-Yin deficiency syndrome in type 2 diabetes mellitus (T2DM) patients with macroangiopathy and the apolipo-protein E (APOE) E4 and E3 genotypes [26], the correlation between the serum level of eosinophil cationic protein (ECP) in asthmatic patients and the development of “heat” syndrome manifestations [27], and the close relationship between decreases in skin electrical conductance and the severity of qi vacuity [28]. Therefore, the classification of different TCM patterns for specific diseases is important in both clinical practice and related basic research.

RA patients can be classified into two main TCM patterns: the cold pattern and the heat pattern. Manifestations of the cold pattern are described as severe fixed pain in a joint or muscle, pain relief upon warming, severe pain when cooling, and a white tongue coating. In contrast, manifestations of the heat pattern are characterized by severe pain, hot, red, swollen and inflamed joints, pain relief upon cooling, severe pain upon warming, fever, thirst, restlessness, deep-colored urine, and a red tongue with a yellow coating [13, 29]. Our previous study showed that the effective rate of biomedical combination therapy was higher in patients with a cold pattern than in patients with a heat pattern (P < 0.01) [29]. These different responses to treatment led to the hypothesis that these patterns in TCM have their own specific markers.

An important part of systems biology, microarrays provide a powerful method for cataloging genes that are affected or mediated by a given disease [30, 31], such as RA [32, 33]. Increasingly informative software programs that aid in the translation of these catalogues of genes from microarray experiments into an understanding of disease pathology are available. Biological systems are both redundant and highly networked. As a consequence, many functionally interrelated genes tend to be affected when the pathological effects on a given pathway are significant. These pathways can be characterized through a secondary analysis of differentially regulated gene sets using clustering and networking algorithms and a visual analysis of gene function. Ma T coupled the classical TCM cold syndrome with methods from systems biology to explore the macro-and microsystems biology frameworks of TCM syndromes. The results of this study indicate that genes related to the cold syndrome play an essential role in energy metabolism [34]. In a previous report, we used microarray technology to reveal gene expression profiles in CD4+ T cells and showed that 29 genes were differentially regulated in RA patients with cold and heat patterns [35]. We identified the differentially expressed genes in the cold and heat pattern RA patients. However, the comparisons between the cold-pattern RA patients and healthy persons and those between the heat-pattern RA patients and healthy persons may be more important in identifying the biomarkers of the cold and heat patterns. In the present study, we applied gene expression profiling of peripheral CD4+ T cells from cold-pattern RA patients, heat-pattern RA patients, and healthy volunteers to identify the differentially expressed genes and related networks for cold- and heat-pattern RA patients based on the differences from healthy persons and to further reveal the network-based biomarkers for the cold and heat patterns using the Database for Annotation, Visualization and Integrated Discovery (DAVID), the GeneSpring Software, and a PPI (protein-protein interactions) analysis.

2. Materials and Methods

2.1. Patients

A total of 33 female RA patients from the China-Japan Friendship Hospital and 12 healthy female volunteers from the China Academy of Chinese Medical Science in Beijing, aged 18 to 70 years old, participated in the study. RA patients were eligible to participate if they had met the American College of Rheumatology (ACR) criteria for RA for at least one year, with a functional class of I, II, or III [36]. Patients were diagnosed and classified into either the heat pattern group or the cold pattern group according to TCM theory using a questionnaire, a tongue examination, and pulse diagnosis by 3 appointed TCM practitioners [13, 29]. Patients were included in the study only if the 3 practitioners reported consistent results. This ensured that all of the selected patients had typical manifestations of the heat or cold patterns according to TCM theory. Healthy females without any diagnosed diseases were included as controls.

Patients who had continuously received nonsteroidal anti-inflammatory corticosteroid drugs for more than 6 months or who had received these medicines within one month were excluded from the study. Patients with severe cardiovascular, lung, liver, kidney, mental, or blood system diseases and women who were pregnant, breast-feeding, or planning pregnancy in the next 8 months were excluded from the study. A complete joint function and biochemical function evaluation was available for all participants in the study.

2.2. Sample Preparation

For the microarray, 8 mL of venous blood was collected in anticoagulation tubes from each of the 45 participants (12 patients with the heat pattern, 21 patients with the cold pattern, and 12 healthy volunteers) before breakfast. CD4+ T cells were extracted and purified using the RosetteSep Human CD4+ T Cell Enrichment Cocktail (StemCell Technologies, Inc., Vancouver, Canada). Total RNA was isolated from the CD4+ T cells using the Trizol extraction method (Invitrogen, Carlsbad, Canada), as described by the manufacturer. mRNA was amplified linearly using the MessageAmp™ aRNA Kit (Ambion, Inc., Austin, USA), in accordance with the instructions of the manufacturer. cRNA was purified using an RNeasy Mini Kit (QIAGEN, Hilden, Germany) based on a standard procedure.

2.3. Microarray Assay

A two-color whole Human Genome Microarray Kit,  K (Agilent Technologies) was used in this study. Microarray hybridizations were performed on labeled cRNAs. All arrays had the same labels: Cy3 for samples and Cy5 for controls. The arrays were incubated at 65°C for 17 h in Agilent’s microarray hybridization chambers and subsequently washed according to the Agilent protocol. The arrays were scanned at a 5 μm resolution using a Genepix 4000B scanner (Molecular Devices Corporation, Sunnyvale, CA). Auto-photomultiplier tube (PMT) gains were adjusted to obtain a ratio of Cy3 and Cy5 channel intensities. Scanned picture information was transformed into data using the GenePix Pro Microarray Image Analysis Software. One assay was completed for each sample, and biological replication was adopted to reduce the systematic sources of variation common in macroarray studies [37].

2.4. Statistics and Functional Analysis

2.4.1. Microarray Statistics

All of the data were analyzed using the SAS9.1.3 statistical package (order no. 195557). The data were normalized to correct for technical variations among individual microarray hybridizations using the two-step procedure described in detail by Jarvis and colleagues. The signal intensity of each expressed gene was globally normalized (LOWESS) using the R statistics program [38]. Any ratio between two groups of more or less than 1 : 1.4 was taken as the differential gene expression criteria. Statistical significance was tested using the Student’s t-test. Changes greater than 1.4-fold (cold pattern to control group or heat pattern to control group) were recorded as upregulations, and those less than 1.4-fold (cold pattern to control group or heat pattern to control group) were recorded as downregulations. Other fold changes for gene expression were recorded as normal expression. Changes in gene expression (1.4-fold change) were required in more than 50% of the patients. A chi-squared test was used for these comparisons and to identify similar and different genes in the cold and heat pattern groups of differentially expressed genes.

Gene assemblages were obtained using a principal components analysis and an iterated principal factor analysis. A cut-off value of 1.5 for the foreground and background signal ratios was used. The fluorescence ratio for each spot was log-transformed for normalization. A cluster analysis was performed using Cluster 3.0 and Tree View software.

2.4.2. GeneSpring Analysis

A global comparison of all cell lines was performed using GeneSpring GX v 7.3.1 (Agilent Technologies) and the gene annotations available in March 2009 to find differentially expressed genes in the majority of resistant cell lines. Triplicate samples for the two conditions (sensitive and resistant) in each of the seven cell lines (45 samples in total) were imported into one single experiment. The expression of each gene was calculated as the ratio of the value obtained for each condition relative to the control condition after data normalization. The data were filtered using the control strength, and a control value was calculated using the Cross-Gene Error Model on replicates based on the average base/proportional value. Measurements with higher control strengths are relatively more precise than measurements with lower control strengths. Genes that did not reach this value were discarded. Additional filtering was performed to identify the differentially expressed genes. We selected the genes that displayed a P value corrected by a false discovery rate (Benjamin and Hochberg FDR) of less than 0.05.

Of the significantly differentially expressed RNA, only genes with a greater than 1.4-fold increase or 1.4-fold decrease in expression compared to the controls were used for further analysis. All microarray data in this study are in accordance with MIAME guidelines and have been deposited in the NCBI GEO database.

2.4.3. Gene Ontology and KEGG Pathway Analysis Using DAVID

We used DAVID Bioinformatics Resources 6.7 (the Database for Annotation, Visualization, and Integrated Discovery) [39], “http://david.abcc.ncifcrf.gov/”, a comprehensive set of functional annotation tools for understanding the biological meaning behind large lists of genes, to obtain gene ontology and KEGG pathway information for differential genes between the cold pattern and control samples and between the heat pattern and control samples. Differentially expressed genes that were similar and different between the cold pattern and the control and between the heat pattern and the control were compared. In the analysis, the control was set at the 0.01 level for the number of false positives using two statistics: False Detection Rate (FDR) and Bonferroni correction (FWER).

2.4.4. Protein-Protein Interaction Analysis

PPIS are the basic skeleton for the self-organization and homeostasis of living organisms [40]. In this study, information on human PPI networks from significant genes was obtained from databases, including the BIND, BIOGRID, DIP, HPRD, IntAct, and MINT, and complemented with relationships that were parsed from the literature using Agilent Literature Search. The PPI network was visualized using the Cytoscape software [41].

We integrated the databases and networks and used an IPCA to analyze the characteristics of the networks. The IPCA algorithm can detect highly connected regions or clusters in the interactome network [42]. Interactomes with a score greater than 2.0 and at least four nodes were taken as significant predictions.

Gene ontology categories were further analyzed to identify the function of each highly connected region generated by the IPCA. The latest version of the Biological Network Gene Ontology (BINGO) tool was used to statistically evaluate groups of proteins with respect to the existing annotations of the Gene Ontology Consortium. The degree of functional enrichment for a given cluster was quantitatively assessed (P value) using a hypergeometric distribution implemented in the BINGO tool.

3. Results

3.1. Differential Gene Expression : RA with the Cold Pattern versus the Healthy Control and RA with the Heat Pattern versus the Healthy Control

A comparison between RA patients with the cold pattern and the healthy control revealed 35 differentially regulated genes (Figure 1, Table 1). Among these differentially regulated genes, 16 were upregulated (fold ≥ 1.4) and 19 were downregulated (fold ≤ 1.4).

Table 1: Differentially expressed genes in cold-pattern RA patient versus healthy control.

Gene name
Definition
Accession no.
Cold/control ratio

Upregulated gene (>1.4)

    MMGT1
Membrane magnesium transporter 1
AL157477
1.70

    TDRD7
Tudor domain containing 7
AB025254
1.56

    SLC4A7
Solute carrier family 4, sodium bicarbonate Cotransporter, member 7
AF047033
1.54

    MRPL28
Mitochondrial ribosomal protein L28
NM_006428
1.54

    CCNT1
Cyclin T1
NM_001240
1.52

    GTF3C6
General transcription factor IIIC, polypeptide 6, alpha 35 kDa
BC011593
1.49

    BCL2A1
BCL2-related protein A1
NM_004049
1.45

    TRIM22
Tripartite motif-containing 22
NM_006074
1.45

    CTLA4
Cytotoxic T-lymphocyte-associated protein 4
AF414120
1.45

    KLRAQ1
KLRAQ motif containing 1
BC011978
1.45

    GBP3
Guanylate binding protein 3
AL136680
1.43

    MIS12
MIS12, MIND kinetochore complex component, homolog (S. pombe)
NM_024039
1.43

    PDCD2
Programmed cell death 2
AJ420535
1.43

    ATOX1
ATX1 antioxidant protein 1 homolog (yeast)
NM_004045
1.43

    PIGB
Phosphatidylinositol glycan, class B
NM_004855
1.41

    NPEPPS
Aminopeptidase puromycin sensitive
NM_006310
1.41

Downregulated gene (<1 : 1.4)

    Unknown
Homo sapiens cDNA FLJ11459 fis, clone HEMBA1001566
AK021521
0.59

    Unknown
Homo sapiens mRNA; cDNA DKFZp434I138 (from clone DKFZp434I138)
AL122121
0.66

    Unknown
Homo sapiens clone 24629 mRNA sequence
AF052160
0.66

    ZNF330
Zinc finger protein 330
NM_014487
0.67

    PLP2
Proteolipid protein 2 (colonic epithelium-enriched)
NM_002668
0.68

    PTGES2
Prostaglandin E synthase 2
NM_025072
0.68

    PSMD8
Proteasome (prosome, macropain) 26S subunit, non-ATPase, 8
NM_002812
0.68

    IL7R
Interleukin 7 receptor
NM_002185
0.68

    STMN3
Stathmin-like 3
NM_015894
0.68

    ABCC10
ATP-binding cassette, sub-family C (CFTR/MRP), member 10
AK000002
0.69

    ELAC2
ElaC homolog 2 (E. coli)
NM_018127
0.69

    PLEKHA1
Pleckstrin homology domain-containing, family A (phosphoinositide binding specific) member 1
NM_021622
0.69

    DCXR
Dicarbonyl/L-xylulose reductase
NM_016286
0.70

    EIF4A2
Eukaryotic translation initiation factor 4A2
NM_001967
0.70

    IL16
Interleukin 16 (lymphocyte chemoattractant factor)
M90391
0.70

    GRAMD1B
GRAM domain containing 1B
AB033027
0.71

    LAGE3
L antigen family, member 3
NM_006014
0.71

    TMEM80
Transmembrane protein 80
BC008671
0.71

    SLIT3
Slit homolog 3 (Drosophila)
AL122074
0.71

(a)

(b)

Figure 1: Cluster diagram of the expression of significantly expressed genes in 21 RA patients with the cold pattern, 12 RA patients with the heat pattern, and 12 healthy volunteers. Patients are indicated as vertical column headings, and gene symbols of transcripts are given in horizontal rows. Red represents a relative expression greater than the median expression level across all samples, and green represents an expression level lower than the median. (a) Significantly expressed genes between TCM cold-pattern RA patients and healthy control. (b) Significantly expressed genes between TCM heat-pattern RA patients and healthy control.

An analysis of RA patients with the heat pattern versus the healthy control showed 21 differentially expressed genes (Figure 1, Table 2); 15 were upregulated (fold ≥ 1.4), and 6 were downregulated (fold ≤ 1.4) (Figure 1, Table 2).

Table 2: Differentially expressed genes in heat-pattern RA patient versus healthy control.

Gene name
Definition
Accession no.
Heat/control ratio

Up-regulated gene (>1.4)

    CAMP
Cathelicidin antimicrobial peptide
NM_004345
2.33

    MMGT1
Membrane magnesium transporter 1
AL157477
1.69

    TDRD7
Tudor domain containing 7
AB025254
1.56

    LSM6
LSM6 homolog, U6 small nuclear RNA associated (S. cerevisiae)
NM_007080
1.54

    CCDC55
Coiled-coil domain containing 55
NM_032141
1.54

    GTF3C6
General transcription factor IIIC, polypeptide 6, alpha 35 kDa
BC011593
1.54

    CD59
CD59 molecule, complement regulatory protein
NM_000611
1.52

    CTLA4
Cytotoxic T-lymphocyte-associated protein 4
AF414120
1.52

    ASCC3
Activating signal cointegrator 1 complex subunit 3
AY013288
1.47

    PRKAA1
Protein kinase, AMP-activated, alpha 1 catalytic subunit
AK024252
1.45

    BCL2A1
BCL2-related protein A1
NM_004049
1.45

    TMEM60
Transmembrane protein 60
NM_032936
1.43

    HSPA1A
Heat shock 70 kD protein 1A
NM_005345
1.43

    PARP9
Poly (ADP-ribose) polymerase family, member 9
NM_031458
1.43

    HIST1H1C
histone cluster 1, H1c
NM_005319
1.43

Down-regulated gene (<1 : 1.4)

    HSPA8
Heat shock 70 kD protein 8
NM_006597
0.63

    SLC43A3
Solute carrier family 43, member 3
NM_014096
0.67

    Unknown
Homo sapiens clone 24629 mRNA sequence
AF052160
0.67

    PSMD8
Proteasome (prosome, macropain) 26S subunit, non-ATPase, 8
NM_002812
0.69

    NCL
Nucleolin
NM_005381
0.70

    LBH
Limb bud and heart development homolog (mouse)
NM_030915
0.71

Six genes showed the same pattern of differential expression between the cold and control groups and between the heat and control groups (Table 3): MMGT1, TDRD7, GTF3C6, BCL2A1, CTLA4, and PSMD8. Except for PSMD8, which was downregulated, these genes were upregulated in both the cold pattern and the heat pattern comparisons.

Table 3: The same significantly expressed genes found in the differences between the cold and control groups and between the heat and control groups.

Gene name
Definition
Accession no.
Ratio in cold versus control
Ratio in heat versus control

MMGT1
Membrane magnesium transporter 1
AL157477
1.70
1.69

TDRD7
Tudor domain containing 7
AB025254
1.56
1.56

GTF3C6
General transcription factor IIIC, polypeptide 6, alpha 35 kDa
BC011593
1.49
1.54

BCL2A1
BCL2-related protein A1
NM_004049
1.45
1.45

CTLA4
Cytotoxic T-lymphocyte-associated protein 4
AF414120
1.45
1.52

PSMD8
Proteasome (prosome, macropain) 26S subunit, non-ATPase, 8
NM_002812
0.68
0.69

3.2. Pathway Analysis Using DAVID

The KEGG pathway analyses of the significantly expressed genes using DAVID are shown in Tables 4, 5, and 6. The six genes shared between the cold pattern and the heat pattern comparisons were related to the following pathways (Table 4): CAMs, T-cell receptor signaling pathway, and the proteasome. CTLA4 participates in the pathways of autoimmune thyroid disease, CAMs, and the T-cell receptor signaling and rheumatoid arthritis pathways. PSMD8 is one subunit of the proteasome. These results showed that the RA patients with either the cold or heat pattern had the same five pathways, which were disordered compared to healthy persons. These were the common points in the pathways of the cold and heat pattern comparisons.

Table 4: The pathways shared by the differently expressed genes of the cold pattern versus control and heat pattern versus control comparisons.

KEGG pathway
Related genes
Ratio in cold versus control
Ratio in heat versus control

hsa05320: Autoimmune thyroid disease
CTLA4
1.45
1.52

hsa04514: Cell adhesion molecules (CAMs)
CTLA4
1.45
1.52

hsa04660: T cell receptor signaling pathway
CTLA4
1.45
1.52

hsa05323: Rheumatoid arthritis
CTLA4
1.45
1.52

hsa03050: Proteasome
PSMD8
0.68
0.69

Table 5: The pathways related to the differently expressed genes between the cold pattern and control comparison.

KEGG pathway
Related genes
Ratio

hsa00563: Glycosylphosphatidylinositol (GPI) anchor biosynthesis
PIGB
1.41

hsa00590: Arachidonic acid metabolism
PTGES2
0.68

hsa04630: Jak-STAT signaling pathway
IL7R
0.68

hsa04640: Hematopoietic cell lineage
IL7R
0.68

hsa05340: Primary immunodeficiency
IL7R
0.68

hsa04060: Cytokine-cytokine receptor interaction
IL7R
0.68

hsa02010: ABC transporters
ABCC10
0.69

hsa00040: Pentose and glucuronate interconversions
DCXR
0.70

hsa04360: Axon guidance
SLIT3
0.71

Table 6: The pathways related to the differently expressed genes between the heat pattern and control comparison.

KEGG pathway
Related genes
Ratio

hsa03040: Spliceosome
LSM6, HSPA1A, HSPA8
1.54, 1.43, 0.63

hsa04612: Antigen processing and presentation
HSPA1A, HSPA8
1.43, 0.63

hsa04144: Endocytosis
HSPA1A, HSPA8
1.43, 0.63

hsa04010: MAPK signaling pathway
HSPA1A, HSPA8
1.43, 0.63

hsa03018: RNA degradation
LSM6
1.54

hsa04640: Hematopoietic cell lineage
CD59
1.52

hsa04610: Complement and coagulation cascades
CD59
1.52

hsa04150: mTOR signaling pathway
PRKAA1
1.45

hsa04920: Adipocytokine signaling pathway
PRKAA1
1.45

hsa04140: Regulation of autophagy
PRKAA1
1.45

hsa05410: Hypertrophic cardiomyopathy (HCM)
PRKAA1
1.45

hsa04910: Insulin signaling pathway
PRKAA1
1.45

hsa05020: Prion diseases
HSPA1A
1.43

hsa05130: Pathogenic Escherichia coli infection
NCL
0.70

Compared to the heat pattern versus control analysis, the cold pattern versus control analysis revealed that different pathways were related to the following differentially expressed genes (Table 5): glycosylphosphatidylinositol (GPI) anchor biosynthesis, arachidonic acid metabolism, Jak-STAT signaling, hematopoietic cell lineage, primary immunodeficiency, cytokine-cytokine receptor interactions, ABC transporters, pentose and glucuronate interconversions, and axon guidance. These pathways, which were related to the differentially expressed genes, may be characteristic of the cold pattern.

The differentially expressed genes identified in the heat pattern versus control group analysis were related to the following pathways (Table 6): antigen processing and presentation, endocytosis, MAPK signaling, RNA degradation, hematopoietic cell lineage, complement and coagulation cascades, mTOR signaling, adipocytokine signaling, regulation of autophagy, and insulin signaling. These results may reveal characteristics of the heat pattern.

3.3. Results of Gene Ontology Analysis Using DAVID

The GO-discovered categories using DAVID analysis for similar and different significantly expressed genes between the cold pattern and the heat pattern comparisons are shown in Tables 7, 8, and 9. As shown in Table 7, the six shared differentially expressed genes of the cold and heat pattern comparisons were predominantly grouped into functional classes of protein binding and binding. BCL2A1, CTLA4, and PSMD8 belong to biological processes of negative cellular regulation, and TDRD7, GTF3C6, and PSMD8 are cellular constituents of intracellular organelle lumens. The results of these gene ontology analyses were common to the cold and heat patterns.

Table 7: GO-discovered categories for the same expressed genes in the cold versus control and heat versus control groups.

GO ID
GO name
Number of genes
P value
Gene name

Biological process

    GO:0048523
Negative regulation of cellular process
3
0.1085
BCL2A1, CTLA4, PSMD8

    GO:0048519
Negative regulation of biological process
3
0.1264
BCL2A1, CTLA4, PSMD8

Cellular constituent

    GO:0070013
Intracellular organelle lumen
3
0.0993
TDRD7, GTF3C6, PSMD8

Molecular function

    GO:0005515
Protein binding
5
0.2392
TDRD7, BCL2A1, GTF3C6, CTLA4, PSMD8

    GO:0005488
Binding
6
0.3880
MMGT1, TDRD7, BCL2A1, GTF3C6, CTLA4, PSMD8

Table 8: GO-discovered categories for the genes that differed between the cold pattern and the control.

GO ID
GO name
Number of genes
P value
Gene name

Biological process

    GO:0044419
Interspecies interaction between organisms
3
0.0600
IL16, EIF4A2, CCNT1

    GO:0051704
Multiorganism process
4
0.0692
IL16, EIF4A2, CCNT1, TRIM22

    GO:0006955
Immune response
4
0.0714
IL16, IL7R, TRIM22, GBP3

Cellular constituent

    GO:0005902
Ribonucleoprotein complex
2
0.0508
SLC4A7, DCXR

Molecular function

    GO:0005515
Protein binding
18
0.0210
PLP2, ELAC2, PTGES2, IL16, STMN3, ATOX1, CCNT1, NPEPPS, IL7R, TRIM22, MIS12, ZNF330, PDCD2, SLIT3, LAGE3, EIF4A2, SLC4A7, DCXR

    GO:0005488
Binding
22
0.0866
PLP2, ELAC2, PTGES2, IL16, STMN3, ATOX1, CCNT1, ABCC10, NPEPPS, IL7R, TRIM22, MIS12, ZNF330, PDCD2, SLIT3, MRPL28, LAGE3, EIF4A2, SLC4A7, GBP3, DCXR, PLEKHA1

Table 9: GO-discovered categories for the genes that differed between the heat pattern and the control.

GO ID
GO name
Number of genes
P value
Gene name

Biological process

    GO:0006950
Response to stress
5
0.0235
CAMP, CD59, PRKAA1, HSPA1A, HSPA8

    GO:0009607
Response to biotic stimulus
3
0.0287
CAMP, HSPA1A, HSPA8

    GO:0006986
Response to unfolded protein
2
0.0492
HSPA1A, HSPA8

    GO:0051789
Response to protein stimulus
2
0.0733
HSPA1A, HSPA8

Cellular constituent

    GO:0030529
Ribonucleoprotein complex
4
0.0076
LSM6, HSPA1A, NCL, HSPA8

Molecular function

    GO:0000166
Nucleotide binding
5
0.0205
ASCC3, PRKAA1, HSPA1A, NCL, HSPA8

    GO:0005524
ATP binding
4
0.0357
ASCC3, PRKAA1, HSPA1A, HSPA8

    GO:0032559
Adenyl ribonucleotide binding
4
0.0369
ASCC3, PRKAA1, HSPA1A, HSPA8

    GO:0030554
Adenyl nucleotide binding
4
0.0423
ASCC3, PRKAA1, HSPA1A, HSPA8

    GO:0001883
Purine nucleoside binding
4
0.0440
ASCC3, PRKAA1, HSPA1A, HSPA8

    GO:0001882
Nucleoside binding
4
0.0447
ASCC3, PRKAA1, HSPA1A, HSPA8

    GO:0051082
Unfolded protein binding
2
0.0592
HSPA1A, HSPA8

    GO:0032553
Ribonucleotide binding
4
0.0623
ASCC3, PRKAA1, HSPA1A, HSPA8

    GO:0032555
Purine ribonucleotide binding
4
0.0623
ASCC3, PRKAA1, HSPA1A, HSPA8

    GO:0017076
Purine nucleotide binding
4
0.0695
ASCC3, PRKAA1, HSPA1A, HSPA8

Genes that were uniquely differentially expressed between the cold pattern and the control were predominantly grouped into functional classes of protein binding (18 genes) and binding (22 genes) (Table 8). A few genes belonged to interspecies interactions between organisms (IL16, EIF4A2, and CCNT1), multiorganism processes (IL16, EIF4A2, CCNT1, TRIM22), and the immune response (IL16, IL7R, TRIM22, GBP3). SLC4A7 and DCXR are cellular constituents of the ribonucleoprotein complex.

As shown in Table 9, the genes that were uniquely differentially expressed between the heat pattern and the control were predominantly grouped into functional binding classes, including nucleotide binding, ATP binding, adenyl ribonucleotide binding, adenyl nucleotide binding, purine nucleoside binding, nucleoside binding, unfolded protein binding, ribonucleotide binding, purine ribonucleotide binding, and purine nucleotide binding. A few genes belonged to the biological stress response process, the biotic stimuli response, the unfolded protein response, and the protein stimulus and cellular constituent response of the ribonucleoprotein complex.

The common point of GO-discovered categories in the cold and heat pattern comparisons was the cellular constituent of the ribonucleoprotein complex. The characteristic molecular function is protein binding in the cold pattern, whereas it is nucleoside and ribonucleotide binding in the heat pattern.

3.4. Gene Function Using GeneSpring Analysis

The networks in Figure 2 were constructed using the GeneSpring GX v 7.3.1 software, as described in the Methods Section, starting with the lists of genes that were differently expressed between the cold pattern and the control (network A) and between the heat pattern and the control (network B). These networks revealed the functions and biological processes of the significantly expressed genes.

(a)

(b)

(c)

(d)

(e)

Figure 2: Networks of cold- and heat-pattern comparisons based on GeneSpring GX v 7.3.1 software. (a) Differently expressed genes between the cold pattern and the control. (b) Differently expressed genes between the heat pattern and the control. (c) Differently expressed genes between the cold pattern and the control in T-cell regulation. (d) Differently expressed genes between the heat pattern and the control in T-cell regulation.

CTLA4, PSMD8, and TDRD7 were in both the cold and heat pattern comparisons, similar to the DAVID analysis results. CTLA4 and CCNT1 are important centers of the cold pattern networks, and they cooperated with the IL7R, IL16, and EIF4A2 proteins and participated in the negative regulation of T cells, T cell homeostasis, intracellular tyrosine, Janus Kinase 3, RNA elongation, and transcription.

In the network of differently expressed genes between the heat pattern and the control, CAMP, HSPA8, and HSPA1A were highly connected with other differently expressed genes and were related to the negative regulation of T cells, protein refolding, keratinocyte migration, and neutrophil apoptosis.

Networks C and D revealed that the differently expressed genes in the cold and heat pattern comparisons were related to the regulation of T cells but in different ways. In the cold pattern, IL7R, CD80, and IL-6 processed the negative regulation of T-cell activity and T-cell homeostasis through CTLA4. However, in the heat pattern, except for CTLA4, which was related to the negative regulation of T cell activation, CAMP, SOCS1, HSPA1A, and TLR7 participated in T-cell-mediated immunity, macrophage activation, and keratinocyte migration.

3.5. Results of PPI Analysis

The PPI networks of the significantly expressed genes between the cold pattern and the control and between the heat pattern and the control are shown in Figure 3. There were four subnetworks for the cold-pattern comparison and four subnetworks for the heat patten comparison. The gene ontology analyses of each subnetwork are shown in Table 10.

Table 10: The gene ontology analysis of the cold pattern versus control and the heat pattern versus control comparisons.

GO ID
Description
P value

Cold pattern versus control

Subnetworks A

    6455
Translational elongation
4.26E-168

    6416
Translation
6.96E-146

    44249
Cell biosynthetic process
1.87E-104

    9059
Macromolecule biosynthetic
8.699E-70

    10467
Gene expression
1.156E-68

Subnetworks B

    31145
Anaphase-promoting complex-dependent proteasomal ubiquitin-depended protein catabolic process
4.326E-81

    51436
Negative regulation of ubiquitin protein ligase activity during mitotic cell cycle
4.326E-81

    51352
Negative regulation of ligase activity
1.235E-80

    51444
Positive regulation of ligase activity during mitotic cell cycle
1.235E-80

    51437
Positive regulation of ligase activity
9.253E-80

Subnetworks C

    6374
Nuclear mRNA splicing, via spliceosome
9.582E-17

    375
RNA splicing, via transesterification reactions
9.582E-17

    377
RNA splicing, via transesterification reactions with bulged adenosine as nucleophile
9.582E-17

    6395
RNA splicing
1.284E-14

    6397
mRNA process
2.703E-14

Subnetworks D

    7259
JAK-STAT cascade
2.231E-09

    19221
Cytokine and chemokine mediated signaling pathway
1.603E-08

    1553
Luteinization
5.211E-07

    1779
Natural killer cell differentiation
5.211E-07

    46544
Development of secondary male sexual characteristics
5.211E-07

Heat pattern versus control

Subnetworks E

    31145
Anaphase-promoting complex-dependent proteasomal ubiquitin-depended protein catabolic process
1.45E-105

    51436
Negative regulation of ubiquitin protein ligase activity during mitotic cell cycle
1.45E-105

    51352
Negative regulation of ligase activity
4.16E-105

    51444
Positive regulation of ligase activity during mitotic cell cycle
4.16E-105

    51437
Positive regulation of ligase activity
3.14E-104

Subnetworks F

    51101
Regulation of DNA binding
1.508E-05

    51098
Regulation of binding
3.286E-05

    43123
Positive regulation of I-κB kinase/NF-κB cascade
5.187E-05

    9607
Response to biotic stimulus
5.593E-05

    43122
Regulation of I-κB kinase/NF-κB cascade
6.457E-05

Subnetworks G

    6395
RNA splicing
1.284E-14

    6397
mRNA splicing
2.703E-14

    16071
mRNA metabolic process
7.908E-14

    6394
RNA processing
1.458E-12

    16070
RNA metabolic process
1.175E-10

Subnetworks H

    6633
Fatty acid biosynthetic process
1.498E-09

    46394
Carboxylic acid biosynthetic process
2.761E-09

    16053
Organic acid biosynthetic process
2.761E-09

    6631
Fatty acid metabolic process
7.384E-08

    32787
Monocarboxylic acid metabolic process
2.767E-07

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 3: The PPI subnetworks based on the differently expressed genes made up of highly connected regions. (a–d): The subnetworks of the cold pattern versus control comparison. (e–h): The subnetworks of the heat pattern versus control comparison. Diamonds represent seed nodes. Cycles represent neighbor nodes. All edges represent interactions.

The PPI subnetworks of the cold-pattern versus control comparison (subnetworks A-D) revealed the seed nodes EIFA2, CCNT1, and PSMD8 and their related proteins, which were similar to the differently expressed genes in the aformentioned analysis. The PPI subnetworks of the heat pattern versus control comparison (subnetworks E-H) revealed the seed nodes PSMD8, HSPA8, LSM6, and PRKAB2, which were also similar to the significantly differently expressed genes in the heat pattern and control comparison. The seed proteins and subnetworks related to the gene ontology analysis results were matched with the KEGG pathways, the gene ontology, and the functional analysis results using DAVID and GeneSpring based on the differently expressed genes in the cold- and heat-pattern comparisons (Table 10).

3.6. The Intersection Network between the Cold Pattern and the Heat Pattern

Figure 4 shows the intersection network between the cold pattern and the heat pattern, which includes the majority of differently expressed genes in the cold and heat pattern comparisons. All subnetworks are shown in Figures 3 and 4.

Figure 4: The PPI network of functional relationships between the cold and heat patterns. Diamonds represent seed nodes. Cycles represent neighbor nodes. All edges represent interactions. (a) Is the network specific to the cold pattern. (c) Is the network specific to the heat pattern. (b) Is the commonly shared network of both patterns, which is likely the common network of RA. The parts between (a) and (b) and between (b) and (c), as shown by arrows, are the connections between the highly condensed TCM pattern-related networks and the common shared network.

Based on the gene ontology analysis results (Table 10) and the other results, Figure 5 was drawn to outline the networks of the cold and heat patterns and to reveal the relationships between the two patterns in the functional and biological processes of the network.

Figure 5: The gene ontology relationship between the TCM cold and heat patterns in RA. Cycles on the left represent cold-pattern-related biological processes. Cycles on the right represent heat-pattern-related biological processes. Cycles in the middle represent the biological processes common to the cold and heat patterns. Arrows represent the direction of the regulation of differently expressed genes: up signifies up-regulation, and down signifies down-regulation.

In these networks, protein ubiquitination and RNA splicing were the common biological processes in the TCM cold and heat patterns of RA. The different biomarkers for the TCM cold and heat patterns of RA are obvious: the cold pattern was related to the regulation of translation and the Jak-STAT cascade, while the heat pattern was related to fatty acid metabolism and the I-κB kinase/NF-κB cascade. In addition, protein ubiquitination, proliferation, and apoptosis related to the cell cycle are the biological process connections between the cold pattern and the heat pattern. CTLA4 and PSMD8 were the significant biomarkers in both the cold and heat patterns. The significant biomarkers of the cold pattern were EIF4A2, CCNT1, and IL7R, while the significant biomarkers of the heat pattern were PRKAA1, HSPA8, and LSM6.

4. Discussion

In this study, we utilized software to identify network-based gene expression biomarkers and biomarkers that were organized by sets of differentially expressed genes that were members of established functional networks. Our major findings were the network-based gene expression biomarker pathways that were similar and different between the TCM cold-pattern RA patients and the TCM heat-pattern RA patients.

The CAMs, T-cell receptor signaling pathway, and proteasome could be related to both the TCM cold and heat patterns in RA patients. Specifically, CTLA4 (cytotoxic T lymphocyte antigen 4), which was a seed gene in these pathways and was up-regulated in both patterns, participates in the pathways of CAMs and T-cell receptor signaling. The CTLA4 molecule is expressed on activated T lymphocytes and has recently been identified as an important negative regulator in autoimmune diseases. Quantitative alterations of CTLA4 contribute to autoimmune tissue destruction [43], and the expression of CTLA4 plays a downregulatory role in rheumatoid articular inflammation. CTLA4 was also up-regulated on T cells in patients with RA, and the increase in CTLA4 expression might exert a downregulation effect on tumor necrosis factor alpha and interleukin 1 beta production [44]. CTLA4-Ig in activated macrophages induces significant down-regulation in the cellular production of IL-6, TNF-alpha, IL1-beta, and TGF-beta for the treatment of RA [45]. In the GeneSpring analysis results, the functions of CTLA4 and its relationships with other immune molecules, such as IL-16, IL-7R, FOXP3, and CD80, which participates in T-cell homeostasis and the negative regulation of T-cell activation in cold-pattern RA patients, were identified. In the heat pattern RA patients, except for the negative regulation of T-cell activation, CTLA4 participates in T-cell-mediated immunity, macrophage activation, and keratinocyte migration with the other seed genes, such as CAMP, SOCS1, HSPA1A, and TLR9. Our study showed that CTLA4 is an important common negative regulator in both cold- and heat-pattern RA patients compared to healthy controls.

PSMD8 (Proteasome (prosome, macropain) 26S subunit, non-ATPase, 8), a significantly expressed gene in both cold- and heat-pattern RA patients, is one subunit of a protein-destroying apparatus that is involved in many essential cellular functions, such as the regulation of the cell cycle, cell differentiation, signal transduction pathways, antigen processing for appropriate immune responses, stress signaling, inflammatory responses, and apoptosis [46]. PSMD8 was down-regulated in both cold- and heat-pattern RA patients compared to healthy controls. The PPI analysis (Figure 3) showed that in the subnetwork B of cold-pattern RA patients and the subnetwork E of heat pattern RA patients, the PSMD8-related family showed similar biological functions; it was involved in the regulation of protein ubiquitination in the cell cycle. Therefore, in RA patients, the regulation of protein ubiquitination in the cell cycle is down-regulated in both cold- and heat-pattern patients. In PPI sub-networks C and G, a similar biological process, RNA splicing, was observed in both cold-pattern and heat-pattern RA patients.

In TCM cold-pattern RA patients, pathways related to GPI anchor biosynthesis, arachidonic acid metabolism, Jak-STAT signaling, hematopoietic cell lineage, primary immunodeficiency, cytokine-cytokine receptor interaction, ABC transporters, pentose and glucuronate interconversions, and axon guidance were found. In these pathways, CCNT1, IL7R, IL16, and EIF4A2 genes were included as the seeds.

CCNT1, or Cyclin T1, is a protein in the highly conserved cyclin family, whose members are characterized by a dramatic periodicity in protein abundance throughout the cell cycle. Cyclins function as regulators of CDK kinases. Different cyclins exhibit distinct expression and degradation patterns, which contribute to the temporal coordination of each mitotic event. Cyclin T1 is closely associated with CDK9 kinase and is a major subunit of the transcription elongation factor p-TEFb. This cyclin and its kinase partner are involved in the phosphorylation and regulation of the carboxy-terminal domain (CTD) of the largest RNA polymerase II subunit [47]. Cyclin T1 protein expression is highly regulated in CD4+ T cells and macrophages. Cyclin T1 expression is low in resting CD4+ T cells that have been isolated from healthy donors, but upon T-cell activation, it is induced by a mechanism that involves posttranscriptional regulation [48, 49]. Cyclin T1 expression is also low in freshly isolated monocytes, and it is up-regulated by a posttranscriptional mechanism within one to two days after the cells are cultured under conditions that allow for macrophage differentiation [50]. However, after one to two weeks in culture, Cyclin T1 protein expression is shut off in macrophages by proteasome-mediated proteolysis. Treatment of macrophages with the immunosuppressive cytokine IL-10 accelerates this proteasome-mediated shut-off of Cyclin T1 [51]. Cyclin T1 can be reinduced by activation with agents, such as LPS, which indicates that the induction of Cyclin T1 is a component of an innate immune response in mature macrophages [52, 53]. In our study, CCNT1, which was up-regulated, participated in RNA transcription and protein ubiquitination in cold-pattern RA patients, but no change was observed in heat-pattern RA patients.

IL7R (Interleukin 7 receptor) participated in the Jak-STAT signaling pathway, hematopoietic cell lineage, primary immunodeficiency, and cytokine-cytokine receptor interactions in the KEGG pathway database for cold-pattern RA patients. The protein encoded by this gene is a receptor for interleukin 7 (IL7). The function of this receptor requires the interleukin 2 receptor gamma chain (IL2RG), a common gamma chain shared by the receptors of various cytokines, including interleukins 2, 4, 7, 9, and 15. This protein plays a critical role in the V (D) J recombination during lymphocyte development. This protein also controls the accessibility of the TCR gamma locus by STAT5 and histone acetylation. Knockout studies in mice suggest that apoptosis blockade is an essential function of this protein during the differentiation and activation of T lymphocytes. The functional defects in this protein may be associated with the pathogenesis of severe combined immunodeficiency (SCID) [54]. NF-κB-dependent gene expression in peripheral leukocytes is highly correlated with RA activity, as measured by the Disease Activity Score and C-reactive protein (DAS28-CRP). IL7R was one of the notably expressed associated with DAS28-CRP in the evaluation of the peripheral blood expression of genes regulated by NF-κB, a key mediator of tumor necrosis factor-alpha (TNF-alpha) signaling, in patients with RA before and during treatment with anti-TNF-alpha or methotrexate (MTX). These results identify candidate markers, such as IL7R, that could lead to the development of a simple, minimally invasive pharmacodynamic assay for RA treatments directed at the NF-κB pathway [55]. In our study, according to the GeneSpring analysis, IL7R was down-regulated and participated in T-cell regulation in cold-pattern RA patients but not in heat-pattern RA patients. In this case, IL7R can block apoptosis and promote the proliferation of CD4+ T cells in cold-pattern RA patients.

IL16 (Interleukin 16, a lymphocyte chemoattractant factor) is a pleiotropic cytokine that functions as a chemoattractant, a modulator of T-cell activation, and an inhibitor of HIV replication. The signaling process of this cytokine is mediated by CD4. The product of this gene undergoes proteolytic processing, which yields two functional proteins. The cytokine function is exclusively attributed to the secreted C-terminal peptide, and the N-terminal product may play a role in cell cycle control. Caspase 3 is involved in the proteolytic processing of this protein. Alternate splicing results in multiple transcript variants [56]. Many studies have shown that IL16 plays a role in the disease process underlying RA and joint destruction [57–59]. In our study, IL16 was down-regulated and participated in T-cell regulation with IL7R in cold-pattern RA patients but not in heat pattern RA patients.

EIF4A2 (eukaryotic translation initiation factor 4A2) is a gene for one of the protein-synthesis initiation factors involved in the binding of mRNA to the ribosome [60]. One study indicated that EIF4A2 controls mRNA-specific translation and the protein synthesis rate in pancreatic beta-cells and that EIF4A2 is down-regulated by high glucose levels in rat beta-INS832/13 cells [61]. In our study, according to the PPI analysis, EIF4A2 as a seed gene was down-regulated in cold-pattern RA patients and was related to the regulation of translation and cell biosynthetic processes, but it was not altered in heat-pattern RA patients.

For TCM heat-pattern RA patients, closely related pathways, including antigen processing and presentation, endocytosis, MAPK signaling, RNA degradation, hematopoietic cell lineage, complement and coagulation cascades, mTOR signaling, adipocytokine signaling, regulation of autophagy, hypertrophic cardiomyopathy (HCM), and insulin signaling, were found. The differentially expressed major genes were CAMP, PRKAA1, HSPA1A, HSPA8, and LSM6 in heat-pattern RA patients.

CAMP (cathelicidin antimicrobial peptide) is an antimicrobial protein found in specific granules of polymorphonuclear leukocytes (PMNs) [62]. The cathelicidin family is characterized by a conserved N-terminal cathelin domain and a variable C-terminal antimicrobial domain that can be released from the precursor protein after cleavage by proteinases [63]. In our study, CAMP was up-regulated 2.33-fold and participated in T-cell regulation and cell proliferation in TCM heat-pattern RA patients but not in cold-pattern RA patients.

PRKAA1 (protein kinase, AMP-activated, alpha 1 catalytic subunit) belongs to the ser/thr protein kinase family. It is the catalytic subunit of the 5’-prime-AMP-activated protein kinase (AMPK). AMPK regulates the activities of a number of key metabolic enzymes through phosphorylation and protects cells from stresses that cause ATP depletion by switching off ATP-consuming biosynthetic pathways [64]. It participated in the mTOR signaling pathway, the adipocytokine signaling pathway, the regulation of autophagy, hypertrophic cardiomyopathy (HCM), and the insulin signaling pathway of the KEGG pathway in heat-pattern, RA patients but not in cold-pattern RA patients. According to the PPI analysis, PRKAA1, as a seed gene, was up-regulated and participated in fatty acid metabolism in TCM heat-pattern RA patients.

HSPA1A (heat shock 70 kDa protein 1A) is a 70 kDa heat shock protein and a member of the heat shock protein 70 family. In conjunction with other heat shock proteins, this protein stabilizes existing proteins against aggregation and mediates the folding of newly translated proteins in the cytosol and in organelles. It is also involved in the ubiquitin-proteasome pathway through interaction with the AU-rich element RNA-binding protein 1. It is located in the major histocompatibility complex class III region in a cluster with two closely related genes that encode similar proteins [65]. According to the GeneSpring analysis results, HSPA1A was up-regulated and participated in spliceosome, antigen processing and presentation, endocytosis, MAPK signaling pathway, and prion diseases using the KEGG pathway and in T-cell regulation in heat pattern RA patients. Based on the PPI analysis results, it also participated in the I-κB kinase/NF-κB cascade in heat pattern RA patients. HSPA8 (heat shock 70 kDa protein 8) belongs to the heat shock protein 70 family, which contains both heat-inducible and constitutively expressed members. The latter members are called heat-shock cognate proteins. This gene encodes a heat-shock cognate protein, which binds to nascent polypeptides to facilitate correct folding. It also functions as an ATPase in the disassembly of clathrin-coated vesicles during the transport of membrane components through the cell [66]. Similar to HSPA1A, HSPA8 also participated in the spliceosome, antigen processing and presentation, endocytosis, and MAPK signaling pathway in heat-pattern RA patients. According to the PPI analysis results, HSPA8, as a seed gene, was down-regulated and also participated in the I-κB kinase/NF-κB cascade in TCM heat-pattern RA patients. LSM6 is related to Sm-like proteins, which form a stable heteromer that is present in tri-snRNP particles, and is important for pre-mRNA splicing [67]. According to the KEGG pathway, mRNA splicing, and PPI analysis results, LSM6, as a seed gene, was up-regulated and participated in the spliceosome and RNA degradation in TCM heat-pattern RA patients.

We outlined the gene ontology relationship between the TCM cold and heat patterns in RA patients (Figure 5) by first combining the PPI results of interactions between TCM cold and heat patterns (Figure 4) and then combining them with all other data. In this network, common genes and biological processes and significantly different gene-based pathways and biological processes between the TCM cold- and heat-pattern patients were identified, and these can be used as biomarkers for the classification of the two TCM RA patterns. In the network, protein ubiquitination and RNA splicing were the common biological processes in the TCM cold and heat patterns. The different biomarkers for the TCM cold and heat patterns were obvious: the TCM cold pattern was related to the regulation of translation and the Jak-STAT cascade, and the TCM heat pattern was related to fatty acid metabolism and the I-κB kinase/NF-κB cascade. In addition, protein ubiquitination, proliferation and apoptosis related to the cell cycle may be the connecting biological processes between the TCM cold and heat patterns. CTLA4 and PSMD8 were the same significant biomarkers in both the TCM cold and heat patterns. The significant biomarkers for the TCM cold pattern could be EIF4A2, CCNT1, and IL7R, and the significant biomarkers for the TCM heat pattern could be PRKAA1, HSPA8, and LSM6.

This study has some limitations, including the lack of validation by PCR or Western blotting or other clinical analysis approaches for the microarray data. Because TCM pattern classification is based on an assemblage of multiple clinical manifestations, including clinical signs and symptoms and tongue and pulse appearances, the relevant biomedical networks are likely equally complex. Thus, various genes and complicated interactions were included in the results, which makes short-term validation of the results difficult. However, in this study, we use gene interaction analysis software (IPA) to represent functional clusters of different genes. Based on these findings, additional Omic (metabolics and proteomics) studies have been designed to validate the functional clusters. In addition, a more optimal study design, which includes detailed aspects such as female menstruation, diet, and activity controls, is warranted.

5. Conclusions

EIF4A2, CCNT1, and IL7R, which were related to the up-regulation of cell proliferation and the Jak-STAT cascade, were significant gene biomarkers of the TCM cold pattern of RA. PRKAA1, HSPA8, and LSM6, which were related to fatty acid metabolism and the I-κB kinase/NF-κB cascade, were significant biomarkers of the TCM heat pattern of RA. The network-based gene expression biomarkers for the TCM cold and heat patterns may be helpful for the further stratification of RA patients when deciding on interventions or clinical trials.

Authors’ Contribution

C. Lu, X. Niu, and C. Xiao have contributed equally to this work.

Acknowledgment

This paper is supported in part by the projects from the National Science Foundation of China (nos. 90709007, 30825047, and 30902003).

References

H. E. Jasin, “Mechanisms of tissue damage in rheumatoid arthritis,” in Arthritis and Allied Conditions: A Textbook of Rheumatology, W. J. Koopman, Ed., pp. 1128–1152, Lippincott Williams & Wilkins, Philadelphia, PA, USA, 2001.

M. Chang, C. M. Rowland, V. E. Garcia et al., “A large-scale rheumatoid arthritis genetic study identifies association at chromosome 9q33.2,” PLoS Genetics, vol. 4, no. 6, Article ID e1000107, 2008.

C. A. Seemayer, O. Distler, S. Kuchen, et al., “Rheumatoid arthritis: new developments in the pathogenesis with special reference to synovial fibroblasts,” Zeitschrift für Rheumatologie, vol. 60, no. 5, pp. 309–318, 2001.

M. C. Boissier, E. Assier, J. Biton, A. Denys, G. Falgarone, and N. Bessis, “Regulatory T cells (Treg) in rheumatoid arthritis,” Joint Bone Spine, vol. 76, no. 1, pp. 10–14, 2009.

A. P. Cope, H. Schulze-Koops, and M. Aringer, “The central role of T cells in rheumatoid arthritis,” Clinical and Experimental Rheumatology, vol. 25, no. 5, pp. S4–S11, 2007.

U. Schurigt, C. Pfirschke, I. M. Irmler, et al., “Interactions of T helper cells with fibroblast-like synoviocytes: up-regulation of matrix metalloproteinases by macrophage migration inhibitory factor from both Th1 and Th2 cells,” Arthritis and Rheumatism, vol. 58, no. 10, pp. 3030–3040, 2008.

G. S. Panayi, J. S. Lanchbury, and G. H. Kingsley, “The importance of the T cell in initiating and maintaining the chronic synovitis of rheumatoid arthritis,” Arthritis and Rheumatism, vol. 35, no. 7, pp. 729–735, 1992.

P. K. Gregersen, J. Silver, and R. J. Winchester, “The shared epitope hypothesis. An approach to understanding the molecular genetics of susceptibility to rheumatoid arthritis,” Arthritis and Rheumatism, vol. 30, no. 11, pp. 1205–1213, 1987.

P. Stastny, “Mixed lymphocyte cultures in rheumatoid arthritis,” Journal of Clinical Investigation, vol. 57, no. 5, pp. 1148–1157, 1976.

A. P. Lu, H. W. Jia, C. Xiao, and Q. P. Lu, “Theory of traditional chinese medicine and therapeutic method of diseases,” World Journal of Gastroenterology, vol. 10, no. 13, pp. 1854–1856, 2004.

R. Goldbach-Mansky, M. Wilson, R. Fleischmann et al., “Comparison of Tripterygium wilfordii Hook F versus sulfasalazine in the treatment of rheumatoid arthritis: a randomized trial,” Annals of Internal Medicine, vol. 151, no. 4, pp. 229–240, 2009.

C. Chai, J. Kou, D. Zhu, Y. Yan, and B. Yu, “Mice exposed to chronic intermittent hypoxia simulate clinical features of deficiency of both Qi and Yin syndrome in traditional Chinese medicine,” Evidence-Based Complementary and Alternative Medicine, vol. 2011, Article ID 356252, 7 pages, 2011.

W. Y. Jiang, “Therapeutic wisdom in traditional Chinese medicine: a perspective from modern science,” Discovery Medicine, vol. 5, pp. 455–461, 2005.

A. P. Lu, X. R. Ding, and K. J. Chen, “Current situation and progress in integrative medicine in China,” Chinese Journal of Integrative Medicine, vol. 14, no. 3, pp. 234–240, 2008.

Y. Gong, Q. Zha, L. Li, et al., “Efficacy and safety of Fufangkushen colon-coated capsule in the treatment of ulcerative colitis compared with mesalazine: a double-blinded and randomized study,” Journal of Ethnopharmacology. In press

F. S. Zhou, L. J. Hu, R. J. Wang, Z. X. Huang, and Q. Luo, “Study on clinical effect and mechanism of jianpi qingre huayu recipe,” Chinese Journal of Integrative Medicine, vol. 13, no. 1, pp. 22–26, 2007.

W. Li, B. Mao, G. Wang et al., “Effect of Tanreqing injection on treatment of acute exacerbation of chronic obstructive pulmonary disease with Chinese medicine syndrome of retention of phlegm and heat in Fei,” Chinese Journal of Integrative Medicine, vol. 16, no. 2, pp. 131–137, 2010.

Y. Yang, S. C. Wang, W. J. Bai, R. L. Li, and J. Ai, “Evaluation by survival analysis on effect of traditional Chinese medicine in treating children with respiratory syncytial viral pneumonia of phlegm-heat blocking Fei syndrome,” Chinese Journal of Integrative Medicine, vol. 15, no. 2, pp. 95–100, 2009.

Z. H. Min, Y. Zhou, and H. M. Zhang, “Effect of treatment based on syndrome differentiation by Chinese medicine on post-traumatic elbow arthritis,” Chinese Journal of Integrative Medicine, vol. 16, no. 3, pp. 264–269, 2010.

N. Zhou, Y. P. Bai, X. H. Man et al., “Effect of new pulian ointment in treating psoriasis of blood-heat syndrome: a randomized controlled trial,” Chinese Journal of Integrative Medicine, vol. 15, no. 6, pp. 409–414, 2009.

H. B. Shen, Y. J. Bai, Z. J. Huo, W. N. Li, and X. P. Tang, “Assessment of clinical effect of therapy combining disease with syndrome on rheumatoid arthritis,” Journal of Traditional Chinese Medicine, vol. 31, no. 1, pp. 39–43, 2011.

Y. F. Yang, J. Z. Ge, Y. Wu et al., “Cohort study on the effect of a combined treatment of traditional Chinese medicine and Western medicine on the relapse and metastasis of 222 patients with stage II and III colorectal cancer after radical operation,” Chinese Journal of Integrative Medicine, vol. 14, no. 4, pp. 251–256, 2008.

L. Q. Qian, B. Wang, J. Y. Niu, S. Gao, and D. Y. Zhao, “Assessment of the clinical effect of Chinese medicine therapy combined with psychological intervention for treatment of patients of peri-menopausal syndrome complicated with hyperlipidemia,” Chinese Journal of Integrative Medicine, vol. 16, no. 2, pp. 124–130, 2010.

Q. Zhang, Z. Wu, Y. Feng, and J. Shi, “Levels of sexual hormones in relation with syndrome-differentiation of TCM in patients of chronic renal failure,” Journal of Traditional Chinese Medicine, vol. 10, no. 2, pp. 132–135, 1990.

Z. Linhua, Y. Xiaoping, W. Jianming, et al., “Correlation between pathophysiology features and Chinese medical four diagnostic information in rheumatoid arthritis,” Guangzhou Zhongyiyao Daxue Xuebao, vol. 23, pp. 112–117, 2006.

Y. Q. Xie, H. Wang, Y. P. Wu, D. H. Yin, Z. S. Wang, and Y. H. Huang, “Association of APOE polymorphisms and insulin resistance with TCM syndromes in type 2 diabetes patients with macroangiopathy,” Molecular Medicine Reports, vol. 4, no. 6, pp. 1219–1223, 2011.

C. H. Hsu, M. C. Yu, C. H. Lee, T. C. Lee, and S. Y. Yang, “High eosinophil cationic protein level in asthmatic patients with “heat” Zheng,” American Journal of Chinese Medicine, vol. 31, no. 2, pp. 277–283, 2003.

H. M. Yu, H. H. Chang, S. Y. Liou, S. F. Li, M. M. Hou, and M. F. Chen, “The correlation between skin electrical conductance and the score of Qi vacuity,” American Journal of Chinese Medicine, vol. 26, no. 3-4, pp. 283–290, 1998.

C. Lu, Q. Zha, A. Chang, Y. He, and A. Lu, “Pattern differentiation in traditional Chinese medicine can help define specific indications for biomedical therapy in the treatment of rheumatoid arthritis,” Journal of Alternative and Complementary Medicine, vol. 15, no. 9, pp. 1021–1025, 2009.

S. J. Wheelan, F. Martinez Murillo, and J. D. Boeke, “The incredible shrinking world of DNA microarrays,” Molecular BioSystems, vol. 4, no. 7, pp. 726–732, 2008.

D. Diez, A. M. Wheelock, S. Goto et al., “The use of network analyses for elucidating mechanisms in cardiovascular disease,” Molecular BioSystems, vol. 6, no. 2, pp. 289–304, 2010.

E. D. Zanders, M. G. Goulden, T. C. Kennedy, and K. E. Kempsell, “Analysis of immune system gene expression in small rheumatoid arthritis biopsies using a combination of subtractive hybridization and high-density cDNA arrays,” Journal of Immunological Methods, vol. 233, no. 1-2, pp. 131–140, 2000.

A. A. Alizadeh, M. B. Elsen, R. E. Davis et al., “Distinct types of diffuse large B-cell lymphoma identified by gene expression profiling,” Nature, vol. 403, no. 6769, pp. 503–511, 2000.

T. Ma, C. Tan, H. Zhang, M. Wang, W. Ding, and S. Li, “Bridging the gap between traditional Chinese medicine and systems biology: the connection of Cold syndrome and NEI network,” Molecular BioSystems, vol. 6, no. 4, pp. 613–619, 2010.

C. Lu, C. Xiao, G. Chen et al., “Cold and heat pattern of rheumatoid arthritis in traditional Chinese medicine: distinct molecular signatures indentified by microarray expression profiles in CD4-positive T cell,” Rheumatology International, vol. 32, no. 1, pp. 61–8, 2012.

F. C. Arnett, S. M. Edworthy, D. A. Bloch et al., “The American Rheumatism Association 1987 revised criteria for the classification of rheumatoid arthritis,” Arthritis and Rheumatism, vol. 31, no. 3, pp. 315–324, 1988.

G. A. Churchill, “Fundamentals of experimental design for cDNA microarrays,” Nature Genetics, vol. 32, no. 5, pp. 490–495, 2002.

Y. H. Yang, S. Dudoit, P. Luu et al., “Normalization for cDNA microarray data: a robust composite method addressing single and multiple slide systematic variation,” Nucleic Acids Research, vol. 30, no. 4, p. e15, 2002.

W. da Huang, B. T. Sherman, and R. A. Lempicki, “Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources,” Nature Protocols, vol. 4, no. 1, pp. 44–57, 2009.

A. Real-Chicharro, I. Ruiz-Mostazo, I. Navas-Delgado et al., “Protopia: a protein-protein interaction tool,” BMC Bioinformatics, vol. 10, supplement 12, article S17, 2009.

P. Shannon, A. Markiel, O. Ozier et al., “Cytoscape: a software environment for integrated models of biomolecular interaction networks,” Genome Research, vol. 13, no. 11, pp. 2498–2504, 2003.

M. Li, J. E. Chen, J. X. Wang, B. Hu, and G. Chen, “Modifying the DPClus algorithm for identifying protein complexes based on new topological structures,” BMC Bioinformatics, vol. 9, article 398, 2008.

H. Ueda, J. M. M. Howson, L. Esposito et al., “Association of the T-cell regulatory gene CTLA4 with susceptibility to autoimmune disease,” Nature, vol. 423, no. 6939, pp. 506–511, 2003.

M. F. Liu, C. Y. Yang, J. S. Li, K. A. Lai, S. C. Chao, and H. Y. Lei, “Increased expression of down-regulatory CTLA-4 molecule on T lymphocytes from rheumatoid synovial compartment,” Scandinavian Journal of Immunology, vol. 50, no. 1, pp. 68–72, 1999.

M. Cutolo, S. Soldano, P. Montagna et al., “CTLA4-Ig interacts with cultured synovial macrophages from rheumatoid arthritis patients and downregulates cytokine production,” Arthritis Research & Therapy, vol. 11, no. 6, article R176, 2009.

http://www.genome.jp/dbget-bin/www_bget?pathway:map03050

http://www.ncbi.nlm.nih.gov/gene/904

T. L. Sung and A. P. Rice, “Effects of prostratin on Cyclin TI/P-TEFb function and the gene expression profile in primary resting CD4⁺ T cells,” Retrovirology, vol. 3, article 66, 2006.

C. H. Herrmann, R. G. Carroll, P. Wei, K. A. Jones, and A. P. Rice, “Tat-associated kinase, TAK, activity is regulated by distinct mechanisms in peripheral blood lymphocytes and promonocytic cell lines,” Journal of Virology, vol. 72, no. 12, pp. 9881–9888, 1998.

L. Y. Liou, C. H. Herrmann, and A. P. Rice, “Transient induction of cyclin T1 during human macrophage differentiation regulates human immunodeficiency virus type 1 Tat transactivation function,” Journal of Virology, vol. 76, no. 21, pp. 10579–10587, 2002.

Y. Wang and A. P. Rice, “Interleukin-10 inhibits HIV-1 LTR-directed gene expression in human macrophages through the induction of cyclin T1 proteolysis,” Virology, vol. 352, no. 2, pp. 485–492, 2006.

L. Y. Liou, C. H. Herrmann, and A. P. Rice, “Human immunodeficiency virus type 1 infection induces cyclin T1 expression in macrophages,” Journal of Virology, vol. 78, no. 15, pp. 8114–8119, 2004.

L. Y. Liou, R. E. Haaland, C. H. Herrmann, and A. P. Rice, “Cyclin T1 but not cyclin T2a is induced by a post-transcriptional mechanism in PAMP-activated monocyte-derived macrophages,” Journal of Leukocyte Biology, vol. 79, no. 2, pp. 388–396, 2006.

http://www.ncbi.nlm.nih.gov/gene/3575

A. Parker, E. S. Izmailova, J. Narang et al., “Peripheral blood expression of nuclear factor-kappab-regulated genes is associated with rheumatoid arthritis disease activity and responds differentially to anti-tumor necrosis factor-alpha versus methotrexate,” Journal of Rheumatology, vol. 34, no. 9, pp. 1817–1822, 2007.

http://www.ncbi.nlm.nih.gov/pubmed?Db=gene&Cmd=retrieve&dopt=full_report&list_uids=3603

L. R. Lard, B. O. Roep, R. E. Toes, and T. W. Huizinga, “Enhanced concentrations of interleukin 16 are associated with joint destruction in Patients with heumatoid arthritis,” Journal of Rheumatology, vol. 31, no. 1, pp. 35–39, 2004.

S. Blaschke, H. Schulz, G. Schwarz, V. Blaschke, G. A. Müller, and M. Reuss-Borst, “Interleukin 16 expression in relation to disease activity in rheumatoid arthritis,” Journal of Rheumatology, vol. 28, no. 1, pp. 12–21, 2001.

M. L. Cho, Y. O. Jung, K. W. Kim, et al., “IL-17 induces the production of IL-16 in rheumatoid arthritis,” Experimental and Molecular Medicine, vol. 40, no. 2, pp. 237–245, 2008.

K. Sudo, E. Takahashi, and Y. Nakamura, “Isolation and mapping of the human EIF4A2 gene homologous to the murine protein synthesis initiation factor 4A-II gene Eif4a2,” Cytogenetics and Cell Genetics, vol. 71, no. 4, pp. 385–388, 1995.

C. Cheyssac, C. Dina, F. Lepretre et al., “EIF4A2 is a positional candidate gene at the 3q27 locus linked to type 2 diabetes in French families,” Diabetes, vol. 55, no. 4, pp. 1171–1176, 2006.

http://www.ncbi.nlm.nih.gov/gene/820

P. Mendez-Samperio, “The human cathelicidin hCAP18/LL-37: a multifunctional peptide involved in mycobacterial infections,” Peptides, vol. 31, no. 9, pp. 1791–1798, 2010.

http://www.ncbi.nlm.nih.gov/gene/5562

http://www.ncbi.nlm.nih.gov/gene/3303

http://www.ncbi.nlm.nih.gov/gene/3312

http://www.ncbi.nlm.nih.gov/gene/11157