Integrated analysis of Transcriptomic Data Reveal AP-1 as a Potential Regulation Hub in Inammation Induced Hyperalgesia Rat

Background Inammation associated chronic pain is a global clinical problem, affecting millions of people worldwide. However, the underlying mechanism of inammation associated chronic pain remains unclear. Complete Freund's adjuvant (CFA) induced cutaneous inammation rat has been widely used as an inammation-induced pain hypersensitivity model. Methods Herein, we presented a transcriptomic prole of CFA-induced rat dorsal root ganglion via an approach targeting gene expression, DNA methylation and post-transcriptional regulation. Results We identied 418 differentially expressed mRNAs, 120 differentially expressed microRNAs and 2670 differentially methylated regions, which are all highly associated with multiple inammation related pathways, including NF-κB signaling pathways and IFN signaling pathways. Integrated analysis further demonstrated that AP-1 network, which may work as a regulator of inammation response, is sophisticatedly regulated at both transcriptomic and epigenetic level. We believe our data will not only provide drug screening targets for treating chronic pain and inammation, but also shed a light on the molecular network inammation induced hyperalgesia. Our study provides the rst comprehensive molecular prole of the CFA-induced inammatory pain model using the integrated omics approach. Our results revealed AP-1 signal may play a central role in regulating pain-related signaling.


Introduction
Chronic pain, which affects more than 17% of the all Canadians, as well as approximately 100 million people in the United States [22; 23], has been one of the most common and distressing issues, in uencing the society and individuals signi cantly [5]. Chronic pain is usually accompanied with in ammation, tissue damage, hyperalgesia as a separate condition [18]. Risks factors of chronic pain include sociodemographic, mental, clinical as well as the biological factors [29]. To date, in ammations and tissue damages are widely believed as the leading biological factors to induce chronic pain [7; 24]. However, more studies on the pathogenesis of chronic pain are imperative at molecular level.
Animal in ammation models are widely used to assess the production of in ammatory mediators.
Complete Freund's adjuvant (CFA)-induced cutaneous in ammation rat model has been widely characterized in literature and has been used for nding novel anti-in ammatory molecules for years [6].
Swelling, hyperalgesia and allodynia can be used to evaluate the in ammation response in the CFAinduced model [6]. Long-term hyperexcitability of the primary sensory neurons in the dorsal root ganglion (DRG), and the secondary sensory neurons in the spinal cord dorsal horn (SDH), are considered as the cause of the alteration in pain perception [30]. CFA-induced chronic in ammation causes chronic pain.
CFA-induced rat has been used as one of the most common chronic pain animal models [14; 15; 36; 37].
However, the mechanism of how CFA-induced chronic pain remains not fully understood.
Recently, integrated analysis of multi-omics data, which relies on next generation sequencing, has revealed numerous complex molecular networks of basic biological process as well as human diseases [26; 32]. DNA, RNA and protein work synergistically to perform certain biological functions. Focused on chronic pain studies for years, we are extremely interested in understanding the underlying biological network of how CFA induce chronic pain in animals. In fact, there exists a pioneer study published by Tan et al. in 2017. Although only an array-based analysis focused on microRNA (miRNA) and mRNA network was used, they were still able to reveal that several miRNA-mRNA interactions contributed to the in ammatory pain process. For instance, they found that miR-124-3p could attenuate in ammatory pain and decrease IL-6R expression in the spinal cord [17]. Aside from this study, miRNAs have already been widely associated with in ammatory models [1; 12; 27; 33]. Inhibiting some of these miRNAs has been considered as novel treatment for chronic pain [13; 20; 21]. However, only focusing on post-transcriptional regulation is not enough to illustrate how CFA-induced in ammatory process is regulated comprehensively. In fact, our previous work has demonstrated that dysregulation of DNA methylation is involved in the CFA-induced in ammation model [14]. A recent study also described that hypomethylation of nerve growth factors affects In ammatory hyperalgesia in rats [39]. These studies suggested that epigenetic regulation, especially DNA methylation could contribute to in ammatory pain. In addition, although Tan et al. provided rat SDH tissue expression pro les days after CFA treatment using microarray technology [17], the transcriptomic pro le of the dorsal root ganglion (DRG) tissues from the chronic pain animal models is still lacking. This suggested that how the gene expression is affected by CFA treatment in the primary sensory neurons is still not clear.
In order to provide a comprehensive pro le of DRG tissues' transcriptomic response to CFA injection, herein, we used a an approach that target DNA methylation, gene transcription and post-transcriptional regulation, which includes mRNA-seq, small-RNA-seq and targeted bisul te sequencing, to detect the transcription level alteration. We identi ed 418 differentially expressed genes, 120 differentially expressed miRNAs and 2670 differentially methylated regions, which are associated with multiple in ammation related pathways. Integrated analysis further indicates that AP-1, a key regulator of in ammation response, is sophisticatedly regulated at both transcriptomic and epigenetic level.

