A human coronavirus OC43‐derived polypeptide causes neuropathic pain

Abstract Human coronaviruses have been recently implicated in neurological sequelae by insufficiently understood mechanisms. We here identify an amino acid sequence within the HCoV‐OC43 p65‐like protein homologous to the evolutionarily conserved motif of myelin basic protein (MBP). Because MBP‐derived peptide exposure in the sciatic nerve produces pronociceptive activity in female rodents, we examined whether a synthetic peptide derived from the homologous region of HCoV‐OC43 (OC43p) acts by molecular mimicry to promote neuropathic pain. OC43p, but not scrambled peptides, induces mechanical hypersensitivity in rats following intrasciatic injections. Transcriptome analyses of the corresponding spinal cords reveal upregulation of genes and signaling pathways with known nociception‐, immune‐, and cellular energy‐related activities. Affinity capture shows the association of OC43p with an Na+/K+‐transporting ATPase, providing a potential direct target and mechanistic insight into virus‐induced effects on energy homeostasis and the sensory neuraxis. We propose that HCoV‐OC43 polypeptides released during infection dysregulate normal nervous system functions through molecular mimicry of MBP, leading to mechanical hypersensitivity. Our findings might provide a new paradigm for virus‐induced neuropathic pain.

In this report, we identified a coronavirus HCoV-OC43-encoded amino acid sequence with a striking similarity to MBP  . Past research provided us with biochemical evidence that proprotein convertase furin and/or matrix metalloproteinase (MMP) inflammatory proteolysis releases cryptic MBP fragments implemented in multiple sclerosis (Shiryaev et al, 2009) and PNS injury (Kobayashi et al, 2008;Kim et al, 2012;Liu et al, 2012;Hong et al, 2017). We propose that MBP-like polyproteins generated during HCoV infection, similarly, mediate biological activities in PNS/CNS that promote neuropathic pain. Using a synthetic peptide specific to the HCoV-OC43/MBP 84-104 homologous region, we tested its activity in mechanosensitivity behavior upon sciatic nerve injection followed by RNA-seq, bioinformatics, and proteomic analyses of the unique DRG and spinal cord molecular signatures relative to the scrambled peptide.

Results
Identification of HCoV-OC43 fragment with high amino acid homology to MBP We have identified the nociceptive activity of the PNS injuryreleased MBP-derived peptides with highly conserved sequence motifs Ko et al, 2016;Hong et al, 2017;Chernov et al, 2018). The MBP 84-104 amino acid sequence, corresponding to positions 84-104 of the human classic MBP isoform 4, was used to search for homologous sequences in public depositories of clinical and environmental coronavirus isolates. Position-specific iterative PSI-BLAST algorithm (www.ncbi.nlm.nih.gov/blast) (Altschul et al, 1997;Altschul & Koonin, 1998) was applied with a standard threshold of 0.005. As a result, MBP 84-104 homologous regions were detected within six individual polyproteins from HCoV-OC43, canine respiratory coronavirus, bovine coronavirus E-AH65, Bat SARS-CoV Rm1/2004 and CoV279/2005, and nucleocapsid protein from porcine epidemic diarrhea (corona)virus ( Fig 1A). Because of its relevance to human disease, we selected the amino acid sequence of HCoV-OC43 ORF1ab polyprotein for further analysis. A highly homologous 12 amino acid sequence IVHFFKTFTTST (OC43 656-668 ) is localized at positions 656-668 of ORF1ab (at positions 409-421 of p65-like, also known as mouse hepatitis virus (MHV) p65-like protein) ( Fig 1B).

