Meta-SCI app: a comprehensive transcriptional reference for understanding the severity and progression of spinal cord injury

doi:10.21203/rs.3.rs-4358401/v1

Download PDF

Research Article

Meta-SCI app: a comprehensive transcriptional reference for understanding the severity and progression of spinal cord injury

https://doi.org/10.21203/rs.3.rs-4358401/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background

Spinal cord injury (SCI) is a devastating condition that leads to motor, sensory, and autonomic dysfunction. Current therapeutic options remain limited, emphasizing the need for a comprehensive understanding of the underlying SCI-associated molecular mechanisms.

Results

This study characterized distinct SCI phases and severities at the gene and functional levels, focusing on biomarker gene identification. Our approach involved a systematic review, individual transcriptomic analysis, gene meta-analysis, and functional characterization. We compiled a total of fourteen studies with 273 samples, leading to the identification of severity-specific biomarker genes for injury prognosis (e.g., Srpx2, Hoxb8, Acap1, Snai1, Aadat) and phase-specific genes for the precise classification of the injury profile (e.g., Il6, Fosl1, Cfp, C1qc, Cp). We investigated the potential transferability of severity-associated biomarkers and identified a twelve-gene signature that predicted injury prognosis from human blood samples. We also report the development of MetaSCI-app - an interactive web application designed for researchers - that allows the exploration and visualization of all generated results (https://metasci-cbl.shinyapps.io/metaSCI).

Conclusions

Overall, we present a transcriptomic reference and provide a comprehensive framework for assessing SCI considering severity and time perspectives, all integrated into a user-friendly tool.

Spinal cord injury (SCI) leads to alterations in the motor, sensory, and autonomic systems immediately below the affected spinal segment, which represent a devastating blow to patients' quality of life due to the loss of voluntary movement and associated physiological dysfunctions [1–3]. The World Health Organization (WHO) estimates that 250,000-500,000 people suffer from an SCI each year. Sadly, there is no cure for SCI nowadays, although substantial efforts currently support a wide range of preclinical research aims [4, 5]. Recent reports from the Courtine lab have defined a neuronal subpopulation responsible for functional regeneration after epidural electrical stimulation in mice [6]; nevertheless, we still require concerted efforts in basic research to elucidate the dynamic, complex molecular mechanisms that underlie the distinct stages of SCI. Likewise, an improved understanding of the mechanisms (and their biological relevance) at play at early stages after SCI is needed to control evolution in acute cases and achieve functional recovery in those patients suffering from chronic injuries.

Transcriptomics represents a widely used technology for the study of pathological states such as SCI; this technique provides a comprehensive view of the cellular response to a given condition from a systems biology perspective and allows the identification of differentially expressed genes (DEGs). Unfortunately, statistically-significant DEGs from studies can display distressingly little overlap when distinct research groups study the same biological system [7]. This lack of robustness and reproducibility can lead to conflicting or inconclusive results. Although many transcriptomic studies on SCI exist, we still lack a consensus regarding the description of the injury-induced transcriptomic profile and its functional correlation, which hinders the identification of therapeutic targets with sufficient depth over time and based on injury severity due to the dispersion among injuries/models.

To solve such issues, we propose a novel characterization of SCI through a meta-analysis of transcriptomic studies. A meta-analysis represents an integrative approach that allows the measurement of the effect of interest more precisely compared to individual studies, thus increasing statistical power and considering the variability of individual studies to provide more consistent results [8]. We based our strategy on the systematic review and selection of public transcriptomic studies performed in rats since 2004. We grouped samples based on severity and time after injury, allowing the characterization of SCI in two dimensions. We subsequently processed and analyzed all datasets similarly to avoid introducing biases associated with bioinformatics pipelines [9]. We performed a gene expression meta-analysis to integrate individual DGE analysis results, providing a consensus gene expression signature for each group. We identified specific severity- and phase-associated biomarker genes from the generated gene expression patterns. We then functionally characterized these transcriptomic profiles to estimate the dysregulated pathways, infer transcription factor activity, and construct gene co-expression networks. We validated the results bioinformatically with two external datasets and experimentally by quantitative (q-)PCR. We also explored the potential transferability of severity-associated biomarker genes found in rats to predict injury prognosis using human blood samples. Finally, we created the Meta-SCI app, a platform that allows the research community to access and deeply explore all generated results.

Identifying gene consensus signatures based on severity and time

We first performed a systematic review and selection of studies in GEO [10] and ArrayExpress [11] databases following the Preferred Reporting Items for Systematic Reviews and Meta-Analyses (PRISMA) guidelines [12]. After applying the inclusion and exclusion criteria we finally selected 14 studies and 273 samples for our meta-analysis (Supplementary Fig. 1, Additional File 1). Given the studies diversity, we standardized metadata to create experimental groups with uniform nomenclature. We categorized injury severity into moderate (M) and severe (S) injuries (Additional File 2). We established four injury phase groups based on time after injury and sample collection: acute phase (0–3 days, T1), subacute phase (4–14 days, T2), early chronic phase (15–35 days, T3), and late chronic phase (> 35 days, T4) [3, 13, 14]. Controls consisted of samples from sham-operated or naïve rats (non-operated). Additional File 3 contains sample identifiers and annotations, while Additional File 4 and Supplementary Fig. 2 provides details on sample distribution into the distinct experimental groups.

We followed a common pipeline that included data preprocessing, exploratory analysis, and DGE analysis for the individual analysis of the selected datasets. Since we aimed to characterize SCI in all phases and analyze the influence of injury severity, we compared each injury group against the control group within their study to identify DEGs. The eight comparisons performed and the abbreviations used for simplicity are: M1 (Moderate Acute [T1] vs. Control), M2 (Moderate Subacute [T2] vs. Control), M3 (Moderate Early chronic [T3] vs. Control), M4 (Moderate Late chronic [T4] vs. Control), S1 (Severe Acute [T1] vs. Control), S2 (Severe Subacute [T2] vs. Control), S3 (Severe Early chronic [T3] vs. Control), and S4 (Severe Late chronic [T4] vs. Control).

By conducting a meta-analysis, we integrated the DGE analysis results for each comparison to obtain a transcriptional consensus signature for each experimental group. Supplementary Fig. 3A reports the number of studies per group while Supplementary Fig. 3B shows a high number of significant genes (false discovery rate [FDR] < 0.05) in each comparison in our study. The “Individual analysis” and "Meta-analysis" modules of the Meta-SCI app contains a detailed list of all results.

In order to achieve an overview of the similarities between groups we performed a principal component analysis (PCA) and hierarchical clustering analysis of transcriptional consensus signatures. Together, we observe a separation between severe and moderate injury, with similarities between M1 and S1, suggesting the acute phase as the most distinctive temporal phase irrespective of severity (Fig. 1A and B).

Identification of predictive biomarker genes for severity assessment

To identify severity-specific biomarker genes, we compared transcriptional consensus signatures of the four severe injury groups against the four moderate injury groups (Additional File 5). We identified 282 genes with an FDR < 0.1 (Additional File 6). Interestingly, the hierarchical heatmap revealed a higher gene dysregulation (both up-regulated or down-regulated) in severe compared to moderate injuries (Supplementary Fig. 4), resulting in a clear separation into two distinct heatmap branches based on severity. Figure 2A and Fig. 2B illustrate that the top 10 most significantly altered genes (Srpx2, Hoxb8, Acap1, Snai1, Aadat, Ppic, Lrrc17, Map7, Actg2, and H19) enabled a clear separation between moderate and severe injury, with PC1 explaining 93% of the variance in the PCA. The expression patterns of the top 5 severity-associated biomarker genes demonstrated consistent and more heightened dysregulation over time following severe injury (Fig. 2C), which represents a suitable pattern for injury prognosis.

We divided the list of 282 potential severity-specific genes into up-regulated and down-regulated genes and conducted a protein-protein interaction (PPI) analysis to evaluate functional relationships. The PPI network of down-regulated genes (152) revealed enrichments in functions related to nervous system development and the regulation of trans-synaptic signaling (Supplementary Fig. 5). Meanwhile, the PPI network of up-regulated genes (130) exhibited enrichments in functions related to extracellular matrix organization, animal organ development, and anatomical structure morphogenesis (Supplementary Fig. 6).

Effective group classification using phase-specific biomarker genes

To systematically identify phase-specific biomarker genes (independent of severity), we compared each phase against the remaining phases. Additionally, we investigated potential phase-specific biomarker genes for each severity (Additional File 5). Table 1 reports the top 10 genes for each comparison. The acute phase displayed the most distinct transcriptional signature, with 770 genes at FDR < 0.1, while other comparisons did not yield significant genes at this threshold. Nevertheless, using the top 10 genes with the lowest p-value in each case allowed clear group separation (Supplementary Fig. 7–10 and Additional File 7–10).

We selected the top 10 genes from the T1, T2, and T4 groups to create a phase-specific consensus signature, resulting in robust group separation (Fig. 2D and E). In the PCA plot (Fig. 2D), PC1 separated M1/S1 from the remaining groups, while PC2 classified the remaining groups into pairs, with the more similarity between M3/S3 and M4/S4 indicating that the selected biomarkers allow the identification of the injury phase independently of severity. Figure 2D depicts a similar segregation. While we identified the top 10 potential biomarker genes for T3 (demonstrating the segregation of M3/S3 from remaining groups), we do not consider them appropriate biomarkers (fig. S8) as their characteristic expression pattern depicts a decrease between T2 and T4. This could be explained by a lack of consensus in results due to study heterogeneity rather than a genuine drop in gene expression during this phase. T3 likely represents a transitional phase between T2 and T4 without any expected expression biological process that peaks in T3, as one might anticipate in T2.

Table 1

**Top 10 specific biomarker genes for each severity- and time-associated comparison. Severity**, (S1, S2, S3, S4) vs. (M1, M2, M3, M4); T1, (M1, S1) vs. (M2, M3, M4, S2, S3, S4); T2, (M2, S2) vs. (M1, M3, M4, S1, S3, S4); T3, (M3, S3) vs. (M1, M2, M4, S1, S2, S4); T4, (M4, S4) vs. (M1, M2, M3, S1, S3, S4); M1, M1 vs. (M2, M3, M4); M2, M2 vs. (M1, M3, M4); M3, M3 vs. (M1, M2, M4); M4, M4 vs. (M1, M2, M3); S1, S1 vs. (S2, S3, S4); S2, S2 vs. (S1, S3, S4); S3, S3 vs. (S1, S2, S4); S4, S4 vs. (S1, S2, S3). Genes in bold are upregulated in each comparison, while the others are down-regulated.
Comparison	Genes
Severity	Srpx2, Hoxb8, Acap1, Snai1, Aadat, Ppic, Lrrc17, Map7, Actg2, H19
T1	Fosl1, Il33, Vcam1, Mx2, Mdk, Ldlr, Hmgcr, Alpk3, Dusp5, Serpina3n
T2	Cfp, Slc46a3, Scn8a, Ly86, Casr, C1qc, Slc2a5, Fgd2, Aif1, Slc7a7
T3	Slc28a2, Iqgap1, Cd14, Nt5dc2, Xdh, Msn, Lyve1, Tmem37, Rdh10, Enpp3
T4	Cp, Rab27a, Rab3b, Cybrd1, Thrsp, Plcb1, Pmp22, Acot2, Cyp27a1, RT1-Db1
M1	Hmgcs1, Fdps, Cd68, Gpx2, Alpk3, Crym, Apln, Slc15a2, Pltp, Ldlr
M2	Mmp8, Cxcl13, Cd8a, Casr, C1qb, RT1-CE4, Fap, Cfp, Mmp7, Slc7a7
M3	Sost, Rbp1, Cd14, Esm1, Hck, Plag1, Fam102b, Atf3, Slc36a3, Itgb8
M4	Thrsp, Ccl3, Tmem204, Calml4, Cpxm2, Ckap2, Ttk, Csf3r, G0s2, Mpz
S1	Fosl1, Il6, Kng2, Tnfrsf12a, Upp1, Myc, Gprc5a, C3, Exo1, Chodl
S2	Apoc1, Cfp, Tnr, Slc28a3, C1qa, Slc46a3, Casr, Postn, Slc40a1, Pnlip
S3	Enpp3, Ccl4, Cd36, Lox, Igsf6, Xdh, Folr2, Msr1, Slc28a2, Rnf4
S4	Mmrn1, Esm1, Anxa8, Crabp2, Bnc2, Igfbp6, Slco2a1, Ptgis, Tspan8, Pde5a

The PPI analysis with up-regulated T1-specific genes revealed enrichment in responses to stress, cell cycle, and various metabolic processes (Supplementary Fig. 11). Notably, the up-regulated genes demonstrated an increase in expression after injury that declined over time (as evidenced by gene expression patterns of Fosl1 and Il33 in Fig. 2F, for instance). We also detected genes such as Lcat, Vcam1, or Fxyd1, which become down-regulated in T1 but up-regulated in the following phases. Down-regulated genes displayed enrichments for processes related to nervous system development and lipid metabolism (Supplementary Fig. 12). Overall, these two divergent gene expression profiles represent valuable tools to distinguish T1 from the other stages. By utilizing the top 10 candidate T1-specific genes, PCA1 effectively captures 95% of the variance, resulting in a distinct separation between M1/S1 and the remaining groups (Supplementary Fig. 7). T2-specific genes display gene upregulation peaking in T2 and declining in subsequent phases; the Cfp and C1qc complement genes associated with secondary immune response [15, 16] exemplify this pattern (Fig. 2F). Potential T4-specific biomarker genes displayed an increasing expression pattern from T1 that peaked in T4, as seen for Cp (ceruloplasmin) (Fig. 2F). For T2 and T4, the enrichment of the top 25 genes (arbitrary threshold) again revealed immune response involvement in T2 (Supplementary Fig. 13) and lipid metabolism enrichment in T4 (Supplementary Fig. 14). Overall, our meta-analysis allowed the successful identification of biomarker genes for each injury phase after SCI in independent studies with different rat models.

Functional characterization of gene consensus signatures

We next conducted three different analyses in order to characterize the gene consensus signatures obtained in the meta-analysis. First we performed a functional enrichment analysis for each group via the gene set analysis (GSA). Supplementary Fig. 15A reports the number of significant functions (FDR < 0.05) grouped by database, severity, injury phase, and the direction of enrichment. Venn diagrams for moderate and severe injuries indicated a high number of shared functions between the four time phases for up-regulated and dow-nregulated functions (Supplementary Fig. 15B), suggesting the existence of dysregulated mechanisms induced by the SCI that persist into chronic phases. We used Reactome annotation to classify significantly affected pathways to achieve an overview of the biological categories most affected after SCI (Fig. 2A). The immune system represents one of the most prevalent categories for significantly up-regulated pathways in all comparisons independently of severity. These pathways form part of the core of permanent dysregulated pathways, which can be explained by the massive infiltration of inflammatory cells during the acute phase, the intrinsic microglial activation at the injury site, and a failure to efficiently resolve inflammation at the chronic stage [17]. In contrast, the elevated presence of pathways related to cell cycle, DNA repair, or DNA replication in the acute phase decreases over time, which could be explained by the proliferation of spinal cord resident cells such as microglia and macroglia (astrocytes and oligodendrocyte precursor cells) and infiltrated cells such as pericytes, fibroblast or Schwann cells [18–21].

The permanent loss of neuronal system functions and metabolism characterizes significantly down-regulated pathways. Regarding metabolism-associated functions, the dotplot in Supplementary Fig. 16 reports the down-regulation of functions related to fatty acid synthesis, cholesterol synthesis, and the Krebs cycle in all groups, as previously described [22, 23].

Supplementary Fig. 17 reports in more detail the down-regulation of genes involved in cholesterol synthesis in a hierarchical heatmap, demonstrating the more pronounced level of down-regulation in the chronic compared to the acute phase. Our results agree with the findings of Spann et al. [22], who reported the significant down-regulation of the cholesterol biosynthesis pathway and the down-regulated expression of Hmgcs1, Hmgcr, Cyp51, Idi1, and Fdft1, genes involved in the production of cholesterol. We also observed permanent mitochondrial dysfunction, suggested by the down-regulation of functions/pathways related to mitochondria and ATP synthesis [23] (Supplementary Fig. 18). The totality of the described events acts in detriment of the correct functioning of the nervous system and any regenerative attempts, as indicated by the large number of down-regulated functions related to synapse regulation, neurotransmitter regulation, or nervous system development in all groups (Supplementary Fig. 19–21).

We next estimated the transcription factor (TF) activity levels based on the expression levels of their target genes, with a methodology that considers whether TF-target interactions activate or repress said target gene (Additional File 11). Thus, the up-regulation (respectively down-regulation) of target genes of an activating (respectively inhibiting) TF becomes interpreted as TF activation following injury [normalized enrichment score (NES) > 0]. Conversely, activating (resp. inhibiting) TFs with down-regulated (resp. up-regulated) target genes, are interpreted as lost TF function or deactivation following injury (NES < 0). We observed a higher prevalence of activated vs. deactivated TFs in all scenarios (Supplementary Fig. 22A). These findings underscore the more substantial impact on TF activation during the early stages after SCI, which gradually decreases in the chronic phases. Figure 2B (left panel) depicts consistently dysregulated TFs shared between moderate and severe injury. TFs such as Nfkb1, Stat1, and Stat3 become activated during injury progression, contributing to processes such as inflammatory response modulation [24–26]. This activation is consistent with the persistently up-regulated pathways related to the immune system identified in GSA, highlighting their role in sustaining the inflammatory response post-injury. Interestingly, Srebf1 exhibits an increasing up-regulated gene expression pattern (Fig. 2C); however, this TF becomes differentially-deactivated in all comparisons except M1. Both Srebf1 and Srebf2 (also deactivated) represent activators of lipid metabolism genes [27], which could explain the previously described down-regulation of fatty acid and cholesterol synthesis.

We use the gene consensus signatures matrix (Additional File 16) from moderate and severe injuries separately as input to find genes with similar expression patterns. We identified 19 and 18 clusters of co-expressed genes in moderate and severe injuries, respectively, and constructed a PPI network for each cluster. The extensive constructed networks displayed high levels of interconnection, indicating the functional association of genes with similar expression patterns. To gain deeper insights, we generated 261 subnetworks, conducted a functional enrichment analysis for each subnetwork, and then classified them into manually predefined biological categories of interest. Among the 261 generated subnetworks, we successfully classified 153 (Fig. 5A). The remaining 108 subnetworks comprised only a few proteins without enriched functions. As observed in the GSA results, we found the prominent representation of categories such as nervous system, immune system, and metabolism; moreover, we observed notable functional similarity across severity levels. Figure 2D depicts cluster 11, subnetwork 1 for moderate injury, comprising genes related to the cell cycle (proteins in green, red, and purple also form a physical part of the same complexes). These genes possess an expression pattern characterized by an initial up-regulation that decreases over time (Fig. 2E); as discussed previously, this explains the presence of enriched functions related to the cell cycle in T1 that decrease in subsequent phases.

Consistency in temporal gene expression patterns between the mouse and rat

To validate and compare gene expression patterns obtained in our meta-analysis, we selected the mouse RNA-seq dataset from Li et al., 2022 [28], which possessed an experimental design similar to the studies included in our work. They induced injuries at the T9 level and collected samples from a spinal cord segment at different days post-injury (dpi) that covered the four times phases established in our work. We then followed the same pipeline applied to individual studies to generate gene expression patterns, ensuring their comparability with those generated in our meta-analysis. Additional File 12-S16 reports the DGE analysis results of each comparison and the gene expression patterns matrix.

We then conducted a correlation analysis to compare gene expression patterns between the mouse and rat models (rat moderate injury vs. mouse; rat severe injury vs. mouse; rat moderate injury vs. rat severe injury), which demonstrated a generally positive relationship between them. The density plot of correlation values indicates a median close to 0.5 in three comparisons (Fig. 3A and Additional File 17). A total of 892 genes displayed a strong positive relationship (r > 0.9) between mouse and moderate rat injury, another 892 genes demonstrated a positive relationship between mouse and severe rat injury, and 922 genes presented a positive relationship between moderate rat injury and severe rat injury. Among these, 198 genes exhibited a positive correlation across all three comparisons (Fig. 6A). We next performed a clustering analysis using a curated set of potential phase-specific genes (Fosl1, Il6, Cfp, Sipa1, Cp, and Dbp) that display similar gene expression patterns across the three comparisons (Fig. 3B). Hierarchical clustering in Fig. 3C revealed a clear separation of groups based on time post-injury similar to observed in Fig. 1G.

In this case, the mouse T3 and T4 signatures become grouped within the same branch as M3 and S3 of the rat, which could be explained due to slight differences in timings. Mouse T3 is defined as 28 dpi, and T4 as 42 dpi, which rat T3 covers 15 to 35 dpi, and T4 extends from 56 to 168 dpi. Longer intervals between phases could lead to more distinct differences in expression profiles, as the injury has additional time to stabilize.

Our comparative analysis with an SCI mouse model validated the gene expression patterns established in our meta-analysis and demonstrated a positive correlation between mouse and rat injury responses, supporting the reliability of the identified biomarker genes.

Our proposed severity- and phase-specific biomarker genes effectively classify novel rat transcriptomic profiles

We used the GSE218088 dataset [29] published after the study search and selection period as a validation set, given that it meets the inclusion criteria established in the meta-analysis. We aimed to evaluate whether the selected phase- and severity-specific biomarker genes could classify a new transcriptomic profile correctly. In this new study, the authors indicated that the free fall of a 10 g hammer with a diameter of 2.5 mm produced the injury, which is considered severe [30, 31]. The authors indicated that they extracted spinal cord tissue three days after SCI, corresponding to the acute phase (T1) defined in our work.

Using the top 10 severity-specific genes (Table 1), we observed the clustering of the transcriptional signature of the GSE218088 dataset with the severe injury groups in our analysis (Fig. 3D). Additionally, using phase-specific biomarkers for T1, T2, and T4 (Table 1), we observed the clustering of the transcriptional profile of GSE218088 with the T1 signatures from our meta-analysis (Fig. 3E). Therefore, this comparison demonstrates the robustness of the results obtained and the capacity of the biomarker genes to classify a new transcriptional profile.

Concordance between meta-analysis gene expression patterns and qPCR analysis

We aimed to confirm through experimental validation the expression patterns identified in our meta-analysis. Following severe injury induction, we extracted spinal cord RNA from injured rats at T1, T2, T3, and T4 post-injury. We selected eight genes to validate their expression profiles by qPCR. We chose four TFs (Cebpd, Hes1, Jun and Onecut1) due to their functional relevance, whereas the other four genes (Il6, Vcam1, Vegfa, and Vim) represent targets of these TFs and potential phase-specific biomarker genes (Fig. 4).

In particular, Il6 exhibited general overexpression with notably higher expression in T1 than in other phases, in which it remained similar to pre-injury levels. Vcam1 presented down-regulation in T1 followed by a gradual up-regulation expression increase over time, except in the mouse dataset, where we only appreciate the gradual expression increase over time. Oppositely, Vegfa showed up-regulation in T1 and down-regulation in the other phases. Vim confirmed the constant up-regulation, although qPCR did not validate its potential role as S4-specific biomarker.

With respect to TFs, Cebpd, Hes1 and Jun presented similar up-regulated patterns in qPCR, meta-analysis and mouse dataset. Although there is lower consistency in the expression patterns across all four systems, we observe a general agreement. Contrarily, Onecut1 presented different patterns of change but confirmed its constant down-regulation over time.

In addition, we conducted a correlation analysis to compare the logFC values obtained from qPCR with those derived from the meta-analysis and mouse dataset (Supplementary Fig. 23). For target genes, 10 of 12 qPCR comparisons had strong Pearson correlations above 0.85. In contrast, for TFs only 5 of 11 comparisons exceeded 0.7.

Taken together, although reproducing the exact expression pattern remained challenging, this validation confirmed general dysregulation patterns or at least the direction of dysregulation at specific phases after injury.

Severity-specific rat biomarker genes predict injury prognosis from blood samples from human SCI patients

To investigate the potential clinical translation of severity-specific biomarker genes identified in rats to the prediction of injury prognosis from blood samples of human SCI patients, we analyzed the GSE151371 dataset [32], which contains gene expression profiles of human blood samples collected in the first 10 days after SCI. Human SCIs are categorized based on severity and prognosis using the American Spinal Cord Injury Association (ASIA) with a 5-grade scale- the American Impairment Scale (AIS) - ranging from A-E, moving from more to less severe. Grade A denotes complete sensory and motor loss; Grade B signifies complete motor loss but preserved sensation; Grade C and D represent various degrees of motor function preservation; while Grade E indicates normal sensory and motor function [33].

To achieve a similar scenario to our rat-based meta-analysis, we selected samples from AIS A and AIS D to ensure similarity in severity between the two biological systems. After DGE analysis between groups, we identified 931 genes (FDR < 0.1; Additional File 18); we compared this list with the 282 severity-specific genes previously identified in rat models (Fig. 5A). We found 12 intersecting genes (EXT1, FBN1, FOSB, GNAO1, GRB10, MMP9, NFE2, PRF1, SLC25A23, SNAI1, ST6GALNAC3 and STX11) (Fig. 5B), which expression pattern in the GSE151371 dataset is characterized by an increasing dysregulation as SCI severity increases (Supplementary Fig. 24).

We then used this 12-gene signature to analyze its ability to stratify groups according to severity by clustering analysis in rats and humans. Human blood sample analysis revealed the presence of three clusters: severe injury, moderate injury, and healthy controls (Fig. 5C), while the gene consensus signatures demonstrated stratification based on the severity level in rats (Fig. 5D). These findings indicate that expression levels of the 12 identified biomarker genes can delineate severity-specific subgroups in human blood and rat spinal cord samples.

Meta-SCI app - an interactive, user-friendly platform to explore results

We have developed the Meta-SCI app (https://metasci-cbl.shinyapps.io/metaSCI/), an interactive web application that provides a user-friendly interface to consult and visualize the results obtained in this work, allowing users to gain a deeper understanding of the underlying mechanisms of SCI. This application consists of eight main modules (Fig. 6), allowing users not only to explore all the generated results but also perform further analysis and research with the available tools. For further information on the modules of Meta-SCI app consult the “Help” section.

SCI presents a complex challenge that demands novel approaches for a better understanding of its underlying mechanisms. Here, we present a transcriptomic meta-analysis consolidating the most significant number of studies/samples regarding rat SCI datasets to our knowledge. Furthermore, we integrated diverse experimental models based on severity and phase post-injury. As a pioneering effort, we also established a freely accessible transcriptomic reference for SCI research.

The strength of our work emerges from the generation of transcriptomic consensus signatures within established groups. The application of uniform processing and integration through meta-analysis enhances the statistical robustness in determining the magnitude of gene alteration, mitigating inherent dataset heterogeneity and biases stemming from disparate pipelines – a significant feat, given that studies with identical experimental designs can yield contradictory expression changes for the same gene [34]. Experimental validation represents another strength of our study, which we demonstrated via i) qPCR of eight gene signatures, confirming in all cases the direction of dysregulation; ii) bioinformatic analysis, finding a positive correlation in temporal transcriptional changes between rat and mouse; and iii) the correct classification of an external transcriptional signature using the proposed severity and

