Twenty-seven adult, lab male Sprague–Dawley (SD) rats (180–220 g) were randomly divided into three groups (Sham, CFA, and CFA+BoNT/A group) of nine rats each. All rats were housed in a standard laboratory environment (22 ± 2°C room temperature, 12-h light/dark cycle with lights) with six animals per cage, fed according to a standard rat diet and free water. The housing conditions and experimental procedures were in accordance with the protocols approved by the Bioethics Committee of Qingdao University. The procedures were approved by the Institutional Animal Care and Use Committee and performed in accordance with the NIH guidelines for the care and use of laboratory animals and the International Association for the Study of Pain.

Monoarthritis (ankle arthritis) was induced as previously reported (Butler et al., 1992). After being anesthetized with 5% isoflurane in oxygen, the rats were placed in a prone position and injected with 50 μl of CFA (Sigma, United States) into the left ankle joint cavity as the CFA group. The control group was injected with saline as Sham group. After 21 days, the CFA model was set up. After the rats were anesthetized with isoflurane, 20 μl of BoNT/A (Allergan Pharmaceuticals Ireland) diluted to 3 U with saline (3 U BoNT/A/20 μl) or 20 μl of saline was also injected into the left ankle joint cavity in the Sham, CFA, and CFA+BoNT/A groups.

On the fifth day after BoNT/A injection, the three groups of rats (n = 3/group) were deeply anesthetized. Then, they were perfused with 0.9% saline and 4% paraformaldehyde (pH 7.4) in sequence. The left DRG was removed and placed in perfusion fixative (4°C) for 24 h. After fixation, paraffin embedding was performed, and 8-μm-thick sections were cut along the long axis of the DRG and mounted on 3-aminopropyl-triethoxysilane-coated glass slide. Sections were sealed with phosphate-buffered saline (PBS) containing 3% donkey serum albumin and 0.3% Triton X-100, and then incubated overnight with primary antibody, which was mouse anti-rat cl-SNAP-25 (1:500, MyBioSource, United States). After washing with PBS, DRG sections were incubated with donkey anti-mouse IgG (1:1,000 dilution Invitrogen) labeled with fluorescein isothiocyanate at 37°C for 40 min to determine the immunoreactivity of Cl-SNAP-25. Images were captured on a Leica DMIRB fluorescence microscope at a magnification of 20x. The immunofluorescent staining of protein expression was quantified and analyzed by ImageJ.

Five days after BoNT/A treatment, the left L4 DRG of the rats were extracted from the Sham group, the CFA group, and the CFA+BoNT/A group (n = 3/group). Next, the total RNA of each sample was extracted with TriPure isolation reagent according to the instructions of the kit (Roche), and then genomic DNA was digested in the column in accordance with specific protocols. Thereafter, NanoDrop ND-1000 (Thermo Scientific, Wilmington, DE, United States) was used to measure RNA quality and quantity, and RNA integrity was assessed using the RNA Nano 6000 Assay kit and the Bioanalyzer 2100 System (Agilent, Santa Clara, CA, United States). Total RNA was used as input material for the RNA sample preparations. Briefly, mRNA was purified from total RNA by using poly-T oligo-attached magnetic beads. Fragmentation was carried out using divalent cations under elevated temperature in First Strand Synthesis Reaction Buffer (5×). First-strand cDNA was synthesized using a random hexamer primer and M-MuLV Reverse Transcriptase, and then RNaseH was used to degrade the RNA. Second-strand cDNA synthesis was subsequently performed using DNA Polymerase I and dNTP. Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. After adenylation of 3’ ends of DNA fragments, Adaptor with a hairpin loop structure was ligated to prepare for hybridization. In order to select cDNA fragments of preferentially 200 bp in length, the library fragments were purified with AMPure XP system (Beckman Coulter, Beverly, United States). After PCR amplification, the PCR product was purified by AMPure XP beads, and the library was finally obtained. In order to ensure the quality of the library, the library needs to be tested. After the construction of the library, the library was initially quantified by Qubit2.0 Fluorometer and then diluted to 1.5 ng/μl, and the insert size of the library is detected by Agilent 2100 Bioanalyzer. After insert size meets the expectation, qRT-PCR is used to accurately quantify the effective concentration of the library (the effective concentration of the library is higher than 2 nM) to ensure the quality of the library. After the library is qualified, the different libraries are pooled according to the effective concentration and the target amount of data off the machine and then sequenced by the Illumina NovaSeq 6000. The end reading of 100 bp pairing is generated. The basic principle of sequencing is to synthesize and sequence at the same time (Sequencing by Synthesis). Four fluorescently labeled dNTP, DNA polymerase, and splice primers were added to the sequenced flowcell and amplified. When the sequence cluster extends the complementary chain, each dNTP labeled by fluorescence can release the corresponding fluorescence. The sequencer captures the fluorescence signal and converts the optical signal into the sequencing peak by computer software, so as to obtain the sequence information of the fragment to be tested.

