Orthodontic force induces nerve injury-like transcriptomic changes driven by TRPV1-expressing afferents in mouse trigeminal ganglia.

Orthodontic force produces mechanical irritation and localized inflammation in the periodontium, which causes pain in most patients. Nocifensive behaviors resulting from orthodontic force in mice can be substantially attenuated by intraganglionic injection of resiniferatoxin (RTX), a neurotoxin that specifically ablates a subset of neurons expressing transient receptor potential vanilloid 1 (TRPV1). In the current study, we determined changes in the transcriptomic profiles in the trigeminal ganglia (TG) following the application of orthodontic force, and assessed the roles of TRPV1-expressing afferents in these transcriptomic changes. RTX or vehicle was injected into the TG of mice a week before the placement of an orthodontic spring exerting 10 g of force. After 2 days, the TG were collected for RNA sequencing. The application of orthodontic force resulted in 1279 differentially expressed genes (DEGs) in the TG. Gene ontology analysis showed downregulation of gliogenesis and ion channel activities, especially of voltage-gated potassium channels. DEGs produced by orthodontic force correlated more strongly with DEGs resulting from nerve injury than from inflammation. Orthodontic force resulted in the differential expression of multiple genes involved in pain regulation, including upregulation of Atf3, Adcyap1, Bdnf, and Csf1, and downregulation of Scn10a, Kcna2, Kcnj10, and P2ry1. Orthodontic force-induced DEGs correlated with DEGs specific to multiple neuronal and non-neuronal subtypes following nerve injury. These transcriptomic changes were abolished in the mice that received the RTX injection. These results suggest that orthodontic force produces transcriptomic changes resembling nerve injury in the TG and that nociceptive inputs through TRPV1-expressing afferents leads to subsequent changes in gene expression not only in TRPV1-positive neurons, but also in TRPV1-negative neurons and non-neuronal cells throughout the ganglia. Orthodontic force-induced transcriptomic changes might be an active regenerative program of trigeminal ganglia in response to axonal injury following orthodontic force.

Introduction

Malocclusion of teeth is treated orthodontically by realigning the malpositioned teeth using physical force exerted by various orthodontic appliances. Orthodontic tooth movement produces pain and soreness in almost all patients with fixed orthodontic appliances.[1] However, methods for alleviating orthodontic pain are not well developed and a better understanding of the mechanisms involved should facilitate the development of novel treatments for orthodontic pain.[2]

Orthodontic force produces compression and tension of the periodontal ligament around teeth and generates mechanical irritation to the primary afferent terminals, which is projected to the periodontal ligament.[3,4] Such compression of the periodontal ligament leads to localized inflammatory responses, which is highly likely to produce peripheral sensitization in the periodontal nociceptors.[5] Consequently, orthodontic pain involves spontaneous pain and pain around the teeth upon mastication.[1] Our group and others have developed methods for reproducing orthodontic pain by applying an experimental orthodontic force in rodents.[3,6-8] By placing a well-calibrated coil spring between the maxillary first molar and incisors in mice, precise mechanical force can be applied to the localized periodontium around the affected teeth. We assessed pain-like behaviors such as changes in facial grimace or bite force in mice,[3] which mimic spontaneous pain and bite-evoked pain in patients. Such pain-like behaviors peak at 1 day and resolve after approximately 1 week in rodents,[3] which is consistent with the time course in patients fitted with orthodontic appliances.[1] This mouse model is therefore a clinically relevant pain model that produces reliable pain-like behaviors and allows us to study the mechanisms of tooth-related pain. Using the experimental tooth movement model, TRPV1 and TRPV1-expressing (TRPV1+) afferents were shown to mediate orthodontic pain-like behaviors.[3,9,10] Several other molecules in the primary afferents, such as calcitonin gene-related peptides (CGRP), acid-sensing ion channel 3 (ASIC3), or P2X3, are also implicated in orthodontic pain.[11-13] Orthodontic force leads to changes in the expression of these molecules, suggesting that the neural plasticity of the primary afferents contributes to orthodontic pain.

Peripheral injury induces widespread changes within the sensory ganglia. Such intraganglionic changes could lead to peripheral sensitization, primary hyperalgesia, and extraterritorial ectopic pain.[14] In response to peripheral injury or inflammation, infiltration of inflammatory cells, proliferation of satellite glia, and release of cytokines and chemokines occur in the trigeminal ganglia (TG), which eventually produces alterations in gene expression in neuronal and non-neuronal cells. For example, approximately 17% of all genes are differentially regulated in the TG following inflammation of the masseter muscle.[15] Transcriptomic analysis has shown that masseter inflammation upregulates a group of pronociceptive genes, allowing us to infer events that are relevant to the hyperalgesia associated with craniofacial muscle inflammation. Injury of the sciatic nerve also produces differential regulation of approximately 6.5% of all genes in the dorsal root ganglia (DRG) within 3 days.[16] Nerve injury differentially regulates the expression of genes related to chronic pain development and axonal regeneration. Since orthodontic force produces persistent mechanical irritation to the afferent terminals within the periodontal ligaments accompanied by localized inflammation, orthodontic force could lead to transcriptomic changes similar to those observed in response to both peripheral inflammation and nerve injury.

Although it is not clearly understood how peripheral injury or inflammation leads to changes in widespread gene expression within sensory ganglia, it is highly likely that nociceptive signals transmitted from the periphery to the sensory ganglia through nociceptive afferents drive the changes. Since TRPV1+ afferents contribute to orthodontic force-induced nocifensive behaviors,[3] persistent activation of TRPV1+ nociceptors could, in principle, deliver neural inputs and drive changes in gene expression within the TG. However, the contribution of TRPV1+ afferents to transcriptomic changes in the sensory ganglia following peripheral injury has not been determined.

The objectives of this study were to investigate orthodontic force-induced changes in transcriptomic profiles within the TG and to determine the role of TRPV1+ afferents in such changes. To do so, we performed transcriptomic analysis in the TG following the application of orthodontic force in mice. To determine the contribution of nociceptors to transcriptomic changes, TRPV1+ trigeminal afferents were chemically ablated. We also compared the transcriptomic profiles and differentially expressed genes following the application of orthodontic force with datasets obtained from four previous studies employing masseter inflammation and nerve injury.[15-18]

