Global Analysis of mRNA Alternative Splicing in the Trigeminal Ganglion at Different Stages of Trigeminal Neuropathic Pain in Mice

Introduction

Trigeminal neuropathic pain (TNP) is a debilitating chronic pain disorder that arises from injury or dysfunction of the trigeminal nerve. It occurs more frequently in females than in males and significantly impairs quality of life, with an estimated rate of 5.5 cases per 100,000 individuals.¹ Clinically, TNP is challenging to manage because of its refractory nature and limited responsiveness to conventional analgesics. Furthermore, current therapeutic strategies often yield suboptimal efficacy and are frequently associated with adverse side effects. These limitations underscore the urgent need to elucidate the molecular mechanisms underlying TNP and to identify more precise therapeutic targets.

Emerging evidence indicates that gene expression changes in the trigeminal ganglion (TG) are essential for the development and persistence of neuropathic pain.² Although trigeminal and spinal sensory systems share many nociceptive features, they are not fully equivalent in their molecular organization and injury responses. These distinctions suggest that splicing-regulatory mechanisms identified in spinal neuropathic pain may not be fully generalizable to trigeminal pain and support the need for TG-specific investigation of post-injury molecular regulation.^3,4 Alternative splicing (AS) is a ubiquitous and tightly regulated post-transcriptional process that enables a single gene to produce multiple mRNA isoforms, thereby enhancing transcriptomic and proteomic diversity.^5,6 It is estimated that more than 95% of multi-exon genes in the human genome undergo AS. This mechanism influences protein structure and function and plays vital roles in tissue-specific gene expression, developmental timing, and disease pathogenesis.^7,8 In the nervous system, AS can influence neuronal excitability, ion channel function, synaptic plasticity, and neuroinflammatory signaling.⁹ Notably, alternative splicing of NaV1.7 has been shown to alter the functional consequences of a painful channel mutation, supporting a direct link between AS and pain-related ion channel regulation, while voltage-gated sodium channels more broadly are central determinants of sensory neuron excitability and pain signaling.^10,11 These observations suggest that stage-dependent AS may influence trigeminal nociceptor excitability and downstream pain signaling after nerve injury.

Although AS has been widely implicated in neurological disease, its role in pathological pain, particularly within the trigeminal system, remains poorly defined.^12–17 To date, no study has comprehensively characterized TG-specific, stage-resolved AS dynamics during the progression of TNP. Recent behavioral studies have shown that pain responses in the pIONT model vary significantly across postoperative time points, providing a temporal framework for linking pain phenotypes to underlying molecular changes.¹⁸ These findings raise the possibility that AS is dynamically regulated during the onset and maintenance phases of TNP rather than remaining static throughout disease progression.

In this study, we addressed this gap by profiling genome-wide AS events in a well-established partial infraorbital nerve transection (pIONT) mouse model of TNP. Poly(A)-enriched RNA sequencing was employed to analyze TG samples collected at two critical stages: day 3 (onset phase) and day 10 (maintenance phase) following surgery. These stages correspond to distinct phases of pain behavior observed in the von Frey test, with significant pain responses emerging at day 10, as documented previously.¹⁸ Because AS-related changes may not be fully captured at a single molecular level, transcriptomic findings were further integrated with TMT-based quantitative mass spectrometry proteomic analysis, RNA-binding protein (RBP) analysis, functional enrichment, and single-cell RNA-seq-based cell-type mapping to define the temporal, molecular, and cellular context of TNP-associated splicing alterations. Collectively, this integrative strategy provides a TG-focused and stage-resolved framework for understanding post-transcriptional regulation in trigeminal neuropathic pain.

Materials and Methods

Animals and Surgery

Eight-week-old male ICR mice were obtained from the Experimental Animal Center of Nantong University. Animals were housed in a temperature-controlled environment (22 ± 1°C) under a 12-h light/dark cycle, with ad libitum access to standard chow and water. All animal procedures were reviewed and approved by the Institutional Animal Care and Use Committee of Nantong University Experimental Animal Center (Approval No. S20251027-007) and were conducted in accordance with the Standard Operating Procedures for Laboratory Animal Center of Nantong University. The principles of Replacement, Reduction, and Refinement (3Rs) were followed wherever possible. Animals were monitored daily for signs of pain or distress, and humane endpoints were applied when necessary. For animals not used for terminal perfusion-based tissue collection, euthanasia was performed by CO₂ inhalation according to the approved institutional protocol.

The pIONT model was established as previously described.^18,19 Anesthesia was induced by intraperitoneal injection of pentobarbital sodium (60 mg/kg). With the animal in the supine position, the tongue was gently retracted to expose the upper left inner gingiva. A 5-mm incision was made anterior to the white tendon to expose the infraorbital nerve (ION), which was carefully isolated using a glass electrode. The ION was ligated with an 8–0 nylon suture and transected distal to the ligation site. In sham-operated controls, the ION was exposed but neither ligated nor transected. After surgery, animals were monitored closely during recovery in a temperature-controlled chamber for at least 24 h.

RNA Extraction and Sequencing

For RNA-seq experiments, forty mice were subjected to pIONT and forty mice underwent sham surgery. TGs were collected at 3 and 10 days after surgery. For terminal tissue collection, mice were deeply anesthetized with 5% isoflurane until loss of reflexes, followed by transcardial perfusion with 0.9% saline for approximately 3 min until the liver turned pale. TGs were then rapidly dissected and collected for RNA extraction.

