A human coronavirus OC43-derived polypeptide causes neuropathic pain.

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.

Introduction

Human coronaviruses (HCoV) causing mild and severe respiratory distress syndromes show evidence for the peripheral and central nervous systems (PNS/CNS) involvement (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é et al, 2018), potentially contributing to neurological conditions (Boziki et al, 2020; Gutiérrez‐Ortiz et al, 2020; Koralnik & Tyler, 2020; Manji et al, 2020; Montalvan et al, 2020; Romoli et al, 2020; Troyer et al, 2020; Ermis et al, 2021), including Guillain–Barré syndrome (Kilinc et al, 2020; Koralnik & Tyler, 2020; Montalvan et al, 2020; Sancho‐Saldaña et al, 2020; Zhao et al, 2020; Koike et al, 2021), multiple sclerosis (Burks et al, 1980; Cook & Dowling, 1980; Talbot et al, 1993; Edwards et al, 2000; Boziki et al, 2020), and states of neuropathic pain (Kemp et al, 2020; Mao et al, 2020; Widyadharma et al, 2020; Attal et al, 2021; McFarland et al, 2021; Şahin et al, 2021). The virus‐mediated pathologies can be accompanied by damage to the myelin sheath of the nervous system and cause rapid‐onset demyelination (Croxford et al, 2005).

Cationic myelin basic protein (MBP) controls myelin compaction, cytoskeletal interactions, and calcium homeostasis through electrostatic interactions with anionic lipids and proteins (Boggs & Moscarello, 1978; Boggs, 2006). MBP is also a major autoantigen contributing to autoimmune demyelinating disorders, including Guillain–Barré syndrome and multiple sclerosis (Kadlubowski & Hughes, 1979; Musse et al, 2006). Molecular mimicry between host and viral proteins (Roos, 1983; Weise & Carnegie, 1988; Adelmann & Linington, 1992; Stohlman & Hinton, 2001; Getts et al, 2013), including myelin sheath and HCoV proteins (Wege et al, 1983; Talbot et al, 2001; Savarin & Bergmann, 2017), are thought to contribute to the etiology of these conditions.

Our earlier work (Kobayashi et al, 2008; Kim et al, 2012; Liu et al, 2012; Ko et al, 2016; Shubayev et al, 2016, 2018; Hong et al, 2017; Chernov et al, 2018, 2020; Remacle et al, 2018b) implicated immunodominant MBP regions, proteolytically released after PNS damage, in initiating mechanical hypersensitivity through autoreactivity targeted at myelin on mechanosensory neurons. MBP84–104 peptide injection into an intact sciatic nerve was sufficient to induce sustained pain (Liu et al, 2012; Ko et al, 2016) via transcriptional reprogramming of metabolic, pronociceptive, and inflammatory signaling in the segmental dorsal root ganglia (DRG) and spinal cord in sex‐specific manner (Chernov et al, 2020). MBP84–104 amino acid sequence conservation is critical for its interactions, trafficking, and pronociceptive activity (Chernov et al, 2018).

In this report, we identified a coronavirus HCoV‐OC43‐encoded amino acid sequence with a striking similarity to MBP84–104. 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/MBP84–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 injury‐released MBP‐derived peptides with highly conserved sequence motifs (Liu et al, 2012; Ko et al, 2016; Hong et al, 2017; Chernov et al, 2018). The MBP84–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, MBP84–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 (OC43656–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).

Figure 1

Amino acid sequence homology and structural similarity of OC43p and MBP84–104

MBP84–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.

Schematic organization of HCoV‐OC43 coronavirus ORF1ab polyprotein. Triangles indicate the canonical cleavage sites. Amino acid sequence alignment of homologous MBP84–104 and OC43p is plotted using Clustal Omega.

3D structure of MBP84–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.

Structural similarity revealed by predictive analysis

Although MBP is an intrinsically disordered protein, several regions exhibit secondary structure, including the conserved MBP84–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 OC43656–668 and MBP84–104 using AlphaFold2 using identical computational parameters. 3D structure models with the highest scores were compared by alignment. A remarkable structural similarity was observed between MBP84–104 and OC43656–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 OC43656–668 exhibited high amino acid sequence homology and structural similarity to the pronociceptive MBP84–104 peptide.

HCoV‐OC43‐derived peptide causes persistent mechanical hypersensitivity in female rats

