Transcriptome Analysis of Dorsal Root Ganglion in Rats with Knee Joint Inflammation

Introduction

Patients with rheumatoid arthritis (RA) suffer from chronic widespread pain and chronic regional pain. These symptoms significantly decrease their quality of life.¹ However, current therapies for chronic knee pain are either not effective or with serious side effects.² Thus, it is urgent to understand the underlying mechanisms for developing a new therapy for such pain.

Long non-coding RNAs (lncRNAs) are transcribed RNA molecules with no or very low protein-coding potential.³ Numerous lncRNAs have been identified to exert vital parts in regulating the expression of specific genes and are involved in different biological functions.^4,5 The abnormal expression of lncRNAs will result in different neurological disorders. Previous studies have shown that changes in the expression of lncRNAs in dorsal root ganglia (DRG) neurons have been implicated in the pathogenesis of neuropathic pain.^6–11 Gene expression alteration in knee cartilage, synovium or peripheral blood from RA patients have been identified using microarray,^12–15 but it remains unclear whether lncRNAs-mediated gene regulation occurs in primary sensory neurons of DRG during RA-like joint inflammation.

In this study, we analyzed lncRNAs and related mRNAs in the DRG after knee joint inflammation induction. A unilateral intra-articular injection of complete Freund’s adjuvant (CFA) in rats was performed to prepare a rodent model of joint inflammation. Using microarray RNA sequencing, we identified 69 lncRNAs and 102 related mRNAs that displayed significant expression alteration in the DRG after knee joint inflammation. These lncRNAs and mRNAs could be potential targets for developing novel therapies for inflammatory joint pain. As indicated by the technical validation in this study, high-quality data were obtained, and the expression levels of genes took place within the desired regions of genes. Our findings provide precious data for the scientific community to more intensively investigate the alterations of DRG at molecular level in the CFA-induced inflammatory joint pain model.

Materials and Methods

Animals

Male inbred Wistar rats (200 g) were provided by Animal Laboratories Center of Zhengzhou University in China and were housed in the central facility at Zhengzhou University School of Medicine under a standard 12 h light/12 h dark cycle, with water and food available adlibitum. The study protocol gained approval from the Institutional Ethics Review Committee at Zhengzhou University. Each experiment was performed according to the guidelines for the care and use of laboratory animals released by the National Institutes of Health (NIH). The animal number utilized together with their sufferings was minimized as far as possible. The reviewers were blinded to animal treatment.

CFA-Induced Inflammatory Joint Pain Rat Model

Three percent isoflurane was injected into rats for anesthesia, then 150 μL CFA (1 mg/mL Mycobacterium tuberculosis, Sigma, USA) was also injected into their right knees via the patella tendon. We prepared the animal model as described previously.¹⁶ As a control, 150 μL of saline (0.9% sodium chloride) was injected into the right knee joint. Rats were allowed to recover from anesthesia before being returned to their cages.

Determination of Rat Knee Diameter

For animals under anesthesia, their knee joint diameters just below the level of the patella were determined using a digital vernier caliper (ID:110–601, Langir Electric CO., Ltd.) as described previously.¹⁶ On the 7th day following intra-articular injection with CFA, the edema in knee was quantified.

Pain Behavioral Test

Mechanical behavioral testing (6 rats per group from a different set of rats) was performed according to our prior articles.^17–19 The up–down testing paradigm was employed to measure the 50% paw withdrawal threshold (PWT) responding to the mechanical stimulus by von Frey filaments (North Coast Medical Inc., Gilroy, CA, USA), according to the previous description.^18,19

RNA Extraction

Right L4 and L5 DRGs from CFA group and saline control group were harvested on day7 after CFA or saline injection (3 rats per group from a different set of rats). Next, the mini RNA easy kit (QIAGEN, Valencia, CA, USA) was used to extract total RNA, and then genomic DNA was digested in the column in accordance with specific protocols. Thereafter, NanoDrop ND-1000 (Thermo Scientific, Wilmington, DE, USA) was used to measure RNA quality and quantity, whereas the standardized denatured AGE was applied in assessing the RNA integrity.

