Next-generation RNA sequencing elucidates transcriptomic signatures of pathophysiologic nerve regeneration

The cellular and molecular underpinnings of Wallerian degeneration have been robustly explored in laboratory models of successful nerve regeneration. In contrast, there is limited interrogation of failed regeneration, which is the challenge facing clinical practice. Specifically, we lack insight on the pathophysiologic mechanisms that lead to the formation of neuromas-in-continuity (NIC). To address this knowledge gap, we have developed and validated a novel basic science model of rapid-stretch nerve injury, which provides a biofidelic injury with NIC development and incomplete neurologic recovery. In this study, we applied next-generation RNA sequencing to elucidate the temporal transcriptional landscape of pathophysiologic nerve regeneration. To corroborate genetic analysis, nerves were subject to immunofluorescent staining for transcripts representative of the prominent biological pathways identified. Pathophysiologic nerve regeneration produces substantially altered genetic profiles both temporally and in the mature neuroma microenvironment, in contrast to the coordinated genetic signatures of Wallerian degeneration and successful regeneration. To our knowledge, this study presents as the first transcriptional study of NIC pathophysiology and has identified cellular death, fibrosis, neurodegeneration, metabolism, and unresolved inflammatory signatures that diverge from pathways elaborated by traditional models of successful nerve regeneration.

signatures to guide prognostic biomarker discovery or identify mechanistic targets to prevent neuroma formation and imminent neurologic consequences, such as chronic pain and limited functional recovery.
To address the translational gap from basic science models to clinical injuries, we have created and validated a laboratory model of rapid-stretch nerve injury, which is the predominant mechanism associated with traumatic nerve injury and neuroma formation in the clinical setting 2 . Biomechanical 31 , histologic 32 , and behavioral 33 outcomes faithfully recapitulate the neuroma phenotype and neurologic deficits found in the clinical population. Expanded comparison to traditional models of crush, transection, and delayed transection repair has delineated recovery profiles specific to the rapid-stretch mechanism of injury and associated neuroma pathophysiology 34 . To further this understanding at the molecular level, we sought to elucidate genetic signatures associated with NIC formation and identify critical pathways associated with the transition from favorable regeneration towards pathologic remodeling. Here, we present insights into the aberrant regenerative program through a temporal examination of genetic activity and application to downstream biological and functional pathways.