Methods

Experimental animals

We used sixteen 8-week-old adult male C57BL/6J mice (Jackson Laboratory, Bar Harbor, ME, USA). All animals were housed in a temperature-controlled room under a 12:12 light–dark cycle with ad libitum access to food and water. All animal experimental studies and procedures were conducted in accordance with the NIH Guide for the Care and Use of Laboratory Animals (Publication No. 80-23) and under a University of Maryland-approved Institutional Animal Care and Use Committee protocol.

Microinjection into the TG

To determine the contribution of TRPV1+ afferents to orthodontic force-induced transcriptomic changes, resiniferatoxin (RTX) or vehicle was directly injected into the TG as described previously.[3] RTX is a highly efficacious TRPV1 agonist. The activation of TRPV1 by localized injection of RTX leads to ablation of nociceptor terminals or soma.[3,19,20] Direct injection of RTX into TG ablates approximately half of TRPV1+ neurons in ophthalmic/maxillary regions, but not in mandibular regions of ganglia, and substantially attenuates nocifensive behaviors induced by the ocular application of capsaicin or experimental orthodontic force.[3] The animals were anesthetized using ketamine (100–150 mg/kg) and xylazine (10–16 mg/kg) and then placed in a Kopf stereotaxic apparatus. A midline incision and an opening to the skull were made. A 0.5-μl Hamilton microsyringe was used for microinjection. The microsyringe needle was placed in the left TG region according to the stereotaxic coordinates of the mouse brain (0.2 mm posterior to bregma, 1.3 mm lateral to the midline, and 6.5 mm deep), targeting the ophthalmic/maxillary (V1/V2) region. RTX (50 ng/0.5 μl; Sigma-Aldrich, St. Louis, MO, USA) was dissolved in phosphate-buffered saline (PBS) containing 1% dimethyl sulfoxide and 10% Tween-80. Mice injected with vehicle (0.5 μl) served as a control group. Injection was performed at a rate of 0.5 μl/min and the injection needle was held in the tissue for 2 minutes before removal to allow diffusion.

Experimental orthodontic force

To apply the experimental orthodontic force in mice, a coil spring was placed between the maxillary first molar and maxillary incisors as described previously.[3] The animals were anesthetized with ketamine (100–150 mg/kg) and xylazine (10–16 mg/kg). A 0.010-in stainless steel ligature wire was looped around the first molar and a second ligature wire was looped around the maxillary incisors. We used two nickel–titanium orthodontic coil springs (Xu Jia Chuang Spring, Guangdong, China) exerting a 10 g force (wire diameter, 0.15 mm; outer diameter 1.8 mm; length, 2.2 mm) upon activation by 1 mm. A 10g force produces reliable nocifensive behaviors and tooth movement in mice.[3] In the orthodontic force group, the coil spring was extended mesially and ligated to the incisors. In the sham group, the orthodontic spring was irreversibly deformed by extension beyond the elastic limit and ligated so that the spring delivered no force. To secure the ligature wires, a self-etching primer and a light-cured adhesive resin cement (Transbond; 3M Unitek, Monrovia, CA, USA) were applied to the palatal surfaces of the maxillary incisors and first molars. After insertion of the spring, the animals were fed a soft diet (DietGel Recovery; ClearH2O; Portland, ME, USA) ad libitum. The appliances were inspected daily, and additional bonding material was applied when necessary.

RNA isolation

Two days following the placement of coil springs, the mice were anesthetized using ketamine and xylazine and the TG in the ipsilateral side were harvested from the RTX- or vehicle-injected mice. We chose day 2 since the maximum pain-like behavior was observed on day 1 through day 3 in our previous study,[3] and we presumed that changes of gene expression in TG could occur during days 1 to 3. We also assumed that gene expression changes on day 2 represented changes occurring during day 1 to 3 and might contribute to pain-like behaviors not only on day 2 but also on day 1 and 3. A total of sixteen dissected TG were stored in RNAlater (Invitrogen, Carlsbad, CA, USA). Total RNA was extracted using Trizol (Invitrogen, Carlsbad, CA, USA) and purified using an RNeasy kit (Qiagen, Germantown, MD, USA), which included a DNase treatment for removing genomic DNA. RNA integrity was evaluated by Agilent 2100 Bioanalyzer analysis (Agilent Technologies, Palo Alto, CA, USA). The RNA integrity number of all samples was greater than 8.5.

Library construction

A total amount of 3 µg of RNA per sample was used for library construction. Ribosomal RNA (rRNA) was removed using an Epicentre Ribo-zero rRNA Removal Kit (Epicentre, USA). Subsequently, strand-specific sequencing libraries were generated according to the dUTP method using the RNA produced using the NEBNext Ultra Directional RNA Library Prep Kit for Illumina (New England Biolabs, Ipswich, MA, USA). The library quality and quantity were assessed using the Agilent Bioanalyzer 2100 system and real-time PCR.

RNA sequencing and data analysis

RNA sequencing (RNA-seq) was performed using an Illumina Hiseq 2000 platform and 150 bp paired-end reads were generated. Reference genome and gene model annotation files were downloaded from NCBI, UCSC, and Ensembl. The index of the reference genome was built, and paired-end clean reads were aligned to the reference genome using STAR (v2.5).[21] The maximal mappable prefix method was used to generate a precise mapping result for junction reads.

HTSeq v0.6.1 was used for quantification of gene expression levels. The quantitative gene expression levels were calculated for each sample based on the number of fragments per kilobase of exon per million fragments mapped (FPKM), which is a normalized value for the length of the gene and the sequencing depth.[22] Differential expression analyses between two groups (four biological replicates per condition) were performed using the DESeq2 R package (2_1.6.3).[23] To cluster the samples based upon the similarity of their patterns of gene expression, we performed principal component analysis (PCA). Clustering of differentially expressed genes (DEGs) was performed using pheatmap. Enrichment analysis was performed using the ClusterProfiler R package (v2.4.3).[24] Gene ontology (GO) is a major bioinformatics classification system that is used to unify the presentation of gene properties across all species. This includes three main areas: cellular components, molecular functions, and biological processes. Kyoto Encyclopedia of Genes and Genomes (KEGG) is a database for pathway analysis. The Reactome is a database of reactions, pathways, and biological processes. In these analyses, the terms with q (adjusted p) < 0.05 were considered as significant enrichment.

