A total of 25,788 genes were detected by RNA-seq in the proximal nerve stumps, including 2,553 differentially expressed genes (DEGs, Supplementary Table 1) that conformed to the significance criteria (adjusted p-values (P_adj) < 0.1) and the expression fold change (FC) filter (log₂FC > 1 or log₂FC < –1, Figure 1A). In summary, more unique, highly expressed up-regulated DEGs (n = 232) were identified in males relative to females as presented on MA and scatter plots (Figures 1B,C). In addition, 65 and 159 DEGs were downregulated in female and male nerves, respectively.

Heatmaps (Figure 2A) and volcano plots (Figure 2B) present the identity of the most significant protein-coding DEGs in each sex. Principal component analysis (PCA) identified their contribution to variance in axotomized samples from males and females based on normalized count ranking. Sample-to-sample comparisons using PCA demonstrated that female and male samples aggregated as distinct clusters in the PC1/PC2 dimensions (Figure 2C).

The most significant DEGs that contributed to gene expression variance in females encoded a major acute-phase serum amyloid A-1 protein (Saa1), leucine-rich repeat-containing protein 15 (Lrcc15), neuronal pentraxin-2 (Nptx2), heterodimeric neurotropic cytokine receptor-like factor 1 (Crlf1), cardiotrophin-like cytokine factor 1 (Clcf1), and sonic hedgehog protein (Shh). In males, DEGs encoding the potassium channel subfamily T member 1 (Kcnt1), apoptosis-promoting dendrin (Ddn), inactive serine/threonine-protein kinase (Plk5), Toll-like receptor 9, and myelin-associated oligodendrocyte basic protein (Mobp) were the most significant drivers of variance. These and other DEGs identified by PCA (Figures 2C,D) could present potentially important sex-specific markers of early PNS injury response and nerve regeneration.

Many sexually dimorphic and monomorphic DEGs were recognized as important components of cell signaling pathways in the nervous and immune systems. To gain further mechanistic insight into their role in nerve regeneration, gene ontology (GO) enrichment analysis was conducted. Molecular-level biochemical characteristics of DEGs were used to identify the most relevant GO clusters of molecular function (Figures 3A,B). GO clusters that included sexually dimorphic DEGs relevant to nerve injury response are illustrated in Figure 3C and detailed below.

Signaling pathways, upstream regulation, and interactive networks were predicted by Ingenuity Pathway Analysis (IPA) software using a highly stringent subset of DEGs (P_adj < 0.001 and log₂FC > 1 or log₂FC < –1). Significant canonical pathways were further ranked based on Fisher’s exact test (p < 0.05). Positive or negative z-scores determined pathways’ activation and deactivation.

Axotomy-induced significant changes in cellular signaling processes in both sexes (Figure 4), including pathways with female- or male-specificity. Accordingly, GDNF Family Ligand Receptor Interactions, Inhibition of Matrix Metalloproteinases, Synaptogenesis Signaling Pathway, LXR/RXR Activation, PD1/PDL1 cancer immunotherapy pathway, and Coronavirus Pathogenesis Pathway demonstrated activation in both sexes. PTEN Signaling, CHK Proteins in Cell Cycle Checkpoint Control, Cell Cycle G2/M DNA Damage Checkpoint Regulation, and Endocannabinoid Neuronal Synapse pathways were activated in males. In females, Senescence Pathway and Pancreatic Adenocarcinoma Signaling were activated. Many pathways exhibited downregulation in both sexes, including Neuroinflammation Signaling, Pattern Recognition Receptors, Phagosome Formation, Pyroptosis Signaling, Leukocyte Extravasation Signaling, HOTAIR Regulatory Pathway, G Protein-Coupled Receptor Signaling, and Natural Killer Cell Signaling. In addition, males specifically reduced PI3K/AKT Signaling, Kinetochore Metaphase Signaling, IL-8 Signaling, and TREM1 Signaling pathways. In both sexes, regenerating nerves displayed activation of Macrophages in Rheumatoid Arthritis, Axonal Guidance, Thyroid Hormone Metabolism I via Deiodination, Atherosclerosis Signaling, and Agranulocyte/Granulocyte Adhesion and Diapedesis canonical pathways.