Results
Survey of genomic expression profiles. After rupture injury, there was a dramatic increase in genetic activity, with upregulated genes peaking at day 2 (D2) and declining linearly over time; this pattern was mirrored by downregulated genetic activity patterns (Fig. 1a). Despite the decrease in differentially expressed genes (DEGs), at the final timepoint (D48) there remained a large number of DEGs versus baseline. In contrast, although sham injury also demonstrated marked increase in DEGs in the acute phase after surgery, DEGs returned to baseline expression patterns.
Global changes in genetic expression were further delineated through principal component analysis (PCA), which aims to deconstruct multidimensional differences in gene expression to elucidate similarities in genetic activity profiles (Fig. 1b). Collectively, samples for each timepoint for the different injury severities demonstrated reliable clustering patterns, except that the sham injury was indistinct from control by D48. Interestingly, sham and rupture demonstrated relatively tight clustering for D0 and D2, then decoupled into distinct profiles for the remaining timepoints, indicative of unique genetic programs starting by D7.
Cluster analysis through Euclidian distance matrix mapping was performed to identify gene clusters demonstrating similar temporal activity profiles (Fig. 2a). Five clusters were distinguished that broadly bifurcated into 2 temporal activation patterns: (1) immediate up-or downregulation (Clusters 1 and 4) and (2) delayed up-or downregulation (Clusters 2, 3, 5). Gene clusters were then subject to Ingenuity Pathway Analysis (IPA) to identify metabolic and canonical signature pathways (Fig. 2b).
The immediate downregulated genes (Cluster 1; decreased D0 with partial return to control expression at D48) were grossly associated with neuronal signaling pathways, such as Axonal Guidance Signaling and Netrin Signaling (Fig. 2b). Immediate upregulated genes (Cluster 4; increased D0 with partial decrease to control expression at D48) demonstrated immunologic signatures, with the top statistically significant pathways implicating cytokine signaling (IL-10, IL-6), T cell activity (Th1 and Th2 Activation Pathway), and the innate immune response (Granulocyte Adhesion and Diapedesis, HMGB1 Signaling). Delayed downregulation genes (Cluster 2; decreased D2 with partial return to control expression at D48) aligned with pathways implicated in metabolic regulation. Cluster 3 (increased D2 with partial return to control expression at D48) genes were associated with cell-cycle pathways and a second wave of the immune response. Finally, Cluster 5 (decreased D2 and no substantial return to control expression at D48) genes were related to pathways primarily involved in neuronal processes, such as Synaptogenesis and Glutamate Receptor Signaling, although the most significant path was associated with metabolic dysregulation, i.e., Superpathway of Cholesterol Biosynthesis. Importantly, in contrast to the other clusters, which chiefly resolved towards the control activity state, Cluster 5 genes had persistently decreased expression and involved axonal pathways in IPA canonical pathways analysis.
Pathway and functional analysis of differentially expressed genes at distinct timepoints. To develop a framework for understanding pathways involved in neuroma formation, DEGs were subject to Gene Set Enrichment Analysis (GSEA), IPA canonical pathways analysis, and IPA diseases and functions analysis. Utilization of two different analysis platforms, each with unique methods of attributing genetic activity to pathways, permits a corroborative framework for results.
GSEA. DEGs between the rupture and control groups were subjected to GSEA analysis, broadly demonstrating early and prolonged upregulation of inflammatory gene sets such as Complement, IFNγ Response, and TNFα Signaling via NFκβ in the rupture group. A modest increase in angiogenesis occurred at early timepoints, while a notable increase in genes associated with myogenesis and oxidative phosphorylation occurred from D7 to D48. Consistent upregulation of Apoptosis and Epithelial-Mesenchymal Transition was also present. In contrast, pathways associated with metabolism were downregulated across all experimental timepoints (Fig. 3a).
Ingenuity canonical pathways analysis of DEGs. IPA canonical pathway analysis applies the experimentally derived genetic expression profiles to well-defined canonical cell signaling and metabolic pathways to identify significant associations and predicted activation states (Fig. 3b,c). In the acute phase after injury (D0 and D2), the majority of upregulated canonical pathways were significantly associated with the inflammatory response, such as TREM1 Signaling and IL-17 Signaling. In comparison, D7 was hallmarked by increased paths associated with the cell cycle and paths associated with neuronal activity, such as CREB Signaling in Neurons. Integration of neuronal, inflammatory, and cell-cycle pathways presented for significant upregulation for D14. Finally, persistent inflammatory response on D48 was highlighted by significant upregulation of pathways such as Neuroinflammation Signaling (Fig. 3b). Downregulated pathways were predominately associated with metabolism www.nature.com/scientificreports/ (Fig. 3c). Inhibition of extracellular matrix degradation was also implicated, with significant downregulation of the Inhibition of Matrix Metalloproteases for multiple timepoints (D0, D7, D14), while axonal regeneration signatures, e.g. Synaptogenesis Signaling Pathway, were decreased at D48 (Fig. 3c).
Diseases and functions analysis of DEGs. In addition to IPA pathway analysis, diseases and functions analysis was utilized to provide broad biological and pathophysiologic context, with the goal of identifying relative  www.nature.com/scientificreports/ influence or change over time from injury (Fig. 3d). Broad categories are defined by IPA, with further stratification into sub-annotations of specific diseases and functions. In our data set, six categories were identified to have robust longitudinal representation: Cell Death, Inflammatory Response, Cellular Growth and Proliferation, Neurological Disease, Skeletal and Muscular System Development and Function ("Muscle"), and Fibrosis. To deconstruct these specific categories and assimilate into tangible patterns, a part-to-whole approach was employed; the counts of sub-annotations for each broad category, out of sub-annotations for all categories, was calculated per day and represented as a percent. Inflammatory Response predominated the transcriptional landscape across all timepoints until tapering at D48 (Fig. 3d). Representation of Neurologic Disease increased over time, with the highest proportion of representation on D48. Cellular Growth and Proliferation demonstrated an inverse relationship compared with Neurologic Disease, with the highest proportions in the acute and early timepoints (D0-D7) and decreases for the later timepoints. Cell Death was parabolic, with a higher proportion at D0 and then decrease until D7, and then increase for the later timepoints with peak representation by D48.
D48 overlapping pathways and critical genes. Because many genes overlap or coordinate in function, we investigated the frequency of overlap of genes considered critical in IPA to provide an assessment of which genes are involved in the interplay of the predominant pathways in our dataset. Described simply, we sought to identify specific genes suspected of participating in the mature neuroma pathology. Collectively, the categories of Inflammatory Response, Neurologic Disease, and Cell Death displayed the highest overlap in gene products at the terminal endpoint (D48) where we have previously demonstrated mature neuroma formation. Therefore, we further examined the degree to which these gene products were shared with multiple categories or were unique to the category in the mature pathophysiologic neuroma (Fig. 3e) Interaction model to identify genes distinctly significant to rupture injury. To define the cardinal genes associated with neuroma formation, we employed an interaction model to identify the statistically significant genes distinct or influential to the rupture-injury state as compared with the sham injury. Sham injury is an acute inflammatory response associated with the surgical procedure but not neurologic injury, as clearly demonstrated by resolution of gene expression between D7 and D48 and subjection of DEGs to IPA diseases and functions analysis ( Supplementary Fig. S1). This approach permits specific identification of pathways associated with the neuroma microenvironment to identify distinctions that might be otherwise too nuanced to appreciate with the overwhelming inflammatory response. The significant genes determined by the interaction model were then subject to IPA network analysis, which constructs de-novo networks based on the interconnectedness of genes in the dataset and the premise that highly interconnected genes reflect germane biological networks. These networks were then independently overlaid with IPA canonical pathways ( Supplementary Fig. S2) and IPA diseases and functions ( Fig. 4) to provide context of how the de novo-constructed genetic networks are related to traditionally established biologic and molecular pathways. Identified pathways from both overlays predominately aligned with 5 motifs: Vascularity, Neuronal, Fibrosis, Cell Death, and Inflammation. Overlaying both IPA canonical pathways and IPA diseases and functions to the network analysis permits independent corroboration of the de novo-constructed networks to downstream biological function.
We focused on IPA diseases and functions analysis as this analysis employs a predication algorithm to determine activation/inhibition of downstream biological function states based on relative contribution of DEGs in our dataset to the genes in the established disease/function pathways (Fig. 4). Only pathways with a p-value < 1 × 10 -5 are presented and annotated to the identified motifs. Furthermore, we selected the transition between D0 and D2 ( Fig. 4a,b) and D2 to D7 (Fig. 4c,d) , when the pathways diverged on PCA, thereby focusing on the early events that discriminated outcomes. The top network for D0 vs. D2 (Fig. 4a) largely predicted activation of pathways associated with Vascularity, Neuronal, Fibrosis, and Inflammation. Cell Death paradoxically predicted global activation yet specific inhibition of Neuronal Cell Death. Critical hub genes with significant upregulation included CCL11, GDNF, VEGFA, INHBA, and GFRA1. Significantly downregulated hub genes included FGF4, CNMD, SLCO1C1, and FGFR2. Critical hub genes with predicted inhibition included P38 MAPK, FBXO32, and MAL while predicted activation included JNK, P2RY6, HS6ST1, CHR1, and Histone H3. The second highest network for D0 vs. D2 (Fig. 4b) presented activation of cell death and inflammation pathways and predicted inhibition of fibrosis pathways. Although neuronal pathways presented conflicting activation/inhibition patterns, pathways were largely associated with neuronal inhibition, such as activation of a negative pathway. Critical hub genes with significant upregulation included SLC6A1, SCN4A, WNT16, MOBP, TPPP, BFSP2, and CYP4F12, while significantly downregulated hub genes included GDNF, UTP14C, and ALDH1A3. Critical hub genes with predicted inhibition included SOX2, NFASC, MAG, RPTOR, and CTNNB1, while only MYH6 demonstrated predicted activation. For the top network of D2 vs. D7 (Fig. 4c), we observed reversal of the patterns described on D0 vs. D2; pathways aligning to neuronal, fibrotic, and cell death motifs largely predicted inhibition, while there was limited significant annotation to vascular or inflammatory paths. Critical hub genes with significant upregulation included CAV3, DES, and SGCA. Significantly downregulated hub genes included KCNJ10, PMP22, Immunofluorescence for protein corroboration. Structural changes. Longitudinal and multiple independent genetic analyses identified repeating motifs of Neurologic Disease, Inflammation, Fibrosis, Muscle, Cell Proliferation, and Cell Death. Because these categories were globally represented in GSEA and other IPA analyses, key molecules were assessed through immunofluorescent histology for protein corroboration. Molecules that are often used clinically stood out among the implicated genes in the pathways and thus were used to provide a translational comparison (Fig. 5). The microarchitectural features of Neurologic Disease demonstrated progressive decline in axons (anti-NF200) and myelin (anti-MBP). Thus, inverse correlation between protein expression and representation of neurologic disease in the diseases and functions analysis in Fig. 2 was observed. MYH4 was frequently associated with the muscle pathways yet is also present in the myotubules of many cellular species; thus, the identification of muscle-specific pathways is at risk of misattribution. MYH4 demonstrated apparent overlap with NF200 at D0, D2, and D14 after injury as well as in control samples. Limited expression was observed from D7 to D48 within the nerve. Thus, although the increase in transcript expression could be attributed to small quantities of muscle surrounding the neuroma, which cannot be resected in a gross manner without risking compromising the neuroma itself, immunofluorescent staining suggests MYH4 is common within other cellular origins. Genetic patterns of fibrosis were identified throughout multiple analyses, and aberrant fibrotic deposition is a hallmark of clinical neuroma pathology. In particular, collagen is a well-known component of the endoneurial tubules; thus, we sought to corroborate multiple types in the pathologic progression. Col1a is the predominant collagen of the collagenous fiber type while Col3a is associated with reticular fiber types. By D2, we observed alterations in collagen expression, hallmarked by disorganized matrix deposition that persists longitudinally throughout our timepoints (Fig. 5).
Immunologic changes. Inflammatory pathways ubiquitously presented across analysis platforms with significant upregulation across timepoints (Fig. 6). CD11b is a pan-granulocyte marker with robust expression after injury-corroborative of the persistent activation of inflammatory paths in the genetic analysis. Pockets of these cells were observed on D2 and D7 that clustered in locations on the exterior of the nerve with some internal infiltration. Surprisingly, robust expression was still observed on D14; although topographically, these cells may exist interior of the formed neuroma, they remained localized to the exterior boundary and were segregated from the apparent dead ends of injured nerve ( Supplementary Fig. S3). CD68 presents as a marker for activated M1-like (pro-inflammatory) macrophages, which we observed transition into the nerve by D2, increase in size by D7, develop a "foamy" phenotype at D14, and then decrease in size yet remain in the neuromatous tissue through D48. Interestingly, at this terminal endpoint, CD68 staining appeared fragmented. CD4 is a marker for T-helper cells and dendritic cells in the mouse and thus suggests crosstalk to adaptive immunity. Interestingly, we observed high expression on D2 surrounding the CD11b pockets. Although mRNA expression levels increased by D7 and D14, this did not translate to a functional protein per immunofluorescent staining. At D48, we visualized clusters of CD4 within the neuroma, yet morphologically distinct cells from D2, thus potentially representing a different cell type at this time. DAPI expression appeared to correlate with increased immune cell infiltration. Interestingly, we observed denuclearized DAPI deposited along the exterior of the nerve and dramatic intensity. By D2, robust cellularity was observed on the exterior of the nerve, yet the DAPI largely delineated distinct nuclei when compared with D0. At D7, we observed the infiltration of this cellularity, which persisted through D48.  www.nature.com/scientificreports/  www.nature.com/scientificreports/

