Ketogenic diet-mediated steroid metabolism reprogramming improves the immune microenvironment and myelin growth in spinal cord injury rats according to gene and co-expression network analyses

Experimental procedures and animal conditions

The entire experimental process is shown in Figure 1A. To establish a C7 spinal cord hemi-contusion model, unilateral C7 laminectomies were performed using a Precision Systems and Instrumentation-IH Impactor according to the manufacturer’s instructions (Supplementary Figure 1A). Sprague-Dawley rats were anaesthetized with sodium pentobarbital to expose the spinal cord. We fixed the C6-T1 vertebrae rigidly in a frame tilted at a 25.0° angle and then used the Infinite Horizon impactor to deliver a set contusion force of 150 kdyne (Supplementary Figure 1B); the contusion forces were monitored by a computer (Supplementary Figure 1C). An average of 150 kdyne contusion force was obtained for each sample. These results indicated that the C7 spinal cord hemi-contusion model was successful. Based on myelin staining, the myelin sheath area in the SCI_KD group was significantly larger than that in the spinal cord injury with standard diet (SCI_SD) group (Figure 1B–1C), indicating that SCI rats exhibited protection or promotion of myelin sheath growth at the injured site after feeding with a KD.

Figure 1. Experimental protocol and animal myelin area. (A) Detailed experimental protocol. (B) Representative images showing EC staining after SCI. Scale bar = 1 mm. (C) Quantitative analysis of the results in panel (A). All data are presented as the mean ± standard deviation, n = 3. Rostral of the injury site, p = 0.0255; Epicentre, p = 0.0004; Caudal, p = 0.0015.

Differentially expressed mRNAs

To further understand the effects and mechanism of the KD on SCI, RNA-seq was performed on each group with 5 biological replicate samples, in which 5 samples were degraded and eliminated. To explore transcription factor expression in Sprague-Dawley rats treated with the KD under spinal cord injury, we established 15 RNA-seq cDNA libraries from the rats divided into 4 groups (SCI_KD, SCI_SD, Sham_KD, and Sham_SD). A summary of the sequence assembly after Illumina sequencing is presented in Supplementary Table 1. The average of raw reads was 23295214.7, the Q30 basic mass fraction of the sample was 100%, and the error rate (%) was 0, meeting the requirements for subsequent analysis. Using a comparison of the transcriptome and reference genome, greater than 82% of reads were uniquely mapped; these results met the subsequent analysis requirements (Supplementary Table 2). In short, the sequencing data were suitable for subsequent data analysis.

A total of 32,883 genes were detected by RNA-seq, and Venn analysis between samples showed the number of expressed genes in each group (Figure 2A). We first used an unsupervised classification method-principal component analysis (PCA) to characterize the differences among the gene expression profiles of the 15 samples. Figure 2C shows that the biologically replicated samples could be divided into the following 4 groups: Sham_SD1-3, Sham_KD1-4, SCI_SD1-3 and SCI_KD1-5. A PCA model showed that the maximum variation in PC1 could explain 16.92% of the variation, while principal component 2 explained 9.85% of the variation [11]. In Figure 2B, the scatter plot of differential genes between the SCI_KD and SCI_SD groups shows that 463 genes were differentially expressed between the two groups, among which 184 DEGs were upregulated, and 279 DEGs were downregulated. The top 50 upregulated and downregulated genes are listed in Supplementary Tables 3 and 4, respectively.

Figure 2. Changes in mRNA expression profiles between groups. (A) The mRNA expression of the four groups in the entire transcriptome. (B) Volcano plots showing the up- and downregulated mRNA transcripts in the SCI_KD vs. the SCI_SD group. (C) Principal component analysis (PCA) score plot with the two principal components (PC1 = 16.92%; PC2 = 9.85%).

Reprogramming steroid biosynthesis and metabolism in SCI rats on KD

Compared with the SCI_SD group, the KD group exhibited upregulated DEGs in many enrichment pathways, such as the steroid synthesis pathway in lipid metabolism (path ID: map00100, q-value = 5.82043E-08), biosynthesis of unsaturated fatty acids, terpenoid backbone biosynthesis, etc. A bubble chart of the DEGs in the top 20 most enriched pathways is provided (Figure 3A). The downregulated DEGs from the KD diet group were significantly enriched in ribosomes, lysosomes, osteoclast differentiation, proteasome, tuberculosis, and antigen processing and display. Figure 3B shows the KEGG enrichment distribution of the DEGs.

