Introduction

Trigeminal neuropathic pain is a severe health problem and is difficult to treat because of its complex mechanisms. It often occurs after injuries to the trigeminal sensory fibers induced by trauma and tooth iatrogenic injuries from dental treatments, such as local anesthetic injections, root canal therapies, extractions, oral surgery, dental implants, orthognathic surgery, and other invasive procedures.1,2 Although lots of potential intervention molecular targets in the trigeminal ganglion (TG) for treating trigeminal neuropathic pain have been found,3 for instance, multiple groups of molecules (ATP, glutamate, SP, IL-1β, and TNF-α), receptors (P2X, NMDAR, NK1, and GPCRs), ion channels (voltage-gated potassium channels and voltage-gated calcium channels), and intracellular signaling pathways (MAPKs, PKC, PKA),3–5 consistently effective treatment options are still not available.3,6 Therefore, more preclinical research is expected to elucidate the molecular mechanism underlying trigeminal neuropathic pain and identify novel therapeutic targets to develop a more effective treatment for trigeminal neuropathic pain.

In recent years, emerging evidence suggests that non-coding ribose nucleic acids (ncRNAs) play a vital role in the pathophysiological process of neuropathic pain,7,8 which brings novel landscapes for the diagnosis and treatment of trigeminal neuropathic pain. According to their size, ncRNAs mainly consist of micro RNAs (miRNAs), long ncRNAs (lncRNAs), circular RNAs (circRNAs), and small nucleolar RNAs (snoRNAs), which typically do not encode proteins but functionally regulate protein expression.9,10 miRNAs in the TG have been proved to play an essential role in oral cancer pain and orofacial inflammatory pain.11,12 In addition, emerging data indicate that lncRNAs and circRNAs interact with miRNAs to regulate mRNA expression in the progression of neuropathic pain.8,13 Notably, recent studies have reported that lncRNAs in the TG also contribute to trigeminal neuropathic pain and chronic inflammatory temporomandibular joint pain.14,15 These data support that ncRNAs in the TG may be associated with trigeminal neuropathic pain. However, the regulatory functions of ncRNAs in the TG under trigeminal neuropathic pain and their functional mechanisms have not been systematically described.

Chronic constriction injury of the infraorbital nerve (CCI-ION) is a classic animal model of trigeminal neuropathic pain,16,17 resulting in long last orofacial pain status. In this study, a whole-transcriptome sequencing technique was performed to systematically analyze ncRNAs and mRNA expression signatures in the TG of mice with or without trigeminal neuropathic pain induced by CCI-ION. Furthermore, bioinformatics analysis interpreted the potential biological functions of differentially expressed (DE) ncRNAs, DE mRNAs, and DE genes (DEGs) related to pain and analyzed the possible regulatory mechanisms of ncRNAs in trigeminal neuropathic pain.

Methods Animals

Male and female C57BL/6 mice (7 weeks) were used in the present study. All animal procedures performed in this study were reviewed and approved by the Research Ethics Committee of West China Hospital of Stomatology and were conducted according to the International Association for the Study of Pain guidelines. Animals were habituated for one week before experiments and maintained on a 12:12 light-dark cycle at a room temperature of 22 ± 1 °C with free access to food and water.

Trigeminal Neuropathic Pain Model Induced by CCI-ION

Twenty-four male mice were randomly divided into CCI-ION (n = 12) and sham (n = 12) groups. In addition, six female mice were randomly divided into CCI-ION (n = 3) and sham groups (n = 3). Chronic constriction injury of the infraorbital nerve (CCI-ION) was carried out to establish orofacial neuropathic pain.17 Animals were anesthetized with 4% pentobarbital sodium (50 mg/kg, i.p.) and laid on their back. Intraoral surgery was performed under aseptic conditions. Briefly, a 3-mm-long incision was made along the gingivobuccal margin in the buccal mucosa, beginning next to the maxillary first left molar. The distal segment of ION outside the infraorbital foramen was isolated from surrounding connective tissue and clearly visualized from the infraorbital foramen. Then, it was loosely tied with two chromic gut (4-0) ligatures at an interval of 2 mm. Finally, the mucosa incision was closed by 5-0 silk sutures. Sham-operated mice received only unilateral nerve exposure without ligature.

Behavioral Test

Nociceptive responses of Sham (n = 6) and CCI-ION (n = 6) male mice to mechanical stimuli were measured on 1 day before and 3, 7, 10, 14, 21, and 28 days after surgery. A series of 8 von Frey filaments (0.04 g, 0.16 g, 0.4 g, 0.6 g, 1 g, 1.4 g, 2 g, and 4 g, Stoelting Company, Wood Dale, USA) were applied to the left orofacial skin to examine the head withdrawal threshold (HWT) of mice in sham and CCI-ION groups in the up-down method described previously.18 The mouse was placed in a metal mesh cage (8 cm × 8 cm × 8 cm). After it calmed down, the hairs were presented at the left whisker pad, of which the hair was cut before the test. The starting filament was the middle of the series. Each filament was tested at an interval of several seconds. A positive response was noted if the head was sharply withdrawn. Once a negative response followed by a positive response or vice versa was approached, another four von Frey presentations were done according to the up-down rules. The mechanical pain threshold was quantified as EF50 (the response frequency), which was calculated by combining the value of the factor based on the response pattern.18 When continuous positive or negative responses were observed to the exhaustion of the stimulus set, values of 4 g and 0.04 g were assigned, respectively. The data of behavioral tests were presented as the means ± SD. The results of the Sham and CCI-ION groups on the same day were statistically analyzed using a t-test. The time tendency of behavioral testing data was detected with repeated measure ANOVA. Significance was established at P < 0.05.

RNA Extraction, Library Preparation, and RNA Sequencing