phase biomarker genes. Our work facilitates the comparison of gene expression levels, pathways and TF activity across time/injury severities. Furthermore, direct identification of gene expression levels associated with specific pathways or TF targets is possible. While one of our study's strengths lies in assembling 14 studies and 273 samples, unequal representation across experimental groups represents the main limitation, highlighting an emphasis on acute and a lack of samples in chronic phases. As far as our knowledge extends, Squier et al. [35] previously integrated the larger number of rat SCI transcriptomic datasets. Nevertheless, our research represents a progress as we incorporated a more comprehensive array of diverse datasets (14 vs. 3) and employed more sophisticated integration methods using meta-analysis instead of a basic intersection of results. This underscores the necessity to publish data following FAIR principles (findability, accessibility, interoperability, and reusability) [36], as certain studies do not provide data or adequate sample descriptions. Another potential limitation arises from including datasets derived from diverse technologies and platforms (10 microarrays and 4 RNA-seq experiments). However, the standardized reanalysis and subsequent integration after DGE analysis effectively captured changes in transcriptional regulation [9].

This study generated eight consensus transcriptional profiles that comprehensively characterized each predefined group. These transcriptional profiles serve a dual purpose: first, they provide insight into phase- and severity-specific expression patterns of genes, and second, they are used for the identification of potential biomarkers and functional characterization. The robustness of these expression profiles is reinforced by different sources. Previously, Tica et al. [37] combined four transcriptomic and proteomic datasets to identify genes/proteins persistently dysregulated at 7 days and 8 weeks post injury. They found 40 up-regulated and 48 down-regulated genes/proteins, which consistently matched our meta-analysis results. Notably, our research identified one of these genes - Vim - as consistently up-regulated across all time points after SCI, corroborating results from the meta-analysis, the mouse dataset, and qPCR validation. We also validated the expression profiles of the meta-analysis by comparing them with the expression profiles of a mouse dataset [28], which shared a very similar experimental design with the studies included in our work. We observed an overall positive correlation between mouse and rat gene expression patterns for moderate and severe injury. While the exact replication of expression patterns posed challenges, qPCR-based validation reaffirmed the direction of dysregulation (up-regulation or down-regulation) at specific time points. Overall, these transcriptomic consensus signatures provide insight into how gene expression patterns evolve, offering valuable information regarding the development of potential therapeutic strategies suitable for each SCI phase.