We have shown repeatedly that MBP84–104 peptide produces mechanical allodynia with no effect on thermal sensitivity (Liu et al, 2012; Ko et al, 2016) likely due to its myelin‐dependent pronociceptive activity on myelinated A‐afferents, sparing unmyelinated heat‐sensitive nociceptors (Shubayev et al, 2016). Thus, to test the ability of the OC43‐derived peptide to regulate mechanical hypersensitivity characteristics to the homologous MBP84–104, we used a synthetic 20‐amino acid OC43653–673 peptide (OC43p, VSK ALAFA), and scrambled peptides OC43p‐SCR1 and OC43p‐S2 designed to mismatch the MBP84–104 amino acid sequence. Female rodents display robust mechanical hypersensitivity to intrasciatic MBP84–104 relative to males (Chernov et al, 2020). Female rats received a single bolus intrasciatic injection (Fig 2A) of OC43653–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 OC43653–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 MBP84–104 peptide sequences and PBS vehicle in the equivalent experimental designs (Liu et al, 2012; 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 pain‐like behavioral indices, although differences from the control animals were not statistically significant (Fig 2D). We concluded that OC43p, like MBP84–104, induced a robust and sustained pain mechanical hypersensitivity in female rats.

Figure 2

OC43p is a prospective determinant of pronociceptive activity

A schematic of the injections into the sciatic nerve followed by ipsilateral DRG and dorsal spinal cord tissue analysis.

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) ± standard deviation; two‐way analysis of variance (ANOVA) with Bonferroni post hoc test: **P ≤ 0.005; ***P ≤ 0.0005; and ****P ≤ 0.00005.

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.

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.

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 (|Log2(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.

Figure 3

Transcriptome changes induced by OC43p

Volcano plots of most significant DEGs in (A) DRG and (B) spinal cord. The size of each circle is proportional to Log2FC. Red and green colors indicate up‐ and downregulated DEGs, respectively, relative to thresholds (|Log2FC| > 0.58, P adj < 0.1, n = 3/group) displayed by dashed lines. Selected DEGs are labeled.

Figure 4

Transcriptome changes in the spinal cord

Principal component (PC) analysis of DEGs in OC43p (red) and OC43p‐SCR (blue) groups (n = 3/group).

Hierarchical clustering plot of 548 significant upregulated DEGs (Log2FC > 1, P adj < 0.1, n = 3/group). Heatmap color scheme corresponds to logarithms of variance stabilized counts.

Hierarchical clustering plot of X‐linked upregulated DEGs (Log2FC > 1, P adj < 0.1, n = 3/group).

Hierarchical clustering (Fig 4B) showed the most significant DEGs ranked by Log2FC between the OC43p and OC43p‐SCR groups in the spinal cord. A remarkable enrichment of transcripts encoding many voltage‐gated ion channels was observed in the OC43p group. Calcium (Cacna2d1, Cacnb4, Cacna1c, Cacna1d, Cacng4, Cacna1e, Cacna1a, and Cacna1b) and sodium (Scn7a, Scn3a, Scn1a, and Scn2a) voltage‐gated channels exhibited an increase. A set of potassium voltage‐gated channels was upregulated, including Kcnq3, Kcna3, Kcnj13, Kcnk9, Kcnma1, Kcnh7, Kcnj3, and Kcnh5. Transient receptor channels (Trpm3, Trps1, Trpc5, and Trpm7), nicotinic receptors (Chrnb3, Chrna7, Chrna5, and Chrnb2), muscarinic Chrm2 receptors, glutamate ionotropic NMDA‐type receptor Grin2A, glutamate ionotropic AMPA‐type receptors (Gria2 and Gria1), and GABA receptors (Gabrg3, Gabra2, Gabra3, and Gabrg2) were significantly upregulated.

Crucial innate immune system genes encoding pattern‐recognition receptors (PRPs), including Toll‐like receptors (TLRs) Tlr4, Tlr7, Tlr8, Tlr12, and Tlr13, exhibited upregulation. Interleukin receptors Il17, Il1, Il20, Il7, Il2, and Il6, and chemokine receptors Cx3cr and Ccr5 increased, although no increase in cytokine ligands was detected.

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

OC43p‐regulated signaling pathways were predicted using the Ingenuity Pathway Analysis (IPA) knowledgebase and DEGs with |Log2FC| > 1 and P adj < 0.1. Pathogen response‐specific pathways, FXR/RXR, LXR/RXR, T‐cell receptor signaling, and immune response‐specific 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 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.

Figure 5

Canonical pathways affected by OC43p

Canonical pathways (A, DRG and B, spinal cord) identified by IPA are ranked by P (‐log10 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.

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 channel Kcnma1, PI3K family member Pik3c2a, protein kinase C epsilon type (PKCε), Irs1, and Socs4 (Fig 6A).

Figure 6

Pronociceptive signaling induced by OC43p in rat spinal cord

Prolactin signaling pathway plotted in IPA. Red and green colors indicate up‐ and downregulated DEGs, respectively, relative to thresholds (|Log2FC | > 0.58, P adj < 0.1). Orange arrows indicate activation of signaling.

Gene ontology (GO) molecular function clusters are schematically shown on a two‐dimensional scaling plot. Cluster groups 1–11 are marked by colors and explained on panel 6C. Clusters of potentially pronociceptive molecular functions are circled.

Molecular function clusters are summarized on the best‐match average (BMA) distance heatmap. Heatmap colors correspond to a number of GO terms in each cluster.

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.

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).

Figure 7

OC43p affinity to the Na+/K+‐transporting ATPase complex

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 log10P – confidence score and peptide sequence coverage (%) were calculated in PEAKS Studio™ 8.5 (BSI).

Relationship of ATPase subunits identified in string‐db.org. Lines indicate evidence‐based interactions.

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 Log2(TPM).

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.

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é 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 al, 2000; Boziki et al, 2020; Gutiérrez‐Ortiz et al, 2020; Kemp et al, 2020; Kilinc et al, 2020; Koralnik & Tyler, 2020; Manji et al, 2020; Mao et al, 2020; Montalvan et al, 2020; Romoli et al, 2020; Sancho‐Saldaña et al, 2020; Troyer et al, 2020; Widyadharma et al, 2020; Zhao et al, 2020; Attal et al, 2021; Koike et al, 2021; McFarland et al, 2021).

As a proof of concept, the perspective synthetic peptide used in this study was derived based on strong sequence and structure homology to MBP84–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 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 infection‐related 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, 2015, 2018a).