Figure 3. Enrichment analysis of DEGs in the SCI_KD vs. the SCI_SD group. (A) KEGG enrichment bubble diagram of the upregulated DEGs from the SCI_KD and SCI_SD groups. (B) KEGG enrichment bubble diagram of the downregulated DEGs from the SCI_KD and SCI_SD groups. (C) The GSEA analysis diagram of steroid biosynthesis between SCI_KD and SCI_SD. (D) Cluster analysis of KEGG-enriched DEGs. (E) Protein-protein interaction of DEGs between the SCI_KD and SCI_SD groups.

To further understand the biological functions of the SCI_SD and SCI_KD DEGs, we performed gene set enrichment analysis (GSEA) of the sample sets based on the C2 MSigDB database, which includes the Broad Institute recombined gene sets. The selected gene set was tested, and it was found that SCI_KD was significantly negatively correlated with inflammatory cytokines and their receptors (the NES values were ˗1.44 and ˗1.46, respectively; Supplementary Figure 2), suggesting that a KD reduces the lethality of SCI, which echoes the results of our previous studies. Similarly, based on the KEGG functional annotation, we found that SCI_KD was significantly positively correlated with the steroid synthesis pathway (the NES is 1.62; Figure 3C) by GSEA analysis, which further shows that the KD is involved in the metabolic reprogramming of steroids in SCI rats, which is a positive effect of the KD. An important mechanism of SCI function may be related to reducing metabolic responses. iPath3.0 was used to visually analyse the enriched steroid metabolism-related genes, as shown in Supplementary Figure 3, which shows that the steroid metabolism DEGs in the biological system transmit information.

Cluster analysis and PPI analysis of DEGs in the steroid synthesis pathway

Figure 3D shows that the Sqle (squalene epoxidase), Sc5d (sterol-C5-desaturase), Cyp51 (cytochrome P450, family 51), Dhcr24 (24-dehydrocholesterol reductase), Msmo1 (methylsterol monooxygenase 1), Hsd17b7 (hydroxysteroid (17-β) dehydrogenase 7), and Fdft1 (farnesyl diphosphate farnesyltransferase 1) genes were obviously highly expressed in the SCI_KD group. DEGs corresponded to interaction relationships for network integration and allowed us to identify key pathways in the interaction network according to indicators such as the connectivity between genes (Figure 3E).

WGCNA and hub module identification and analysis

According to the KEGG enrichment analysis of differential gene expression and the GSEA, the mean values showed that SCI rats in the KD group exhibited significantly increased steroid anabolism, which indicates that the KD may be protective by mediating steroid metabolism reprogramming. By changing certain physiological functions or pathways, we could analyse the SCI_KD, SCI_SD, Sham_KD, and Sham_SD results in a gene network (WGCNA) and fully understand the roles of a KD and standard diet (SD) on genetic changes in Sprague-Dawley rats. Genes that were commonly expressed in the four groups were placed in the same gene network through the WGCNA algorithm. A total of 32,883 mRNAs were screened by WGCNA, and a total of 5,163 mRNAs were used for subsequent analysis after pre-processing and filtering, such as excluding deletions and outliers. Figure 4A and 4B show the scale-free fitness curve and the average connectivity curve, respectively. The scale-free fitness curve is in a smooth position, and the power exponent is a weighted β value; thus, the soft threshold value is set to 6, and the scale-free topological index is 0.9. Therefore, the network conforms to the power law distribution and is close to the actual biological network state. The resulting gene tree diagram and the colours in the corresponding modules are shown in Supplementary Figure 4A–4B, with the number of mRNAs for each module shown in Supplementary Table 5.

Figure 4. Data pre-processing and module identification for WGCNA analysis. (A–B) Analysis of the network topology with different soft threshold powers. The left panel shows the influence of soft threshold power on the scale-free topology fitting index. The right panel shows the influence of soft threshold power on the average connectivity. (C) Analysis of the correlation between the module and the diet phenotype. Red indicates that the module has a higher correlation with the phenotype, while blue means the module has a lower correlation with the phenotype. (D) WGCNA module correlation; the colour indicates the correlation between modules.

