Abstract
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 (; ). 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 (; ). 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) (). 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 (). 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 (). 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) (). Numerous eRNAs have been identified in Homo sapiens so far, and many of them were suggested to play crucial parts in the transcriptional circuitry ().
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 (). Large amounts of cytokines are quickly generated and released from various cells in the damaged structures after SCI (; ). Cytokines activate apoptosis of pathological nerve cells and recruitment of active leukocytes in myelopathic lesions (; ). 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 (). 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 () and E-GEOD-69901 (ArrayExpress)2 (), 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 (). Hallmark signaling pathways were collected from the Molecular Signatures Database (MSigDB, v7.2)6 (). 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) ().
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 (; ). 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 (). 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 (). 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 () and GSEA () 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; ). 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) ().
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 (). 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 (), 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 (). 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 (; ).
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 (). 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 (). 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 (). Role of H3K27ac in eRNA transcription was evaluated by analyzing ChIP-seq data (accession number: GSE134744) in peripheral blood ().
The ChIP-seq data of Splicing factor proline and glutamine rich (SFPQ) (accession number: GSM2827312, GSM1411215, GSM2825596, GSM1097497, GSM2827311, and GSM1097496) (; ; ) and H3K27 (accession number: GSM732912, GSM575294, and GSM663427) (; ; ) 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 (; ).
Assay for Transposase-Accessible Chromatin Using Sequencing Validation
The ATAC-seq () 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) () and performed the correlation analysis using ggpubr.15 In order to explore the potential super-enhancers, we performed the peak calling with macs14 ().16 Furthermore, we ranked enhancers and determined super-enhancers using ROSE ().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 ().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
FIGURE 2
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
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
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
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
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
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
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
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 (). SCI is a pathological process with multi-gene mutation and multi-pathway dysfunction that causes serious neurological dysfunction (). 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 (). 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 (). Currently, studies of SCI are mainly concentrated in mice, whereas research on gene expression changes in peripheral whole blood is relatively inadequate (; ).
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 () and EGFR-coamplified and overexpressed protein (ECOP) (), is upregulated in multiple human tumors, such as gastric cancer and squamous cell carcinoma (; ). 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 (). 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 (). 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 (). Inhibition of VOPP1 expression in PBMCs alleviated the inflammatory response in advanced sepsis patients (). 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 ().
SFPQ refers to a multifunctional protein, which contributes to substantial biological processes, including DNA synthesis, gene expression regulation, DNA repair, RNA splicing, and apoptosis (). Importantly, SFPQ plays a critical role in neuronal differentiation and survival (). Furthermore, central nervous system (CNS) injury stimulates SFPQ expression, hindering neurogenesis and nerve regeneration, which may inhibit the repair of injured spinal cord (). SFPQ is rich in paraspeckles, RNA-protein structures identified in the interchromatin space of nucleus, which increase in response to pro-inflammatory stimuli (). Additionally, SFPQ were significantly upregulated in chronic inflammatory diseases, including Crohn’s disease and ulcerative colitis (). 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 (; ). 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 (). Being a neurotrophin receptor to initiate cellular signaling and regulate neuronal processes, EGFR has many important roles in the CNS (). 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 (). Furthermore, EGFR is crucial in regulating the differentiation of precursors to astrocytes, and increased EGFR expression levels determine their differentiation to astrocytes (). 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 (). Modulation of EGFR expression may be propitious to activate regeneration and counteract neurodegeneration (). 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 (). Importantly, SCI is associated with immune depression syndrome, owing to the dysregulated hypothalamic–pituitary–adrenal (HPA) axis and dysfunctioned sympathetic nervous system (). 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 (). 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 (). 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 (). 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 ().
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 (; ). 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 (). 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 (). Furthermore, a recent treatment showed TSA was able to decrease nitrosylation of spinal cord tissues, protecting nerve cells from free radical attack after SCI (). 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.
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.
Statements
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.
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 1In SFPQ Chromatin immunoprecipitation sequence (ChIP-seq) data, multiple binding peaks were found in VOPP1 sequences.
Supplementary Figure 2In H3K27 Chromatin immunoprecipitation sequence (ChIP-seq) data, multiple binding peaks were identified in VOPP1 sequences.
Supplementary Figure 3In SFPQ Assay for Transposase-Accessible Chromatin with high-throughput sequencing (ATAC-seq) data, multiple binding peaks were identified in VOPP1 sequences.
Supplementary Figure 4Single-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 1Baseline information of 27 normal peripheral blood mononuclear cell (PBMC) samples and 25 spinal cord injury PBMC samples.
Footnotes
1.^https://www.ncbi.nlm.nih.gov/geo/
2.^https://www.ebi.ac.uk/arrayexpress/
3.^https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL21975
6.^http://software.broadinstitute.org/gsea/msigdb/
7.^https://cibersort.stanford.edu/
8.^https://portals.broadinstitute.org/cmap/
10.^https://www.10xgenomics.com/instruments/chromium-x-series
12.^http://biocc.hrbmu.edu.cn/CellMarker/
13.^https://github.com/Coolgenome/iTALK
15.^https://github.com/kassambara/ggpubr
16.^https://liulab-dfci.github.io/
17.^http://younglab.wi.mit.edu/super_enhancer_code.html
18.^https://www.ebi.ac.uk/gxa/sc/experiments/E-MTAB-7320/results/tsne
References
1
AndersonM. A.O’sheaT. M.BurdaJ. E.AoY.BarlateyS. L.BernsteinA. M.et al (2018). Required growth facilitators propel axon regeneration across complete spinal cord injury.Nature561396–400. 10.1038/s41586-018-0467-6
2
AtharA.FullgrabeA.GeorgeN.IqbalH.HuertaL.AliA.et al (2019). ArrayExpress update - from bulk to single-cell expression data.Nucleic Acids Res.47D711–D715. 10.1093/nar/gky964
3
BadhiwalaJ. H.WilsonJ. R.FehlingsM. G. (2019). Global burden of traumatic brain and spinal cord injury.Lancet Neurol.1824–25. 10.1016/S1474-4422(18)30444-7
4
BarasA. S.SolomonA.DavidsonR.MoskalukC. A. (2011). Loss of VOPP1 overexpression in squamous carcinoma cells induces apoptosis through oxidative cellular injury.Lab. Investig.911170–1180. 10.1038/labinvest.2011.70
5
BarrettT.TroupD. B.WilhiteS. E.LedouxP.EvangelistaC.KimI. F.et al (2011). NCBI GEO: archive for functional genomics data sets–10 years on.Nucleic Acids Res.39D1005–D1010. 10.1093/nar/gkq1184
6
BeckK. D.NguyenH. X.GalvanM. D.SalazarD. L.WoodruffT. M.AndersonA. J. (2010). Quantitative analysis of cellular inflammation after traumatic spinal cord injury: evidence for a multiphasic inflammatory response in the acute to chronic environment.Brain133433–447. 10.1093/brain/awp322
7
BhattacharyaS.AndorfS.GomesL.DunnP.SchaeferH.PontiusJ.et al (2014). ImmPort: disseminating data to the public for the future of immunology.Immunol. Res.58234–239. 10.1007/s12026-014-8516-1
8
BlackwoodE. M.KadonagaJ. T. (1998). Going the distance: a current view of enhancer action.Science28160–63. 10.1126/science.281.5373.60
9
BortonD.BonizzatoM.BeauparlantJ.DigiovannaJ.MoraudE. M.WengerN.et al (2014). Corticospinal neuroprostheses to restore locomotion after spinal cord injury.Neurosci. Res.7821–29. 10.1016/j.neures.2013.10.001
10
BuenrostroJ. D.GiresiP. G.ZabaL. C.ChangH. Y.GreenleafW. J. (2013). Transposition of native chromatin for fast and sensitive epigenomic profiling of open chromatin, DNA-binding proteins and nucleosome position.Nat. Methods101213–1218. 10.1038/nmeth.2688
11
BurrowsR. C.WancioD.LevittP.LillienL. (1997). Response diversity and the timing of progenitor cell maturation are regulated by developmental changes in EGFR expression in the cortex.Neuron19251–267. 10.1016/S0896-6273(00)80937-X
12
ButlerA.HoffmanP.SmibertP.PapalexiE.SatijaR. (2018). Integrating single-cell transcriptomic data across different conditions, technologies, and species.Nat. Biotechnol.36411–420. 10.1038/nbt.4096
13
CeresoliG. L. (2012). Role of EGFR inhibitors in the treatment of central nervous system metastases from non-small cell lung cancer.Curr. Cancer Drug. Targets12237–246. 10.2174/156800912799277430
14
ChenK. B.UchidaK.NakajimaH.YayamaT.HiraiT.WatanabeS.et al (2011). Tumor necrosis factor-alpha antagonist reduces apoptosis of neurons and oligodendroglia in rat spinal cord injury.Spine361350–1358. 10.1097/BRS.0b013e3181f014ec
15
ChenW. T.LuA.CraessaertsK.PavieB.Sala FrigerioC.CorthoutN.et al (2020). Spatial Transcriptomics and In Situ Sequencing to Study Alzheimer’s Disease.Cell182976–991.e919. 10.1016/j.cell.2020.06.038
16
DangondF.HwangD.CameloS.PasinelliP.FroschM. P.StephanopoulosG.et al (2004). Molecular signature of late-stage human ALS revealed by expression profiling of postmortem spinal cord gray matter.Physiol. Genomics16229–239. 10.1152/physiolgenomics.00087.2001
17
DasguptaS.ZhouY.JanaM.BanikN. L.PahanK. (2003). Sodium phenylacetate inhibits adoptive transfer of experimental allergic encephalomyelitis in SJL/J mice at multiple steps.J. Immunol.1703874–3882. 10.4049/jimmunol.170.7.3874
18
DelileJ.RayonT.MelchiondaM.EdwardsA.BriscoeJ.SagnerA. (2019). Single cell transcriptomics reveals spatial and temporal dynamics of gene expression in the developing mouse spinal cord.Development146:dev173807. 10.1242/dev.173807
19
DunhamI.KundajeA.AldredS. F.CollinsP. J.DavisC.DoyleF.et al (2012). An integrated encyclopedia of DNA elements in the human genome.Nature48957–74. 10.1038/nature11247
20
DurinckS.SpellmanP. T.BirneyE.HuberW. (2009). Mapping identifiers for the integration of genomic datasets with the R/Bioconductor package biomaRt.Nat. Protoc.41184–1191. 10.1038/nprot.2009.97
21
EleyG. D.ReiterJ. L.PanditaA.ParkS.JenkinsR. B.MaihleN. 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.486–94. 10.1093/neuonc/4.2.86
22
ElsaeidiF.BembenM. A.ZhaoX. F.GoldmanD. (2014). Jak/Stat signaling stimulates zebrafish optic nerve regeneration and overcomes the inhibitory actions of Socs3 and Sfpq.J. Neurosci.342632–2644. 10.1523/JNEUROSCI.3898-13.2014
23
ENCODE Project Consortium (2012). An integrated encyclopedia of DNA elements in the human genome.Nature48957–74. 10.1038/nature11247
24
ErsozG.FicicilarH.PasinM.YorganciogluR.YavuzerS. (1999). Platelet aggregation in traumatic spinal cord injury.Spinal Cord37644–647. 10.1038/sj.sc.3100903
25
FangZ.WuL.DaiH.HuP.WangB.HanQ.et al (2020). The role of vesicular overexpressed in cancer pro-survival protein 1 in hepatocellular carcinoma proliferation.Cancer Biomark289–20. 10.3233/CBM-190574
26
FengJ.LiuT.QinB.ZhangY.LiuX. S. (2012). Identifying ChIP-seq enrichment using MACS.Nat. Protoc.71728–1740. 10.1038/nprot.2012.101
27
FerreiraM. R.SantosG. A.BiagiC. A.Silva JuniorW. A.ZambuzziW. F. (2021). GSVA score reveals molecular signatures from transcriptomes for biomaterials comparison.J. Biomed. Mater. Res. A.1091004–1014. 10.1002/jbm.a.37090
28
FornerS.MartiniA. C.AndradeE. L.RaeG. A. (2016). Neuropathic pain induced by spinal cord injury: Role of endothelin ETA and ETB receptors.Neurosci. Lett.61714–21. 10.1016/j.neulet.2016.02.005
29
FoxA. H.LamondA. I. (2010). Paraspeckles.Cold Spring Harb. Perspect. Biol.2:a000687. 10.1101/cshperspect.a000687
30
GaoC.PangM.ZhouZ.LongS.DongD.YangJ.et al (2015). Epidermal growth factor receptor-coamplified and overexpressed protein (VOPP1) is a putative oncogene in gastric cancer.Clin. Exp. Med.15469–475. 10.1007/s10238-014-0320-7
31
GaoT. Y.HuangF. F.XieY. Y.WangW. Q.WangL. D.MuD.et al (2021). Dynamic changes in the systemic immune responses of spinal cord injury model mice.Neural. Regen. Res.16382–387. 10.4103/1673-5374.290910
32
GautierL.CopeL.BolstadB. M.IrizarryR. A. (2004). affy–analysis of Affymetrix GeneChip data at the probe level.Bioinformatics20307–315. 10.1093/bioinformatics/btg405
33
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.18459–480.
34
GhlissiZ.KrichenF.KallelR.AmorI. B.BoudawaraT.GargouriJ.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.123335–342. 10.1016/j.ijbiomac.2018.11.044
35
GrubertF.SrivasR.SpacekD. V.KasowskiM.Ruiz-VelascoM.Sinnott-ArmstrongN.et al (2020). Landscape of cohesin-mediated chromatin loops in the human genome.Nature583737–743. 10.1038/s41586-020-2151-x
36
GuentherM. G.FramptonG. M.SoldnerF.HockemeyerD.MitalipovaM.JaenischR.et al (2010). Chromatin structure and gene expression programs of human embryonic and induced pluripotent stem cells.Cell Stem Cell7249–257. 10.1016/j.stem.2010.06.015
37
HakimiM. A.BocharD. A.ChenowethJ.LaneW. S.MandelG.ShiekhattarR. (2002). A core-BRAF35 complex containing histone deacetylase mediates repression of neuronal-specific genes.Proc. Natl. Acad. Sci. U.S.A.997420–7425. 10.1073/pnas.112008599
38
HanzelmannS.CasteloR.GuinneyJ. (2013). GSVA: gene set variation analysis for microarray and RNA-seq data.BMC Bioinform.14:7. 10.1186/1471-2105-14-7
39
HaslerR.KerickM.MahN.HultschigC.RichterG.BretzF.et al (2011). Alterations of pre-mRNA splicing in human inflammatory bowel disease.Eur. J. Cell Biol.90603–611. 10.1016/j.ejcb.2010.11.010
40
HayashiM.UeyamaT.NemotoK.TamakiT.SenbaE. (2000). Sequential mRNA expression for immediate early genes, cytokines, and neurotrophins in spinal cord injury.J. Neurotrauma.17203–218. 10.1089/neu.2000.17.203
41
HenniganA.O’callaghanR. M.KellyA. M. (2007). Neurotrophins and their receptors: roles in plasticity, neurodegeneration and neuroprotection.Biochem. Soc. Trans.35424–427. 10.1042/BST0350424
42
IrizarryR. A.HobbsB.CollinF.Beazer-BarclayY. D.AntonellisK. J.ScherfU.et al (2003). Exploration, normalization, and summaries of high density oligonucleotide array probe level data.Biostatistics4249–264. 10.1093/biostatistics/4.2.249
43
JinY.BouyerJ.HaasC.FischerI. (2015). Evaluation of the anatomical and functional consequences of repetitive mild cervical contusion using a model of spinal concussion.Exp. Neurol.271175–188. 10.1016/j.expneurol.2015.06.001
44
JonesT. B. (2014). Lymphocytes and autoimmunity after spinal cord injury.Exp. Neurol.25878–90. 10.1016/j.expneurol.2014.03.003
45
KangY.KimY. W.KangJ.KimA. (2021). Histone H3K4me1 and H3K27ac play roles in nucleosome eviction and eRNA transcription, respectively, at enhancers.FASEB J.35:e21781. 10.1096/fj.202100488R
46
KarinM.Ben-NeriahY. (2000). Phosphorylation meets ubiquitination: the control of NF-[kappa]B activity.Annu. Rev. Immunol.18621–663. 10.1146/annurev.immunol.18.1.621
47
KarinM.CaoY.GretenF. R.LiZ. W. (2002). NF-kappaB in cancer: from innocent bystander to major culprit.Nat. Rev. Cancer2301–310. 10.1038/nrc780
48
KimT. K.HembergM.GrayJ. M.CostaA. M.BearD. M.WuJ.et al (2010). Widespread transcription at neuronal activity-regulated enhancers.Nature465182–U165. 10.1038/nature09033
49
KohlM.WieseS.WarscheidB. (2011). Cytoscape: software for visualization and analysis of biological networks.Methods Mol. Biol.696291–303. 10.1007/978-1-60761-987-1_18
50
KonishiY.LehtinenM.DonovanN.BonniA. (2002). Cdc2 phosphorylation of BAD links the cell cycle to the cell death machinery.Mol. Cell91005–1016. 10.1016/S1097-2765(02)00524-5
51
LambJ.CrawfordE. D.PeckD.ModellJ. W.BlatI. C.WrobelM. J.et al (2006). The Connectivity Map: using gene-expression signatures to connect small molecules, genes, and disease.Science3131929–1935. 10.1126/science.1132939
52
LiJ. M.ZhangH.ZuoY. J. (2018). MicroRNA-218 alleviates sepsis inflammation by negatively regulating VOPP1 via JAK/STAT pathway.Eur. Rev. Med. Pharmacol. Sci.225620–5626.
53
LiW.HuY.OhS.MaQ.MerkurjevD.SongX.et al (2015). Condensin I and II Complexes License Full Estrogen Receptor alpha-Dependent Enhancer Activation.Mol. Cell59188–202. 10.1016/j.molcel.2015.06.002
54
LiW.NotaniD.RosenfeldM. G. (2016). Enhancers as non-coding RNA transcription units: recent insights and future perspectives.Nat. Rev. Genet.17207–223. 10.1038/nrg.2016.4
55
LiY. J.ZhangW.XiaH.ZhangB. S.ChenP.ZhaoY. L.et al (2017). miR-218 suppresses epithelial-to-mesenchymal transition by targeting Robo1 and Ecop in lung adenocarcinoma cells.Fut. Oncol.132571–2582. 10.2217/fon-2017-0398
56
LiberzonA.BirgerC.ThorvaldsdottirH.GhandiM.MesirovJ. P.TamayoP. (2015). The Molecular Signatures Database (MSigDB) hallmark gene set collection.Cell Syst.1417–425. 10.1016/j.cels.2015.12.004
57
ListerR.PelizzolaM.DowenR. H.HawkinsR. D.HonG.Tonti-FilippiniJ.et al (2009). Human DNA methylomes at base resolution show widespread epigenomic differences.Nature462315–322. 10.1038/nature08514
58
LiuB.NeufeldA. H. (2007). Activation of epidermal growth factor receptors in astrocytes: from development to neural injury.J. Neurosci. Res.853523–3529. 10.1002/jnr.21364
59
LoweryL. A.RubinJ.SiveH. (2007). Whitesnake/sfpq is required for cell survival and neuronal development in the zebrafish.Dev. DYN2361347–1357. 10.1002/dvdy.21132
60
MaynardF. M.BrackenM. B.CreaseyG.DitunnoJ. F.DonovanW. H.DuckerT. B.et al (1997). International standards for neurological and functional classification of spinal cord injury.Spinal. Cord35266–274. 10.1038/sj.sc.3100432
61
MurphyK. M.ReinerS. L. (2002). The lineage decisions of helper T cells.Nat. Rev. Immunol.2933–944. 10.1038/nri954
62
NakaeA.NakaiK.YanoK.HosokawaK.ShibataM.MashimoT. (2011). The animal model of spinal cord injury as an experimental pain model.J. Biomed. Biotechnol.2011:939023. 10.1155/2011/939023
63
NewmanA. M.LiuC. L.GreenM. R.GentlesA. J.FengW.XuY.et al (2015). Robust enumeration of cell subsets from tissue expression profiles.Nat. Methods12453–457. 10.1038/nmeth.3337
64
ParkS.JamesC. 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.Oncogene242495–2502. 10.1038/sj.onc.1208496
65
PhillipsA.MatinN.WestC.ZhengA.GaleaL.DorranceA.et al (2015). Autonomic Dysreflexia Impairs Cerebrovascular Health and Cognition in Experimental Spinal Cord Injury.FASEB J.2959–70. 10.1096/fasebj.29.1_supplement.800.10
66
PineauI.LacroixS. (2007). Proinflammatory cytokine synthesis in the injured mouse spinal cord: multiphasic expression pattern and identification of the cell types involved.J. Comp. Neurol.500267–285. 10.1002/cne.21149
67
PouwM. H.HosmanA. J.Van KampenA.HirschfeldS.ThietjeR.Van De MeentH. (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. Cord49307–312. 10.1038/sc.2010.114
68
RobinsonJ. T.ThorvaldsdottirH.WincklerW.GuttmanM.LanderE. S.GetzG.et al (2011). Integrative genomics viewer.Nat. Biotechnol.2924–26. 10.1038/nbt.1754
69
RobsonJ. P.WagnerB.GlitznerE.HeppnerF. L.SteinkellnerT.KhanD.et al (2018). Impaired neural stem cell expansion and hypersensitivity to epileptic seizures in mice lacking the EGFR in the brain.FEBS J.2853175–3196. 10.1111/febs.14603
70
RomanoR.BucciC. (2020). Role of EGFR in the Nervous System.Cells9:1887. 10.3390/cells9081887
71
RopperA. E.ZengX.AndersonJ. E.YuD.HanI.HaragopalH.et al (2015). An efficient device to experimentally model compression injury of mammalian spinal cord.Exp. Neurol.271515–523. 10.1016/j.expneurol.2015.07.012
72
RudmanM. D.ChoiJ. S.LeeH. E.TanS. K.AyadN. G.LeeJ. K. (2018). Bromodomain and extraterminal domain-containing protein inhibition attenuates acute inflammation after spinal cord injury.Exp. Neurol.309181–192. 10.1016/j.expneurol.2018.08.005
73
RyuH.LeeJ.OlofssonB. A.MwidauA.DedeogluA.EscuderoM.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.1004281–4286. 10.1073/pnas.0737363100
74
SelvarajahS.HammondE. R.SchneiderE. B. (2015). Trends in Traumatic Spinal Cord Injury.JAMA314:1643. 10.1001/jama.2015.11194
75
Shav-TalY.ZiporiD. (2002). PSF and p54(nrb)/NonO–multi-functional nuclear proteins.FEBS Lett.531109–114. 10.1016/S0014-5793(02)03447-6
76
SmithE.O’reillyA.SynnottK.MorrisS.TimlinM. (2015). Review of Time to Surgical Decompression in Traumatic Spinal Cord Injured Patients.Ir. Med. J.108265–267.
77
SmythG. K. (2004). Linear models and empirical bayes methods for assessing differential expression in microarray experiments.Stat. Appl. Genet. Mol. Biol.3:Article3. 10.2202/1544-6115.1027
78
SubramanianA.KuehnH.GouldJ.TamayoP.MesirovJ. P. (2007). GSEA-P: a desktop application for Gene Set Enrichment Analysis.Bioinformatics233251–3253. 10.1093/bioinformatics/btm369
79
WangH.ZouJ.ZhaoB.JohannsenE.AshworthT.WongH.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.10814908–14913. 10.1073/pnas.1109023108
80
WangX. Y.SmithD. I.FrederickL.JamesC. D. (1998). Analysis of EGF receptor amplicons reveals amplification of multiple expressed sequences.Oncogene16191–195. 10.1038/sj.onc.1201476
81
WangY.WangR.ZhangS.SongS.JiangC.HanG.et al (2019). iTALK: an R Package to Characterize and Illustrate Intercellular Communication.bioRxiv[preprint]. 10.1101/507871
82
WestJ. A.DavisC. P.SunwooH.SimonM. D.SadreyevR. I.WangP. I.et al (2014). The long noncoding RNAs NEAT1 and MALAT1 bind active chromatin sites.Mol. Cell.55791–802. 10.1016/j.molcel.2014.07.012
83
WhyteW. A.OrlandoD. A.HniszD.AbrahamB. J.LinC. Y.KageyM. H.et al (2013). Master transcription factors and mediator establish super-enhancers at key cell identity genes.Cell153307–319. 10.1016/j.cell.2013.03.035
84
WolfF. A.AngererP.TheisF. J. (2018). SCANPY: large-scale single-cell gene expression data analysis.Genome Biol.19:15. 10.1186/s13059-017-1382-0
85
XiaoB.LiuL.LiA.XiangC.WangP.LiH.et al (2020). Identification and Verification of Immune-Related Gene Prognostic Signature Based on ssGSEA for Osteosarcoma.Front. Oncol.10:607622. 10.3389/fonc.2020.607622
86
YaoJ. S. (1992). Deep vein thrombosis in spinal cord-injured patients. Evaluation and assessment.Chest102645S–648S. 10.1378/chest.102.6_Supplement.645S
87
YuG.WangL. G.HanY.HeQ. Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters.OMICS16284–287. 10.1089/omi.2011.0118
88
ZhangT.JiangM.ChenL.NiuB.CaiY. (2013). Prediction of gene phenotypes based on GO and KEGG pathway enrichment scores.Biomed. Res. Int.2013:870795. 10.1155/2013/870795
89
ZhangX.LanY.XuJ.QuanF.ZhaoE.DengC.et al (2019). CellMarker: a manually curated resource of cell markers in human and mouse.Nucleic Acids Res.47D721–D728. 10.1093/nar/gky900
90
ZhangY.GuanZ.ReaderB.ShawlerT.Mandrekar-ColucciS.HuangK.et al (2013). Autonomic dysreflexia causes chronic immune suppression after spinal cord injury.J. Neurosci.3312970–12981. 10.1523/JNEUROSCI.1974-13.2013
91
ZhengR.WanC.MeiS.QinQ.WuQ.SunH.et al (2019). Cistrome Data Browser: expanded datasets and new tools for gene regulatory analysis.Nucleic Acids Res.47D729–D735. 10.1093/nar/gky1094
92
ZhuY.SoderblomC.KrishnanV.AshbaughJ.BetheaJ. R.LeeJ. K. (2015). Hematogenous macrophage depletion reduces the fibrotic scar and increases axonal growth after spinal cord injury.Neurobiol. Dis.74114–125. 10.1016/j.nbd.2014.10.024
Summary
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
Volume
9 - 2021
Edited by
Chuan-Xing Li, Karolinska Institutet (KI), Sweden
Reviewed by
Zhi-an Cheng, Guangdong Provincial Hospital of Chinese Medicine, China; Lixin Cheng, Jinan University, China
Updates
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, gzhuangzq@163.comLiming Cheng, limingcheng@tongji.edu.cnJie Zhang, jiezhang@tongji.edu.cn
†These authors have contributed equally to this work
This article was submitted to Molecular and Cellular Pathology, a section of the journal Frontiers in Cell and Developmental Biology
Disclaimer
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.