Systemic bioinformatics analysis of skeletal muscle gene expression profiles of sepsis.

Sepsis is a type of systemic inflammatory response syndrome with high morbidity and mortality. Skeletal muscle dysfunction is one of the major complications of sepsis that may also influence the outcome of sepsis. The aim of the present study was to explore and identify potential mechanisms and therapeutic targets of sepsis. Systemic bioinformatics analysis of skeletal muscle gene expression profiles from the Gene Expression Omnibus was performed. Differentially expressed genes (DEGs) in samples from patients with sepsis and control samples were screened out using the limma package. Differential co-expression and coregulation (DCE and DCR, respectively) analysis was performed based on the Differential Co-expression Analysis package to identify differences in gene co-expression and coregulation patterns between the control and sepsis groups. Gene Ontology terms and Kyoto Encyclopedia of Genes and Genomes pathways of DEGs were identified using the Database for Annotation, Visualization and Integrated Discovery, and inflammatory, cancer and skeletal muscle development-associated biological processes and pathways were identified. DCE and DCR analysis revealed several potential therapeutic targets for sepsis, including genes and transcription factors. The results of the present study may provide a basis for the development of novel therapeutic targets and treatment methods for sepsis.

Introduction

Sepsis is a type of systemic inflammatory response syndrome with high morbidity and mortality due primarily to multiple organ failure (MOF), including lung injury, renal dysfunction and skeletal problems, particularly in elderly people (1–3). Despite recent advances, sepsis remains an economic and humanistic burden to society that is increasing in developing and developed countries due to delayed diagnoses (4–6). Identifying the biomarkers of MOF may increase our understanding of the underlying mechanisms and allow for earlier diagnosis of sepsis.

Bacterial culturing of blood samples, which requires several days to complete, is a traditional method used to diagnose sepsis (4). Polymerase chain reaction (PCR)-based molecular quantification has emerged as a novel diagnostic method for several diseases, such as sepsis (7). For example, 16S ribosomal RNA gene amplification by quantitative PCR and sequencing has been reported to improve the sensitivity, specificity, and positive and negative predictive value of bacteria detection in neonatal sepsis (7). Furthermore, the development of ‘omics’ technologies, including high-throughput sequencing and mass spectrometry, has accelerated the exploration of biomarkers for early diagnosis and treatment of sepsis (8). A study by Scherag et al (9) identified vacuolar protein sorting 13 homolog A, cysteine rich secretory protein LCCL domain containing 2 and chromosome 13 loci as potential hot points for sepsis via a genome-wide association study. A study by Gosiewski et al (10) demonstrated that bacterial DNA in the blood was a valuable biomarker for sepsis via high-throughput cationic trypsinogen sequencing. It has also been reported that elevated levels of serum trypsin and the PRSS1 mutation may also contribute to susceptibility of sepsis (11).

Skeletal muscle dysfunction is a major risk factor of sepsis, as well as a complication, and contributes to sepsis-associated mortality (12). The perturbation of several hormones, proteins and gene expression in skeletal muscle have been demonstrated to be associated with the outcome of sepsis, including glucocorticoids (13), lipopolysaccharide (14), regulated in development and DNA damage response 1 (15) and glucose transporter type 4 (16) expression. Gene co-expression and coregulation are important indicators of their functions in some biological processes or diseases (17). Altered co-expression or coregulation patterns under different conditions may indicate the role of different genes in specific diseases. To the best of our knowledge, there have been no previous studies investigating differential co-expression (DCE) or differential coregulation (DCR) analysis in the skeletal muscles of patients with sepsis. In the present study, DCE and DCR analysis of gene expression was performed using the skeletal muscles of patients with sepsis. Several genes and transcription factors (TFs) that may contribute to the progression of sepsis were identified. These results may be helpful to improve our understanding of the underlying mechanisms and potential effective treatments for sepsis.

Materials and methods

Microarray datasets

The skeletal muscle transcriptome-wide expression profiles of 13 patients with sepsis and 8 controls were downloaded from the Gene Expression Omnibus (ncbi.nlm.nih.gov/geo/; accession number, GSE13205) deposited by Fredriksson et al (18). No significant differences were observed in age (P=0.115; Student's t-test) or sex (P=0.336; Fisher's exact test) between the groups.

Gene expression profile analysis