After the data were pre-processed, the genes were classified, and genes with similar expression patterns were divided into categories called modules. We performed module clustering and module correlation analysis for the obtained modules, as shown in Figure 4D. Interaction relationships were analysed in the six different modules, and a network heat map of the correlation between the modules and the phenotype was drawn (Figure 4C). The results show that the modules are independent of one another, which indicates the relative independence of gene expression in different modules. Here, we mainly analysed the correlation modules between the KD- and SD-fed rats. Among them, the MEgreen module is highly related to a KD, and the MEblue module is highly related to a SD. We then calculated the correlation coefficient between each gene in the module and specific phenotype data to obtain the gene significance (GS) value; this value was combined with the MM value of each gene in the module (i.e., kME, eigengene-based connectivity) to obtain the MM-GS scatter plot. The stronger the correlation, the closer the relationship between the module and the phenotype; moreover, correlations with a lower p-value are more reliable. This analysis shows that MEgreen is an important module in the KD, while MEblue is an important module in the SD (Supplementary Figure 4C–4F).

Functional annotation analysis of hub modular genes

To explore the biological functions of the KD and SD modules (green and blue, respectively) and compare the differences in their biological functions, we performed GO and KEGG enrichment analysis with annotation and visualization. GO function annotation analysis classifies DEGs using different perspectives, such as the biological processes involved, the components that constitute cells, and the molecular functions achieved. Figure 5A shows the GO annotation functions of the MEgreen module, and Figure 5B shows the GO annotation functions of the MEblue module. Different from other groups, there was a significant positive correlation between SCI rats that received a SD and the MEblue module. GO functional analysis of immune system processes is expressed as biological processes, which indicated that the immune system significantly affected SCI rats, and immune system processes included all processes of inflammation-related pathways and congenital and acquired immune expression changes.

Figure 5. GO function annotation analysis of WGCNA module genes. (A) The GO function annotation analysis histogram of the MEgreen module. (B) The GO function annotation analysis histogram of the MEblue module, which includes biological processes, cell components, and molecular functions.

Similarly, KEGG annotation analysis was also performed on the green and blue modules. The results showed that the KEGG annotation molecule in the MEblue module significantly identified immune system, immune system and infectious diseases: viruses and other immune disease pathways among the subtypes. This module mainly focused on rats with SCI that received a SD, which is consistent with the GO function analysis results. Compared with the MEblue module, in terms of expression ratio, immune system, immune system and infectious diseases with the sub-categories of viruses and other immune disease pathways were significantly downregulated, and the module is related to the KD (Supplementary Figure 5).

Multi-gene set KEGG enrichment analysis and gene interaction network analysis

The immune system and infectious diseases: viral and immune disease pathways were used as pathways of interest for multi-gene KEGG enrichment analysis. Figure 6A shows the KEGG enrichment results of the MEgreen module, and Figure 6B shows the KEGG enrichment results of the MEblue module. The MEgreen module is primarily a module closely related to the KD rats, while the MEblue module is mainly a module closely linked to the SD rats, especially the SD injury rats. The results of the KEGG enrichment signal pathways and intersection of the top 20 pathways for the two modules show that only the viral myocarditis and phagosome pathways are present in both modules. The MEgreen module mainly activates innate immune-related pathways, such as the Toll-like receptor signalling pathway and the NOD-like receptor signalling pathway. These results indicate that MEgreen and MEblue act on different immune pathways. The WGCNA and functional analysis showed that rats fed a KD showed a reduction in immune system pathway processes, including immune diseases, infectious diseases and other immune microenvironments.

Figure 6. Multi-gene set KEGG enrichment analysis and modular gene visualization. (A) The bubble diagram of the intersection of the KEGG pathway in the MEblue module, including the immune system and infectious diseases: viral and immune diseases gene sets. (B) The bubble diagram of the intersection of the KEGG pathway in the MEgreen module, including the immune system and infectious diseases: viral gene sets.

Animals and diets