The image data measured by the high-throughput sequencer are converted into sequence data (reads) by CASAVA base recognition. Raw data (raw reads) of fastq format were firstly processed through in-house perl scripts. In this step, clean data (clean reads) were obtained by removing reads containing adapter, reads containing N base, and low-quality reads from raw data. At the same time, Q20, Q30, and GC content of the clean data were calculated. All the downstream analyses were based on the clean data with high quality. HISAT2 (v2.0.4; parameters: -p 8 –phred64 –sensitive -I 1 -X 1000) was employed to map the clean reads to the rat reference sequences. Then, clean reads were mapped to the reference transcript using Bowtie2 (v2.2.5; parameters: -q –phred64 –sensitive –dpad 0 –gbar 99999999 –mp 1 ,1 –np 1 –score-min L,0, -0.1 - p 16 -k 200), and the gene expression values were calculated by RSEM (v1.3.1; parameters: default). Differential expression analysis of three groups was performed using the DESeq2 R package (1.20.0). DESeq2 provides statistical routines for determining differential expression in digital gene expression data using a model based on the negative binomial distribution. The resulting p-values were adjusted using Benjamini and Hochberg’s approach for controlling the false discovery rate. |log2(foldchange)| >2 and adjusted p-value <0.05 were set as the threshold for significantly differential expression.

Functional enrichment analysis of DEGs was performed using KEGG and Metascape. KEGG analysis was obtained from KEGG, and biological processes were performed by Metascape.

Protein–protein interaction (PPI) network analysis was performed using the online database STRING to predict PPIs. In addition, the PPI network was constructed and visualized using Cytoscape V3.8.2 software. MCODE was used to conduct cluster analysis on the gene network to determine key PPI network modules.

Five days after BoNT/A treatment, total RNA was extracted from the DRG of the rats in each group (n = 3/group) using the TriPure Isolation Reagent according to the instructions of the kit (Roche). RNA purity was determined by Quantus Fluorometer. First, the total template RNA was reverse transcribed into complementary DNA (cDNA). Then, the SYBR PrimeScript PLUS RTPCR kit was mixed with the cDNA samples for qPCR cycling in the AriaMX HRM analysis. GAPDH was used as the reference gene to normalize cycle threshold (Ct) data, and the comparative Ct method was used to calculate each sample.

Throughout the research process, a statistically significant difference analysis (p < 0.05) was performed using GraphPad Prism 8 (GraphPad Software, CA, United States). All of the data are shown as mean ± standard deviation (SD). Two-way repeated measure analysis of variance (RMANOVA) was used to compare the results of pain behavior test (PWT and PWL) in each group. Immunofluorescence staining and RT-PCR analysis were performed to analyze the differences between the groups by one-way ANOVA.

All animals developed prolonged mechanical allodynia after CFA injection. During the entire observation period, PWT of mechanical stimulation was significantly lower than that of the Sham group (Figure 1A). Similarly, PWL of thermal stimulation was distinctly shorter than that of the Sham group (Figure 1B). The PWT of CFA+BoNT/A group on mechanical stimulation increased significantly at 5 days after BoNT/A injection. Similarly, the PWL to thermal stimulation was also significantly reversed by the injection of BoNT/A.