Predictive analysis of upstream regulators was conducted in IPA to identify regulatory molecules, including transcription factors (TFs), any gene or small molecule that, with high probability, could affect the expression of their target DEGs. The activation or inhibition efficiency of each upstream regulator on target DEGs could be defined by positive or negative z-scores, respectively, as illustrated in Figure 5 and summarized in Supplementary Table 2. Worth noting that upstream regulator z-scores do not directly reflect the expression of the regulator itself (Krämer et al., 2014), as different organs or tissues can contribute regulatory molecules. IPA identified many sexually dimorphic upstream regulators, including the brain-derived neurotrophic factor (Bdnf), a known stimulator of the nerve axon growth (Zhou and Shine, 2003). Bdnf was predicted to strongly affect gene expression in males but not females. Bdnf level was low in the proximal nerve stumps in both sexes. However, Bdnf mRNA was highly abundant in DRG and demonstrated a male-dominant increase post-axotomy (Chernov and Shubayev, 2021). Other male-specific upstream positive regulators included Sox2, Aire, Egr2, Cdkn1a, Map2k1/2, and other molecules. Ikzf3, Mek, Erk, and Akt1 were predicted as prospective positive regulators in females. Remarkably, IRFs, STATs, and cytokines were identified among significant negative regulators in both sexes.

Protein kinases and their regulatory molecules were identified as significant sexually dimorphic upstream regulators (Figure 6A and Supplementary Table 2). Male-specific upregulation of Kndc1 (also known as Very-KIND) may modulate Hras, Myc, Neu1, Ras, and Map2 (Figure 6B) during neuronal growth (Huang et al., 2007). Highly expressed Cdk5 was moderately activated by axotomy in male nerves. Cdk5 activator molecule that also controls cytoplasmic or nuclear localization of the Cdk5 complex demonstrated male-dominant upregulation. Co-upregulation of Cdk5-interacting and target molecules in males was significantly stronger than in females (Figure 6C). While female nerves upregulated fewer protein-coding DEGs than males, two female-specific interactive networks were determined using the STRING database (Szklarczyk et al., 2021). The first proposed network potentially conveys Gdnf signaling via the neuron-specific adaptor Shc3, the Gdnf receptor α1 (Gfra1), and neural cell adhesion proteins Chl1 and Cadm1 to regulate cell-cell adhesion and neuronal plasticity (Figure 6D). Another network is centered around the regulatory component of the Cyclin D1 (Ccnd1)-Cdk4 complex (Figure 6E). It is known to respond to Timp1 that activates cyclic AMP-induced regenerative program gene post-axotomy (Liu et al., 2015) partly via Cd63/Pi3k/Akt signaling (Rossi et al., 2015).

Regulatory signaling events reflect distinct molecular programs involved in regenerating nerves. Programs exhibiting sexually dimorphic regulation are as follows:

After nerve injury, Schwann cells assume a pro-regenerative function (Jessen and Mirsky, 2016, 2021; Merrell and Stanger, 2016; Milichko and Dyachuk, 2020), exhibiting an exceptional phenotypic plasticity (Jessen and Arthur-Farraj, 2019; Stierli et al., 2019; Nocera and Jacob, 2020) that supports their de-differentiation, proliferation, and re-differentiation into myelinating and non-myelinating phenotype to facilitate repair of the respective axons. Schwann cell reprogramming is controlled by a set of TFs (Balakrishnan et al., 2021), including Myt1l, Pou3f1/Oct6, Myrf, Olig1/2, Jun, Sox-, Hox-, and Fox-family members TFs upregulated at 24 h post-axotomy predominantly in males. Transcripts of TF genes, such as Fos, Smad3, Creb1, Rac2, Klf6, and Atf3, that function as RAGs (reviewed in [Van Kesteren et al., 2011; Ma and Willis, 2015)] in the corresponding DRGs after PNS injury (Chernov and Shubayev, 2021) were detected in proximal nerve segments. In addition to transcriptional activity in cells of the PNS, this accumulation may reflect an anterograde transport from DRG in coordinated axonal regeneration processes post-axotomy. Single-cell RNA-seq and spatial transcriptomics analysis could provide crucial information on the identity of cell lineages activating these sexually dimorphic TFs.

The GDNF family of ligands (GDNF, neurturin, artemin, and persephin), known survival factors for neurons, bind to GFRA receptors, trigger the phosphorylation of RET tyrosine kinase receptor (Durbec et al., 1996) that is known to regulate Pi3k/Akt and Plcγ/Ip3r dependent Ca²⁺ signaling during neurogenesis (Fukuda et al., 2002; Lundgren et al., 2012) and, sometimes, hyperalgesia (Bogen et al., 2008). Furthermore, GDNF/Gfra1 could interact with neural cell adhesion molecules such as NCAM to induce an axonal expansion (Nielsen et al., 2009). GDNF family of ligands can act in a synergetic manner with other growth factors, including transforming growth factor-β (Tgf-β) and sonic hedgehog (Shh) [reviewed in Jessen and Mirsky (2021)]. GDNF signaling can stimulate the migration of neuronal precursors and Schwann cells (Iwase et al., 2005) to the injury site in both sexes.