The analysis of gene expression profiles consisted of two steps: Preprocessing and differential expression (DE) analysis. Briefly, the raw microarray datasets were imported into R version 3.4.3 (https://www.r-project.org/), a free open-source statistical software package, and the affy version 1.56.0 (http://bioconductor.org/packages/release/bioc/html/affy.html) (19) package was used for background correction and normalization in order to compare expression profiles of different samples. DE genes (DEGs) with the thresholds of fold change >2 and Benjamin adjusted P<0.05 were screened out using the limma version 3.34.8 (http://www.bioconductor.org/packages/release/bioc/html/limma.html) (20) package. Additionally supervised two-way clustering, including cluster analysis of samples and DEGs, was conducted through the Euclidean distance method (21) based on pheatmap 1.0.8 package (https://cran.r-project.org/web/packages/pheatmap/index.html).

Functional enrichment analysis

To explore the functions associated with DEGs, functional enrichment analysis was performed using the Database for Annotation, Visualization and Integrated Discovery (DAVID; david.ncifcrf.gov/) (22). Gene Ontology (GO) terms (http://www.geneontology.org/) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways (http://www.kegg.jp/), which had P<0.05 and a gene count ≥2, were identified in the present study.

DCE and DCR analysis

As a type of polygenic inheritance, the progression of sepsis is associated with multiple gene interactions; similarly, the co-expression of two genes indicates that they may serve similar functions in sepsis (17). Differences in co-expression patterns in the control and sepsis groups may reveal the roles of particular genes in a specific biological process or disease. This type of link is a DCE link (DCL). Genes involved in DCLs are considered to be DCE genes if the links surrounding them are significantly enriched in a binomial probability model (P<0.05); these genes are differential co-expression genes (DCGs).

TFs serve important roles in the regulation of gene expression. Multiple genes are regulated by common TFs, the differential expression of which may result in DCE among those genes. This is referred to as DCR (17). In the present study, DCE and DCR analyses were performed based on the DCGL version 2.1.2 (https://cran.r-project.org/web/packages/DCGL/index.html) (17) package. DCGs were screened out when the P-value of the differential co-expression enrichment (DCe) test and Q-value of the differential co-expression profile (DCp) test in DCGL were <0.05. Two types of DCL were identified: TF_bridged_DCLs, in which the two DCGs share ≥1 TF, and TF2target_DCLs, in which TF-gene regulation relationships were identified. The DCL network was visualized using Cytoscape 3.6.0 software (http://www.cytoscape.org/) with nodes and edges, representing genes and interactions, respectively.

Results

Differential expression genes

The preprocessing step produced comparable expression profiles among all the samples (Fig. 1A), which were subsequently used for DE analysis. A total of 1,052 DEGs were screened out in sepsis samples, including 441 downregulated and 611 upregulated genes (Fig. 1B). Supervised two-way hierarchical clustering of DEGs and samples was performed based on the pheatmap package (Fig. 1C).

Figure 1.

Analysis of gene expression microarray dataset. (A) Relative mRNA level following the preprocessing step. Red and blue boxes represent patients with sepsis and controls, respectively. (B) Scatter plot of DEGs. (C) Supervised clustering of samples and DEGs. DEGs, differentially expressed genes.

Enriched functions of DEGs

A total of 129 significantly enriched GO terms were identified using DAVID according to the thresholds of P<0.05 and gene count ≥2. The 10 most enriched significant GO terms are presented in Table I. The majority of GO terms were associated with biological processes in the skeletal muscle, RNA binding and cell responses to stimuli. A total of 12 KEGG pathways were identified and these were primarily associated with substance metabolism and synthesis, and inflammatory and cancer-associated processes (Fig. 2). The number of upregulated genes was markedly greater than the number of downregulated genes in the four most enriched pathways, suggesting that upregulated genes serve an important role in these pathways.

Table I.

Top 10 most significantly enriched GO terms of differentially expressed genes.

Category	GO name	Hits	P-value
CC	Z disc	19	<0.001
BP	Skeletal muscle tissue development	12	<0.001
CC	Cytoplasm	318	<0.001
CC	T-tubule	10	<0.001
BP	Muscle contraction	17	<0.001
BP	Cellular response to cadmium ion	7	<0.001
BP	Skeletal muscle contraction	8	<0.001
CC	Cytosol	212	<0.001
CC	Nucleolus	69	<0.001
BP	Response to denervation involved in regulation of muscle adaptation	5	<0.001

GO, gene ontology; CC, cellular component; BP, biological process.

Figure 2.

Significantly enriched Kyoto Encyclopedia of Genes and Genomes pathways of differentially expressed genes. Blue and red bars represent downregulated and upregulated genes contained in the corresponding pathways, respectively, and the green line is -log10 (P-value). PKG, protein kinase G; AMPK, AMP-activated protein kinase; HTLV-1, human T-cell lymphotropic virus type 1; HIF-1, hypoxia-inducible factor 1.

DCLs and DCRs

A total of four DCGs were identified using the DCGL package, including GTP binding protein 1 (GTPBP1), chaperonin containing T-complex 1 complex (CCT2), meiosis regulator of oocyte development (MIOS) and cytochrome P450 family 2 subfamily J member 2 (CYP2J2), and 93 DCLs. The network composed by the four DCGs and 90 non-DCGs was visualized using Cytoscape software (Fig. 3A) (23). No TF2target_DCLs were identified; however, 73 TF_bridged_DCLs were revealed using DCR analysis. The DCR network is presented in Fig. 3B. A total of eight TFs, including TATA box binding protein, melanocyte inhibiting factor-1 (MIF-1), signal transducer and activator of transcription 5B (STAT5B), POZ-, AT hook-, and zinc finger-containing protein 1, LIM homeobox 3 (LHX3)a, LHX3b, AP-2γ and AP-2αA, were considered to be associated with the development of sepsis due to their target enrichment in the DCGs.

Figure 3.

Differential co-expression and coregulation analysis. (A) Differential co-expression network. Circles and triangles represent DCGs and non-DCGs, respectively. (B) Differential coregulation network composed of TF_bridged_DCLs. Squares represent TFs, circles represent non-TFs, the pink circle represents DCGs, the blue circle represents non-DCGs and the gray square indicates TFs that were not tested in expression microarray data. Solid lines represent DCLs and dashed lines represent non-DCLs. Arrows indicate TF-to-target associations. DCG, differentially co-expressed gene; TF, transcription factor; DCL, differential co-expression link.

Discussion

Severe sepsis induces skeletal muscle dysfunction (24,25) and, in turn, skeletal muscle dysfunction is associated with the outcome of sepsis (12). Identifying biomarkers associated with this process may increase our understanding of the underlying mechanisms and allow for the development of therapeutic targets and novel treatment methods. In the present study, DE, DCE and DCR analysis were performed for this purpose. Functional enrichment analysis identified several inflammatory and skeletal muscle development-associated GO terms and KEGG pathways, which was consistent with previous studies (26,27).

Several skeletal muscle development-associated biological processes were identified in the present study, including skeletal muscle tissue development, muscle contraction and skeletal muscle atrophy. Several KEGG pathways associated with substance biosynthesis and metabolism were identified as significantly enriched in skeletal muscle samples from patients with sepsis compared with control samples. The majority of genes in those pathways were upregulated in patients with sepsis and all of the genes associated with ribosome biogenesis in eukaryotes were upregulated. Ribosome biogenesis dysfunction is associated with skeletal muscle inflammation and results in a decrease in muscle mass, which accounts for ~60% of sepsis-associated mortality (28,29). It may be a previously unidentified pathway that is associated with the dysfunction of skeletal muscle in sepsis. Several cancer-associated pathways were also identified in the present study, which may be because cancer and sepsis are inflammation-inducing diseases.

The regulation of gene expression serves a key role in several diseases (17). In the present study, 73 TF_bridged_DCLs, which contained four DCGs (CCT2, CYP2J2, GTPBP1 and MIOS) were identified. A total of 29, 17, 28 and 19 co-expression genes were identified for CCT2, CYP2J2, GTPBP1 and MIOS, respectively, in the DCE network. The direct association between these genes and sepsis progression or skeletal muscle dysfunction was not investigated in the present study. Previous studies have focused on their association with inflammation and cancer; for instance, CYP2J2 has been reported to regulate metabolic dysfunction via peroxisome proliferator-activated receptor-γ by reducing hepatic inflammation (30) and to protect against lung ischemia/reperfusion injury based on its anti-inflammatory effects (31). CCT2 and TCP1 have also been demonstrated to be valuable potential biomarkers for the prognosis of breast cancer (32). Such studies support the diverse functions of these genes and their role as potential therapeutic targets for sepsis or skeletal muscle dysfunction. The differential regulated rank (DRrank) method, which is based on DCGs and their co-expressed genes surrounding a specific TF, identified AP-2α and AP-2γ as the top two highest ranked TFs significantly associated with the development of sepsis, and they are considered to be associated with the progression of inflammatory diseases, such as sepsis (33–35). DNA binding of STAT5B, another potential target of sepsis, was also identified to serve a role in sepsis in a study by Chen et al (36), along with some novel potential sepsis-associated TFs, including MIF-1, ETS domain-containing protein Elk-1.

In conclusion, the present study identified aspects of the mechanisms of sepsis progression and skeletal muscle dysfunction using systemic bioinformatics analysis. DCE and DCR analysis provide an insight into gene regulation loops in sepsis and may be utilized to develop therapeutic targets and treatment methods for patients with sepsis.