Identification of the TF-miRNA-mRNA co-regulatory networks involved in sepsis

Sepsis is a life-threatening medical condition caused by a dysregulated host response to infection. Recent studies have found that the expression of miRNAs is associated with the pathogenesis of sepsis and septic shock. Our study aimed to reveal which miRNAs may be involved in the dysregulated immune response in sepsis and how these miRNAs interact with transcription factors (TFs) using a computational approach with in vitro validation studies. To determine the network of TFs, miRNAs, and target genes involved in sepsis, GEO datasets GSE94717 and GSE131761 were used to identify differentially expressed miRNAs and DEGs. TargetScan and miRWalk databases were used to predict biological targets that overlap with the identified DEGs of differentially expressed miRNAs. The TransmiR database was used to predict the differential miRNA TFs that overlap with the identified DEGs. The TF-miRNA-mRNA network was constructed and visualized. Finally, qRT-PCR was used to verify the expression of TFs and miRNA in HUVECs. Between the healthy and sepsis groups, there were 146 upregulated and 98 downregulated DEGs in the GSE131761 dataset, and there were 1 upregulated and 183 downregulated DEMs in the GSE94717 dataset. A regulatory network of the TF-miRna target genes was established. According to the experimental results, RUNX3 was found to be downregulated while MAPK14 was upregulated, which corroborates the result of the computational expression analysis. In a HUVECs model, miR-19b-1-5p and miR-5009-5p were found to be significantly downregulated. Other TFs and miRNAs did not correlate with our bioinformatics expression analysis. We constructed a TF-miRNA-target gene regulatory network and identified potential treatment targets RUNX3, MAPK14, miR-19b-1-5p, and miR-5009-5p. This information provides an initial basis for understanding the complex sepsis regulatory mechanisms.


Introduction
Sepsis is defined as a life-threatening organ dysfunction caused by a dysregulated host immune response to an infection (Singer et al. 2016). It is estimated that approximately 30 million people worldwide are diagnosed with sepsis every year (Walters 2018). The leading causes of sepsis in intensive care units are pneumonia, abdominal cavity infection, bloodstream infection, and urinary tract infections (Vincent et al. 2009). Although the research into sepsis and thus our understanding of the condition is increasing, the underlying pathological mechanisms of sepsis still need further study. The capacity to treat sepsis by supporting organ function and hemodynamics in intensive care units (ICUs) has increasingly strengthened, but sepsisrelated organ damage and its associated mortality remain high (Angus 2010). Therefore, because there is a paucity in understanding the precise molecular mechanisms underlying sepsis, it is important to find biological markers for early diagnosis and novel treatment targets of sepsis.
MicroRNAs (miRNAs) are short, noncoding RNAs of 18-25 nucleotides in length that regulate the translation of mRNAs (Doench and Sharp 2004). Mature miRNAs can bind to the complementary site of the 3′ untranslated region (UTR) of their specific target mRNAs and regulate translation at the post-transcriptional level (Ebert and Sharp 2012). Recent studies have found that the expression of miRNAs is associated with the pathogenesis of sepsis and septic shock (Benz et al. 2016;Reithmair et al. 2017). Transcription factors are DNA-binding proteins that can play significant roles in regulating gene expression and sepsis-induced organ dysfunction (Kumar et al. 2005). Zhang et al. have found that MYC and STAT3 may be the key regulatory genes in the underlying pathology of sepsis-induced ARDS (Zhang et al. 2019). In a previous study conducted by Mussbacher et al., the transcription factor NF-κB is critically involved in these pathophysiological processes as it induces both inflammatory and thrombotic responses (Mussbacher et al. 2019).
With the development and application of gene chip technology in the past decades, bioinformatics has received increasing attention as a research tool. Bioinformatics analysis can reveal molecular regulatory mechanisms in the progression of sepsis. In this study, we aimed to identify and construct the comprehensive miRNA-TF-mRNA coregulatory network in sepsis. Firstly, the potential targets of sepsis-related TFs and miRNAs were identified through the datasets from GEO. Secondly, further prediction of miRNA based on transcription factors that target differentially expressed genes (DEGs) and prediction of miRNAtranscription factor interactions using online tools. Based on these results, a TF-miRNA-mRNA network was then constructed, which can help us to reveal the complicated regulatory mechanisms underlying sepsis. Finally, RT-PCR was used to verify the expression of the identified miRNAs and TFs using bioinformatics.

