Systematic analysis of inflammation and pain pathways in a mouse model of gout.

Gout is a prevalent and painful inflammatory arthritis, and its global burden continues to rise. Intense pain induced by gout attacks is a major complication of gout. However, systematic studies of gout inflammation and pain are lacking. Using a monosodium urate (MSU) crystal-induced gout model, we performed genome-wide transcriptome analysis of the inflamed ankle joint, dorsal root ganglion (DRG), and spinal cord of gouty mice. Our results revealed important transcriptional changes, including highly elevated inflammation and broad activation of immune pathways in both the joint and the nervous system, in gouty mice. Integrated analysis showed that there was a remarkable overlap between our RNAseq and human genome-wide association study (GWAS) of gout; for example, the risk gene, stanniocalcin-1 (STC1) showed significant upregulation in all three tissues. Interestingly, when compared to the transcriptomes of human osteoarthritis (OA) and rheumatoid arthritis (RA) joint tissues, we identified significant upregulation of cAMP/cyclic nucleotide-mediated signaling shared between gouty mice and human OA with high knee pain, which may provide excellent drug targets to relieve gout pain. Furthermore, we investigated the common and distinct transcriptomic features of gouty, inflammatory pain, and neuropathic pain mouse models in their DRG and spinal cord tissues. Moreover, we discovered distinct sets of genes with significant differential alternative splicing or differential transcript usage in each tissue, which were largely not detected by conventional differential gene expression analysis approaches. Based on these results, our study provided a more accurate and comprehensive depiction of transcriptomic alterations related to gout inflammation and pain.

Introduction

Gout is the most common form of inflammatory arthritis with an approximately 2–4% prevalence worldwide, and it is characterized by the deposition of monosodium urate (MSU) crystals in and around joints. The natural history of gout is typically divided into four stages as follows: asymptomatic hyperuricemia (high uric acid levels in the blood), acute gout, interval period, and chronic tophaceous gout. Gout flares can start suddenly and be extremely painful for days or weeks, followed by a long period of time with no symptoms, which is known as remission. Etiological studies of gout can be traced back to the eighteenth century. It has been well established that hyperuricemia and other influencing factors, including reduced temperatures, lower pH, high concentration of sodium ions, and abnormal connective tissues, lead to crystals of uric acid to form and settle in the joints, eventually causing gouty arthritis.[2-4]

Half a century ago, multiple studies have shown that the injection of MSU crystals into the joints of normal men and gouty patients successfully induced acute inflammatory responses that closely resembled gouty attacks. Since then, several animal models have been generated to mimic human gouty arthritis by direct injection of MSU-crystals into different anatomical structures, including ankles, knees, and paws.[2,5-7]

During the past few decades, gout-associated immune pathways and inflammatory responses have been gradually elucidated. A previous study has shown that MSU crystals are proinflammatory by activating the NALP3 inflammasome, resulting in the production of active IL-1β and IL-18. Following the activation of the NLRP3 inflammasome and the release of cytokines (especially IL-1β) and chemokines, innate immune cells, such as monocytes and neutrophils, are recruited to the inflammatory sites. Although MSU crystals activate NLRP3, which cleaves pro-IL-1β into mature IL-1β, it remained unclear how pro-IL-1β was induced in gout. Subsequent studies have shown that myeloid-related proteins S100A8 and S100A9, which are endogenous TLR4 ligands, secreted by MSU crystal-activated phagocytes, induced pro-IL-1β in a TLR4-dependent manner. However, large deposits of tophi alone may not be sufficient to elicit gouty attacks. Dietary factors, such as high-purine foods and alcohol, can trigger gout flares. Joosten et al. reported that the interaction of free fatty acids with TLR2 synergized with MSU crystals to induce the release of a large amount of IL-1β. Another study has suggested that through an NALP3 inflammasome independent pathway, the AIM2 inflammasome is involved in the full development of acute gouty inflammation. Additionally, serum levels of IL-33 in gout patients are significantly higher than those in healthy controls. In mouse gout model, IL-33 and its receptor, ST2, promote neutrophil-dependent ROS production and activate TRPA1 in the DRG.

To date, most studies have focused on the acute inflammatory responses induced by exogenous MSU crystals. However, little is known about the cellular and molecular mechanisms of gout pain. As an excruciatingly painful form of arthritis, gout attacks can prevent one from walking or going to work. If left untreated, gout can develop into chronic arthritis with chronic pain. Furthermore, previous transcriptomic analyses of gout animal models are limited to inflammatory tissues or joints.[5,7] Therefore, it remains unclear how the DRG and spinal cord are affected to produce highly intensified and excruciating pain in gout. Taking all these factors into consideration, we decided to investigate gout from the perspectives of inflammation and pain. In the present work, we sought to systematically investigate the transcriptional changes by genome-wide expression profiling of the ankle joints, DRG, and spinal cord after MSU-induced acute inflammation in the ankle joints of mice. We then integrated our RNAseq data with large-scale GWAS of gout to identify the top candidate genes. To gain new insights into the molecular mechanisms of different types of arthritis, we performed comparative RNAseq analysis of joint tissue samples among the mouse gout model, human osteoarthritis (OA), and human rheumatoid arthritis (RA). Moreover, we integrated our RNAseq data of the DRG and spinal cord with public RNAseq datasets of complete Freund's adjuvant (CFA) and spared nerve injury (SNI) pain models of mice to discover common dysregulated genes and pathways in pain. In addition, by comparison with RNAseq of human OA with high knee pain, we aimed to identify genes for arthritis-related severe joint pain.