On 7 days after the operation, the left TG was collected only from six male mice (3 mice in sham group and 3 mice in CCI-ION group) after being deeply anesthetized with isoflurane. Total RNA was extracted from the isolated TG using TRIZOL (Invitrogen, Carlsbad) according to the manufacturer’s instructions. RNA degradation and contamination were monitored on 1% agarose gels. RNA purity and concentration were examined using the NanoPhotomet ® spectrophotometer (IMPLEN, CA, USA) and Qubit® RNA Assay Kit in Qubit® 2.0 Fluorometer (Life Technologies, CA, USA). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA). Library construction and sequencing were conducted by Novogene (Beijing, China). For lncRNA, circRNA, and mRNA sequencing, about 3 µg of total RNA from each sample was used to deplete ribosomal RNA according to the Epicenter Ribo-zeroTM rRNA Removal Kit (Epicenter, USA) following the manufacturer’s instruction. After the rRNA-free residue was cleaned up by ethanol precipitation, RNA was then fragmented into 250~300 bp fragments. Fragmented RNAs were reverse-transcribed to create the first-strand cDNA, which was used to synthesize U-labeled second strand DNAs with DNA polymerase I, RNase H, and dNTPs followed by adapter ligation and PCR amplification. For small RNA sequencing, about 3 μg total RNA per sample was used to generate a small RNA library using NEBNext® Multiplex Small RNA Library Prep Set for Illumina® (NEB, USA) following the manufacturer’s recommendations. Briefly, after 3’ end of small RNAs was ligated by NEB 3’ SR Adaptor, the SR RT Primers were used to transform the single-stranded DNA into a double-stranded DNA. Then, PCR amplification was performed using LongAmp Taq 2X Master Mix, SR Primer for illumina and index (X) primer. The RNA fragments from the above two libraries were finally sequenced on an Illumina NovaSeq 6000 platform, and 50 bp single-end reads for the small RNAs were generated, and 150 bp paired-end reads for the library of lncRNAs, circRNAs, and mRNAs were generated. The raw and clean reads number of the library for miRNAs for the samples were 13,077,884–16,706,212 and 12,891,853–16,477,632, respectively. The raw and clean reads number of the library for lncRNAs, circRNAs, and mRNAs for the samples were 79,399,034–105,087,144 and 77,718,888–103,727,040, respectively. The raw and clean bases of the library for miRNAs for the samples were, respectively, 0.654G–0.835G and 0.645G–0.824 G. The raw and clean bases of the library for lncRNAs, circRNAs, and mRNAs for the samples were, respectively, 11.91G–15.76G and 11.66G–15.56G. The error rate of the library for miRNAs and library of lncRNAs, circRNAs, and mRNAs for the samples were, respectively, <0.01 and <0.03%. The raw data of RNA-seq for miRNAs (SRA ID: SUB10791226, BioProject ID: PRJNA788681) and lncRNAs, circRNAs, and mRNAs (SRA ID: SUB10816303, BioProject: PRJNA788734) have been published in the NCBI SRA database.

Differential Analysis of ncRNAs and mRNAs

The small RNA tags were mapped to the reference sequence by Bowtie with 1 mismatch to analyze their expression and distribution. Clean reads with a length of 18–35 nt were screened as miRNA. miRBase 22 was used as a reference to look for the known miRNA. Modified software mirdeep2 and srna-tools-cli were used to obtain the potential miRNA and draw the secondary structures.19 Only miRNAs with Fold Change 1.4 and P value < 0.05 were classified as DE miRNAs. All transcriptomes from all samples were merged to reconstruct a comprehensive transcriptome using Perl scripts. Hisat220 was used to map reads to the genome of mice. LncRNAs were identified from the assembled transcripts following steps: (1) removal of short transcripts < 200 bp and < 2 exons; (2) removal of the transcripts with protein-coding capability determined using CNCI (Coding-Non-Coding-Index),21 Pfam and CPC (Coding Potential Calculator)22 database; Those remaining were considered to be lncRNAs. Novel lncRNAs were named following the rules of HGNC (The HUGO Gene Nomenclature Committee). StringTie23 was used to perform the expression level of mRNAs and lncRNA by calculating FPKM (fragments per kilobase of transcript per million mapped reads). The differential expression analysis was performed using the edgeR package (1.10.1). The resulting P values were adjusted using Benjamini and Hochberg’s approach for controlling the false discovery rate. The adjusted P value < 0.05 and |log2 (fold change)| > 1 were defined as the cutoff criteria for differentially expressed mRNAs (DE mRNAs) and DE Genes (DEGs) encoding proteins. The DE lncRNAs were selected with adjusted P-value < 0.05 and |log2 (fold change)| > 1. Find_circ and CIRI224,25 were used to de novo assemble the unmapped reads to circRNAs. The raw counts were first normalized using TPM (Read Count*1,000,000)/libsize). DE circRNAs were extracted with |log2 (fold change)| >1 and with statistical significance (P value < 0.05).

ncRNA Target Analysis

The target gene of miRNAs was performed by miRanda and RNAhybrid.26,27 The intersection result of the above two software was considered the target of miRNAs. The target mRNA analysis of DE lncRNAs was performed with two methods: co-location and co-expression analysis by bedtools and R package cor.test, respectively. Co-location analysis is based on the Cis-acting theory that lncRNAs are acting on the nearby protein-coding genes. The threshold of co-location was set within 100 kb upstream and downstream of lncRNAs. For co-expression analysis (Trans role of target gene prediction), target gene prediction is based on the expression correlation of lncRNAs and mRNAs. The absolute value of correlation was >0.95. For circRNA target analysis, some source genes of circRNAs can also be normally transcribed into linear mRNAs.28 Therefore, the expression level of circRNA may also affect the expression level of mRNA from the same source gene. The correlation analysis between the expression level of DE circRNA source genes and DE mRNA source genes was performed in this study. In order to predict the role of DE ncRNAs and DEGs with trigeminal neuropathic pain, the list of pain genes (Supplementary Table 1) was constructed based on PainNetworks,29 PainGene databases,30 and literature curation.31 The visualization of the regulatory network of DE ncRNAs and targeted DE mRNAs were represented by Cytoscape 3.9.0.