Dataset collection and identification of DEGs
The GEO (http:// www. ncbi. nlm. nih. gov/ geo) database is a public functional genomics data repository of high-throughput gene expression data, chips, and microarrays (Edgar et al. 2002). Two datasets, GSE131761 (Martinez-Paz et al. 2020) and GSE94717 (Ge et al. 2017), met our retrieval requirement for this study. The differentially expressed mRNAs and miR-NAs between healthy controls and samples collected from patients who experienced sepsis were identified using the web tool named GEO2R (http:// www. ncbi. nlm. nih. gov/ geo/ geo2r). Adjusted P-values (adj. P) and the Benjamini and Hochberg procedure were applied to increase the reliability of the results. LogFC (fold change) values greater than 2 and adj. P-values less than 0.05 were considered statistically significant.

Functional analysis and pathway enrichment analysis of DEGs
The biological functions and signaling pathways of the identified differentially expressed genes were analyzed using the Database for Annotation, Visualization, and Integrated Discovery (DAVID). DAVID (http:// david. ncifc rf. gov) is a common tool used for GO (gene ontology) annotation and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways enrichment analysis (Huang et al. 2007). It provides a comprehensive set of tools for functional annotation of genes and proteins (Ashburner et al. 2000). P < 0.05 was considered statistically significant.

PPI network construction and identification of hub genes
In order to identify hub genes among the identified DEGs, the STRING database (Search Tool for the Retrieval of Interacting Genes; version 10.0; http:// string-db. org) was used to obtain the predicted interactions (Franceschini et al. 2013). Cytoscape (version 3.6.1) is an open source bioinformatics software platform that visualizes a PPI (protein-protein interaction) network (Smoot et al. 2011). The MCODE (Molecular Complex Detection) plugin in Cytoscape can identify highly connected proteins and significant modules from the network analyzed by STRING. Genes with degree connectivity > 8 were defined as hub genes.
Prediction and construction network of the TF-miRNA-mRNA network TargetScan (http:// www. targe tscan. org/ vert_ 72/) and miRWalk (http:// mirwa lk. umm. uni-heide lberg. de) were used to predict the biological targets of differentially expressed miRNAs (DEMs), and those that had overlap with the targets of DEGs were identified. The TransmiR database (http:// www. cuilab. cn/ trans mir) was used to predict TF-miRNA relations. Based on the regulation relationship between TFs, miRNAs, and mRNAs, the TF-miRNA-mRNA network was constructed and visualized using Cytoscape (version 3.6.1).
Cell culture and treatment Human umbilical vein endothelial cells (HUVECs) were obtained from the Cell Resource Center of the Shanghai Institute of Life Science (Shanghai, China). The cells were cultured in endothelial cell medium (ECM) (ScienCell, San Diego, CA) and were maintained at 37 °C in 5% CO2 using standard cell culture methods. After reaching a confluence of approximately 80%, HUVECs were detached using 0.25% trypsin in ethylenediaminetetraacetic acid (EDTA). HUVECs were treated with 1 μg/mL lipopolysaccharide (LPS) for 24 h.
Cell viability assay Cell Count Kit-8 (CCK-8, Dojindo, Kumamoto, Japan) was performed to measure cell viability. Approximately 5 × 10 3 transfected cells were seeded into each well of a 96-well plate. Ten microliters of CCK-8 was added to each well every 24 h for 30-min incubation. Subsequently, the absorbance values of the experimental wells were analyzed by a microplate reader at 490 nm.
Total RNA extraction, reverse transcription, and quantitative real-time PCR Total RNA was isolated from cells using TRIzol reagent (Invitrogen, Carlsbad, CA) according to the manufacturer's instructions. A reverse transcription reagent was added to the RNA sample according to the commercial kit instructions, and the resulting cDNA was used for quantitative real-time PCR (qRT-PCR). QRT-PCR was accomplished using the FastStart Universal SYBR Green Master (Rox) (Roche) in the ABI PRISM® 7300 real-time PCR system (Applied Biosystems, Foster City, CA, USA). We used melting curves to monitor nonspecific amplifications. Relative expression level was computed using 2 − ΔΔCt method.

Statistical analysis
Results were reported as mean ± standard deviation (Mean ± SD). When comparing two groups of data, the unpaired t-test or Student's t-test was used to analyze real-time gene expression analysis. One-way ANOVA was used when comparing multiple groups of data. P < 0.05 was considered statistically significant. All data and graphs were constructed and analyzed using SPSS software (version 19.0, USA) and GraphPad Prism software (version 8.0, USA).

Identification of DEGs and DEMs in patients with sepsis
A total of 146 upregulated and 62 downregulated DEGs were identified in the GSE131761 dataset, while 1 upregulated and 183 downregulated DEMs were identified in the GSE94717 dataset. The volcano plots and heatmap are shown in Fig. 1.
Gene ontology enrichment analysis and KEGG pathway analysis GO functional enrichment analysis of the identified DEGs showed that the biological processes of these DEGs were significantly enriched in immune response, innate immune response, defense response to fungi, defense response to bacteria, and response to lipopolysaccharide ( Fig. 2A). The cellular component was significantly enriched in the extracellular exosome, T-cell receptor complex, extracellular space, extracellular region, and plasma membrane (Fig. 2B). The molecular function was significantly enriched in MHC class II receptor activity, RAGE receptor binding, and peptide antigen binding (Fig. 2C). These DEGs were also significantly associated with the pathway involved in asthma, graft-versus-host disease, systemic lupus erythematosus, and allograft rejection (Fig. 2D).
PPI network construction and hub gene selection The PPI network was obtained and visualized using Cytoscape. Genes with connectivity degree > 8 were identified as hub genes. A total of 48 genes were identified as hub genes (Fig. 3).
Predicted target mRNA analysis and miRNA-mRNA network construction The top 16 significantly differentially expressed miRNAs of GSE94717 were selected for further analysis. The target genes of these miRNAs were intersected with the differentially expressed mRNAs of GSE131761. The miRNA-mRNA network was constructed and visualized using Cytoscape software (version 3.6.1) based on the negative regulation relationship between miRNAs and mRNAs. The constructed miRNA-mRNA network is shown in Fig. 4A. The intersection of 48 hub genes obtained from the above analysis and 16 target genes of differentially expressed miRNAs was taken to obtain 25 hub genes. Additionally for the dataset GSE131761 differentially expressed 208 genes and 454 potential transcription factors predicted by the TransmiR database for differentially expressed 16 miRNAs, the intersection set was taken to obtain 4 transcription factors including LEF1, MAPK14, RUNX3, and TBX21 (Fig. 4B). There are six miRNAs (miR-150-3p, miR-19b-1-5p, miR-5009-5p, miR-2362, miR4659a-5p, and miR-3681-5p) corresponding to these four transcription factors. In Fig. 4A, there are 17 hub genes that have regulatory relationships with these six miRNAs (ELANE, CCR7, MPO, CHI3L1, STOM, CEACAM8, FCAR, MMP8, CZMK, LCK, CD3G, CD8B, CEACAM1, KLRB1, MMP9, CD5, and CXCR3).
To confirm the reliability of our findings, we validated the identified TFs and miRNAs in our well-established The analysis was conducted using the DAVID (https:// david. ncifc rf. gov/ summa ry. jsp) database. GO, Gene Ontology; KEGG, the Kyoto Encyclopedia of Genes and Genome Fig. 1 Expression profiles of distinct mRNAs. A The volcano diagram of GSE131761 showed the differential expression of mRNAs between the normal and sepsis groups. B The volcano diagram of GSE94717 showed the differential expression of mRNAs between the normal and sepsis groups. C The heatmap based on GSE131761 showed the differential expression patterns of several mRNAs LPS-HUVECs model by qRT-PCR. HUVECs were administrated with different concentrations of LPS (0, 0.5, or 1 μg/ mL) at different time points (6, 12, and 24 h). As shown in Fig. 5A, incubation of HUVECs with LPS induced a dosedependent reduction in cell viability (P < 0.05). Figure 5B shows that the expression levels of inflammatory cytokines were markedly increased upon LPS stimulation (P < 0.05). According to our experimental results, RUNX3 expression was downregulated and MAPK14 expression was upregulated, consistent with our computational expression analysis (Fig. 6A). In the HUVECs model, miR-19b-1-5p and miR-5009-5p were significantly downregulated (Fig. 6B). Other transcription factors and mirna expression profiles were inconsistent with our bioinformatic analysis.

Discussion
Sepsis is a complex disorder that develops as a dysregulated host response to an infection, and is associated with acute organ dysfunction and a high risk of death (Cecconi et al. 2018). Despite advances in medical treatment, there is still no treatment method containing the dysregulated host response. Therefore, the molecular mechanisms that regulate the progression of sepsis are important discoveries that might lead to the prevention of this potentially life-threatening disease.
Transcription factors (TFs) are DNA-binding proteins that can play significant roles in regulating gene expression and in sepsis-induced organ dysfunction by inducing the expression of pro-inflammatory genes. Meanwhile, miRNAs play an important role in regulating the immune response and molecular pathology in sepsis (Taheri et al. 2020;Xu et al. 2020). Our study aimed to identify the potential network of TFs, miRNA, and mRNAs involved in the dysregulated immune response using computational bioinformatics.
RUNX3 is one of the three mammalian Runt-domain transcription factors (Levanon and Groner 2004). A previous study TFs from TransmiR database. The hub genes were intersected with the target genes of the top 20 miRNAs from TargetScan and miR-Walk databases. C TF-miRNA-mRNA regulatory networks were constructed for sepsis based on the above analysis and their inside interaction relationships. The red color represents TFs. The orange color represents miRNAs. The green color represents mRNAs has shown that RUNX3 might be involved in the development of hematopoietic malignancies and epithelial malignancies (Elinav et al. 2013). Other studies have shown that RUNX3 plays an important function in controlling immunity and inflammation. Dicken et al. and Yin et al. found that RUNX3 and its downstream target genes in immune and inflammatory cells act as protectors against a range of immune-related diseases (Dicken et al. 2013;Yin et al. 2018). Currently, there are a few reports on RUNX3 in the field of sepsis research. The downregulated expression of RUNX3 is consistent with the findings of these previous studies, suggesting that RUNX3 may protect against the development of sepsis.
The other differentially expressed TF, MAPK14, plays an important role in initiating numerous conditions, including inflammation, cardiovascular diseases, and cancer (Madkour et al. 2021). MAPK14 can be activated by environmental stress and by various proinflammatory cytokines. Therefore, with the progression of sepsis and septic shock, MAPK14 would be upregulated. Li et al. reported that miR-128-3p could enhance the protective effect of dexmedetomidine on acute lung injury (ALI) in a septic mouse model by inhibiting MAPK14 expression (Ding et al. 2020). Pan et al. have also found that miR-124 could alleviate acute lung injury symptoms of by inhibiting the activation of the MAPK signaling pathway via inhibiting MAPK14 expression (Pan et al. 2019). A previous bioinformatics study has already speculated that MAPK14 may be a hub gene of sepsis (Godini et al. 2018). In the TransmiR v2.0 database, MAPK14 is a transcription factor that can regulate several miRNAs, but this needs to be further verified experimentally.
At present, there is not much information available on miR-19b-1-5p. miR-19b-1-5p is one of the microRNAs from the miR-17-92 cluster, which can regulate the development of cardiovascular diseases (Singh et al. 2021). Additionally, it was found that the low expression of miR-19b-1-5p in platelets can lead to sustained platelet aggregation and increase the risk of future adverse cardio-cerebrovascular events. Kok et al. found that the expression of miR-19b-1-5p was also associated with aspirin insensitivity after aspirin use (Kok et al. 2016). In the TransmiR v2.0 database, the relationships between MAPK14 and miR-19b-1-5p are listed as regulation and feedback, which means the action type is unclear, and TF may be a target gene for miRNA. In our PCR assay, the expression of MAPK14 was high, while the expression of miR-19b-1-5p is low in the LPS group. This may indicate a negative regulatory relationship between miR-19b-1-5p and HUVECs were administrated with different concentrations of LPS (0, 0.5, or 1 μg/ mL) at different time points (6, 12, and 24 h). A HUVECs with LPS induced a dose-dependent reduction in cell viability (P < 0.05). B The expression level of inflammatory cytokines (IL-6, IFN-γ, IL-1β, and TNFα) by qRT-PCR. The results showed that the expression levels of inflammatory cytokines were markedly increased upon LPS stimulation (P < 0.05) MAPK14, though this would require further experimental proof. Previous reports of miR-5009-5p could not be found after an extensive literature search. However, its downregulation in the LPS-treated group suggests that it may play a role in regulating the sepsis response.
While these results are promising, they are still inconclusive. Verification experiments confirming the relevance of the hub genes were also not performed. Therefore, future in vitro and in vivo experiments are required to investigate the function and pathway of these genes in sepsis pathology. Studies with larger cohorts of patients with sepsis are also required to confirm the diagnostic and therapeutic value of the identified genes, transcription factors, and miRNAs.

Conclusions
A TF-miRNA-mRNA regulation network was constructed through the bioinformatics analysis of mRNA and miRNA expression data. After experimental validation, it was found that the TFs MAPK14 and RUNX3 and the miRNAs miR-19b-1-5p and miR-5009-5p potentially play important roles in the development of sepsis. These results contribute to a deeper understanding of the molecular mechanisms of sepsis and present potential targets for clinical treatment.
Author contribution Qinghui Fu conceived this study. Qinghui Fu designed the study. Xiaoqian Luo and Jun Hu acquired and analyzed the data. Enjiang Chen and Weina Lu contributed analysis tools. Xiaoqian Luo wrote the paper. Shi Fu and Jianfeng Zhao were of immense help in the preparation of the manuscript. All authors read and approved the final manuscript.

Declarations
Ethics approval and consent to participate Not applicable.

Consent for publication
All authors have consented this study to publication.

Competing interests
The authors declare no competing interests.