The outcome of nerve regeneration after sciatic nerve injury is partly dependent on the integrity of ECM-rich Schwann cell basal lamina and endoneurial tubes (Fawcett and Keynes, 1990). Proteolytic cleavage mediated by ECM proteinases and their endogenous inhibitors is indispensable for morphological remodeling of regenerating axons (Rivera et al., 2010; Fujioka et al., 2012; Bajor and Kaczmarek, 2013), after sciatic nerve injury (Liu et al., 2010). In addition, MMPs control the activity of critical growth factors, their ligands, and receptors in peripheral nerve, including IGF and neuregulin/ErbB (Chattopadhyay and Shubayev, 2009). Mmp17 was upregulated in both sexes and can convey inflammatory responses due to the TNFα-converting activity (English et al., 2000).

Intriguingly, a group of MMPs co-localized on the mouse chromosome 9 (Mmp3, Mmp8, Mmp12, Mmp13, and Mmp27) have decreased in both sexes. Because MMPs in the orthologous human chromosome 11q22.3 region could carry different histone modification marks (Chernov et al., 2010), epigenetic control of this set of metzincins may be involved. Consistent with our previous report, MMPs, and many disintegrin ADAM/ADAM-TS families of the ECM-remodeling metzincins demonstrate regulated expression in injured nerves (Chernov et al., 2015).

Notably, protease expression levels were high in intact and axotomized nerves. Processes of post-translational protease activation and binding to intrinsic proteinase inhibitors control potentially cytotoxic proteolytic activities. Aberrant cleavage of mediators of nociceptive signaling during neurogenesis potentially causes neuropathic pain (Remacle et al., 2015). Previously we identified Timps as key genes in the normal and damaged nerves (Kim et al., 2012; Nishihara et al., 2015) and post-axotomy DRG (Liu et al., 2015; Chernov and Shubayev, 2021). Timp1 binds and inhibits all MMPs, except Mmp14; Timp2 inhibits Mmp14 (Brew and Nagase, 2010). In the proximal nerve, Timp1 demonstrated strong upregulation in both sexes. Other upregulated endogenous inhibitors of proteinases Wfdc5, Fetuin B, Spink1, Itih3, and Serpineb1b could serve a not well-understood role in nerve damage response.

An intact connection to DRG allows for dynamic axonal transport of protein and mRNA content between neuronal soma and the proximal segment of the axotomized nerve (Cavalli et al., 2005). The anterogradely transported mRNAs from DRG are directly used for the local protein synthesis at the site of the injury (Sahoo et al., 2018) and complemented by the transcriptional repertoire of the resident non-neuronal cells in the damaged nerve (e.g., Schwann cells, macrophages) that establish a microenvironment generally permissive for axon regeneration (Avraham et al., 2021).

Kinesins facilitate anterograde axonal transport of neurotransmitter receptors, mitochondria [reviewed in (Fan and Lai, 2021)], and also mRNAs [reviewed in (Kanai et al., 2004; Turner-Bridger et al., 2020)]. Remarkably, most kinesin mRNAs were depleted in proximal nerve segments at 24 h post-axotomy, except for members of the conventional kinesin-1 family Kif5a/b/c increased in males.

Retrograde transport of transcription factors from the site of nerve injury to DRG has been shown to facilitate the neuronal survival (Cox et al., 2008) and axon regeneration (Ito and Enomoto, 2016). In the corresponding DRGs of the same animal cohort, at 24 h post-axotomy, we observed male-specific reduction of mRNAs encoding calcium, sodium, and potassium ion channels as well as ionotropic AMPA, NMDA, and GABA receptors (Chernov and Shubayev, 2021) that correlated to their increase in proximal segments in the present study and presumably, anterograde axonal transport. We conclude that axonal transport of mRNAs contributes to sex-specific protein synthesis and downstream signaling programs both in DRG neurons and regenerating axons at the early stages of response to sensory nerve damage.

Sex chromosome-encoded genes were activated in regenerating sciatic nerves. Several X-linked genes were differentially regulated in male and female nerves. Among X-linked genes induced in nerves of both sexes were Timp1, known to promote regeneration (Kim et al., 2012; Liu et al., 2015; Nishihara et al., 2015), and Map7d2, involved in kinesin-mediated cargo transport into the axon (Pan et al., 2019). It is worth noting that genes associated with the female-specific epigenetic X-chromosome inactivation (Xi) process were differentially expressed in DRGs within 24 h after post-axotomy in the same animal cohort (Chernov and Shubayev, 2021).