Microarray Analysis

A lncRNA array of rats was prepared to profile the lncRNAs as well as the protein-coding genes. Meanwhile, gene-coding probes were printed for once, whereas the lncRNAs-coding ones were printed for thrice, and negative ones together with housekeeping genes-coding ones were printed several times to ensure the quality of hybridization. The slightly modified Agilent One-Color Microarray-Based Gene Expression Analysis protocol (Agilent Technology) was used to label samples and hybridize the array. In brief, rRNA was removed from total RNA using the mRNA-ONLY™ Eukaryotic mRNA Isolation Kit (Epicentre), followed by mRNA purification. Subsequently, all samples were subjected to amplification and transcription to the fluorescence cRNA down the full-length transcripts in the absence of 3ʹ bias according to the random priming approach (Arraystar Flash RNA Labeling Kit, Arraystar). Then, the RNeasy Mini Kit (Qiagen) was used to purify those labeled cRNAs. The NanoDrop ND-1000 was used to measure labeled cRNAs content along with the activity. Afterward, 10 × blocking reagent (5 μL) together with 25 × fragmentation buffer (1 μL) were added to fragment the labeled cRNAs (1 μg each), followed by 30 min heating under 60 °C. At last, those labeled cRNAs were diluted by 2 × GE hybridization buffer (25 μL). Later, the hybridization solution (50 μL) was added to the gasket slide for the assembly of lncRNA expression microarray slides. Then, each slide was incubated within the Agilent Hybridization Oven under 65 °C for 17 h. Afterward, an Agilent DNA Microarray Scanner (Agilent, part # G2505C) was utilized to wash, fix, and scan those hybridized arrays.

Data Processing

Altogether six samples were used for sequencing, multiplexing, analysis for differentially expressed (DE) genes, as well as analysis of expression of transcripts and lncRNAs. In brief, the array images were analyzed by Agilent Feature Extraction (version 11.0.1.1). The GeneSpring GX v12.1 (Agilent Technologies) was employed to normalize quantiles and process data. When the raw data were normalized for their quantiles, mRNAs and lncRNAs with flags in ≥3 among the six samples in Marginal or Present (“All Targets Value”) were selected to analyze data. Then, the significant DE lncRNAs and DE mRNAs between CFA injection and normal control groups were discovered according to the same threshold of P-value as that in screening DE lncRNAs and DE mRNAs. The data were submitted to NCBI Gene Expression Omnibus (GEO: GSE146881). This GEO project includes raw data for all samples.

Enrichment and Pathway Analysis

GO analysis facilitates to interpret the properties of genes and their products within an organism (http://www.geneontology.org). We used top GO for GO analysis of differential mRNAs to infer their biological functions: 1) Biological Processes (BPs); 2) Molecular Functions (MFs); 3) Cellular Components (CCs). The KEGG pathway analysis was performed for identifying those biological pathways related to DE mRNAs.²⁰

Quantitative Real-Time PCR

The ThermoScript reverse transcriptase (Invitrogen) was utilized to reversely transcribe total RNA in accordance with specific protocols. Table 1 lists the primers used to amplify cDNA through qRT-PCR. All samples were examined for three times by the BIO-RAD CFX96 RT-PCR system. The PCR conditions were shown below: 3 min of initial denaturation under 95°C; 10 s under 95°C, 30 s under 60°C, 30 s under 72°C for 40 cycles. The ipsilateral to contralateral lncRNA expression ratio was determined by the standardized 2^−∆∆Ct method. GAPDH was used as the internal reference.

Table 1 Primers for RT-qPCR

Establishment of the Co-Expression Network

Pearson correlation calculation (PCC) was utilized to construct the coding-noncoding gene co-expression (CNC) network to normalize the DE lncRNAs signal intensity relative to that of DE mRNAs. The results of PCC > 0.9 were screened as interrelationship pairs. We used the Cytoscape (v2.8.1) to draw a CNC network.

Statistical Analysis

All values were presented in the manner of means ± SEM. SPSS 17.0 was employed for all statistical analyses. Differences in specific lncRNAs or mRNAs levels between CFA injection and saline control groups discovered from microarray as well as PCR were analyzed by Student’s t-test. A difference of two-sided p < 0.05 indicated statistical significance.

Results

CFA-Induced Inflammatory Joint Pain

In the CFA-induced inflammatory joint pain rat model, we observed that unilateral injection of CFA into the knee joint not only significantly increased knee joint diameter (Figure 1A) but also markedly decreased PWT (Figure 1B) at the ipsilateral side on day 7 post-injection, which indicates intra-articular injection of CFA induces inflammatory joint pain. Our observation is consistent with previous studies.^16,21

Figure 1 Injection of CFA into the rat knee joint causes knee joint swelling and induces inflammatory joint pain. (A) Intra-articular injection of CFA significantly increased knee joint diameter on day 7 post-injection. (B) Intra-articular injection of CFA significantly decreased 50% PWT at the ipsilateral, but not contralateral, side on day 7 post-injection. *p < 0.05, **p < 0.01 vs the saline control group (n = 6 for each group).

Expression Profiles of lncRNAs and mRNAs in the DRGs After CFA Injection

On the 7th day following CFA treatment, a graphical overview of the expression signatures of lncRNAs and mRNAs in the L4 and L5 DRGs were obtained by using scatter plot and hierarchical cluster analyses of the microarray results. As observed from the scatter plot, numerous mRNAs and lncRNAs showed differential expression in CFA treatment group compared with saline control group (Figure 2A and B). Besides, it was illustrated from the hierarchical clustering analysis on mRNAs and lncRNAs that, those respective three samples under saline or CFA treatment were clustered, with identical signal intensity in CFA or saline control group (Figure 2C and D). The heatmap showed that DE lncRNAs or mRNAs that were either up-regulated or down-regulated twofold (Figure 2E and F), indicating the high level of concordance in both saline-treated and CFA-treated samples.

Figure 2 Intra-articular injection of CFA alters expression profiles of lncRNAs and mRNAs in the L4 and L5 DRGs on day 7 post-injection. (A and B) The scatter plots showed that a large number of lncRNAs (A) and mRNAs (B) were differentially expressed between CFA and saline groups. (C and D) The hierarchical cluster analysis of all lncRNAs (C) or mRNAs (D) showed that the 3 saline-treated and 3 CFA-treated samples were clustered together, respectively, and signal intensity was consistent in the saline or CFA group. (E and F) The heatmap showed that DE lncRNAs (E) or mRNAs (F) with either up-regulated or down-regulated twofold were magnified, indicating the high level of concordance in both saline-treated and CFA-treated samples.

DE lncRNAs and mRNAs Following RA Induction

We further analyzed DE lncRNAs and mRNAs using the significance analysis of microarray. Totally, 69 DE lncRNAs, including 42 up-regulated and 27 down-regulated, were significantly altered in the CFA-treated group compared with the saline control group. lncRNAs with the highest up-regulation level included MRAK139263, uc.247+, MRAK163594, AY539883, MRAK081329, and MRAK081329, of which MRAK139263 showed the largest up-regulation (log2 fold change = 4.3602434). lncRNAs with the highest down-regulation level included XR_005599, uc.288+, AY321352, uc.83+ and uc.344-, of which XR_005599 showed the largest down-regulation (log2 fold change = −3.6773075). The detailed information about the top 20 up-regulated and 20 down-regulated lncRNAs is listed in Table 2. For the DE mRNAs, 51 up-regulated mRNAs and 51 down-regulated mRNAs were identified in the CFA-induced RA group compared with the saline control group. Those DE mRNAs contained multiple known genes related to the processing of pain, such as 1) Scn2a1 (encodes NaV1.2, a voltage-dependent sodium channel), 2) C1qtnf3 (encodes C1q/TNF-related protein 3, and 3) Mpo (encodes myeloperoxidase). In addition, several genes with unknown roles in pain also showed significant alteration following RA induction, including Ampd1 (encodes adenosine monophosphate deaminase 1) and S100a8 (encodes S100 calcium-binding protein A8). The detailed information about the top 20 up-regulated and 20 down-regulated mRNAs are listed in Table 3.