We successfully identified biomarker genes that clearly distinguish moderate and severe injury, displaying significant and persistent dysregulation in severe injuries. The down-regulated genes play roles in the functions and development of the nervous system and synapses. These findings mirror the increased loss of function associated with severe injury due to more extensive damaged tissue and complete axonal rupture. The up-regulated genes are associated with tissue remodeling and extracellular matrix organization. Notably, 45 of these genes are translated into proteins that will be located in the extracellular region. In severe injuries these genes might become activated to regulate mechanisms involved in wound healing and fibrosis [38]. By employing the top 10 severity-specific genes, we effectively classified an external dataset involving a severe rat injury, demonstrating their validity. The expression profiles of the identified severity-specific genes make them suitable as potential prognostic biomarkers. Since no significant expression changes occurred in moderate damage, but in severe injuries there is considerable and persistent dysregulation, these genes can be employed at the early stages of injury to predict its progression. Among potential severity biomarkers, we identified matrix metalloproteinases (MMPs), Mmp9 and Mmp23, which exhibited significantly higher up-regulation in severe compared to moderate injury. Previous studies suggested MMPs as potential predictive biomarkers for worse neurological outcomes in human patients as well as in rat and canine models [38–42]. This consistency across species underlines the reliability and potential application in clinical research as prognosis biomarkers. In addition, Microtubule-Associated Protein 2 (Map2) has previously been proposed as a severity biomarker [40–43]. Our results indicate a greater down-regulation of the expression of Map2 (and Map7) in severe injury, which could represent a new candidate biomarker gene for future studies.