Structural similarity revealed by predictive analysis
Although MBP is an intrinsically disordered protein, several regions exhibit secondary structure, including the conserved MBP 84-104 (Ahmed et al, 2012). Predictions conducted with MBP isoform 4 sequence using AlphaFold2 (Jumper et al, 2021;Tunyasuvunakool et al, 2021) supported these conclusions. Because experimentally resolved 3D structures of HCoV-OC43 p65-like were not available, we predicted structures of both OC43 656-668 and MBP 84-104 using Alpha-Fold2 using identical computational parameters. 3D structure models with the highest scores were compared by alignment. A remarkable structural similarity was observed between MBP 84-104 and OC43 656-668 ( Fig 1C). Each peptide is folded into a characteristic N-terminal α-helix and unstructured C-terminal tail. The aligned α-helices comprised the VHFFK motif, including the invariable histidine-89 (Chernov et al, 2018). Less conserved C-terminal tails consisted of multiple threonine/serine residues in close proximity to α-helix. When AlphaFold2 predictions were performed with the full-length proteins (human classic MBP, isoform 4, NP_001020263; HCoV-OC43 p65-like, YP_009555238.1), predicted structures of the corresponding protein regions and the respective peptides were highly similar. We concluded that OC43 656-668 exhibited high amino acid sequence homology and structural similarity to the pronociceptive MBP 84-104 peptide.
HCoV-OC43-derived peptide causes persistent mechanical hypersensitivity in female rats We have shown repeatedly that MBP 84-104 peptide produces mechanical allodynia with no effect on thermal sensitivity Ko et al, 2016) likely due to its myelin-dependent pronociceptive activity on myelinated A-afferents, sparing unmyelinated heatsensitive nociceptors . Thus, to test the ability of the OC43-derived peptide to regulate mechanical hypersensitivity characteristics to the homologous MBP 84-104 , we used a synthetic 20-amino acid OC43 653-673 peptide (OC43p, VSKIVHFFKTFTTSTALAFA), and scrambled peptides OC43p-SCR1 and OC43p-S2 designed to mismatch the MBP 84-104 amino acid sequence. Female rodents display robust mechanical hypersensitivity to intrasciatic MBP 84-104 relative to males (Chernov et al, 2020). Female rats received a single bolus intrasciatic injection (Fig 2A) of OC43 653-673 , scrambled peptides (10 µg in 5 µl, each), or PBS vehicle (5 µl) (n = 6/group), followed by von Frey testing. The rats displayed a significant reduction in the mechanical force required to evoke hind paw withdrawal after OC43 653-673 injection, and the effect was sustained during the 3-week observation period (Fig 2B and C). In contrast, the withdrawal thresholds remained significantly higher in rats injected with OC43p-SCR and alternative OC43p-S2 peptides, or PBS. In agreement, our prior studies found no hypersensitivity arising from scramble MBP 84-104 peptide sequences and PBS vehicle in the equivalent experimental designs Ko et al, 2016;Hong et al, 2017;Chernov et al, 2020). No significant contralateral hypersensitivity was observed in response to either peptide. Unstimulated pain-like behavior was measured according to the method described by Attal et al (1990) with modifications. Rats injected with OC43p exhibited slightly higher unstimulated painlike behavioral indices, although differences from the control animals were not statistically significant ( Fig 2D). We concluded that OC43p, like MBP 84-104 , induced a robust and sustained pain mechanical hypersensitivity in female rats.