Go Annotations, KEGG Pathways, and PPI Analysis

Gene Ontology (GO) enrichment analysis and KEGG pathway analysis were performed to predict the function of DE ncRNAs and DEGs encoding proteins in the condition of trigeminal neuropathic pain. GO analysis was used to elucidate genetic regulatory networks of the DE ncRNAs and DEGs according to the biological process (BP), cellular component (CC), and molecular function (MF). GOseq R package-based Wallenius non-central hypergeometric distribution was implemented for GO enrichment analysis.32 The number of genes was statistically analyzed with significant enrichment of each GO term and displayed in a histogram. KEGG (Kyoto Encyclopedia of Genes and Genomes) is a database resource for understanding high-level functions and utilities of the biological system from molecular-level information.33 KOBAS software was used to test the statistical enrichment of the target gene candidates in KEGG pathways.34 A protein–protein interaction (PPI) network of DEGs was drawn using Cytoscape 3.9.0 based on the reported 1002 PPIs involved in pain.35 Single edges not connected to the main network were removed. The STRING was used to analyze the PPI of the DEGs related to the known pain genes. Nociceptor-related DEGs were identified in the whole TG, referencing the reported results of single-cell RNA-seq analysis of DEGs in DRG sensory neurons during SNT injury in mice.36

Construction and Analysis of ceRNAs Regulatory Network

To reveal the role and interactions among ncRNAs and mRNAs in trigeminal neuropathic pain pathogenesis, competing endogenous RNA (ceRNA) regulatory networks were constructed for DE lncRNAs, DE circRNAs, DE miRNAs, and DE mRNAs with miRanda software. The interaction network was built and visually displayed using Cytoscape 3.9.0 based on the screening of DE lncRNA-DE miRNA-DE mRNA pairs and DE circRNA-DE miRNA-DE mRNA pairs. In the ceRNA regulation networks, DE lncRNAs/DE circRNAs were determined as a decoy, DE miRNAs as the center, DE mRNAs as the target. The size of the node was in direct proportion to the degree.

Real-Time Quantitative Polymerase Chain Reaction (RT-qPCR)

Totally, six male mice (3 mice in sham group and 3 mice in CCI-ION group) and six female mice (3 mice in sham group and 3 mice in CCI-ION group) were used for testing the randomly chosen DE ncRNAs and DE mRNA in the TG by RT-qPCR. RNA Extraction Kit (Takara, China) was used to extract the total RNA from TGs. The RNA was reversely transcribed into the first-strand cDNA with Prime Script™ RT reagent Kit (Takara, China). The Mir-X miRNA First-Strand Synthesis Kit (Cat. No. 638313, Takara, China) was used for converting miRNAs. The sequences of polymerase chain reaction (PCR) primers (Tsingke Biotechnology, China) are shown in Supplementary Table 2. With an ABI QuantStudio 7 (Applied Biosystems, Waltham, MA), the PCR amplification was performed at parameters: Step 1 (95°C for 30s); Step 2 (40 cycles: 95°C for 5s, 60°C for 34s); Step 3 (Melt Curve). Calculation of the relative expression of lncRNAs and circRNAs was consistent with mRNA using the 2−ΔΔCt method. The internal control was Gapdh. The expression levels of miRNAs were normalized to the U6 level in each sample. The data of the behavioral test and RT-qPCR were presented as the means ± SD. The RT-qPCR results were analyzed using a t-test. Significance was established at P < 0.05.

Results CCI-ION Mice Suffer from Trigeminal Neuropathic Pain

The HWT of male Sham and CCI-ION mice were recorded on day 1 before and day 3, 7, 10, 14, 21, and 28 after surgery using the von Frey filament test. The behavioral test results showed that CCI-ION mice suffered from facial allodynia from day 3 after surgery compared with the Sham group (Figure 1). In addition, from day 7 after surgery, the HWT of CCI-ION mice reached a low level and lasted to at least day 28 after surgery. Moreover, compared with the baseline, there were no significant changes in the HWT of the Sham group. However, the HWT of CCI-ION mice began to decrease significantly from day 3 after surgery. These results indicated that CCI-ION treated mice suffered from trigeminal neuropathic pain from day 3 to day 28 after surgery.