We also identified phase-specific biomarker genes that allowed for phase-specific stratification of transcriptomic profiles regardless of severity. The acute phase (T1) - the most characteristic and readily distinguishable phase - provided the highest number of potential biomarker genes. This may be explained by The cascade of events triggered by the initial trauma, including the activation of stress and inflammatory response mechanisms to contain the injury [1, 17]. Additionally, we also observe a great activation of genes related to the cell cycle, driven by the proliferation/activation of mitotic cells (such as microglia and astrocytes), as well as the initiation of apoptotic processes in damaged cells [18, 19]. In fact, it has been reported that inhibition of cell cycle gene expression may improve motor and cognitive recovery [44–47]. Subacute phase-specific genes, including complement-related genes (e.g., Cfp, C1qa, C1qb, and C1qc), exhibited a distinct gene expression pattern peaking at T2. This pattern aligns with the dynamic behavior of cells such as macrophages and microglia that express and secrete complement proteins [15, 16]. Monocytes migrate to injury sites and differentiate into macrophages at approximately 3 dpi, culminating in a peak response at around 7 dpi while resident microglia reach an activation peak 1 week post-injury [48, 49]. Utilizing a combination of phase-specific genes supports the accurate classification of transcriptomic profiles for meta-analysis data and new transcriptomic profiles, such as those in the GSE218088 dataset. Additionally, while some of these genes may display species-specificity, common genes (such as Fosl1, Il6, Cfp, Sipa1, Cp, and Dbp) can be used for the temporal classification of rat and mouse samples. Experimental validation of Il6, Vcam1, and Vegfa confirmed the expression patterns observed in the meta-analysis and validated their similarity in the mouse dataset. Thus, our study identified phase-specific biomarkers that may aid in selecting therapies based on transcriptomic profiles and improve patient stratification.

Our work provides a functional framework for assessing the degree of dysregulation in biological functions of interest. Our findings align with the findings of previous SCI research, characterized by a maintained up-regulation of functions related to the immune system and down-regulation of processes associated with the nervous system function and development. This phenomenon is influenced by other processes that sustain an unfavorable environment for functional recoveries, such as disturbed ionic homeostasis, vascular injury, ischemia, free radical stress, cell death [13, 17, 50], mitochondrial dysfunction [23], metabolic alterations (including the down-regulation of lipid and cholesterol metabolism) [22], or extracellular matrix remodeling [51]. Considering this complex context, any strategy for developing effective therapies should include a combination of multiple targets, considering the interconnected nature of these systems. At the functional level, we observed differences between phases, with the decreasing significance of pathways related to cell cycle and RNA transcription as time progresses post-injury particularly noteworthy. However, in contrast to gene level results, we did not identify substantial differences between severities. Likewise, the constructed co-expression networks, comprising different genes with altered expression based on severity, demonstrated similarities in the enriched functions' quantity and nature. This analysis corroborates the functional interconnection of genes with similar expression patterns [52, 53] and agrees with the study of De Biase et al. [54], which failed to observe consistent alterations in functional patterns associated with injury severity, even given the apparent variations occurring at the gene expression level. This phenomenon may be explained by the design of functional enrichment techniques that capture the coordinated action of gene groups. Consequently, even though specific genes may display more pronounced dysregulation in severe injuries, their collective direction points toward the same outcome.