Nerve injections of OC43p induced vast transcriptional changes in the spinal cord
Next, genome-wide transcriptomes were compared in animal groups injected with either OC43p or OC43p-SCR control. Total RNAs from ipsilateral L4-5 DRG and lumbar spinal cord (dorsal quarter) were collected at day 21 post-injection and analyzed by RNA-seq. We detected 17 up-and 21 downregulated differentially expressed genes (DEGs, adjusted P (P adj ) < 0.1) in DRG in the OC43p group relative to the scrambled-injected group (Fig 3A). In the spinal cord, 724 up-and 160 downregulated DEGs (|Log 2 (fold change (FC))| > 1, P adj < 0.1) were recorded ( Fig 3B, Dataset EV1). The principal component analysis (PCA) (Fig 4A) attributed 87.8% of the variance (PC1) to the effect of peptide injections, highlighting LOC108348215, Col8a1, Six1, Slc26a7, Kcnj13, and Tlr12 genes as the most potent drivers of variation.
In female rats, which are susceptible to MBP peptide-induced pain (Chernov et al, 2020), we detected a large number of X-linked upregulated DEGs in response to OC43p (Fig 4C). The expression of LOC100911498 (a homolog of XIST non-coding RNA in rats), a marker for Xi in females, exhibited upregulation. Another epigenetic factor crucial for the Xi state, the X-linked chromatin remodeling helicase II (Atrx), also demonstrated robust upregulation. Taken together, we concluded that OC43p induced multifaceted transcriptional responses in the PNS/CNS consistent with pronociceptive and proinflammatory signaling.
Gene ontology (GO) analysis identified pronociceptive signaling pathways activated by OC43p  A MBP 84-104 homology search in Coronaviridae (taxid: 11118) datasets by position-specific iterative (PSI)-BLAST is displayed on a distance tree with a maximum sequence difference of 0.85. B Schematic organization of HCoV-OC43 coronavirus ORF1ab polyprotein. Triangles indicate the canonical cleavage sites. Amino acid sequence alignment of homologous MBP 84-104 and OC43p is plotted using Clustal Omega. C 3D structure of MBP 84-104 (blue) and OC43p (yellow) according to AlphaFold2 predictions. Key amino acid residues are labeled by single-letter codes and numbers corresponding to the relative position in the human classic MBP, isoform 4 (NP_001020263), and HCoV-OC43 p65-like (YP_009555238.1). Structures are visualized in PyMOL.
ª 2022 The Authors EMBO reports 23: e54069 | 2022 |Log 2 FC| > 1 and P adj < 0.1. Pathogen response-specific pathways, FXR/RXR, LXR/RXR, T-cell receptor signaling, and immune responsespecific pathways were affected in DRG ( Fig 5A). In the spinal cord, due to a robust upregulation of multiple voltage-and ligand-gated ion channels, signaling pathways involving neural signal transduction were predictably activated ( Fig 5B). Notably, the activation of A A schematic of the injections into the sciatic nerve followed by ipsilateral DRG and dorsal spinal cord tissue analysis. B von Frey testing in female rats (n = 6/group) at 1-21 days after injections of OC43p peptide, respective control scrambled OC43p-SCR, OC43p-S2 (10 μg in 5 μl, each) peptides, and PBS vehicle. Responses were recorded in ipsilateral hind paws. Mean tactile withdrawal thresholds are in gram force (g) AE standard deviation; two-way analysis of variance (ANOVA) with Bonferroni post hoc test: **P ≤ 0.005; ***P ≤ 0.0005; and ****P ≤ 0.00005. C Areas under the curve (AUC) were calculated for days 1-21 (n = 6 animals/group). Bars show the mean AUC and standard deviations (error bars) for each injection group. Data were analyzed by two-way ANOVA with Tukey's post hoc test: ***P ≤ 0.0005; ****P ≤ 0.00005. D Observational assays of unstimulated pain-like behavior. Assays were conducted in female rats (n = 6/group) after injection of OC43p (red), OC43p-SCR (purple), or PBS (green) on days 4, 14, and 20 post-injection. Each animal was video-recorded for 2 min three times within a 2 h period of time. Hind paw positions were scored to calculate unstimulated pain-like behavioral indices. Index means and standard deviation are shown; two-way ANOVA with Tukey's post hoc test was used for group comparisons.  synaptogenesis signaling, CREB signaling in neurons, neuropathic pain signaling in dorsal horn neurons, glutamate receptor signaling, calcium signaling, opioid signaling, and endocannabinoid neuronal synapse pathway was expected in connection with persistent pain hypersensitivity demonstrated by behavioral tests. The activation of estrogen receptor signaling and androgen signaling pathways due to robust Esr2 and Ar upregulation was recorded. The activation of TLR signaling, neuroinflammation signaling pathway, IL-2/IL-6/IL-8, and PI3K signaling in B lymphocytes was low to moderate relative to other pathways. It is worth noting the activation of the long-coding RNA HOTAIR regulatory pathway was previously not associated with pain signaling. Pathways related to mitochondrial function, metabolic pathways, and protein synthesis demonstrated a decline. Remarkably, prolactin signaling was elevated (Fig 5B) due to more than twofold elevation of genes encoding prolactin receptor dimer (Prlp/Prlr), Irs1, potassium A Principal component (PC) analysis of DEGs in OC43p (red) and OC43p-SCR (blue) groups (n = 3/group). B Hierarchical clustering plot of 548 significant upregulated DEGs (Log 2 FC > 1, P adj < 0.1, n = 3/group). Heatmap color scheme corresponds to logarithms of variance stabilized counts. C Hierarchical clustering plot of X-linked upregulated DEGs (Log 2 FC > 1, P adj < 0.1, n = 3/group). channel Kcnma1, PI3K family member Pik3c2a, protein kinase C epsilon type (PKCε), Irs1, and Socs4 (Fig 6A). The most significant DEGs were attributed to molecular functions (MFs) using the GO database to enrich the biological interpretation. The best-match average (BMA) distance plot represents the proximity of the eleven MF clusters (Fig 6B) summarized on the distance heatmap ( Fig 6C). A broad range of cellular functions was affected, including transcriptional regulation, ribosome activity, binding to fatty acids, hormones, and ECM (collagen and fibronectin). We identified signaling pathways, including voltage-gated ion channel activity, ligand-gated channel activity, and cytokine receptor activity directly relevant to pain hypersensitivity and nociception.