All animal experiments conformed to the Guide for the Care and Use of Laboratory Animals (NIH Publications No. 8023, revised 1978). All procedures and euthanasia performed in this study were approved by the Ethics Committee for Animal Experiments of the Department of Laboratory Animal Science, Peking University Health Science Center (No. LA201456, 24 February 2014).

We randomly divided 32 8-week-old male Sprague-Dawley rats (300–360 g) into the following four groups (n = 8 in each group): the control group was undamaged and divided into the SD group (Sham_SD) and the KD group (Sham_KD), while the SCI rats received either the SD (SCI_SD group) or KD (SCI_KD group); 20 rats were used for RNA-seq, and 12 rats were used for histopathological examination of fixed specimens. All rats were reared at 21°C with a light/dark cycle of 12/12 hours, and food and water were provided freely. Rats in the SCI_KD group were fed a KD (70 g/kg body weight/day), and the ratio of fat to carbohydrate and protein was 4:1 (fat: 50.5 g/100 g KD; carbohydrate: 23.6 g/100; four hours after SCI, protein content: 5.5 g/100 g KD; dietary fibre: 16.1 g/100 g KD (Shenzhen Zeneca Inc., Shenzhen, China)). The rats in the Sham_SD group and SCI _SD group were fed normal carbohydrates. All groups had free-access water supply.

C7 spinal cord semi-contusion

A C7 hemi-contusion was performed in the two SCI groups as described previously for a different site of contusion [20]. Briefly, rats were anaesthetized with sodium pentobarbital after a unilateral (left or right side, chosen randomly) C7 laminectomy, and the C6-T1 vertebrae were rigidly fixed in a frame tilted at a 25.0° angle. We then used an Infinite Horizon impactor to deliver a set contusion force of 150 kdyne [21]. The C7 hemi-contusion is hereafter described as ‘SCI.’ Rats were fed for four weeks according to the above grouping and diet plan, and spinal cord tissue was sampled after four weeks.

EC myelin staining

Spinal cord tissue was placed in a fully automatic tissue dehydrator for an overnight dehydration programme (xylene, alcohol isocratic dehydration), embedded in paraffin, and cut it into 5-μm continuous sections with a microtome after cooling overnight at 40°C. The tissues were soaked in xylene I and II for 30 minutes each, after which they were soaked in 100%, 100%, 95%, 95%, and 80% alcohol for 5 minutes each and then rinsed twice with double-distilled water. The tissues were then stained with eriochrome cyanine (EC) for 15–20 minutes before undergoing ferric chloride (5.6%) differentiation for a specific time under a microscope. We made sure not to over-differentiate, as the staining would be too light. Afterwards, the samples were dehydrated and mounted with neutral gum.

RNA isolation, RNA-seq library construction and sequencing

We randomly divided 20 Sprague-Dawley rats into 4 groups: SCI_KD; SCI_SD; Sham_KD; and Sham_SD. Total RNA was extracted from a 1.6-cm spinal cord segment at the centre of the impact site with the TRIzol reagent (Invitrogen, Thermo Fisher Scientific Inc., MA, USA) and subjected to RNA-seq analysis to assess gene expression (n = 5). Genomic DNA was removed using DNase I. Agarose gels (1%) were used to monitor RNA degradation and contamination. A NanoPhotometer^® spectrophotometer (IMPLEN, CA, USA) was used to check RNA purity. A Qubit^® RNA Assay Kit was used to measure the RNA concentration on a Qubit^® 2.0 Fluorometer (Life Technologies, CA, USA). An RNA Nano 6000 Assay Kit was used to assess the RNA integrity on a Bioanalyser 2100 system (Agilent Technologies, CA, USA). RNA sequencing of 5 biological replicate samples was performed for each group, and 5 samples were degraded and eliminated. Therefore, 15 RNA-seq cDNA libraries were established.