RNAseq based alternative splicing analysis provides robust detection of disease-associated splice events. Alternative splicing is known to be involved in regulating normal physiological functions as well as diseases. Aberrant alternative splicing events are associated with cancer, neurological and psychiatric diseases, inflammation, immune and metabolic disorders. Numerous alternatively spliced isoforms of ion channels and GPCRs in nociceptors have been recognized. A recent study has shown that alternative splicing of the Nrcam gene in the DRG contributed to neuropathic pain genesis and targeting Nrcam alternative splicing effectively attenuated pain hypersensitivity. We observed widespread and distinct alternative splicing events and differential transcript usage (DTU) in the ankle joints, DRG, and spinal cord, suggesting that tissue-specific alternative splicing may contribute to the tissue-specific manifestation of gout.

Materials and methods

Animals

C57BL/6 mice were purchased from the Laboratory Animal Center of Wenzhou Medical University. Mice were housed with a 12-h light/12-h dark cycle and have free access to rodent lab chow and water. Twelve-to sixteen-week-old mice were selected for the assay. The handling of mice and experimental procedures in this study were performed in accordance with the requirements of the Institutional Animal Care and Use Committee of Wenzhou Medical University.

MSU-induced gout mouse model

For the preparation of MSU crystals: 1.0 g uric acid (Sigma-Aldrich) was dissolved in 200 mL boiling distilled water containing 6.0 mL of 1 M NaOH. The pH of the solution was adjusted to 7.2 with HCl until the solution temperature dropped to room temperature. The solution was stirred on a magnetic stirrer for 24 h and then crystals formed. Crystals were then washed with ddH20, dry at 55°C for 24 h and autoclaved. Sterile MSU crystals were suspended in phosphate-buffered saline (PBS, 1x) at a concentration of 20 mg/mL. All reagents were prepared under pyrogen-free conditions.

In this model, the right ankle joint of mice in the MSU group received an intra-articular injection of 100 μL MSU suspension (2 mg MSU crystals in total). For the control group, 100 μL PBS solution was injected into the right ankle joint. Eight hours after injections, the degree of joint swelling was measured with an electronic vernier caliper and photographed. Mechanical hyperalgesia was evaluated using Von Frey filament tests.

RNAseq library preparation and RNAseq data analysis

Eight hours post-injections of MSU crystals or 1 x PBS, mice were anesthetized with isoflurane and decapitated. Ankle joints, DRG, and spinal cord were collected and stored in RNAlater (Thermo Fisher Scientific). Total RNAs were extracted using TRIzol reagent (Thermo Fisher Scientific). RNAseq libraries were prepared for sequencing using a VAHTS Universal V6 RNAseq Library Prep Kit (Vazyme) for Illumina. The library was sequenced on an Illumina HiSeq platform. Clean sequencing reads were aligned to the mouse genome (mm39) using the HISAT2 alignment program, and a >95% alignment rate was achieved for each sample. Differential expression analyses were performed with the DESeq2 R package. For each type of tissue, genes with total counts less than 50 summarized for all 6 samples were filtered out first. For our own RNAseq and public RNAseq datasets, differential expression analysis was performed with the DESeq2 R package. Genes were considered differentially expressed when they had a | log2(fold change) | > log21.5, and p-value < 0.05. Functional enrichment analysis was conducted using the clusterProfiler package. Venn diagrams were generated by the VennDiagram package.

Gene regulatory network construction and visualization

Transcriptional regulatory interactions between transcription factors and their target genes were determined according to the TRRUST database. The gene regulatory network was inferred using GENIE3 and visualized using Gephi (0.9.2).

Alternative splicing analysis

Differential alternative splicing and differential transcript usage (DTU) were analyzed by SUPPA2. Thresholds of statistical significance for |ΔPSI| (PSI: percent spliced-in) and p-values were set at 0.1 and 0.05, respectively. Gene Ontology (GO) analysis for genes with significantly differential alternative splicing and DTU was performed using the clusterProfiler package. Venn diagrams were generated by the VennDiagram package.

Statistical analysis

All data are expressed as the mean ± SEM. Data were analyzed by a two-tailed Student’s t-test between two groups (experimental and control). The criterion for statistical significance was p < 0.05.

Data availability statement

