Introduction
Rheumatoid arthritis (RA) is an autoimmune inflammatory joint disease/disorder which affects the global population. Without treatment, RA can lead to joint destruction and severe disability. Although disease-modifying anti-rheumatic drugs and biological treatments are available for managing the condition, there is currently no cure for RA. Acupuncture is now widely accepted as a treatment for RA patients [1], however, knowledge about the underlying mechanism of the effect of acupuncture on RA is sparse.
Gene expression microarray is a platform for investigating gene activities under different physiological conditions e.g. disease status, treatment, and patient prognosis. Using an animal model of RA, Nardini et al [2] analyzed single gene expression patterns in association with acupuncture treatment of RA in rats and reported a reduced inflammatory response. In this study, a network-based analysis was used to detect and examine gene cluster profiles that are co-expressed during acupuncture treatment.
Materials and Methods
Samples
The Gene Expression Omnibus database (accession number GSE48025) was used in this study. The data for gene expression microarray was derived from 6 female Wistar rats with collagen induced arthritis (CIA) and a group of 6 female control Wistar rats. All rats were anesthetized with ether and shaved on acupunctural areas. The shaving was performed on all animals to minimize bias. Whole blood samples (5 mL) were extracted prior to acupuncture treatment in 2 CIA and 2 control rats (Time 0). Acupuncture was performed on the remaining animals by punctuation at 4 points; bilateral on the lower spine, between 2nd and 3rd lumbar vertebra (BL 23) and between the tibia and fibula on the hind legs (ST 36). The acupuncture was performed by inserting a thin needle and twirling for 20 seconds clockwise every other day. The insertion points used were similar to earlier studies on classical Chinese medicine as described by Nardini et al [2].
At the Time 1 timepoint (1 hour after acupuncture treatment), whole blood samples (5 mL) were extracted from 2 CIA rats and 2 control rats. At the Time 2 timepoint (34 days after acupuncture treatment), the 2 remaining CIA rats, and the 2 remaining control rats were anesthetized, then whole blood samples were drawn. All 12 rats (4 rats each from Time 0 to Time 2), were sacrificed after blood sampling at the end of the study.
The whole blood samples were analyzed using Affymetrix rat230 2.0 arrays which comprised of more than 31,000 probe sets, analyzing over 30,000 transcripts and variants from over 28,000 well-substantiated rat genes, which were hybridized and scanned according to the manufacturer’s recommendations. The output files were analyzed using the mas5 function from the R package affy to generate expression intensities for each probe. Details on raw microarray data pre-processing (QC and normalization) are described in Nardini et al [2].
Gene co-expression network analysis
Since multiple probes were mapped to the same gene, the expression for each gene was first summarized by calculating the mean. This resulted in 15,397 unique genes. The WGCNA R package provides functions for weighted correlation network analysis which were applied to identify the co-expression modules. Fig. 1 is the workflow for WGCNA where it is assumed that the genetic network in the system is scale free, converting their adjacency matrix into a matrix of connection of strength by raising it to a power of beta, which is chosen to satisfy the assumption of a scale free network.
In this study, the power of 10 was chosen based on a scale-free connectivity plot (Fig. 2). Pearson correlation was used for statistical analysis by assuming a linear relationship between gene expression profiles, as well as between gene expression and experimental variables. To identify the modules, the adjacency matrix was used to construct the Topology Overlap Matrix (TOM) and from this the dissimilarity matrix 1-TOM was constructed to perform hierarchical clustering. This returns a number of modules corresponding to clusters in other methods, indicated by colors (Fig. 3). Modules highly correlated by their eigengenes (1st principal component of all gene expression in a module; > 0.75) were merged to reduce the number of clusters.
Linear regression was used to explore interactions between early and late responses to treatment of RA. Each trait was focused on one by one (disease/condition: CIA = 1, control = 0; acupuncture time) and the absolute correlation between a module eigengene and the specific trait was established. The association of the individual genes to each module was then quantified.
In WGCNA, the gene significance (GS) is defined as: GSi = |cor (xi, Clinical Trait)|
Here xi is the expression of gene i. As we perform a linear regression, GS was defined as: GSi = −log10 (p value of Trait for gene i). The interpretation of GS is that the bigger the GSi, the more biologically significant the gene. As a measurement of relative importance for individual genes in a module, a module-membership (MM) was defined and calculated MM(i) = cor (xi, ME). The interpretation of this function was that the higher the absolute value of MM(i), the more important gene i was to the module. A high MM would indicate that the gene was a possible hub-gene in the module, and therefore, is of special interest. By examining the relationship between GS and MM, important hub-genes with high gene significance for the specific traits, as well as high module-membership can be identified.
Functional annotation
To investigate the biological processes for the genes in the modules associated significantly with traits, gene ontology (GO)5 and KEGG by the R package ClusterProfiler [3] was used. For the GO results, the command EnrichGO was used. All 3 ontologies were chosen to achieve a better overview of the different processes, and a p value cut off at 0.05 was selected. For better understanding of the biological processes, the GO results were submitted to AmiGO [4], which is an open access database for ontology and annotation data. For the KEGG analysis, enrich KEGG command was used, with the p value cut-off at 0.05. For further information about the pathways, KEGG results can be viewed on GenomeNet (www.genome.jp) in the KEGG PATHWAY database [5]. Hub-genes are deducted from GO results and they are visualized using the CNet plot in ClusterProfiler, the lowest p value shown in blue. Selected possible hub-genes are inputted online to retrieve possible associations with RA.
Results
RA associated gene expression networks
A cluster of 140 genes were significantly differentially expressed in the disease/condition trait (p < 0.001). The main biological process of the GO terms (Table 1) for this module was antigen processing and presentation of exogenous peptide antigen (r = −0.26, p < 0.001). In the molecular function part of the GO results, the top result was lytic vacuole and lysosome, both with p < 0.001. The cellular component GO terms showed activation in the lysosome, endocytic vesicles, and vacuoles. The top KEGG pathways (Table 2) were lysosomal, cAMP signaling, salivary secretion, and adrenergic signaling in cardiomyocytes, all with p < 0.001. Fig. 4 illustrates the gene expression networks and their linked genes in relation to CIA.
Time-dependent RA associated gene expression networks in acupuncture
WGCNA identified a module consisting of 347 genes as the best representation for the interaction between acupuncture time and disease/condition, or in other words, CIA-specific time effect of acupuncture treatment. The module had a highly significant marginal association with CIA (correlation −0.44 and p < 0.001) but also displays a time-dependent pattern with a module-RA correlation at Time 1 (−0.27; p < 0.001) and at Time 2 (−0.33; p < 0.001). The top GO (Table 3) biological processes included the hemoglobin system2 involving porphyrin-containing compound biosynthetic processes, tetrapyrrole biosynthetic processes, porphyrin-containing compound metabolic processes, and heme biosynthetic processes, all with p < 0.001. The GO results for the molecular functions were peroxidase activity, oxidoreductase activity, acting on peroxide as acceptor, oxygen carrier activity and antioxidant activity, all with p < 0.001.
The GO results were consistent with the KEGG pathway analysis results (Table 4), as porphyrin and chlorophyll metabolism (p < 0.001) were the most significant pathway. The remaining KEGG results were malaria, collecting duct acid secretion and ferroptosis, all with p < 0.001.
In Fig. 5, there were many possible hub-genes of interest for the interaction between time and disease/condition or treatment effect. The most significant hub-genes for the interaction include Slc25a39, Alas2 and Urod related to the biological processes, Hba-a1, Hbe2, Cat, Prxl2a and Prdx2 related to the molecular function, and genes Prdx2, Hba-a1, Hbe2 and Cat related to both biological processes and molecular functions.
Discussion
Network-based analysis identified a gene cluster significantly associated with RA disease/condition status. In the gene cluster, antigen presenting cells (APC) were highly implicated in the biological process part of the GO-terms correlated with RA [6]. The APC in the disease/condition pattern may come from the CIA process. It is unknown if the presence of APC is due to collagen induction of arthritis or T cell activation initiating RA and involved in the joint destruction [7]. In the cellular component part of the GO results, an activation in lysosome, endocytic vesicles and vacuoles was observed which was connected to several genes that have high significance and which could be possible hub-genes. The GO results are consistent with the KEGG results, where the top pathway is lysosomal, which is an organelle in animal cells involved in degradation of macromolecules requiring a series of processes such as endocytosis, phagocytosis, and autophagy. The molecular function part of the GO results suggested an increased activity in the transport through membranes by calcium gated channels, which correlates with the increased activity by the lysosome in various cellular components.
Although every one of the significant genes in the results could be a hub-gene to focus on, a few genes were selected. The genes Unc93b1, Ifi30 and LOC100909593 are connected to both biological processes and cellular components. The genes Anxa11 and Anxa2 are related to molecular processes and cellular components, and the genes Clec4a3 and Clec4a were related to both molecular functions and biological processes.
The Unc93b1 gene encodes a protein placed in the endoplasmic reticulum and is connected to the control of trafficking toll-like receptor involved in autoimmune diseases [8]. Listed in the rat genome database the gene LOC100909593 is a Major Histocompatibility Complex (MHC) Class II histocompatibility antigen, M alpha chain-like [9]. Its encoded protein is predicted to be part of antigen processing and presentation of antigens via MHC Class II and positive regulation of the immune response [10]. MCH Class II has a strong association with the development of CIA in animal models of RA, where MCH Class II presents the immunodominant Collagen II 256-270 peptide to specific T cells which mediate the inflammatory response [11]. This corresponds well with the module’s relation to the disease/condition status. The gene Ifi30 is also called gamma-interferon-inducible lysosomal thiol reductase. It is expressed in APC and is associated with MHC Class II by contributing to the presentation of peptides loaded on the MHC [12]. Anxa11 has been associated with the inflammatory disease sarcoidosis, which is a chronic inflammatory disorder dominantly manifested in the lungs [13]. Interestingly, the disease has been shown to co-occur with RA [14].
The interaction between time and disease i.e. the effect of acupuncture over time is conditional on or specific to the disease. This interaction effect on gene expression was of high interest in this study because it reflected changes in gene activities associated with acupuncture treatment over time, specifically in the disease group. The module significant for this interaction consisted of 347 genes. The top GO biological processes were linked to the hemoglobin system (Table 3). Among the top GO processes associated with porphyrin, the component in the red blood cells which contains the heme group and involved in heme biosynthesis. This could indicate blood formation and improved oxygen transport in the body. On the other hand, this could also be an effect of the CIA-treatment which is causing increased blood flow in the skin in the CIA rats [2]. The GO results are consistent with the KEGG pathway analysis results, where porphyrin metabolism is the most significant KEGG pathway. Porphyrin derivatives have also been suggested as a photodynamic treatment where the porphyrin derivatives are used as a photosensitizer together with oxygen molecules (which produces a single oxygen and possible reactive oxygen species) and external light to interfere with cellular signaling pathways, This leads to necrosis or apoptosis in the target tissues reducing inflammation of the tissue [15].
In the KEGG results, ferroptosis, a regulated form of cell death that is present in some neurodegenerative diseases and tissue damage, may be associated with RA. Ferroptosis can be induced therapeutically by sulfasalazine, which is a drug used to treat RA [16], but as the CIA rats were not given that drug, the acupuncture treatment may have activated some of the same pathways induced by the drug. As a reaction of the metabolism of porphyrin, it could be a very interesting target for RA researchers [17]. Ferroptosis requires iron, while the latter is associated with activated hemes. Ferroptosis is the accumulation of reactive oxygen species [16], which could be interacting with the porphyrin and be activated by an external light source. This could lead to a decrease in inflammation, but further studies need to be conducted for verification.
This study was limited by the sample size and design. As an association study in an animal model, the results presented here represent correlations rather than causal relationships. Well-defined randomized controlled trials are needed to establish deterministic effects of the treatment.
Conclusion
This network-based analysis detected significant gene clusters enriched for functional pathways and biological processes implicated in an RA model in association with acupuncture treatment. The identified gene expression networks and their hub-genes could help with the understanding of mechanisms involved in the pathogenesis of RA, as well as in the effects of acupuncture treatment of RA.
PDF Links
PubReader
ePub Link
XML Download
Full text via DOI
Download Citation
CrossRef TDM
Print
