A transcriptomic study of probenecid on injured spinal cords in mice

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™ 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. 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 “ECM-receptor 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).