Characterization of lncRNA and mRNA profiles in the process of repairing peripheral nerve defects with cell-matrixed nerve grafts

Abstract

Background

Decellularized extracellular matrix (dECM) is an intriguing natural biomaterial that has garnered significant attention due to its remarkable biological properties. In our study, we employed a cell-matrixed nerve graft for the repair of sciatic nerve defects in rats. The efficacy of this approach was assessed, and concurrently, the underlying molecular regulatory mechanisms were explored to elucidate how such grafts facilitate nerve regeneration. Long noncoding RNAs (lncRNAs) regulate mRNA expression via multiple mechanisms, including post-transcriptional regulation, transcription factor effects, and competitive binding with miRNAs. These interactions between lncRNAs and mRNAs facilitate precise control of gene expression, allowing organisms to adapt to varying biological environments and physiological states. By investigating the expression profiles and interaction dynamics of mRNAs and lncRNAs, we can enhance our understanding of the molecular mechanisms through which cell-matrixed nerve grafts influence neural repair. Such studies are pivotal in uncovering the intricate networks of gene regulation that underpin this process.

Results

Weighted gene co-expression network analysis (WGCNA) utilizes clustering algorithms, such as hierarchical clustering, to aggregate genes with similar expression profiles into modules. These modules, which potentially correspond to distinct biological functions or processes, can subsequently be analyzed for their association with external sample traits. By correlating gene modules with specific conditions, such as disease states or responses to treatments, WGCNA enables a deeper understanding of the genetic architecture underlying various phenotypic traits and their functional implications. We identified seven mRNA modules and five lncRNA modules that exhibited associations with treatment or time-related events by WGCNA. We found the blue (mRNAs) module which displayed a remarkable enrichment in “axonal guidance” and “metabolic pathways”, exhibited strong co-expression with multiple lncRNA modules, including blue (related to “GnRH secretion” and “pyrimidine metabolism”), green (related to “arginine and proline metabolism”), black (related to “nitrogen metabolism”), grey60 (related to “PPAR signaling pathway”), and greenyellow (related to “steroid hormone biosynthesis”). All of the top 50 mRNAs and lncRNAs exhibiting the strongest correlation were derived from the blue module. Validation of key molecules were performed using immunohistochemistry and qRT-PCR.

Conclusion

Revealing the principles and molecular regulatory mechanisms of the interaction between materials and biological entities, such as cells and tissues, is a direction for the development of biomimetic tissue engineering technologies and clinically effective products.

Background

Peripheral nerve injury (PNI) is a common clinical condition that affects millions of people worldwide each year [1]. It is primarily caused by accidents, trauma, and other factors, leading to loss of sensory and motor functions, accompanied by neuropathic pain and target organ atrophy [2]. Despite extensive research and understanding of the regenerative mechanisms of peripheral nerve injury, reliable treatment methods to ensure complete structural and functional recovery are still lacking. In recent years, many biodegradable polymers such as polylactic acid and poly (lactic-co-glycolic acid) (PLGA) have been widely used in the research of tissue-engineered nerve scaffold materials [3, 4]. Although these materials have made some progress in the repair of peripheral nerve defects, they still have deficiencies in cell interaction sites and cell adhesion [5].

Decellularized extracellular matrix (dECM) has been widely used in the field of tissue engineering due to its tissue-specific composition of extracellular matrix (ECM) [6]. ECM is a three-dimensional polymer network composed of proteins and polysaccharides, which possesses complex morphological and structural signals that finely regulate cellular processes, signaling pathways, gene expression profiles, and injury repair [7]. Through the decellularization process, dECM retains the morphology that mimics the natural microenvironment and eliminates most antigenic substances, effectively preventing pathogen transmission [8, 9].

Following established protocols, we harvested a dECM from human bone marrow mesenchymal stem cells (hBMSCs) that retained intact conformation and biological functionality. This material was then utilized to encapsulate a multilayered PLGA fiber scaffold, successfully creating a three-dimensional, matrix-oriented nerve graft. This engineered graft was deployed for the reconstructive bridging of a 10 mm defect in the rat sciatic nerve, subsequently undergoing a comprehensive evaluation to determine the functional recovery of the regenerated nerve. To assess pain sensitivity, which encompasses responses to thermal and mechanical stimuli, both thermal and mechanical pain tests were implemented. Gait and motor function were scrutinized using CatWalk gait analysis, which measured gait parameters and symmetry. Furthermore, neuroelectrophysiological tests were conducted to evaluate nerve conduction function. The extent of muscle atrophy was determined by comparing the wet weight ratio of the impaired muscle to the normal muscle. Key morphological features of nerve regeneration, such as the cross-sectional area of the regenerated nerve axons and the thickness of the myelin sheath, were quantitatively assessed. Additionally, we conducted a thorough statistical evaluation of the muscle fiber cross-sectional area to ascertain the degree of muscle atrophy and performed motor endplate maturity analysis to determine the morphological and functional maturity of the neuromuscular junction, as depicted in Fig. 1. This holistic approach allowed a nuanced exploration of the therapeutic efficacy of the novel nerve graft across multiple physiological dimensions. The cell-matrixed nerve graft provides a favorable microenvironment for axonal growth and achieves functional recovery similar to autografts in nerve regeneration [10]. During the process of nerve regeneration, significant changes in various genes [11,12,13,14], as well as in microRNAs (miRNAs), long non-coding RNAs (lncRNAs), and circular RNAs (circRNAs) are often accompanied [15, 16].