Total RNA was extracted from TG tissues using TRIzol reagent (Tiangen, Beijing, China) according to the manufacturer’s instructions. RNA concentration and purity were measured using a Qubit fluorometer (Invitrogen, Carlsbad, CA, USA), and RNA integrity was assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Only samples with an RNA integrity number (RIN) > 7.0 and a 28S:18S rRNA ratio > 1.8 were used for subsequent analyses.

TG tissues from ten mice under the same experimental condition were pooled to generate one biological replicate. Thus, each condition at each time point (pIONT day 3, sham day 3, pIONT day 10, and sham day 10) yielded two pooled biological replicates, resulting in a total of eight RNA-seq samples. Because individual TG samples yielded limited amounts of RNA, pooling was used to ensure sufficient input for library preparation. Accordingly, the RNA-seq variance estimates reflect variation between pooled biological replicates rather than inter-individual biological variability. The day 3 and day 10 time points were selected to represent the onset and maintenance phases of neuropathic pain, respectively, based on prior behavioral characterization of the pIONT model. This pooled design was adopted for discovery-oriented transcriptomic profiling. Behavioral outcomes were not re-assessed in the current molecular cohort; therefore, associations between molecular alterations and pain-related outcomes were interpreted in the context of this previously established behavioral staging of the pIONT model.¹⁸

Poly(A)-enriched RNA libraries were prepared and sequenced on an Illumina platform using 150-bp paired-end reads. Raw sequencing data have been deposited in the Sequence Read Archive under accession number PRJNA814329.

Identification of Alternative Splicing Events

AS events were identified using an integrated computational pipeline combining ASprofile and rMATS (v4.1.0). ASprofile was used for global AS event classification, whereas rMATS was used for differential AS event detection. Transcriptome assembly was first performed using StringTie to reconstruct transcript annotations from the RNA-seq data.²⁰ rMATS was then used to detect and quantify five major types of AS events: skipped exon (SE), alternative 5′ splice site (A5SS), alternative 3′ splice site (A3SS), mutually exclusive exons (MXE), and retained intron (RI).²¹ The inclusion level was calculated by rMATS based on the proportion of sequencing reads mapped to transcript isoforms with or without the alternatively spliced exon, as schematically illustrated in Figure 1A. Differences in mean inclusion levels between pooled biological replicates in the experimental and sham groups were evaluated using a likelihood ratio test to identify differential splicing events. Differential AS events were defined as those with an inclusion level difference (|IncLevel difference|) > 0.05 and a false discovery rate (FDR) < 0.01.

Figure 1 Overview of alternative splicing (AS) events during trigeminal neuropathic pain (TNP). (A) Schematic illustration of the five major rMATS-defined AS event types analyzed in this study: skipped exon (SE), alternative 3′ splice site (A3SS), alternative 5′ splice site (A5SS), mutually exclusive exon (MXE), and retained intron (RI). (B) Distribution of AS events in individual TG RNA-seq samples and in the merged dataset based on ASprofile classification, including transcription start site-associated events (TSS), transcription termination site-associated events (TTS), alternative exon events (AE), exon skipping events (SKIP), and intron retention events (IR). (C) Number of differential alternative splicing events at days 3 and 10 after pIONT, defined using the criteria |IncLevel difference| > 0.05 and FDR < 0.01. (D–F) Representative sashimi plots for Kcnq2 (D), Fgf1 (E), and Kcnc3 (F), showing exon inclusion or exclusion patterns in the pIONT group (red) and sham group (orange). (G–I) qPCR validation of splice variant-related expression changes in Kcnq2 (G), Fgf1 (H), and Kcnc3 (I) in TG samples from sham and pIONT mice. (J–L) Agarose gel electrophoresis of PCR products for Kcnq2 (J), Fgf1 (K), and Kcnc3 (L), showing bands of the expected sizes corresponding to distinct transcript isoforms. Data are presented as mean ± SEM. *P < 0.05, **P < 0.01 versus sham, n = 8, two-tailed unpaired Student’s t-test.

Sequencing Data Processing and Analysis

