Runzhi Huang1,2,3†
Siqiao Wang3†
Rui Zhu1,2,3†Shuyuan Xian3
Zongqiang Huang4*
Liming Cheng1,2*
Jie Zhang1,2*- 1Key Laboratory of Spine and Spinal Cord Injury Repair and Regeneration (Tongji University), Ministry of Education, Shanghai, China
- 2Division of Spine Surgery, Department of Orthopedics, Tongji Hospital, Tongji University School of Medicine, Shanghai, China
- 3Tongji University School of Medicine, Shanghai, China
- 4Department of Orthopedics, The First Affiliated Hospital of Zhengzhou University, Zhengzhou, China
Background: Spinal cord injury (SCI) is a severe neurological deficit affecting both young and older people worldwide. The potential role of key enhancer RNAs (eRNAs) in SCI remains elusive, which is a prominent challenge in the trauma repair process. This study aims to investigate the roles of key eRNAs, transcription factors (TFs), signaling pathways, and small-molecule inhibitors in SCI using multi-omics bioinformatics analysis.
Methods: Microarray data of peripheral blood mononuclear cell (PBMC) samples from 27 healthy volunteers and 25 chronic-phase SCI patients were retrieved from the Gene Expression Omnibus database. Differentially expressed transcription factors (DETFs), differentially expressed enhancer RNAs (DEeRNAs), and differentially expressed target genes (DETGs) were identified using the Linear Models for Microarray Data (limma) package. Fraction of immune cells was estimated using CIBERSORT algorithm. Gene Set Variation Analysis (GSVA) was applied to identify the downstream signaling pathways. The eRNA regulatory network was constructed based on the correlation results. Connectivity Map (CMap) database was used to find potential drugs for SCI patients. The cellular communication analysis was performed to explore the molecular regulation mechanism of SCI based on single-cell RNA sequencing (scRNA-seq) data. Chromatin immunoprecipitation sequencing (ChIP-seq) and Assay for Transposase-Accessible Chromatin using sequencing (ATAC-seq) data were used to validate the key regulatory mechanisms. scRNA-seq dataset was used to validate the cell subtype localization of the key eRNAs.
Results: In total, 21 DETFs, 24 DEeRNAs, and 829 DETGs were identified. A regulatory network of 13 DETFs, six DEeRNAs, seven DETGs, two hallmark pathways, two immune cells, and six immune pathways was constructed. The link of Splicing factor proline and glutamine rich (SFPQ) (TF) and vesicular overexpressed in cancer prosurvival protein 1 (VOPP1) (eRNA) (R = 0.990, p < 0.001, positive), VOPP1 (eRNA) and epidermal growth factor receptor (EGFR) (target gene) (R = 0.974, p < 0.001, positive), VOPP1, and T helper (Th) cells (R = −0.987, p < 0.001, negative), and VOPP1 and hallmark coagulation (R = 0.937, p < 0.001, positive) was selected. Trichostatin A was considered the best compound target to SCI-related eRNAs (specificity = 0.471, p < 0.001).
Conclusion: VOPP1, upregulated by SFPQ, strengthened the transient expression of EGFR. Th cells and coagulation were the potential downstream pathways of VOPP1. This regulatory network and potential inhibitors provide novel diagnostic biomarkers and therapeutic targets for SCI.
Introduction
Spinal cord injury (SCI) refers to functional or structural injury of the spinal cord causing total or partial loss of motor, sensory, and sphincter function below the injured segment (Nakae et al., 2011; Anderson et al., 2018). Whether the pathogenesis is disease or trauma, SCI exhibits high disability rates. It affects approximately 347,000 individuals in the United States, with about 17,500 new cases diagnosed every year (Badhiwala et al., 2019; GBD 2016 Neurology Collaborators, 2019). According to the illness course, SCI can be divided in three phases: the acute phase (0–15 days), the subacute phase (3–5 months), and the chronic phase (6–12 months) (Pouw et al., 2011). The most common pathological characteristics of SCI are a cascade of molecular and cellular events triggered by inflammation, and excitotoxicity impairs endogenous regeneration, namely, axonal outgrowth and remyelination.
It is arduous to repair injured neurons and restore conducting function of axons, so treatments of SCI have become worldwide problems (Borton et al., 2014). Importantly, there are no efficacious drugs or therapeutic approaches for SCI. Thus, many patients suffer substantial physical and psychological consequences. To date, the molecular mechanisms of SCI remain unclear, so it is difficult to develop novel drugs or treatments. Thus, it is urgent to determine specific molecular mechanisms that underlie the pathogenesis of SCI.
The gene mutation and aberrant expression are implicated in the development of numerous diseases and may involve various types of regulatory molecules. Enhancer refers to a kind of distal regulatory DNA element that is able to enhance the transcription of corresponding target genes via coordinating with target gene promoters (Blackwood and Kadonaga, 1998). Widely acknowledged as DNA elements that nucleate transcription factor (TF) binding, enhancers were identified to also transcribe non-coding RNAs recently, which are referred to as enhancer RNAs (eRNAs) (Kim et al., 2010). Numerous eRNAs have been identified in Homo sapiens so far, and many of them were suggested to play crucial parts in the transcriptional circuitry (Li et al., 2016).
Inflammation after SCI is a fundamental basis of secondary damage. Inflammatory response is coordinated by various signaling modalities that include the epigenetic modification of promoters and eRNAs (Rudman et al., 2018). Large amounts of cytokines are quickly generated and released from various cells in the damaged structures after SCI (Hayashi et al., 2000; Pineau and Lacroix, 2007). Cytokines activate apoptosis of pathological nerve cells and recruitment of active leukocytes in myelopathic lesions (Beck et al., 2010; Chen et al., 2011). Leukocyte infiltration contributes to cascaded amplification of inflammation cytokine signaling, increasing neurotoxicity, and thus promoting fibrotic scar formation via recruiting fibroblasts in damaged regions (Zhu et al., 2015). In addition, damaged areas chronically maintain pro-inflammatory phenotype after injury, hindering regeneration and healing of the injured spinal cord (Beck et al., 2010). Because of the detrimental inflammatory effects after SCI, reducing inflammatory response is a major demand for improving clinical outcomes of SCI patients. Nevertheless, eRNA-targeted anti-inflammation agents that effectively inhibit inflammation and improve the prognosis of SCI patients are still insufficient, and novel anti-inflammation interventions are urgently needed.
In this study, transcription factors (TFs), differentially expressed enhancer RNAs (DEeRNAs), and target genes in patients with SCI were identified. Significant immune cells and immune-related pathways were identified by cell type identification by estimating relative subsets of RNA transcripts (CIBERSORT) and single-sample Gene Set Enrichment Analysis (ssGSEA) algorithms, respectively. In addition, TFs, eRNAs, target genes, immune cells, immune-related gene sets, hallmark gene sets (signaling pathways) acquired by the Gene Set Variation Analysis (GSVA) were merged in correlation analysis. Only remarkable interactions were extracted for construction of eRNA regulatory network. Moreover, we performed Connectivity Map (CMap) analysis to explore the potential inhibitors specific to SCI-related eRNAs. Furthermore, chromatin immunoprecipitation sequencing (ChIP-seq) and Assay for Transposase-Accessible Chromatin using Sequencing (ATAC-seq) data were utilized for validation of the key regulation mechanisms in this network.
Materials and Methods
Data Acquisition
The present study was endorsed by the Ethics Committee of the Shanghai Tongji Hospital affiliated to Tongji University. Microarray data and clinical information of 27 peripheral blood mononuclear cell (PBMC) samples from healthy volunteers and 25 PBMC samples from SCI patients were obtained from the Gene Expression Omnibus (GEO) database (accession number: GSE82152)1 (Barrett et al., 2011) and E-GEOD-69901 (ArrayExpress)2 (Athar et al., 2019), respectively. In addition, these two batches of microarray data were both retrieved from a platform called GPL21975 [PrimeView] Affymetrix Human Gene Expression Array [Brainarray ENTREZG Version 17],3 compared with data merging based on multiple different platforms, and resulted in less error. TFs were obtained from the Cistrome Cancer database4 (Zheng et al., 2019). Immunologically related genomic expression profiles were downloaded from the ImmPort database5 (Bhattacharya et al., 2014). Hallmark signaling pathways were collected from the Molecular Signatures Database (MSigDB, v7.2)6 (Liberzon et al., 2015). ChIP-seq data of H3K27ac (accession number: GSE134744) and ATAC-seq data of key eRNAs (accession number: GSE139099) were obtained from the GEO database (see text foot note 1) (Barrett et al., 2011).
Data Preprocessing and Differential Expression Analysis
Samples and patients with incomplete clinical information were excluded. All original microarray data were read using affy package, followed by robust multi-array average (RMA) background correction, standardization, probe-specific background correction, and summarizing probe set values in one expression measure (Irizarry et al., 2003; Gautier et al., 2004). Then, two batches of microarray data aforementioned were corrected using normalizeBetweenArrays function in Linear Models for Microarray Analysis (limma) package, which were then merged for differential expression analysis. Differential expression analysis of eRNAs (DEeRNAs), target genes (DETGs) and TFs (DETFs) between SCI and normal blood was conducted using the Linear Models for Microarray Data (limma) package (Smyth, 2004). For p-value, false discovery rate (FDR) was utilized for multiple testing correction. Absolute log2[Fold Change (FC)] ≥ 1.0 and the FDR < 0.05 were cutoff criteria.
Functional Enrichment Analysis
We conducted Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses to identify biological processes and pathways that were most related to DEeRNAs with FDR p-value < 0.05 as cutoff value (Zhang T. et al., 2013).
Identification of Potential Immune Cells, Immune Pathways, and Hallmark Pathways
Immune cell proportions in the SCI and normal samples were analyzed utilizing the cell type identification by estimating relative subsets of RNA transcripts (CIBERSORT) algorithm (Newman et al., 2015). To identify the association between eRNA signature and immune cell infiltration in SCI tissues, we uploaded gene expression matrix data to CIBERSORT database7 to purify cellular subtype-specific gene expression. Infiltrating immune cells were all extracted for further analysis. Furthermore, non-parametric tests were used to identify correlations between these immune cells/pathways and various clinical phenotypes.
In addition, ssGSEA algorithm was conducted to evaluate and quantify the enrichment level of 10 immune-related pathways in SCI samples (Xiao et al., 2020). Furthermore, Pearson correlation analysis was carried out to determine the correlations between key eRNAs and immune-related signaling pathways or immune genes, in which a p-value < 0.05 was considered statistically significant.
Identification of Potential Downstream Hallmark Pathways
GSVA (Ferreira et al., 2021) and GSEA (Subramanian et al., 2007) were both utilized for exploring potential downstream hallmark pathways of DEeRNAs. Absolute quantification of 50 hallmark signaling pathways was evaluated to extract differentially expressed hallmark pathways between SCI and normal blood using ClusterProfiler package and GSVA package (Yu et al., 2012; Hanzelmann et al., 2013). GSEA was carried out to reveal the significant enrichment of upregulated and downregulated hallmark pathways in SCI and normal blood. Furthermore, correlations of hallmark pathways of GSVA and GSEA were extracted, and the interactional pathways were suggested as key pathways.
Construction of Regulation Network for Transcription Factors, Enhancer RNAs, Immune Cells, Immune-Related Pathways, Target Genes, and Hallmark Gene Sets
Pearson correlation analysis was conducted based on key TFs, eRNAs, immune genes, immune cells, target genes, and hallmark pathways aforementioned. An eRNA regulatory network was then reconstructed by Cytoscape (3.7.1) (Kohl et al., 2011).
Molecular complex detection (MCODE) plugin of Cytoscape software was utilized to identify significant molecules within this network in a recognition standard MCODE score ≥ 4 to extract the modules of hub genes. Interaction relationships between key eRNAs and other components were controlled based on p-value < 0.05 and | correlation coefficient| > 0.40.
Connectivity Map Analysis
Here, we used CMap (build 02) to find potential inhibitors that may target SCI-related eRNAs. In total, 6,100 gene expression cases covering 1,309 drugs were obtained from the CMap database8 (Lamb et al., 2006). That is, a candidate drug might correspond to various gene expression cases. Genes in each case were ranked via taking differential expression values between drug-untreated and drug-treated cell lines, and 6,100 gene lists related to drugs were then generated. Based on identified key DEeRNAs involved in SCI and 6,100 drug-related cases, we conducted a non-parametric test to explore the relationship between drugs and SCI.
Information on targeting compounds is available in the mechanism of actions (MoA)9 that includes transcriptional responses of various human cell lines to perturbagens, structural formulas, and protein targets. On the basis of MoA, compounds that may target SCI-related eRNAs/enhancers in this study were extracted.
Single-Cell RNA Sequencing Data Processing
Following the procedure of 10x Genomics Chromium10 (Chen et al., 2020), the preprocessing of samples and scRNA-seq data was conducted. After demultiplexing, the sequencing results were divided into two pair-ended reads fastq files that were then trimmed to eliminate template switch oligo (TSO) sequence and polyA tail sequence. Additionally, clean reads were aligned with the GRCh38 (Version: 100) genome assembly, which were quantified using the Cell Ranger Software (Version 1.0.0).11
The quantitative gene expression matrices (The row names of matrices were genes and column names were barcodes) acquired from seven libraries that included seven vertical section samples and 14 cross-section samples were analyzed with Seurat pipeline (Version: 3.2.2) for further analysis (Butler et al., 2018). Only cells with less than 10% mitochondrial gene mapped and more than 100,000 transcripts expressing were extracted for subsequent analysis. In addition, genes that expressed in over three single cells were included in follow-up analyses. After the completion of quality control (QC), all samples were merged in one Seurat object with the function of “IntegrateData” that were then scaled and standardized with the function of “ScaleData.” Top 1,500 variable genes were identified using “vst” method. To reduce model dimensionality, principal component (PC) analysis (PCA) was initially carried out, and the top 20 PCs were incorporated as input file for Uniform Manifold Approximation and Projection (UMAP) for dimension reduction analysis. The UMAP plots that illustrated cell subclusters were constructed using the “DimPlot” and “RunUMAP” function.
Differential Expression Analysis of Single-Cell RNA Sequencing
Genes with remarkable differential expression from the top 1,500 variable genes were defined as differentially expressed genes (DEGs) with “wilcox” method using “FindAllMarkers” function.
Cell Type Annotation
To identify the cell type of each unsupervised cluster, DEGs of all subclusters were utilized as potential references that were combined with known specific cell surface biomarkers obtained from CellMarker12 (Zhang et al., 2019) for a comprehensive annotation of cell type. Given the variable gene expression patterns in SCI, a specific cellular annotation method was utilized in the present study. Firstly, known cell surface biomarkers of neuron (Rph3a, Tubb3, Gnal), neuron precursor cell (NPC) (Sox2, Prom1, Sox9), astrocyte (Gfap, S100b, Vim), fibroblast (Col1a1, Col1a2, Col4a1), oligodendrocyte (Mbp, Olig2, Gjb1), oligodendrocyte precursor cell (OPC) (Olig1, Apoc4, Epn2), microglia (Cx3cr1, Aif1, Itgam), and macrophage (Cd68, Itgam, Fcgr3) were used to annotate the cell type of the eight cells aforementioned. Furthermore, because several biomarkers of macrophage and microglia were overlapping, cells with transcripts per million (TPM) of Cd68 < 2 and Aif1 > 1.38 were identified as microglia, while cells with TPM of Cd68 > 2 were identified as macrophages. In addition, the cellular feature plots, dot plots, violin plots, and heat maps were constructed to show the marker genes of each cell type using SCANPY (Version: 1.7.1) and Seurat R package (Version: 3.2.2) under the environment of Python 3.6 (Butler et al., 2018; Wolf et al., 2018).
Cellular Communication Analysis
In order to elucidate the significant cellular communication patterns and ligand–receptor pairs among various different cell types in the spinal cord, cellular communication analysis was carried out using iTALK R package (Version: 0.1.0)13 (Wang et al., 2019). Firstly, because of the scarcity of musculus resources in present mainstream cellular communication algorithms, the top 1,500 variable genes were transformed to human genes using biomaRt package (Version: 2.46.0) in a certain homologous degree with “getLDS” function (Durinck et al., 2009). Secondly, the normalized expression matrix of these genes was incorporated into the iTALK object with the “rawParse” function. Finally, the top 200 ligand–receptor pairs were shown by ligand–receptor plots and iTALK networks.
Chromatin Immunoprecipitation Sequencing Validation
Histone H3K27ac was implicated in enhancer-specific modifications, which were essential for enhancers to activate the transcription of relevant target genes (Kang et al., 2021). Role of H3K27ac in eRNA transcription was evaluated by analyzing ChIP-seq data (accession number: GSE134744) in peripheral blood (Grubert et al., 2020).
The ChIP-seq data of Splicing factor proline and glutamine rich (SFPQ) (accession number: GSM2827312, GSM1411215, GSM2825596, GSM1097497, GSM2827311, and GSM1097496) (Dunham et al., 2012; ENCODE Project Consortium, 2012; West et al., 2014) and H3K27 (accession number: GSM732912, GSM575294, and GSM663427) (Lister et al., 2009; Guenther et al., 2010; Wang et al., 2011) were obtained from the Cistrome database for validating binding relationships between key eRNAs and other significant biomarkers in this study. Furthermore, we determined the eRNA binding relationships using the University of California Santa Cruz (UCSC) Genome Browser14 based on the original ChIP-seq data (Robinson et al., 2011; Li et al., 2015).
Assay for Transposase-Accessible Chromatin Using Sequencing Validation
The ATAC-seq (Buenrostro et al., 2013) refers to an impressively flexible, simple, and powerful technique to profile chromatin regions genome-wide compared to traditional methods like functional assays or sequence conservation analyses.
We downloaded peak data of key TFs and eRNAs (accession number: GSE139099) (Grubert et al., 2020) and performed the correlation analysis using ggpubr.15 In order to explore the potential super-enhancers, we performed the peak calling with macs14 (Feng et al., 2012).16 Furthermore, we ranked enhancers and determined super-enhancers using ROSE (Whyte et al., 2013).17 Multiple SCI-enriched ATAC-seq peaks that represented candidate regulation elements were located near established eRNAs and were enriched in distinct sets of TF binding sites.
Online Single-Cell RNA Sequencing Validation
To identify the expression of the key biomarkers in single-cell level, scRNA-seq data of developing mouse spinal cord (accession number: GO0006836 and GO0098609) were downloaded from Single Cell Expression Atlas (Delile et al., 2019).18 Then, we validated the cellular location of these biomarkers’ expression.
Statistics Analysis
All statistical analyses were put into effect using R version 3.6.1 (Institute for Statistics and Mathematics, Vienna, Austria).19 In descriptive statistics, mean ± standard deviation was utilized for continuous variables in normal distribution.
In addition, when encountering continuous variables in abnormal distribution, the median was utilized. Two-sided p-value < 0.05 was suggested to be necessary for statistics.
Results
Identification of Differentially Expressed Enhancer RNAs
Analysis process of the present study was illustrated in Figure 1. Samples and patients with incomplete clinical information were excluded, and conformers were shown in Supplementary Table 1. Based on the threshold, a total of 3,979 eRNAs were identified as DEeRNAs between 27 healthy volunteers and 25 SCI patients from 5,100 eRNAs, which were shown in the Heat map (Figure 2A) and volcano plot (Figure 2B).
Figure 1. The analysis flowchart of all analysis processes.
Figure 2. Identification of differentially expressed enhancer RNAs (DEeRNAs) and functional annotation. (A) Heat map showing the expression level of DEeRNAs. (B) Volcano plot showing the p-value and log| FC| value of DEeRNAs by the differential expression analysis above. (C) The output of Gene Ontology (GO) analysis. (D) The output of Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enrichment analysis.
Functional Enrichment Analysis
GO and KEGG enrichment analyses were both conducted to explore the potential mechanism of identified DEeRNAs. In GO analysis, the most important terms in biological process (BP), cellular component (CC), and molecular function (MF) were neutrophil activation, endosomal part, and cell adhesion molecule binding (Figure 2C). The most significant KEGG pathway was mitogen-activated protein kinase (MAPK) signaling pathway (Figure 2D). In the heat map (Figure 3A) and volcano plot (Figure 3B), 24 differentially expressed SCI-related eRNAs were identified, which were considered to be key eRNAs in patients with SCI.
Figure 3. Identification of differentially expressed spinal cord injury (SCI)-related enhancer RNAs (eRNAs) and transcription factors (TFs). (A) Heat map showing the expression level of differentially expressed SCI-related eRNAs. (B) Volcano plot showing the p-value and log| FC| value of differentially expressed SCI-related eRNAs. (C) Heat map showing the expression level of differentially expressed TFs. (D) Volcano plot showing the p-value and log| FC| value of differentially expressed TFs.
Correlation Analysis of Key Differentially Expressed Enhancer RNAs and Differentially Expressed Transcription Factors
Twenty-one DETFs (log2 FC > 1 or < −1 and FDR < 0.05) were identified using limma package that were shown in the heat map (Figure 3C) and volcano plot (Figure 3D). In addition, to explore the relationship between DETFs and DEeRNAs, Pearson correlation analysis was carried out, and 27 regulatory relationships were identified (correlation coefficient < −0.300 or > 0.300, and p-value < 0.001). Based on the correlation analysis results for regulation relationships of key DETFs and DEeRNAs, link of SFPQ (TF) and vesicular overexpressed in cancer prosurvival protein 1 (VOPP1, eRNA) was extracted (R = 0.990, p < 0.001, positive).
Correlation Analysis of Key Differentially Expressed Enhancer RNAs and Differentially Expressed Target Genes
In total, 829 DETGs (log2 FC > 1 or < −1 and FDR < 0.05) were identified using limma package and illustrated in the heat map (Figure 4A) and volcano plot (Figure 4B). Furthermore, to determine the relationships between the identified DETGs and DEeRNAs, eight regulation relationships were identified using Pearson correlation analysis (correlation coefficient < −0.300 or > 0.300 and p-value < 0.001). Additionally, based on significant regulation relationships of key DEeRNAs and DETGs acquired by Pearson analysis, the link of VOPP1 (eRNA) and epidermal growth factor receptor (EGFR, target gene) was extracted (R = 0.974, p < 0.001, positive).
Figure 4. Identification of differentially expressed target genes (DETGs) and evaluation of immune cells/gene sets. (A) Heat map showing the expression level of DETGs. (B) Volcano plot showing the p-value and log| FC| value of DETGs. (C) Bar chart showing the percentage of 22 kinds of immune cells from healthy volunteers and spinal cord injury (SCI) samples by cell type identification by estimating relative subsets of RNA transcripts (CIBERSORT). (D) Heat map showing the expression level of 22 immune cells by CIBERSORT. (E) Box plot showing the fraction of 22 immune cells from healthy volunteers and SCI samples by non-parametric tests. (F) Heat map showing the expression level of 28 immune cells/pathways via single-sample Gene Set Enrichment Analysis (ssGSEA).
Correlation Analysis of Key Differentially Expressed Enhancer RNAs, Immune Cells/Pathways
Twenty-two types of immune cells/pathways were identified using CIBERSORT. The bar chart showed the percentage of 22 kinds of immune cells in 27 healthy volunteers and 25 SCI patients (Figure 4C). The heat map showed the expression level of 22 immune cells by CIBERSORT (Figure 4D). Results of non-parametric tests demonstrated significant correlation between immune cells/pathways and SCI (Figure 4E). Furthermore, ssGSEA was conducted to identify 28 immune cells/pathways that were significantly correlated with differentially expressed SCI-related eRNAs. Results of ssGSEA were shown in the heat map (Figure 4F).
Then, by co-analyzing VOPP1 and 22 types of immune cells or pathways, the top 3 immune cells or immune-related pathways were extracted for further analysis, which included Antigen-presenting cells (APC) co-inhibition (R = 0.936, p < 0.001, positive), MHC class I (R = −0.973, p < 0.001, negative), and T helper (Th) cells (R = −0.987, p < 0.001, negative). Th cells were finally extracted in further analyses.
Identification of Potential Downstream Hallmark Pathways of Key Differentially Expressed Enhancer RNAs
Results of GSVA (| log2 FC| > 0.1, p < 0.05) were shown in the heat map (Figure 5A). The volcano plot showed 28 differentially expressed hallmark pathways between 27 healthy volunteers and 25 SCI patients (Figure 5B). The t-value of GSVA score of these pathways was shown in the t-test bar chart (Figure 5C). Furthermore, 38 hallmark pathways were extracted from the 50 hallmark pathways for further analysis, controlled by cutoff values of | log2 FC| > 0.1 and p < 0.05. Importantly, VOPP1 and hallmark coagulation (R = 0.937, p < 0.001, positive) showed the strongest interaction.
Figure 5. Identification of key hallmark pathways. (A) Heat map showing the co-expression level of hallmark sets in spinal cord injury (SCI) samples by Gene Set Variation Analysis (GSVA). (B) Volcano plot showing the co-expression level of hallmark sets in SCI samples by GSVA. (C) The bar plot revealing the t-value of GSVA score. (D) Overview of the protein–protein interaction network of OA-related enhancer RNAs (eRNAs), translational factors (TFs), immune cells/pathways, and hallmark gene sets. The diamonds represented key eRNAs. The arrows represented TFs. The ellipses represented immune cells. The triangles represented immune pathways. The hexagons represented target genes. The rectangle represented hallmark gene sets.
Construction of the Enhancer RNA Regulation Network
A co-expression regulation network of key DETFs, DEeRNAs, DETGs, immune cells/pathways, and hallmark pathways was constructed, which elaborated on the regulatory relationships among the aforementioned components (Figure 5D). Furthermore, to quantify the interaction coefficients among them, co-expression analysis was performed at the transcriptional level (Figure 6A).
Figure 6. Correlation analysis of the key biomarkers and identification of potential small-molecule drugs for spinal cord injury (SCI). (A) Heat map showing the results of correlation analysis (Pearson analysis) of enhancer RNAs (eRNAs), transcription factors (TFs), target genes, immune cells/pathways, and hallmark pathways. (B) Heat map showing perturbagens from the Connectivity Map (CMap) that might be capable of targeting SCI-related eRNAs. (C) Structural formula of trichostatin A. (D) Structural formula of 15-delta prostaglandin J2.
Connectivity Map Analysis
Because eRNAs were implicated in the pathological processes of SCI and traditional long-term treatment with drugs may result in severe side effects and/or insufficient inflammation and pain relief, it is urgent to find potential compounds that target SCI-related eRNAs. Through exploring potential compounds in CMap database, we identified 19 compounds that have been validated to target SCI-related eRNAs in multiple clinical trials, including 15-delta prostaglandin J2, anisomycin, azaperone, bumetanide, dosulepin, felodipine, lycorine, metrizamide, MG-262, mitoxantrone, naftifine, pirinixic acid, Prestwick-864, probenecid, puromycin, rifampicin, trichostatin A (TSA), troglitazone, and valinomycin (Figure 6B). TSA (specificity = 0.471, p < 0.001) and 15-delta prostaglandin J2 (specificity = 0.128, p < 0.001) with the highest specificity were considered the best compounds to target SCI-related eRNAs (Figures 6C,D), and TSA was extracted for further analysis.
Integrated Analysis of Single-Cell RNA Sequencing
The UMAP scatter plots and cellular feature plots showed the cell types and marker genes of different cell types (Figure 7A). Furthermore, the top 1,500 variable genes (| log2(FC)| > 0.5 and FDR < 0.05) were extracted by DEG analysis. Gene expression levels of the most significant DEGs of multiple inflammation cells (astrocyte, dendritic, Div-myeloid, endothelial, ependymal, fibroblast, lymphocyte, macrophage, microglia, monocyte, neutrophil, oligodendrocyte, OPC, and pericyte) and neuron were illustrated in Figure 7B. In addition, the proportion of these cells in different SCI samples was also shown in the bar plot (Figure 7B). Specifically, average number and cell proportion of these cells were illustrated in Figure 7C. The feature plots of VOPP1, EGFR, and SFPQ were shown in Figure 7D. Importantly, VOPP1 was significantly highly expressed in lymphocytes, which was consistent with the previous analysis. In addition, EGFR and SFPQ were highly expressed in fibroblast and microglia, respectively. Finally, intersected cellular communication network and ligand–receptor plot showed the mechanisms of intercellular signal transduction, and fibroblast was the core cellular component among these intersected ligand–receptor pairs (Figures 7E,F). It showed that lymphocytes received signals from myeloid cells and various neurocytes including microglia and then played a role in the downstream pathological process of SCI.
Figure 7. Integrated analysis of single-cell RNA sequencing (scRNA-seq) and cellular communication. (A) Uniform Manifold Approximation and Projection (UMAP) scatter plots and cellular feature plots showing the cell types and marker genes of each cell type. (B) Dot plot showing gene expression levels of the most significant differentially expressed genes (DEGs) of multiple inflammation cells and neuron. In addition, the proportion of these cells in different spinal cord injury (SCI) samples was also shown in the bar plot. (C) Bar plots showing the average number and cell proportion of these cells. (D) Feature plots of vesicular overexpressed in cancer prosurvival protein 1 (VOPP1), epidermal growth factor receptor (EGFR), and Splicing factor proline and glutamine rich (SFPQ). (E) Network showing the intersected cellular communication genes. (F) Circle plot showing the intersected cellular communication genes.
The UMAP scatter plots and cellular feature plots showed various cell types and marker genes for corresponding cell types, respectively (Figure 8A). Cellular feature plots of marker genes in different cell types were shown in Figure 8B. The most significant DEGs of neuron and inflammation cells were annotated on the DEG heat map (Figure 8C).
Figure 8. Cell type annotation and cellular communication analysis. (A) Uniform Manifold Approximation and Projection (UMAP) scatter plots and cellular feature plots showing the cell types and marker genes for each cell type. (B) Cellular feature plots showing marker genes of different cell types. (C) Heat map showing the most significant differentially expressed genes (DEGs) of neuron and inflammation cells.
Chromatin Immunoprecipitation Sequencing and Assay for Transposase-Accessible Chromatin Using Sequencing Validation
In order to explore the role of enhancer-specific histone in modifications of eRNA transcription, ChIP-seq data of H3727ac were obtained and analyzed. The UCSC Genome Browser tracks illustrated enrichment of H3K27ac on multiple loci in various key eRNAs identified in this study (CPT1A, MAPK6, PCCA, PRRC2C, USP3, and VOPP1) (Figure 9).
Figure 9. Chromatin immunoprecipitation sequence (ChIP-seq) analysis of key differentially expressed enhancer RNAs (DEeRNAs). (A) In ChIP-seq data of H3K27ac, multiple binding peaks were found in CPT1A sequences. (B) In ChIP-seq data of H3K27ac, multiple binding peaks were found in MAPK6 sequences. (C) In ChIP-seq data of H3K27ac, multiple binding peaks were found in PCCA sequences. (D) In ChIP-seq data of H3K27ac, multiple binding peaks were found in PRRC2C sequences. (E) In ChIP-seq data of H3K27ac, multiple binding peaks were found in USP3 sequences. (F) In ChIP-seq data of H3K27ac, multiple binding peaks were found in vesicular overexpressed in cancer prosurvival protein 1 (VOPP1) sequences.
Regulatory relationships between VOPP1 and other key markers (SFPQ and H3K27) were further explored using ChIP-seq and ATAC-seq methods. We retrieved relative data from the Cistrome database to detect DNA fragments binding with SFPQ and H3K27ac in various SCI-related and inflammation-related cases. The results showed that in ChIP-seq data of SFPQ (Supplementary Figure 1) and H3K27 (Supplementary Figure 2), there were strong binding peaks in the position of VOPP1, which basically coincided with the open chromosome region shown in the ATAC-seq results (Supplementary Figure 3). Furthermore, we downloaded the original ChIP-seq data of SFPQ and H3K2 to analyze the results with UCSC Genome Browser. Binding peaks in SFPQ and H3K27 groups were obviously stronger than those in control groups in the chromosomal position of VOPP1, further validating the binding relationship.
Online Single-Cell RNA Sequencing Validation
To explore the cell subtype localization of the key eRNAs in the speculative regulation mechanisms, transcriptome combined with scRNA-seq data of SCI samples was analyzed. Sixteen clusters were identified using t-distributed stochastic neighbor embedding (t-SNE) (Supplementary Figure 4A). The heat map plot showed the top 5 marker genes of each cluster in SCI samples (Supplementary Figure 4B). The feature plots illustrated the distribution and expression level of SFPQ (Supplementary Figure 4C), VOPP1 (Supplementary Figure 4D), and EGFR (Supplementary Figure 4E) in SCI tissues.
Discussion
Multiple pathogenic factors contribute to the damage of the spinal cord and lead to SCI, which is a common orthopedic disease (Phillips et al., 2015). SCI is a pathological process with multi-gene mutation and multi-pathway dysfunction that causes serious neurological dysfunction (Forner et al., 2016). Primary injuries of SCI during the early stages are usually accompanied by secondary ischemia, tissue edema, and ischemia–reperfusion injury. Because of the non-renewable features of nerve cells and extremely high deformity rates, SCI has brought economic and psychological burdens (Smith et al., 2015). Research reported that primary damages induced by early stage of SCI caused massive neuronal apoptosis and keratinocyte regeneration. Moreover, formation of glial scars also inhibited nerve fiber growth, as well as exogenous and endogenous factors (Selvarajah et al., 2015). Currently, studies of SCI are mainly concentrated in mice, whereas research on gene expression changes in peripheral whole blood is relatively inadequate (Jin et al., 2015; Ropper et al., 2015).
In this study, on the basis of comprehensive bioinformatics analysis, DETFs, SCI-related eRNAs, and target genes were all identified. Furthermore, according to 27 SCI-related DEeRNAs and significant DETF links, link of SFPQ (TF) and VOPP1 (eRNA) was extracted (R = 0.990, p < 0.001, positive). Based on eight regulation relationships of key DEeRNAs and DETGs, the link of VOPP1 (eRNA) and EGFR (target gene) was extracted (R = 0.974, p < 0.001, positive). And by co-analyzing VOPP1 and 22 types of immune cells or immune-related pathways, Th cells (R = −0.987, p < 0.001, negative) were postulated to be the most significant immune cells. On the basis of comprehensive consideration of GSVA and correlation results of VOPP1 and 50 hallmark gene sets, VOPP1 and hallmark coagulation (R = 0.937, p < 0.001, positive) showed the strongest interaction. Eventually, positively regulated by SFPQ, VOPP1 strengthened the transient expression of EGFR. Th cells and hallmark coagulation were the downstream immune cells/pathways of VOPP1 in SCI.
Neurocytes with enhanced expression of SFPQ and EGFR signal the lymphocytes overexpressing VOPP1 to infiltrate peripheral blood, which may play important roles in apoptosis and promote the development of inflammation in patients with SCI based on scRNA-seq and cellular communication analysis. In addition, based on ChIP-seq and ATAC-seq analysis, multiple binding peaks in SFPQ and H3K27 groups were identified in the chromosomal position of VOPP1, further validating the binding relationship.
VOPP1, also known as glioblastoma-amplified secreted protein (Dunham et al., 2012) and EGFR-coamplified and overexpressed protein (ECOP) (Fang et al., 2020), is upregulated in multiple human tumors, such as gastric cancer and squamous cell carcinoma (Baras et al., 2011; Gao et al., 2015). Overwhelming evidence demonstrates that VOPP1 acts as a critical regulator of nuclear factor kappa B (NF-κB) signaling pathway, and it could be implicated in apoptosis resistance (Li et al., 2017). The nuclear factor of kappa light polypeptide gene enhancer in B-cells inhibitor alpha (IκBα) is correlated with NF-κB that prevents NF-κB from binding cognate promoters and also leads to steady-state cytoplasm localization. When stimulated by activation signals including tumor necrosis factor alpha (TNFα), IκBα is quickly degraded by phosphorylation-dependent ubiquitination (Karin and Ben-Neriah, 2000). Degradation of IκBα causes NF-κB to translocate to the nucleus and leads to the activation of multiple genes that are significant in various processes, including apoptosis, immune-related and inflammation response, and cellular proliferation (Park and James, 2005). Inhibition of VOPP1 expression in PBMCs alleviated the inflammatory response in advanced sepsis patients (Li et al., 2018). Taken together, increased VOPP1 expression could confer resistance to apoptosis and promote the development of inflammation. Further study of VOPP1 in SCI pathological process may bring insights into refining therapeutic regimens, as well as other disease states that are associated with elevated VOPP1 activity, such as asthma, arthritis, chronic/acute inflammation, and diverse tumors (Karin et al., 2002).
SFPQ refers to a multifunctional protein, which contributes to substantial biological processes, including DNA synthesis, gene expression regulation, DNA repair, RNA splicing, and apoptosis (Shav-Tal and Zipori, 2002). Importantly, SFPQ plays a critical role in neuronal differentiation and survival (Lowery et al., 2007). Furthermore, central nervous system (CNS) injury stimulates SFPQ expression, hindering neurogenesis and nerve regeneration, which may inhibit the repair of injured spinal cord (Elsaeidi et al., 2014). SFPQ is rich in paraspeckles, RNA-protein structures identified in the interchromatin space of nucleus, which increase in response to pro-inflammatory stimuli (Fox and Lamond, 2010). Additionally, SFPQ were significantly upregulated in chronic inflammatory diseases, including Crohn’s disease and ulcerative colitis (Hasler et al., 2011). Hence, this provides novel options for researchers in the field of SCI because SFPQ not only shows SCI-specific effects of inflammation that serves as a diagnostic biomarker but also can be a potential target for therapeutic intervention of SCI.
Based on integrated multinomial bioinformatics analysis and other studies, EGFR was suggested to be a target gene of VOPP1 (Wang et al., 1998; Eley et al., 2002). EGFR belongs to the receptor tyrosine kinase (RTK) superfamily, consisting of three other members, ErbB2/Neu/HER-2, ErbB3/HER-3, and ErbB4/HER-4 (Romano and Bucci, 2020). Being a neurotrophin receptor to initiate cellular signaling and regulate neuronal processes, EGFR has many important roles in the CNS (Hennigan et al., 2007). Moreover, EGFR is important in neural stem cells self-renewal, and loss of EGFR signaling induces these cells to differentiate preferentially into glia, which may promote glia scar formation in SCI and influence functional recovery (Robson et al., 2018). Furthermore, EGFR is crucial in regulating the differentiation of precursors to astrocytes, and increased EGFR expression levels determine their differentiation to astrocytes (Burrows et al., 1997). EGFR also plays an important role in astrocytes’ morphology. Blockade of its expression causes disorganization of astrocytes in CNS development, losing their processes surrounding neurons, which induces degeneration of abundant axons of nerve (Liu and Neufeld, 2007). Modulation of EGFR expression may be propitious to activate regeneration and counteract neurodegeneration (Ceresoli, 2012). Therefore, EGFR is critical for differentiation, growth, and repair of the injured tissue in the spinal cord, which could be a potential target for therapy of patients with SCI.
Immune response in peripheral blood plays an important role in maintaining spinal cord homeostasis. Dysfunction or disorder caused by SCI in vegetative innervation of lymphatic and endocrine systems leads to long-term abnormal inflammation responses (Zhang Y. et al., 2013). Differentiating to functionally distinct Th subsets, Th cells are crucial in normal immune surveillance and proper immune regulation (Murphy and Reiner, 2002). Importantly, SCI is associated with immune depression syndrome, owing to the dysregulated hypothalamic–pituitary–adrenal (HPA) axis and dysfunctioned sympathetic nervous system (Jones, 2014). Based on bioinformatics analysis and other studies, a rapid decrease in Th cells was identified in patients with SCI, contributing to immunosuppression in the acute phase of SCI (Gao et al., 2021). We postulated that overexpression of VOPP1 decreased the counts of Th cells in peripheral blood of SCI patients. It is an important reason that acute SCI patient becoming chronic, revealing novel targets for future SCI immunotherapy.
SCI is correlated with a remarkable thrombophilia and aberrant coagulation. In addition, the relation between deep venous thrombosis (DVT) and SCI has been well established (Yao, 1992). Numerous experimental research with fibrinogen scanning has shown the presence of DVT in nearly 100% of patients with acute SCI (Maynard et al., 1997). Sequelae from SCI affect the three components of Virchow triad, which are responsible for DVT, blood flow deceleration, vascular endothelial damage, release of procoagulants, and coagulation cascade activation, resulting in hypercoagulability state. Furthermore, patients with SCI often show dehydration and secondary increase of blood viscosity, aggravating the stasis and hypoxia in injured spinal cord (Ersoz et al., 1999). On the basis of bioinformatics analysis and other studies, coagulation was postulated to contribute to the pathogenesis and progress of SCI via the positive regulation of VOPP1, which is a promising entry point of therapeutic interventions for SCI (Ghlissi et al., 2019).
TSA refers to a highly specific hydroxamic acid that can inhibit histone deacetylase (HDAC) enzymes. Importantly, HDACs repress repressor element 1 silencing transcription factor (REST) that is able to counteract neuronal differentiation traits and Sp1, a TF that can mediate neuronal antioxidant signaling pathways. TSA leads to histone hyperacetylation, accompanied by activation of antioxidant gene expression and neuronal maturation, which indicate TF derepression of REST and Sp1 (Hakimi et al., 2002; Ryu et al., 2003). Moreover, TSA can block immune cell proliferation and suppress pro-Th1 factor IFN-γ, which could cause the transformation of Th1 to Th2 phenotype, indicating neuroimmunoprotective effects (Dangond et al., 2004). Importantly, TSA can also inhibit CNS immune cell infiltration and lipid breakdown products absorbing of microglial cells, which are consistent with anti-microglial activation of TSA recently reported (Konishi et al., 2002). Furthermore, a recent treatment showed TSA was able to decrease nitrosylation of spinal cord tissues, protecting nerve cells from free radical attack after SCI (Dasgupta et al., 2003). TSA increases histone acetylation and E2F-dependent transcripts in injured spinal cords, indicating its significant neuroprotective effect. Taken together, the general effect of TSA is to adjust dysregulated homeostatic processes, promoting histone acetylation and harnessing much more favorable genes than pathogenicity gene such as chemokine and pro-Th1, which may be effective for SCI treatment.
There were still several limitations to this study. First, the data were obtained from public sources statistically imperfect with limited samples. Second, information on other confounding variables in this study, such as smoking, was not available. Third, a prospective study is needed to evaluate the significance of these key biomarkers in terms of long-term clinical outcomes and possible applications of molecular drugs for SCI therapy. Finally, further experiment is an absolute necessity in demonstrating the regulatory mechanisms of key eRNAs implicated in SCI. Therefore, ChIP-seq data of H3K27ac from online databases were obtained and analyzed, which broadened the scope of validation and supplemented the specific regulatory mechanisms of eRNA action involved in the pathogenesis of SCI. ATAC-seq data of VOPP1 (key eRNA) were also utilized to validate the eRNA regulatory mechanisms. Moreover, the cell subtype localizations of the key eRNAs and TFs were identified by scRNA-seq validation fluorescence immunohistochemistry. Additionally, a comprehensive transcriptome bioinformatics analysis of spatial transcriptome and scRNA-seq, fluorescence immunohistochemistry, and eRNA-related direct mechanism experiments would be the further research directions.
Conclusion
Based on integrated multinomial bioinformatics analysis, we found that SFPQ was the most significant TF and VOPP1 was the most significant key eRNA in the progression of SCI patients. In addition, during this pathological process, VOPP1 upregulated the transient expression of EGFR. Th cells and hallmark coagulation were the potential downstream pathways of VOPP1. Moreover, this study provided candidate small-molecule compounds as potential targets for the treatment of SCI patients.
Data Availability Statement
Publicly available datasets were analyzed in this study. This data can be found here: Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/), ArrayExpress (https://www.ebi.ac.uk/arrayexpress/), and Cistrome Cancer database (http://cistrome.org/).
Author Contributions
RH, SW, RZ, SX, ZH, LC, and JZ: conception and design, collection and assembly of data, data analysis and interpretation, the manuscript writing, and final approval of the manuscript. All authors contributed to the article and approved the submitted version.
Funding
This study was supported in part by the National Natural Science Foundation of China (Nos. 81501203 and 81820108013), the Shanghai Municipal Health Commission (No. 201940306), the Henan Medical Science and Technology research project (No. 201602031), and the Key Project of Provincial and Ministerial Co-construction of Henan Medical Science and Technology (No. SBGJ202002031). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2021.728242/full#supplementary-material
Supplementary Figure 1 | In SFPQ Chromatin immunoprecipitation sequence (ChIP-seq) data, multiple binding peaks were found in VOPP1 sequences.
Supplementary Figure 2 | In H3K27 Chromatin immunoprecipitation sequence (ChIP-seq) data, multiple binding peaks were identified in VOPP1 sequences.
Supplementary Figure 3 | In SFPQ Assay for Transposase-Accessible Chromatin with high-throughput sequencing (ATAC-seq) data, multiple binding peaks were identified in VOPP1 sequences.
Supplementary Figure 4 | Single-cell RNA-seq validation. (A) A total of 25 numbered clusters were identified by t-NSE. (B) Feature plots illustrating the distribution and expression of SFPQ. (C) Feature plots illustrating the distribution and expression of VOPP1. (D) Feature plots illustrating the distribution and expression of EGFR. (E) Heat map showing the expression level of the top 5 marker genes of each cluster.
Supplementary Table 1 | Baseline information of 27 normal peripheral blood mononuclear cell (PBMC) samples and 25 spinal cord injury PBMC samples.
Footnotes
- ^ https://www.ncbi.nlm.nih.gov/geo/
- ^ https://www.ebi.ac.uk/arrayexpress/
- ^ https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL21975
- ^ http://cistrome.org/
- ^ http://www.immport.org/
- ^ http://software.broadinstitute.org/gsea/msigdb/
- ^ https://cibersort.stanford.edu/
- ^ https://portals.broadinstitute.org/cmap/
- ^ http://clue.io/
- ^ https://www.10xgenomics.com/instruments/chromium-x-series
- ^ http://10xgenomics.com/
- ^ http://biocc.hrbmu.edu.cn/CellMarker/
- ^ https://github.com/Coolgenome/iTALK
- ^ http://genome.ucsc.edu
- ^ https://github.com/kassambara/ggpubr
- ^ https://liulab-dfci.github.io/
- ^ http://younglab.wi.mit.edu/super_enhancer_code.html
- ^ https://www.ebi.ac.uk/gxa/sc/experiments/E-MTAB-7320/results/tsne
- ^ www.r-project.org
References
Anderson, M. A., O’shea, T. M., Burda, J. E., Ao, Y., Barlatey, S. L., Bernstein, A. M., et al. (2018). Required growth facilitators propel axon regeneration across complete spinal cord injury. Nature 561, 396–400. doi: 10.1038/s41586-018-0467-6
Athar, A., Fullgrabe, A., George, N., Iqbal, H., Huerta, L., Ali, A., et al. (2019). ArrayExpress update - from bulk to single-cell expression data. Nucleic Acids Res. 47, D711–D715. doi: 10.1093/nar/gky964
Badhiwala, J. H., Wilson, J. R., and Fehlings, M. G. (2019). Global burden of traumatic brain and spinal cord injury. Lancet Neurol. 18, 24–25. doi: 10.1016/S1474-4422(18)30444-7
Baras, A. S., Solomon, A., Davidson, R., and Moskaluk, C. A. (2011). Loss of VOPP1 overexpression in squamous carcinoma cells induces apoptosis through oxidative cellular injury. Lab. Investig. 91, 1170–1180. doi: 10.1038/labinvest.2011.70
Barrett, T., Troup, D. B., Wilhite, S. E., Ledoux, P., Evangelista, C., Kim, I. F., et al. (2011). NCBI GEO: archive for functional genomics data sets–10 years on. Nucleic Acids Res. 39, D1005–D1010. doi: 10.1093/nar/gkq1184
Beck, K. D., Nguyen, H. X., Galvan, M. D., Salazar, D. L., Woodruff, T. M., and Anderson, A. J. (2010). Quantitative analysis of cellular inflammation after traumatic spinal cord injury: evidence for a multiphasic inflammatory response in the acute to chronic environment. Brain 133, 433–447. doi: 10.1093/brain/awp322
Bhattacharya, S., Andorf, S., Gomes, L., Dunn, P., Schaefer, H., Pontius, J., et al. (2014). ImmPort: disseminating data to the public for the future of immunology. Immunol. Res. 58, 234–239. doi: 10.1007/s12026-014-8516-1
Blackwood, E. M., and Kadonaga, J. T. (1998). Going the distance: a current view of enhancer action. Science 281, 60–63. doi: 10.1126/science.281.5373.60
Borton, D., Bonizzato, M., Beauparlant, J., Digiovanna, J., Moraud, E. M., Wenger, N., et al. (2014). Corticospinal neuroprostheses to restore locomotion after spinal cord injury. Neurosci. Res. 78, 21–29. doi: 10.1016/j.neures.2013.10.001
Buenrostro, J. D., Giresi, P. G., Zaba, L. C., Chang, H. Y., and Greenleaf, W. J. (2013). Transposition of native chromatin for fast and sensitive epigenomic profiling of open chromatin, DNA-binding proteins and nucleosome position. Nat. Methods 10, 1213–1218. doi: 10.1038/nmeth.2688
Burrows, R. C., Wancio, D., Levitt, P., and Lillien, L. (1997). Response diversity and the timing of progenitor cell maturation are regulated by developmental changes in EGFR expression in the cortex. Neuron 19, 251–267. doi: 10.1016/S0896-6273(00)80937-X
Butler, A., Hoffman, P., Smibert, P., Papalexi, E., and Satija, R. (2018). Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat. Biotechnol. 36, 411–420. doi: 10.1038/nbt.4096
Ceresoli, G. L. (2012). Role of EGFR inhibitors in the treatment of central nervous system metastases from non-small cell lung cancer. Curr. Cancer Drug. Targets 12, 237–246. doi: 10.2174/156800912799277430
Chen, K. B., Uchida, K., Nakajima, H., Yayama, T., Hirai, T., Watanabe, S., et al. (2011). Tumor necrosis factor-alpha antagonist reduces apoptosis of neurons and oligodendroglia in rat spinal cord injury. Spine 36, 1350–1358. doi: 10.1097/BRS.0b013e3181f014ec
Chen, W. T., Lu, A., Craessaerts, K., Pavie, B., Sala Frigerio, C., Corthout, N., et al. (2020). Spatial Transcriptomics and In Situ Sequencing to Study Alzheimer’s Disease. Cell 182, 976–991.e919. doi: 10.1016/j.cell.2020.06.038
Dangond, F., Hwang, D., Camelo, S., Pasinelli, P., Frosch, M. P., Stephanopoulos, G., et al. (2004). Molecular signature of late-stage human ALS revealed by expression profiling of postmortem spinal cord gray matter. Physiol. Genomics 16, 229–239. doi: 10.1152/physiolgenomics.00087.2001
Dasgupta, S., Zhou, Y., Jana, M., Banik, N. L., and Pahan, K. (2003). Sodium phenylacetate inhibits adoptive transfer of experimental allergic encephalomyelitis in SJL/J mice at multiple steps. J. Immunol. 170, 3874–3882. doi: 10.4049/jimmunol.170.7.3874
Delile, J., Rayon, T., Melchionda, M., Edwards, A., Briscoe, J., and Sagner, A. (2019). Single cell transcriptomics reveals spatial and temporal dynamics of gene expression in the developing mouse spinal cord. Development 146:dev173807. doi: 10.1242/dev.173807
Dunham, I., Kundaje, A., Aldred, S. F., Collins, P. J., Davis, C., Doyle, F., et al. (2012). An integrated encyclopedia of DNA elements in the human genome. Nature 489, 57–74. doi: 10.1038/nature11247
Durinck, S., Spellman, P. T., Birney, E., and Huber, W. (2009). Mapping identifiers for the integration of genomic datasets with the R/Bioconductor package biomaRt. Nat. Protoc. 4, 1184–1191. doi: 10.1038/nprot.2009.97
Eley, G. D., Reiter, J. L., Pandita, A., Park, S., Jenkins, R. B., Maihle, N. J., et al. (2002). A chromosomal region 7p11.2 transcript map: its development and application to the study of EGFR amplicons in glioblastoma. Neuron. Oncol. 4, 86–94. doi: 10.1093/neuonc/4.2.86
Elsaeidi, F., Bemben, M. A., Zhao, X. F., and Goldman, D. (2014). Jak/Stat signaling stimulates zebrafish optic nerve regeneration and overcomes the inhibitory actions of Socs3 and Sfpq. J. Neurosci. 34, 2632–2644. doi: 10.1523/JNEUROSCI.3898-13.2014
ENCODE Project Consortium (2012). An integrated encyclopedia of DNA elements in the human genome. Nature 489, 57–74. doi: 10.1038/nature11247
Ersoz, G., Ficicilar, H., Pasin, M., Yorgancioglu, R., and Yavuzer, S. (1999). Platelet aggregation in traumatic spinal cord injury. Spinal Cord 37, 644–647. doi: 10.1038/sj.sc.3100903
Fang, Z., Wu, L., Dai, H., Hu, P., Wang, B., Han, Q., et al. (2020). The role of vesicular overexpressed in cancer pro-survival protein 1 in hepatocellular carcinoma proliferation. Cancer Biomark 28, 9–20. doi: 10.3233/CBM-190574
Feng, J., Liu, T., Qin, B., Zhang, Y., and Liu, X. S. (2012). Identifying ChIP-seq enrichment using MACS. Nat. Protoc. 7, 1728–1740. doi: 10.1038/nprot.2012.101
Ferreira, M. R., Santos, G. A., Biagi, C. A., Silva Junior, W. A., and Zambuzzi, W. F. (2021). GSVA score reveals molecular signatures from transcriptomes for biomaterials comparison. J. Biomed. Mater. Res. A. 109, 1004–1014. doi: 10.1002/jbm.a.37090
Forner, S., Martini, A. C., Andrade, E. L., and Rae, G. A. (2016). Neuropathic pain induced by spinal cord injury: Role of endothelin ETA and ETB receptors. Neurosci. Lett. 617, 14–21. doi: 10.1016/j.neulet.2016.02.005
Fox, A. H., and Lamond, A. I. (2010). Paraspeckles. Cold Spring Harb. Perspect. Biol. 2:a000687. doi: 10.1101/cshperspect.a000687
Gao, C., Pang, M., Zhou, Z., Long, S., Dong, D., Yang, J., et al. (2015). Epidermal growth factor receptor-coamplified and overexpressed protein (VOPP1) is a putative oncogene in gastric cancer. Clin. Exp. Med. 15, 469–475. doi: 10.1007/s10238-014-0320-7
Gao, T. Y., Huang, F. F., Xie, Y. Y., Wang, W. Q., Wang, L. D., Mu, D., et al. (2021). Dynamic changes in the systemic immune responses of spinal cord injury model mice. Neural. Regen. Res. 16, 382–387. doi: 10.4103/1673-5374.290910
Gautier, L., Cope, L., Bolstad, B. M., and Irizarry, R. A. (2004). affy–analysis of Affymetrix GeneChip data at the probe level. Bioinformatics 20, 307–315. doi: 10.1093/bioinformatics/btg405
GBD 2016 Neurology Collaborators (2019). Global, regional, and national burden of neurological disorders, 1990-2016: a systematic analysis for the Global Burden of Disease Study 2016. Lancet Neurol. 18, 459–480.
Ghlissi, Z., Krichen, F., Kallel, R., Amor, I. B., Boudawara, T., Gargouri, J., et al. (2019). Sulfated polysaccharide isolated from Globularia alypum L.: Structural characterization, in vivo and in vitro anticoagulant activity, and toxicological profile. Int. J. Biol. Macromol. 123, 335–342. doi: 10.1016/j.ijbiomac.2018.11.044
Grubert, F., Srivas, R., Spacek, D. V., Kasowski, M., Ruiz-Velasco, M., Sinnott-Armstrong, N., et al. (2020). Landscape of cohesin-mediated chromatin loops in the human genome. Nature 583, 737–743. doi: 10.1038/s41586-020-2151-x
Guenther, M. G., Frampton, G. M., Soldner, F., Hockemeyer, D., Mitalipova, M., Jaenisch, R., et al. (2010). Chromatin structure and gene expression programs of human embryonic and induced pluripotent stem cells. Cell Stem Cell 7, 249–257. doi: 10.1016/j.stem.2010.06.015
Hakimi, M. A., Bochar, D. A., Chenoweth, J., Lane, W. S., Mandel, G., and Shiekhattar, R. (2002). A core-BRAF35 complex containing histone deacetylase mediates repression of neuronal-specific genes. Proc. Natl. Acad. Sci. U.S.A. 99, 7420–7425. doi: 10.1073/pnas.112008599
Hanzelmann, S., Castelo, R., and Guinney, J. (2013). GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinform. 14:7. doi: 10.1186/1471-2105-14-7
Hasler, R., Kerick, M., Mah, N., Hultschig, C., Richter, G., Bretz, F., et al. (2011). Alterations of pre-mRNA splicing in human inflammatory bowel disease. Eur. J. Cell Biol. 90, 603–611. doi: 10.1016/j.ejcb.2010.11.010
Hayashi, M., Ueyama, T., Nemoto, K., Tamaki, T., and Senba, E. (2000). Sequential mRNA expression for immediate early genes, cytokines, and neurotrophins in spinal cord injury. J. Neurotrauma. 17, 203–218. doi: 10.1089/neu.2000.17.203
Hennigan, A., O’callaghan, R. M., and Kelly, A. M. (2007). Neurotrophins and their receptors: roles in plasticity, neurodegeneration and neuroprotection. Biochem. Soc. Trans. 35, 424–427. doi: 10.1042/BST0350424
Irizarry, R. A., Hobbs, B., Collin, F., Beazer-Barclay, Y. D., Antonellis, K. J., Scherf, U., et al. (2003). Exploration, normalization, and summaries of high density oligonucleotide array probe level data. Biostatistics 4, 249–264. doi: 10.1093/biostatistics/4.2.249
Jin, Y., Bouyer, J., Haas, C., and Fischer, I. (2015). Evaluation of the anatomical and functional consequences of repetitive mild cervical contusion using a model of spinal concussion. Exp. Neurol. 271, 175–188. doi: 10.1016/j.expneurol.2015.06.001
Jones, T. B. (2014). Lymphocytes and autoimmunity after spinal cord injury. Exp. Neurol. 258, 78–90. doi: 10.1016/j.expneurol.2014.03.003
Kang, Y., Kim, Y. W., Kang, J., and Kim, A. (2021). Histone H3K4me1 and H3K27ac play roles in nucleosome eviction and eRNA transcription, respectively, at enhancers. FASEB J. 35:e21781. doi: 10.1096/fj.202100488R
Karin, M., and Ben-Neriah, Y. (2000). Phosphorylation meets ubiquitination: the control of NF-[kappa]B activity. Annu. Rev. Immunol. 18, 621–663. doi: 10.1146/annurev.immunol.18.1.621
Karin, M., Cao, Y., Greten, F. R., and Li, Z. W. (2002). NF-kappaB in cancer: from innocent bystander to major culprit. Nat. Rev. Cancer 2, 301–310. doi: 10.1038/nrc780
Kim, T. K., Hemberg, M., Gray, J. M., Costa, A. M., Bear, D. M., Wu, J., et al. (2010). Widespread transcription at neuronal activity-regulated enhancers. Nature 465, 182–U165. doi: 10.1038/nature09033
Kohl, M., Wiese, S., and Warscheid, B. (2011). Cytoscape: software for visualization and analysis of biological networks. Methods Mol. Biol. 696, 291–303. doi: 10.1007/978-1-60761-987-1_18
Konishi, Y., Lehtinen, M., Donovan, N., and Bonni, A. (2002). Cdc2 phosphorylation of BAD links the cell cycle to the cell death machinery. Mol. Cell 9, 1005–1016. doi: 10.1016/S1097-2765(02)00524-5
Lamb, J., Crawford, E. D., Peck, D., Modell, J. W., Blat, I. C., Wrobel, M. J., et al. (2006). The Connectivity Map: using gene-expression signatures to connect small molecules, genes, and disease. Science 313, 1929–1935. doi: 10.1126/science.1132939
Li, J. M., Zhang, H., and Zuo, Y. J. (2018). MicroRNA-218 alleviates sepsis inflammation by negatively regulating VOPP1 via JAK/STAT pathway. Eur. Rev. Med. Pharmacol. Sci. 22, 5620–5626.
Li, W., Hu, Y., Oh, S., Ma, Q., Merkurjev, D., Song, X., et al. (2015). Condensin I and II Complexes License Full Estrogen Receptor alpha-Dependent Enhancer Activation. Mol. Cell 59, 188–202. doi: 10.1016/j.molcel.2015.06.002
Li, W., Notani, D., and Rosenfeld, M. G. (2016). Enhancers as non-coding RNA transcription units: recent insights and future perspectives. Nat. Rev. Genet. 17, 207–223. doi: 10.1038/nrg.2016.4
Li, Y. J., Zhang, W., Xia, H., Zhang, B. S., Chen, P., Zhao, Y. L., et al. (2017). miR-218 suppresses epithelial-to-mesenchymal transition by targeting Robo1 and Ecop in lung adenocarcinoma cells. Fut. Oncol. 13, 2571–2582. doi: 10.2217/fon-2017-0398
Liberzon, A., Birger, C., Thorvaldsdottir, H., Ghandi, M., Mesirov, J. P., and Tamayo, P. (2015). The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst. 1, 417–425. doi: 10.1016/j.cels.2015.12.004
Lister, R., Pelizzola, M., Dowen, R. H., Hawkins, R. D., Hon, G., Tonti-Filippini, J., et al. (2009). Human DNA methylomes at base resolution show widespread epigenomic differences. Nature 462, 315–322. doi: 10.1038/nature08514
Liu, B., and Neufeld, A. H. (2007). Activation of epidermal growth factor receptors in astrocytes: from development to neural injury. J. Neurosci. Res. 85, 3523–3529. doi: 10.1002/jnr.21364
Lowery, L. A., Rubin, J., and Sive, H. (2007). Whitesnake/sfpq is required for cell survival and neuronal development in the zebrafish. Dev. DYN 236, 1347–1357. doi: 10.1002/dvdy.21132
Maynard, F. M., Bracken, M. B., Creasey, G., Ditunno, J. F., Donovan, W. H., Ducker, T. B., et al. (1997). International standards for neurological and functional classification of spinal cord injury. Spinal. Cord 35, 266–274. doi: 10.1038/sj.sc.3100432
Murphy, K. M., and Reiner, S. L. (2002). The lineage decisions of helper T cells. Nat. Rev. Immunol. 2, 933–944. doi: 10.1038/nri954
Nakae, A., Nakai, K., Yano, K., Hosokawa, K., Shibata, M., and Mashimo, T. (2011). The animal model of spinal cord injury as an experimental pain model. J. Biomed. Biotechnol. 2011:939023. doi: 10.1155/2011/939023
Newman, A. M., Liu, C. L., Green, M. R., Gentles, A. J., Feng, W., Xu, Y., et al. (2015). Robust enumeration of cell subsets from tissue expression profiles. Nat. Methods 12, 453–457. doi: 10.1038/nmeth.3337
Park, S., and James, C. D. (2005). ECop (EGFR-coamplified and overexpressed protein), a novel protein, regulates NF-kappaB transcriptional activity and associated apoptotic response in an IkappaBalpha-dependent manner. Oncogene 24, 2495–2502. doi: 10.1038/sj.onc.1208496
Phillips, A., Matin, N., West, C., Zheng, A., Galea, L., Dorrance, A., et al. (2015). Autonomic Dysreflexia Impairs Cerebrovascular Health and Cognition in Experimental Spinal Cord Injury. FASEB J. 29, 59–70. doi: 10.1096/fasebj.29.1_supplement.800.10
Pineau, I., and Lacroix, S. (2007). Proinflammatory cytokine synthesis in the injured mouse spinal cord: multiphasic expression pattern and identification of the cell types involved. J. Comp. Neurol. 500, 267–285. doi: 10.1002/cne.21149
Pouw, M. H., Hosman, A. J., Van Kampen, A., Hirschfeld, S., Thietje, R., and Van De Meent, H. (2011). Is the outcome in acute spinal cord ischaemia different from that in traumatic spinal cord injury? A cross-sectional analysis of the neurological and functional outcome in a cohort of 93 paraplegics. Spinal. Cord 49, 307–312. doi: 10.1038/sc.2010.114
Robinson, J. T., Thorvaldsdottir, H., Winckler, W., Guttman, M., Lander, E. S., Getz, G., et al. (2011). Integrative genomics viewer. Nat. Biotechnol. 29, 24–26. doi: 10.1038/nbt.1754
Robson, J. P., Wagner, B., Glitzner, E., Heppner, F. L., Steinkellner, T., Khan, D., et al. (2018). Impaired neural stem cell expansion and hypersensitivity to epileptic seizures in mice lacking the EGFR in the brain. FEBS J. 285, 3175–3196. doi: 10.1111/febs.14603
Romano, R., and Bucci, C. (2020). Role of EGFR in the Nervous System. Cells 9:1887. doi: 10.3390/cells9081887
Ropper, A. E., Zeng, X., Anderson, J. E., Yu, D., Han, I., Haragopal, H., et al. (2015). An efficient device to experimentally model compression injury of mammalian spinal cord. Exp. Neurol. 271, 515–523. doi: 10.1016/j.expneurol.2015.07.012
Rudman, M. D., Choi, J. S., Lee, H. E., Tan, S. K., Ayad, N. G., and Lee, J. K. (2018). Bromodomain and extraterminal domain-containing protein inhibition attenuates acute inflammation after spinal cord injury. Exp. Neurol. 309, 181–192. doi: 10.1016/j.expneurol.2018.08.005
Ryu, H., Lee, J., Olofsson, B. A., Mwidau, A., Dedeoglu, A., Escudero, M., et al. (2003). Histone deacetylase inhibitors prevent oxidative neuronal death independent of expanded polyglutamine repeats via an Sp1-dependent pathway. Proc. Natl. Acad. Sci. U.S.A. 100, 4281–4286. doi: 10.1073/pnas.0737363100
Selvarajah, S., Hammond, E. R., and Schneider, E. B. (2015). Trends in Traumatic Spinal Cord Injury. JAMA 314:1643. doi: 10.1001/jama.2015.11194
Shav-Tal, Y., and Zipori, D. (2002). PSF and p54(nrb)/NonO–multi-functional nuclear proteins. FEBS Lett. 531, 109–114. doi: 10.1016/S0014-5793(02)03447-6
Smith, E., O’reilly, A., Synnott, K., Morris, S., and Timlin, M. (2015). Review of Time to Surgical Decompression in Traumatic Spinal Cord Injured Patients. Ir. Med. J. 108, 265–267.
Smyth, G. K. (2004). Linear models and empirical bayes methods for assessing differential expression in microarray experiments. Stat. Appl. Genet. Mol. Biol. 3:Article3. doi: 10.2202/1544-6115.1027
Subramanian, A., Kuehn, H., Gould, J., Tamayo, P., and Mesirov, J. P. (2007). GSEA-P: a desktop application for Gene Set Enrichment Analysis. Bioinformatics 23, 3251–3253. doi: 10.1093/bioinformatics/btm369
Wang, H., Zou, J., Zhao, B., Johannsen, E., Ashworth, T., Wong, H., et al. (2011). Genome-wide analysis reveals conserved and divergent features of Notch1/RBPJ binding in human and murine T-lymphoblastic leukemia cells. Proc. Natl. Acad. Sci. U.S.A. 108, 14908–14913. doi: 10.1073/pnas.1109023108
Wang, X. Y., Smith, D. I., Frederick, L., and James, C. D. (1998). Analysis of EGF receptor amplicons reveals amplification of multiple expressed sequences. Oncogene 16, 191–195. doi: 10.1038/sj.onc.1201476
Wang, Y., Wang, R., Zhang, S., Song, S., Jiang, C., Han, G., et al. (2019). iTALK: an R Package to Characterize and Illustrate Intercellular Communication. bioRxiv [preprint]. doi: 10.1101/507871
West, J. A., Davis, C. P., Sunwoo, H., Simon, M. D., Sadreyev, R. I., Wang, P. I., et al. (2014). The long noncoding RNAs NEAT1 and MALAT1 bind active chromatin sites. Mol. Cell. 55, 791–802. doi: 10.1016/j.molcel.2014.07.012
Whyte, W. A., Orlando, D. A., Hnisz, D., Abraham, B. J., Lin, C. Y., Kagey, M. H., et al. (2013). Master transcription factors and mediator establish super-enhancers at key cell identity genes. Cell 153, 307–319. doi: 10.1016/j.cell.2013.03.035
Wolf, F. A., Angerer, P., and Theis, F. J. (2018). SCANPY: large-scale single-cell gene expression data analysis. Genome Biol. 19:15. doi: 10.1186/s13059-017-1382-0
Xiao, B., Liu, L., Li, A., Xiang, C., Wang, P., Li, H., et al. (2020). Identification and Verification of Immune-Related Gene Prognostic Signature Based on ssGSEA for Osteosarcoma. Front. Oncol. 10:607622. doi: 10.3389/fonc.2020.607622
Yao, J. S. (1992). Deep vein thrombosis in spinal cord-injured patients. Evaluation and assessment. Chest 102, 645S–648S. doi: 10.1378/chest.102.6_Supplement.645S
Yu, G., Wang, L. G., Han, Y., and He, Q. Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 16, 284–287. doi: 10.1089/omi.2011.0118
Zhang, T., Jiang, M., Chen, L., Niu, B., and Cai, Y. (2013). Prediction of gene phenotypes based on GO and KEGG pathway enrichment scores. Biomed. Res. Int. 2013:870795. doi: 10.1155/2013/870795
Zhang, X., Lan, Y., Xu, J., Quan, F., Zhao, E., Deng, C., et al. (2019). CellMarker: a manually curated resource of cell markers in human and mouse. Nucleic Acids Res. 47, D721–D728. doi: 10.1093/nar/gky900
Zhang, Y., Guan, Z., Reader, B., Shawler, T., Mandrekar-Colucci, S., Huang, K., et al. (2013). Autonomic dysreflexia causes chronic immune suppression after spinal cord injury. J. Neurosci. 33, 12970–12981. doi: 10.1523/JNEUROSCI.1974-13.2013
Zheng, R., Wan, C., Mei, S., Qin, Q., Wu, Q., Sun, H., et al. (2019). Cistrome Data Browser: expanded datasets and new tools for gene regulatory analysis. Nucleic Acids Res. 47, D729–D735. doi: 10.1093/nar/gky1094
Keywords: spinal cord injury, multi-omics bioinformatics analysis, biomarkers, eRNA regulatory network, therapeutic targets
Citation: Huang R, Wang S, Zhu R, Xian S, Huang Z, Cheng L and Zhang J (2021) Identification of Key eRNAs for Spinal Cord Injury by Integrated Multinomial Bioinformatics Analysis. Front. Cell Dev. Biol. 9:728242. doi: 10.3389/fcell.2021.728242
Received: 21 June 2021; Accepted: 25 August 2021;
Published: 11 October 2021.
Edited by:
Chuan-Xing Li, Karolinska Institutet (KI), SwedenReviewed by:
Zhi-an Cheng, Guangdong Provincial Hospital of Chinese Medicine, ChinaLixin Cheng, Jinan University, China
Copyright © 2021 Huang, Wang, Zhu, Xian, Huang, Cheng and Zhang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Zongqiang Huang, Z3podWFuZ3pxQDE2My5jb20=; Liming Cheng, bGltaW5nY2hlbmdAdG9uZ2ppLmVkdS5jbg==; Jie Zhang, amllemhhbmdAdG9uZ2ppLmVkdS5jbg==
†These authors have contributed equally to this work