The noted remodeling represents baseline sexual dimorphism in uninjured peripheral nerves and an immediate regenerative response to nerve injury. The functional significance of the reported sexual dimorphism remains to be determined. It is conceivable that, at least partially, sex differences were equalized over the time-course of nerve injury. Single-cell RNA-seq and spatial transcriptomics analysis could provide crucial information on the identity of cell lineages activating these sexually dimorphic coding and non-coding RNAs.

Detailed descriptions of reagents and resources are included in Supplementary Table 3.

Female and male C57BL6/J mice (6–8 weeks old, Jackson Labs, n = 24) were housed in a temperature-controlled room (22°C) on a 12 h light/dark cycle with ad libitum access to food and water. The mice were randomly assigned to axotomy and sham groups by sex (n = 6/group). Under isoflurane anesthesia, the left sciatic nerve was exposed at the mid-thigh level, followed by a complete transection using sterile microsurgery scissors. In the sham group, anatomically equivalent nerves were subjected to exposure without transection. The muscle was sutured, and the skin stapled. At 24 h after surgery, proximal segments of the sciatic nerve, corresponding to previously investigated lumbar (L)4 and L5 DRG tissues (Chernov and Shubayev, 2021), were collected for RNA isolation. 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.

All surgical and tissue harvesting instruments were sterilized and repeatedly treated with RNase Away reagent followed by RNase-free water rinse. Tissues were immediately submerged in 500 μl RNAlater Stabilization Solution to preserve RNA integrity, placed at 4°C overnight, then transferred for storage at –20°C. All sample groups were processed in parallel to minimize batch effects. Proximal stumps of sciatic nerve tissues were pooled from 2 mice per group for RNA purification.

Nerve tissues were transferred in Trizol solution (Invitrogen) and disrupted by mechanical homogenization. Total RNAs were purified using RNeasy RNA purification reagents. RNA concentrations and quality were determined using Nanodrop absorbance ratios at 260/280 nm and 260/230 nm. RNA integrity indices were determined using the Agilent Bioanalyzer Nano RNA chip. 500 ng of total RNA samples (n = 6 mice/group, 2 mice/sample (pooled), 3 sample/group) with RIN ≥ 7.0 were used for RNA-seq.

mRNA libraries were generated following the TruSeq Stranded mRNA library preparation protocol (Illumina). In brief, the 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 magnetic 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 Agilent Bioanalyzer DNA high sensitivity chip and normalized. RNA-seq was performed using the paired-end 100 cycle program on the NovaSeq 6000 system at the Genomics High Throughput Facility (University of California Irvine). Base calls were recorded and converted to FASTQ files containing sequencing reads and the corresponding quality scores using Illumina software. Sequencing was conducted until we acquired at least 50 million paired-end reads per sample.

FASTQ files were uploaded to the Amazon S3 server and processed using Elastic Compute Cloud (EC2) [Amazon Web Services (AWS)] running Ubuntu Server 20.04 LTS (64-bit ARM). Data analysis steps are summarized in Supplementary Figure 1 and Supplementary Table 3. 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 and mapping were performed using Salmon (Patro et al., 2017) and the Gencode M29 mouse genome. Mapping coverage was estimated in MultiQC (Ewels et al., 2016). Transcript- to gene-level quantifications were done using Tximeta (Love et al., 2020). Gene count matrices were processed in the DESeq2 (Love et al., 2014). Log₂FC values were calculated using the Wald test and adjusted using the adaptive t-prior apeglm method (Love et al., 2014). Significant DEGs were identified by P_adj values below a false discovery rate (FDR) cutoff (P_adj < 0.1) (Supplementary Table 1). P_adj < 0.1 were used in downstream analyses unless otherwise noted. Batch effects were controlled using the removeBatchEffect tool (Ritchie et al., 2015) and RUVseq (Risso et al., 2014) functions. DEGs were visualized using ggVennDiagramm, PCAtools, ComplexHeatmap, and EnhancedVolcano R packages.

Prediction and biological interpretation of the regulated signaling pathways and upstream regulators in female and male animals was made in IPA using the default parameters (Krämer et al., 2014). Significance criteria of P_adj < 0.1, log₂FC > 1 or log₂FC < -1 were applied. Regulation directionality was estimated in IPA based on z-scores. Gene ontology terms were determined using Geneontology¹ and the STRING database using the default parameters.