The genome-wide distribution of AS events was visualized using Circos.²² Differentially expressed alternatively spliced genes (ASGs) and proteins were analyzed using Enrichr (https://maayanlab.cloud/Enrichr/) for pathway annotation and Gene Ontology (GO) enrichment. For GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses, P values were adjusted for multiple testing using the Benjamini–Hochberg method, and adjusted P values/FDR were used to determine significance. GO biological process and pathway enrichment scores were represented as –log10-transformed adjusted P values. Hierarchical and k-means clustering were performed using Cluster 3.0 (Stanford University), and the results were visualized with Java Treeview. Protein–protein interaction analysis was performed using STRING v11.5 with a medium confidence threshold (score > 0.4) using the Mus musculus protein interaction background. In STRING networks, edge thickness indicates interaction confidence, with thicker lines denoting higher confidence. Single-cell RNA-seq-based cell-type mapping was performed using the online portal tg.painseq.com, which is based on R Shiny and ShinyCell.²³ Ligand–receptor interaction analysis was conducted using CellLinker and the R package SingleCellSignalR (v0.0.1.8) to investigate potential signaling pathways involved in TNP. Multi-omics integration was performed as a descriptive cross-platform comparison by intersecting ASGs, differentially expressed genes (DEGs), differentially expressed proteins (DEPs), and RBP genes harboring differential AS events based on shared gene symbols and time-point-specific changes. No formal joint modeling framework was applied. Cross-platform overlaps were evaluated separately at day 3 and day 10 by matching ASGs, DEGs, DEPs, and RBP genes harboring differential AS events using shared official gene symbols.

Quantitative Real-Time PCR

The expression of splice variants from selected genes was validated by quantitative real-time PCR (qPCR). Total RNA was extracted from TG samples according to the manufacturer’s instructions, and 1 μg of RNA was reverse-transcribed into complementary DNA (cDNA) using a Takara cDNA synthesis kit. For each selected gene, two primer pairs were designed to distinguish transcript isoforms with or without the skipped exon. qPCR was performed to quantify splice variants of Kcnq2, Fgf1, and Kcnc3, with Gapdh used as the internal control. Primer sequences are listed in Table 1. Amplification reactions were carried out on a StepOne system (Corbett Life Science) using SYBR Premix Ex Taq II (Takara). Melting curve analysis was performed to confirm amplification specificity. Relative expression levels were calculated using the 2^−ΔΔCt method. For qPCR validation, n = 8 independent biological samples per group were analyzed.

Table 1 Primer Sequences Used in Real-Time PCR

Agarose Gel Electrophoresis Analysis

To further confirm the presence of distinct splice variants, PCR products of Kcnq2, Fgf1 and Kcnc3 were separated by 2% agarose gel electrophoresis in 1× TAE buffer. Prior to electrophoresis, PCR products were mixed with 6× loading buffer (Takara) to a final concentration of 1×. A DL1000 DNA marker (Takara) was used as the molecular size standard. Electrophoresis was performed at 140 V for 20 min, and gels were visualized under UV illumination using a Bio-Rad imaging system. The presence of multiple bands was used to confirm distinct splice variants of the target genes.

Quantification and Statistical Analysis

All data are presented as mean ± standard error of the mean (SEM). qPCR data were analyzed using a two-tailed unpaired Student’s t-test in GraphPad Prism v9.0. A P value < 0.05 was considered statistically significant.

TMT-Based Quantitative Proteomic Analysis of TG Tissues Following pIONT

TG samples from the pIONT10d and sham10d groups were subjected to TMT-based quantitative proteomic analysis to identify protein-level alterations associated with trigeminal neuropathic pain. Three pooled biological replicates were analyzed per group, and each replicate consisted of two TGs.

Proteins were extracted from TG tissues using an SDS-containing lysis buffer, and protein concentrations were determined using a BCA assay. Equal amounts of protein from each sample were reduced with dithiothreitol, alkylated with iodoacetamide, and digested with trypsin. The resulting peptide mixtures were desalted, dried, reconstituted in TEAB buffer, and labeled with tandem mass tag (TMT) reagents according to the manufacturer’s protocol. After labeling, peptide samples were combined, desalted, and dried before chromatographic separation.

The labeled peptides were fractionated by high-pH reversed-phase HPLC and subsequently analyzed by liquid chromatography–tandem mass spectrometry (LC-MS/MS). Raw mass spectrometry data were processed using Proteome Discoverer software (v2.4.1.15), and peptide identification was performed with the SEQUEST search engine against the Mus musculus UniProt database. Trypsin was specified as the digestion enzyme with up to two missed cleavages. Methionine oxidation was set as a variable modification, whereas carbamidomethylation of cysteine and TMT labeling of peptide N-termini and lysine residues were set as fixed modifications. Peptide identifications were filtered at a false discovery rate of 1%, and protein abundance was calculated based on normalized unique-peptide intensities.

DEPs were defined using a fold-change cutoff of ≥1.2 or ≤0.83 together with P value < 0.05. P values for proteomic comparisons were calculated using a two-tailed unpaired Student’s t-test, and the DEP screening was used for exploratory integration with transcriptomic AS data.

Results

Dynamic Changes of AS Events During pIONT-Induced TNP

To investigate dynamic changes in AS during pIONT-induced TNP, TG samples were collected at days 3 and 10 after surgery and subjected to poly(A)-enriched RNA sequencing. AS events were first summarized using ASprofile-based classification at both the individual-sample level and in the merged dataset, including transcription start site-associated events (TSS), transcription termination site-associated events (TTS), alternative exon events (AE), exon skipping events (SKIP), and intron retention events (IR) (Figure 1B). This analysis showed that SKIP-, TSS-, and AE-related events represented the major AS categories in the merged dataset, whereas IR and TTS events accounted for relatively smaller proportions, indicating broad transcript-structure diversity in TG samples during TNP progression.

Differential AS events between pIONT and sham groups were then identified using rMATS. The five rMATS-defined AS event types are schematically illustrated in Figure 1A, and differential AS events between pIONT and sham groups were summarized in Figure 1C. Among these, SE was the predominant differential AS pattern at both time points, accounting for approximately 67% and 51% of differential AS events at days 3 and 10, respectively. A progressive increase in differential AS events was observed over the course of TNP progression (Figure 1C). On day 3, 176 differential SE events were identified, and this number increased to 588 on day 10. Similar increases were also observed for other AS event types, including RI (11 vs. 44), MXE (16 vs. 388), A5SS (17 vs. 58), and A3SS (44 vs. 72) at days 3 and 10, respectively. These findings indicate stage-dependent changes in differential AS events between day 3 and day 10, two time points that were previously associated with distinct pain-related behavioral states in the pIONT model.¹⁸

To validate representative AS events, Kcnq2, Fgf1 and Kcnc3 were selected based on their known roles in neuronal excitability, synaptic plasticity, and pain modulation. Two primer sets were designed for each gene to target isoforms with or without the skipped exon (Figure 1D–1F). qPCR analysis confirmed significant differences in splice variant expression between the sham and pIONT groups (Figure 1G–1I). Agarose gel electrophoresis further showed PCR products with the expected sizes for the corresponding splice variants of Kcnq2, Fgf1 and Kcnc3 (Figure 1J–1L), confirming the successful amplification of distinct transcript isoforms.

Genomic distribution analysis showed enrichment of SE events on chromosomes 4 and 7 (Figure 2). At day 3, 66 and 44 SE events were detected on chromosomes 7 and 4, respectively; these numbers increased to 191 and 184 by day 10. Together, these results indicate a marked increase in differential AS events, particularly SE events, during the maintenance phase of TNP.

Figure 2 Chromosomal distribution of alternative splicing (AS) events in the mouse genome during trigeminal neuropathic pain (TNP). Circos plots showing the genomic distribution of AS events across mouse chromosomes at two stages after partial infraorbital nerve transection (pIONT). (A) Alternative 3′ splice site (A3SS) events at day 3 post-pIONT. (B) Alternative 5′ splice site (A5SS) events at day 3 post-pIONT. (C) Mutually exclusive exon (MXE) events at day 3 post-pIONT. (D) Retained intron (RI) events at day 3 post-pIONT. (E) Skipped exon (SE) events at day 3 post-pIONT. (F) A3SS events at day 10 post-pIONT. (G) A5SS events at day 10 post-pIONT. (H) MXE events at day 10 post-pIONT. (I) RI events at day 10 post-pIONT. (J) SE events at day 10 post-pIONT. Among the detected events, skipped exon (SE) events showed notable enrichment on chromosomes 4 and 7.

Functional Analysis of Alternatively Spliced Genes

To further investigate the functional significance of ASGs in pIONT-induced TNP, we analyzed their stage-specific expression abundance patterns and biological annotations. Hierarchical clustering of ASG expression abundance showed more prominent changes on day 10 than on day 3, suggesting stronger transcript-level alterations among ASGs during the maintenance phase (Figure 3A and B). Representative pain-related genes also exhibited temporal specificity. For example, Trpa1 was increased at day 3, whereas Trpv1 was decreased at day 10. These findings suggest stage-dependent expression changes among pain-related ASGs.²⁴

Figure 3 Functional characterization of alternatively spliced genes (ASGs) during trigeminal neuropathic pain (TNP) progression. (A) Hierarchical clustering heatmap showing transcript abundance patterns of ASGs across experimental groups at day 3 after pIONT. (B) Hierarchical clustering heatmap showing transcript abundance patterns of ASGs across experimental groups at day 10 after pIONT. Red and blue indicate relatively higher and lower abundance levels, respectively. Samples are ordered according to hierarchical clustering. (C and D) Gene Ontology (GO) enrichment analysis of ASGs at day 3 (C) and day 10 (D), highlighting enrichment of cellular structure- and polarity-related terms at day 3 and synaptic remodeling- and neurotransmission-related terms at day 10. (E) KEGG pathway enrichment analysis of day 10 ASGs, showing significant enrichment in pain-related pathways, including MAPK signaling and Ras signaling. For GO and KEGG enrichment analyses, P values were adjusted using the Benjamini–Hochberg method, and enrichment scores are shown as −log10-adjusted P values.

GO enrichment analysis revealed distinct functional profiles at the two time points. On day 3, ASGs were primarily associated with biological processes related to cellular structure and polarity (Figure 3C), including establishment or maintenance of cell polarity (GO:0007163), actin filament organization (GO:0007015), and cell junction assembly (GO:0034329), suggesting early structural adaptation after injury. In contrast, ASGs on day 10 were enriched in synaptic and neurotransmission-related functions (Figure 3D), including synapse organization (GO:0050808) and postsynapse organization (GO:0099173), indicating a shift toward synaptic remodeling during the maintenance phase.

KEGG pathway analysis showed that no significantly enriched pathways were detected for day 3 ASGs under the preset significance criteria. In contrast, day 10 ASGs were significantly enriched in several pain-related pathways (Figure 3E), including the MAPK signaling pathway (adjusted P = 0.03710), Ras signaling pathway (adjusted P = 0.03842), pathways of neurodegeneration - multiple diseases (adjusted P = 0.03842), and axon guidance (adjusted P = 0.03710). Among these, the MAPK pathway, particularly the ERK branch, has been widely implicated in synaptic plasticity and neuroinflammatory signaling in pain pathophysiology.²⁴

Collectively, these findings indicate that ASGs exhibit distinct temporal functional profiles during pIONT-induced TNP, with day 3 changes being more closely associated with structural adaptation and day 10 changes being more strongly linked to synaptic remodeling and pain-related signaling pathways.

Functional Overlap Between Alternative Splicing and Differential Gene Expression

To investigate stage-specific relationships between alternative splicing and transcriptional regulation during TNP progression, overlaps between differentially expressed alternatively spliced genes (DEASGs) and DEGs were analyzed. A total of 19 overlapping genes were identified on day 3, whereas this number increased markedly to 127 on day 10 (Figure 4A). Heatmap analysis showed that overlapping genes at day 3 exhibited relatively moderate expression changes (Figure 4B), whereas overlapping genes at day 10 showed more pronounced expression changes (Figure 4C), suggesting a stronger association between AS-related changes and transcriptional alterations during the maintenance phase of TNP.

Figure 4 Functional overlap between alternative splicing and transcriptional regulation during trigeminal neuropathic pain (TNP). (A) Venn diagrams showing the overlap between differentially expressed alternatively spliced genes (DEASGs) and differentially expressed genes (DEGs) at day 3 and day 10 after partial infraorbital nerve transection (pIONT). (B and C) Heatmaps showing the expression profiles of overlapping genes at day 3 (B) and day 10 (C). Samples are ordered according to hierarchical clustering. (D and E) Gene Ontology (GO) enrichment analysis of overlapping genes at day 3 (D) and day 10 (E). (F) KEGG pathway enrichment analysis of overlapping genes at day 10, showing significant enrichment in neural activity- and synaptic plasticity-related pathways, including circadian entrainment, long-term potentiation, cholinergic synapse, and glutamatergic synapse, as well as several hormone-associated signaling pathways. For GO and KEGG enrichment analyses, P values were adjusted using the Benjamini–Hochberg method, and enrichment scores are shown as −log10-adjusted P values.

To further examine the functional relevance of these overlapping genes, GO enrichment and KEGG pathway analyses were performed. On day 3, overlapping genes were enriched in biological processes related to cellular morphology, synaptic signaling, and intercellular communication (Figure 4D), including positive regulation of homotypic cell-cell adhesion (GO:0034112), regulation of integrin-mediated signaling pathway (GO:2001044), regulation of receptor-mediated endocytosis (GO:0048259), and signal release from synapse (GO:0099643). These findings suggest that early overlap between AS and differential expression is associated with neuronal adaptation and intercellular communication after injury. By day 10, overlapping genes were significantly enriched in biological processes associated with neuronal survival, axonal remodeling, basement membrane organization, and cell junction regulation (Figure 4E), including neuron death (GO:0070997), regulation of neuron death (GO:1901214), basement membrane organization (GO:0071711), axonogenesis (GO:0007409), negative regulation of neuron death (GO:1901215), basement membrane assembly (GO:0070831), and regulation of cell junction assembly (GO:1901888).

KEGG pathway analysis further showed that day 10 overlapping genes were significantly enriched in neural activity-, synaptic plasticity-, and hormone-associated signaling pathways, including circadian entrainment (adjusted P = 0.00667), long-term potentiation (adjusted P = 0.00667), cholinergic synapse (adjusted P = 0.01720), and glutamatergic synapse (adjusted P = 0.01720) (Figure 4F). No significantly enriched KEGG pathways were detected for day 3 overlapping genes under the preset significance criteria.

Collectively, these results indicate that the overlap between DEASGs and DEGs exhibits stage-specific functional characteristics during TNP progression, with day 3 changes being more closely associated with early neuronal adaptation and intercellular communication, whereas day 10 changes are more strongly linked to neuronal remodeling, extracellular matrix organization, synaptic plasticity-related pathways, and hormone-associated signaling.

Cellular Localization and Intercellular Communication of ASGs

To characterize the cellular context of ASGs, we mapped genes harboring AS events onto published single-cell RNA-seq-defined TG cell populations. TG cells were classified into eight neuronal and seven non-neuronal subtypes based on marker gene expression profiles (Supplementary Figure 1A and B).²⁵ This analysis was used to infer the cell-type distribution of ASGs at the gene-expression level, rather than to directly quantify cell-type-specific splice isoforms. Many ASGs showed marked cell-type-enriched expression patterns across TG cell populations (Supplementary Figure 1C–S).

Among neuronal subtypes, several ASGs were associated with excitability, nociceptive signaling, or neuronal remodeling. Nav2 (Supplementary Figure 1D), a neuron navigator gene associated with neuronal development and axonal processes, showed cell-type-enriched expression patterns.²⁶ Kcnq2 (Supplementary Figure 1E), a regulator of neuronal membrane excitability, was detected in neuronal populations and may be relevant to altered excitability after nerve injury.²⁷ Trpv1 and Trpa1 (Supplementary Figure 1F and G), which mediate thermal and chemical nociception, respectively, showed distinct cellular expression patterns, suggesting potential involvement of nociceptive neuronal subsets in ASG-associated changes.^28,29 Ccl25 (Supplementary Figure 1J) was detected in immune-responsive neuronal populations, whereas Hgf (Supplementary Figure 1K), a neurotrophic regulator, may be involved in nerve repair and pain-related signaling.^30,31 In addition, Shc4 (Supplementary Figure 1M), an adaptor protein linked to MAPK signaling, and Gal (Supplementary Figure 1N), a neuropeptide-related gene, were detected in neuronal subsets, suggesting possible involvement in synaptic plasticity and pain maintenance.³²

Several ASGs also showed selective expression in non-neuronal subtypes, suggesting possible roles in neuroimmune modulation and glial support. For example, P2ry14 (Supplementary Figure 1R), which is linked to inflammatory signaling, was enriched in immune-like glial populations.³³ Some ASGs, such as Adgrl2 (Supplementary Figure 1H), were expressed across both neuronal and non-neuronal populations, suggesting broader relevance to cell adhesion and intercellular signaling.³⁴ Genes related to calcium signaling, including Cacnb4 (Supplementary Figure 1C), Cacna1g (Supplementary Figure 1I), Cacna2d3 (Supplementary Figure 1L), and Camk1 (Supplementary Figure 1O), were also detected, indicating potential involvement in synaptic regulation. In addition, Pde2a (Supplementary Figure 1P), which is involved in cAMP/cGMP signaling,³⁵ Gde1 (Supplementary Figure 1Q), which is related to phospholipid metabolism, and Adgrg6 (Supplementary Figure 1S), which is associated with myelination and axonal support, may contribute to signal transduction and maintenance of TG homeostasis.³⁶

To further investigate potential intercellular communication associated with ASGs, ligand–receptor interaction networks were constructed for days 3 and 10 post-pIONT using CellLinker (Figure 5A and B). Receptor partners were computationally predicted, revealing complex many-to-many ligand–receptor associations. On day 3, predicted interactions such as Hgf–Met and Slit2–Robo2 involved ligand–receptor families previously implicated in nerve repair, axon guidance, neural remodeling, and synaptic connectivity, suggesting possible involvement in early injury-associated cellular communication.^37,38 By day 10, the predicted interaction network became more complex, with interactions involving Fgf–Fgfr2, Slit2–Robo, and Col1a2–Ddr1. These ligand–receptor families have been previously implicated in post-injury synapse formation, axon guidance, synaptic connectivity, and collagen-associated extracellular matrix remodeling, suggesting their possible involvement in persistent pain-related intercellular communication.^38–40

Figure 5 ASG-associated ligand–receptor interaction networks during trigeminal neuropathic pain (TNP) progression. (A and B) Predicted ligand–receptor interaction networks associated with representative alternatively spliced genes (ASGs) at day 3 (A) and day 10 (B) after partial infraorbital nerve transection (pIONT), generated using CellLinker. Edge width indicates predicted interaction strength. (C–F) Representative expression patterns of selected ligand–receptor-related genes mapped onto trigeminal ganglion (TG) cell clusters: Eng and Zyx (C), Gpr151 and Gal (D), Tspan5 and Adam10 (E), and Dag1 and Agrn (F). Color intensity indicates relative expression levels across cell populations.

Representative expression patterns of selected ligand–receptor-related genes are shown in Figure 5C–F. Eng and Zyx (Figure 5C) were detected in vascular-associated and neural cell populations, suggesting potential relevance to vascular–neural communication. Gpr151 and Gal (Figure 5D) were expressed across multiple cell populations; given the known involvement of GPR151 in neuropathic pain and the broader role of neuron–glia signaling in TG pain states, these patterns may reflect cellular contexts relevant to pain-associated intercellular communication.⁴¹ Tspan5 and Adam10 (Figure 5E) showed broad expression patterns, consistent with the known role of Tspan5 in ADAM10 regulation and the involvement of ADAM10 in cell-surface protein processing and cell adhesion.⁴² Dag1 and Agrn (Figure 5F) may be associated with extracellular matrix-linked cell communication and synaptic organization.⁴³

Collectively, these findings suggest that genes harboring AS events show distinct gene-level expression and localization patterns across diverse TG cell types. Together with predicted ligand–receptor interactions, these results provide cellular-context support for a possible association between ASG-related changes and intercellular communication, calcium signaling, neuroimmune regulation, and extracellular matrix remodeling during TNP progression.

Identification of RBP Genes Harboring Differential AS Events During TNP

To investigate post-transcriptional regulators potentially involved in AS during TNP, RNA-seq data were compared with a curated mouse RBP database (https://rbp2go.dkfz.de/). RBPs are key regulators of post-transcriptional processes and play central roles in splicing regulation. A total of 222 and 912 DEASGs were identified at days 3 and 10 post-pIONT, respectively. Intersecting these genes with 2,887 annotated mouse RBP genes identified 44 and 229 RBP genes harboring differential AS events at the two respective time points (Figure 6A and B), indicating a marked expansion of putative splicing regulators during the maintenance phase.

Figure 6 Identification of RNA-binding protein (RBP)-associated alternative splicing changes during trigeminal neuropathic pain (TNP). (A and B) Venn diagrams showing the overlap between differentially expressed alternatively spliced genes (DEASGs) and annotated mouse RNA-binding protein (RBP) genes at day 3 (A) and day 10 (B) after partial infraorbital nerve transection (pIONT). (C) Heatmap showing transcript abundance patterns of RBP genes identified at both time points. (D) Heatmap showing the expression profiles of selected pain-related or RBP-associated genes. Color scale indicates normalized transcript abundance/expression level. (E) STRING-based interaction network of selected RBP genes harboring differential AS events and associated molecules. Edge thickness indicates interaction confidence.

Overlap analysis revealed five RBP genes that harbored differential AS events at both time points: Rbm34, Srrt, Atxn2, Rbfox1, and Sltm. Among them, Atxn2 was consistently downregulated throughout the course of TNP (Figure 6C), whereas Rbm34 and Rbfox1 showed dynamic expression changes across stages, suggesting possible involvement in stage-dependent regulation of neuron-related splicing events.

Several RBP-associated molecules in the regulatory analysis were also linked to established neuropathic pain mechanisms. In particular, SRSF1 has previously been implicated in pathological pain-related splicing regulation.^9,44 Additional network-associated genes, including Atf4, Asic1, Nrg1, and Trpv1, also showed differential expression at day 10 post-pIONT (Figure 6D), supporting links between the RBP-associated network and pathways related to neuronal excitability, synaptic transmission, and chronic pain progression.

To further delineate the regulatory architecture, an interaction network involving RBP genes harboring differential AS events and associated molecules was constructed (Figure 6E). The resulting network showed extensive connectivity, with Trpv1 and Bdnf connected to multiple network-associated molecules. Members of the ASIC family, including Asic1 and Asic2, occupied central positions in the network, while Fus and Atxn2 displayed distinct topological features, suggesting that they may act as important regulatory hubs in the AS-associated network.

Collectively, these findings support the involvement of RBPs in the modulation of AS during TNP and identify candidate post-transcriptional regulators associated with pain-related molecular networks.

Identification of Differentially Expressed Proteins During TNP

To evaluate protein-level alterations associated with TNP, TMT-based quantitative proteomic analysis was performed on TG tissues collected at day 10 post-pIONT. By integrating the proteomic data with transcriptome-based splicing analysis, we identified 45 DEPs that overlapped with protein products encoded by genes harboring differential AS events (Figure 7A).

Figure 7 TMT-based proteomic profiling reveals functional convergence with alternative splicing in trigeminal neuropathic pain. (A) Venn diagram showing the overlap between differentially expressed proteins (DEPs) identified by TMT-based quantitative proteomic analysis and protein products encoded by genes harboring differential AS events at day 10 after partial infraorbital nerve transection (pIONT). (B) Gene Ontology (GO) enrichment analysis of the overlapping proteins, highlighting biological processes related to synaptic vesicle endocytosis, synaptic vesicle recycling, and regulation of postsynaptic neurotransmitter receptor levels. Representative protein abundance profiles of overlapping proteins are provided in Supplementary Figure 2. For GO enrichment analyses, P values were adjusted using the Benjamini–Hochberg method, and enrichment scores are shown as −log10-adjusted P values.

GO enrichment analysis indicated that these overlapping DEPs were significantly associated with biological processes related to synaptic function (Figure 7B), including synaptic vesicle endocytosis (GO:0048488), presynaptic endocytosis (GO:0140238), synaptic vesicle recycling (GO:0036465), positive regulation of synaptic vesicle recycling (GO:1903423), regulation of postsynaptic neurotransmitter receptor levels (GO:0099072), and regulation of synaptic vesicle cycling (GO:0099504). These findings suggest that protein-level alterations at day 10 are closely linked to processes required for sustained synaptic transmission and neuronal signaling.⁴⁵

The 45 overlapping proteins were associated with synaptic transmission, neuronal viability, intracellular signaling, and metabolic regulation. Representative examples included GAL (Supplementary Figure 2A), C8G (Supplementary Figure 2B), and MED17 (Supplementary Figure 2C), which are related to pain modulation, neuroinflammatory responses, and transcriptional regulation, respectively.^46–49 Additional overlapping proteins, including PRNP, DLG3, FGF12, and DNM1 (Supplementary Figure 2D, E, H and I), are associated with synaptic plasticity, glutamatergic signaling, voltage-gated ion channel regulation, and synaptic vesicle cycling.^50–53 Other candidates, such as FGF1, OSF1, CAMK1, SHC4, and GDE1 (Supplementary Figure 2F, G, J–L), are linked to nerve repair, intracellular signaling, calcium regulation, and phospholipid metabolism.^32,54–57

Collectively, these results show that protein-level alterations during TNP converge with AS-related transcriptomic changes in pathways associated with synaptic function, neuroinflammation, and neuronal maintenance, supporting a potential multilayered association between AS-related transcriptomic alterations and protein-level changes during persistent pain.

A concise summary of the key validated ASGs, RBPs, and overlapping DEPs identified across stages is provided in Table 2.

Table 2 Condensed Summary of Key Validated ASGs, Candidate RBPs, and Overlapping DEPs in pIONT-Induced TNP

Discussion

This study provides a temporal analysis of AS dynamics in the TG during the onset and maintenance phases of TNP. AS is a major post-transcriptional mechanism that expands transcriptomic complexity by generating multiple mRNA isoforms from single genes.⁵⁸ In the present study, exon skipping (SE) was the predominant AS event at both time points and showed a marked increase in complexity at day 10 after pIONT. These findings suggest that stage-specific splicing programs are associated with TNP progression and are consistent with temporal changes in neuronal signaling during pain development.

An additional notable finding was the enrichment of SE events on chromosomes 4 and 7. Because these chromosomes contain pain-related genes, including Asic1,⁵⁹ this non-random distribution raises the possibility that chromosomal context may influence AS regulation after nerve injury. Although the mechanisms remain to be clarified, this pattern is consistent with the idea that nerve injury triggers coordinated transcriptomic remodeling rather than isolated gene-level changes.^60,61

At the molecular level, several AS-related changes involved genes associated with neuronal excitability and synaptic regulation, including Kcnq2, Kcnc3, and Fgf1, which were validated by qPCR and agarose gel electrophoresis. In addition, enrichment analyses highlighted stage-dependent functional differences. Day 3 ASGs were more closely associated with structural adaptation, whereas day 10 ASGs were enriched in synaptic and signaling-related pathways, including MAPK signaling. These findings are consistent with a transition from early injury-related cellular remodeling to later changes associated with synaptic reorganization and persistent pain-related signaling.^62,63

Our results also extend previous observations of AS dysregulation in neuropathic and inflammatory pain models. Prior studies in spinal and peripheral nerve injury models have identified stage-specific splicing changes and implicated splicing regulators such as SRSF1 in pain-related plasticity.⁴⁴ In this context, the present study adds a TG-focused, stage-resolved, and multi-omic perspective to a comparatively underexplored area. The integration of AS profiling with TMT-based quantitative proteomic analysis and single-cell RNA-seq-based cell-type mapping further suggests that AS-related changes in TNP are linked not only to transcript-level variation but also to protein-level alterations and cell-type-specific localization patterns.

The overlap analyses further support this interpretation. During the onset phase, AS-related changes were mainly associated with processes such as cytoskeletal remodeling, cell junction organization, and neurotransmitter-related communication, whereas during the maintenance phase they were more strongly linked to synaptic remodeling, vesicle cycling, and neural circuit-related pathways. At the protein level, overlapping candidates such as DNM1, DLG3, FGF1, FGF12, and CAMK1 point to potential convergence between splicing-related regulation and pathways involved in synaptic function, intracellular signaling, and neuronal maintenance. These findings do not establish direct causal effects of individual splice variants, but they support the view that AS is associated with broader molecular reorganization during TNP progression.

Cell-type mapping and ligand–receptor analyses provide additional context for these stage-dependent changes. Several ASGs were preferentially detected in nociceptive neuronal subtypes, including genes related to ion channel function and excitability, whereas others were enriched in non-neuronal populations associated with immune or glial support. These observations suggest that AS in TNP is not restricted to neurons, but may also involve cell-type-specific contributions from non-neuronal compartments. Similarly, predicted ligand–receptor interactions, including Hgf–Met, Slit2–Robo, and Fgfr2–Fgf, support the possibility that AS-related molecular changes are linked to altered intercellular communication during pain maintenance.^37,64 In parallel, RBP genes harboring differential AS events, including Atxn2, Rbfox1, and Rbm34, were identified as candidate post-transcriptional regulators associated with these AS programs. Together, these results support a model in which AS-related regulation, cell-type specificity, and intercellular signaling are linked within a broader network of TNP-associated molecular adaptation.^65,66 Network analysis further suggests that these RBPs may function as regulatory hubs linking extracellular cues with coordinated transcriptomic reprogramming, ultimately shaping the functional architecture of pain circuits. These findings parallel recent studies highlighting the role of splicing factors in pain pathogenesis and support further investigation of their therapeutic relevance.^9,44

Despite these findings, several limitations should be acknowledged. First, both transcriptomic and proteomic analyses were performed using pooled trigeminal ganglion (TG) samples, and the number of pooled biological replicates was limited for each condition. Although this design was adopted to obtain sufficient material from small TG tissues for high-throughput analyses, it reduced the ability to assess mouse-to-mouse biological variability and limited statistical power. The use of ICR mice matched for strain, age, sex, source, and experimental condition was intended to reduce non-experimental variability among animals, and inter-individual variability was expected to be relatively limited under these controlled conditions. Nevertheless, because individual TG tissues were pooled before sequencing, the present study was not designed to characterize AS variability among individual animals. Accordingly, the primary aim of this study was to identify AS changes associated with the pIONT pain-model condition compared with sham-operated control mice without nerve transection, rather than to evaluate AS variability attributable to inter-individual differences. Therefore, the present findings should be interpreted as exploratory and descriptive rather than as definitive evidence of causality. Second, many of the regulatory relationships described here were derived from predictive bioinformatics analyses and therefore require independent experimental validation. Third, although the pIONT model reproduces important features of trigeminal neuropathic pain, differences between this mouse model and human TNP should be considered when evaluating the generalizability of the present findings.

In summary, our data support the presence of dynamic, stage-specific AS programs in the TG during TNP progression, with stronger AS complexity and broader multi-omic convergence during the maintenance phase. These findings highlight AS, RBP-associated regulation, and cell-type-specific signaling as potentially important components of TNP pathophysiology. From a translational perspective, these observations raise the possibility that AS-related mechanisms may represent future therapeutic targets. However, this implication remains speculative at present and will require validation in independent cohorts, mechanistic experiments, and studies designed to determine whether specific splice variants or their regulators can be therapeutically modulated in trigeminal pain states.

Conclusion

This study demonstrates dynamic and stage-specific alterations in AS in the trigeminal ganglion during TNP. Exon skipping was the predominant AS event during both the onset and maintenance phases, with greater AS complexity observed at the maintenance stage. Integrative analyses further suggest that AS-related changes are associated with cell-type-specific expression patterns and with multilayered regulation involving RNA-binding proteins and protein-level alterations. Together, these findings provide a descriptive and multi-omic framework for understanding AS-related molecular changes in TNP and support further investigation of AS-associated mechanisms in trigeminal pain states.