Our study also provides a way to infer the activity of master transcriptional regulators represented by pioneering evaluation of TF activity based on the expression levels of target genes, giving a new perspective for more efficient therapeutic interventions. In addition, we experimentally validated four TFs. Hes1 - a transcriptional repressor involved in neurogenesis [55] – displayed overexpression in our meta-analysis and the experimental validation and differential activation in all experimental groups. Previous studies support these findings, demonstrating Hes1 overexpression after SCI [56] and emphasizing the role of Hes1 in promoting neurogenesis through inhibition of its target genes, consequently enhancing functional recovery [57]. Onecut1 – a transcriptional activator regulating the production, diversification, distribution, and maintenance of various neuronal populations [58–60] – became down-regulated following SCI; while we experimentally validated this result, the results were not wholly consistent with the expression pattern obtained from the meta-analysis, possible by an inter-experimental deviations. Activity inference analysis indicated deactivation, especially in the early phases of injury. These findings suggest overexpression or increased activity levels could benefit neuronal development and function. Our bioinformatic and experimental results also indicated an up-regulation of Cebpd - a transcriptional activator that regulates genes involved in immune responses [61] - after injury, accompanied by differential activation in acute phases. Previous studies reported Cebpd overexpression following injury [61, 62], along with enhanced expression during functional recovery in mice deficient in this gene [62]. Jun - a transcriptional activator highly induced following neuronal damage [63, 64] – became consistently overexpressed and activated in all experimental groups; in a related study, Zhang et al. proposed that miR-152 overexpression inhibited inflammatory responses and promoted functional recovery by inhibiting Jun [64].

We found of a great interest the correlation within the expression profiles in human samples and the rat spinal cord transcriptional data. We identified a 12-gene signature predicting injury severity in human blood and rat spinal cord data. Among these 12 genes, MMP9 has been proposed as a potential severity biomarker in species including humans, rats, and dogs [39, 41, 65]. The clinical application of this signature could be valuable, enabling patient stratification through simple techniques such as blood extraction and transcriptional analysis. Although the ASIA system, accompanied by the anatomical evidence of the extent of the affected and preserved tissue by MRI, represents the standard for predicting SCI severity, there exists a need to find new operative and non-invasive methods for its prognosis. The ASIA system relies on neurological status evaluations through sensory and motor exams; however, this method suffers from limitations [32]. Further investigation will be required to explore the prognostic capacity of these genes in detail with larger sample sizes in humans and additional relevant models.

The development of the Meta-SCI app during this study offers an additional dimension. This platform extends beyond a mere results repository; Meta-SCI serves as an analysis tool for the research community. Also, Meta-SCI includes two injury models, allowing researchers to choose that which best suits their needs. Consequently, the Meta-SCI app emerges as a valuable resource for hypothesis generation and validation and constitutes a transcriptional reference database for studying SCI.

Our study offers a comprehensive systems biology frame of SCI, providing a holistic view of the gene and functional levels and the dimensions of severity and phase. We identified potential biomarker genes that effectively stratified transcriptomic profiles based on time and severity, offering critical insights into SCI progression. Thanks to our user-friendly Meta-SCI web application, researchers can explore this extensive results repository, tailoring their investigations to their fields of interest. Our study sheds light on the intricate mechanisms underlying SCI and equips researchers with the means to translate these insights into diagnostic and therapeutic interventions.

Study search and selection

A systematic review and selection of studies were conducted in September of 2022. The search was performed in two public databases for transcriptomic data - GEO [10] and ArrayExpress [11] - according to PRISMA [12] statement guidelines and carried out with the keywords "spinal cord injury." The inclusion criteria were: (1) organism, Rattus norvegicus; (2) tissue, thoracic spinal cord; and (3) RNA sequencing or microarray data from Affymetrix or Agilent used as gene expression platforms. The exclusion criteria were: (1) studies without control samples (sham or naive rats); (2) individual samples from rats that had received any type of treatment; and (3) samples rostral or caudal to the epicenter of the injury.

We classified injury severity into moderate and severe injuries. The severe group included lesions resulting from a spinal cord impact of ≥ 200 kdyn, a 10 g weight drop from 50 mm, or a complete transection of the spinal cord; meanwhile, we considered injuries not meeting these criteria as moderate (Additional File 2). Additionally, we defined four phases of injury based on the timeline post-injury and when samples were collected: the acute phase (0–3 days, T1), the subacute phase (4–14 days, T2), the early chronic phase (15–35 days, T3), and the late chronic phase (over 35 days, T4). Control groups included samples from rats either sham-operated or naïve rats.

Data preprocessing

For microarray studies, normalized data were downloaded with the GEOquery package and probes were annotated to their corresponding gene symbol using platform-specific annotation packages. For RNA-seq studies, raw count matrices were manually downloaded from GEO and identifiers were converted to gene symbols using the org.Rn.eg.db package. In all cases, the median of the expression values of duplicated identifiers was calculated.

Exploratory analysis

An exploratory analysis was performed (including boxplot, PCA, and clustering analysis) to observe sample distribution in the groups of interest and evaluate anomalous behavior and possible batch effects. Samples from the GSE183591 study were taken in two different series, separated into two blocks in the PCA. This batch effect was taken into consideration for DGE analysis. The negative values detected in boxplots for the GSE2599 and GSE464 datasets were normalized by adding the minimum of all values followed by a logarithmic transformation of the data.

Differential gene expression analysis

Each injury group was compared with its corresponding control within the study to identify DEGs. DGE analysis of microarray platforms was carried out using the limma package [66]. To ensure comparable results between technologies, RNA-seq datasets were analyzed using the limma-voom pipeline, which involves the removal of genes with low counts and a logarithmic transformation of the expression matrix. p-values were corrected with the Benjamini–Hochberg (BH) method [67].

Gene expression meta-analysis

DGE results were integrated to obtain a transcriptional consensus signature for each comparison. A meta-analysis was performed for each gene in at least two studies in each comparison. The random-effects model proposed by Der Simonian & Laird [68], implemented in the metafor package [69], was used to evaluate the combined effect. This method weighs the heterogeneity of each study and incorporates their inherent variability into the overall estimate of the measure of effect. Consequently, studies with greater variability will have lower weight in this measure. In this case, logFC was used to measure the effect, while the variance was used to measure variability. Thus, values for the p-value, logFC, and 95% confidence interval (CI) were calculated for each gene evaluated in the meta-analysis. Since multiple meta-analyses were performed, p-values were adjusted using the BH method. Adjusted p-values lower than 0.1 were considered significant. Funnel and forest plots were used to assess the variability and measure each study's contribution to the meta-analysis.

Clustering analysis

All clustering analyses were performed on the gene consensus signature matrix. We constructed this matrix with consensus logFC values from the meta-analysis, selecting genes evaluated in all groups and significant in at least one group. The mixOmics package for PCA was used [70]. For hierarchical clustering analysis, the Euclidean distance was first calculated as input for the hclust function.

Biomarker gene identification

In this context, biomarker genes were defined as those displaying more significant deregulation in certain groups of interest than others, allowing the stratification of different groups based on their expression pattern. A limma test was applied to identify possible biomarker genes that could classify samples according to severity/phase. The four severe were compared against the four moderate groups of samples to elucidate severity-specific biomarker genes. Phase-specific biomarker genes common to both severities and specific to moderate or severe injury were also explored. Detailed comparisons are listed in Additional File 4. Subsequently, the top 10 genes with the lowest p-values from each comparison were selected to assess their classification ability via clustering analysis comprising PCA and hierarchical clustering.

Gene set analysis

A functional enrichment analysis was performed for each gene consensus signature with the GSA method implemented in the mdgsa R package [71]. A gene ranking and functional annotation are required as input for GSA. The ranking is made by ordering all genes according to the p-value obtained in meta-analysis and the logFC value sign. For functional enrichment annotation, three annotation databases have been included: the biological processes of Gene Ontology [72], KEGG pathways [73], and Reactome pathways [74]. Gene sets with less than 10 or more than 500 genes were excluded. P-values were corrected with the BH method.

Transcription factors activities

The identification of differentially active TFs used the msviper function of the viper package [75] with each of the gene consensus signatures as input. Mouse regulons from the DoRothEA package [76] with a confidence level of A, B, C, or D were selected, excluding those with less than 10 genes (Additional File 11). The p-values were corrected with the BH method.