All data used to support the findings of this study are available within the paper and supplementary materials. Our own 18 RNAseq datasets have been deposited in the Gene Expression Omnibus (GEO) repository under access number GSE190138. RNAseq datasets for 18 normal and 20 OA human knee cartilage tissues were downloaded from the GSE114007 dataset. RNAseq datasets of 5 OA patients with high knee pain and 5 OA patients with low knee pain were downloaded from the GSE99662 dataset. RNAseq datasets of CFA and SNI pain mouse models were obtained from the GSE111216 dataset. Human RA single-cell RNAseq data were downloaded from ImmPort (https://www.immport.org/shared/study/SDY998).

Results

Transcriptomics of the mouse gout model

To study the pathophysiology of gouty inflammation and pain at the molecular level, we used an established mouse model of gout in which MSU crystals were injected into the ankle joints of C57BL/6 male mice as described. The control group received 100 μL of 1 x PBS injection. Eight hours after MSU injection, mice developed obvious joint edema and mechanical allodynia (Figure 1(a)-(c)). We next performed RNAseq from the MSU-induced acute inflammation model. For each injected mouse, total RNAs were extracted from three different types of tissues, namely, the ankle joint, DRG (L3-L5), and spinal cord. The quality of RNAseq was assessed by mapping rates on the mouse genome and unsupervised clustering of samples based on gene expression. On average, the alignment rate was >95% for each sample when the reads were mapped to the reference genome. Principal component analysis (PCA) effectively distinguished the samples from different tissues under control and gout conditions (Figure 1(d)). In particular, the three types of tissues presented clear separations with only two principal components. For each type of tissue, samples were well separated between the control and gout models (Figure 1(d)).

Figure 1.

MSU-induced mouse gout model and overview of transcriptomes of the injected ankle joint, DRG, and spinal cord. a. Representative photos showing mouse ankle joints injected with 1 x PBS or MSU. The pictures were taken 8 h after injections. b. Von Frey filament tests showing mechanic allodynia of MSU-injected mice (n = 6 mice/group). c. Swelling degree of MSU-injected mice as indicated by increased ankle diameter (n = 6 mice/group). d. Tissue-based principal component analyses (PCA) of mouse RNAseq under gout condition showed clustering and clear separation for three different types of tissues and between case and control samples of each tissue (n = 3 mice/group).

Differential gene expression was analyzed by DESeq2. Genes with a | log2(fold change) | > log21.5 and p-value < 0.05 were selected as DEGs. In total, we identified 5439 DEGs (2397 upregulated and 3042 downregulated) from the ankle joint, 2757 DEGs (1621 upregulated and 1136 downregulated) from the DRG, and 1157 DEGs (771 upregulated and 386 downregulated) from the spinal cord (Figure 2(a) and 2(b)). We then performed functional enrichment analysis for each tissue, which demonstrated that upregulated DEGs in all three tissues were mainly associated with the inflammatory response and immune pathways (Figure 2(c)). However, the top-ranked GO terms for downregulated DEGs were different (Figure 2(c)). Furthermore, the Venn diagrams in Figure 2(a) and 2(b) show the overlapping upregulated and downregulated DEGs from the ankle joint, DRG, and spinal cord. There were 343 shared upregulated DEGs and only 37 shared downregulated DEGs (Supplemental Table 1). Importantly, the 343 shared upregulated DEGs were also significantly enriched in the inflammatory response and immune pathways (Figure 2(c)), which suggested that severe inflammation was triggered in both the ankle joint and the nervous system.

Figure 2.

Transcriptomic analysis of the injected ankle joint, DRG, and spinal cord. a. Intersection of upregulated DEGs between the ankle joint, DRG, and spinal cord of gouty mice. b. Intersection of downregulated DEGs between the ankle joint, DRG, and spinal cord of gouty mice. c. GO functional enrichment analysis of up- and downregulated DEGs from the three tissues. d. Gene regulatory network for upregulated inflammatory and immune DEGs in the ankle joint.

To better understand the transcriptional regulation of gout inflammation, we further constructed a gene regulatory network for the upregulated DEGs of inflammatory and immune responses in the ankle joints. Consequently, 28 TFs were first extracted and their targeted genes were inferred based on the TRRUST database. In total, 107 target genes were selected by intersecting inferred targeted genes with the upregulated DEGs related to inflammatory and immune responses in the ankle joints. We applied GENIE3, a top-ranked network inference algorithm, to build the gene regulatory network, and we visualized the network by Gephi (Figure 2(d)). The size of the circle represents the number of target genes. Many key inflammatory transcription factors, such as Irf1/2/7, Stat1/3/5a, Myc, and Nr4a1, are the most important regulators in the network that directly regulate a broad range of inflammatory mediators including cytokines, chemokines, and matrix metalloproteinases (MMPs) (Figure 2(d)).

Table 1 shows the top 25 DEGs for the ankle joint, DRG, and spinal cord in the mouse gout model based on log2FC. In addition to significantly upregulated chemokines and cytokines, serum amyloid A (SAA) family proteins, a major type of acute phase proteins and the most sensitive indicators in response to inflammation and tissue injury, were highly expressed in all three tissues. In the ankle joint, the expression of three SAA family proteins, SAA1/2/3 expression was increased approximately 500-fold. In the DRG, SAA1/2/3 were the most upregulated DEGs and their levels increased as much as in the ankle joint. Similarly, the gene expression levels of SAA1/2/3 significantly increased in the spinal cord. SAA proteins have strong proinflammatory properties and induce a wide range of cytokines and chemokines. Such increases were not observed in the DRG and spinal cord of mouse CFA and SNI pain models, or in human OA and RA joint tissues analyzed in this study, indicating that there are important differences between gout and other related diseases.

Table 1.

List of top 25 DEGs for the ankle joint, DRG and spinal cord in the mouse gout model based on log2FC.

Ankle joint	log2FC	DRG	log2FC	SC	log2FC
Csf3	11.5211	Saa1	10.3678	Csf3	9.6181
Cxcl3	10.2888	Saa2	9.0135	Cxcl1	9.1630
Il6	9.8986	Saa3	8.4723	Selp	7.9485
Saa3	9.8221	Csf3	8.1977	Saa3	7.7465
Cxcl2	9.7280	Odf3l2	8.1241	Serpina3f	7.0725
Saa2	9.6969	Tnfsf11	8.0714	Mefv	7.0175
Saa1	8.7028	Elobl	7.1844	Sele	6.3845
Acod1	8.3274	Sele	7.1192	Il6	6.2670
Ccl4	8.1948	Fgb	7.0158	Slc10a6	6.2058
Nos2	7.2523	Cxcl1	6.8378	Cxcl2	6.1489
Cxcl1	7.2176	Catsper4	6.6992	Il1rn	5.6162
Klk1b1	7.1834	Tmem52	6.3292	Saa1	5.5035
Ccl3	7.0601	Myhas	6.3137	Plin4	5.3598
Sult6b2	6.9770	Prok2	6.2708	F10	5.1351
Ptx3	6.7073	Dupd1	6.2434	Ms4a8a	5.0921
Mbl1	6.6097	Lrrc30	6.0803	Hcar2	5.0827
Cxcl5	6.5312	Selp	6.0472	Lcn2	4.8363
Slc22a6	6.3333	Ptx3	5.9484	Saa2	4.7556
Il1b	6.2989	Myl1	5.7829	Cd244a	4.5848
Mir7678	6.2929	Plin4	5.6349	Ms4a6d	4.3855
Slfn4	6.1824	Ampd1	5.5826	Plaur	4.3416
Ms4a20	6.0637	Abra	5.5299	Ms4a4a	4.0995
Dmbt1	5.9288	Cox6a2	5.4676	Ms4a4c	4.0984
Ptgs2	5.8963	Myf6	5.4250	Ptgs2	4.0929
Clec4e	5.8540	Mylpf	5.3839	Clec4d	4.0369

Relevance of DEGs in the DRG and spinal cord of gouty mice to mouse pain models

Pain is first detected by DRG neurons and then transmitted to the dorsal horn of the spinal cord. The spinal dorsal horn is a major site of nociceptive signal processing and relay. From there, pain information is sent to the higher brain centers including the thalamus and cortex along the ascending pathway. Therefore, to better understand how excruciating pain is produced and transmitted with gout flares, it is crucial to investigate the transcriptional changes in the DRG and spinal cord in the gout model.

Complete Freund’s adjuvant (CFA)-induced inflammatory pain and spared nerve injury (SNI) induced neuropathic pain rodent models have been well established and widely used. Several transcriptomic studies have revealed transcriptional alterations in different pain models.[22,24,25] Here, we conducted transcriptomic comparisons of the DRG and spinal cord among the gout, CFA and SNI models. Transcriptome profiles of the CFA and SNI models were retrieved from a previous article. Overall, there were 155 overlapping DEGs in the DRG of the gout, CFA and SNI models (126 upregulated and 29 downregulated, Supplemental Table 2) (Figure 3(a)). The most enriched GO biological processes for the 129 upregulated DEGs included regulation of the inflammatory response, positive regulation of TNF and cytokine production (Figure 3(c)). For the spinal cord, 136 DEGs were shared among the three models (114 upregulated and 22 downregulated, Supplemental Table 3) (Figure 3(b)). The 111 upregulated DEGs were strongly associated with granulocyte, leukocyte, and neutrophil chemotaxis as well as, leukocyte migration and activation (Figure 3(d)). Our results indicated that upregulation of inflammatory and immune responses is a common feature in the DRG and spinal cord in the gout, CFA and SNI mouse models.

Figure 3.

Comparative transcriptomic analysis of the DRG and spinal cord between gout, CFA and SNI pain mouse models. a. Intersection of DEGs in the DRG between gout, CFA and SNI pain mouse models. b. Intersection of DEGs in the spinal cord between gout, CFA and SNI pain mouse models. c. GO enrichment analysis for shared DEGs of the DRG. d. GO enrichment analysis for shared DEGs of the spinal cord.

Integrative analysis of gouty mouse RNAseq and human gout GWAS

Genome-wide association studies (GWAS) have identified many susceptibility loci associated with serum urate concentrations and gout.[26-28] In 2012, the Global Urate Genetics Consortium reported a large-scale GWAS from >140,000 individuals of European ancestry and identified 28 genome-wide significant loci. Additionally, a GWAS conducted on a Japanese population has discovered 10 loci. We compared our DEGs with the human GWAS top hits and revealed many overlaps, which mostly showed downregulation in our RNAseq data, including ATXN2, CNIH2, CUX2, HLF, KCNQ1, NIPAL1, NRXN2, PRKAG2, and SLC2A9. Several overlaps were upregulated in our study, including IGF1R, INHBB, STC1, UBE2Q2, and VEGFA. In total, among these candidate genes identified by GWAS, 14 of 37 (37.84%) were differentially expressed in at least one type of tissue of the mouse gout model. Notably, only stanniocalcin-1 (STC1) was significantly upregulated in all three types of tissues, namely, the ankle joint, DRG, and spinal cord (Table 2). STC1 is a secreted glycoprotein, implicated in several pathologies including cancer, retinal degeneration, and inflammation. It has been shown that tumor STC1 inhibits phagocytosis and promotes tumor immune evasion and resistance. It will be interesting to determine its role during gout development.

Table 2.

Overlaps between the gouty mice RNAseq and human GWAS of gout, corresponding statistics of DEGs are bold, underlined and italic.

	Ankle joint		DRG		SC
Symbol	log2FC	p-value	log2FC	p-value	log2FC	p-value
ATXN2	−0.03389	0.81245	−0.63983	0.00000	−0.02006	0.82734
CNIH2	−0.96389	0.02244	−0.21701	0.12839	0.13911	0.31001
CUX2	0.01525	0.97320	−0.70206	0.00000	−0.08247	0.44834
HLF	−1.86339	0.00000	−0.01119	0.98194	−0.08137	0.45980
IGF1R	0.77819	0.00000	−0.07880	0.44310	0.12685	0.23950
INHBB	0.94277	0.00000	−0.14309	0.36845	−0.17877	0.28079
KCNQ1	−1.42708	0.00067	−0.07852	0.80479	−0.31812	0.37203
NIPAL1	−1.52108	0.00000	0.27202	0.35190	−0.19287	0.65037
NRXN2	−1.71093	0.00002	−1.07030	0.00000	0.16269	0.03693
PRKAG2	−0.58726	0.00001	−0.29779	0.00042	0.08028	0.42038
SLC2A9	−1.05144	0.00012	−0.56304	0.02822	0.06642	0.88348
STC1	1.52280	0.00000	2.75347	0.00002	0.83847	0.00560
UBE2Q2	1.04378	0.00000	−0.10197	0.36463	−0.05931	0.53433
VEGFA	0.85048	0.00000	−0.06485	0.64870	−0.56125	0.00000

The GWAS loci not shared with DEGs of the gouty mice RNAseq are listed as below.

A1CF, ABCG2, ACVR1BACVRL1, B3GNT4, BAZ1B, BCAS3, FAM35A, GCKR, HIST1H2BF-HIST1H4E, HNF4G, INHBC, MAF, NFAT5, PDZK1, RFX3, RREB1, SFMBT1, SLC16A9, SLC17A1, SLC22A11, SLC22A12, TMEM171, TRIM46.

The above comparisons between the DEGs identified in our study and human GWAS support the contention that some of the DEGs uncovered in the mouse gout model have biological relevance to human gout pathogenesis. The 14 genes simultaneously detected as DEGs and by GWAS have a higher probability of being involved in gout pathogenesis.

Comparative analysis between gouty mouse RNAseq and human OA transcriptomes

Osteoarthritis (OA) is the most common form of arthritis that can cause severe joint pain, swelling, and stiffness, affecting more than 10% of the population aged 60 years or older. OA is a chronic “wear and tear” disease characterized by the deterioration of cartilage in joints. A recent RNAseq study has been performed on 18 normal and 20 OA human articular cartilage samples in the knee joint to identify dysregulated genes and pathways. We retrieved the RNAseq data from NCBI GEO (access no. GSE114007) and analyzed them with DESeq2. A Venn diagram showed 1450 overlapping DEGs between the ankle joint of the gout model and human OA (Figure 4(a)). Overrepresentation analysis was performed using GO databases for overlapping DEGs. The results showed that upregulated overlapping DEGs were mainly associated with extracellular matrix organization, ossification, and neutrophil activation (Figure 4(d)).

Figure 4.

Comparative transcriptomic analysis of the joint tissues between gouty mice, human OA and RA. a. Intersection of DEGs in the joint tissues between gouty mice and human OA. b. Intersection of DEGs in the joint tissues between gouty mice and human OA with high knee pain. c. Intersection of DEGs in the joint tissues between gouty mice and human RA. d. GO enrichment analysis for shared DEGs of the joint tissues between gouty mice and human OA. e. GO enrichment analysis for shared DEGs of the joint tissues between gouty mice and human OA with high knee pain. f. GO enrichment analysis for shared DEGs of the joint tissues between gouty mice and human RA. g. Cell-type specific GO enrichment analysis for shared DEGs of the joint tissues between gouty mice and human RA.

Severe joint pain is one of the leading causes of disability and morbidity in arthritis. The etiology of joint pain in OA and RA is still unclear. Transcriptomic comparison of OA patients with high knee pain and low knee pain has identified many candidate genes associated with joint pain in OA knee synovial tissues. A total of 246 of these OA joint pain genes were overlapped with the ankle joint DEGs in the present study, and the overlapping DEGs are mostly upregulated (Figure 4(b)). Importantly, the upregulated DEGs are involved in cAMP/cyclic nucleotide-mediated signaling and GPCR signaling pathways (Figure 4(e) and Table 3). There is growing evidence indicating that pain processing is dependent upon activation of the cAMP signaling pathway in a dose-dependent manner.[32-34] Inflammatory mediators reduce the threshold for nociceptor neurons to fire action potentials and lead to pain sensitivity. Inflammatory pain behaviors were greatly reduced or abolished in adenylyl cyclase knock-out mice. Therefore, cAMP-mediated signaling and GPCR signaling play a key role in nociceptor sensitization caused by inflammatory and immune mediators. The overlaps identified in the present study are promising targets of drug screening for effective pain relief (Table 3).

Table 3.

List of shared upregulated DEGs in the ankle joint between gouty mice and human OA with high knee pain, involved in cAMP/cyclic nucleotide-mediated signaling and GPCR signaling pathways.

Symbol	Gene name	log2FC	p-value
GAL	Galanin and GMAP prepropeptide	2.7442	0.0010
MRAP	Melanocortin 2 receptor accessory protein	1.9223	0.0177
RASD1	Ras related dexamethasone induced 1	1.5736	0.0000
RAMP1	Receptor activity modifying protein 1	1.4898	0.0315
PRKAR2B	Protein kinase cAMP-dependent type II regulatory subunit beta	1.3688	0.0382
ATP1A2	ATPase Na+/K+ transporting subunit alpha 2	1.2874	0.0036
CNR1	Cannabinoid receptor 1	1.2838	0.0103
PTGER3	Prostaglandin E receptor 3	1.2799	0.0121
ADCY5	Adenylate cyclase 5	1.1300	0.0055
PDE4D	Phosphodiesterase 4D	1.0775	0.0000
PTGFR	Prostaglandin F receptor	1.0622	0.0005
CASQ2	Calsequestrin 2	1.0067	0.0251
NPR1	Natriuretic peptide receptor 1	1.0042	0.0169
PLN	Phospholamban	0.9696	0.0056
GNAZ	G protein subunit alpha z	0.9567	0.0290
ADRA1A	Adrenoceptor alpha 1A	0.9532	0.0117
DMD	Dystrophin	0.9209	0.0000
ADCY3	Adenylate cyclase 3	0.8855	0.0000
KSR1	Kinase suppressor of ras 1	0.8558	0.0000
S1PR3	Sphingosine-1-phosphate receptor 3	0.7286	0.0309
GNAL	G protein subunit alpha L	0.7149	0.0002
ADGRL3	Adhesion G protein-coupled receptor L3	0.7018	0.0413
PALM	Paralemmin	0.6707	0.0123
NFATC4	Nuclear factor of activated T cells 4	0.6698	0.0019
ACKR3	Atypical chemokine receptor 3	0.6416	0.0150
PTGIR	Prostaglandin I2 receptor	0.6050	0.0071

Comparative analysis between gouty mouse RNAseq and human RA transcriptomes

Rheumatoid arthritis (RA) is a chronic autoimmune and inflammatory disease that causes swelling, stiffness, pain, and loss of function primarily in the joints. The prevalence of RA is relatively constant, accounting for approximately 1% of the world population. RA is characterized by systemic inflammation and persistent synovitis, leading to serious complications in many other tissues as well. GWASs and transcriptomic studies have significantly improved the understanding of RA pathogenesis. Recently, synovial T cells, B cells, monocytes, and fibroblasts from 36 RA and 15 OA patients have been analyzed by utilizing an integrated approach with single-cell RNAseq, bulk RNAseq, flow cytometry, and mass cytometry. Characterized by a high abundance of CD4+ T cells, CD8+ T cells, and B cells, leukocyte-rich RA had significantly higher Krenn inflammation scores than leukocyte-poor RA or OA samples, and subsequently, single-cell RNAseq analysis defined 18 distinct cell subpopulations for the four types of cells. We selected enriched genes for leukocyte-rich RA samples, which had 143 overlapping genes with the ankle joint DEGs of gouty mice (Figure 4(c)). GO analysis of all 143 genes demonstrated that the top enriched biological processes included T-cell activation and positive regulation of cell-cell adhesion for upregulated genes; however, downregulated genes are significantly associated with extracellular matrix organization (Figure 4(f)). Cell type enrichment analysis demonstrated that distinct GO pathways were statistically overrepresented in each cell type (Figure 4(g)), as follow: the “immune response-activating cell surface receptor signaling pathway” and “immune response-activating signal transduction” GO terms were enriched in B cells; the “extracellular matrix organization” GO term was enriched in fibroblasts; the “negative regulation of immune system process” GO term was enriched in monocytes; and the “positive regulation of cell-cell adhesion” and “positive regulation of cell activation” GO terms were enriched in T cells. Together, the comparative studies between gouty mouse RNAseq and human RA single-cell RNAseq revealed shared inflammatory/immune response genes and pathways for specific cell types in the joint tissues.

Differential alternative splicing under gout conditions

Alternative splicing (AS) of pre-mRNA is an essential mechanism to generate genomic and proteomic diversity. In humans, more than 95% of multiexon genes undergo alternative splicing to encode two or more splice isoforms. Compared to gene expression profiling alone, alternative splicing analysis of RNAseq provided robust detection of disease-associated splice events and allowed us to analyze genomic diversity in an unprecedented depth.

In the present study, we utilized our RNAseq data to incorporate splice variant analysis. In total, 49,270, 51,976, and 51,521 alternative splicing events were identified in the ankle joint, DRG, and spinal cord, respectively (Figure 5(a)). The distribution of seven different types of alternative splicing events in the three tissues was similar with alternative first exons (AF) at the highest frequency (>32%), followed by skipping exon (SE) (>26%), alternative 3′ splice sites (A3) (∼13%), and alternative 5′ splice site (A5) (∼11%) (Figure 5(a)). To define genes differentially alternatively spliced in response to gout conditions, differential alternative splicing events were detected by SUPPA2 using cutoffs of | ΔPSI | > 0.1 and p-value < 0.05. We identified 613, 653, and 371 differentially alternatively spliced genes in the ankle joint, DRG, and spinal cord, respectively (Figure 5(b)-5(d) and Supplemental Table 6). We then examined the relationships between DEGs and differentially alternatively spliced genes. The number of DEGs was greater than that of differentially alternatively spliced genes. In the ankle joint, 239 out of 613 (38.99%) differentially alternatively spliced genes are also DEGs. However, for the DRG and spinal cord, only 111 out of 653 (17.00%) and 33 out of 371 (8.89%) differentially alternatively spliced genes are DEGs, respectively (Figure 5(b)-5(d)).

Figure 5.

Alternative splicing analysis for gouty mice. a. Overview of the seven different types of alternative splicing and their frequencies in the three types of tissues. b. Intersections between DEGs, differentially alternatively spliced (DAS) genes, and differential transcript usage (DTU) genes for the ankle joint of gouty mice. c. Intersections between DEGs, differentially alternatively spliced (DAS) genes, and differential transcript usage (DTU) genes for the DRG of gouty mice. d. Intersections between DEGs, differentially alternatively spliced (DAS) genes, and differential transcript usage (DTU) genes for the spinal cord of gouty mice. e. GO enrichment analysis for differentially alternatively spliced (DAS) genes in the ankle joint, DRG, and spinal cord of gouty mice. f. GO enrichment analysis for differential transcript usage (DTU) genes in the ankle joint, DRG, and spinal cord of gouty mice.

Functional enrichment analysis showed that the major overrepresented GO terms for differentially alternatively spliced genes in the ankle joint were “mRNA processing” and “myeloid cell differentiation”. The enriched GO biological processes for differentially alternatively spliced genes in the DRG included “mitochondrion organization”, “vesicle organization”, and “organelle localization”. However, no significant GO terms were identified for differentially alternatively spliced genes in the spinal cord (Figure 5(e)). Overall, these results showed that the ankle joint, DRG, and spinal cord underwent different alternative splicing changes under gout conditions, which may be related to different aspects of gout pathogenesis.

Differential transcript usage under gout conditions

Differential gene expression has been the primary focus of most transcriptomic studies. However, isoform diversity information is lost in conventional DEG analysis. Recent studies have demonstrated that specific transcript usage profiles are associated with cancer and neurological diseases.[40,41] Differential transcript usage (DTU) measures and compares how much each transcript contributes to the total expression of the gene between conditions. In the present study, we performed a genome-wide analysis of DTU in gouty mice. A gene exhibiting at least one DTU event with statistical significance was defined as a DTU gene. As a result, 484, 527, and 255 DTU genes were detected in the ankle joint, DRG, and spinal cord, respectively. Interestingly, the majority of DTU genes were not differentially expressed, and DTU genes largely did not overlap with differentially alternatively spliced genes as well (Figure 5(b)-5(d)). Importantly, DTU genes from different tissues were enriched in different functional pathways. The most enriched biological processes for DTU genes from the ankle joint are “extracellular matrix organization”, “cell chemotaxis”, “positive regulation of cytokine production”, and “myeloid leukocyte activation”. In contrast, the most enriched GO terms for DTU genes from DRG are “mitochondrion organization”, “positive regulation of protein polymerization”, and “purine ribonucleotide metabolic process”. Moreover, no enriched GO terms were found for DTU genes from the spinal cord (Figure 5(f)). Thus, the DTU analysis provided additional insights into the transcriptomic landscape in gout.

Discussion

The present study is the first genome-wide transcriptomic analysis of the inflamed ankle joint, DRG, and spinal cord in a mouse model of gout. Sudden onset of severe inflammation and intense pain are characteristic features of gout arthritis. Although gout inflammation has been extensively studied, systematic analysis of transcriptomic alterations in gout is lacking and the mechanisms involved in gout pain are unclear.

According to The Global Burden of Disease Study 2017, there are more than 41 million adults with gout worldwide.[42-44] Maintaining patient adherence to urate-lowering medications is an ongoing challenge.[42,45] As a result, most gout patients continue to experience gout incidence and flares. In addition, increasing evidence supports that gout is strongly associated with metabolic disorders. The rise in obesity, hypertension, coronary artery disease, and diabetes is an important contributor to the marked increase in gout prevalence globally.[43,46] Gout adds a considerably burden to the public healthcare costs. Currently, non-steroidal anti-inflammatory drugs (NSAIDs) and Colchicine are the most widely used pharmaceuticals for the treatment of gout inflammation and pain. However, accumulating clinical data suggested that benefits of these drugs were limited for gout patients and statistically significantly greater number of adverse events were reported for long term and/or high-dose usages.[47,48] Management and outcomes in people with chronic gout are poor. Novel treatments with high efficacy and minimal side effects are urgently needed for gout patients. Therefore, a renewed interest and a deeper understanding of the pathogenesis of gout will help us develop novel therapies for better gout care.

Systematic transcriptional profiling of the inflamed ankle joint, DRG, and spinal cord in gouty mice

Even though it is well established that an excessively increased serum urate concentration is the major risk factor for the development of gout, the molecular mechanisms driving inflammatory manifestations and severe pain of gout are still poorly studied. Here, we performed subcutaneous injections of MSU crystals in the ankle joints of mice, a model of human gout, and then conducted systematic transcriptome analysis of the ankle joints, DRG and spinal cord, to investigate how intensive pain is triggered and transmitted under acute gouty inflammation. Our results demonstrated that the inflammatory responses and immune pathways were activated and significantly upregulated in all three tissues. For example, the transcriptional levels of the SAA family genes, SAA1/2/3, increased several hundred-fold in the gout model, which may be an indicator of the severity of acute gout inflammation and pain. In addition, such observation might be one of the unique features of gout that are not shared in other pain models or human arthritis. Furthermore, the gene regulatory network provided an overall picture of the interactions of upregulated inflammatory TFs and their target genes during MSU-induced acute inflammation.

Comparative studies of gouty mice and human GWAS of gouty

To assess whether the DEGs identified in our mouse study implicated similar genes as revealed in human GWAS data, we retrieved the genetic association summary statistics of two recent large-scale GWAS cohorts.[27,28] Comparison of our DEGs with the human GWAS revealed 14 overlapping genes. Among these genes, SLC9A2 inhibition reduced serum urate, and activation of KCNQ1 relieved gout pain. Overall, 12 out of the 14 genes were identified among the ankle joint DEGs. Two of the 14 genes, ATXN2 and CUX2 were only identified among the DRG DEGs. In particular, STC1, which is immune modulator and phagocytosis checkpoint that promotes tumor immunotherapy resistance, was the only molecule present in all three types of tissues . Thus, the overlapping genes between our RNAseq DEGs and human GWAS may be excellent candidate genes worthy of further functional characterization.

Comparative gene expression analysis of gouty mice, human OA and human RA

OA and RA are the two most common forms of arthritis. OA is a degenerative disease of the joints due to the wear and tear of the cartilage. RA is an autoimmune disease in which the immune system attacks the joints. Gout is an inflammatory type of arthritis that mainly affects men over the age of 40. Multiple transcriptomic studies have examined the gene expression changes in OA and RA patients. However, comparative studies among OA, RA and gout are rare. We systematically compared our RNAseq data of the ankle joints from gouty mice with RNAseq data from human OA and RA knee joints. For human OA and our gout RNAseq, the upregulated overlapping genes are significantly associated with extracellular matrix organization, ECM-receptor interaction, ossification, and neutrophil activation. In contrast, when comparing our gout RNAseq data with human RA single-cell RNAseq data, enriched biological processes for the overlapping genes include immune response, T-cell activation, leukocyte activation, cell-cell adhesion, and extracellular matrix organization. The results highlight the key differences between human OA and RA, suggesting that gout arthritis may share distinctive features with human OA and RA.

Comparative gene expression analysis of gouty mice, CFA pain model and SNI pain model, and human OA with high knee pain

Sudden and intense pain is one of the characteristics of gout. To determine genes and pathways associated with high pain in gout attacks, we investigated the transcriptomes of DRG and the spinal cord of the gout mouse. For the first time, we revealed that cytokine, cell chemotaxis, leukocyte, and T cell activation related genes were all dramatically upregulated in both the DRG and spinal cord in the gouty mice. We further identified shared DEGs in the DRG and spinal cord among the gouty, CFA, and SNI mouse models. Our results showed that upregulation of cytokine and TNF production in the DRG and increasing cell chemotaxis and leukocyte activation in the spinal cord were shared by the three models. These results suggest that a higher level of inflammatory responses and broad immune activation in the DRG and spinal cord might be one of the underlying causes for severe pain with gout flares.

A previous transcriptomic study of OA patients with low and high pain has identified many DEGs in knee synovial tissues, which are mostly upregulated and likely associated with high pain. We further compared these DEGs with our RNAseq data and then performed functional enrichment analysis. We found that cyclic-nucleotide-mediated signaling, cAMP-mediated signaling, and GPCR signaling pathway coupled to cyclic nucleotide second messengers were significantly enriched for the overlapping upregulated DEGs. Inflammatory mediators abnormally elevate intracellular cAMP levels and activate cyclic nucleotide signaling which in turn causes either sensitization or increased expression of pain channels in sensory neurons.[32,34,50] More importantly, the cAMP-dependent protein kinase (PKA) transduction cascade mediates pain sensation in a dose-dependent manner. Together, our analysis suggests that markedly increased inflammatory responses, immune responses, cyclic nucleotide signaling (including cAMP and cGMP), and GPCR signaling pathways in the DRG and spinal cord may be the primary reasons for significant gout pain.

Alternative splicing analysis of gouty mice

RNAseq data have provided researchers unprecedented insight into the transcriptional landscape of cells. The current sequencing depth is sufficient to reveal the complexity of genomes and transcriptomes, such as gene fusions, mutations/indels, alternative splicing, posttranslational modification, and novel RNA species, in addition to differential gene expression. Alternative splicing is an important mechanism used to increase the diversity of transcripts and proteins in the genome. Alternative splicing is crucial in the regulation of normal cellular function. Aberrant splicing events are commonly reported in human diseases, particularly cancer and neurological diseases. Alternative splicing plays an important role in immune cell development and regulation of inflammatory responses. Targeting alternative splicing has led to the development of novel therapeutics. In the present study, we discovered differential alternative splicing and differential transcript usage in the inflamed ankle joint, DRG and spinal cord of the mouse gout model, which would remain concealed in traditional transcriptomic studies. We found that distinct sets of genes were differentially spliced under gout conditions for different tissues, and these genes were enriched in different biological processes and may be relevant to gout pathogenesis.

Conclusions

In the present study, we systematically investigated the transcriptional profiles of the ankle joint, DRG, and spinal cord in gouty mice. The prominent result was the broad upregulation of immune and inflammatory responses in all three tissues. Integrative analysis with human GWAS data of gout showed that 14 of 37 known susceptibility loci were also DEGs. In particular, STC1 was the only locus upregulated in all three tissues. In addition, differentially alternatively spliced genes and genes with differential transcript usage were detected under gout conditions, which showed tissue specific patterns. Interestingly, for each tissue, DEGs, differentially spliced genes, and genes with differential transcript usage are largely not overlapped. Comparison with human OA and RA transcriptomic datasets indicated that significant upregulation of cAMP/cyclic nucleotide-mediated signaling was likely associated with high pain. Our findings provide key insight into how severe gouty inflammation and pain are triggered by whole-genome gene expression analysis, which not only provides candidate genes to further study gout pathogenesis but also aids in the development of novel therapies for better gout care.

Click here for additional data file.

Supplemental Material for Systematic analysis of inflammation and pain pathways in a mouse model of gout by Yuanlan Fan, Jiajun Yang, Song Xie, Jiangshuang He, Shenghui Huang, Jiawang Chen, Shangying Jiang, Lei Yu, Yujuan Zhou, Xuxia Cao, Xiang Ji and , Yi Zhang in Molecular Pain

Click here for additional data file.

Supplemental Material for Systematic analysis of inflammation and pain pathways in a mouse model of gout by Yuanlan Fan, Jiajun Yang, Song Xie, Jiangshuang He, Shenghui Huang, Jiawang Chen, Shangying Jiang, Lei Yu, Yujuan Zhou, Xuxia Cao, Xiang Ji and , Yi Zhang in Molecular Pain

Click here for additional data file.

Supplemental Material for Systematic analysis of inflammation and pain pathways in a mouse model of gout by Yuanlan Fan, Jiajun Yang, Song Xie, Jiangshuang He, Shenghui Huang, Jiawang Chen, Shangying Jiang, Lei Yu, Yujuan Zhou, Xuxia Cao, Xiang Ji and , Yi Zhang in Molecular Pain

Click here for additional data file.

Supplemental Material for Systematic analysis of inflammation and pain pathways in a mouse model of gout by Yuanlan Fan, Jiajun Yang, Song Xie, Jiangshuang He, Shenghui Huang, Jiawang Chen, Shangying Jiang, Lei Yu, Yujuan Zhou, Xuxia Cao, Xiang Ji and , Yi Zhang in Molecular Pain

Click here for additional data file.

Supplemental Material for Systematic analysis of inflammation and pain pathways in a mouse model of gout by Yuanlan Fan, Jiajun Yang, Song Xie, Jiangshuang He, Shenghui Huang, Jiawang Chen, Shangying Jiang, Lei Yu, Yujuan Zhou, Xuxia Cao, Xiang Ji and , Yi Zhang in Molecular Pain

Click here for additional data file.

Supplemental Material for Systematic analysis of inflammation and pain pathways in a mouse model of gout by Yuanlan Fan, Jiajun Yang, Song Xie, Jiangshuang He, Shenghui Huang, Jiawang Chen, Shangying Jiang, Lei Yu, Yujuan Zhou, Xuxia Cao, Xiang Ji and , Yi Zhang in Molecular Pain

Click here for additional data file.

Supplemental Material for Systematic analysis of inflammation and pain pathways in a mouse model of gout by Yuanlan Fan, Jiajun Yang, Song Xie, Jiangshuang He, Shenghui Huang, Jiawang Chen, Shangying Jiang, Lei Yu, Yujuan Zhou, Xuxia Cao, Xiang Ji and , Yi Zhang in Molecular Pain