The expression of cl-SNAP-25 in DRG was analyzed to determine whether retrograde transport of BoNT/A could partially explain the observed effect of anti-inflammatory pain. As reported by Restani et al., it is reasonable to use the cl-SNAP-25 as a marker for BoNT/A. Therefore, the anti-cleaved SNAP-25 antibody was used to mark the BoNT/A action at distant sites. The results showed typical examples of immunofluorescence image taken from ipsilateral dorsal root ganglion (Figure 2). The expression of cl-SNAP-25 could hardly be discovered in the Sham group injected with saline and the CFA group. Conversely, the expression of cl-SNAP-25 was obviously detected after injection of BoNT/A.

To determine the potential mechanism of BoNT/A’s anti-inflammatory effect, mRNA was respectively extracted from the Sham, CFA, and CFA+BoNT/A treatment groups, and then transcriptome analysis was performed using RNA-seq. The analysis results showed that 97 DEGs were found between the CFA group and the Sham group. Among them, CFA induced the upregulation of 29 DEGs and downregulation of 68 DEGs. Similarly, a total of 71 DEGs were found between the CFA+ BoNT/A group and the CFA group, of which 52 genes were upregulated and 19 genes were downregulated.

To further analyze the genes associated with CFA-induced arthritis pain, 16 different GO terms enrichment analyses were performed on 97 DEGs between the CFA group and the Sham group. All terms were listed, several of which were associated with response to lipopolysaccharide and inflammatory response (Figures 3A,B). In addition, the study also performed pathway analysis on the screened DEGs and matched them to the physiological processes involved in CFA-induced arthritis pain in rats using the KEGG pathway analysis. CFA affected some signaling pathways, for example, metabolic pathways and IL-17 signaling pathway (Figure 3C). We obtained a PPI network that contains 78 nodes and 262 edges. Among 97 DEGs between the CFA group and the Sham group, 21 did not form molecular networks with other molecules. The nodes were used to represent genes and edges to indicate interactions between genes. Upregulated genes were marked in red, and downregulated genes were labeled in cyan (Figure 3D). A cluster including 11 genes was obtained by MCODE in Cytoscape. We discovered these genes (SERPINE1, MMP9, MMP8, PF4, S100A8, S100A9, MPO, CAMP, HP, FN1, and FCNB) encoding proteins about inflammatory response, and these genes might be the key genes of CFA-induced arthritis pain in rats (Figure 3E).

To identify the key genes involved in the analgesic effect of BoNT/A on arthritis, this study analyzed the CFA group vs. the Sham group and the CFA+ BoNT/A group vs. the CFA group. Twenty-two overlapping genes were found (Figure 4A). Furthermore, among the 21 DEGs, 19 genes were upregulated and two genes were downregulated, which were compared between the Sham group and the CFA group. However, the expression levels of these 21 different genes were reversed between the CFA+ BoNT/A group and the CFA group (Figure 4B). The study also performed GO and KEGG analysis on DEGs between the CFA+ BoNT/A group and the CFA group. The results indicated that the mechanism of the anti-inflammatory effect of BoNT/A against CFA-induced arthritis pain involves the reactions to lipopolysaccharides and peptides (Figures 5A,B) and their related signaling pathways such as PI3K-Akt signaling pathway and IL-17 signaling pathway (Figure 5C). A PPI network was obtained that contains 53 nodes and 163 edges. Among the 71 genes, 18 were not involved in the formation of the molecular network. The nodes were used to represent genes and edges to indicate interactions between genes. Upregulated genes are marked in red, and downregulated genes are labeled in cyan (Figure 5D). We gained a cluster including 11 genes by MCODE in Cytoscape. We discovered these genes (MMP9, MPO, S100A9, MMP8, FN1, HP, S100A8, ELANE, CAMP, NGP, and FCNB) encoding proteins mostly about chronic inflammatory response, and these genes might be involved in the anti-inflammatory effect of BoNT/A against arthritis pain (Figure 5E).

Three genes with high reliability and associated with inflammatory response were selected for verification. MCODE identified S100A9, S100A8, and MMP8 as the genes with higher score in cluster. We detected the expression of these genes in DRG by qRT-PCR. The results showed that the expression levels of S100A9, S100A8, and MMP8 in rats with arthritis pain were visibly lower than the Sham group. After injection of BoNT/A, S100A9, S100A8, and MMP8 showed high expression levels compared with the CFA rat model (Figure 6).