Fig. 1

Implantation of cell-matrixed nerve grafts in vivo and corresponding evaluation. A The schematic diagram of the cell-matrixed nerve graft. The cell-matrixed nerve grafts, as described in Chinese Patent ZL202210982152.0, were fabricated from chitosan conduits with an internal diameter of 2 mm and PLGA scaffolds constructed from 140 to 160 monofilaments (provided by Jiangsu Yitong Biotechnology Company, Nanjing, China). A gelatinous membrane-like hBMSC-dECM, with a thickness of approximately 0.12 mm and PLGA serving as the axis, was rolled into 6 to 8 layers and threaded into the chitosan conduits to create the cell-matrixed nerve grafts prior to surgery. B Schematic illustration of the transplantation of cell-matrixed nerve grafts into rat sciatic nerve defects. C Behavioral testing. The gait analysis was used to assess the recovery of motor function, while thermal and mechanical sensitivity tests evaluated the recovery of sensory function. D Electrical stimuli with an intensity of 10 mA were applied to the sciatic nerve trunk at the proximal and distal ends of the graft using a portable digital MYTO electrophysiological apparatus. The compound muscle action potential (CMAP) of the gastrocnemius muscle was recorded during this process (left part). Morphological assessment of the regenerative nerve and target muscles (right part)

LncRNAs, characterized by their length exceeding 200 nucleotides and lack of significant protein-coding capacity, have been shown to play a role in neural regeneration, including the regulation of gene expression, modulation of signaling pathways, and interaction with other molecules [17,18,19,20,21]. However, approaches adopted in such articles usually examine individual genes or lncRNAs, whereas they exert their functions through the co-expression network showing consistent biological functions.

Weighted gene co-expression network analysis (WGCNA) has been widely employed to uncover the relationships between gene-based connections and phenotypes using microarray data or RNA sequencing in various samples [22,23,24]. This comprehensive approach encompasses the identification of modules comprising highly correlated genes, summarizing these modules through an eigengene network, and evaluating the correlations between modules and external sample traits using eigengene network methodology. Moreover, WGCNA allows for the identification of intramodular hub genes that play key roles within modules [25]. This approach provides a more comprehensive and integrated understanding of gene expression patterns, allowing for a more nuanced exploration of the underlying biological mechanisms.

The present study collected RNA sequencing data from the sciatic nerve graft segments in rats to elucidate significant co-expression modules between the human bone marrow mesenchymal stem cell-derived decellularized extracellular matrix (hBMSC-dECM) group and the sham group. Next, we performed module-trait relationship analysis to select important modules. Subsequently, we conducted functional and pathway enrichment analysis of the modules using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG). In the significant modules, we also constructed mRNA-lncRNA co-expression networks. The blue module was further analyzed because it was not only associated with treatment and time features but also located at the core of the co-expression network with lncRNAs. Currently, there is no systematic study using WGCNA to construct mRNA-lncRNA co-expression networks related to sciatic nerve repair in nerve grafts.

These findings provide the scientific evidence and resource for better understanding of the molecular mechanisms that orchestrate the regenerative processes in peripheral nerves through the cell-matrixed nerve graft.

Methods

Preparation of the rat sciatic nerve defect model

Healthy male adult Sprague-Dawley (SD) rats, aged 8–10 weeks, weighing 200–220 g, which were purchased from the Experimental Animal Center of Nantong University, were used in this study. Rats were randomly divided into two groups: the dECM group undergoing repair with cell-matrixed nerve grafts; and the sham group with exposure of the sciatic nerve without defect. Each experimental group (n = 9) for evaluation at eight time points (day 1, day 3, day 7, week 2, week 3, week 4, week 8, week 12). The surgical approach was performed as described previously [10]. During the surgery, all animals were subjected to deep anesthesia (3 ml/kg) by intraperitoneal injection of a composite anesthesia solution, which included sodium chlorate, magnesium sulfate, pentobarbital sodium, ethanol, and propylene glycol. Before exposing the sciatic nerve in the middle of the left thigh, the skin was shaved off and disinfected with iodine. Then, approximately 8 mm of the nerve was excised near the bifurcation of the sciatic nerve, causing it to retract and form a 10 mm defect. Then, the prepared nerve graft was sutured to the defect using 8 − 0 micro sutures, followed by muscle and skin sutures. Iodine disinfection was performed after surgery and rats were transferred to a specific pathogen-free (SPF) facility. All experimental procedures followed the Guide for the Care and Use of Laboratory Animals issued by the US National Research Council and approved by the Laboratory Animal Management Committee of Jiangsu Province, China.