Table 2 The Detail Information of the Top 20 Up-Regulated and 20 Down-Regulated lncRNAs

Table 3 The Detail Information of the Top 20 uUp-Regulated and 20 Down-Regulated mRNAs

Class Distribution Analysis of DE lncRNAs

LncRNAs are shown previously to modulate overlapping or adjacent mRNAs expression within the genome.^22–24 In this study, all DE lncRNAs were classified into five categories according to their genomic proximity to protein-coding genes: 1) intergenic lncRNAs without a coding transcript in the 30 kb range of lncRNAs; 2) sense lncRNAs whose exons were overlapping with the exons of coding transcripts or introns in the coding transcripts in an identical genomic strand; 3) antisense lncRNAs transcribed based on an antisense strand, which overlapped with the coding transcripts or did not share the overlapping exons; 4) bidirectional lncRNAs directed head-to-head to the coding transcripts in the range of 1000 bp; 5) other lncRNAs. In the 69 DE lncRNAs modulated by CFA-induced joint inflammation, we identified that 20 Intergenic lncRNAs, 17 Sense lncRNAs, 10 Antisense lncRNAs, 1 Bidirectional lncRNAs, and 22 other types of lncRNAs (Figure 3).

Figure 3 Distribution of various types of DE lncRNAs. Five classes (Intergenic lncRNAs, Sense lncRNAs, Antisense lncRNAs, Bidirectional lncRNAs, and other lncRNAs) were identified. In the 69 DE lncRNAs modulated by CFA-induced joint inflammation, we identified 20 Intergenic lncRNAs, 17 Sense lncRNAs, 10 Antisense lncRNAs, 1 Bidirectional lncRNAs, and 22 other type of lncRNAs.