Construction of gene co-expression networks

First, in order to find genes with similar expression patterns we used the clust software [77] with the following arguments: -n 0 -t 100 -cs 4. Gene consensus signatures matrix from moderate and severe injuries were used separately as input. A PPI was then constructed for each cluster of co-expressed genes obtained from clust software. The STRINGdb package with the interactions of STRING version 11.5 [78] and medium confidence of 0.400 were used to build PPI networks. The PPI network of each cluster was divided into different subnetworks of more closely related genes using the fastgreedy algorithm implemented in the get_clusters function.

Each subnetwork was then classified into manually pre-defined functional categories of interest. The get_enrichment function was first used to obtain overrepresented functional terms in each subnetwork. Subsequently, the biological processes were selected from the Gene Ontology, KEGG, and Reactome pathways, ad hoc classifying them into pre-established general categories. Finally, the frequency distribution of each category was calculated for each subnetwork, assigning the category with the highest frequency as that which best defines the network.

Bioinformatic validation datasets

For the mouse transcriptomic dataset, a normalized count matrix provided from the supplementary material provided by Li et al. [28] (https://doi.org/10.6084/m9.figshare.17702045) was downloaded. Samples representative of the four temporal phases established in our study (1, 7, 28, and 42 dpi) and controls were selected (Additional File 19). For the GSE218088 dataset, the normalized expression matrix was downloaded using GEOquery, and the probes were annotated to the gene symbol. Those genes that displayed significance in at least one comparison and were present in all groups of our meta-analysis were selected for both datasets. A DGE analysis using the limma package was performed, comparing each injured group against uninjured controls.

Experimental validation (qPCR)

A severe traumatic SCI model in adult female Sprague Dawley rats was employed for RNA isolation and qPCR in silico data validation. The animals were housed at the Animal Experimentation Unit of the Research Institute Príncipe Felipe (Valencia, Spain) under standard conditions. All experimental procedures adhered to the guidelines established by the European Communities Council Directive (86/609/ECC), the Spanish Royal Decree 53/2013, and the Animal Care and Use Committee of the Research Institute Príncipe Felipe (2021/ VSC/PEA/0032).

Rats were subcutaneously pre-medicated with morphine (2.5 mg/kg) and anesthetized with 2% isoflurane in a continuous oxygen flow of 1 L/min. Laminectomy was performed on thoracic vertebrae 8–9 to expose the spinal cord, and severe SCI was induced at the thoracic vertebrae 8 level by contusion, applying a force of 250 kdyn using the Infinite Horizon Impactor, as previously described (13, 94). Post-surgery care included manual bladder drainage twice a day until vesical reflex recovery and subcutaneous administration of 5 mg/kg of Enrofloxacin (Alsir) for 7 days, as well as 0.1 mg/kg of Buprenorphine twice a day for 4 days after each intervention.

At each endpoint: 1 dpi (n = 3), 7 dpi (n = 3), 28 dpi (n = 4), and 56 dpi (n = 4) after injury to represent (T1, T2, T4, and T8, respectively), animals were overdosed of sodium pentobarbital (100 mg/kg) and transcardially perfused with a 0.9% saline solution. The spinal cord tissue was extracted, and the injury epicenter was immediately frozen in liquid nitrogen and stored at – 80 ºC until use. RNA extraction was performed using the TriZol standard method, followed by an additional cleanup step using RNeasy MinElute Cleanup (Qiagen, Germany) to ensure sample quality (A260/280 ≈ 2 and A260/230 ≥ 1.8). Reverse transcription was performed using the high-capacity RNA-to-cDNA™ kit (Applied Biosystems, Massachusetts, USA).

Specific primers (Additional File 20) for each gene of interest were designed using primer-BLAST (NCBI, Maryland, USA) and validated by efficiency curve performance. qPCR was conducted in triplicate using AceQ SYBR qPCR Master Mix (ThermoFisher) in the Light-Cycler 480 detection System (Roche, Basel, Switzerland). Ct data were calculated with the LightCycler 480 relative quantification software (Roche, Basel, Switzerland). GAPDH mRNA levels served as an internal control for normalization.

Analysis of human blood samples

The GSE151371 normalized gene expression matrix was downloaded from GEO, and the expression values were log-transformed. The samples corresponding to the AIS A, AIS D, and the healthy control groups were then selected for further analysis (Additional File 21). DGE analysis was then performed to compare AIS A vs. AIS D with the limma package following the same pipeline as the previous analyses. Finally, significantly altered genes (FDR < 0.1) were selected to compare with the list of rat biomarker genes. The list of intersecting genes was used for the clustering analysis of both biological systems. The healthy controls from the GSE151371 dataset were also included for visualization and clustering analysis.

Meta-SCI app

Meta-SCI app is powered by RStudio Shiny package and deployed on a shinyapps.io server, available at https://metasci-cbl.shinyapps.io/metaSCI. Plots are generated using ggplot, plotly, and heatmaply. All data processing and analysis are carried out using R.

Acknowledgments

The authors thank the Principe Felipe Research Center (CIPF) for providing access to the cluster, co-funded by European Regional Development Funds (FEDER) in Valencian Community 2014-2020. The authors also thank Stuart P. Atkinson for reviewing the manuscript.

Funding

This research was supported and partially funded by the Institute of Health Carlos III (project IMPaCT-Data, IMP/00019), co-funded by ERDF, "A way to make Europe", PID2021-124430OA-I00 funded by MCIN/AEI/10.13039/501100011033 and by “ERDF A way of making Europe".

Author contributions:

Conceptualization: FGG

Methodology: RGR, MRH, BMR, VMM, FGG

Investigation: FGG, MRH, BMR, VMM, RGR

Visualization: RGR, MRH

Supervision: MRH, VMM, FGG

Writing—original draft: RGR, MRH, BMR, VMM, FGG

Writing—review & editing: RGR, MRH, BMR, VMM, FGG

Availability of data and materials

All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Material. All data and results generated in this study are freely available through the Meta-SCI app: https://metasci-cbl.shinyapps.io/metaSCI. Transcriptomic data are available at GEO under accession numbers: GSE183591, GSE182796, GSE138637, GSE133093, GSE115067, GSE102964, GSE104317, GSE93249, GSE52763, GSE45006, GSE46988, GSE29488, GSE2599, GSE464, GSE151371 and GSE218088. The code developed for the analyses described in this work is publicly available at GitLab (https://gitlab.com/Biorubens/meta-sci-app).

Competing interests

Authors declare that they have no competing interests.

Silva NA, Sousa N, Reis RL, Salgado AJ. From basics to clinical: A comprehensive review on spinal cord injury. Prog Neurobiol. 2014 Mar;114:25–57.
Kang Y, Ding H, Zhou H, Wei Z, Liu L, Pan D, et al. Epidemiology of worldwide spinal cord injury: a literature review. J Neurorestoratology. 2017 Dec;Volume 6:1–9.
Ahuja CS, Wilson JR, Nori S, Kotter MRN, Druschel C, Curt A, et al. Traumatic spinal cord injury. Nat Rev Dis Primer. 2017 Apr 27;3(1):17018.
Courtine G, Sofroniew MV. Spinal cord repair: advances in biology and technology. Nat Med. 2019 Jun;25(6):898–908.
Nas K. Rehabilitation of spinal cord injuries. World J Orthop. 2015;6(1):8.
Kathe C, Skinnider MA, Hutson TH, Regazzi N, Gautier M, Demesmaeker R, et al. The neurons that restore walking after paralysis. Nature. 2022 Nov 17;611(7936):540–7.
Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci. 2005 Oct 25;102(43):15545–50.
Normand SLT. Meta-analysis: formulating, evaluating, combining, and reporting. Stat Med. 1999 Feb 15;18(3):321–59.
Ramirez Flores RO, Lanzer JD, Holland CH, Leuschner F, Most P, Schultz J, et al. Consensus Transcriptional Landscape of Human End‐Stage Heart Failure. J Am Heart Assoc. 2021 Apr 6;10(7):e019667.
Barrett T, Wilhite SE, Ledoux P, Evangelista C, Kim IF, Tomashevsky M, et al. NCBI GEO: archive for functional genomics data sets—update. Nucleic Acids Res. 2012 Nov 26;41(D1):D991–5.
Kolesnikov N, Hastings E, Keays M, Melnichuk O, Tang YA, Williams E, et al. ArrayExpress update—simplifying data submissions. Nucleic Acids Res. 2015 Jan 28;43(D1):D1113–6.
Liberati A, Altman DG, Tetzlaff J, Mulrow C, Gøtzsche PC, Ioannidis JPA, et al. The PRISMA Statement for Reporting Systematic Reviews and Meta-Analyses of Studies That Evaluate Health Care Interventions: Explanation and Elaboration. PLoS Med. 2009 Jul 21;6(7):e1000100.
Mataliotakis GI, Tsirikos AI. Spinal cord trauma: pathophysiology, classification of spinal cord injury syndromes, treatment principles and controversies. Orthop Trauma. 2016 Oct;30(5):440–9.
Kjell J, Olson L. Rat models of spinal cord injury: from pathology to potential therapies. Dis Model Mech. 2016 Oct 1;9(10):1125–37.
Merle NS, Noe R, Halbwachs-Mecarelli L, Fremeaux-Bacchi V, Roumenina LT. Complement System Part II: Role in Immunity. Front Immunol [Internet]. 2015 May 26 [cited 2023 Oct 10];6. Available from: http://www.frontiersin.org/Molecular_Innate_Immunity/10.3389/fimmu.2015.00257/abstract
Lubbers R, Van Essen MF, Van Kooten C, Trouw LA. Production of complement components by cells of the immune system. Clin Exp Immunol. 2017 Apr 6;188(2):183–94.
Alizadeh A, Dyck SM, Karimi-Abdolrezaee S. Traumatic Spinal Cord Injury: An Overview of Pathophysiology, Models and Acute Injury Mechanisms. Front Neurol. 2019 Mar 22;10:282.
Di Giovanni S, Knoblach SM, Brandoli C, Aden SA, Hoffman EP, Faden AI. Gene profiling in spinal cord injury shows role of cell cycle in neuronal death. Ann Neurol. 2003 Apr;53(4):454–68.
Wu J, Stoica BA, Faden AI. Cell Cycle Activation and Spinal Cord Injury. Neurotherapeutics. 2011 Apr;8(2):221–8.
Zhou Y, Wang Z, Li J, Li X, Xiao J. Fibroblast growth factors in the management of spinal cord injury. J Cell Mol Med. 2018 Jan;22(1):25–37.
Kanno H, Pearse DD, Ozawa H, Itoi E, Bunge MB. Schwann cell transplantation for spinal cord injury repair: its significant therapeutic potential and prospectus. Rev Neurosci [Internet]. 2015 Jan 1 [cited 2023 Sep 27];26(2). Available from: https://www.degruyter.com/document/doi/10.1515/revneuro-2014-0068/html
Spann RA, Lawson WJ, Grill RJ, Garrett MR, Grayson BE. Chronic spinal cord changes in a high-fat diet-fed male rat model of thoracic spinal contusion. Physiol Genomics. 2017 Sep 1;49(9):519–29.
Rabchevsky AG, Michael FM, Patel SP. Mitochondria focused neurotherapeutics for spinal cord injury. Exp Neurol. 2020 Aug;330:113332.
Karova K, Wainwright JV, Machova-Urdzikova L, Pisal RV, Schmidt M, Jendelova P, et al. Transplantation of neural precursors generated from spinal progenitor cells reduces inflammation in spinal cord injury via NF-κB pathway inhibition. J Neuroinflammation. 2019 Dec;16(1):12.
Wang T, Yuan W, Liu Y, Zhang Y, Wang Z, Zhou X, et al. The role of the JAK-STAT pathway in neural stem cells, neural progenitor cells and reactive astrocytes after spinal cord injury. Biomed Rep. 2015 Mar;3(2):141–6.
Satriotomo I, Bowen KK, Vemuganti R. JAK2 and STAT3 activation contributes to neuronal damage following transient focal cerebral ischemia. J Neurochem. 2006 Sep;98(5):1353–68.
Bengoechea-Alonso MT, Ericsson J. The phosphorylation-dependent regulation of nuclear SREBP1 during mitosis links lipid metabolism and cell growth. Cell Cycle. 2016 Oct 17;15(20):2753–65.
Li C, Wu Z, Zhou L, Shao J, Hu X, Xu W, et al. Temporal and spatial cellular and molecular pathological alterations with single-cell resolution in the adult spinal cord after injury. Signal Transduct Target Ther. 2022 Mar 2;7(1):65.
Cao J, Tang T, Tan J, Chen Q, Yuan J, Li T, et al. Expression Profiles of Long Noncoding RNAs and Messenger RNAs in a Rat Model of Spinal Cord Injury. V E S, editor. Comput Math Methods Med. 2023 Jan 19;2023:1–18.
Verma R, Virdi JK, Singh N, Jaggi AS. Animals models of spinal cord contusion injury. Korean J Pain. 2019;32(1):12.
Khuyagbaatar B, Kim K, Kim YH. Conversion Equation between the Drop Height in the New York University Impactor and the Impact Force in the Infinite Horizon Impactor in the Contusion Spinal Cord Injury Model. J Neurotrauma. 2015 Dec 15;32(24):1987–93.
Kyritsis N, Torres-Espín A, Schupp PG, Huie JR, Chou A, Duong-Fernandez X, et al. Diagnostic blood RNA profiles for human acute spinal cord injury. J Exp Med. 2021 Mar 1;218(3):e20201795.
Roberts TT, Leonard GR, Cepela DJ. Classifications In Brief: American Spinal Injury Association (ASIA) Impairment Scale. Clin Orthop. 2017 May;475(5):1499–504.
Merikangas AK, Shelly M, Knighton A, Kotler N, Tanenbaum N, Almasy L. What genes are differentially expressed in individuals with schizophrenia? A systematic review. Mol Psychiatry. 2022 Mar;27(3):1373–83.
Squair JW, Tigchelaar S, Moon KM, Liu J, Tetzlaff W, Kwon BK, et al. Integrated systems analysis reveals conserved gene networks underlying response to spinal cord injury. eLife. 2018 Oct 2;7:e39188.
Wilkinson MD, Dumontier M, Aalbersberg IjJ, Appleton G, Axton M, Baak A, et al. The FAIR Guiding Principles for scientific data management and stewardship. Sci Data. 2016 Mar 15;3(1):160018.
Tica J, Bradbury E, Didangelos A. Combined Transcriptomics, Proteomics and Bioinformatics Identify Drug Targets in Spinal Cord Injury. Int J Mol Sci. 2018 May 14;19(5):1461.
Haggerty AE, Marlow MM, Oudega M. Extracellular matrix components as therapeutics for spinal cord injury. Neurosci Lett. 2017 Jun;652:50–5.
De Castro R, Burns CL, McAdoo DJ, Romanic AM. Metalloproteinase increases in the injured rat spinal cord: NeuroReport. 2000 Nov;11(16):3551–4.
Albayar AA, Roche A, Swiatkowski P, Antar S, Ouda N, Emara E, et al. Biomarkers in Spinal Cord Injury: Prognostic Insights and Future Potentials. Front Neurol. 2019 Jan 29;10:27.
Kwon BK, Casha S, Hurlbert RJ, Yong VW. Inflammatory and structural biomarkers in acute traumatic spinal cord injury. cclm. 2011 Mar 1;49(3):425–33.
Rosell A, Cuadrado E, Ortega-Aznar A, Hernández-Guillamon M, Lo EH, Montaner J. MMP-9–Positive Neutrophil Infiltration Is Associated to Blood–Brain Barrier Breakdown and Basal Lamina Type IV Collagen Degradation During Hemorrhagic Transformation After Human Ischemic Stroke. Stroke. 2008 Apr;39(4):1121–6.
Zhang SX, Underwood M, Landfield A, Huang F fa, Gibson S, Geddes JW. Cytoskeletal Disruption Following Contusion Injury to the Rat Spinal Cord. J Neuropathol Exp Neurol. 2000 Apr;59(4):287–96.
Wu J, Stoica BA, Dinizo M, Pajoohesh-Ganji A, Piao C, Faden AI. Delayed cell cycle pathway modulation facilitates recovery after spinal cord injury. Cell Cycle. 2012 May;11(9):1782–95.
Di Giovanni S, Movsesyan V, Ahmed F, Cernak I, Schinelli S, Stoica B, et al. Cell cycle inhibition provides neuroprotection and reduces glial proliferation and scar formation after traumatic brain injury. Proc Natl Acad Sci. 2005 Jun 7;102(23):8333–8.
Kabadi SV, Stoica BA, Byrnes KR, Hanscom M, Loane DJ, Faden AI. Selective CDK Inhibitor Limits Neuroinflammation and Progressive Neurodegeneration after Brain Trauma. J Cereb Blood Flow Metab. 2012 Jan;32(1):137–49.
Skovira JW, Wu J, Matyas JJ, Kumar A, Hanscom M, Kabadi SV, et al. Cell cycle inhibition reduces inflammatory responses, neuronal loss, and cognitive deficits induced by hypobaria exposure following traumatic brain injury. J Neuroinflammation. 2016 Dec;13(1):299.
Milich LM, Ryan CB, Lee JK. The origin, fate, and contribution of macrophages to spinal cord injury pathology. Acta Neuropathol (Berl). 2019 May;137(5):785–97.
Kopper TJ, Gensel JC. Myelin as an inflammatory mediator: Myelin interactions with complement, macrophages, and microglia in spinal cord injury. J Neurosci Res. 2018 Jun;96(6):969–77.
Chamankhah M, Eftekharpour E, Karimi-Abdolrezaee S, Boutros PC, San-Marina S, Fehlings MG. Genome-wide gene expression profiling of stress response in a spinal cord clip compression injury model. BMC Genomics. 2013 Dec;14(1):583.
Gaudet AD, Popovich PG. Extracellular matrix regulation of inflammation in the healthy and injured spinal cord. Exp Neurol. 2014 Aug;258:24–34.
Piro RM, Ala U, Molineris I, Grassi E, Bracco C, Perego GP, et al. An atlas of tissue-specific conserved coexpression for functional annotation and disease gene prediction. Eur J Hum Genet. 2011 Nov;19(11):1173–80.
Pierson E, the GTEx Consortium, Koller D, Battle A, Mostafavi S. Sharing and Specificity of Co-expression Networks across 35 Human Tissues. Rigoutsos I, editor. PLOS Comput Biol. 2015 May 13;11(5):e1004220.
De Biase A, Knoblach SM, Di Giovanni S, Fan C, Molon A, Hoffman EP, et al. Gene expression profiling of experimental traumatic spinal cord injury as a function of distance from impact site and injury severity. Physiol Genomics. 2005 Aug 11;22(3):368–81.
Sueda R, Imayoshi I, Harima Y, Kageyama R. High Hes1 expression and resultant Ascl1 suppression regulate quiescent vs. active neural stem cells in the adult mouse brain. Genes Dev. 2019 May 1;33(9–10):511–23.
Wang X, Wang Q, Tian H, Lv W, Song L, Li Z, et al. Electroacupuncture in promoting neural repair after spinal cord injury: Inhibiting the Notch signaling pathway and regulating downstream proteins expression. Anat Rec. 2021 Nov;304(11):2494–505.
Ruan W, Ning G, Feng S, Gao S, Hao Y. MicroRNA‑381/Hes1 is a potential therapeutic target for spinal cord injury. Int J Mol Med [Internet]. 2018 May 8 [cited 2023 Oct 10]; Available from: http://www.spandidos-publications.com/10.3892/ijmm.2018.3658
Kabayiza KU, Masgutova G, Harris A, Rucchin V, Jacob B, Clotman F. The Onecut Transcription Factors Regulate Differentiation and Distribution of Dorsal Interneurons during Spinal Cord Development. Front Mol Neurosci. 2017 May 26;10:157.
Roy A, Francius C, Rousso DL, Seuntjens E, Debruyn J, Luxenhofer G, et al. Onecut transcription factors act upstream of Isl1 to regulate spinal motoneuron diversification. Development. 2012 Sep 1;139(17):3109–19.
Toch M, Harris A, Schakman O, Kondratskaya E, Boulland JL, Dauguet N, et al. Onecut-dependent Nkx6.2 transcription factor expression is required for proper formation and activity of spinal locomotor circuits. Sci Rep. 2020 Jan 22;10(1):996.
Chevreau R, Ghazale H, Ripoll C, Chalfouh C, Delarue Q, Hemonnot-Girard AL, et al. RNA Profiling of Mouse Ependymal Cells after Spinal Cord Injury Identifies the Oncostatin Pathway as a Potential Key Regulator of Spinal Cord Stem Cell Fate. Cells. 2021 Nov 27;10(12):3332.
Wang SM, Hsu JYC, Ko CY, Chiu NE, Kan WM, Lai MD, et al. Astrocytic CCAAT/Enhancer-Binding Protein Delta Contributes to Glial Scar Formation and Impairs Functional Recovery After Spinal Cord Injury. Mol Neurobiol. 2016 Nov;53(9):5912–27.
Raivich G, Bohatschek M, Da Costa C, Iwata O, Galiano M, Hristova M, et al. The AP-1 Transcription Factor c-Jun Is Required for Efficient Axonal Regeneration. Neuron. 2004 Jul;43(1):57–67.
Zhang T, Gao G, Chang F. miR-152 promotes spinal cord injury recovery via c-jun amino terminal kinase pathway. Eur Rev Med Pharmacol Sci. 2019 Jan;23(1):44–51.
Nagano S, Kim SH, Tokunaga S, Arai K, Fujiki M, Misumi K. Matrix metalloprotease-9 activity in the cerebrospinal fluid and spinal injury severity in dogs with intervertebral disc herniation. Res Vet Sci. 2011 Dec;91(3):482–5.
Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015 Apr 20;43(7):e47–e47.
Benjamini Y, Hochberg Y. Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing. J R Stat Soc Ser B Methodol. 1995 Jan;57(1):289–300.
DerSimonian R, Laird N. Meta-analysis in clinical trials. Control Clin Trials. 1986 Sep;7(3):177–88.
Viechtbauer W. Conducting Meta-Analyses in R with the metafor Package. J Stat Softw [Internet]. 2010 [cited 2023 Oct 18];36(3). Available from: http://www.jstatsoft.org/v36/i03/
Rohart F, Gautier B, Singh A, Lê Cao KA. mixOmics: An R package for ‘omics feature selection and multiple data integration. Schneidman D, editor. PLOS Comput Biol. 2017 Nov 3;13(11):e1005752.
Montaner D, Dopazo J. Multidimensional Gene Set Analysis of Genomic Data. Hoheisel J, editor. PLoS ONE. 2010 Apr 27;5(4):e10348.
Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene Ontology: tool for the unification of biology. Nat Genet. 2000 May;25(1):25–9.
Kanehisa M. The KEGG resource for deciphering the genome. Nucleic Acids Res. 2004 Jan 1;32(90001):277D – 280.
Vastrik I, D’Eustachio P, Schmidt E, Joshi-Tope G, Gopinath G, Croft D, et al. Reactome: a knowledge base of biologic pathways and processes. Genome Biol. 2007;8(3):R39.
Alvarez MJ, Shen Y, Giorgi FM, Lachmann A, Ding BB, Ye BH, et al. Functional characterization of somatic mutations in cancer using network-based inference of protein activity. Nat Genet. 2016 Aug;48(8):838–47.
Garcia-Alonso L, Holland CH, Ibrahim MM, Turei D, Saez-Rodriguez J. Benchmark and integration of resources for the estimation of human transcription factor activities. Genome Res. 2019 Aug;29(8):1363–75.
Abu-Jamous B, Kelly S. Clust: automatic extraction of optimal co-expressed gene clusters from gene expression data. Genome Biol. 2018 Dec;19(1):172.
Szklarczyk D, Gable AL, Nastou KC, Lyon D, Kirsch R, Pyysalo S, et al. The STRING database in 2021: customizable protein–protein networks, and functional characterization of user-uploaded gene/measurement sets. Nucleic Acids Res. 2021 Jan 8;49(D1):D605–12.

Supplementary Tables are not available with this version

No competing interests reported.

SupplementaryFigures.docx

Download PDF

Version 1

posted

You are reading this latest preprint version

Meta-SCI app: a comprehensive transcriptional reference for understanding the severity and progression of spinal cord injury

Status:

Version 1

Abstract

Background

Results

Conclusions

Figures

Background

Results

Identifying gene consensus signatures based on severity and time

Identification of predictive biomarker genes for severity assessment

Effective group classification using phase-specific biomarker genes

Functional characterization of gene consensus signatures

Consistency in temporal gene expression patterns between the mouse and rat

Our proposed severity- and phase-specific biomarker genes effectively classify novel rat transcriptomic profiles

Concordance between meta-analysis gene expression patterns and qPCR analysis

Severity-specific rat biomarker genes predict injury prognosis from blood samples from human SCI patients

Meta-SCI app - an interactive, user-friendly platform to explore results

Discussion

Conclusions

Methods

Study search and selection

Data preprocessing

Exploratory analysis

Differential gene expression analysis

Gene expression meta-analysis

Clustering analysis

Biomarker gene identification

Gene set analysis

Transcription factors activities

Construction of gene co-expression networks

Bioinformatic validation datasets

Experimental validation (qPCR)

Analysis of human blood samples

Meta-SCI app

Declarations

References

Supplementary Tables

Additional Declarations

Supplementary Files

Status:

Version 1