The N‐terminal invariable 87‐(V/I)VHFFK‐92 motif of MBP84–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 MBP84–104 and OC43p peptides may be crucial to promoting pronociceptive and 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 inflammation‐induced 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 MBP‐induced 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 MBP84–104 peptide (Liu et al, 2012; Ko et al, 2016) consistent with the model of a myelin‐dependent pronociceptive activity on myelinated A‐afferents, sparing unmyelinated heat‐sensitive nociceptors (Shubayev et al, 2016).

The proallodynic MBP84–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 female‐specific 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.

Materials and Methods

Peptides

Peptides OC43p (VSKIVHFFKTFTTSTALAFA), OC43p‐SCR (VFIAHSVKFTKSFTLATTFA), and OC43p‐S2 (DNPVLHYFASTEKSN) were synthesized with > 95% purity, N‐terminal acetyl, and C‐terminal amide groups. Trifluoroacetic acid was removed after synthesis, and counterions were exchanged for acetates. Biotin‐tagged peptides were synthesized with N‐terminal biotin modifications. Peptides were dissolved in sterile PBS (vehicle). Key reagents and resources are described in Table EV1.

Antibodies

Antibodies used for protein immunoblotting were as follows: anti‐ Na+/K+‐transporting ATPase α‐1 antibody, clone C464.6 ZooMAb® mouse monoclonal IgG (Millipore Sigma, Cat. no. ZMS1029, used at 1:5,000 dilution); rabbit polyclonal anti‐Na+/K+‐transporting ATPase α‐2 antibody (Millipore Sigma, Cat. no. 07‐674, at 1:2,500 dilution); cross‐adsorbed donkey anti‐rabbit, horseradish peroxidase (HRP)‐conjugate (Thermo Fisher, Cat. no. 0031458, at 1:5,000 dilution); and goat anti‐mouse IgG (H + L) HRP‐conjugate (Bio‐Rad, Cat. no. 1706516, at 1:5,000 dilution).

Amino acid homology search

Amino acid sequence homology search was conducted using position‐specific iterative PSI‐BLAST tool (www.ncbi.nlm.nih.gov/blast) (Altschul et al, 1997; Altschul & Koonin, 1998) within NCBI (www.ncbi.nlm.nih.gov), VIPR (www.viprbrc.org), and GISAID SARS‐CoV‐2 mutant variants (www.gisaid.org) depositories.

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 33‐gauge 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 mid‐hind 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 pain‐like 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 high‐sensitivity 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 dual‐index 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).

Figure EV1

Steps of RNA‐seq analysis workflow

Unprocessed sequence data generated by NovaSeq 6000 were provided in a fastq format. Purple rectangles indicate software used. Blue rectangles indicate analysis results and intermediate data. Green rectangles show reference data.

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. Log2FC was calculated using the Wald test. The adjusted (shrunken) Log2FC 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 removeBatchEffect (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ämer 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.

Affinity capture assay

Peptide‐bound beads were prepared by incubating biotin‐tagged OC43p and OC43p‐SCR peptides (9 nmole) with 100 µl Dynabeads MyOne Streptavidin T1 (10 mg/ml) for 2 h at 25°C in 500 µl of TBST (20 mM Tris, 150 mM NaCl, 0.1% Tween‐20, and pH 7.4). Beads were washed using a magnetic separation rack 6 times with 750 µl TBST to remove unbound peptides, and resuspended in 100 µl of affinity capture buffer (TBST supplemented with 1 mM CaCl2, 1 mM MgCl2, 1 mM Na3VO4, and EDTA‐free protease inhibitors (Roche)).

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,4‐dithiothreitol. 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™ Lumos Tribrid (Thermo Fisher Scientific). Proteomics data were analyzed using PEAKS Studio™ 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.

Author contributions

Veronica I Shubayev: Conceptualization; Resources; Funding acquisition; Writing—review & editing. Jennifer Dolkas: Investigation; Methodology. Glaucilene Ferreira Catroli: Visualization; Methodology. Andrei V Chernov: Conceptualization; Data curation; Software; Formal analysis; Supervision; Validation; Investigation; Visualization; Methodology; Writing ‐ original draft; Project administration; Writing—review & editing.

In addition to the CRediT author contributions listed above, the contributions in detail are:

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.

Table EV1

Click here for additional data file.

Dataset EV1

Click here for additional data file.