1) We utilized two methods to treat the total RNA. Oligo(dT) magnetic beads were used to select mRNAs with polyA tails. rRNAs were hybridized with DNA probes and then digested to generate a DNA/RNA hybrid strand, followed by DNase I digestion to remove DNA probes. The target RNAs were then obtained by purification. 2) The target RNA was fragmented and then reverse transcribed into double-stranded cDNA (dscDNA) with N6 random primers. 3) The dscDNA was end-repaired with phosphates at the 5′ end and sticky ′A′s at the 3' end, after which the sequence was ligated to an adaptor with sticky 'T's at the 3' end of the dscDNA. 4) Two specific primers were used to amplify the ligation product. 5) The PCR products were denatured with heat, and the single strand DNA was cyclized with a splint oligo and DNA ligase. 6) Finally, the prepared library was sequenced with the Illumina HiSeq xten (2 × 150-bp read length, Illumina, San Diego, CA, USA). The sequencing programmes were performed by BGI Company (China, Shenzhen).

Quality control and expression level

The raw paired-end reads were trimmed and subjected to quality control by SeqPrep (https://github.com/jstjohn/SeqPrep) and Sickle (https://github.com/najoshi/sickle) with default parameters. The paired-end clean reads were aligned to the Rattus_norvegicus reference genome (Rnor_6.0) using the default parameters in TopHat (http://tophat.cbcb.umd.edu/, version 2.1.1) [47]. The reference genome and gene model annotation files were directly downloaded from genome website (http://www.ensembl.org/Rattus_norvegicus/Info/Index). RSEM (http://deweylab.biostat.wisc.Edu/rsem/) [48] was used to quantify the gene abundances. To identify DEGs between two different samples, the expression levels of each transcript were calculated based on the number of transcript reads per million reads (TPM). The reads were then normalized using DESeq2 version 1.24.0 (http://bioconductor.org/packages/stats/bioc/DESeq2/) [49], and the statistical significance of the DEGs was assessed. The following default parameters and criteria were used here: Benjamini & Hochberg (BH) p-adjust < 0.05 and |log2FC|≥1 after multiple comparisons. In the initial data exploration, we performed principal component analysis (PCA) to directly calculate the coefficient of variation between the groups in R ‘prcomp’.

Gene function annotation analysis

The genes and transcripts are aligned to six databases (NR, Swiss-Prot, Pfam, EggNOG, GO, and KEGG) to obtain comprehensive annotation information. DIAMOND software (https://github.com/bbuchfink/diamond) [50] was used to sequence the genes and transcripts with the NR (ftp://ftp.ncbi.nlm.nih.gov/blast/db/), Swiss-Prot (ftp://ftp.uniprot.org/pub/databases), and EggNOG (http://eggnogdb.embl.de/#/app/home) databases [51]. Sequence alignment was performed with the BLAST2GO and GO database, followed by alignment with the Pfam database using HMMER software [52]. GO (http://www.geneontology.org/) is a database established by the Gene Ontology Federation. Its purpose is to standardize the biological terminology of genes and gene products in different databases and to determine gene and protein functions. Using the GO database, the genes in the selected gene set can be classified according to the biological process (BP) they participate in, the cell component (CC), and the molecular function (MF) they achieve. The results of KEGG Orthology (Kyoto Encyclopaedia of Genes and Genomes, http://www.genome.jp/kegg/) were obtained using KOBAS2.1. KEGG is a knowledge base that systematically analyses gene functions and links genomic information and functional information. Using the KEGG database, the genes in the gene set can be classified according to the pathways involved or the functions they perform [53].

KEGG pathway enrichment analysis

The KEGG PATHWAY enrichment analysis was performed on the gene set using an R script in KOBAS (http://kobas.cbi.pku.edu.cn/home.do) [53]. The BH method was used to calculate the corrected p value. The threshold of the corrected p value (corrected P-value) was 0.05. The KEGG pathway that met the conditions showed that the pathway was significantly enriched for the gene set.

iPath metabolic pathway analysis

iPath3.0 (http://pathways.embl.de) was used to visually analyse the metabolic pathways and facilitate observation of the metabolic pathway information of DEGs in the entire biological system. Nodes represent various biochemical molecules, and lines represent biochemical reactions. The overall metabolism of the biological system includes 146 metabolism-related KEGG pathways, 22 regulatory-related KEGG pathways, and 58 KEGG pathways involved in the biosynthesis of secondary metabolites. The nodes in the figure represent different compounds, and the boundaries represent different enzymatic reactions. iPath3.0 can outline the biosynthesis and important regulatory pathways of secondary metabolites and easily identify complex metabolic pathways [54].

Gene set enrichment analysis (GSEA)

GSEA is performed on a predefined set of genes, usually from functional annotations or the results of previous experiments (included in MSigDB, version 6.2, http://software.broadinstitute.org/gsea/downloads.jsp), and organizes the genes according to two sample types. The GSEA sorts the degree of differential expression and then checks whether the gene set clusters at the top or bottom of the ranking list. The first round of enrichment is based on C2, which includes curated gene sets (from online pathway databases, publications in PubMed, and knowledge of domain experts). The second enrichment is based on the results of KEGG functional annotation in this study. For the analysis results, we generally considered that gene sets with the following characteristics were meaningful: |NES|>1, NOM p-val<0.05, and FDR q-val<0.25.

Protein-protein interaction (PPI) analysis

We used an online search tool (STRING database, version 10.5; http://string-db.org/) to search for interacting genes to construct a protein-protein interaction (PPI) network based on the extensive functions of the genome-wide data [55].

Weighted gene co-expression network analysis (WGCNA)

A co-expression network was constructed using the WGCNA (version 1.47) package in R [21]. WGCNA is an undirected, weighted network. "Unweighted network" means that the correlation between genes can only be 0 or 1. 0 means that the two genes are not related, while 1 means they are related. A "weighted network" means not only that genes are related or not but also that their correlation value is recorded. The value is the weight (correlation) of the connection between genes. The soft threshold (β power) in the co-expression network is the expression of each gene at a specific time or space and is regarded as a node. To obtain the correlation between genes, it is necessary to calculate the correlation coefficient (Pearson coefficient) any two genes; this allows us to know whether the expression profiles of two genes are similar. In this study, the β power = 6.

1) The module refers to the grouping of genes with similar expression patterns. Genes of a specific type below to a module.

2) An eigengene is defined as the first principal component of a given module. It is the most representative expression pattern of the module. It is an indicator that approximately represents the expression status of all genes in a module but is not a real gene.

3) A hub module refers to a module that is highly related to phenotypic data or closely related to research.

4) Gene significance (GS), or a measure of gene importance, is the correlation between gene expression patterns and the characteristics of a sample. Generally, the larger the correlation value, the more important the gene, which can be evaluated as the log conversion of the p-value. The correlation index between a gene and a specific pathway can also be the correlation between a gene and a specific phenotypic feature.

5) Module membership (MM), also known as the eigengene-based connectivity (kME), is the correlation between the expression pattern of a gene and a characteristic gene in the module. The larger the value, the greater the possibility that the gene belongs to the module.

Genes with the highest MM and the highest GS are genes with high significance (hub genes). All these analyses were performed using commands implemented in the WGCNA software package with the following parameters: minModuleSize = 30, minKMEtoStay = 0.3, and mergeCutHeight = 0.2.

Modular gene expression correlation analysis

After a gene of interest is obtained based on the correlation of gene expression, the correlation coefficient between genes is determined using the Spearman correlation algorithm. The correlation coefficient threshold is 0.5, the significance level is the p value after BH correction (corrected p-value) ≤0.05, and these values were drawn into a visual network diagram. In the figure, the nodes represent genes, and the connections between nodes represent the correlation between gene expression. The larger the node is, the greater the expression correlation between the genes.

Multi-gene set KEGG pathway enrichment analysis

After performing KEGG analysis on the different gene sets, we utilized the intersection of the selected gene sets for the display and adjusted the significance level to the corrected p-value ≤0.5 as the threshold to display the top 20 pathways.

Statistical analysis

The data are expressed as the mean ± standard deviation. The statistical analysis was performed with GraphPad Prism 7.0 (GraphPad Software Inc., San Diego, CA). Student's t-test was used for comparisons between two groups. A p value < 0.05 was considered statistically significant. ^*p < 0.05, ^**p < 0.01, ^***p < 0.001, and ^****p < 0.0001.

Availability of data and materials

RNA-seq raw data have been deposited in the NCBI Sequence Read Archive (SRA, https://submit.ncbi.nlm.nih.gov/subs/sra/SUB7836712/overview). Accession IDs for Rattus norvegicus BioProject = PRJNA648739; BioSample  = SAMN15646629 – SAMN15646643 (15 objects); SRA = SRR12327118 – SRR12327132 (15 objects).