Discussion
The regenerative capacity of peripheral nerves is well defined in basic science models, yet failed neurologic recovery continues to plague clinical outcomes. This translational gap may in part be explained by the predominance of animal models investigating mechanisms of successful regeneration. Probing the mechanisms underpinning pathologic nerve regeneration may reveal previously unidentified and underappreciated pathways governing failed regeneration and thereby present as novel targets for therapeutic intervention. The current investigation sought to further interrogate this pathophysiologic focus by applying RNAseq to identify significantly implicated transcripts and annotate their functional capacity to downstream signaling pathways. To our knowledge, this study is the first transcriptional study of NIC pathophysiology. Here, we have identified inflammatory, cellular death, fibrotic, and neurodegeneration signatures that diverge from pathways identified by traditional models of successful nerve regeneration. These insights afford numerous downstream investigations for future interrogation, with the potential to guide prognostic biomarker heuristics or targeted therapeutic intervention.

Inflammation.
Recently, molecular-level assessments of nerve regeneration and WD have ubiquitously revealed substantial enrichment of pathways associated with inflammation [24][25][26][27] . These results corroborate the cellular profiling of seminal WD studies [35][36][37] ; however, the transcriptional profiles depict attenuation of the inflammatory response, as genes specific to inflammatory pathways peak from 4 to 7 days after injury and decline thereafter through day 14, which coincides with onset of signaling pathways associated with axonal guidance, nervous system signaling/plasticity, and nervous system development and function 24,25,27,29 . Most of these studies focus on the timeframe of WD, and only 2 studies to our knowledge have examined timepoints out Figure 5. Immunofluorescent corroboration of GSEA and IPA analysis. Longitudinal slices of sciatic nerve were obtained at 10 μm and images acquired at × 200. Log2fold change values from the raw data are presented for transcript to protein comparison. Repeating motifs of neurologic disease, fibrosis, and muscular pathways, were identified through multiple pathway analyses. Neurologic disease was corroborated with axonal (anti-NF200) and myelin sheath (anti-MBP) staining. Expression between these markers is complementary: both demonstrate a dramatic decrease in expression over time. MYH4 is represented in multiple muscle categories, yet as a myosin, has functional multiplicity in the cytoskeleton, which is exemplified by apparent overlap with NF200 at day 0 (D0), D2, and D14, along with control specimens. Fibrosis presents with pathological remodeling as early as 2 days post-injury. Col3a has minute representation in the control nerve, with fluctuations in expression and aberrant deposition over the time course of injury progression. Col1a demonstrates decreased protein expression D2 and D7, followed by increased intensity and disordered remodeling. At D48, both collagens depict intense staining and dearth of microtubule delineation. Scale bar, 50 μm. www.nature.com/scientificreports/ to 4 weeks after injury, reporting a single pathway significantly associated with inflammation at later timepoints, consistent with the interpretation that inflammation is successfully resolving 25,38 . By contrast, our neuroma model provides evidence of persistently activated inflammatory genetic programs. Multiple IPA analyses and GSEA indicate a peak in inflammatory activity at D2 after injury and sustained longitudinal activation though D48. Inflammatory pathway analysis was corroborated by markers of both innate and adaptive immunity, such as CD11b, CD68, and CD4. CD11b, a pan-myeloid marker, demonstrated persistent presence of both mRNA transcripts and histologic evidence at 14 days post-injury-a timepoint traditionally associated with the culmination of WD and inflammatory resolution in recoverable crush injury models 39,40 . The cellular phenotype potentially indicates chronic myeloid cell recruitment to or de-novo proliferation within the neuroma microenvironment, as it is unlikely for these cells to be long-lived 41,42 . CD68 is predominately associated with activated, phagocytic macrophages, while CD4+ cells traditionally annotate T helper cells; both demonstrated prolonged presence in the mature neuroma. Furthermore, transcriptional activity failed to resolve towards the homeostatic levels of control nerves, while DAPI revealed profound hypercellularity of the neuroma.
To our knowledge, no other studies have examined a protracted timepoint (7 weeks post-injury) and associated nerve-specific genetic activity profile. This terminal endpoint is more faithful to the timeline of the clinical disease course and surgical repair, and clinical neuroma specimens similarly note inflammatory cell presence at prolonged timepoints after injury [9][10][11] . Persistent activation of inflammatory pathways is also a common transcriptional correlate of neuropathic pain in higher-order structures (brain 17,43 , spinal cord 44,45 , and dorsal root ganglia 46,47 ) after peripheral nerve injury. Importantly, although the adaptive immune response to nerve injury is incompletely characterized, previous studies suggest the adaptive immune cells have little impact on the kinetics of WD, and suppression may even may facilitate peripheral nerve regeneration 48,49 . Conversely, aberrant activation of adaptive immunity is implicated in chronic pain states, potentially governing steady recruitment of inflammatory species thereby influencing axonal sensitization 50 . Therefore, the persistent activity of the adaptive Figure 6. Immunofluorescent corroboration of temporal inflammatory signatures associated with neuroma formation. Immunofluorescence images acquired at × 400 for inflammatory markers and × 200 for DAPI. Cd11b is a pan-granulocyte marker; robust increase in transcript expression is observed acutely after injury, as corroborated by immunofluorescent staining. Expression persists through day 14 (D14), with gradual decline and minute expression by D48. CD68 is a marker for activated macrophages, expression compared with cd11b is delayed with a peak at D2, and gradual decline in transcript expression through D48. An expanded and "foamy" phenotype is observed for D7 and D14 that attenuates with time. Notably, persistence of CD68 is observed in the pathologic neuroma state at the D48 endpoint. CD4 is commonly associated with T-helper cells yet also is a marker for dendritic cells. Pockets of CD4 labeled cells converge with those of cd11b at the early timepoints. In contrast to transcript expression levels, protein corroboration was not observed for D7 and D14. Importantly, we see the reemergence of positive staining at D48, although the phenotype is markedly divergent. DAPI is a nuclear marker and displays a dynamic expression pattern. At our earliest timepoint, intense staining is observed on the exterior of the nerve, notably with undefined nucleic staining. D2 through D7 demonstrate cellular accumulation on the exterior of the nerve, reflected by increased cellularity internally. By D14, cellular infiltration is profound and this hypercellularity persists through D48. Scale bar, 50 μm. www.nature.com/scientificreports/ immune system may be a critical event in neuroma formation and neuropathic pain pathogenesis in the nerve microenvironment.
Resolution of the inflammatory response is widely accepted as the precursor to proper wound healing to permit coordinated architectural remodeling 51 . This process is largely attributed to phagocytosis of dead and dying cells in a process known as efferocytosis 52 . Clearance of dead and dying cells removes cells secreting Damage-associated Molecular Patterns (DAMPs), cytokines, and chemokines-molecules that propagate the inflammatory response and are implicated in contributing to neuropathic pain 53 . In SCI, the self-limited nature of the inflammatory response is corrupted, and sustained activation is associated with fibrotic scar development and secondary pathogenesis 13,54 . This environment is paralleled in the pathophysiologic peripheral neuroma 39,55 , suggesting that aberrant inflammatory activation or failed resolution may similarly contribute to a microenvironment inhibitory to successful regeneration 56,57 . While synergistic cooperation between immune and nervous systems is critical for successful nerve regeneration, our data indicate an unmitigated presence of inflammation, consistent with histology from human specimens 11 , suggesting that the nature of the inflammatory response may be critical in pathophysiologic neuroma formation. Cellular death. Although interrogation of neuronal cellular death profiles in the spinal cord and dorsal root ganglia has been widely studied, investigation of cellular death cascades at the level of the nerve has been limited 39,58,59 . Notably, there is scant annotation of cellular death pathways in transcriptional studies of favorable regeneration 28,[60][61][62] . In our study, cellular death pathways demonstrated bimodal activation across multiple timepoints and persisted in the mature neuroma microenvironment 48 days after injury. Given that this study is at the level of the nerve, our results capture a death turnover of cells in the neuroma microenvironment, which could include resident immune, glial, or fibrotic cell lineages in addition to hematogenous immune cell influx. This evidence is corroborated by the interaction model, which specifically predicts inhibition of neuronal cell death, despite implicating activation of apoptotic and necrotic pathways acutely after injury (D0 and D2). Thereafter, overall cell death is suppressed at D7, which coincides with additional inflammatory recruitment and the downregulation of neuronal function. Finally, increased annotation to cell death pathways persists through the remaining timepoints and overlaps with sustained activation of neurodegenerative and inflammatory pathways in the mature neuroma environment.
Cellular death pathways are dually implicated in antagonizing or abrogating inflammatory cascades; apoptosis and efferocytosis promote inflammatory resolution while necroptosis and pyroptosis prime a proinflammatory environment 63,64 . Recent work in a model of sciatic nerve transection has demonstrated dysregulated macrophage cell death in the form of pyroptosis biases the cytokine profile towards a more pro-inflammatory microenvironment. When pyroptosis is specifically ablated in macrophages, clearance of apoptotic cells is enhanced and the pro-inflammatory microenvironment is reversed, which translated to increased axonal regeneration and improved functional recovery 58 . Similarly, in SCI, secondary pathogenesis is associated with activation of numerous cell death pathways and exacerbation of the inflammatory response. Consequently, targeted inhibition attenuates the inflammatory response and improves neurologic recovery profiles; thus, a proposed neuroprotective strategy involves targeted inhibition of cellular death pathways to abrogate pathophysiologic sequelae [64][65][66] . Although further mechanistic study is warranted to elaborate specific cellular death pathways in peripheral nerve trauma (e.g., necroptosis vs. apoptosis), our data suggest coordinated signatures of WD may be perturbed by cell death dysregulation and contribute to neuroma formation in the pathophysiologic scenario.
Fibrosis. The canonical hallmark of neuroma is the profound abundance of fibrosis that enlarges and swells the nerve; however, it remains in doubt whether the fibrosis is causal or sequela to the failure for regeneration 67 . In the studies examining transcriptional correlates of successful regeneration, identification of fibrotic transcripts has been limited, meaning annotation to downstream biological pathways is also limited 26,27 . Single-cell RNAseq has identified fibroblasts and phenotype derivatives within the homeostatic and regenerating peripheral nerve, but, because these models do not develop neuromas, enriched fibrosis is not expected [68][69][70] .
In our pathophysiologic neuroma model, we have identified contribution to fibrogenesis in multiple independent transcriptional analyses, which were corroborated by robust chaotic collagen deposition in our histologic assays. We suspect that an even more pronounced attribution to these pathways may be present yet many fibrotic pathways are latent within broad disease states (e.g., cancer, cystic fibrosis). Although extracellular matrix proteins are necessary for reestablishing the endoneurial tubules that provide the nerve mechanical stability and continuity to the distal target organs, aberrant and uncoordinated fibrosis is implicated in inhibiting axonal growth 71 . These events are recapitulated in SCI, where scar-free healing is associated with improved axonal regeneration 72 ; this provides a rationale that the prevention of scar formation may be a viable mechanistic target for therapeutic intervention.
Neurodegeneration and metabolism. Transcriptional analyses of peripheral nerve injury in successful models of regeneration have highlighted pathways of neuronal regeneration and Schwann cell activation [25][26][27] . In our model, we conversely observed downregulation of similar paths, which contrasted against the immense activation of the inflammatory, cell death, and fibrotic pathways. Metabolic pathways were significantly downregulated across all timepoints, which may correlate to dysregulated Schwann cell activity, as prior research has identified that glycolytic shifts support axonal regeneration 73 . Interestingly, the enriched transcripts specific to the rupture injury state predict activation of neuronal and Schwann cell activity acutely after injury. Yet by day 7, these efforts are reversed and inhibited, suggesting a critical timepoint at which genetic programs may switch towards pathophysiologic regeneration. This neuronal inhibition coincides with enriched inflammatory signatures distinct from those of surgical manipulation, as indicated by the bifurcation of transcriptional clusters www.nature.com/scientificreports/ between sham and rupture injury. Collectively, these insights suggest that manipulation of the microenvironment at this critical timepoint could counteract the transition towards pathophysiologic neuroma formation. The role of metabolic processes in regeneration and immunologic reprograming is gaining intellectual traction 74 . Metabolic pathways such as fatty acid biosynthesis and lipid metabolism are notably downregulated after spinal cord injury, and targeted inhibition is associated with decreased axonal regeneration in vitro. Similar downregulation was not observed in a peripheral nerve transection (recoverable) nerve injury model 75 . A recent body of work has used single-cell RNAseq to elaborate the reprogramming of macrophage metabolism from predominately glycolytic to oxidative phosphorylation in a recovering nerve injury model (crush). This metabolic transition is associated with a proinflammatory to inflammatory resolution phenotype in macrophages and is observed about day 7 after injury 76 . Interestingly, metabolic cycles have similarly been implicated in inflammatory pain phenotypes; reduced oxidative phosphorylation is associated with the apogee of inflammatory-associated hyperalgesia whereas activation of oxidative phosphorylation resolves this pain pathogenesis 77 . Collectively, there are interesting metabolic parallels between our pathophysiologic peripheral nerve injury model, comorbidities such as chronic pain, and the recalcitrant regeneration programs of the central nervous system after injury, which contrasts against the metabolic activity of recovering nerve injury models.
Translational impact: targets for therapeutic intervention. As there are currently no therapeutic interventions to prevent NIC formation, the pathways identified may afford further insight for scientific interrogation and opportunities for clinical translation. Given the intricate coordination among the inflammatory, cellular death, and fibrotic pathways, targeting their genesis may be a most promising avenue. For example, we know efferocytosis is critical for inflammatory resolution and that failure to coordinate the immune response can lead to excessive fibrotic deposition 51,52 . Therefore, efforts to explore how the immune response resolves in pathophysiologic nerve injury, such as identifying whether local neutrophil apoptosis occurs or detailing the events governing immune cell egress, are promising next steps.
However, a common clinical challenge is intervention at an appropriate timepoint after injury, which may not lend itself to immunomodulation efforts at the genesis of the inflammatory response. Therefore, a multipronged approach of acute, middle, or late interventions may provide more clinically feasible strategies. For example, a mid-term approach may center on shifting immune cells towards inflammatory resolution by promoting oxidative phosphorylation. Such techniques have already been applied to preclinical models of CNS injury, where studies manipulating the bioenergetic cascades demonstrate improved axonal regeneration and functional recovery in SCI 78,79 . Efforts to understand metabolic pathways in successful models of peripheral nerve regeneration have similarly demonstrated favorable outcomes, yet metabolism remains unexplored in pathophysiologic peripheral nerve regeneration 73,80 . Therefore, metabolic reprogramming to dually promote inflammatory resolution and redirect the microenvironment towards a more proregenerative state is a feasible strategy to overcome temporal restrictions in the clinical population.
Fibrosis is a cardinal feature of neuroma formation and, as such, lends itself as a viable strategy for more protracted injuries where manipulation of the inflammatory response may be a belated intervention. For example, HSP47 is critical for collagen production and is upregulated in a variety of tissues (including nerve) after injury 81,82 . HSP47 inhibition, through both siRNA knockdown and knockout models, reduces fibrosis and scar formation in liver and peritoneum [83][84][85] . Therefore, translation to peripheral nerve trauma may similarly decrease the fibrotic pathogenesis of neuroma formation and improve regenerative outcomes. Alternatively, application of low-concentration enzymes, which break down collagen (e.g., collagenase), may locally help degrade collagen of the neuroma while preserving the functionality of other cells critical to regeneration, such as perineurial glia and Schwann cells 86,87 . This strategy has already been approved by the US Food and Drug Administration for Dupuytren contracture, another disease hallmarked by unmitigated collagen deposition, with recent efforts focused on creating sustained delivery strategies for the gradual release of collagenase 88 . Collectively, a singular intervention may not suffice for the treatment of traumatic peripheral nerve injuries, and a multidimensional, temporal approach may provide the most translatable intervention from the transcriptional signatures identified in our discovery science model.