Figure 1 Orofacial mechanical allodynia induced by CCI-ION. The results of von Frey filaments indicated that the head withdrawal threshold (HWT) of CCI-ION mice began to decrease from 3 days and lasted at least 28 days after surgery compared to the baseline (n = 6, *P < 0.05). CCI-ION mice almost reached the significantly lowest HWT on day 7 after surgery compared to the Sham group (n = 6, Sham vs CCI, 2.32 ± 0.47 g vs 0.41 ± 0.16 g, #P < 0.05).

CCI-ION Regulates the Expression of ncRNAs and mRNAs in the TG of Mice

Information on the DE miRNAs, DE lncRNAs, DE circRNAs, DE mRNAs, and DEGs in the TG of male CCI-ION mice compared with the male Sham group on 7 days after surgery are listed in Supplementary Tables 3–7, respectively. Totally, 1325 miRNAs were found in the TG of Sham and CCI-ION groups (Figure 2A). And 67 miRNAs (66 annotated miRNAs and 1 novel miRNAs) were differentially expressed compared with the Sham group (Figure 2B and C, Supplementary Table 3), including 36 up-regulation and 31 down-regulation (P value < 0.05, absolute log2FC > 0.485). For lncRNA, 43586 lncRNAs were detected in the TG of Sham and CCI-ION groups (Figure 3A). And 216 lncRNAs (165 annotated lncRNAs and 51 novel lncRNAs) were differentially expressed (adjusted P value < 0.05, absolute log2FC >1), including 144 up-regulation and 72 down-regulation (Figure 3B and C, Supplementary Table 4). For circRNAs, 2627 circRNAs were detected in Sham and CCI-ION groups (Figure 4A). And 14 circRNAs (5 annotated circRNAs and 9 novel circRNAs) were significantly regulated, including 9 up-regulation and 5 down-regulation (P value < 0.05, absolute log2FC >1) (Figure 4B and C, Supplementary Table 5). For mRNA analysis at the transcription level, 43695 mRNAs were detected in Sham and CCI-ION groups (Figure 5A). And 595 DE mRNA transcripts (adjusted P value < 0.05, absolute log2FC > 1) from 539 host genes were determined in the TG of CCI-ION mice compared with the Sham group, including 402 up-regulation and 193 downregulation (Figure 5B and C, Supplementary Table 6). For the analysis of DEGs encoding proteins, 18303 genes encoding proteins were determined, and 421 DEGs (adjusted P value < 0.05, absolute log2FC > 1, 388 up-regulation and 33 down-regulation) were determined in the TG of CCI-ION mice (Supplementary Figure 1A–C, Supplementary Table 7). Figure 6A shows the summary histogram of DE miRNAs, DE lncRNAs, DE circRNAs, DE mRNAs, and DEGs. To validate the reliability of the sequencing results, the change of ncRNAs and mRNAs expression 7 days after surgery were analyzed. One of DE miRNAs (mmu-miR-134-3p), DE lncRNAs (A630023A22Rik-202), DE circRNAs (mmu_circ_0000668), and DE mRNAs (Anxa10) were analyzed by RT-qPCR. All of the validated ncRNAs and mRNA were consistent with the results of RNA sequencing (Figure 6B and C). In addition, to investigate the changes in the expression of the above randomly chosen DE miRNA, DE lncRNA, DE circRNA, and DE mRNA in the TG after trigeminal nerve injury, RT-qPCR was also performed in female mice. The results showed that the changes in the expression of A630023A22Rik-202, mmu_circ_0000668, andAnxa10 mRNA were the same as that of male mice, while mmu-miR-134-3p has no significant change (Supplementary Figure 2A–D).

Figure 2 The expression profiling changes of miRNAs in the TG of CCI-ION mice (n = 3). (A) The Venn diagram shows the number of overlap miRNAs in the CCI-ION group compared with the Sham group. (B) The Volcano plot displays the DE miRNAs in the CCI-ION group compared with the Sham group. (C) The clustering map shows the top 20 up and down-regulated DE miRNAs in each sample TG of CCI-ION and Sham groups.

Figure 3 The expression profiling changes of lncRNAs in the TG of CCI-ION mice (n = 3). (A) The Venn diagram shows the number of overlap lncRNAs in the CCI-ION group compared with the Sham group. (B) The Volcano plot displays the DE lncRNAs in the CCI-ION group compared with the Sham group. (C) The clustering map shows the DE lncRNAs in each sample TG of CCI-ION and Sham groups.

Figure 4 The expression profiling changes of circRNAs in the TG of CCI-ION mice (n = 3). (A) The Venn diagram shows the number of overlap circRNAs in the CCI-ION group compared with the Sham group. (B) The Volcano plot displays the DE circRNAs in the CCI-ION group compared with the Sham group. (C) The clustering map shows the DE circRNAs in each sample TG of CCI-ION and Sham groups.

Figure 5 The expression profiling changes of mRNAs in the TG of CCI-ION mice (n = 3). (A) The Venn diagram shows the number of overlap mRNAs in the CCI-ION group compared with the Sham group. (B) The Volcano plot displays the DE mRNAs in the CCI-ION group compared with the Sham group. (C) The clustering map shows the DE mRNAs in each sample TG of CCI-ION and Sham groups.

Figure 6 Summary of DE ncRNAs, DE mRNAs, and DEGs in the TG of CCI-ION mice and RT-qPCR validation of RNA-seq results. (A) The histogram shows the number of down and up-regulated lncRNAs, circRNAs, miRNAs, and genes encoding mRNAs. (B) RT-qPCR validation of the expression level of randomly chosen DE miRNA, DE lncRNA, DE circRNA, and DE mRNA in the TG of male mice, *P< 0.05, n = 3. (C) The expression level of chosen DE miRNA, DE lncRNA, DE circRNA, and DE mRNA by RNA sequencing.

Some Pain-Related DE mRNAs Could Be Targeted by DE miRNAs and DE lncRNAs

ncRNAs can regulate the expression of mRNAs to exact biological functions. In order to predict DE mRNAs regulated by ncRNAs, Venn diagram analysis was performed between the targeted mRNAs of DE ncRNAs (miRNAs, lncRNAs, and circRNAs) and DE mRNAs (Figure 7A–C). Firstly, we searched the pain-related studies of the 67 DE miRNAs in PubMed and found that 18 DE miRNAs have been reported to be involved in pain (Supplementary Table 8), including neuropathic pain, inflammatory pain, diabetic neuropathic pain, and cancer pain. MiRNAs can degrade mRNAs or inhibit translation by acting on specific sequences in the 3` untranslated region (UTR) of target mRNAs.7 The intersection of the targeted mRNAs of DE miRNAs and DE mRNAs indicated that 125 DE mRNAs were the same as the targeted mRNAs of 37 DE miRNAs (Supplementary Table 9, Figure 7A). Among the 37 DE miRNAs, 9 DE miRNAs have been reported to be associated with pain (Supplementary Table 9). Meanwhile, among the 125 targeted DE mRNAs, 5 DE mRNAs belonging to the known pain genes (Cacna1a, Grm1, Il1rl1, Lmx1b, and Ctsd) were targeted by 6 DE miRNAs (Supplementary Table 9). Moreover, 5 proteins encoded by DE mRNAs (Cacna1a, Pik3cd, Il1rl1, Cga, and Spr) targeted by 5 DE miRNAs were involved with the 8 pain-specific PPI relations (Supplementary Table 10). Finally, the negative regulation network of DE miRNA and DE mRNA showed that 30 DE miRNAs were predicted to negatively target 80 DE mRNAs (Figure 8, Supplementary Table 11). In the network, 7 DE miRNAs have been reported to be associated with pain, and 4 DE mRNAs (Cacna1a, Il1rl1, Lmx1b, and Ctsd) were related to the known pain genes.

Figure 7 Intersection analysis of DE mRNAs and DE ncRNAs targeting mRNAs. (A) Venn diagram showing the overlap number of the targeted mRNAs of DE miRNAs and DE mRNAs. (B) Venn diagram displaying the overlap number of co-expression and co-location targeted mRNAs of DE lncRNAs with DE mRNAs. (C) The Venn diagram shows the overlap number of targeted mRNAs of DE circRNAs and DE mRNAs.

Figure 8 The regulatory network of DE miRNAs negatively targeting DE mRNAs. Totally, 30 DE miRNAs were predicted to target 80 DE mRNAs negatively. In the network, 7 DE miRNAs (circles with black border) have been associated with pain, and 4 DE mRNAs (triangles with black border) were related to the known pain genes. The node size is directly proportional to the number of their corresponding partners.

The Venn diagram shows the interaction between DE lncRNA targeted mRNAs by co-location or co-expression with DE mRNAs (Figure 7B). Co-expression analysis demonstrated that 139 DE lncRNAs were predicted to target 488 DE mRNAs (82% DE mRNAs). Figure 9 shows the regulatory network of 107 DE lncRNAs targeting 35 DE mRNAs related to known pain genes according to the list of pain genes constructed based on PainNetworks, PainGene databases, and previous literature.29–31 The detailed information on the 107 DE lncRNAs and 35 pain-related DE mRNAs is listed in Supplementary Table 12. Co-location analysis suggested that 24 DE lncRNAs were predicted to target 29 DE mRNAs (Supplementary Table 13, Supplementary Figure 3), but no DE mRNA was related to the known pain genes. Furthermore, 22 DE mRNAs were predicted to be regulated by DE lncRNAs both through co-expression and co-location. However, only 3 DE mRNAs (Ly6a-202, Lrp8-201, and Esd-202) were predicted to be regulated by the same DE lncRNAs (Gm20275-201, 4933424M12Rik-202, and TCONS_00087827) through co-expression and co-location, respectively (Supplementary Figure 4), which are more associated with trigeminal neuropathic pain. CircRNAs may regulate the expression of mRNAs from the same host gene.28 Intersection analysis of DE circRNA host genes and DE mRNA host genes indicated the potential regulatory relationship between DE circRNAs and DE mRNAs (Figure 7C). Notably, 2 DE circRNAs (mmu_circ_0000016 and mmu_circ_0000041) were from the same source genes as 2 DE mRNAs (Rims1 and Hecw2), respectively.

Figure 9 The regulatory network of DE lncRNAs targeting DE mRNAs related to the known pain genes by co-expression. Totally, 107 DE lncRNAs were predicted to target 35 DE mRNAs associated with pain by co-expression. The size of the node is in direct proportion to the degree.

Nociception-Related Signaling Pathways are Enriched for DE ncRNAs Targeting Genes

GO and KEGG pathway analyses were performed to predict the gene function of DE ncRNAs in the condition of trigeminal neuropathic pain. According to the distribution of predicted target genes in the Gene Ontology, the target gene function of DE ncRNAs was clarified. According to the GO analysis of target genes of DE miRNAs, the most significant biological processes (BP) were cellular component organization or biogenesis, metabolic process, organelle organization, etc. The most significantly enriched cellular components (CC) were related to the intracellular part, cytoplasm, organelle, etc. The most significantly enriched molecular function (MF) were protein binding, binding, catalytic activity, etc. (Figure 10A). Based on the GO analysis of the co-localization and co-expression genes of DE lncRNAs, the most significantly enriched GO terms were similar to the enriched GO terms of target genes of DE miRNAs. Significantly activated BP of these genes included cellular component organization, cellular component organization or biogenesis, positive regulation of the biological process, etc. The most noteworthy enriched CC was related to the intracellular part, cytoplasm, organelle, etc. The most significantly enriched MF were protein binding, binding, ion binding, etc. (Figure 10B). Based on the GO analysis of predicted target genes of DE circRNAs, the most significantly enriched BP were positive regulation of peptidyl-threonine phosphorylation, regulation of peptidyl-threonine phosphorylation, and intracellular protein transport, etc. The most significantly enriched CC was mitotic spindle, adherens junction, anchoring junction, etc. The most significantly enriched MF were protein binding, ankyrin repeat binding, acetyltransferase activator activity, etc. (Figure 10C).

Figure 10 GO analysis of the target genes of DE miRNAs, DE lncRNAs, and DE circRNAs. The top 20 significant items in the biological process (BP), cellular component (CC), molecular function (MF) by GO analysis of the target genes of DE miRNAs (A), DE lncRNAs (B), and DE circRNAs (C).

KEGG analysis of predicted genes of DE miRNAs in the TG of CCI-ION mice demonstrated the Rap1 signaling pathway, oxytocin signaling pathway, cAMP signaling pathway, MAPK signaling pathway, and lysosome were significantly enriched (Figure 11A), which were related to hyperalgesia.37,38 The KEGG analysis of predicted genes of DE lncRNAs in the TG of CCI-ION mice indicated that the most significantly enriched pathways were focal adhesion. It was noteworthy that MAPK signaling pathway, PI3K-Akt signaling pathway, ECM-receptor interaction, and gap junction were also significantly enriched (Figure 11B), which were involved in the process of nociceptive transmission. For the host genes of 9 DE circRNAs determined in the TG of CCI-ION mice, only retrograde endocannabinoid signaling pathways were significantly enriched (Supplementary Figure 5).

Figure 11 KEGG analysis of the target genes of DE lncRNAs and DE miRNAs. (A) The scatter dot graph shows the top 20 enriched pathways related to the target genes of DE miRNAs; (B) The scatter dot graph shows the top 20 enriched pathways related to the target genes of DE lncRNAs (P < 0.05).

Some DE Genes Encoding Proteins are Related to Pain

Among the DEGs encoding proteins (421 genes from Supplementary Table 7) in the TG of CCI-ION mice, Barx1 was expressed specifically in the TG compared to the DRG.39 In order to characterize the changes in the functional properties of DEGs, GO, KEGG pathway, and PPI analysis were performed. Based on the GO analysis, the most noteworthy enriched BP were defense response, immune system process, response to external stimulus, and response to stress. The extracellular region, extracellular space, extracellular region part, and plasma membrane part in the category of the CC were enriched. The most noteworthy enriched MF were binding, ion binding, anion binding, and carbohydrate derivative binding (Supplementary Figure 6A–D). In the KEGG analysis of DEGs, the most significantly involved pathways in CCI-ION pathogenesis were PI3K−Akt signaling pathway, focal adhesion, ECM−receptor interaction, and protein digestion and absorption (Supplementary Figure 7).

To explore the DEGs in the TG induced by CCI-ION involved in the pain-related PPI network, the previously published pain interactome, a PPI network of 611 interconnected pain-related proteins, was used.35 Totally, 30 up-regulated DEGs and 66 partner genes were involved in the pain-specific PPI network without down-regulated DEGs, including 107 PPIs (Figure 12, Supplementary Table 14). Several up-regulated hub genes (Prkcg, Grin2b, Gal, Npy, Vip, and Mmp9) were obviously centered around by other genes. These genes were related to nerve injury and regeneration (Gal and Npy),40,41 pro-nociceptive signaling transmission (Prckg, Grin2b, Vip, and Mmp9),42–45 and anti-nociceptive signaling transmission (Npy).46 To further determine the DEGs related to the known pain genes in the TG of CCI-ION mice, an intersection analysis of DEGs and the known 804 pain genes (Supplementary Table 1) according to previous databases and literature was performed.29–31 The results showed that 39 DEGs (38 up-regulation and 1 down-regulation) were known pain genes (eg, top 10 up regulated genes: Prok2, Sprr1a, Vip, Npy, Slc6a4, Mpo, Kcnj5, Col9a1, Lcn2, Grin2b, and the down-regulated Nmu) (Supplementary Table 15). The STRING network analysis of the 39 DEGs related to the known pain genes showed an interconnected network of 23 DEGs and 5 partner genes (Tnf, Cd79a, Cyba, Ncf1, and Dlg4) (Figure 13). The cluster of Atf3, Lcn2, Mmp9, Mpo, Serpinb1a, Spp1, Tnf, and Tnfrsf1b was significantly enriched in the IL-17 signaling pathway, TNF signaling pathway, and immune system. The cluster of Gal, Npy, Nts, and Vip was enriched in neuropeptide signaling pathway and peptide hormone binding. The cluster of Dlg4, Grin2a, Grin2b, and Shisa9 was enriched in glutamate-gated calcium ion channel activity and long-term potentiation. The cluster of Cd79a, Itpr1, Lyn, and Prkcg was significantly enriched in long-term depression and potentiation. The cluster of Cckar and Cckbr was implicated in peptide hormone binding and calcium signaling pathway. Moreover, the network included nerve injury-related genes Atf3, Tnf, Mmp9, and Anxa1. Finally, according to the lists of genes related to small nonpeptidergic (NP) nociceptors, medium peptidergic (PEP) nociceptors, and large myelinated (LM) sensory neurons,36 73 DEGs (17% DEGs) were related to nociceptor (NP or PEP), 72 DEGs were related to NP nociceptors, 26 DEGs were related to medium PEP nociceptors, and 18 DEGs were related to LM sensory neurons (Supplementary Table 16).

Figure 12 The pain-related protein–protein intersection (PPI) network of the DEGs encoding proteins in the TG of CCI-ION mice. The network was built based on the pain-specific PPI network. Using DEGs in the TG of CCI-ION mice and partner genes from the pain-specific PPI network (DEGs and non-DEGs). The colored arrows represent the type and direction of interactions; The colored nodes mark the expression changes. Totally, 30 up-regulated genes were involved in the pain-specific PPI network, and no down-regulated DEGs were involved in the pain-specific PPI network. DEG, differentially expressed gene; PPI, protein–protein interaction. The size of the node is in direct proportion to the degree.

Figure 13 PPI analysis of biological interactions of the pain-related DEGs. The analysis was performed with STRING, using high confidence of interactions (score 0.7) and not more than 5 interactors in the first shell. The nodes were clustered using the MCL algorithm. The disconnected genes in the network are not shown.

CeRNA Regulatory Networks Involved with Pain Exist in the TG of CCI-ION Mice

It is well known that lncRNAs and circRNAs can sponge miRNAs to regulate the expression of mRNAs, which is proven to be involved in the development and maintenance of neuropathic pain.7,8 Therefore, to further study the molecular mechanisms in ncRNAs involved in trigeminal neuropathic pain, the lncRNA/circRNA-miRNA-mRNA regulatory networks based on the ceRNA regulatory theory were analyzed. The ceRNA regulation network of DE lncRNAs, DE miRNAs, and DE mRNAs demonstrated that 42 down-regulated lncRNAs were predicted to regulate the expression of 21 down-regulated mRNAs by sponging 15 up-regulated miRNAs in the TG of CCI-ION mice (Figure 14A). Among the DE mRNAs in the above down-up-down regulatory network, Ctsd was related to the known pain gene. The network in Figure 14B displayed that 96 up-regulated lncRNAs were predicted to sponge 15 down-regulated miRNA to regulate the expression of 59 up-regulated mRNAs in the TG of CCI-ION mice. Among the DE mRNAs in the above up-down-up regulatory network, 3 up-regulated mRNAs (Cacna1a, Il1rl1, and Lmx1b) were related to the known pain genes.29–31 In Figure 15, the ceRNA regulatory network of DE circRNAs, DE miRNAs, and DE mRNAs showed that 3 down-regulated circRNAs were predicted to regulate the expression of 45 up-regulated mRNAs by sponging 5 down-regulated miRNAs. Besides, 6 up-regulated circRNAs were predicted to regulate the expression of 20 down-regulated mRNAs through sponging 11 up-regulated miRNAs. In the above ceRNA regulation network of DE circRNA-DE miRNA-DE mRNA, the up-regulated Ctsd and the down-regulated Cacna1a and Lmx1b were associated with the known pain genes.29–31

Figure 14 Regulatory network analysis of DE lncRNA-DE miRNA-DE mRNA in the TG of CCI-ION mice based on the ceRNA regulation theory. (A) The down-regulated DE lncRNAs were predicted to regulate the expression of the down-regulated DE mRNAs by sponging the up-regulated DE miRNAs in the TG of CCI-ION mice. The triangle with a red border (Ctsd) represents the DE mRNA related to pain. (B) The up-regulated lncRNAs were predicted to regulate the expression of the up-regulated mRNAs by sponging the down-regulated miRNAs in the TG of CCI-ION mice, the triangle with a green border (Cacna1a, Il1rl1, and Lmx1b) represents the DE mRNA related to pain. The size of the node is in direct proportion to the degree.

Figure 15 Regulatory network analysis of circRNA-miRNA-mRNA in the TG of CCI-ION mice according to the ceRNA mechanism. The up-regulated circRNAs and down-regulated circRNAs were predicted to regulate the expression of DE mRNAs by sponging DE miRNAs in the TG of CCI-ION mice. The triangles with a black border represent the DE mRNAs related to the known pain genes. The size of the node is in direct proportion to the degree.

Based on the GO analysis of target mRNAs in the ceRNA regulatory network of DE lncRNAs, DE miRNAs, and DE mRNAs, the noteworthy enriched BP were macromolecule localization, secretion, and regulation of cellular localization, etc. The most noteworthy enriched CC were cell projection, membrane region, and dendrite. The most noteworthy enriched MF were binding, protein binding, ion binding, etc. (Supplementary Figure 8A). Based on the GO analysis of target mRNAs of the regulatory network of DE circRNAs, DE miRNAs, and DE mRNAs according to the theory of ceRNA, the most significantly enriched BP were cell–cell junction organization, regulation of neuron differentiation, and regulation of neurogenesis, etc. The most noteworthy enriched CC were cell projection, cell projection part, and dendrite, which were associated with neurons. The most noteworthy MF were binding, protein binding, and ion binding (Supplementary Figure 8B). When the data were KEGG-analyzed with the predicted mRNAs in the DE lncRNA-DE miRNA-DE mRNA regulatory network based on the theory of ceRNA under the pathogenesis of CCI-ION induced neuropathic pain, the most noteworthy involved pathways were retrograde endocannabinoid signaling, progesterone-mediated oocyte maturation, and aldosterone-regulated sodium reabsorption (Supplementary Figure 9A). The KEGG analysis of the DE circRNA-DE miRNA-DE mRNA regulatory network indicated that the most noteworthy KEGG pathways were retrograde endocannabinoid signaling, cholinergic synapse, carbohydrate digestion and absorption, and synaptic vesicle cycle (Supplementary Figure 9B). These results illustrate the regulatory relationship between ncRNAs and mRNAs in the mechanism of trigeminal neuropathic pain. However, the regulatory role of ncRNAs in trigeminal neuropathic pain was complicated. It needs studies furthermore to elucidate their role in the development and maintenance of trigeminal neuropathic pain.

Discussion

Trigeminal neuropathic pain belongs to orofacial pain that is attributed to lesions or diseases of cranial nerves according to the International Classification of Orofacial Pain (ICOP), 1st edition.1 However, highly effective treatments or drugs for trigeminal neuropathic pain are still unavailable because of the lack of a complete understanding of its mechanisms.3,6 Although increasing data demonstrate that ncRNAs (miRNAs, lncRNAs, and circRNAs) play crucial roles in neuropathic pain, little is known about the expression profile and roles of ncRNAs in the TG under the condition of trigeminal neuropathic pain. As far as we know, this study was the first to study the expression profile of miRNAs, lncRNAs, and circRNAs in the TG of CCI-ION mice (a classical trigeminal neuropathic pain model) by whole-transcriptome sequencing and bioinformatics analysis. The results of this study support that miRNA, lncRNA, and circRNA contribute to the development of trigeminal neuropathic pain.

miRNA is a group of small non-coding RNAs consisting of 18~23 nucleotide sequences. Since miRNA was proved to be involved in orofacial inflammatory pain through regulating the cytokine mRNA expression and had an analgesic action at the spinal cord level towards an anti-inflammatory phenotype,47–49 the role of miRNAs in pain regulation has attracted increasing attention. In this study, 18 DE miRNAs have been reported to be involved in pain (Supplementary Table 8). miRNA with the “activated” RNA-induced silencing complex (RISC) attaches to the 3’ UTR of a target mRNA, resulting in the degradation of target mRNAs (full complementarity) or translational repression (incomplete complementarity).7 According to the above mechanism, 37 DE miRNAs were predicted to target 125 DE mRNAs negatively, and 4 DE target mRNAs (Ctsd, Lmx1b, Il1rl1, and Cacnala) of DE miRNAs were related to the known pain genes. These results suggest that miRNAs may play a role in trigeminal neuropathic pain by targeting the known pain genes. Furthermore, several nociception-related signaling pathways were significantly enriched in the KEGG analysis, including Rap1 signaling pathway, oxytocin signaling pathway, cAMP signaling pathway, MAPK signaling pathway, lysosome, etc.37,50,51 For instance, protein kinase A is involved in neuropathic pain by activating the p38 MAPK pathway to mediate spinal cord cell apoptosis.52 The above results provide a theoretical basis for using miRNAs as therapeutic targets in trigeminal neuropathic pain.

LncRNAs with lengths exceeding 200 nucleotides are involved in many physiological and pathological processes including embryonic development, cancer, inflammation, cardiovascular and neurobiological diseases, etc.9,53,54 Recently, the role of lncRNAs has been quickly attracted in neuropathic pain.55 For instance, Zhou et al discovered 126 DE lncRNAs in the spinal cord of rats using RAN-sequencing 7 days after SNI surgery.56 Most recently, evidence also has shown that lncRNAs (Gm14461 and uc.48) in the TG play crucial roles in trigeminal neuralgia.14,57 In the present study, 216 DE lncRNAs were determined in the TG of mice treated with CCI-ION. LncRNAs can modulate gene expression through trans (co-expression) and cis (co-location) regulation.58 Intersection analysis of the DE mRNAs and the target mRNAs of DE lncRNA by co-expression demonstrated that 139 DE lncRNAs were predicted to target 488 DE mRNAs, including 35 DE mRNAs transcribed from the known pain genes.29–31 Furthermore, 3 DE mRNAs (Ly6a, Lrp8, and Esd) were predicted to be regulated by the same DE lncRNAs (Gm20275-201, 4933424M12Rik-202, and TCONS_00087827) through co-expression and co-location, respectively. Notably, Lrp8 is involved in the reelin signaling pathways that contribute to nociceptive processing.59 KEGG analysis suggested MAPK signaling pathway, PI3K-Akt signaling pathway, ECM-receptor interaction, and gap junction were enriched, which were related to neuropathic pain.52,60,61

CircRNAs are non-coding RNAs with covalently closed circular structures with 5’ and 3’ ends. The sequences of circRNAs are conserved among mammalian species and have higher stability than liner mRNAs because the particular ring structure enables them to escape from RNase degradation.28,62 Dysregulation of circRNAs is associated with tumors, cardiovascular diseases, endocrine diseases, immune system diseases, and neurodegenerative diseases.28,62 Recently, emerging studies found that circRNAs play important roles in developing and maintaining neuropathic pain.13,63 For instance, Zhang et al found 22 DE circRNAs in rat spinal cord on days 7 and 14 following SNL. And they proved that the up-regulated CircAnks1a in the spinal cord of rats contributes to the hypersensitivity in the paw induced by SNL injury.13 Besides, Zhou et al determined 188 DE circRNAs in the rat spinal cord on day 14 after SNI was controlled to the sham group with RNA-seq and bioinformatic analysis.10 In the present study, 9 DE circRNAs were determined in the TG of CCI-ION mice. Intersectional analysis showed that 2 DE circRNAs (mmu_circ_0000016 and mmu_circ_0000041) were from the same source genes as 2 DE mRNAs (Rims1 and Hecw2), respectively. RIMS1 can regulate neurotransmitter release at the presynaptic active zones.64,65 Besides, mmu_circ_0000041 (circRNA HECW2) engaged in the endothelial-mesenchymal transition related to the damage of the blood–brain barrier.66 KEGG analysis showed that only retrograde endocannabinoid signaling pathways were significantly enriched for the target genes of DE circRNAs. Tissue injury can activate retrograde endocannabinoid signaling,67,68 which is involved in multiple pathophysiological functions, such as anti-nociception, cognition and memory, endocrine function, nausea, and vomiting, inflammation, immune recognition, and modulation of synaptic transmission and plasticity.69,70

The best-known mechanism underlying the functional regulation of ncRNAs is that lncRNA/circRNA acts as a ceRNA to sponge miRNAs, reducing the number of available miRNAs and helping to improve the translation of target mRNAs.8 For example, increased circAnks1a sponged miR-324-3p to promote the VEGFB expression in the spinal cord of rats after SNL, resulting in neuropathic pain.13 Among the target DE mRNAs in the ceRNA regulatory network of DE lncRNAs, Ctsd, Cacna1a, Il1rl1, and Lmx1b are related to the known pain genes.29–31 Among the target DE mRNAs in the ceRNA regulatory network of DE circRNA, the up-regulated Ctsd and the down-regulated Cacna1a and Lmx1b were associated with the known pain-specific genes.29–31 These results support that lncRNAs/circRNAs in the TG of CCI-ION mice may play an essential role in the development of trigeminal neuropathic pain through ceRNA regulation.

Trigeminal nerve damage or compression of the trigeminal nerve root could cause apparent changes in the expression of various genes in the TG, which contribute to the initiation and sustainability of trigeminal neuropathic pain.3,6 In the present study, a total of 421 DEGs encoding proteins were found in the TG of mice on 7 days following CCI-ION surgery. KEGG pathway analysis showed PI3K−Akt signaling pathway, focal adhesion, ECM−receptor interaction, protein digestion and absorption and so on were significantly enriched, which was similar to the results of pathway enrichment of DEGs in the spinal cord of rats treated with SNI.10 Based on a previously reported pain-related PPI network,35 PPI analysis indicated several up-regulated hub genes (Prkcg, Grin2b, Gal, Npy, Vip, and Mmp9) were involved in the pain-specific PPI network. These genes were related to nerve injury and regeneration and nociception regulation (Gal and Npy),40,41 pro-nociceptive signaling transmission (Prckg, Grin2b, Vip, and Mmp9),42–45 and anti-nociceptive signaling transmission (Npy).46 In addition, 39 DEGs were known pain-specific genes.29–31 TG mainly consists of trigeminal ganglion neurons (TGNs), satellite glial cells (SGCs), Schwann cells, and macrophages.71 In the present study, we examined the expression profile of ncRNAs and mRNAs in the whole TG while not determining the changes in the expression of ncRNAs and mRNAs in different types of cells in the TG. In order to predict the cell-type origin of the identified DEGs, previously published data of single-cell RNA-seq analysis of DEGs in the mouse nonpeptidergic nociceptors (NP), peptidergic nociceptors (PEP), and large myelinated (LM) sensory neurons under the st