RNA extraction and cDNA library analysis

Bridged segments connecting the proximal and distal junctions were sampled thrice from each experimental group (dECM and sham groups) on day 1, 3, and 7, as well as on weeks 2, 3, 4, 8, and 12. Following the manufacturer’s protocol, a total of 48 ribonucleic acid (RNA) samples were extracted using Trizol (Life technologies, Carlsbad, CA, USA). The purity of the RNA samples was evaluated using NanoDrop 2000 (Thermo Fisher, USA), while the integrity was assessed using Agilent 2100 (Agilent Technologies, Santa Clara, CA, USA), respectively. The samples with RNA Integrity Number (RIN) of greater than or equal to 7 were analyzed subsequently. Strand-specific RNA-seq libraries were generated, and subsequently, the samples were subjected to sequencing using the Illumina Novaseq 6000 platform (Gene Denovo, Guangzhou, China). All subsequent analyses were conducted using clean reads. The raw sequencing data from this study have been deposited in the Genome Sequence Archive in BIG Data Center (https://bigd.big.ac.cn), Beijing Institute of Genomics (BIG), Chinese Academy of Sciences, under the accession number CRA010000.

Principal component analysis

Principal component analysis (PCA) was performed on the two groups using all genes expressed. With R package gmodels (https://www.rproject.org/), PCA was performed to determine the principal components and to identify the differences in gene expression patterns among these two groups. Each sample was differentially colored.

Differential expression analysis

Differential expression analysis was performed using the edgeR package (http://www.rproject.org/) to identify genes that were differentially expressed between the two groups. Significantly differentially expressed genes (DEGs) were defined as those with a false discovery rate (FDR) < 0.05 and |log₂FC|>1, as determined from the analysis results. Venn diagrams were plotted utilizing the VennDiagram R package (https://bioinformatics.psb.ugent.be/webtools/Venn/).

Construction of co-expression networks

Gene co-expression network analysis was specifically performed on sciatic nerve bridging segments using the R package WGCNA. In short, the correlation network was first constructed by creating a Pearson correlation matrix between all pairs of genes. Subsequently, when the soft threshold power (β) value (range 1 to 20) was at least 5 for mRNAs and at least 8 for lncRNAs, the corresponding scale-free R² value was 0.8 for both mRNAs and lncRNAs to obtain a good scale-free topology model. Then, to capture the network connectivity, topological overlap matrix (TOM) was calculated [26]. Using various measures of TOM, we employed average linkage hierarchical clustering to group genes with similar expression patterns into modules. To ensure robustness, each module contained a minimum of 50 genes. The dynamic tree cutting algorithm was utilized to determine the initial module assignments. Subsequently, we employed a combined threshold function with a threshold of 0.6 for mRNA and 0.75 for lncRNA.

Calculating module eigengenes, module membership and intramodular connectivity

The module eigengene (ME), represented by the first principal component of each module, was used to assess the gene’s contribution within the module. The module membership (MM) of a gene was determined by calculating the Pearson correlation between the gene and its corresponding ME. Genes with higher MM values are considered to have a more significant role within the module. This rigorous selection criterion allowed us to prioritize genes that are likely to have a substantial impact on the module’s functionality. Modules exhibiting both high independence and distinct connectivity effects compared to other modules were identified for subsequent analysis of module-trait relationships. The intra-module connectivity of a gene within a specific module was quantified as the sum of correlation coefficients between the module and other nodes. Genes displaying substantial module connectivity are considered to have a pivotal role within their respective modules.

Enrichment analysis

The DEGs were annotated with GO terms (https://www.geneontology.org/), and the number of genes associated with each term was calculated. Pathway enrichment analysis was carried out with KEGG. A hypergeometric P-value was calculated to assess the significance of enrichment compared to the background genes, and the FDR was controlled at a threshold of ≤ 0.05 for multiple testing correction.

Visualization of co-expression of mRNA-lncRNA modules

The co-expression relationship between mRNA and lncRNA: Pearson correlation coefficient method was used to analyze the expression correlation between lncRNA and mRNA among samples. The relationship pairs with an absolute correlation value greater than 0.95 were considered as mRNA-lncRNA co-expression relationships. Positive correlation results and negative correlation results were plotted separately. We Selected modules related to time or treatment in mRNA and lncRNA, and used circos (version 0.69-3) to draw a positive correlation trans interaction network between mRNA and lncRNA, used R (version 3.6.0) package circlize (version 0.4.13) to draw a negative correlation trans interaction network between mRNA and lncRNA. Interaction networks were draw by Cytoscape (3.7.0) from the top 10/50 lncRNAs and mRNAs with the most trans-target relationships. Blue module network: We conducted a correlation analysis between the gene expression levels in the blue module of the mRNA WGCNA results and the first principal component, which represents the module eigengene value. We calculated the MM for the top 50 genes and used Cytoscape to draw an interaction network, with the top 5 MM genes highlighted in red. Additionally, we utilized the ClueGO plugin in Cytoscape, which allows for KEGG enrichment analysis and visualization of functionally related genes as a network.

Immunohistochemical staining

At 12 weeks after surgery, bridging segments of the sciatic nerve in the rats of hBMSC-dECM group were removed before euthanasia, and cut into 12 μm thick longitudinal sections. To validate the expression and tissue localization of proteins encoded by genes, the immunofluorescent triple staining was performed as described previously [11]. Briefly, the sections were allowed to incubate with primary antibodies: mouse anti-NF200 antibody (1:200, Sigma) with rabbit anti-Lgals3 antibody (1:200, SAB); rabbit anti-Atp1a3 antibody (1:200, Abcam), and mouse anti-Hmox1 antibody (1:200, Immunoway) at 4℃ for 24 h, followed by further reaction with the secondary antibody (Goat anti-Mouse IgG-Alex-488, 1:500; Goat anti-Rabbit IgG-Cy3, 1: 1000) at room temperature for 1 h. The nuclei were stained by DAPI (SouthernBiotech). The resultant images were observed under fluorescence microscopy (AxioImager M2, Zeiss).

Quantitative real-time polymerase chain reaction (qRT-PCR)

Total RNA was reverse transcribed using the PrimeScript RT Reagent Kit (TaKaRa Bio). Gene expression was quantified through qRT-PCR utilizing the TB Green Premix Ex Taq (TaKaRa) on a StepOnePlus RT-PCR system (Applied Biosystems). Relative mRNA levels were determined using the comparative 2^−∆∆Ct method and normalized to Gapdh levels. Each experiment was conducted three times for accuracy. Primer sequences can be found in Table S1.

Results

Differential expression of mRNAs and lncRNAs

Both the dECM and sham groups demonstrated a significant overlap in the expression of mRNAs and lncRNAs, with 17,396 and 5,020 shared transcripts, respectively. Both at the mRNA and lncRNA levels, the differentially expressed genes shared by the dECM and sham groups far exceed those unique to each group (Fig. 2A). Through PCA analysis, we observed that the overall mRNA expression of dECM group and sham group had their own characteristics. Whereas for lncRNA expression, most of the time points in the dECM group clustered together with the sham group, except for weeks 2, 3 and 12 (Fig. 2B).

We assessed the differential expression of lncRNAs and mRNAs between the two groups at day 1 and seven other time points. The differences in both mRNAs and lncRNAs expression at various time points were more pronounced in the dECM group compared to the sham group. In the sham group, differentially expressed genes emerged rapidly, peaking at a minor crest around the 2-week mark before stabilizing in subsequent weeks. In contrast, the peak occurrence of differentially expressed genes in the dECM group was observed at 3 weeks. Over time, both the sham and dECM groups exhibited a proportional increase in differentially expressed genes, converging at the 12-week time point where the count of differentially expressed genes was approximately double that observed at 3 days (Fig. 2C). In the dECM group, the number of up-regulated genes reached its peak at week 2 and remained relatively stable thereafter. As for lncRNAs, both up-regulated and down-regulated lncRNAs in the dECM group exhibited their highest expression levels on day 1 and week 2 (Fig. 2D).

Fig. 2

The discrete expression patterns of mRNAs (left column) and lncRNAs (right column). A Venn diagram of detected mRNA and lncRNA genes in two groups. B PCA of two groups across the eight time points based on normalized mRNAs and lncRNAs expression level. C Bar plot presentation of differentially expressed mRNAs and lncRNAs at different time points in two groups (compared with that at day 1). D The number of up- or down-regulated differentially expressed mRNAs and lncRNAs at different time points in the hBMSC-dECM group compared with the sham group

Identification of RNA and LncRNA co-expression modules

Before constructing the network, the samples underwent clustering analysis to identify any notable outliers, revealing the absence of any outliers. First, the soft thresholding power β was chosen separately for the mRNA and lncRNA groups. For the mRNA group, β = 5 was selected, while for the lncRNA group, β = 8 was chosen. These values were selected to achieve a higher average connectivity degree (scale R² = 0.8) in both groups (Fig. 3A). The cluster dendrogram of all selected genes was generated using the adjacency matrix, resulting in the identification of co-expression modules. These modules are displayed in the analysis (Fig. 3B). The co-expression modules were observed to vary in size, ranging from small to large, reflecting the number of genes included. As a result, a total of 21 distinct modules were identified in the mRNA group, while 19 distinct modules were identified in the lncRNA group (Fig. 3C). In the mRNA group, the blue module had 6348 genes, the saddlebrown module had 4262 genes, and the brown module had 2648 genes, while in the lncRNA group, the blue and brown modules had 2480 and 860 genes, respectively (Fig. 3D).

Fig. 3

Construction and identification of weighted gene co-expression based on mRNA data and lncRNA data. A Analysis of network topology for various soft-threshold powers, including the scale-free fit index (y-axis) and the mean connectivity (degree, y-axis). B Hierarchical cluster dendrograms of all differentially expressed mRNAs or lncRNAs modules. C Hierarchical clustering heat map of all differentially expressed 21 mRNAs and 19 lncRNAs modules with different colors. D Number of genes within each module of mRNAs and lncRNAs. E Filter key modules related to treatment or point in time

The results from hierarchical clustering of sample groups reveal distinct expression patterns for mRNAs and lncRNAs. Specifically, the variation in mRNA expression appears to be more strongly influenced by treatment effects, with samples from the sham and dECM groups clustering into their respective branches, whereas samples are distributed across different branches based on time points within each group. In contrast, the expression changes in lncRNAs clustered into distinct groups based on either sham and dECM treatments or time points, suggesting that the effects of treatment and time on lncRNA expression are comparable and have not resulted in significant changes. This differential response highlights the unique regulatory mechanisms affecting mRNA and lncRNA expression under varying experimental conditions (Fig. S1A). To ascertain whether the identified modules are influenced predominantly by temporal dynamics or treatment effects, we applied Pearson correlation analysis to evaluate the relationships between module eigengenes and sequential time points. A correlation coefficient (Cor) greater than 0, significant at a p-value less than 0.05, signifies that the module eigengenes exhibit variation corresponding to the progression through time points, thus indicating a time-dependent nature. Conversely, for modules whose eigengenes show no significant variation with time, we conducted paired t-tests to determine differences between the sham and dECM groups at identical time points. A significant difference, denoted by a p-value less than 0.05, suggests that the variations in these modules are predominantly driven by treatment effects, rather than temporal changes. We discovered that the saddlebrown module (t-test:0.004206329) was specifically associated with the treatment process. Additionally, we found that the lightgreen (cor:0.883197516), steelblue (cor:-0.80204246), lightyellow (cor:0.78710782), and orangered4 (cor:-0.711325) modules were significantly linked to different time points. The blue (cor:0.874914822; t-test:3.28E-06) and brown (cor:0.78813544; t-test:0.04810722) modules were found to be associated not only with time but also with the treatment process. In lncRNAs, the grey60, greenyellow, and black modules were linked to time, whereas the blue and green modules were specifically related to the treatment process (Fig. S1B, Fig. 3E).

Functional analysis of the key modules

Selecting seven key modules, we plotted GO and KEGG enrichment bar charts for the top 5 enriched pathways (Fig. 4). In the saddlebrown module, we observed significant enrichment in the following biological process (BP) terms: “organic substance metabolic process” and “metabolic process”. In terms of cellular component (CC), the term “membrane-bounded organelle” was enriched. For molecular function (MF), the terms “protein binding” and “catalytic activity” were found to be enriched. Additionally, the KEGG pathway analysis revealed significant enrichment in the “metabolic pathways” (Fig. 4A). Moving on to the blue and brown modules (Fig. 4B). In the blue module, the terms “developmental process” and “nervous systerm development” were found to be significantly enriched in the BP category. Regarding CC, the term “intracellular part” was identified. The MF term “protein binding” was enriched. Furthermore, the blue module also exhibited enrichment in the “metabolic pathways” and “axon guidance”. In the brown module, the terms “positive regulation of biological process” and “regulation of response to stimulus” were found to be enriched in the BP category. In terms of MF and CC, the brown module displayed the similar with the blue module. As for time related modules (Fig. 4C), such as the lightgreen module, the the BP terms “mRNA processing” and “RNA processing” were found to be significantly enriched. In terms of CC, the term “ribonucleoprotein complex” was identified, in the MF terms “nucleic acid binding” and “RNA binding” were exhibited. Moreover, the “spliceosome” pathway was significantly enriched in the lightgreen module. In the steelblue, “immune response” and “phagocytosis, recognition” were enriched in BP terms, while “immunoglobulin receptor binding” and “extracellular space” were significantly enriched in MF and CC terms, respectively. The two most enriched BP, CC, MF terms and KEGG pathways for other modules are shown in Fig S2.

Fig. 4

Go and KEGG enrichment analyses for the 7 modules. A-C The top 5 enriched BP, CC, MF terms and KEGG pathways for each module are shown (with p-values shown as lines and bar lengths representing the number of genes corresponding to each pathway)

GO analysis of the 7 modules revealed several crucial biological processes, with the top 4 processes exhibiting significant P-values, while the time-course expression curves in these 7 modules were generated by plotting the average expression levels of the top 30 genes with the highest MM, and the gene expression pattern of seven modules were displayed (Fig. 5A-C), the average fold change of these 30 genes between the two groups were calculated (Fig. S3). Meanwhile, we plotted the expression levels of genes within 7 modules at eight time points between two groups (after normalization) (Fig. 5D).

The saddlebrown module exhibits clear correlations with treatment outcomes. Within this module, gene expression patterns are characteristically low in the sham group, whereas high expression levels are observed in the dECM group at all time points except day 1. The top 30 genes with the highest MM values within this module demonstrate an M-shaped expression trajectory over time, with slight decreases at 2 weeks and 8 weeks, but elevated expression at other time points. Prominent among these are genes such as Fbxw4, Tie1, Fbxo31, Scarf1, and Bahd1 (Fig. 5A). The blue and brown modules exhibit associations with both time and treatment efficacy. In blue, the top 30 genes with the highest MM values in this module show a gradual increase in expression levels over time, such as Arhgap32, Dock9, Tgfbr3, Olfml1, Ppara. Notably, the leading pathways in this module are connected to the development of the nervous system, highlighting their relevance to neural regeneration and recovery processes (Fig. 5B).

Expression of mRNA genes in time related modules were specific to regenerative time but were less associated with treatment. GO analysis of the lightgreen module revealed that the related genes (Ranbp1, Ebna1bp2, N6amt1, Cckar) showed significantly up regulated expression immediately after transplantation, and a small peak appeared at 2 w, then the expression remained low even at 12 w. The expression of genes (Igll1, Foxl2, Fscb) in steelblue module maintained low at the early stage, but rose slightly at 7 d, peaked at the 4 w, and then fell to the normal levels (Fig. 5C). Such a treatment or time related modules were also shown by lncRNAs (Fig. 6). In the identified grey60, greenyellow, and black modules, we observed a generally low level of lncRNA expression in the dECM group, whereas, in the sham group, expression levels peaked at 2 weeks, 3 weeks, and 4 weeks, respectively, as depicted in Fig. 6A. Conversely, the blue and green modules demonstrated associations with treatment effects. For the green module, lncRNA expression was low in the sham group across the board, but in the dECM group, there was a notable increment over time, particularly evident at 2 weeks, 3 weeks, and 12 weeks, as shown in Fig. 6B and further detailed in Fig. 6C.

Fig. 5

The discrete expression modules of mRNA expression by the WGCNA analysis. A Treatment related module, B Time & treatment related modules and C Time related modules. D The line diagrams with different colors represented the expression levels of genes within the 7 modules during the regeneration process between two groups (after homogenization)

Fig. 6

The discrete expression modules of lncRNA expression by the WGCNA analysis. A Eigengene bar plot of time related modules and B Treatment related modules of lncRNAs. C The line diagrams with different colors represented the expression levels of genes within the 5 modules during the regeneration process between two groups (after homogenization)

Regulatory networks between lncRNA-mRNA and immunohistochemical validation

We performed trans analysis on lncRNAs and mRNAs in modules related to time or treatment, and screened for positively correlated results (correlation greater than 0.95). The modules in lncRNAs exhibited the highest enrichment for KEGG pathways (Fig. 7A).

We analyzed the number of lncRNA-mRNA pairs between any two module pairs (5 lncRNAs and 7 mRNA modules). We found numbers of strong module-module co-expression pairs, which include pairs such as green (lncRNAs)-saddlebrown (mRNAs) and greenyellow (lncRNAs)-lightyellow (mRNAs). Interestingly, the blue (mRNAs) had strong multiple pairs, which include blue, green, black, grey60 and greenyellow modules in lncRNAs (Fig. 7B).

Remarkably, we observed a prevalence of positive co-expression pairs, aligning with the co-expression pattern observed in the majority of the genes under investigation. We performed an analysis to identify lncRNAs within each distinct module that exhibited negative correlations with mRNAs (Fig. 7C). Figure 7D showed that a substantial number of lncRNAs in blue module could negatively regulate the expression of both blue, brown, lightgreen and saddlebrown modules in mRNAs.

Fig. 7

Co-expression modules between lncRNAs and mRNAs. A Heat map presentation of the KEGG pathways for mRNAs associated with each lncRNA modules. B Circular presentation of module-module interaction between lncRNAs and mRNAs. The red and green line represented the trans targeting relationship between lncRNA and mRNA; mRNA and mRNA in modules, respectively. C Heat map presentation of functional clustering by the negatively paired mRNAs of each lncRNA module. D Circular presentation of association between lncRNA and mRNA modules that were negatively coregulated by lncRNAs

Finally, we selected 7 modules of mRNAs and 5 modules of lncRNAs, and merged the mRNAs-lncRNAs relationships among them. Then we selected the top 50 mRNAs with the most associations with lncRNAs, as well as the top 50 lncRNAs with the highest associations with mRNAs, and draw an interaction network diagram between them (Fig. 8A). For each of the seven mRNA modules identified as related to time or treatment, we selected the five most closely associated mRNAs. Correspondingly, we identified the top five associated lncRNAs for each lncRNA module linked to these mRNAs. Similarly, for each of the five lncRNA modules identified as pertinent to either time or treatment, we filtered out the top five associated lncRNAs. From the mRNA modules connected to these lncRNAs, the five most closely associated mRNAs were also selected (Fig. 8B).

Fig. 8

Co-expression network illustration between lncRNAs and mRNAs. Constructions of the A Top50 and B Top 5 lncRNA-mRNA network. mRNAs were represented by circles, while lncRNAs were represented by triangles. The color of the circles corresponded to the color of the module. The thickness of the connecting lines between the circles and triangles represented the strength of their correlation (absolute value). The mRNA associated with the top 5 of each module and the lncRNA associated with the top 5 of each module were marked with red borders around them

Following 12 weeks of bridging sciatic nerve defects with cell-matrix nerve grafts, significant positive expression of ATP1A3, LGALS3, and HMOX1 were observed in the regenerated sciatic nerve tissues. Notably, ATP1A3 and LGALS3 demonstrated elevated expression levels predominantly in axons, whereas Hmox1 exhibited heightened expression in Schwann cells. These findings underscore the differential expression patterns of these genes in specific cellular components within the nerve regeneration environment, highlighting their potential roles in mediating neural repair processes (Fig. 9A-C). We instigated the qRT-PCR targeting the expression levels of five key genes implicated in the regulatory process of nerve regeneration, including Fgf4, Syt11, Lhx2, Tnr and IL6 (Fig. 9D). These findings confirm the accuracy and reliability of our data.

Fig. 9

Validation of key molecules. A-C After 12 weeks, longitudinal sections of the regenerated sciatic nerve bridging segment post cell-matrixed nerve grafts implantation were subjected to immunofluorescent staining of anti-ATP1A3, anti-LGALS3 (red) merged with anti-NF-200 (green) and DAPI (blue), while anti-Hmox1(green) merged with anti-S100 (red) and DAPI (blue). Scale bar, 20 μm (high magnification). D Line charts showed the qRT-PCR of relative mRNA expressions of Fgf4, Syt11, Lhx2, Tnr, and IL6 in the dECM group, which were normalized to Gapdh. Each point represents the average measurement of three biological replicates. Bars correspond to standard error of mean ( mean ± SEM)

Discussion

Our previous research successfully demonstrated the construction of a cell-matrixed nerve graft, which created a favorable microenvironment for nerve regeneration. Furthermore, we elucidated the unique molecular regulation patterns and characteristics associated with different modes of nerve repair. These findings, significantly contribute to our understanding of the intricate mechanisms underlying nerve regeneration and provide valuable insights for the development of novel therapeutic strategies in this field.

However, molecular pattern recognition has emerged as a highly effective strategy in scientific research, catalyzing a paradigm shift in the traditional medical concept from a single-parameter model to a more comprehensive multi-parameter systematic mode [27]. In contrast to the high degree of complexity and dynamic nature of proteoforms within a proteome, RNAs encompassed within a transcriptome exhibit a comparatively simpler and more stable profile [28]. The RNA repertoire comprises various types, including mRNAs, lncRNAs, and miRNAs, each playing distinct roles in gene regulation and cellular processes [29]. In this study, we employed the WGCNA method to identify 21 co-expression modules comprising 18,609 mRNAs and 19 co-expression modules comprising 5,708 lncRNAs from a comprehensive analysis of 48 bridged sciatic nerve samples obtained from the dECM group and sham group. This approach, which prioritized the correlations between traits and co-expression modules, yielded results of exceptional biological significance and reliability [30]. Our analysis revealed disparities in the number of differentially expressed genes between the sham and dECM groups, with these discrepancies being sufficiently pronounced at identical time points. For instance, at comparative time intervals such as 1 day versus 3 days, the dECM group exhibited twice as many differentially expressed genes as the sham group. This significant difference underscores a more pronounced effect of the dECM treatment, indicative of its potent impact in eliciting a broader spectrum of gene expression alterations. We performed pairwise correlation analysis and paired t-tests on the characteristic values of each module at different time points in two groups. This allowed us to identify seven modules in mRNAs and five modules in lncRNAs that were significantly associated with progressing or time points.

While mRNA co-expression networks have been described as important in understanding the nerve regeneration, very few of them appear to reflect the complexity of regeneration process [31, 32]. By utilizing this method, we were able to unravel the complex network of mRNA-mRNA and lncRNA-mRNA interactions, thereby gaining valuable insights into the comprehensive molecular mechanisms underlying the changes induced by dECM in the sciatic nerve. We demonstrated how the data set can be used to profile trajectories of genes associated with specific neurobiological categories or disorders, many of which are not likely to be evident from transcriptomic profiles of commonly studied model systems. Coupled with analysis of co-expressed genes in the data set, these mRNA co-expression networks provide information on specific timing and treating of various genes expressed in the dECM group, which will also offer insights regarding their function. Additionally, we employed pathway enrichment analysis and gene co-expression network analysis to delve deeper into the complex regulatory mechanisms involving genes and lncRNAs that govern these processes. This integrative approach facilitates a comprehensive understanding of the intricate molecular interactions that underpin biological functions. Scarf1 is predominantly expressed in macrophages and endothelial cells, playing a crucial role in the clearance of bodily waste and pathogens. It recognizes and binds to a variety of ligands, including bacteria, viruses, and apoptotic cells. By mediating the endocytosis and degradation of these ligands, Scarf1 exerts a significant influence on the maintenance of tissue homeostasis [33]. Dock9 is identified as a Rho family guanine nucleotide exchange factor (RhoGEF). Its primary function lies in modulating the activity of Rho family GTPases within the cell. Dock9 is integral to the regulation of critical biological processes including cytoskeletal reorganization, cell migration, cell adhesion, and cell polarity, thereby playing a vital role in maintaining cellular structure and function [34]. LncRNAs play a pivotal role in the regeneration of peripheral nerves [17]. Transcriptome-wide studies have revealed the temporal regulation of lncRNA expression in dorsal root ganglia (DRG) during peripheral nerve regeneration, highlighting the association between the downregulation of lncRNA BC089918 and promoted neurite outgrowth in DRG neurons [35]. Beyond DRG neurons, lncRNAs are also linked with phenotypic changes in Schwann cells. A recent identification of lncRNA Loc680254 and BC088259, which is crucial for Schwann cell proliferation and migration [36, 37]. In the co-expression network depicted in Fig. 8, lncRNA Kcnq1 was found to primarily target Tacr2. LncRNA Kcnq1 influences gene expression and various cellular functions, such as proliferation, migration, epithelial-mesenchymal transition, apoptosis, viability, autophagy, and inflammation [17, 38]. It can also suppress microglial activation, reduce the production of pro-inflammatory cytokines, and alleviate neuropathic pain [39]. Tacr2, a G-protein-coupled receptor, primarily functions in the central and peripheral nervous system by binding to the Neurokinin A, which regulates various physiological responses [40, 41]. The potential effects of lncRNA Kcnq1 on peripheral nerve growth and functional recovery warrant further investigation to deepen our understanding of these molecular interactions.

In view of the significance of the blue (mRNAs) module, we explored the interactions among genes within the module. Notably, the hub genes in the blue module, which are involved in anatomical structure morphogenesis and nervous system development, include Septin8, Fgf4, Syt11, Arhgef28, and Ankrd50 (Fig. S4A). A study employed an osmotic ex vivo axoplasm isolation method, using label-free quantitative proteomics to analyze changes in the axoplasm of adult rat sciatic nerves before and after injury. The analysis identified proteins such as netrin- 1, the metalloprotease TIMP-2, calmodulin, nestin, and septin proteins 2, 5, 7, and 9, all of which are associated with axonal response, axonal guidance, and plasticity processes [42]. Septins, characterized by their GTP-binding and membrane-interacting properties, possess highly conserved domain structures and participate in various cellular processes including cytoskeletal organization, cytokinesis, and membrane dynamics. To date, 13 distinct septin genes (SEPT1 to SEPT12 and SEPT14) have been identified in mammals [43]. In sensory neurons, SEPT6 and SEPT7 accumulate at early sites of filopodia formation in axons. Schwann cells lacking Sept7 fail to ensheath DRG neuron axons and exhibit disorganized actin cytoskeletons [44]. Within the blue module of our study, we noted a high expression of Septin 8. As a member of this group, Septin 8 plays crucial roles in regulating cytoskeletal reorganization, maintaining cellular polarity, and mediating intercellular interactions. It also significantly contributes to morphogenesis, signal transduction, and the dynamic of the cytoskeleton in neuronal cells. This aligns with the primary pathways involved in nervous system development identified in the blue module, highlighting the importance of Septin 8 in these processes. ClueGO/CluePedia network analysis was performed on the blue module, selecting the pathway with its enrichment result of top5 and screening the genes involved (Fig. S4B). Such as the neurotrophic factor Sema3d, which can promote the growth and guidance of neuronal axons and participate in the migration and positioning processes of neuronal cells [45, 46]. To validate whether key genes were translated into proteins to exert their biological functions in peripheral nerve regeneration, we performed immunofluorescent staining at the tissue cellular level [11, 47, 48]. ATP1A3, LGALS3, and HMOX1 were validated for their important roles in peripheral nerve regeneration. ATP1A3 plays a critical role in nervous system and early brain development [49, 50]. LGALS3 plays an important role in oligodendrocyte differentiation, myelin functionality [51], and may influence neural regeneration by affecting inflammatory responses [52]. HMOX1 participates in signaling pathways during the acute phase of neural injury, serving as an anti-inflammatory agent [53], and playing a significant role in related metabolic diseases [54].

In our study, we provided a systematic and comprehensive analysis of the gene expression profile related to the repair of peripheral nerve defects using cell-matrixed nerve grafts. This analysis revealed patterns and molecular regulatory mechanisms in terms of the interaction between materials and organisms, as well as between tissue cells and materials. It represents a direction for the development of biomimetic tissue engineering technology and clinically effective products.

Conclusion

Our data enhance genome-wide associations and linkage studies by narrowing the focus to any candidate genes that are specifically expressed during nerve regeneration in the dECM group. We report here on the dynamic changes observed in mRNA co-expression networks that may serve as a regulatory system that truly contributes to the complexity of nerve regeneration by the cell-matrix nerve graft.