Limitations
Applying restrictions to differentially expressed genes selects for those that are mathematically significant. Pragmatically, genes with transcriptional levels less than our fold change or significance thresholds could be biologically germane to neuroma formation, and this analysis may fail to capture these biological nuances. We additionally restricted our analysis to protein coding transcripts. Although recent efforts have revealed pivotal roles for noncoding RNAs in the regulation of transcriptional activity, these were not a focus of the current study given there is limited pathway annotation of this knowledge in the analysis software. Finally, bulk RNAseq cannot distinguish to which cell types the transcripts may belong, and it is also difficult to draw direct comparisons to other studies given the heterogeneity of sequencing platforms and bioinformatic analysis pipelines.

Conclusion
Collectively, through transcriptional profiling and histologic validation, we identified significant genetic activity annotated to longitudinal activation of the inflammatory, cellular death, and fibrotic pathways in our pathophysiologic neuroma model. These three processes are intrinsically entwined; inflammation resolves through coordinated cellular death processes, which orchestrate an environment for concerted regeneration as opposed to scar formation. Antithetically, activation of this triad is coupled with persistent downregulation of neuronal and myelination pathways, which are the focus of recoverable nerve injury models. These insights afford a novel lens through which the scientific and clinical community can approach therapeutic manipulation of the immune, www.nature.com/scientificreports/ metabolic, cellular death, and fibrotic processes to prevent aberrant activity and promote the intrinsic regenerative capacity of the peripheral nervous system.