Quantitative Real-Time PCR Validation of lncRNA Expression

To validate the results of microarray, the up-regulated lncRNAs (MRAK163594 and uc.247+) and the down-regulated lncRNAs (uc.83-, XR_005599, NONRATT021972, and AY321352) were analyzed by quantitative real-time PCR. The DRG tissues were harvested for assessing lncRNA expression on day 7 after intra-articular injection of CFA or saline. We observed that MRAK163594 and uc.247+ were significantly increased in CFA-treated rats in comparison with normal controls (Figure 4A and B) and that uc.83-, XR_005599, NONRATT021972, and AY321352 were remarkably declined in CFA-treated rats in comparison with normal controls (Figure 4C−F). Those markedly changed lncRNAs were used to construct the CNC network (Figure 4G).

Figure 4 Quantitative real-time PCR validation of six regulated lncRNAs in the DRG and CNC network. (A and B) The expressions of MRAK163594 (A) and uc.247+ (B) were significantly increased on day 7 after intra-articular injection of CFA. (C−F) Four lncRNAs including uc.83- (C), XR_05599(D), NONRATT021972 (E), and AY321352 (F) were significantly decreased on day 7 after intra-articular injection of CFA. *p < 0.01, **p < 0.01, ***p < 0.001vs the saline control group (n =3 for each group, the error bars show SEM). (G) In the CNC network drawn with the Cytoscape (v2.8.1), pink circle nodes represent non-coding gene, green hexagonal nodes represent the coding gene, yellow star-marked verified lncRNA and as the predicted potential target, solid lines represent a positive correlation between the two nodes, and dotted lines represent negative correlation between the two nodes.

Functional Prediction of DE mRNAs in Inflammatory Joint Pain

To explore the molecular mechanism of inflammatory joint pain, we conducted GO and KEGG enrichment analysis of regulated genes in the CFA-induced inflammatory joint pain rat model. Using GO, we showed that 1) the most significantly enriched molecular functions of the up-regulated genes included protein complex binding, exonuclease activity, and antioxidant activity (Figure 5A), 2) the most significantly enriched biological processes of the up-regulated genes included defense response, response to external stimulus, and response to stress (Figure 5B), 3) the most markedly enriched cellular components of the up-regulated genes included extracellular space, extracellular region part, and extracellular region (Figure 5C), 4) the most significantly enriched molecular functions of the down-regulated genes included alcohol binding, transporter activity, and calcium channel binding (Figure 5D), 5) the most significantly enriched biological processes of the down-regulated genes included muscle organ development, muscle structure development, and skeletal muscle tissue development (Figure 5E), and 6) the most markedly enriched cellular components of the down-regulated genes included extracellular region, extracellular space, extracellular region part, and cell surface (Figure 5F).