To compare the DEGs from orthodontic force with other injury conditions, four datasets from previous studies were used: 1) 3499 DEGs from TG at 3 days following induction of inflammation of masseter muscle by injecting complete Freund’s adjuvant (CFA) in rats;[15] 2) 1615 DEGs from DRG at 3 days following sciatic nerve transection in mice;[16] 3) 2221 DEGs from TG at 2 days following infraorbital nerve transection determined by single-cell RNAseq;[18] 4) cell type-specific DEGs from DRG at 2 days following spinal nerve transection determined by single-cell RNAseq.[17] To functionally classify DEGs, we used ‘Molecular Function’ and ‘Protein Class’ modules in Gene List Analysis in the PANTHER Classification System (Version 15.0; www.pantherdb.org).[25] To infer the transcriptomic changes associated with pain, we compared the DEGs with the list of pain genes associated with acute and chronic pain. This list includes 681 genes from Pain Research Panels (Algynomics) and other sources.[26-28]

Real-time PCR assay

Real-time PCR was performed and analyzed as described previously.[15] Reverse transcription was carried out using the High Capacity cDNA Reverse Transcription Kit (Applied Biosystems, Foster City, CA). To generate cDNA, 500 ng of RNA was used as template. The sequences of primers were derived from previous studies[29-31] or were newly designed using Primer 3[32] (http://www.ncbi.nlm.nih.gov/tools/primer-blast/index.cgi). All of the primer pairs used in this study and average cycle numbers for each gene are described in Table 1. Dissociation curve analysis was performed to ensure single-product amplification for all primer pairs. Real time PCR was performed on the BioRad CFX384 Real Time System (BioRad, Hercules, CA) using assays specific to the genes of interest. Each reaction well contained 5 µl of PowerUp SYBR Green Master Mix (Applied Biosystems), cDNA equivalent to 10 ng of total RNA and 250 nM each of forward and reverse amplification primers in a final reaction volume of 10 µL. Cycling conditions were as follows: 95°C for 10 minutes for polymerase activation, followed by 40 cycles of 95°C for 15 seconds and 60°C for 1 minute, with a final incubation for a dissociation curve after cycling was complete. Data analysis was performed using CFX Manager software from BioRad, version 3.1. We calculated the ratios of the experimental Cq (cycle quantification) between the genes of interest and the endogenous control product glyceraldehyde-3-phosphate dehydrogenase (Gapdh), and used them to calculate the relative abundance of mRNA in each sample. Relative quantification of the mRNA was achieved using the comparative CT method (2−ΔΔ method).

Table 1.

Sequences of primers for real-time PCR analysis and average CT values.

Gene	Primer sequences (5’ to 3’)	Average CT
Adcyap1	TTAGCTTCTTCTCCCGGTGG	30.3
	CGCTACACATGGTCATTGGTG
Atf3	GTCACCAAGTCTGAGGCGG	27.5
	GTTTCGACACTTGGCAGCAG
Bdnf	TGTAGTCGCCAAGGTGGATG	32.0
	ACCTGGTGGAACATTGTGGC
Calca	TGACAGCATGGTTCTGGCTT	18.2
	GTCCCCAGAAGAGCAAGAGG
Csf1	TGCTAAGTGCTCTAGCCGAG	27.3
	CCCCCAACAGTCAGCAAGAC
Gapdh	GGACCTCATGGCCTACATGG	19.2
	TAGGGCCTCTCTTGCTCAGT
Kcnj10	CACCTTCGAGCCAAGATGAC	26.0
	GGCCACAGCTACCAGATACC
Kcnt2	CTTGAGAGCATGGGCTGTGA	29.1
	ACTGCTGCCCTTCTTGCC
Trpv1	ATCATCAACGAGGACCCAGG	23.6
	TGCTATGCCTATCTCGAGTGC

Statistical analysis

For differential gene expression analysis, the p-values were generated in DESeq and adjusted using the Benjamin–Hochberg method to control for false discovery.[33] The difference was considered to be significant when the q value, which is a p value with corrected false discovery ratio, was less than 0.05. For the real-time PCR assay, we conducted one-way ANOVA followed by Sidak post-hoc test for three pairs of comparisons as indicated in the figure legend. Data are presented as mean ± SEM and p < 0.05 was considered significant.

Results

Orthodontic force alters overall gene expression in TG

As illustrated diagrammatically in Figure 1(a), in adult C67BL/6 mice, we injected RTX or vehicle into TG. After 1 week, coil springs were placed between the maxillary first molar and the incisors with 10 g active orthodontic force. In the control group, coil springs without active force were placed (Sham). Consequently, we compared four groups of TG: vehicle-injected sham group (Veh/Sham; VS), vehicle-injected orthodontic force group (Veh/Force; VF), RTX-injected sham group (RTX/Sham; RS), and RTX-injected orthodontic force group (RTX/Force; RF). After 2 days, the TG were dissected out. To determine the gene expression profiles, RNA-seq assay was performed. We obtained on average 51 ± 15 million clean reads per sample with no significant difference between groups (p > 0.3 in one-way ANOVA). The averaged sequencing quality score Q30 was 94.3 ± 0.06%. Reads were mapped to the reference mouse genome and the total mapping rate was 95.4 ± 0.3%. Gene expression levels were estimated from the abundance of transcripts that mapped to the genome. The distribution of FPKM among the four groups was not different (not shown). To assess the inter- and intra-group variability among the four groups, Pearson correlation coefficients were calculated between all samples (Figure 1(b)), which showed a slight weak correlation between one sample (RS4) and the other samples. Since Pearson’s correlation analysis of RS4 showed r2 > 0.95 across all the other samples, we did not exclude this sample from further analysis. Principal component analysis (PCA) was also performed using the read counts for all genes in each sample normalized by library size (Figure 1(c)). The plot displays all 16 samples along PC1, PC2, and PC3, which describe 33.9%, 24.3%, and 6.85% of the variability, respectively. Overall, all four groups were distinct from one another in the three-dimensional representation, but only the segregation of the Veh/Force group from the other three groups was evident. The RTX/Sham group showed relatively high intragroup variation owing to the variation of RS4.

Figure 1.

Intra- and inter-group variability of transcriptomic data. (a) Time course of experiment. RTX, resiniferatoxin; Veh, vehicle; TG, trigeminal ganglia. (b) Pearson correlation coefficient (R2) matrix between all samples in four groups (n = 4 mice/group): Veh/Sham (VS), Veh/Force (VF), RTX/Sham (RS), and RTX/Force (RF). (c) Principal component analysis (PCA) plot of gene expression levels in four groups (VS, VF, RS, and RF).

Orthodontic force-induced gene expression is abolished in TG injected with RTX

In the differential gene expression analysis, pair-wise comparisons were performed between the groups of interest (Figure 2). Compared to the Veh/Sham, the vehicle-injected orthodontic force group (Veh/Force) exhibited 1279 DEGs (Figure 2(a) and (d); 636 upregulated and 643 downregulated). Among these, 1228 genes showed FPKM > 1 (609 upregulated and 619 downregulated). Compared to the Veh/Sham group, the RTX/Sham group showed 2083 DEGs (Figure 2(b) and (e); 1543 upregulated and 540 downregulated). Among these, 1537 genes showed FPKM > 1 (1016 upregulated and 521 downregulated), with a third of the upregulated genes detected at low frequency. DEGs following RTX injection were largely different from those DEGs following orthodontic force and only 181 genes overlapped (Figure 2(g); 9.5% of RTX-induced DEGs and 16.5% of orthodontic force-induced DEGs). Interestingly, the RTX/Force group showed only six DEGs compared to the RTX/Sham group (Figure 2(c) and (f); two upregulated and four downregulated; FPKM > 1 in all genes), suggesting that RTX injection into TG largely abolished orthodontic force-induced transcriptomic changes in the TG. It is possible that the paucity of DEGs in the RTX/Sham vs. RTX/Force comparison could be derived from large intra-group variation owing to one sample (RS4) showing large variation from the other samples. However, the impact of the sample appeared to not be dominant since differential gene analysis excluding RS4 showed only a modest increase in DEGs to 50 (not shown).

Figure 2.

Transcriptomic changes abolished by intra-TG RTX following orthodontic tooth movement. (a–c) Overall distribution of differentially expressed genes (DEGs) in volcano plots in pair-wise comparisons as indicated. The threshold of differential expression of genes is q (adjusted p value) < 0.05. (d–f) Top 10 DEGs showing upregulation (red) or downregulation (blue) in each comparison. *q (adjusted p) < 0.05, #q < 0.005, +q < 10−5, §q < 10−10, &q < 10−15. (g) The number of overlapping DEGs among three pair-wise comparisons.

Orthodontic force induces unique DEGs in TG

Functional enrichment analysis was performed to identify the biological relevance of the orthodontic force-induced transcriptomic changes (Figure 3). GO analysis of biological processes, cellular components, and molecular functions showed upregulation and downregulation of multiple terms (Figure 3(a) to (c)). Cell growth, sterol and cholesterol biosynthetic and metabolic processes, anion transmembrane transporter activity, and extracellular vesicle-related components were upregulated. In contrast, gliogenesis, ensheathment of axons and neurons, synapse organization, and synaptic membrane-related components were downregulated. Overall, ion channel activities, such as cationic or potassium channels, were also downregulated. Connective tissue development, collagen fibril organization, and extracellular matrix-related components were downregulated. In the KEGG pathway analysis, pathways for mineral absorption and transcriptional misregulation in cancer were upregulated, whereas pathways for antigen processing and presentation, phospholipase D signaling, and glycosaminoglycan biosynthesis were downregulated. In the Reactome enrichment assay, cholesterol biosynthesis was upregulated, whereas extracellular matrix organization was downregulated.

Figure 3.

Enrichment analysis of orthodontic force-induced DEGs. Gene ontology (GO) enrichment analysis in cellular components (a), molecular functions (b), and biological processes (c). (a) Regulation of cell morphogenesis involved in differentiation, (b) transcription factor activity, direct ligand-regulated sequence-specific DNA binding, (c) extracellular matrix structural constituent conferring tensile strength, and (d) glycosaminoglycan biosynthesis – heparan sulfate/heparin. (d) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. (e) Reactome enrichment analysis. *q (adjusted p) < 0.05, #q < 0.005, +q < 0.0005, §q < 10−4, &q < 10−5.

Orthodontic force induces nerve injury-like transcriptomic changes in TG

We performed a more detailed analysis of the differentially expressed potential pain-related genes in the Veh/Sham vs Veh/Force comparisons. To determine if orthodontic force-induced gene changes resemble the transcriptomic signatures of inflammation or nerve injury, we compared the orthodontic force-induced DEGs with masseter inflammation- and sciatic nerve injury-induced DEGs obtained from bulk RNAseq of entire ganglia. We compared the GO of the three groups of DEGs through functional classification using the PANTHER system (Figure 4(a)). The three groups of DEGs from orthodontic force, masseter inflammation, and nerve injury were categorized similarly: Genes encoding molecules for binding and catalytic activity were the top two categorized functions in all three conditions. DEGs enriched in different protein classes were also similar, except that orthodontic force-induced DEGs were apparently more enriched with transporters than the other two groups. Next, we compared the DEGs from the three groups at the individual gene level. Between orthodontic force-induced DEGs and masseter inflammation-induced DEGS, 490 genes exhibited significant changes in common (Figure 4(b)). Among these, only 159 genes showed changes correlated in the same direction (blue dots in Figure 4(c)), either upregulation or downregulation, whereas 331 genes exhibited anti-correlated changes in opposite directions (red dots in Figure 4(c)). Consequently, the fold changes of the DEGs from masseter inflammation and orthodontic force showed a weak correlation with a negative slope (r = −0.26, R2 = 0.04). On the other hand, the orthodontic force and nerve injury groups showed 290 overlapping DEGs (Figure 4(d)). Among these, 259 DEGs showed correlated changes (blue dots in Figure 4(e)), whereas only 31 DEGs showed anti-correlated changes (red dots in Figure 4(e)). Pearson correlation analysis of fold changes between orthodontic force and nerve injury showed apparently stronger correlation (r = 0.69, r2 = 0.47) than the correlation between orthodontic force and masseter inflammation. Orthodontic force, nerve injury and masseter inflammation showed common differential regulation of 116 genes. These genes showed relatively stronger correlation between orthodontic force and nerve injury (r = 0.71, r2 = 0.50). However, these genes showed poor correlation between orthodontic force and masseter inflammation (r = −0.37, r2 = 0.14) or between nerve injury and masseter inflammation (r = -0.33, r2 = 0.11).

Figure 4.

Comparison of DEGs in response to orthodontic force, masseter inflammation, and nerve injury. (a) Proportion of DEGs arising from orthodontic force (OF; red), masseter inflammation (MI; blue), and sciatic nerve injury (NI; gray) obtained from bulk RNAseq of entire ganglia. The results are categorized by molecular function and protein class using the PANTHER classification system. (a) Chromatin/chromatin-binding or -regulatory protein. Venn diagrams showing the number of DEGs overlapping between two conditions: (b) OF vs MI; (d) OF vs sciatic NI from bulk RNAseq; (f) OF vs trigeminal NI from single-cell RNAseq. Fold changes of the overlapping DEGs that are correlated (blue) or anti-correlated (red) between two conditions: (c) OF vs MI; (e) OF vs sciatic NI; (g) OF vs trigeminal NI.

Since our approach was to analyze bulk transcriptomes derived from entire ganglia, the data do not provide any information at the single cell level. To gain further insight on cell-specific gene regulation, we compared our data with those obtained from a single cell RNAseq study using TG from mice with partial infraorbital nerve transection.[16] A group of injured TG neurons was identified to have distinct transcriptome profiles from uninjured neurons. Among 2221 DEGs identified from a group of single TG neurons 2 days following trigeminal nerve injury, 150 genes overlapped with DEGs following orthodontic force. Among these, 132 genes showed correlated changes (blue dots in Figure 4(g)), whereas only 18 DEGs showed anti-correlated changes (red dots in Figure 4(g)). Pearson correlation analysis of fold changes between orthodontic force and trigeminal nerve injury showed high correlation (r = 0.69, r2 = 0.48), which is reminiscent of sciatic nerve injury in Figure 4(e).

Orthodontic force-induced DEGs overlap with DEGs from multiple neuronal and non-neuronal cells following nerve injury

A recent single cell RNAseq study using cells dissociated from DRG from mice with spinal nerve transection showed injury-induced transcriptional reprogramming in distinct subsets of neurons.[17] At two days after nerve injury, nine clusters of neurons showed transcriptional changes: Tac1+/Gpx3+ peptidergic nociceptors (PEP1), Tac1+/Hpca+ peptidergic nociceptors (PEP2), Mrgprd+ non-peptidergic nociceptors (NP), Sst+ pruriceptors (SST), Nefh+ A fibers including Aβ low-threshold mechanoreceptors (LTMR) (NF1), Pvalb+ proprioceptors (NF2), Cadps2+ Aδ-LTMRs (NF3), Fam19a4+/Th+ C-fiber LTMRs (cLTMR1), and Fam19a4+ putative cLTMR2.[17] In addition, six types of non-neuronal cells within ganglia also showed transcriptional changes: Apoe+ satellite glia, Mpz+ myelinating Schwann cells, Mpz-/Scn7a+ non-myelinating (Remak) Schwann cells, Macrophages, fibroblasts, and endothelial cells.[17] Comparison of DEGs following orthodontic force with cell type-specific DEGs following spinal nerve injury suggested that DEGs following orthodontic force involve gene changes in different subsets of sensory neurons as well as non-neuronal cells. Approximately 20-35% of cell type-specific DEGs from nine subsets of neurons induced by nerve injury overlap with DEGs following orthodontic force (Figure 5(a)). Approximately 88-95% of the overlapping genes showed correlated changes with DEGs by orthodontic force; R2 in Pearson correlation was between 0.41 and 0.64 (Figure 5(a) to (g)). DEGs by orthodontic force also overlapped with DEGs from non-neuronal cells following nerve injury. Approximately 11–30% of DEGs from six non-neuronal cell types induced by nerve injury overlaped with DEGs following orthodontic force. Approximately 85–98% of the overlapped genes showed correlated changes with DEGs by orthodontic force; R2 in Pearson correlation was between 0.29 and 0.66 (Figure 5(b) and (h)).

Figure 5.

Correlation of cell-type specific DEGs in response to spinal nerve injury with DEGs in response to orthodontic force. (a) Percent ratio of overlapping orthodontic force-induced DEGs to nerve injury-induced DEGs in nine neuronal subtypes; Tac1+/Gpx3+ peptidergic nociceptors (PEP1), Tac1+/Hpca+ peptidergic nociceptors (PEP2), Mrgprd+ non-peptidergic nociceptors (NP), Nefh+ A fibers including Aβ low-threshold mechanoreceptors (NF1), Pvalb+ proprioceptors (NF2), Cadps2+ Aδ-LTMRs (NF3), Sst+ pruriceptors (SST), Fam19a4+/Th+ C-fiber LTMRs (cLTMR1), and Fam19a4+ putative cLTMR2. The number of overlapping DEGs and R2 value from Pearson correlation in each cell type is shown below each bar. (b) Percent ratio of overlapping orthodontic force-induced DEGs to nerve injury-induced DEGs in six non-neuronal subtypes; Apoe+ satellite glia (Sg), Mpz+ myelinating Schwann cells (Sch_M), Mpz-/Scn7a+ non-myelinating (Remak) Schwann cells (Sch_N), macrophages (Mac), fibroblasts (Fb), and endothelial cells (EN). Fold changes of the overlapping DEGs that are correlated (blue) or anti-correlated (red) between orthodontic force-induced DEGs from the entire TG and nerve injury-induced cell-type specific DEGs in PEP2 (c), NP (d), NF1 (e), NF3 (f), cLTMR1 (g), and Schwann_N (f).

DEGs implicated in pain processing show correlation with DEGs from neuronal and non-neuronal cells induced by nerve injury

At the individual gene level, orthodontic force induced 84 DEGs that are implicated in pain processing (pain genes; Table 2). These genes include neuropeptides (e.g., Adcyap1 and Gal), neurotrophins (e.g., Bdnf) and neurotrophin receptors (e.g., Gfra1), cytokines (e.g., Csf1 and Cx3cl1) and cytokine receptors (e.g., Cxcr4 and Tnfrsf1a), transcription factors (e.g., Atf3 and Sox11), and ion channels (e.g., Trpa1 and Trpv2). Masseter inflammation also induces differential regulation of 41 of the 84 pain genes, among which 15 DEGs show correlative changes (e.g., Adcyap1, Bdnf, Gfra1, Trpa1, and Sox11), whereas 26 DEGs show anti-correlative changes. It is noteworthy that many of the pronociceptive genes that were upregulated following masseter inflammation were not upregulated following orthodontic force (e.g., Calca, Tac1, TRPV1, P2X3, Prkca, and Gfap). In contrast, nerve injury-induced DEGs include 33 of the 84 pain genes, among which 28 DEGs show correlative changes, whereas only five DEGs show anti-correlative changes. The proportions of correlative and anti-correlative changes in the masseter inflammation and nerve injury were significantly different (p < 0.0001 in Fisher’s exact test). It was also remarkable that orthodontic force downregulated multiple ion channel genes. Four voltage-gated sodium channel subunits (Scn11a, Scn1a, Scn2b, and Scn7a) were downregulated. Scn10a encoding Nav1.8, one of the most abundant sodium channel subunits in the TG, was also downregulated but did not reach significance (q = 0.09). In contrast, following masseter inflammation, all five of these sodium channel subunits were upregulated. Orthodontic force downregulated 11 voltage-gated potassium channel genes (Kcna2, Kcna6, Kcnb2, Kcnc4, Kcnd1, Kcnj10, Kcnj12, Kcnk5, Kcns1, Kcns3, and Kcnt2). Among these, 10 genes were also differentially regulated following masseter inflammation. However, only Kcnk5 was correlatively downregulated, whereas the remaining nine genes were anti-correlatively upregulated. In contrast, three of the 10 potassium channel genes were correlatively downregulated in the nerve injury dataset. Multiple DEGs from TG neurons assessed by single cell RNAseq following trigeminal nerve injury (e.g., Adcyap1, Atf3, Cacna2d1, Kcna2, and Kcnb2) also showed correlated changes with DEGs following orthodontic force.

Table 2.

Selected differentially regulated genes potentially implicated in hyperalgesia during experimental orthodontic tooth movement in mice.

Symbol	Gene name	Log2(FC)	Q	MI	SNI	TNI	NI-SC
							Neuronal	Non-neuronal
Adcyap1	Adenylate cyclase activating polypeptide 1	1.00	E	↑	↑	↑	P2,N,F1,F2,F3,S,C1,C2	All
Aqp4	Aquaporin 4	−0.34	B	↑	↓
Atf3	Activating transcription factor 3	1.06	E		↑	↑	All	1,3,4,5
Bdnf	Brain derived neurotrophic factor	0.95	E	↑	↑		F1,F2
Cacna2d1	Calcium channel voltage-dependent alpha2/delta subunit 1	0.27	B	↑	↓	↑	P2,N,F1,F2,F3,S,C1,C2	3,4,5
Cacng2	Calcium channel voltage-dependent gamma subunit 2	−0.22	A	↑	↓		S
Calb1	Calbindin 1	−0.29	D	↑	↓		F2
Casp1	Caspase 1	0.32	A	↓
Cckbr	Cholecystokinin B receptor	0.31	A		↑		P1,N,C1,C2
Csf1	Colony stimulating factor 1 (macrophage)	0.41	B		↑	↑	All
Cx3cl1	Chemokine (C-X3-C motif) ligand 1	−0.31	A
Cxcr4	Chemokine (C-X-C motif) receptor 4	−0.60	E	↓
Esr1	Estrogen receptor 1 (alpha)	−0.30	A
Gabrg2	Gamma-aminobutyric acid (GABA) A receptor subunit gamma 2	−0.22	B	↑	↓		P1,N,F1,F2,F3, C1
Gal	Galanin	0.39	B		↑	↑	All	1,3
Gfra1	Glial cell line derived neurotrophic factor family receptor alpha 1	0.47	E	↑	↑	↑	F1,F2,F3,S	1,3,5
Grik1	Glutamate receptor ionotropic kainite 1	−0.20	B		↓		N,C1,C2	2,6
Gstm1	Glutathione S-transferase mu 1	0.46	D	↑
Hdac1	Histone deacetylase 1	0.18	B	↓
Hdac4	Histone deacetylase 4	−0.21	B	↑
Kcna2	Potassium voltage-gated channel, shaker-related subfamily, member 2 (Kv1.2)	−0.25	B	↑	↓	↓	All	4,5,6
Kcnb2	Potassium voltage-gated channel, shab-related subfamily, member 2 (Kv2.2)	−0.24	B	↑	↓	↓	P1,N,F1,F3,S	2,5,6
Kcnj10	Potassium inwardly-rectifying channel subfamily J member 10 (Kir4.1)	−0.28	B	↑
Kcns1	Potassium voltage-gated channel, delayed-rectifier, subfamily S, member 1 (Kv9.1)	−0.25	B	↑	↓		F2,F3
Kcnt2	Potassium channel, subfamily T, member 2 (KCa4.2)	−0.35	B	↑
P2rx5	Purinergic receptor P2X ligand-gated ion channel 5	0.29	B
P2ry1	Purinergic receptor P2Y G-protein coupled 1	−0.21	A	↑	↓		F3,C1,C2
Plcd1	Phospholipase C delta 1	0.31	B		↑		S
Prkca	Protein kinase C alpha	−0.26	D	↑		↓	P1,P2,N,S	6
Procr	Protein C receptor endothelial	0.89	E	↓	↑		N,F2,F3,S, C1,P2	6
Sema6a	Semaphorin 6a	0.55	C		↑	↑	All	3,4,5
Sox11	SRY (sex determining region Y)-box 11	0.61	D	↑	↑	↑	All	All
Scn11a	Sodium channel voltage-gated type XI alpha (Nav1.9)	−0.16	A	↑	↓		P1,N,C1,C2	5,6
Scn1a	Sodium channel voltage-gated type I alpha (Nav1.1)	−0.15	A	↑	↓	↓
Tnfrsf1a	Tumor necrosis factor receptor superfamily member 1a	0.33	C	↓			F3
Trpa1	Transient receptor potential cation channel subfamily A member 1	0.37	C	↑	↓		N
Trpv2	Transient receptor potential cation channel subfamily V member 2	0.29	C			↑
Vgf	VGF nerve growth factor inducible	0.53	C		↑	↑	P1,P2,N,F1,F2,S,C2	1,3,5,6

FC, fold change; q, false discovery rate adjusted p value; MI, genes differentially regulated following masseter inflammation;[15] SNI, genes differentially regulated following sciatic nerve injury in bulk RNAseq;[16] TNI, genes differentially regulated following trigeminal nerve injury in single cell RNAseq;[18] NI-SC, genes differentially regulated following spinal nerve injury in single-cell RNAseq;[17] ↑, upregulated; ↓, downregulated; a, q<0.05; b, q<0.01; c, q<10−4; d, q<10−6; e, q<10−10; P1, PEP1; P2, PEP2; N, NP; F1, NF1; F2, NF2; F3, NF3; S, SST; C1, cLTMR1; C2, cLTMR2; 1, satellite glia; 2, myelinating Schwann cells; 3, non-myelinating (Remak) Schwann cells; 4, macrophages; 5, fibroblasts; 6, endothelial cells. blue, correlated; red, anti-correlated. Color coding presents the correlation intuitively.

Comparison of orthodontic force-induce DEGs with cell type-specific DEGs following nerve injury showed genes that are commonly or differentially regulated in neurochemically distinct subpopulations of primary afferents and different non-neuronal cells (Table 2; NI-SC). For example, Adcyap1 and sox11 were regulated in almost entire populations of neuronal and non-neuronal cells. Csf1 was upregulated in an entire subpopulation of neurons, but not in any non-neuronal cell. Some genes (e.g., Bdnf, Calb1, Cckbr, Grik1, Kcns1, or Plcd1) were differentially regulated only in specific subtypes of neurons, but not in non-neuronal cells.

Validation of DEGs using qPCR

We performed qPCR to confirm changes in gene expression determined by RNA-seq (Figure 6). Atf3, Adcyap1, Csf1, and Bdnf were significantly upregulated in Veh/Force compared to Veh/Sham. Upregulation of Atf3, Adcyap1, and Csf1 was significantly less in RTX/Force compared to Veh/Force, whereas upregulation of Bdnf in RTX/Force was significantly greater than RTX/Sham, but was not different from Veh/Force. Kcnj10 and Kcnt2 were significantly reduced in Veh/Force compared to Veh/Sham. The expression of TRPV1 was significantly lower in RTX/Sham and RTX/Force than in Veh/Sham and Veh/Force. Although the expression of Calca was substantially lower in RTX/Sham compared to Veh/Sham, the difference was not significant. Calca expression in Veh/Sham and Veh/Force was also not different. The expression of Gfap was not different among all groups. These results are largely consistent with the changes seen in the RNAseq assay results.

Figure 6.

Validation of gene expression using real-time qPCR. Transcripts of Atf3 (a), Adcyap1 (b), Csf1 (c), Bdnf (d), Kcnj10 (e), Kcnt2 (f), TRPV1 (g), Calca (h), and Gfap (i) were assessed using qPCR. Fold changes over the average of the VS group are plotted. *p < 0.05; **p < 0.005; ***p < 0.0005; ****p < 0.0001 in four pairs of Sidak post-hoc comparisons (VS vs. VF, VF vs. RF, RS vs. RF, VS vs. RS) following one-way ANOVA.

Discussion

Orthodontic force produces nocifensive behaviors in mice that last several days, and ablation of TRPV1+ trigeminal afferents substantially abolishes the behaviors[3]; therefore, we determined transcriptomic changes within TG after the application of orthodontic force, and assessed the effects of ablation of TRPV1+ afferents on the transcriptomic changes. Our data showed that application of orthodontic force resulted in differential regulation of approximately 6% of all genes in the TG within 2 days. Orthodontic force induced upregulation and downregulation of multiple terms in GO analysis. Cholesterol biosynthesis processes were upregulated, which is likely associated with axonal growth and regeneration.[34] Genes related to connective tissue organization and extracellular matrix were downregulated, which might suggest adaptation processes by periodontal nerve terminals in response to tissue remodeling. Biological processes relevant to neural functions were also altered: genes related to anion transport activity were upregulated, whereas genes associated with synaptic organization or overall ion channel activities, especially voltage-gated sodium and potassium channels, were downregulated. Although it is difficult to directly implicate these changes in pronociceptive or antinociceptive activities, altered anionic homeostasis potentially affects primary afferent excitability and nociception.[35] Decreased potassium channel activities, such as Kcna2, Kcnj10, Kcns1, and Kcnt2, are associated with increased pathological pain.[36-39] Orthodontic force also altered the expression of genes implicated in pain, such as Adcyap1, Trpa1, Csf1, and Bdnf.[40-43] These changes may be associated with increased nocifensive behaviors upon orthodontic force application, and their contribution to orthodontic pain needs to be determined. Orthodontic force upregulated Gal, which encodes galanin, consistent with a previous report.[44] TRPV1, CGRP, ASIC3, and P2X3 are known to be upregulated in TG following orthodontic force.[9,11-13] In our assay, however, these genes were not significantly upregulated at 2 days following the application of orthodontic force. Although the reasons for this discrepancy are not clear, it is possible that differences in experimental conditions (species, force magnitude, time points of assessment, methods for evaluation of RNA, among others) could have contributed.

When orthodontic force-induced DEGs were compared with DEGs induced by masseter inflammation or nerve injury, we found that the correlation of DEGs between orthodontic force and nerve injury had a larger proportion of correlated changes than the correlation between orthodontic force and masseter inflammation. The transcriptomic signature of three sets of nerve injuries was assessed 2 to 3 days after injury responses,[16-18] suggesting that these transcriptomic changes were early responses to the injury, and may not have been directly relevant to chronic neuropathic pain. Rather, these responses could be considered as immediate responses to nerve injury, including the pro-regenerative transcriptional program.[16-18] These early transcriptomic responses to nerve injury are preserved in sham surgery or minor injuries to facial skin, such as a shallow incision or scratch.[18] This transcriptomic signature includes upregulation (e.g., Atf3, Sox11, Sema6a, Csf1, Gfra1, and Gal) and downregulation (e.g., Grik1, Prkca, Trpc3, Scn10a, Scn1a, Calca, Tac1, and Kcnb2) of groups of genes in TG.[18] In our assay, orthodontic force induced consistent differential regulation of these genes, except for Scn10a, Calca, and Tac1. Orthodontic force produced by a coil spring generates continuous compression for a week or two on the periodontal ligament to which nociceptive afferent terminals are projected. Such mechanical irritation may cause injury to the periodontal ligament and afferent terminals, which can lead to transcriptomic changes similar to those seen in sensory ganglia after sciatic nerve injury. Since gene expression changes within sensory ganglia following nerve injury can cause the development of persistent neuropathic pain (e.g., Csf1[42]), it will be important to determine if these transcriptomic changes are associated with orthodontic pain. Downregulation of a group of voltage-gated potassium channels might contribute to increased pain. However, considering the relatively acute time course (1–5 days) of orthodontic pain in our experiments, it is unlikely that orthodontic force-induced transcriptomic changes in TG are associated with development of persistent neuropathic pain. Rather, the transcriptomic changes more likely represent an active reprogramming of injured neurons which functions to regenerate injured axons, as in neuropathic injury.[16-18] This regeneration program might contribute to preventing the development of nerve injury-induced long-lasting pain or paresthesia following orthodontic treatment, which is a rare complication of orthodontic treatment.[45] In contrast, the DEGs and biological pathways (e.g., satellite cell activation and leukocyte extravasation) that are strongly upregulated following masseter inflammation[15] are not enriched in orthodontic force-induced transcriptomic changes. This difference might be attributed to the different extents (a molar vs. masseter muscle) or sources of inflammation (neurogenic vs. CFA) in the two models.

Since we employed an approach using bulk RNAseq, we could not differentiate cell type-specific changes in gene expression. Comparison of our dataset with cell type-specific transcriptomic changes following spinal nerve injury obtained from single cell RNAseq, however, allowed us to infer cell type-specific changes in our dataset at least for the DEGs showing overlapping changes with DEGs caused by nerve injury. Based on the analysis, it would be reasonable to presume that orthodontic force induces differential regulation of genes in neurochemically diverse subtypes of trigeminal afferents including peptidergic, non-peptidergic afferents, Aβ LTMRs, Aδ LTMRs, and C-fiber LTMRs, as well as non-neuronal cells. We showed that RTX-induced ablation of TRPV1+ afferents in TG decreases pain-like behaviors induced by orthodontic force,[3] suggesting that nociception through TRPV1+ afferents can explain transcriptomic changes at least within TRPV1+ peptidergic afferents, PEP1 and PEP2.

We determined the role of TRPV1+ afferents in orthodontic force-induced transcriptomic changes by directly injecting RTX into the TG a week before the application of orthodontic force. In the mice that received the RTX injection, orthodontic force-induced differential regulation of gene expression was substantially prevented, which suggests that TRPV1+ afferents are necessary for orthodontic force-induced transcriptome changes within the TG. In view of the fact that transcriptomic changes following application of orthodontic force are not confined to TRPV1+ afferents, as discussed above, the results suggest that nociceptive inputs through TRPV1+ afferents lead to subsequent changes in gene expression not only in TRPV1+ neurons, but also in adjacent TRPV1-negative neurons and non-neuronal cells throughout the ganglia. Although 75% of periodontal afferents are TRPV1-negative neurons,[3] the results also suggest that the orthodontic force-induced transcriptomic changes in TRPV1-negative neurons are not likely due to nociceptive inputs through TRPV1-negative afferents projected to periodontium. In our assay, RTX injection alone caused differential expression of approximately 9% of genes, and these preceding changes in expression of a group of genes may occlude further transcriptomic changes by subsequent orthodontic force. This is unlikely, however, since orthodontic force-induced DEGs do not largely overlap with RTX-induced DEGs (Figure 2(g)). Alternatively, we cannot exclude the possibility that RTX injection induces transcriptomic changes in neurochemically diverse subsets of primary afferents, and somehow reduces their responses to the subsequent insult, in this case, orthodontic force. However, this possibility is slim since intraganglionic injection of RTX selectively diminishes the function of a subset of TRPV1+ peptidergic afferents without altering functions of other subsets of afferents involved in touch, proprioception and high-threshold mechanosensitive nociception.[20]

Although our study shows comprehensive changes in gene expression in TG following orthodontic tooth movement, we cannot exclude the possibility that the results include false-positive or -negative outcomes, perhaps because of the small sample size. Trends in differential expression might be affected by the comparison with data obtained from different species (rat vs mouse) or tissues (TG vs DRG) or approaches (Bulk RNAseq vs single cell RNAseq). The functional implications of changes in gene expression within TG following orthodontic force also need to be validated. In this study, we only focused on the role of TRPV1+ afferents, but not TRPV1 itself. Since TRPV1 mediates orthodontic pain,[3] it will be interesting to determine the role of TRPV1 in transcriptomic changes within TG.

In conclusion, we have shown that transcriptomic changes resembling nerve injury occur in the TG upon application of orthodontic force, and that the injury is driven by nociceptive inputs through TRPV1+ afferents. These results expand our understanding of the neuroplasticity in sensory ganglia upon the application of orthodontic force, and will likely lead to identification of targets for better management of pain and sensory disturbances during orthodontic treatment.