Methods
Nerve injury and tissue processing. Sixty-one male C57BL/6 J mice aged 8-10 weeks were randomly assigned to groups: days 0 (6 h post-injury), 2, 7, 14, and 48 for both sham (n = 6) and rupture injuries (n = 5). Control nerves were harvested where no surgical intervention was performed (n = 6). For sham and control samples, 2 biological replicates of 3 nerves were used to obtain sufficient RNA quantity for library preparation. All experiments were approved by the Institutional Animal Care and Use Committee at the University of Utah (Protocol #19-04,007) and performed in accordance with relevant guidelines and regulations. The study is reported per the recommendations detailed by the ARRIVE guidelines. Surgery and injury were performed as previously reported 33,34 . In brief, animals were anesthetized with isoflurane, and using a dorsal approach, the left hindlimb was opened from knee to sciatic notch via blunt dissection to delicately resect the sciatic nerve from surrounding fascia. For the rapid-stretch rupture injury, weight drop applied to a pulley/hook system permitted precise displacement of nerve 31 . After careful placement of nerve on the hook, stretch rupture was obtained with a stretch excursion of > 15 mm while sham-injured nerves remained on the hook for an equivalent time (~ 30 s). Ruptured nerves were repositioned back into the appropriate anatomical area without interventional coaptation, and the incision was closed with nylon suture. Postoperatively, animals were administered bupivacaine (2 mg/kg, topical) and carprofen (5 mg/kg, subcutaneous) for pain control and enrofloxacin (5 mg/kg, subcutaneous) for infection prevention.

Sample preparation.
Nerves were harvested at the predefined timepoints from sciatic notch to distal trifurcation (15 mm), mechanically homogenized in Trizol reagent (Life Technologies, Carlsbad, CA), and stored at − 80 °C until further processing. Control nerves were harvested immediately from euthanized animals that were subject to no experimental intervention. The sample was then mixed with chloroform (1:4 ratio) to produce an aqueous phase containing the total RNA, which was subsequently isolated and used in all further steps. To isolate total RNA for sequencing, the sample was processed through an RNeasy Mini column (QIAGEN, Valencia, CA) per manufacturer instructions. The concentration and quality of isolated RNA were determined by sampling through Agilent ScreenTape Assay (Agilent Technologies, Santa Clara, CA), and only those with an RNA integrity number ≥ 7.0 proceeded for library preparation. Next-generation RNAseq was performed with the Illumina NovaSeq6000 (San Diego, CA) with 2 × 50 base pair sequences and 25 million read-pairs per sample. Sample preparation, sequencing, and further analysis was performed in a blinded fashion.
Bioinformatics analysis. Raw reads derived from sequencing were optimized, trimmed, and aligned to the reference database using STAR software in two-pass mode to output a binary alignment map file sorted by coordinates. Mapped reads were assigned to annotated genes using featureCounts version 1.6.3 (http:// subre ad. sourc eforge. net/) 89 . DEGs for each injury grade were identified using a 5% false discovery rate with DESeq2 version 1.30.0 (https:// bioco nduct or. org/ packa ges/ relea se/ bioc/ html/ DESeq2. html) and log2fold change compared with control expression levels 90 . Global expression profiles were examined through graphing DEG activity, cluster analysis, principal component (PCA) and interaction model analyses. DEG activity was graphed using GraphPad Prism version 6.0.0 for Mac (GraphPad Software, San Diego, California USA, www. graph pad. com). For the cluster analysis, samples were hierarchically clustered using a Euclidean distance matrix of regularized log-transformed counts as input (RStudio 4.2.1: pheatmap 1.0.12. https:// cran.r-proje ct. org/ web/ packa ges/ pheat map/ index. html). PCA was performed with RStudio 4.2.1: highcharter 0.9.4 (https:// cran.r-proje ct. org/ web/ packa ges/ highc harter/ index. html). For the interaction model, we combined injury and day into a single group and identified differentially expressed genes using a 5% false discovery rate with DESeq2 version 1.30.1 90 . In addition, to check whether the injury effect differs across days, we also used an interaction model with injury × day in the design formula. This permitted identification of statistically significant genes distinct to the rupture-injured state for comparison at timepoints D0 vs. D2 and D2 vs. D7 where PCA Euclidian distance mapping suggested similar cluster profiles for rupture and sham injuries.
Critical pathways associated with pathological remodeling were identified using GSEA (version 4.1.0, Broad Institute, San Diego, CA) and IPA (QIAGEN Inc., https:// digit alins ights. qiagen. com/ IPA). Independent analyses were employed to provide an additional layer of validation. GSEA was performed using pre-ranked RNA differential expression data from DESEQ and ENSEMBL hallmarks and chip platforms, resulting in longitudinal comparisons of up-and downregulated gene sets between rupture and control 91 . For IPA, differentially expressed genes were uploaded into the IPA software, and core analysis was run using a p-value < 0.05 and fold change > 2 or <−2 threshold to identify analysis-ready molecules. These molecules were then overlaid with the Ingenuity pathway knowledge base (IPKB), which employs a right-tailed Fischer's exact test to identify statistically significant canonical signaling pathways and biological functional activity associated with the gene set. IPA network analysis was also implemented to identify de novo gene network signatures. This analysis employs the assumption that highly interconnected interactions are likely reflective of considerable biological function and thus may reveal functions that may not be captured in traditional, a priori analyses (such as the canonical pathways or diseases and functions) 92 . Additionally, the networks generated can also be overlaid with either diseases and functions or canonical pathways to provide context of how the de novo-constructed genetic networks are related to traditionally established biologic and molecular pathways. www.nature.com/scientificreports/ identified in the RNAseq analyses. Specifically, we sought to corroborate the repeating motifs of Neurologic Disease, Skeletal and Muscular System Development and Function (abbreviated "Muscle"), Fibrosis, Inflammatory Response, and Cellularity. Timepoints were determined a priori with an n > 5 per timepoint (total n = 34). Upon harvest, nerves were placed in 4% paraformaldehyde for 30 min then transferred to 30% sucrose overnight before embedding in OCT compound, frozen, and stored at − 80 °C. Longitudinal slices were acquired at 10 μm, and corroborative immunofluorescent staining was performed for genes identified in paths of both GSEA and IPA analyses.

Data availability
The datasets used and/or analyzed during the current study available from the corresponding author on reasonable request. www.nature.com/scientificreports/