Figure 5 Biological functions analysis of DE mRNAs using GO. (A−C) The most significantly enriched molecular function (A), biological process (B), and cellular component (C) of up-regulated mRNAs after CFA-induced inflammatory joint pain. (D−F) The most significantly enriched molecular function (D), biological process (E), and cellular component (F) of down-regulated mRNAs after CFA-induced inflammatory joint pain.

Using KEGG, we showed that 1) the up-regulated genes were involved in inflammatory mediators-produced regulation of transient receptor potential channels, nicotine addiction, intestinal immune network for IgA production, and vitamin B6 metabolism (Figure 6A), and 2) the down-regulated genes were involved in hypertrophic cardiomyopathy, peroxisome proliferator-activated receptor (PPARs) signaling pathway, nitrogen metabolism, and drug metabolism-cytochrome P450 (Figure 6B).

Figure 6 Pathway analysis for 51 up-regulated and 51 down-regulated mRNAs using KEGG. (A) The involved significant pathways for up-regulated genes in the CFA group. (B) The involved significant pathways for down-regulated genes in the CFA group.

Discussion

Inflammatory joint pain is associated with gene expression changes within the DRG. Identifying altered genes is critical for defining the mechanisms underlying inflammatory joint pain pathogenesis. In the present study, we conducted the Rat lncRNA and mRNA array to reveal transcriptome alterations after induction of inflammatory joint pain with CFA.

We used quantitative real-time PCR to verify the expression of six randomly selected DE lncRNAs in the DRG on day 7 after CFA injection, including up-regulated lncRNAs MRAK163594 and uc.247+ and down-regulated lncRNAs uc.83-, XR_005599, NONRATT021972, and AY321352. We constructed a CNC network for the most significantly changed lncRNAs. In the CNC network, the up-regulated MRAK163594 and the down-regulated uc.83- were co-expressed with known gene Kcnj16,²⁵ and the up-regulated MRAK163594 and uc.247+ were co-expressed with pain-related gene Scn2a,²⁶ which contributes to many pain conditions, such as inflammatory joint pain.^21,27 In addition, a previous study has shown that lncRNA uc.217 regulates neurite outgrowth in DRG neurons following peripheral nerve injury.²⁸ By co-expressing with Scn9a, the lncRNA uc.217 is likely to play an important role in CFA-induced inflammatory joint pain.²¹

Furthermore, the identified DE mRNAs in our study contain several known genes involved in pain processing. 1) Scn2a (encodes NaV1.2, a voltage-dependent sodium channel): this ion channel has been demonstrated to regulate pain conduction;²⁶ 2) C1qtnf3(encodes C1q/TNF-related protein 3): this protein plays an important role in the development of RA);²⁹ 3) Mpo (encodes myeloperoxidase): the enzyme activity is an index of neutrophilic granulocyte infiltration, which has been recognized as a marker for pain.^30,31

On the other hand, the identified DE mRNAs also contain several genes with unknown roles in pain. 1) Ampd1 (encodes adenosine monophosphate deaminase 1): Patients with Ampd1 deficiency experience muscle pain (myalgia) or weakness after exercise;³² 2) S100a8 (encodes S100 calcium-binding protein A8): the S100a8 protein is positively associated with increased knee symptoms, cartilage defects in patients with knee osteoarthritis.³³ In addition, KEGG pathway analysis for the DE mRNAs indicates some potential targets, such as inflammatory mediator-regulated transient receptor potential channels, which have been demonstrated as important molecules involved in joint inflammation in the DRG.³⁴

Conclusion

Our results suggest that aberrantly expressed lncRNAs contribute to the development of inflammatory joint pain.