1. Introduction
Venous thromboembolism (VTE), including deep venous thrombosis (DVT) and pulmonary embolism (PE), is a common and severe complication of infectious disease. This complication should always be considered in patients who present with a musculoskeletal infection, especially osteomyelitis (OM) in younger ones [Reference Funk and Copley1, Reference Altobelli and Quinonez2]. In recent years, the prevalence of DVT among acute haematogenous osteomyelitis (AHO) cases has been reported to be 6–10% [Reference Ligon, Journeycake, Josephs, Tareen, Lindsay and Copley3, Reference Bouchoucha, Benghachame and Trifa4]. In Taiwan, the risk of developing DVT was 2.49-fold in patients with chronic OM compared with the comparison cohort [Reference Lin, Chen and Huang5].
Previous studies have evaluated the clinical characteristics of OM among children with and without DVT [Reference Ligon, Journeycake, Josephs, Tareen, Lindsay and Copley3, Reference Bouchoucha, Benghachame and Trifa4]. Staphylococcus aureus (S. aureus) is the predominant causative agent for OM [Reference McNeil6], and clinical outcomes were worse in patients caused by methicillin-resistant S. aureus (MRSA) than in those affected by methicillin-susceptible S. aureus (MSSA) [Reference Davis and Gilbert7]. In patients with AHO, Virchow’s triad of hypercoagulability, venous stasis, and injury to the vessel wall also applies and is thought to trigger thrombosis. Animal studies have revealed the pathophysiological roles of various influencing factors, including leukocytes, platelets, and neutrophil extracellular traps (NETs), on thrombosis [Reference Wolberg, Rosendaal and Weitz8]. Besides, there is increasing evidence for an association of Panton–Valentine leukocidin (PVL)-expressing S. aureus strains with AHO severity [Reference Niemann, Bertling and Brodde9]. The possible mechanistic links between DVT and OM may be platelets activated by PVL-damaged neutrophils via neutrophil secretion products [Reference Niemann, Bertling and Brodde9]. However, the mechanisms by which the transcriptomic response of OM contributes to thrombosis development remain incompletely understood, although there is increasing evidence suggesting an extensive “cross-talk” between the inflammatory and coagulation cascades [Reference Branchford and Carpenter10], and these signatures may be useful in the diagnosis of venous thrombosis.
Transcriptional profiles are used increasingly to investigate the severity, subtype, and pathogenesis of a disease, which have implications for diagnosis and therapeutic development [Reference Dunning, Blankley and Hoang11–Reference Berry, Graham and McNab13]. Therefore, blood transcriptional profiles and host response signatures could serve as biomarkers of clinical changes in subjects at risk for or diagnosed with venous thrombosis in OM. To further explore the underlying pathophysiology, we examined the gene expression alterations involved in OM and venous thrombosis. In this study, we downloaded the gene expression profiles for OM and VTE from the Gene Expression Omnibus (GEO) database. We performed a gene expression analysis using the GEO2R web tool to identify the differentially expressed genes (DEGs) of OM and VTE with their respective controls and subsequently developed Gene Ontology and pathway enrichment analysis for screening of DEGs with the g:Profiler toolset. Finally, an integration of the DEGs protein-protein interaction (PPI) network was constructed and the module analysis was performed. The identified hub genes might play important roles in the process of venous thrombosis which developed from OM. Overall, our findings will hopefully deepen the understanding of the relationships between OM and VTE. A flowchart summarizing this study is shown in Figure 1.
FIGURE 1: A flowchart of the present study.
2. Materials and Methods
2.1. Data Acquisition and Processing
Gene expression datasets related to OM and VTE were obtained from the GEO database of NCBI (https://www.ncbi.nlm.nih.gov/geo/) [Reference Barrett, Wilhite and Ledoux14]. The keywords “osteomyelitis” and “venous thromboembolism” with “Homo sapiens” or “human” were employed to mine the dataset. Finally, 3 datasets were selected from published studies and downloaded from the GEO database, including GSE30119 for OM and GSE19151 and GSE48000 for VTE. The related information of these datasets is given in Table 1. These blood samples were further divided into different groups according to the diseases source. DEGs were screened using transcription profile data of whole blood samples. All data were normalized and log (base 2) transformed.
TABLE 1: Published datasets related to OM and VTE used in this study.
2.2. Identification of DEGs
According to the recurrence or degree of risk, VTE patients in the two datasets GSE19151 and GSE48000 were grouped and compared with the healthy controls, respectively. The patients with S. aureus OM in the dataset GSE30119 were compared with healthy controls. GEO2R web tool (https://www.ncbi.nlm.nih.gov/geo/geo2r/) [Reference Davis and Meltzer18] using the GEOquery and limma R packages from the Bioconductor project was utilized to compare the differences in gene expression between two or more groups of samples in a GEO dataset (Table 1) and identify DEGs associated with the diverse experimental conditions. The adjusted P values were used to decrease the false positive rate using the Benjamini and Hochberg false discovery rate method by default. The threshold value for identifying DEGs was set as adjusted P value < 0.05 and |log2(fold change)| ≥ 1. The results of the pairwise comparisons were summarized for subsequent analysis.
2.3. Venn Diagram Analysis of DEGs
We used an online integrative tool Venny (http://bioinfogp.cnb.csic.es/tools/venny/index.html) to analyze the similarities and differences of DEGs in the three datasets, GSE30119, GSE19151, and GSE48000. The identified DEGs of the two datasets associated with VTE were merged, whereas DEGs which were observed in opposite expression directions were discarded. Afterward, two specific clusters of DEGs were defined as OM-related DEGs (OMGs1, obtained from GSE30119) and VTE-related DEGs (VTEGs1, obtained from GSE19151 and GSE48000), respectively. Additionally, the intersections of OMGs1 and VTEGs1 were defined as shared DGEs1.
2.4. Functional and Pathway Enrichment Analyses of DEGs
To discern the implication of specific clusters and shared DEGs on OM and VTE, we used the g:Profiler (https://biit.cs.ut.ee/gprofiler/) [Reference Raudvere, Kolberg and Kuzmin19] toolset to perform Gene Oncology (GO) terms and pathway enrichment analysis. The shared biological process ontology (GO-BP) was identified based on overlapping GO term IDs between functional enrichment analysis results of OMGs1 and VTEGs1. The DEGs enclosed by all GO-BP terms shared between OMGs1 and VTEGs1 were defined as DEGs2 (VTEGs2 and OMGs2, respectively) and used to predict the significant pathways and ensuing crosstalk between these pathways. Additionally, the intersections of OMGs2 and VTEGs2, suggested the potential mechanism of patients with OM in whom VTE developed, were defined as shared DGEs2. The Benjamini and Hochberg false discovery rate method was used to correct the P value. An adjusted P value < 0.05 was considered to have statistical significance and to achieve significant enrichment.
2.5. Identification of Transcription Factors (TFs) That Are Significantly Associated with the DEGs
To link gene expression signatures to upstream cell signaling networks, we used Transcriptional Regulatory Relationships Unraveled by Sentence-based Text mining (TRRUST) (https://www.grnpedia.org/trrust/) [Reference Han, Cho and Lee20] to identify the upstream TFs that regulate VTEGs2 and OMGs2. The shared TF-target gene pairs of VTEGs2 and OMGs2 may be involved in a common pathway.
2.6. PPI Network and Module Analysis
The Search Tool for the Retrieval of Interacting Genes Database (STRING) (https://www.string-db.org/) [Reference Szklarczyk, Franceschini and Wyder21] was used to expose the protein-protein interaction (PPI) information among VTEGs2 and OMGs2 at the protein level. Afterward, the network was recreated and visualized using Cytoscape (http://www.cytoscape.org/) [Reference Shannon, Markiel and Ozier22], based on the coexpression graph of the PPI network. Topological properties (such as degree of distribution, betweenness centrality, and closeness centrality) of the constructed PPI network were calculated by NetworkAnalyzer in Cytoscape. Nodes with a higher degree based on the number of edges (interactions) between various nodes (proteins) were regarded as hub proteins. The Molecular Complex Detection (MCODE) [Reference Bader and Hogue23] plugin of Cytoscape was utilized to perform module analysis and choose hub modules of the PPI network with a degree cutoff = 2, node score cutoff = 0.2, k-core ≥ 2, and max. depth = 100.
3. Results
3.1. Identification of DEGs
Comparisons between differentially gene expression between two or more groups of samples from 1 OM dataset (GSE30119) and 2 VTE datasets (GSE19151 and GSE48000) were performed using the GEO2R web tool, respectively, and mRNAs with an adjusted P value < 0.05 and |log2(fold change)| ≥ 1 were identified as DEGs. Two specific clusters of DEGs are shown in Figure 2. A total of 587 OMGs1 were identified in GSE30119. A total of 382 VTEGs1 were identified in VTE-related datasets, including 186 in GSE19151 and 206 in GSE48000. Furthermore, 27 shared DGEs1 that overlap OMGs1 and VTEGs1 were screened out, with 2 genes (LSP1 and PADI4) found to be observed in opposite expression directions. After applying functional annotation with the Gene Ontology database to find terms associated with the development of VTE in OM patients, 230 VTEGs2 and 333 OMGs2 sharing the common biological process ontology were detected and used for further analysis (Figures 3(d) and 4(a)). Afterward, Venn diagram analysis of VTEGs2 and OMGs2 was performed and 13 shared DGEs2 were detected (Figures 4(a)-4(b)). There were 11 upregulated DEGs and 1 downregulated DEG in shared DGEs2, with 1 gene (LSP1) found to be observed in opposite expression directions (Figure 4(c)).
FIGURE 2: A Venn plot of OMGs1 and VTEGs1. Heart-shaped dotted line represents the shared DGEs1.
FIGURE 3: Results of the GO-BP terms enrichment analyses for (a) OMGs1, (b) VTEGs1, and (c) shared DGEs1 and (d) shared GO-BP enrichment of VTEGs1 and OMGs.
FIGURE 4: The DEGs2 enclosed by all sharing GO-BP terms. (a) A Venn plot of OMGs2 and VTEGs2. (b) PPI network of the shared DGEs2. (c) Fold change of the shared DGEs2.
3.2. Functional Annotation of DEGs
To gain more biological insight, we performed GO enrichment analysis using the g:Profiler toolset. As shown in Figure 3, the activity of neutrophils has attracted widespread attention among biological processes. With respect to OMGs1, the following GO-BP terms were significantly enriched: neutrophil activation, neutrophil degranulation, neutrophil-mediated immune response, and defense response to the bacterium (Figure 3(a)). Highly similar results of enrichment analysis were also presented in shared DGEs1 (Figure 3(c)). There suggested a significant correlation with the biological process of neutrophils and the process of VTE which developed in OM patients. In addition, VTEGs1 were mainly enriched in the viral gene expression, nuclear-transcribed mRNA catabolic process, protein targeting to membrane, and protein localization to the endoplasmic reticulum (Figure 3(b)). The shared GO-BP terms of OMGs1 and VTEGs1 seemed to be focused on cell activation, immune defense, and inflammatory response to stress or biotic stimuli (Figure 3(d)). These results showed that DEGs were mainly enriched in host immune defense processes such as neutrophil activation and degranulation to resist bacterial invasion or other stress.
3.3. Identification of TF-Target Gene Pairs
In the TRRUST database, 36 and 43 upstream TFs were detected to regulate VTEGs2 and OMGs2, respectively. A total of 13 shared upstream TFs were E2F1, GATA1, JUN, NFKB1, RELA, RUNX1, SIRT1, SP1, SPI1, STAT1, STAT3, TP53, and YY1, and these data are given in Table 2. Furthermore, 3 TF-target gene pairs may share a common pathway, consisting of SP1-LSP1, SPI1-FCGR1A, and STAT1-FCGR1A. Interestingly, only one TF-target gene pair has a P value less than 0.05.
TABLE 2: The list of TF-target gene pairs of VTEGs2 and OMGs2.
3.4. Hub Proteins and Module Screening of the PPI Network
A PPI network was utilized to visualize protein-protein interactions involved in the development of VTE in OM patients, and then, potential hub genes between VTEGs2 and OMGs2 were detected. Nodes that did not interact with other nodes were excluded. There were a total of 1611 interactions among 467 nodes (Figure 5(a)). The genes with the highest degree (>30) of interactions within the network were TP53, IL4, MPO, ELANE, FOS, and CD86 (Table 3). To screen significant modules in the PPI network, the MCODE plugin was used in Cytoscape, and the module analysis results are given in Table 4. Top 3 of the 16 modules were illustrated (Figures 5(b), 5(c), 5(d)). The top 5 proteins with higher node degrees in module 1 were ELANE (degree = 32), HP (degree = 29), ORM2 (degree = 25), ARG1 (degree = 25), and CAMP (degree = 24) (Figure 5(b)). The top 5 proteins in module 2 were SOCS3 (degree = 27), RBX1 (degree = 25), UBE2D1 (degree = 24), RPS15A (degree = 22), and RPS17 (degree = 21) (Figure 5(c)). Additionally, the top 5 proteins in module 3 were SNRPG (degree = 26), NDUFA4 (degree = 20), SNRPD2 (degree = 19), COX7C (degree = 19), and HINT1 (degree = 19) (Figure 5(d)). Interestingly, we found that 4 shared DGEs2 were included in these 16 modules, consisting of SLPI (module 1, degree = 17), TCN1 (module 1, degree = 14), CREG1 (module 4, degree = 11), and FCGR1A (module 16, degree = 19) (Table 4).
FIGURE 5: Analysis of the PPI network of DEGs2. (a) PPI network of OMGs2 and VTEGs2. (b) Module 1 of the PPI network of DEGs2. (c) Module 2 of the PPI network of DEGs2. (d) Module 3 of the PPI network of DEGs2. The node sizes correspond to the degree of the node, while the node color denotes betweenness centrality.
TABLE 3: The top 30 hub proteins identified in the topological analysis of the PPI network of DEGs2.
TABLE 4: The module analysis results with MCODE of the PPI network of DEGs2.
4. Discussion
As a fatal disease caused by serious musculoskeletal infection, VTE needs to arouse the attention of clinical staff. It has been reported that children with OM may have increased susceptibility to VTE in previous epidemiologic studies [Reference Ligon, Journeycake, Josephs, Tareen, Lindsay and Copley3–Reference Lin, Chen and Huang5]. Unfortunately, the transcriptional signatures between VTE and OM have not been intensively investigated. The main objective of this study was to focus on revealing the underlying pathophysiology association between VTE and OM.
In the present study, we integrated 2 expression profiles from VTE patients and 1 expression profile from OM patients to identify genes that may play a crucial role in the onset and development of VTE in OM patients. First, 382 VTEGs1 and 587 OMGs1 were identified in VTE and OM patients, respectively. Following these, DEGs1 were evaluated by functional enrichment analysis to get insight into the biological significance in the pathogenesis. It is reported that blood coagulation was overexpressed in patients with osteoarticular infections [Reference Banchereau, Jordan-Villegas and Ardura15, Reference Chen, Yao and Deng24]. Here, our research showed that OMGs1 and shared DGEs1 were mainly enriched in some neutrophils-related pathways, such as neutrophil activation, neutrophil degranulation, neutrophil-mediated immune response, and defense response to the bacterium. Moreover, the shared GO-BP terms of OMGs1 and VTEGs1 seemed to be focused on cell activation, immune defense, and inflammatory response to stress or biotic stimuli. These results indicated that the activity of neutrophils and the immune response to bacterial invasion may explain the molecular mechanisms of VTE developed in OM patients to some extent. After that, 230 VTEGs2 and 333 OMGs2 sharing the common biological process ontology were detected by comparing GO-BP terms of VTEGs1 and OMGs1. Additionally, 13 shared DGEs2 that overlap OMGs2 and VTEGs2 were identified for further analysis.
It is well known that TFs, specific to binding to its target gene, exert facilitative or inhibitory roles in gene expression, showing an important part in the multitude of biological processes involved in diseases [Reference Han, Cho and Lee20]. In this study, based on the Venn diagram analysis method, 13 shared upstream TFs of VTEGs2 and OMGs2 were screened out. It is suggested that 3 TF-target gene pairs, consisting of SP1-LSP1, SPI1-FCGR1A, and STAT1-FCGR1A, may reveal the potential pathogenesis link between VTE and OM. Finally, in the PPI network analysis, TP53, IL4, MPO, ELANE, FOS, CD86, HP, SOCS3, ICAM1, and SNRPG were the top 10 hub proteins with the highest connectivity within the network. Significant modules were further detected in the PPI network, and some specific nodes with the highest connectivity (i.e., ELANE, HP, SOCS3, and SNRPG) were also identified in the top 3 modules. Particularly, 4 shared DGEs2 (SLPI, TCN1, CREG1, and FCGR1A) were found in the 16 modules.
OM, especially haematogenous OM, can cause an excessive inflammatory response in blood vessels and coagulation disorders modulated by various microorganisms [Reference Lichota, Gwozdzinski and Szewczyk25, Reference Chen, Lin and Chen26]. S. aureus is the predominant causative pathogen of OM [Reference Mehkri, Felisma, Panther and Lucke-Wold27]. Surface proteins and exotoxins of S. aureus can trigger thrombosis through effects on the coagulation pathway and on anticoagulation factors. Additionally, S. aureus can activate endothelial surfaces and platelets. It has been reported that some exotoxins such as PVL expressed by S. aureus can contribute to leukocyte lysis and additional damage to endothelial surfaces [Reference Niemann, Bertling and Brodde9, Reference Martin, Cevik and Nugent28]. During S. aureus infection, platelets amplify the recruitment and activation of innate immune cells at the site of infection and help eliminate pathogens. In some cases, these mechanisms can lead to thromboinflammation, leading to severe organ dysfunction. These events can lead to microthrombosis and DVT. Neutrophils act as the first line of defense in the process of acute infection. They can protect the human body from various pathogens through antibacterial mechanisms such as phagocytosis, degranulation, and the production of reactive oxygen species (ROS). According to recent research, neutrophils can also exert their innate immune function by forming a special structure called NETs [Reference Abrams, Morton and Alhamdi29–Reference Eustache, Tohme, Milette, Rayes, Tsung and Spicer31]. Our previous study also suggested that NETs may play a crucial role in the immune response of patients with OM during the S. aureus infection process [Reference Chen, Yao and Deng24]. In this study, some core nodes (such as MPO, ELANE, and CAMP) were essential components of the NETs network structure. Emerging research focused on how NETs affect platelet function, particularly in the setting of infection and inflammation [Reference Elaskalani, Abdol Razak and Metharom32–Reference Yang, Wu and Long34]. In patients with sepsis, systemic inflammation primed neutrophils to release NETs, which promoted thrombin generation through platelet-dependent and platelet-independent mechanisms, with the increased risk of VTE [Reference Gould, Vu and Swystun35, Reference Yang, Qi and Kan36]. There may also be a similar underlying mechanism between S. aureus OM and the development of DVT in the extremities. In addition, early markers detection and targeted intervention may minimize this effect in the inflammatory process of thrombosis in patients with OM.
Nowadays, TFs are becoming attractive biomarkers in complex diseases including infectious diseases. Among the identified regulatory TFs, Sp1 participated in the hypoxia-induced upregulation of VWF [Reference Mojiri, Alavi and Lorenzana Carrillo37]. Increasing evidence indicates that the interaction between NETs and VWF contributes to arterial and venous thrombosis as well as inflammation [Reference Yang, Wu and Long34]. As a target gene of SP1, endothelial LSP1 is recruited into the cytoskeleton during inflammation and plays an important role in the formation of endothelial dome, thereby regulating the transendothelial migration of neutrophils [Reference Petri, Kaur and Long38]. The haematopoietic transcription factor SPI1/PU.1 may be specifically involved in the differentiation or activation of macrophages or B cells. Moreover, SPI1 was able to restrain neutrophil defense by broadly inhibiting the accessibility of enhancers via the recruitment of histone deacetylase 1 [Reference Fischer, Walter and Tonges39]. STAT1 can be phosphorylated by receptor-associated kinases and then shapes’ cell viability in response to different cell stimuli and pathogens. The activation of STAT1 inhibits angiogenesis [Reference Battle, Lynch and Frank40], which may be involved in the inflammatory process of thrombosis in infectious diseases. FCGR1A is the common downstream target gene of SPI1 and STAT1. FCGR1A, also known as CD64, as an important immune-related protein, is closely related to various inflammatory processes. In posttraumatic OM patients, immunocompetent cells, predominantly highly activated polymorphonuclear neutrophils, as characterized by low expression of CD62L and high expression of the adhesion protein CD18, of the high-affinity IgG receptor CD64 and of the LPS receptor CD14, were recruited into the site of infection [Reference Wagner, Kondella, Bernschneider, Heppert, Wentzensen and Hänsch41].
Additionally, other hub proteins may also play a crucial role in the pathogenesis of VTE in patients with OM. TP53, the top hub gene of our study, is a critical tumor suppressor and a key regulator in numerous cellular functions, which maintains critical functions in immunity, inflammation, and tissue repair [Reference Joruiz, Beck, Horikawa and Harris42]. Intriguingly, TP53 may participate in the host-virus interactions that could characterize shared biological mechanisms between acute respiratory distress syndrome (ARDS) and VTE in severe COVID-19 patients [Reference Mishra, Chanchal and Ashraf43]. A similar host-bacterial interaction may also be reflected between OM and VTE. IL4 is directly involved in the bone desorption and osteoclast activity regulation that occur in OM. The association between IL4 gene polymorphisms (-1098-G/T and -590-C/T) and chronic OM was previously observed [Reference Tsezou, Poultsides and Kostopoulou44]. However, IL4 was not significantly elevated in a mouse model of fracture fixation with S. aureus OM [Reference Rochford, Sabate Bresco and Zeiter45]. The possible links between inflammation-related genetic variants (including IL4) and VTE established the fundamental role of genetic background in predisposition to VTE and several inflammation-related conditions [Reference Saghazadeh and Rezaei46]. HP, one of the acute-phase reactants, can stimulate inflammatory bone loss and its phenotype 2-2 has been reported to be a risk factor for spontaneous VTE [Reference Vormittag, Vukovich and Mannhalter47]. SOCS3 was reported to be a key player in bone-associated inflammatory responses, which acted as the cytokine-inducible negative regulator of cytokine signaling [Reference Gao and Van Dyke48]. The FOS proteins have been implicated as regulators of cell proliferation, differentiation, and transformation. The c-Fos gene has been put forward as a new factor in the progression and severity of atherosclerosis [Reference Saliques, Teyssier and Vergely49]. CD86 plays an important role in immune responses as a costimulatory molecule on antigen presenting cells. The gene polymorphisms of CD86 may be related to the risk of sepsis with contradictory results in different studies [Reference Fu, Zhang, Bai, Chen, Liu and Bai50, Reference Zou, Cai, Li and Wu51]. Unfortunately, there is no relevant research report on FOS or CD86 in the pathogenesis of OM or VTE.
This study has several limitations. Although followed the same analysis strategy and threshold for DEGs screening, the number of differential genes shared by the GSE19151 and GSE48000 datasets obtained in the end was small due to the inconsistency of the inclusion criteria and sequencing platform. The DEGs from these two datasets were included in the subsequent analysis, which can provide us with more possibilities, but they were not differentially expressed at the same time. More external datasets were needed to verify the results of the study. In addition, further experimental investigations are warranted to decipher the roles of these hub genes in the development of VTE in OM patients.
5. Conclusions
In this study, we reported for the first time the identification of the DEGs and activated signaling pathways between OM and VTE by the integrated bioinformatic approach. Our results finally identified a new set of novel biomarkers and important molecular targets, including 3 TF-target gene pairs (SP1-LSP1, SPI1-FCGR1A, and STAT1-FCGR1A), 3 structural proteins of NETs (MPO, ELANE, and CAMP), and 6 other hub proteins (TP53, IL4, HP, SOCS3, FOS, and CD86), which might play essential biological roles during the progression of VTE in OM patients. In addition, the enriched neutrophils-related pathways (neutrophil activation, degranulation, and immune response to bacteria) could advance our understanding of the development of OM towards VTE and indicate new avenues to develop therapies for VTE. However, more research, especially experimental research, will be indispensable for the further clinical application of these biomarkers.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This study was supported by the Health Scientific Innovation Platform Program of Fuzhou (2019-S-wp3), Natural Science Foundation of Fujian Province (2020J011186 and 2021J011315), Science and Technology Planning Project of Fuzhou (2020-WS-87), Young and Middle-Aged Scientific Research Project of Fuzhou Municipal Health Commission (2021-S-wq12), Project of Fujian Provincial Clinical Medical Research Center for First Aid and Rehabilitation in Orthopaedic Trauma (2020Y2014), and Project of Fuzhou Trauma Medical Center (2018080303).
Open access