X-linked DEGs
OC43p affinity to the Na + /K + -transporting ATPase complex To identify proteins potentially interacting with OC43p, we conducted the affinity capture in rat protein lysates (spinal cord) using biotinylated OC43p and OC43p-SCR peptides bound to paramagnetic beads. To reduce non-specific binding, in the protein lysates we preincubated with OC43p-SCR beads. Beads were removed and affinity capture was conducted with OC43p-bound beads. Bound proteins were digested with trypsin, digestion products were separated by liquid chromatography, and analyzed by mass spectrometry. Strikingly, Atp1a1, Atp1a2, Atp1a3, and Atp1b1 proteins produced high significance scores and peptide coverages 32% to 57% (Fig 7A). These proteins represent subunits of the Na + /K + -transporting ATPase complexes (Fig 7B) (http:// geneontology.org). According to RNA-seq, the expression of the ATPase subunits is high in the spinal cord and DRG but lower in sciatic nerves (Fig 7C). Protein immunoblotting further confirmed that Atp1a1 and Atp1a2 interaction is specific to OC43p (Fig 7D; Atp1a3 and Atp1b1 were not probed).

Discussion
Neurotropic viruses (Johnson, 1999;Dahm et al, 2016;Maximova et al, 2021), including HCoVs (Burks et al, 1980;Talbot et al, 1993;Arbour & Talbot, 1998;Arbour et al, 1999;Edwards et al, 2000;Glass et al, 2004;St-Jean et al, 2004;Jacomy et al, 2006;Dub e et al, 2018), expose cells to overwhelming quantities of viral proteins and, due to molecular mimicry with host proteins (Roos, 1983;Wege et al, 1983;Weise & Carnegie, 1988;Adelmann & Linington, 1992;Stohlman & Hinton, 2001;Talbot et al, 2001;Getts et al, 2013;Savarin & Bergmann, 2017), may disrupt cellular protein-protein/ RNA/DNA/lipids interactions in the host. Expectedly, HCoV sequence evolution continues to introduce novel amino acid sequence patterns, which may eventually include novel HCoV strains, such as SARS-COV-2. Molecular mimicry can assist the virus in hijacking host-specific functions to (i) mediate immune and neuroimmune responses in the upstream, uninfected regions of the nervous system by axonal trafficking; and (ii) directly affect transcriptional programs in the PNS/CNS neurons in favor of virus survival and immune system evasion. Interference with the PNS/CNS regulatory networks leads to detrimental long-term neurological health outcomes (Burks et al, 1980;Cook & Dowling, 1980;Talbot et al, 1993;Edwards et  As a proof of concept, the perspective synthetic peptide used in this study was derived based on strong sequence and structure homology to MBP 84-104 . The OC43/MBP homologous region is localized at positions 407-422 of the MHV p65-like protein of HCoV-OC43 (GenBank ID YP_009555247), and sequences with such identity are unknown in other (corona)viruses to date. The expression and proteolytic processing of this protein product during viral infections as part of the pp1ab polyprotein were demonstrated in cells infected with coronaviruses and related viruses (reviewed in Weiss et al, 1994). We propose that nociceptive activity can be exhibited by polypeptides of varying lengths. Our investigation centers on the MBP homologous sequence accessible for interaction with respective host protein targets. It is noteworthy that peptide epitopes were used in a conceptually similar study that established high-affinity molecular mimicry based on a short amino acid homologous motif shared by the Epstein-Barr virus-encoded transcription factor EBNA1, and host-encoded GlialCAM protein was implicated in multiple scleroses (Lanz et al, 2022). Future translational and clinical studies in patients with diagnosed infections can ascertain the precise identity of the pronociceptive viral polypeptides.
In MBP, the algesic sequence is buried inside the intact protein and becomes exposed for interaction after proteolytic degradation of MBP by cellular peptidases as we shown previously. If the MBP-like viral sequences are readily exposed to the interface of unprocessed ◀ Figure 5. Canonical pathways affected by OC43p.
A, B Canonical pathways (A, DRG and B, spinal cord) identified by IPA are ranked by P (-log 10 P). Positive and negative z-scores (in the spinal cord) indicate upregulation or downregulation, respectively, according to the color scale. Bar colors correspond to significant activation (orange), deactivation (blue), or no change (gray) in pathway regulation. or partially processed viral proteins, further proteolytic cleavage of the viral polypeptides to shorter peptides is not a prerequisite of molecular mimicry activity. We hypothesize that viral infectionrelated proteolytic mechanisms, such as the inflammatory proprotein convertase/MMP proteolytic pathway (Shiryaev et al, 2009), and viral intrinsic proteinases, can further stimulate a release of cryptic viral peptidic fragments with biological activities. Intriguingly, recent evidence that MMP2/MMP9 can activate the SARS-CoV-2 fusion (preprint: Benlarbi et al, 2022) in cells expressing high levels of MMPs provides a connection with the role of these MMPs in nerve injury and pain (Shubayev et al, 2006;Chattopadhyay et al, 2007;Kobayashi et al, 2008;Kim et al, 2012;Liu et al, 2012;Remacle et al, 2015Remacle et al, , 2018a. The N-terminal invariable 87-(V/I)VHFFK-92 motif of MBP 84-104 and OC43p is folded into structurally similar α-helices. The unstructured C-terminal tails included 3-5 threonine and serine residues, subject to enzymatic phosphorylation by MAPK, CDK5, GSK3 (Pelech, 1995;Chernov et al, 2018), and other kinases. Dynamic phosphorylation/dephosphorylation by cellular kinases can engage a regulatory switch critical for the peptide's biological activity. We propose that the bipartite characteristics of the MBP 84-104 and OC43p peptides may be crucial to promoting pronociceptive and  Figure 7. OC43p affinity to the Na + /K + -transporting ATPase complex.

▸
A Mass spectrometry identification of the Atp1 family of proteins captured by biotinylated OC43p peptide in rat dorsal spinal cords from female rats (n = 2/groups). Negative log 10 Pconfidence score and peptide sequence coverage (%) were calculated in PEAKS Studio TM 8.5 (BSI). B Relationship of ATPase subunits identified in string-db.org. Lines indicate evidence-based interactions. C ATPase gene family expression in sciatic nerves, DRG, and spinal cords of naïve rats. TPM (transcripts per kilobase million) measures the transcription frequency of a specific gene by RNA-seq (n = 3/group). Blue to red color intensity corresponds to Log 2 (TPM). D Affinity capture of proteins by OC43p-and OC43-SCR-bound magnetic beads. Specific reactivity to Atp1a1 and Atp1a2 was detected by protein immunoblotting using specific primary antibodies. OC43p-SCR-coated and uncoated beads (Beads) were used to control non-specific binding. Spinal cords were isolated from male (M) and female (F) rats (n = 2/group). Input protein lysates (right side, in duplicates) served as loading controls. Approximate molecular weights (in kDa) are indicated by arrows.
other neuropathological activity. We identified an intriguing affinity of OC43p to the Na + /K + -transporting ATPase complex responsible for electrochemical cation gradient across the plasma membrane and electrical excitability in the nervous system. The ATPase's multiple subunits have been implicated in Charcot-Marie-Tooth disease, peripheral neuropathies, neuromuscular disorders (reviewed in Clausen et al, 2017), and inflammationinduced mechanical allodynia (Wang et al, 2015). In the context of virus-host interaction, based on our transcriptomics analysis and affinity capture assay we hypothesize that OC43p can directly affect ion transport and aberrant neuroplasticity leading to persistent mechanical allodynia. We demonstrated that OC43p induced persistent pain hypersensitivity in female rats. Females are more susceptible to developing chronic pain states as compared to males, including MBPinduced pain in rodents (Chernov et al, 2020). Bioinformatics analysis of RNA-seq data illuminated pronociceptive transcriptional changes in the dorsal spinal cord established within three weeks after sciatic nerve injection. Supporting the observed pain effects, peripheral terminals of nociceptor neurons and spinal higher-order neurons in the dorsal spinal cord increased expression of a broad spectrum of ion channels in response to OC43p. This observation was consistent with its potential role in neuron excitability in the pain sensation (Suzuki & Dickenson, 2000;Julius & Basbaum, 2001;Kidd & Urban, 2001). Our prior observations of mechanical hypersensitivity were not accompanied by thermal hyperalgesia in response to the homologous MBP 84-104 peptide Ko et al, 2016) consistent with the model of a myelin-dependent pronociceptive activity on myelinated Aafferents, sparing unmyelinated heat-sensitive nociceptors .
The proallodynic MBP 84-104 activity and the downstream signaling are sexually dimorphic (Chernov et al, 2020). In agreement, upregulation of sex hormone receptors Esr2 and Ar in females in response to OC43p predicts activation of the estrogen and androgen signaling pathways, respectively, and potentially virus-induced hypersensitivity. Upregulated cytokine receptor genes and respective signaling pathways outlined the mechanistic link between Esr2, neuroimmune properties of glia, and neuronal excitability as a characteristic of sustained neuropathic states. Accordingly, activation of prolactin signaling (Patil et al, 2019) and the regulatory role of the X chromosome in immunity, (neuro)-inflammation, and neuropathic pain (Syrett et al, 2019;Shenoda et al, 2021;Tang et al, 2021) contributed to the femalespecific pain response. Accordingly, we observed unexpected upregulation of the XIST homolog and other X-linked epigenetic factors. The role of sexual dimorphism in coronavirus-related chronic pain requires focused investigation using both female and male animals.
To summarize, our data strongly support the pronociceptive biological activity of OC43p due to molecular mimicry mechanisms to a neural-specific host protein, MBP. We propose that HCoVs evolve their encoded protein sequences to mimic host proteins in order to hijack cellular programs related to immune, metabolic, and cellular energy functions in the somatosensory nervous system. The HCoV peptide's identification in clinical specimens, their pronociceptive properties, and pathobiochemical processes of their release constitute topics of perspective research.

Protein structure predictions
AlphaFold2 (Jumper et al, 2021) in ColabFold (preprint: Mirdita et al, 2021 with multiple sequence alignments generated by MMseqs2 was used to predict 3D structure. Predicted models were aligned, processed, and visualized in PyMOL (www.pymol.org).

Animal procedures
Sprague Dawley female rats (8-10 weeks old) were obtained from Envigo and housed in a temperature-controlled room (~22°C), on a 12-h light/dark cycle, and with free access to food and water. All procedures were conducted during the daytime. Under isoflurane anesthesia, the common sciatic nerve was exposed unilaterally at the mid-thigh level. A single bolus injection of the peptides (10 µg in 5 μl vehicle) into a nerve fascicle was performed using a 33gauge needle on a Hamilton syringe. All animal procedures were performed according to the Policy on Humane Care and Use of Laboratory Animals and the protocol approved by the Institutional Animal Care and Use Committee at the VA San Diego Healthcare System. Weekly weight measurements and daily assessments of the hydration status and post-operative wound infection were conducted. Animals that developed surgery-related abnormalities were excluded from the study.

Behavioral tests
All behavioral measurements were taken by a tester blinded to the experimental groups. Animal groups were formed randomly. von Frey testing was performed before and at the indicated time points after peptide injections (n = 6/group). Rats were placed in individual compartments with a wire mesh bottom. von Frey filaments (0.41-15.2 g, Stoelting) were applied perpendicularly to the midhind paw for 4-6 s. A withdrawal response was recorded by an experimenter blinded to the groups. The 50% probability of withdrawal threshold was determined by up-down method as described previously (Chernov et al, 2020) using software developed in R (https://github.com/chernov-lab/VonFreyTest). Areas under the curve were calculated using Prism 9 (GraphPad). Unstimulated pain-like behavior was analyzed on days 4, 14, and 20 post-injection as previously described (Attal et al, 1990;Paulson et al, 2002;Chattopadhyay et al, 2007) with modifications. Each animal was video-recorded for 2 min three times within a 2 h period of time. Positions of the injected hind paw were continuously rated according to the scoring system: 0, the paw was placed normally on the floor; 1, the paw was placed lightly on the floor, and the toes were in a ventroflexed position; 2, only the inner edge of the paw was placed on the floor; 3, only the heel was placed on the floor, and the hind paw was inverted; 4, the hind paw was elevated; and 5, animal licked the hind paw. Scoring data were interpreted using custom Java software. Unstimulated painlike behavioral indices were calculated by the time interval the rat spent in each behavior multiplied by weighting factors, and divided by the length of the observational period according to the formula: where t 0 -t 5 is the time duration (s).

Samples
Tissues (DRG and spinal cord, lumbar (L)1-6, quartered) were placed in 500 µl RNAlater, left at 4°C overnight, and then stored at −20°C. All sample groups were processed synchronously to minimize batch effects.

RNA purification
Tissues were homogenized, and total RNAs were purified using RNeasy reagents. RNA concentration and integrity were determined using Qubit 4 and Bioanalyzer, respectively. 500 ng of RNA (3 replicates/group) with RIN ≥ 7.0 was used for RNA-seq.

RNA-seq
RNA-seq library preparations and sequencing were performed at the Genomics High Throughput Facility (University of California, Irvine). In brief, mRNA libraries were generated following the TruSeq Stranded mRNA library preparation protocol (Illumina).
Poly-A-enriched mRNAs were purified using poly-T oligo coupled magnetic beads, followed by mRNA fragmentation, first and second strands synthesis, cleaning on AMPure XP beads, and 3'adenylation. Ligation of TruSeq dual-index adapters was used for barcoding. The quality of RNA-seq libraries was validated using qPCR. Libraries were sized on an Agilent Bioanalyzer DNA highsensitivity chip and normalized. RNA-seq was performed using the paired-end 100 cycle program on the NovaSeq 6000 system. Base calls were recorded and converted to FASTQ files containing sequencing reads and the corresponding quality scores using Illumina software. Sequencing was conducted until at least 25 million paired-end reads per sample were acquired.

Data processing
The data analysis workflow is schematically presented in Fig EV1. FASTQ files were filtered to remove low-quality bases, TruSeq dualindex adapter sequences, and unpaired reads using Trimmomatic (Bolger et al, 2014). Transcript-level quantification was performed using Salmon (Patro et al, 2017) in quasi-mapping mode using the Rat genome version R7. To correct systematic biases commonly present in RNA-seq data, -seqBias and -gcBias options were applied. Transcript-to gene-level conversion was done using Tximeta (Love et al, 2020). RNA-seq coverage and data quality were assessed using MultiQC (Ewels et al, 2016). Gene count matrices were imported into the DESeq2 package (Love et al, 2014). Outliers were identified by Cook's distance method and excluded from further analysis. Dataset's normalization was conducted using trimmed M-values (TMM) included in the DESeq2 package. Log 2 FC was calculated using the Wald test. The adjusted (shrunken) Log 2 FC values were calculated using the adaptive t-prior apeglm method (Zhu et al, 2019). Significant DEGs were identified by P adj values below a false discovery rate cutoff (P adj < 0.05) (Dataset EV1). P adj < 0.05 was used in downstream analyses unless otherwise noted. Batch effects were controlled using remove-BatchEffect (Ritchie et al, 2015) and RUVseq (Risso et al, 2014) functions.

Signaling pathway analysis
Bioinformatics tools used for the processing of RNA-seq data are listed in Table EV1. Ingenuity Pathway Analysis based on the causal network approach (Kr€ amer et al, 2014) was used to predict signaling pathway regulation. The activation directionality was estimated based on z-scores. Gene ontology analysis was performed using ViSEAGO package and other Bioconductor tools.
All affinity capture steps were performed at 4°C. Frozen rat dorsal spinal cord tissues were submerged in 300 µl of lysis buffer (affinity capture buffer supplemented with 50 mM octylthioglucoside (OTG)), and homogenized using BioMasher microhomogenizers for 1 min, followed by centrifugation in QiaShredder units (Qiagen) for 5 min at 21,000 g. Protein lysates were diluted with affinity capture buffer and pre-adsorbed with 100 µl Dynabeads MyOne Streptavidin T1 (10 mg/ml) for 2 h with agitation. Lysates were incubated overnight with 50 µl of respective peptide-bound beads with agitation. Beads were washed seven times with 750 µl of TBST. For mass spectrometry, beads were washed three times with 750 µl PBS. For immunoblotting, proteins were eluted by heating at 70°C for 10 min in 100 µl of NuPAGE LDS Sample Buffer (Thermo Fisher Scientific) supplemented with 50 mM 1,4dithiothreitol. Protein concentrations were measured using a bicinchoninic acid assay.
Liquid chromatography and mass spectrometry LC/MS was performed at the Biomolecular Mass Spectrometry Facility (University of California San Diego). In brief, affinity-captured proteins were trypsin-digested, and peptides were separated by liquid chromatography for 1.5 h using a reverse-phase C18 gradient. Mass spectrometry was performed using Orbitrap Fusion TM Lumos Tribrid (Thermo Fisher Scientific). Proteomics data were analyzed using PEAKS Studio TM 8.5 (BSI).

Immunoblotting
Proteins were separated on Bolt 4-12% Bis-Tris protein gels (Thermo Fisher Scientific) and transferred onto a PVDF membrane (Thermo Fisher Scientific) following the manufacturer's instructions. The membrane was blocked in 5% non-fat milk for 1 h and incubated for 18 h at 4°C with specific primary antibodies. Membranes were washed 6 times with TBST and incubated for 1 h at ambient temperature with respective secondary HRP-conjugated antibodies. Membranes were washed six times with TBST, and chemifluorescence signals were developed using a SuperSignal West Dura Extended Duration Substrate kit (Thermo Fisher Scientific) and documented on X-ray films.
Expanded View for this article is available online.

Acknowledgements
The research was supported by the National Institutes of Health (NIH) R01 DE022757 (to VIS) and the Department of Veterans Affairs Merit Review Award 5I01BX000638 (to VIS). The authors wish to thank Dr. F abio Catroli Andrijauskas, Megh Jariwala, and Mila Angert for their expert technical assistance. RNA-seq was made possible, in part, through access to the Genomics High Throughput Facility Shared Resource of the Cancer Center Support Grant (P30CA-062203) at the University of California, Irvine, and NIH shared instrumentation grants 1S10RR025496-01, 1S10OD010794-01, and 1S10OD021718-01. Mass spectrometry was made possible, in part, through access to the Biomolecular and Proteomics Mass Spectrometry Facility (UCSD) supported by the NIH shared instrumentation grant S10 OD021724. The content is solely the authors' responsibility and does not necessarily represent the official views of the funding agencies. In addition to the CRediT author contributions listed above, the contributions in detail are:

Author contributions
VIS conceptualized the study, provided resources, designed the methodology (animal model), wrote, reviewed, and edited the manuscript, contributed to project administration, and acquired funding. JD investigated animal procedures and behavioral assays. GFC investigated animal procedures, behavioral assays, and analysis of the study. AVC conceptualized the study, designed the methodology (RNA-seq, affinity capture, and bioinformatics), contributed to the software, performed formal analysis and data curation, investigated and visualized the study, wrote the original draft, and wrote, reviewed, and edited the manuscript.

Disclosure and competing interests statement
The authors declare that they have no conflict of interest.