A transcriptomic study of probenecid on injured spinal cords in mice

Background Recent studies have found that probenecid has neuroprotective and reparative effects on central nervous system injuries. However, its effect on genome-wide transcription in acute spinal cord injury (SCI) remains unknown. In the present study, RNA sequencing (RNA-Seq) is used to analyze the effect of probenecid on the local expression of gene transcription 8 h after spinal injury. Methods An Infinite Horizon impactor was used to perform contusive SCI in mice. The SCI model was made by using a rod (1.3 mm diameter) with a force of 50 Kdynes. Sham-operated mice only received a laminectomy without contusive injury. The injured mice were randomly assigned into either the control (SCI_C) or probenecid injection (SCI_P) group. In the latter group, the probenecid drug was intraperitoneally injected (0.5 mg/kg) immediately following injury. Eight hours after the injury or laminectomy, the spinal cords were removed from the mice in both groups. The total RNAs were extracted and purified for library preparation and transcriptome sequencing. Differential gene expressions (DEGs) of the three groups—sham, SCI_C and SCI_P—were analyzed using a DESeq software. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of DEGs were performed using a GOseq R package and KOBAS software. Real-time quantitative reverse-transcriptase polymerase chain reaction was used to validate RNA-Seq results. Results RNA-Seq showed that, compared to the SCI_C group, the number of DEGs was 641 in the SCI_P group (286 upregulated and 355 downregulated). According to GO analysis, DEGs were most enriched in extracellular matrix (ECM), collagen trimer, protein bounding and sequence specific DNA binding. KEGG analysis showed that the most enriched pathways included: cell adhesion molecules, Leukocyte transendothelial migration, ECM-receptor interactions, PI3K-Akt signaling pathways, hematopoietic cell lineages, focal adhesions, the Rap1 signaling pathway, etc. The sequence data have been deposited into the Sequence Read Archive (https://www.ncbi.nlm.nih.gov/sra/PRJNA554464).


INTRODUCTION
The term spinal cord injury (SCI) refers to a variety of injuries to the spinal cord. According to the severity of injury the symptoms may vary, ranging from pain to a complete loss of movement and sensory function. SCI affects millions of people worldwide, typically for life (Friedli et al., 2015). In the United States, there are around 12,000-20,000 new SCI cases a year, with more than 280,000 patients confined to wheelchairs (Singh et al., 2014). In the past decade, the SCI cases in China have risen tenfold, and are now increasing by about 60,000 cases every year (Qiu, 2009). SCI has a high rate of disability and mortality, which creates a heavy burden for patients, their families and society (Krueger et al., 2013). Therefore, it is imperative to explore effective treatment methods for repairing SCI in order to improve the quality of life of patients and reduce the burden of social medical care.
Pathological processes following traumatic SCI can be characterized as primary and secondary injuries (Geisler et al., 2002;McDonald & Sadowsky, 2002). Primary injury refers to the direct injury of the spinal cord by mechanical force, including compression, contusion, laceration and penetration. Secondary injury refers to edema, ischemia, local inflammation and electrolyte changes. These changes can cause an accumulation of lipid peroxides and oxygenfree radicals, as well as a release of inflammatory factors and proteases. This can ultimately lead to a large amount of cell apoptosis or necrosis, which further aggravates the damage to the neurons and axons (Ahuja et al., 2017;Oyinbo, 2011;Tran, Warren & Silver, 2018).
Probenecid is an organic anion-transport protein inhibitor, which has been widely used in clinical settings (Hagos et al., 2017;Töllner et al., 2015). For example, probenecid has been used as a synergist in the treatment of gout and antibiotics (Baranova et al., 2004;Papadopoulos & Verkman, 2008). It can reduce the degree of cognitive impairment in afflicted rats (Mawhinney et al., 2011), as well as reverse cerebral ischemic injury and cellular inflammation (Wei et al., 2015;Xiong et al., 2014). The combination of probenecid and N-Acetylcysteine could potentially both maintain intracellular GSH concentrations and inhibit neuronal death after a traumatic stretch injury (Du et al., 2016). Some studies report that probenecid can also reduce neuropathic pain in the spinal cord (Bravo et al., 2014;Pineda-Farias et al., 2013). Therefore, these reports indicate that probenecid has neuroprotective and reparative effects on central nervous system (CNS) injuries. However, whether the drug can play a role in treating SCI and whether it can affect the gene expression profiles in injured spinal cords remain unknown. To test this, probenecid was injected intraperitoneally into spinal cord-injured mice immediately after injury. Eight hours after the injury or laminectomy, the spinal cords were removed and RNA-Seq was used to analyze the changes in transcriptome expression in the injured area. Then, the key molecules and signal pathways were screened and identified, providing a new theoretical and experimental basis for SCI clinical treatment.

Animals
A total of 27 healthy, clean C57BL/6 female mice (18-20 g, 8 weeks old) were used to model SCI. The Animal Care and Use Committee of Bengbu Medical College provided full approval for this research (037/2017). Animal care following surgery was carried out in compliance with the regulations for the management of experimental animals (revised by the Ministry of Science and Technology of China in June 2004), as well as the guidelines and policies on rodent survival surgery provided by the Animal Care and Use Committee of Bengbu Medical College.

Contusive SCI and drug injection
An Infinite Horizon impactor (Precision Systems & Instrumentation, Lexington, KY, USA) was used to perform contusive SCI on the mice. These mice were first anesthetized with 50 mg/kg pentobarbital, followed by the excision of the T9 lamina. A SCI model of this procedure was created using a rod (1.3 mm diameter) with a force of 50 Kdynes. Sham-operated (sham) mice only received a laminectomy without contusive injury.
The spinal cord-injured mice were randomly assigned to the solvent control (SCI_C) or probenecid injection (SCI_P) group. The solvent or probenecid (0.5 mg/kg) was intraperitoneally injected immediately following injury. The solution (pH 7.3) was prepared as previously described (Hainz et al., 2017).

RNA isolation, quantification and qualification
Eight hours after the injury or laminectomy, the mice were anesthetized and perfused with 10 ml PBS. Their spinal cords (0.5 cm including the injury center) were then removed. The total RNAs from their spinal cords were extracted and purified as previously described (Shi et al., 2017).

Library preparation and transcriptome sequencing
The sequencing libraries were produced by using a NEBNext Ò Ultra TM RNA Library Prep Kit for Illumina Ò (New England Biolabs, Ipswich, MA, USA) as previously described (Shi et al., 2017). Finally, the 125 bp/150 bp paired-end reads were obtained and sequenced on an Illumina Hiseq platform.

Analysis of differentially expressed gene
Prior to DEG analysis, the gene expression statistics were analyzed using RSEM software (http://deweylab.github.io/RSEM/) to convert the read count numbers to Fragments Per Kilobase of transcript per Million fragments mapped (FPKM), and Principal Component Analysis (PCA) analysis was conducted to determine the similarities and differences in the data. DEGs of the three groups of mice were analyzed as previously described (Shi et al., 2017) using DESeq software (http://www.bioconductor.org/). Benjamini and Hochberg's approach was used to control the false discovery rate and adjust the P-values. The adjusted P-value < 0.05 was defined as a standard for significant differences in gene expression. In addition to FPKM hierarchical clustering analysis of DEGs, we further analyzed the subclusters based on log2 (ratios) of their gene expression level relative to that of sham group. The log2 (ratios) in the SCI_C group ≥1 or ≤−1 was used as a cutoff for subcluster analysis. The clustering algorithm divided the DEGs with similar gene expression trends into several subclusters.

Gene ontology and kyoto encyclopedia of genes and genomes enrichment analysis of DEGs
The GO and KEGG analysis were performed using a GOseq R package and KOBAS software as previously described (Shi et al., 2017). In GO analysis, DEGs were implemented using the GOseq R package and gene length bias was corrected. GO terms with corrected P value ≤ 0.05 were considered significantly enriched by DEGs. KEGG is a database resource for understanding the high-level functions and utilities of the biological system (http://www.genome.jp/kegg/). In this study, we used KOBAS software to test the statistical enrichment of DEGs in KEGG pathways.

Real-time quantitative reverse-transcriptase polymerase chain reaction
To validate RNA-Seq results, nine DEGs were randomly selected and verified via Real-time quantitative reverse-transcriptase polymerase chain reaction (RT-qPCR) according to our previous methods (Shi et al., 2017). The analysis was performed in six samples, which included three independent samples and duplicates of these samples to be used in RNA-seq analysis. PCR primer sequences are listed in Table 1. The relative quantitative results of each group of genes were calculated according to the formula " ΔΔ Ct" (Livak & Schmittgen, 2001). The statistical values (n = 6/group) were presented as mean ± standard deviation (SD). The data were analyzed using one-way Analysis of Variance (ANOVA), followed by Student-Newman-Keuls tests. Statistical differences were considered significant at P < 0.05.

Identification of expressed transcripts the mice spinal cords
For the high quality assessment of sequencing data, nine cDNA libraries were established, including sham (sham_1, sham_2 and sham_3), SCI_C (SCI_C1, SCI_C2 and SCI_C3) and SCI_P (SCI_P1, SCI_P2 and SCI_P3). RNA-Seq produced 48,848,744-61,037,096 raw reads for each sample. After filtering out the low-quality reads, there were 48,226,002-60,037,772 clean reads, with the Q30 (%) 93.67-94.31 (Table 2). In order to identify the source of variation within the original data, PCA analysis was conducted. As shown in Fig. 1, PC1, PC2 and PC3 were 54.51%, 12.33% and 7.09%, respectively. Although not too far from one another, the distance between SCI_C (or SCI_P) and sham was apparent and sufficient for the analysis. These distances demonstrated that the data could be used for the next analysis.
Effect of SCI and probenecid treatment on gene expression RPKM and DEGSeq were used to analyze the gene expression level and differential expression profiles, respectively. The results showed that, as compared to the sham group,  there were 4,617 DEGs in the SCI_C group, including 2,904 upregulated and 1,713 downregulated genes ( Fig. 2A; Table S1). Compared to the SCI_C group, there were 641 different genes in the SCI_P group, 286 were upregulated and 355 were downregulated ( Fig. 2B; Table S1). The sequence data have been deposited into Sequence Read Archive (https://www.ncbi.nlm.nih.gov/sra/PRJNA554464).

RT-qPCR identification of DEGs
In order to verify the RNA-Seq results, nine DEGs were randomly selected from the SCI_P group, as compared with the SCI_C group, namely Itga1, Lamb1, Cldn5, Lama2, CD34, Esam, Setdb2, Agrn and Ccnt2. The RNA-Seq and RT-qPCR results indicated that the expression patterns of these DEGs were similar (Fig. 3).

Cluster analysis of DEGs
The DEGs from different groups were analyzed using FPKM hierarchical cluster analysis. As shown in Fig. 4, DEGs were hierarchically clustered and classified into different expression clusters. These clusters contained upregulated or downregulated DEGs. When compared to the sham group, most upregulated DEGs in the SCI_C group were in the middle and upper clusters, while downregulated DEGs were delegated to the lower cluster.
Additionally, compared to the sham group, most upregulated DEGs in the SCI_P group were in the upper cluster, while downregulated DEGs were mainly observed in the lower cluster. When compared to the SCI_C group, some upregulated DEGs in the SCI_P group were observed in upper cluster, while downregulated DEGs were observed in the middle cluster (there were also some clusters in this grouping that showed no significant differences). In addition to FPKM hierarchical clustering analysis of DEGs, the subclusters-which have similar expression trends-were further analyzed. The log2 (ratios) in the SCI_C group ≥1 or ≤−1 was used as a cutoff for subcluster analysis. As shown in Fig. 5, we found several subclusters with similar expression trends. Based on log2 (ratios) of their gene expression levels relative to that of the sham group, the log2 (ratios) of all gene expression levels in the sham group were zero. Figures 5A and 5B show that the two subclusters were strongly upregulated following SCI and then downregulated upon probenecid treatment. The longitudinal coordinates in RT-qPCR were the mRNA expression level calculated using the ΔΔ Ct method and expressed relative to the value in the sham group (designated as oness). All data were calculated with mean ± standard deviation (n = 6, which included three independent samples and the three same samples used for the RNA-seq analysis). ÃÃ P < 0.01 (ANOVA). Figures 5C and 5D show that the two subclusters were strongly downregulated following SCI and then upregulated upon probenecid treatment. In Fig. 5A, six genes (Cybb, Esam, Itgam, Itgb2, Msn and Ncf2) are demonstrated to have been involved in the leukocyte transendothelial migration signaling pathway; six genes (Col4a1, Erbb2, Flt4, Nos3, Syk and Thbs4) were also involved in the PI3K-Akt signaling pathway. Figure 5B displays three genes (Cyba, Ncf1 and Rac2) involved in the NADPH oxidases, two genes (Cflar and Tnfrsf10b) involved in the TRAIL signaling pathway and eight genes (Cd63, Cyba, Ddx58, Fcer1g, Lyn, Myh9, Ncf1 and Psmb8) involved in the innate immune system. Figures 5C and 5D show that no gene can be clustered into valuable signaling pathways.

GO enrichment analysis of DEGs
When compared with the sham group, there were 78 GO terms in upregulated DEGs ( Fig. 6A; Table S2) and nine GO terms in downregulated DEGs ( Fig. 6B; Table S2) in the SCI_C group. The upregulated DEGs were most enriched in: binding, protein binding, chemokine activity, chemokine receptor binding, G-protein coupled receptor binding, anion binding, small GTPase mediated signal transduction, immune system process, immune response. The downregulated DEGs were most enriched in: protein binding, binding, extracellular-glutamate-gated ion channel activity, acid phosphatase activity, transporter activity, mannose metabolic process, excitatory extracellular ligand-gated ion channel activity, transmembrane transporter activity, anion transmembrane and transporter activity. In the SCI_P group, we observed three GO terms in downregulated DEGs ( Fig. 6C; Table S3) and no valuable terms in upregulated DEGs (Table S3) compared to the SCI_C group. The downregulated DEGs were protein binding, binding and sequence-specific DNA binding.

KEGG enrichment analysis of DEGs
Scatter plots were used to express the KEGG enrichment analysis results for the DEGs. When compared to the sham group, the upregulated DEGs in the SCI_C group were most enriched in TNF, NF-kappa B, cytokine-cytokine receptor interaction, Toll-like receptor, Leukocyte transendothelial migration, PI3K-Akt, focal adhesion and apoptosis ( Fig. 7A; Table S4); the downregulated DEGs were most enriched in glutamatergic synapse, basal cell carcinoma, axon guidance, other glycan degradation and nicotine addiction ( Fig. 7B; Table S4). In the SCI_P group vs. SCI_C group, only the "ECMreceptor interaction" was enriched in the upregulated DEGs ( Fig. 7C; Table S5); the downregulated DEGs were enriched in CAMs, malaria, leukocyte transendothelial migration, ECM-receptor interaction, PI3K-Akt signaling pathway, hematopoietic cell lineage, focal adhesion, Rap1 signaling pathway and amebiasis ( Fig. 7D; Table S5).

DISCUSSION
Recent studies have shown that probenecid has neuroprotective and repairing effects on brain disorders (Wei et al., 2015;Xiong et al., 2014). However, its effect on genome-wide transcription in SCI is still unknown. Therefore, in this study, RNA-Seq was used to analyze the effect of probenecid on the local expression of gene transcription 8 h after SCI.
The results showed that when compared with the sham group, there were 4,617 DEGs in the SCI_C group, including 2,904 upregulated and 1,713 downregulated genes. When compared with the SCI_C group, there were 641 DEGs in the SCI_P group, 286 of which were upregulated and 355 downregulated. These are consistent with others' and our previous reports (Chen et al., 2013;Shi et al., 2017). It also shows that the results of this experiment are reliable. As compared to the SCI_C, there were 641 DEGs in the SCI_P group, 286 were upregulated and 355 were downregulated. To further verify the RNA-seq results, we randomly selected nine DEGs (Itga1, Lamb1, Cldn5, Lama2, CD34, Esam, Setdb2, Agrn and Ccnt2) for RT-qPCR. The results showed that the expression patterns of these genes detected by these two methods were similar, indicating that our RNA-seq results are reliable and can be used for subsequent analysis. These also confirmed that probenecid can alter gene transcription after SCI.
To further analyze the DEGs effected by probenecid, we used GO enrichment to reflect the distribution of DEGs on GO term enriched in cell components, molecular functions and biological processes (Huang et al., 2013). In the SCI_P vs. SCI_C group, analysis showed that there were three GO terms in downregulated DEGs (protein binding, binding and sequence-specific DNA binding) and no valuable terms in upregulated DEGs. KEGG analysis showed that the valuable signaling pathways associated with these DEGs included: CAMs, leukocyte transendothelial migration, ECM-receptor interaction, PI3K-Akt signaling pathway, hematopoietic cell lineage, focal adhesion and Rap1 signaling pathway.
Following SCI, probenecid treatment could downregulate some genes, subclusters and signaling pathways. Leukocyte transendothelial migration from the blood into tissues is vital for immune surveillance and inflammation (Cook-Mills, 2006). However there is a large amount of leukocyte infiltration involved in the pathological process of SCI. These infiltrating leukocytes need to bind to endothelial CAM and then migrate between vascular endothelial cells (Wang et al., 2011). Therefore, the inhibition of leukocyte transendothelial migration and CAMs induced by probenecid may play a role in inhibiting inflammation by weakening the infiltration of white blood cells in the injured area. In this study, we clustered six genes (Cybb,Esam,Itgam,Itgb2,Msn and Ncf2) involved in this pathway. Their expression is strongly downregulated following SCI and then upregulated upon probenecid treatment. This demonstrates that probenecid treatment following SCI can play an anti-inflammatory role by inhibiting the infiltration of inflammatory cells.
The ECM plays an important role in tissue and organ morphogenesis (Bonnans, Chou & Werb, 2014;Rabelink et al., 2017) and control of cellular activities such as adhesion, migration, differentiation, proliferation and apoptosis (Yue, 2014). Focal Figure 7 (continued) the P value after multiple hypothesis test correction. The range of Q value is between 0 and 1. The closer the Q value is to 0, the more significant the enrichment is. The KEGG pathways were shown in (A) upregulated DEGs (SCI_C vs. sham); (B) downregulated DEGs; (C) upregulated DEGs (SCI_P vs. SCP_C); (D) downregulated DEGs (SCI_C vs. sham).
Full-size  DOI: 10.7717/peerj.8367/ fig-7 adhesions are specialized, intracellular sites in which aggregated integrin receptors interact with extracellular matrices, while extracellular matrices interact with intracellular actin cytoskeleton (Burridge, 2017;LaFlamme et al., 2018). At the same time, focal adhesions are the result of cell and extracellular matrix (ECM) interactions (Burridge, 2017;De Pascalis & Etienne-Manneville, 2017). ECM and focal adhesions are downregulated after probenecid treatment, indicating that probenecid might improve SCI by inhibiting adhesion, migration, differentiation, proliferation and apoptosis. It has been reported that PI3K-Akt signaling fuzes a variety of extracellular and intracellular signal transduction pathways that regulate macrophage biology, including: the production of pro-inflammatory cytokines, phagocytosis, autophagy and homeostasis (Vergadi et al., 2017). PI3K-Akt signal pathway is downregulated in SCI after probenecid treatment, and the six genes (Col4a1, Erbb2, Flt4, Nos3, Syk and Thbs4) being clustered into this pathway indicate that probenecid might improve SCI by regulating macrophages and inhibiting inflammatory pathways. This is likely to provide important clues towards identifying the mechanism by which probenecid acts.
The relationship between the hematopoietic cell lineage pathway and SCI was referred to in a report on the bioinformatics analysis of SCI (Zhu et al., 2017). Its specific role has yet to be reported and requires further discussion.
Rap1 signal pathway plays an important role in regulating cell-cell and cell-matrix interactions by regulating the function of adhesion molecules (Kim, Ye & Ginsberg, 2011;Pollan et al., 2018). In our study, Rap1 signaling pathways were enriched in downregulated DEGs of SCI after probenecid treatment, suggesting that probenecid may inhibit cell adhesion and polarization by inhibiting the Rap1 signaling pathway, thereby inhibiting inflammation.
In addition, three genes (Cyba, Ncf1 and Rac2) related to the NADPH oxidases, two genes (Cflar and Tnfrsf10b) related to the TRAIL signaling pathway and eight genes (Cd63, Cyba, Ddx58, Fcer1g, Lyn, Myh9, Ncf1 and Psmb8) related to the innate immune system were also strongly downregulated following probenecid treatment. We know that NADPH oxidases are involved in oxidative stress, the TRALL signaling pathway mediates inflammation and apoptosis and the immune system is involved in almost all pathological processes of injury (Chyuan et al., 2018;Ewald, 2018;Tisato et al., 2018). Therefore, probenecid treatment can potentially play a neuroprotective role by inhibiting immune response, oxidative stress, anti-inflammation and anti-apoptosis after SCI.

CONCLUSIONS
Acute SCI can lead to changes in the mRNAs of injured spinal cords. These mRNAs and their related pathways can provide some explanation for the pathological mechanism of acute SCI. More interestingly, we also demonstrated that probenecid can lead to gene expression inhibitions in an acute injured spinal cord. These downregulated DEGs and their associated signaling pathways-such as focal adhesions, leukocyte transendothelial migration, ECM-receptor interaction, PI3K-Akt, Rap1-are mainly related to inflammatory response, local hypoxia, macrophage differentiation, adhesion migration and apoptosis of local cells. This suggests that the application of probenecid in the acute phase can improve the local microenvironment of SCI. However, the application of probenecid as a therapeutic drug for SCI requires further investigation. Next, the detailed research on this subject will be conducted by combining animal models and clinical practice.

ADDITIONAL INFORMATION AND DECLARATIONS Funding
This work was supported by grants from the National Natural Science Foundation of China (Nos. 81571194 and 81772321). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Grant Disclosures
The following grant information was disclosed by the authors: National Natural Science Foundation of China: 81571194 and 81772321.