Methods
Animal Model 15 pairs of adult male SD rats in a SPF grade weighing about 200 g ± 5 g were purchased from Experimental Animal Center of Suzhou University. The SD rats maintained at standardized feeding environment including 12-hours shift light-dark cycle and ambient temperature of 24 ± 1°C. The CFAinduced chronic in ammatory pain model was established by unilateral subcutaneous injection of CFA 0.1 mL at sterile foot skin and was described in our previous research [6] [14]. 10 pairs of the rats (CFA and control) were used for following sequencing experiments, while the rest 5 pairs were used for the validation experiments. All experiments were approved by the Animal Research Ethics Committee from The First People's Hospital of Yancheng, Suzhou Municipal Hospital A liated to Nanjing Medical University and Nantong University. Mechanical and Thermal pain threshold measurement was also described in our previous study in detail [14].
RNA-seq, Small-RNA-seq and Targeted Bisul de sequencing For the NGS-based sequencing, DRG tissues from 10 pairs (Randomly selected) of CFA-treated and untreated rats were extracted. For library construction, RNA-seq library was prepared with NEBNext Ultra RNA with Poly-A selection and was sequenced on an Illumina Hi-Seq 4000. Small RNA-seq library was constructed with NEBNext Multiplex Small RNA Library Prep Kit for Illumina (#7560S). The targeted bisul de sequencing library is prepared with a DNA library construction TruSeq DNA LT Sample Prep Kit v2 for Illumina. The C-T transition is performed via C-T transition EpiTect Bisul te kit from Qiagen (Germany).

RNA sequencing analysis
After running fastQC (0.11.8) for quality control, Cutadapt (2.1.0) was used for trimming adapters and low-quality sequence (Phred score less than 20) from the raw fastq les. The cleaned fastq les were then mapped to rn4 rat genome refence by STAR (2.7.0) with suggested setting. Differential expression analysis was then performed by DEseq2 (1.29.0). Signi cant differentially expressed genes was de ned by P value < = 0.05 and log2 fold change > = 0.5. Further pathway enrichment analysis was performed with EnrichR. The gene set enrichment analysis (GSEA) was performed by using javaGSEA2-3.0.
The RNA-seq wig le generated by STAR were then used as input for DaPars (0.9.1) to identify alternative polyadenylation events. Signi cant events were de ned by P value less than 0.05 and absolute delta PDUI more than 0.1 [35].
Small RNA sequencing data analysis Small RNA-seq data was analyzed by using sRNAnalyzer pipeline [34]. Adapter sequences was provided as input to the pipeline for removing adapter contamination. The cleaned reads were then mapped to rat miRNA database to generate count matrix, which was further analyzed by DEseq2 to identify differentially expressed miRNA. Signi cantly expressed miRNA were de ned based on the P value less than 0.05. The miRNA targets were predicted by using 'Multimir' (1.10.0) R package. Finally, the pathway enrichment analysis was performed by using EnrichR.
DNA methylation data analysis TrimGalore (0.6.4) was used for sequencing data quality control, adapter removing and trimming lowquality sequence (Phred score less than 20). The cleaned data was then mapped to the rn4 rat genome reference by using Bismark (0.22.0). The generated CpG methylation ratios were then used as input for metilene to call differentially methylated regions (DMRs) [10]. The signi cant DMRs were de ned by P value less than 0.05 and absolute methylation ratio change more than 0.1. Signi cant DMRs were annotated to nearest genes by using Homer (4.8). Pathway enrichment analysis and motif enrichment analysis were performed by EnrichR and MEME-suite, respectively.

Statistical analysis
Enrichment analysis was detected by Fisher's Exact test and p-values < 0.05 were considered statistically signi cant. False Discovery Rate (FDR) correction was used to generate p-adj value in all expression matrix data (RNA-seq, small-RNA-seq and DNA methylation). All data (mean ± SE) were analyzed with a Student's t test, and a p-value < 0.05 was considered statistically signi cant.
Quantitative Real-time PCR assays DRG tissues from 4 pairs of the rats (Control and CFA) were separated half portion for qPCR assays.
Total RNAs were extracted from cells by using TRIpure reagents (Bioteke Beijing) and then reversetranscribed by using the rst strand cDNA synthesis kit (TSK302S (RT6 cDNA Synthesis Kit Ver 2).
Reverse-transcribed products were used as templates for real-time PCR using the 2×T5 Fast qPCR Mix (SYBR Green I). The primers used for different target genes were as follows: Reg3a- Western blotting DRG tissues from three pairs of the mice were collected for the protein extraction (Total n = 6). DRG tissues were washed three times before the tissue preparation. Total lysates were fractionated by SDS/PAGE and then transferred for 1 hour to the nitrocellulose lter (NC) membranes. The membranes were then blocked for 45 min with Tris-buffered saline containing Tween 20 (TBST), which contain 3% (mass/vol) nonfat dry milk. Antibodies for Egr1 were purchased from Invitrogen a-actin were purchased from Cell Signaling. After the incubation with primary antibodies, the membranes were washed and then incubated with HRP-conjugated anti-mouse or anti-rabbit secondary antibodies. The membrane was then followed by detection using ECL Western blotting substrate (Bio-Rad). Quantitative analysis of Western blot results were calculated using Image J software.

CFA-induced in ammation rat model
To establish the CFA-induced in ammation rat model, we injected male SD rats (200 ± 5 g of body weight) with 0.1 mL CFA on the plantar skin of the right hind paw. Detail methods were described in the Method section and was published in our previous work. 24 hours after the CFA injection, we measured the rat hind paw withdrawal threshold (PWT) as well as the paw withdrawal latency (PWL) in each group. As expected, we discovered a signi cant downregulation of both PWT and PWL in the CFA group compared to the saline treated group (Figure .1). These results validated our previous experiments as well [14].

RNA-seq pro ling of the CFA-induced in ammation model in rat
To investigate the transcriptomic alternation during the CFA-induced in ammation, we performed total RNA-seq with dorsal root ganglion (DRG) from 10 pairs of rats with or without CFA treatment and generated an RNA-seq dataset with high quality in terms of sequencing quality, duplication levels, adapter contamination and mapping rates. Differential expression analysis revealed that 213 and 205 genes were signi cantly up-regulated or down-regulated, respectively ( Figure. 2a). Among all the signi cant differentially expressed genes (DEGs), Reg3b and Reg3a are the top two genes with more than 6-fold increment of expression level ( Figure. 2b). A quantitative real-time PCR (qRT-PCR) assay was also performed to validate the increment of these mRNA expression (Supplementary Fig. 1). It has been wellknown that regenerating gene (Reg) family plays important role in cell proliferation, migration and in ammation response [40]. In order to examine the transcriptomic alternation globally, we also performed pathway enrichment analysis based on all DEGs and observed the enrichment of many in ammation related pathways, such as Interferon alpha/beta signaling pathway, RIG-1-like receptor signaling pathway, TRAF6 mediated NF-κB activation and AP-1 transcription factor network ( Figure. 2c).
Interestingly, The TRAF6 mediated NF-κB has been previously linked to CFA-induced in ammation already [2; 31; 36]. Altogether, these RNA-seq based results not only further validates our success construction of the CFA-induced in ammation rat model but also reveals the speci c transcriptomic alternations induced by CFA treatment in rat (Supplementary Table 1).
Alternative polyadenylation (APA) has been shown to regulate gene expression by regulating the interaction between miRNA and RNA molecules [4]. In order to investigate the gene regulation mechanism of the CFA-induced in ammation, we therefore rst pro led the APA events between CFA-treated and control rats by DaPars. However, no signi cant APA events were detected (Supplementary Table 2), indicating that APA is irrelevant to the expression change observed in the CFA-treated rats. miRNA pro ling of the CFA-induced in ammation model in rat Although we did not observe APA change in the CFA-treated rats, miRNA might regulate the gene expression directly through the miRNA abundance. In fact, there are several miRNAs were reported to be related to CFA-induced in ammation, including miR-134 in DRG tissue, and miR-124, miR-149, miR-3584 in spinal cord dorsal horn [17; 19]. However, these studies did not include a comprehensive database which should rely on a non-biased next generation sequencing platform. In order to provide a landscape of how miRNA network is regulated right after CFA treatment, we next performed small RNA-seq for the same 10 pairs of rat tissues with or without CFA treatment. As expected, the differential expression analysis showed that 16 miRNAs were up-regulated, such as rno-miR-384-3p, which targets immune related genes like Pik3cd, Nfatc3 and Cblb (Supplementary Table 3). We also identi ed 104 signi cantly down-regulated miRNAs, with the top one being rno-miR-181b, which is predicted to target Tnf and Irs2 ( Figure. 3a). Heatmap of miRNA expression pro le also showed the control and CFA treated rats were well clustered ( Figure. 3b). To study the functions of the down-regulated miRNAs unbiasedly, we then performed pathway enrichment based on their target genes. Interestingly, many in ammation related pathways are highly enriched, such as TGF-beta signaling pathway, PDGFB signaling pathway, SHP2 signaling pathway and MAPK signaling pathway (Figure. 3c). These results suggest that, in response to CFA-treatment, miRNA is involved in the regulation of the in ammation response.

DNA methylation pro ling of the CFA-induced in ammation model in rat
Our previous work demonstrated a promoter demethylation of CXCR4 gene is related to CFA-induced in ammation [14]. To further investigate the gene regulation mechanisms of CFA-induced in ammation response, we pro led the methylation landscape of the CFA-treated rats by performing targeted bisul te sequencing. In the comparison of control and CFA-treated rats, we totally identi ed # differential methylated cytosines, suggesting CFA can indeed cause the change of DNA methylation landscape. We further performed differentially methylated regions (DMRs) calling by using the 'de-novo' mode of metilene and identi ed 1401 hypermethylated DMRs and 1269 hypomethylated DMRs ( Figure. Table 4). The majority of DMRs are enriched to the intergenic and intron regions, with only 19 DMRs are located in the promoter regions, suggesting CFA-induced DNA methylation change might not play the role of gene regulation by direct repressing transcription at the gene promoter regions ( Figure. 4b). We then assigned DMRs to the nearest genes and performed pathway enrichment. For both hypermethylated and hypomethylated DMRs, we observed the enrichment of in ammation related pathways, such as Interleukin-4 signaling pathway, Insulin signaling pathway, TGF-beta receptor activation of SMADs, Signaling by NGF and Toll receptor cascades (Figure. 4c). Finally, we did motif analysis with the DMRs and found the enrichment of the motifs for Tbp and Mtf1 transcription factors, both of which has been reported to linked to stress response and in ammation response previously [16; 28]. In summary, these results show that DNA methylation is involved in the gene regulation of CFAinduced in ammation response, probably through in uencing the binding of in ammation related transcription factors.

AP-1 works as a potential regulation hub for CFA-induced in ammation response
After investigating CFA-induced change of gene expression, miRNA and DNA methylation, we examined their interplay in order to illustrate how these three layers of molecular changes are coordinated. Since majority of miRNAs are downregulated, we then focused on their targets genes which had corresponding increased transcription level (571 genes). These target genes are also enriched for many in ammation related pathways, but with the top one being AP-1 transcription factor network ( Figure. 5a) (Supplementary Table 5). Interestingly, the DEGs are also enriched in this pathway, further indicating that miRNA mediated repression of AP-1 network had to be mitigated during the CFA-induced in ammation response ( Figure. 2c and 5b). A validating experiment based on qRT-PCR also showed that Dusp1, one of the two AP-1 related molecules, was signi cantly highly expressed in the CFA group compared to the control (Supplementary Fig. 2). However, we failed to nd the differential expression for Egr1 gene in the qRT-PCR results (Supplementary Fig. 2). In fact, although the RNA-seq result exhibited a 1.5-fold upregulation of Egr1 expression in the CFA group, most of the reads were enriched in the 3' UTR, indicating a potential post transcriptional regulation event. In order to determine whether Egr1 was regulated by CFA treatment, we detected the protein level of Egr1 and found that Egr1 expression was slightly decreased in the CFA group (n = 3, Supplementary Fig. 3). Compared to the upregulation in mRNA level, this opposite trend of the protein expression may be a consequence of the increased expression of 3'UTR. Furthermore, we also studied the interaction between DNA methylation and gene expression. Since the majority of DMRs are annotated to the intergenic regions, we are not able to observe signi cantly anticorrelation between differential methylation and differential expression. However, based on the motif analysis, we found that CFA-induced DMRs were surprisingly enriched for the AP-1 binding motif ( Figure. 5c) (Supplementary Table 6). Given that AP-1 is a well-known key regulator of in ammation response and has been previously linked to CFA-induced in ammation [3], we speculate that AP-1 may serve an important role as a hub subjected to the regulation of both miRNA and DNA methylation upon CFA treatment.

Discussion
The addition of "omics" to certain molecular term re ects a comprehensive assessment of a set of molecules. Integrated analysis of multiple NGS-based datasets, which highly rely on high-throughput sequencing, has already revolutionized the medical research [8]. Chronic pain is a critical health issue globally, affecting millions of people. However, there exists di culty for directly obtaining certain tissues directly from patients clinically. Thus, we believed that understanding the molecular networks of chronic pain animal models through a comprehensive approach is an alternative. Herein, we provided an integrated analysis based molecular pro le, which focuses on the response of rat DRG tissues to CFA injection, at epigenetic, transcription and post-transcriptional regulation level. Overall, we identi ed 418 differentially expressed mRNAs, 120 differentially expressed miRNAs and recognized more than 2,500 DMRs in CFA-treated groups. Through gene set enrichment analysis, we validated some of the previously reported CFA-response related signaling pathways that are also highly related to CFA induced in ammation in our dataset as well, such as, NF-κB signaling pathways [36] and IFN signaling pathways [38]. Moreover, we also identi ed many new genes/pathways that are potentially highly involved in this pain response model, including Reg family genes (Reg3a, Reg3b) and the AP-1 transcription related pathways. It is interesting that after we adjusted the cell heterogeneity of the methylation sequence data by using CHALM [41], a recent invented software for analyzing methylome, we identi ed 6832 signi cant DMRs, suggesting the CFA's impact on rat's methylome is underestimated by the traditional analysis method. Besides in ammation related pathways or terms, the CHALM identi ed DMRs are also enriched to heart contractions and heart rate, indicating that methylome re-wire is involved in the chronic pain linked heart malfunctions or diseases (Supplementary Table 10). Finally, Based on our multi-omics pro ling of the CFA-induced chronic pain model, we selected top 10 differentially expressed genes, miRNAs and differentially methylated regions as the multi-omics signature for our chronic pain model. CFA treatment group and control group can be clearly separated by these signatures (Fig. 6a-c) (Supplementary Table 7 -9), and is consistent to the idea that miRNA and DNA methylation regulate mRNA transcription.
Pancreatitis-associated proteins, which are from Reg families, have been previously linked to modulation of spinal sensory pathways in pathological pain states [9]. Our results demonstrated that Reg3b and Reg3a were signi cantly upregulated for 12-fold and 6-fold, respectively, in CFA-treated groups, indicating that the Reg family genes may play a crucial role in regulating in ammatory pain response in DRG tissues. AP-1 transcription factor is composed of dimerization of a bZIP (basic region leucine zipper) domain via the Fos and Jun subunits. AP-1 regulation network was previously linked to chronic pain response [3]. However, it was never recognized as a central pathway in pain response before. In our dataset, AP-1 network is considered as the regulation central hub. Not only the CFA-induced miRNA and mRNA interactions are highly enriched for AP-1 networks, but also the CFA-induced DMR motifs are enriched for the AP-1 transcription binding activities. Therefore, we hypothesized that inhibiting certain AP-1 network genes, such as Egr1, which is recognized both as a pleiotropic in ammatory trans-activator [25] and a chronic pain contributor [11], could have a chance to alleviate chronic pain. Interestingly, most of the differentially expressed miRNAs (104/120) were downregulated after CFA treatment. The global downregulation of miRNA itself is an interesting phenomenon, which was missed in the previous arraybased study [17]. Notably, in these downregulated miRNAs, many of them are targeting AP-1 network genes. Thus, we speculate adding back of these miRNAs could inhibit the overactivation of AP-1 network, thus alleviate pain. Overall, we believe AP-1 network plays a central role in regulating in ammatory pain responses through a Methylation-transcription-posttranscription regulation axis.
In conclusion, this study provides a comprehensive transcriptomic pro le of the CFA-induced in ammatory pain rat model via an approach that target DNA methylation, gene expression as well as post-transcriptional regulation. Our study has certain limitations. In the rst place, although we included 10 pairs of rats for the CFA treatment, we only have a single time point post CFA injection, which is 24 hours. Future studies should focus on the 48hours, 72hours or even 7 days post treatment, as CFA could induce a chronic in ammation as well. In addition, although we demonstrated that AP-1 is likely to work as a regulation hub for CFA-induced in ammation response potentially, we did not include the biochemistry assays to further investigate the alteration of the AP-1 signal regulation. Future studies could focus on the effect of individual molecule or gene of AP-1 network on regulating chronic pain.

Conclusion
Our study provides the rst comprehensive molecular pro le of the CFA-induced in ammatory pain model using the integrated omics approach. Our results revealed AP-1 signal may play a central role in regulating pain-related signaling.