Integrated Single Cell and Bulk RNA-Seq Analysis Revealed Immunomodulatory Effects of Ulinastatin in Sepsis: A Multicenter Cohort Study

Sepsis is a leading cause of morbidity and mortality in the intensive care unit, which is caused by unregulated inflammatory response leading to organ injuries. Ulinastatin (UTI), an immunomodulatory agent, is widely used in clinical practice and is associated with improved outcomes in sepsis. But its underlying mechanisms are largely unknown. Our study integrated bulk and single cell RNA-seq data to systematically explore the potential mechanisms of the effects of UTI in sepsis. After adjusting for potential confounders in the negative binomial regression model, there were more genes being downregulated than being upregulated in the UTI group. These down-regulated genes were enriched in the neutrophil involved immunity such as neutrophil activation and degranulation, indicating the immunomodulatory effects of UTI is mediated via regulation of neutrophil activity. By deconvoluting the bulk RNA-seq samples to obtain fractions of cell types, the Myeloid-derived suppressor cells (MDSC) were significantly expanded in the UTI treated samples. Further cell-cell communication analysis revealed some signaling pathways such as ANEEXIN, GRN and RESISTIN that might be involved in the immunomodulatory effects of UTI. The study provides a comprehensive reference map of transcriptional states of sepsis treated with UTI, as well as a general framework for studying UTI-related mechanisms.

Introduction

Sepsis is a leading cause of morbidity and mortality in the intensive care unit (ICU), with an estimated mortality rate of 10 – 50% depending on the severity of the illness (1, 2). Although sepsis is initially caused by infection, subsequent organ dysfunction is the result of uncontrolled inflammatory response. Thus, strenuous effects have been made to develop immunomodulatory drugs, such as corticosteroids, topoisomerase inhibitors, rhodomeroterpene and pterostilbene (3–7). However, few agents have been proven to be clinically effective in reducing sepsis outcomes. The sepsis comprises a heterogenous population, and different subtypes of sepsis may respond differently for a particular treatment (8, 9). Therefore, it is of vital importance to understand the biological mechanisms underlying the pathogenesis of sepsis from the perspective of system biology.

Ulinastatin (UTI) is a broad-spectrum serine protease inhibitor with potent immunomodulatory effects. Clinical evidence showed that UTI had a significant effect on inflammatory biomarkers such as C-reactive protein (CRP), interleukin-6 (IL-6), and tumor necrosis factor-alpha (TNF-α), as well as reduced and prevented MODS and lowered mortality in acute pancreatitis (10, 11). However, the underlying mechanisms of the immunomodulatory effects remain controversial, and a variety of biological pathways have been identified for being associated with UTI . However, none of these studies have explored these potential mechanisms from the perspective of system biology by exploiting a holistic approach to deciphering the complexity of biological systems that starts from the understanding that the networks that form the whole of living organisms are more than the sum of their parts. The high throughput second generation RNA sequencing (RNA-Seq) provides an unbiased approach to exploring the changes in transcriptome under different diseases status or treatment conditions (12). The current study explored potential underlying mechanisms of UTI on sepsis by integrating bulk and single cell RNA-seq data. The study provides a comprehensive reference map of transcriptional states of sepsis treated with UTI, as well as a general framework for studying UTI-related mechanisms.

Methods

Study Population and Setting

The study was conducted in 13 tertiary care hospitals in mainland China from Nov 2020 to August 2021. Patients fulfilling the sepsis-3.0 criteria (suspected or documented infection plus acute increase in SOFA score > 2 points) on admission to ED were potentially eligible for the present study (13). Subjects were excluded if they met one of the following criteria: 1) end-stage cirrhosis with Child-Pugh C; 2) pregnancy; 3) concomitant malignancy, autoimmune disease; 4) patients who signed do-not-resuscitate order; 5) sepsis onset > 48 hours or had been treated in other hospitals when presenting to the participating hospitals; 6) immunosuppression such as long-term use of immunosuppressive agents, chemotherapy, corticosteroids, radiotherapy or HIV infection; and 7) acute myocardial infarction and/or pulmonary embolism. The study was approved by the ethics committee of Sir Run Run Shaw hospital (approval number: 20201014-39).

Demographics and baseline characteristics including age on admission, sex, height, weight, comorbidities and site of infection were obtained on the day of admission. Other clinical and laboratory variables such as lactate, C-reactive protein (CRP) were obtained on day 1, 3 and 5. All patients were followed for mechanical ventilation, continuous renal replacement therapy (CRRT) and vasopressor dependent days during the hospital stay. The use of UTI (Techpool Bio-Pharma Co., Ltd.) and its dosage (300,000 to 1000,000 U/day) were determined by the treating physician. The standard dosage of UTI was 100,000 U once to three times per day. However, there is a few evidence that large dose UTI can have additional benefits (14, 15). Blood samples of 2 ml were drawn on day 1, 3 and 5 after hospital admission and processed within 6 hours. The raw data were deposited in the National Genomic Data Center (https://ngdc.cncb.ac.cn/bioproject/browse/PRJCA006118).

Bulk RNA-Seq Library Preparation and Quantification

Peripheral blood mononuclear cells (PBMC) were isolated by using density-gradient centrifugation as per standard protocol. Total RNA was extracted and purified using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer’s procedure and were then stored at -80°. Then all samples were sent for library preparation and gene expression quantification (LC-BioTechnologies (Hangzhou) Co., LTD.). The remaining RNAs were reverse-transcribed to cDNA following ribosomal RNAs removal. U-labeled double-stranded DNAs was then synthesized with E. coli DNA polymerase I, RNase H and dUTP. The fragments were ligated with single-or dual-index adapters, and a size selection assay was performed with AMPureXP beads. U-labeled double-stranded DNAs were treated with heat-labile UDG enzyme. The ligated products were amplified with PCR under pre-established conditions. The average insert size for the final cDNA library was 300 bp ( ± 50 bp). Finally, we performed paired-end sequencing on an Illumina NovaSeq™ 6000 following the vendor’s recommended protocol.

The reads containing adapter contamination, low-quality bases and undetermined bases were removed using cutadapt-1.9 (cutadapt.readthedocs.io/en/stable/) (16). Then sequence quality was verified using FastQC v0.10.1 (www.bioinformatics.babraham.ac.uk/projects/fastqc/). We used HISAT2-2.0.4 to map reads to the genome of Homo sapiens obtained from the Ensembl v96 database (17). The mapped reads of each sample were assembled using StringTie-1.3.4 with default parameters (18). Then, the Gffcompare tool was used to compare and merge different gene annotations from all samples (19). After the final transcriptome was generated, the expression level of the RNAs were determined using StringTie to obtain FPKM values.

Differential gene expressions between the UTI and control group were tested by using negative binomial generalized linear models in the DESeq2 workflow (20). Clinical variables such as age, gender, SOFA, hospital days and interaction between days and UTI use were included in the model. The batch effects resulted from different hospitals were removed in the model. Weighted correlation network analysis (WGCNA) was used for finding modules of highly correlated genes, for summarizing such clusters using the module eigengene or an intramodular hub gene, for relating modules to one another and to external sample traits (using eigengene network methodology), and for calculating module membership measures (21). Methodological details for WGCNA can be found in electronic digital content. The correlation between module eigengene and UTI group was explored. Gene set enrichment analysis with gene ontology (GO) terms was performed for the modules highly correlated with the UTI group (22). All subontologies including cellular component, biological process, and molecular function were used for enrichment analysis. Regulatory mechanisms of the modules highly correlated with UTI were further explored with transcription factor (TF) enrichment analysis and microRNA correlation analysis (23, 24). Regulatory networks between TF/microRNA and candidate genes were matched using the multiMiR package (version 1.14.0).

Single Cell RNA-Seq Quantifications and Statistical Analysis

Blood samples were processed for Singleron Matrix™ Single cell RNA sequencing. Briefly, Single-cell suspensions with 1×10⁵ cells/mL in concentration in PBS (HyClone) were prepared. Single-cell suspensions were then loaded onto microfluidic devices and scRNA-seq libraries were constructed according to Singleron GEXSCOPE^® protocol by GEXSCOPE^® Single-Cell RNA Library Kit (Singleron Biotechnologies). Individual libraries were diluted to 4nM and pooled for sequencing. Pools were sequenced on Illumina HiSeq X with 150 bp paired end reads.

Raw reads were processed to generate gene expression profiles using an internal pipeline. Briefly, after filtering read one without poly T tails, cell barcode and UMI was extracted. Adapters and poly A tails were trimmed (fastp V1) before aligning read two to GRCh38 with ensemble version 92 gene annotation (fastp v2.5.3a and featureCounts v1.6.2) (25). Reads with the same cell barcode, UMI and gene were grouped together to calculate the number of UMIs per gene per cell. The UMI count tables of each cellular barcode were used for further analysis. Cell type identification and clustering analysis using Seurat (v4.0.4) and clustermole (v1.1.0) (26, 27). Rigorous quality control was performed and poor quality cells with the number of gene features >2500 or < 200, percent of mitochondria genes > 20% were removed. Standard Seurat workflow analysis was performed. Different samples were integrated with the IntegrateData function (26). The parameter resolution was set to 0.15 for FindClusters function to identify clusters. Differentially expressed genes (DEGs) between different clusters were identified with function FindMarkers, and gene set enrichment analysis (GSEA) was performed to identify enriched pathways (28).

Deconvolution of Bulk RNA-Seq Data

The bulk RNA-seq data were deconvoluted by using the CIBERSORT algorithm (29), so that the fraction of each cell component can be estimated for the bulk RNA-seq samples. CIBERSORT exploited nu–support vector regression (v-SVR) to identify cell types by estimating relative subsets of RNA transcripts. The differences in the fraction of cell types between UTI and control samples were compared and visualized with violin plots (30).

Cell-Cell Communication Inference and Intercellular Communication Networks

We inferred cell-cell communication using the CellChat (http://www.cellchat.org/) workflow (31). Specifically, the probability of cell-cell communications was inferred by integrating gene expression profile of each cell and prior knowledge of interactions between receptors, ligands, and their cofactors. Label-based mode was employed with cells labeled by the cell types obtained from the above Seurat workflow. The curated signaling molecule interaction database CellChatDB was used to assign roles for the signaling molecules and their interactions. Significant cell communications were obtained by identifying differentially over-expressed ligands and receptors for each cell group. The global communication patterns and key signals in different cell groups were identified by a pattern recognition method based on non-negative matrix factorization.

Results

Study Population

A total number of 145 sepsis patients were enrolled during the study period, including 22 patients in the UTI group and 123 in the control group. The UTI groups showed higher SOFA score [median (Q1, Q3): 7 (6-13) vs. 7 (5 - 8); p = 0.033] and serum lactate [5.2 (3.0 – 7.3) vs. 2.5 (1.6 – 3.9); p = 0.001] than the control group at baseline. Although the UTI group showed significantly higher serum lactate than the control group on admission, they had the similar vasopressor dependent days (p = 0.422), supporting the clinical and laboratory observations that UTI could mitigate septic shock (32, 33). Other baseline variables such as site of infection, comorbidities and use of MV were comparable between the two groups (Table 1).

TABLE 1

Table 1 Baseline characteristics of the UTI and control groups.

Differential Gene Expression in Bulk RNA-Seq Samples

DESeq2 standard workflow was employed to explore the differential expressed genes between the UTI and control groups. Since there were baseline difference in the cohort, we adjusted for the SOFA, age, sex and days in the negative binomial regression models. There are greater number of genes being suppressed than the number of being activated in UTI versus control groups (Figure 2A). These suppressed genes were significantly enriched in the Fc receptor signaling pathways, phagocytosis, neutrophil degranulation, and neutrophil mediated immunity (Figure 2B). Many genes involved in neutrophil mediated immunity were established to play important roles in the pathogenesis of sepsis, such as MAPK14 (3, 34), TLR2 (35, 36), FPR2 (37) and ITGAM (38). Interestingly, pathways involving neutrophil activation and degranulation were clustered together to form the largest cluster of top enriched pathways (Figure 2D). Many genes showed significant interactions between days and UTI treatment, indicating differential effect of UTI across day 1, 3 and 5 (Figure 2F). For instance, SERPINI2 showed higher expression levels in UTI group versus control group on day 1, but with declining expression levels through day 3 to 5 (Figure 2G), indicating it may take several days for UTI to take the modulatory effect.

Weighted Gene Co-Expression Network Analysis (WGCNA) of Sepsis Gene Expression Profile

WGCNA identified multiple gene co-expression networks (Figure 3A; Figures 1, 2). Most modules showed suppressed expression in the UTI group versus control group. In particular, the eigengene of the turquoise module was significantly correlated to the UTI group (correlation coefficient -0.25, p < 0.001, Figure 3B). The red module was also significantly suppressed by the UTI. The association between turquoise module and UTI group was confirmed by the correlation between module membership and gene significance for UTI (R_pearson=0.73; p < 0.001, Figure 3C). GO term enrichment analysis showed that pathways including neutrophil activation, neutrophil degranulation and leukocyte migration were among the most enriched pathways for the turquoise module (Figure 3D).

FIGURE 1

Figure 1 Schematic illustration of the integrated analysis of bulk and single cell RNA-seq analysis. PBMC, peripheral blood mononuclear cells; UTI, ulinastatin; SOFA, sequential organ failure assessment; GO, gene ontology; PS, propensity score; MDSC, Myeloid-derived suppressor cells; WGCNA, Weighted correlation network analysis; UMAP, Uniform manifold approximation and projection; pDC, Plasmacytoid dendritic cells; DC, dendritic cells; NK, natural killer.

FIGURE 2

Figure 2 Differential gene expression analysis of bulk RNA-seq data. (A) Volcano plot showing the differentially expressed genes between UTI and control samples, while adjusting for SOFA, age, sex and days. P values were adjusted for multiple comparisons. Genes with adjusted p < 0.01 and log2 fold change > 1 were labelled. (B) gene set enrichment analysis of suppressed genes for GO terms. (C) network visualization of enriched terms and associated genes. (D) tree plot showing the hierarchical clustering of enriched terms by using Jaccard’s similarity index. (E) The GSEA enrichment score for several sample pathways. The running sum deviated in the negative direction suggesting suppressed pathways. (F) Heatmap showing the coefficient matrix of the negative binomial regression model, with an interaction term between days and UTI. Top 70 genes of the largest absolute values of coefficients were displayed. (G, H) Example genes with significant interaction effects between days and UTI groups are displayed with box plots. (I) publication trends of the top enriched terms extracted from PubMed search engine.

FIGURE 3

Figure 3 Weighted gene co-expression network analysis (WGCNA) of bulk RNA-seq samples. (A) Dendrogram clustering of the co-expressed genes into modules. Modules were labeled by different colors. (B) module-trait relationship in sepsis samples. The blue color indicates negative correlation and the red color indicates the positive correlation. The numbers in each cell are Pearson’s correlation coefficients between module eigengene and trait values, with p values shown in the parenthesis. (C) Correlation between module membership and gene significance for UTI. Each dot represents a gene. Gene significance was quantified by comparing gene expression between UTI and control groups. Module membership of a given gene was quantified by correlating its expression profile with the module eigengene. (D) Gene set enrichment analysis of the turquoise module with GO terms.

Regulatory Transcription Factor and miRNA Enrichment Analysis

To further explore how the enriched pathways were regulated by upstream molecular events, we explored motif enrichment for the turquoise module (Figures 3, 4). The most significantly enriched motifs included aipale_cyt_meth:SPIB_RWWGRGGAAGTN_eDBD_meth (NES = 5.94), transfac_pro:M05934 (NES = 5.12) and dbcorrdb:SPI1:ENCSR000BGW_1:m1 (NES = 4.89). Transcription factors associated with these motifs included SPIB, ZFP64 and SPI1 (Figure 4D). Our results support previous finding that SPIB plays an important role in maintaining immune homeostasis (39). While the transcription factor Zinc Finger Protein 64 (ZFP64) has been explored in cancer metastasis (40), its role in sepsis has not been fully defined.

FIGURE 4

Figure 4 Regulatory drivers of the turquoise module. (A) Network visualization of the enriched motifs and genes. (B) Cumulative recovery curve for several example motifs. The red line represents the global mean of the number of recovered genes and the green line represents the third standard deviation (3SD). Motifs greater than the 3SD were considered statistically significant. (C) Sequence logo of enriched motifs, the relative width of a position is represented by the information content matrix. (D) Results of the top six enrichment motifs. (E) Network visualization of the top 5 miRNA and their target genes identified from the “mirtarbase” table.

Several microRNA including hsa-miR-335-5p, hsa-miR-16-5p, hsa-miR-26b-5p and hsa-miR-155-5p were identified as the key regulators of the turquoise module involving neutrophil activation and degranulation (Figure 4E).

Single Cell RNA-Seq Showing Expansion of MDSC in UTI Group

Two patients received UTI treatment were matched to corresponding two control subjects by propensity score (Figure 5A). A total of 103,870 PBMCs were harvested and passed quality control (Figures 4, 5), which were clustered into 12 cell types including Myeloid-derived suppressor cells (MDSC), neutrophil, M2 macrophage, Plasmacytoid dendritic cells (pDC), monocytes, CD4+ T cell, dendritic cells (DC), natural killer (NK) cell, B-cell, plasma cell, megakaryocyte, and erythrocyte (Figure 5B and 6). To explore the changes of cell fractions with the treatment of UTI, the bulk RNA-Seq samples were deconvoluted into cellular components using the CIBERSORT method (41). Fifty-three PBMC samples from 22 patients treated with UTI were compared to 365 control samples. While neutrophil fraction was significantly reduced with UTI treatment $(\text{standardized mean difference\hspace{0.17em}}{\widehat{g}}_{hedge}=0.42[0.15,0.69];\text{p}=0.002)$, the MDSC fraction was significantly expanded in the UTI group $[{\widehat{g}}_{hedge}=-0.95(-1.22,-0.66);\text{p}=0.002]$ (-1.22, -0.66); p = 0.002]. MDSCs constitute a heterogeneous population of immature myeloid cells that potently suppress immune responses, and is enriched in sepsis (42–44). However, the regulatory mechanism of MDSC expansion in sepsis is unknown. Our study provides evidence that UTI may take their immunomodulatory effects via regulating MDSC expansion. Transcriptomic biomarkers of each cell type were explored by comparing the differential gene expression for a specific cell type versus all other remaining cell types (Figures 5D, E). Consistent with that reported in the literature (45), transcriptomic biomarkers of MDSC included PI3 and SLPI. However, there are significant overlap of transcriptome profiles between neutrophils and MDSC, raising the question whether human MDSCs and neutrophils are actually different cell types or whether they are one plastic cell type that can functionally polarize from microbial killers to immunosuppressor cells (46).

FIGURE 5

Figure 5 Single cell RNA-seq analysis in UTI and control samples. (A) balance diagnostics of propensity score matching for samples with and without UTI administration. (B) Uniform manifold approximation and projection (UMAP) of single cell gene expression profile, stratified by UTI groups. Twelve cell types were identified for the 103,870 cells. (C) Comparisons of the cell type fractions between UTI and control groups. Cell type fractions of the bulk RNA-seq samples were inferred by deconvolution methods. The effect size of ${\widehat{g}}_{hedge}$ is calculated by dividing the mean difference with pooled standard deviation. (D) Heatmap and (E) dot plot showing the marker genes of each cell type. (F) Gene set enrichment analysis of the differentially expressed genes between UTI and control groups across all cell types.

FIGURE 6

Figure 6 Comparisons of cell-cell communications between UTI and control groups. (A) The number and strength of inferred communication links between UTI and control groups. (B) Circle plot showing differential cell-cell communication networks between UTI and control groups. The width of edges represents the relative number of interactions or interaction strength. Red (or blue) colored edges represent increased (or decreased) signaling in the UTI compared to control. (C) Heatmap showing differential number of interactions or interaction strength in the cell-cell communication network between control and UTI groups; red color indicates increased signaling in the UTI compared to control. The top colored bar plot represents the sum of column of values displayed in the heatmap. The right colored bar plot represents the sum of row of values. (D) Scatter plots showing the dominant senders (sources) and receivers (targets) in a 2D space. x-axis and y-axis are respectively the total outgoing or incoming communication probability associated with each cell group. Dot size is proportional to the number of inferred links (both outgoing and incoming) associated with each cell group. Dot colors indicate different cell groups. (E) 2D visualization of differential outgoing and incoming signaling associated with one cell group. Positive values indicate the increase in the UTI dataset while negative values indicate the increase in the control dataset. (F) Ranking signaling networks based on the information flow. Significant signaling pathways were ranked based on differences in the overall information flow within the inferred networks between sepsis and septic shock. The top signaling pathways colored red are enriched in control, and these colored green were enriched in UTI. (G) 2D visualization of the joint manifold learning of signaling networks from UTI and control datasets. Each dot represents the communication network of one signaling pathway. Dot size is proportional to the overall communication probability. (H) Circle plot showing the inferred signaling network at UTI and control groups. Edge width represents the communication probability, and the edge colors are consistent with the color of the sender cell type. (I) Comparisons of significant interactions (Ligand-Receptor pairs) between UTI and control, which contribute to the signaling from MDSC to M2 macrophage, monocyte, megakaryocytes, and neutrophil subpopulations. Dot color reflects communication probabilities and dot size represents computed p-values. Empty space means the communication probability is zero. p-values are computed from one-sided permutation test.

By comparing gene expression profiles between UTI and control groups across all cell types, some significantly enriched pathways were identified (Figure 5F). For example, the interferon-gamma and alpha response pathways were suppressed in the MDSC by UTI treatment. Many pathways display cell specific activity. While TNF-alpha signaling via NF-_κB was upregulated in monocytes, DC and CD4+ T cells, it was down-regulated in B-cells, indicating the suppression of humoral immunity and activation of cellular immunity. Interestingly, our study showed that the estrogen response in neutrophil is suppressed, consistent with findings from the bulk RNA-seq analysis showing that many neutrophil mediated processes are suppressed in the UTI group. There has been evidence that suppression of estrogen receptor can inhibit the neutrophil extracellular traps formation (47).

Difference of Cell-Cell Communications Between UTI and Control Samples

To further explore the difference of cell-cell communications between UTI and control samples, the CellChat method was employed to infer communications between cell types (31). Overall, there was greater number and strength of cell-cell communications in the UTI samples than the control samples (Figure 6A). However, the signals in neutrophil and MDSC were decreased in both the number and strength (Figure 6B). More specifically, the signals of neutrophils and MDSC as the receiver were decreased. However, signals from neutrophil (sender/source) were not significantly reduced (Figure 6C), supporting the notion that neutrophils are effector cells. In contrast, MDSC is acting as a modulatory cells that send more immunosuppressive signals to other cells in the UTI group than the control group. The specific signaling pathways were explored in more details in MDSC, neutrophil and M2 macrophage (Figure 6E). The signaling pathways such as ANNEXIN, ITGB2, and RESISTIN were found to be reduced in the UTI group. ANNEXIN may play a detrimental role in sepsis, by prolonging neutrophil survival, which is known to contribute to sepsis-mediated organ damage (48, 49). RESISTIN signaling was suppressed in both MDSC and M2 macrophage by the use of UTI, supporting the therapeutic effects of the UTI. Elevated RESISTIN signaling has been associated with adverse outcome in sepsis in the Albumin Italian Outcome Sepsis (ALBIOS) trial (50). The above mentioned signaling pathways were altered between UTI and control groups, which was also supported by the clustering of each signaling in a 2D space (Figure 6G).

Discussions

The study integrated bulk and single cell RNA-seq data to systematically explore the potential mechanisms of the effects of UTI in sepsis. After adjusting for potential confounders, such as age, gender, SOFA, hospital days and interaction between days and UTI, in the negative binomial regression model, there were more genes being downregulated than being upregulated. These down-regulated genes were enriched in the neutrophil involved immunity such as neutrophil activation and degranulation, indicating the immunomodulatory effects of UTI is mediated via regulation of neutrophil activity. Neutrophil is a kind of innate immune cell, acting as forerunners to clear the infection and resolute the inflammation during sepsis. Neutrophil extracellular traps (NETs) have been demonstrated to kill the pathogens by releasing DNA decorated with histone and granular proteins. However, unregulated NETs have a significant influence on the pathogenesis of sepsis-induced multiple organ damage, including arterial hypotension (shock), hypoxemia, renal, neurological, coagulopathy, and hepatic dysfunction (51). In our study cohort, patients in the UTI group showed more severe circulatory shock and organ dysfunction on day 1 (e.g. higher serum lactate and SOFA), but the duration of vasopressor dependent days were similar between the two groups. This results indicate beneficial effects of UTI for ameliorating the severity of organ dysfunctions. Biological process of neutrophil degranulation is closely related to the NET formation and we deduce that the effects of UTI could be mediated via reducing the incidence of NET formation. To the best of our knowledge, there has been no study exploring the effects of UTI on NET formation, and can be explored in further experimental assays.

To further corroborate the results from the bulk RNA-seq data, single cell RNA-seq was performed to identify component cell types of PBMC. The cellular components of the bulk RNA-seq samples can be decomposed with the scRNA-seq matrix as reference. Interestingly, our study found that the MDSC, which shared morphological similarity to neutrophils, were significantly expanded with the UTI treatment. It is well known that MDSC had anti-inflammatory effects that can ameliorates organ injury (45). In our study, the enriched biological processes in MDSC comparing UTI versus control samples were those associated with suppressed inflammatory responses.

Finally, we explored cell-cell communication patterns across cell types in UTI and control samples. Some conventional signaling pathways associated with immunomodulatory potency, such as ANNEXIN, GRN and RESISTIN, were identified to be differentially regulated in UTI versus control groups in our study. For example, many studies have shown that GRN signaling plays an essential role in sepsis immunity, including bacterial clearance, cell growth and survival, tissue repair, and the regulation of inflammation (52–54). Our study showed that the reduced incoming GRN signaling in neutrophils is control specific, indicating that UTI is associated with restoration of GRN signaling.

Several limitations must be acknowledged in the study. Firstly, the study was observational in nature and the treated and control groups were not fully exchangeable. However, we tried to adjust for confounding factors between the two groups by multivariable regression model in the DESeq2 pipeline. The study can be considered as hypothesis-generating and the findings need to be validated in experimental studies. Second, the subjects were enrolled from multiple centers in the study and the batch effects can confounding potential biological effects. We used DESeq2 pipeline to address this issue.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The name of the repository and accession number can be found below: China National Genomics Data Center; PRJCA006118 (https://ngdc.cncb.ac.cn/bioproject/browse/PRJCA006118).

Ethics Statement

The study was approved by the ethics committee of Sir Run Run Shaw hospital (approval number: 20201014-39). The patients/participants provided their written informed consent to participate in this study.

Author Contributions

YH, ZZ, and GZ designed research studies, SJ, LC, MY, CG, YiY, GD, and WZ conducted experiments, JZ, GH, LQ, JW, YX, JS, NW, MW, YaY, YiY, YT, JH, QB, and LX acquired data, HC and ZZ analyzed data, and ZZ and YH wrote the manuscript. All authors contributed to the article and approved the submitted version.

Funding

The study was supported by Health Science and Technology Plan of Zhejiang Province (2021KY745), National Natural Science Foundation of China (81901929) and Key Research & Development project of Zhejiang Province (2021C03071). LC received funding from Zhejiang Provincial Project of Medical and Health Technology (2022PY099), Youth Foundation of Jinhua Municipal Central Hospital (JY2020-2-10) and Key Project of Jinhua City (2021-3-037) and Key Laboratory of Emergency and Trauma (Hainan Medical University), Ministry of Education (Grant No. KLET-202118). Ulinastatin was provided by the Techpool Bio-Pharma Co., Ltd.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2022.882774/full#supplementary-material

References

1. Bauer M, Gerlach H, Vogelmann T, Preissing F, Stiefel J, Adam D. Mortality in Sepsis and Septic Shock in Europe, North America and Australia Between 2009 and 2019- Results From a Systematic Review and Meta-Analysis. Crit Care (2020) 24(1):239. doi: 10.1186/s13054-020-02950-2

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Evans L, Rhodes A, Alhazzani W, Antonelli M, Coopersmith CM, French C, et al. Surviving Sepsis Campaign: International Guidelines for Management of Sepsis and Septic Shock 2021. Crit Care Med (2021) 49(11):e1063–143. doi: 10.1097/CCM.0000000000005337

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Fang M, Zou T, Yang X, Zhang Z, Cao P, Han J, et al. Discovery of Novel Pterostilbene Derivatives That Might Treat Sepsis by Attenuating Oxidative Stress and Inflammation Through Modulation of MAPKs/NF-κb Signaling Pathways. Antioxidant (Basel) (2021) 10(9):1333. doi: 10.3390/antiox10091333

CrossRef Full Text | Google Scholar

4. Wu YM, Shi Q, Zhu PF, Ma HJ, Cui SC, Li J, et al. Rhodomeroterpene Alleviates Macrophage Infiltration and the Inflammatory Response in Renal Tissue to Improve Acute Kidney Injury. FASEB J (2021) 35(11):e21985. doi: 10.1096/fj.202100981RR

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Zhang J, Zhao H, Feng Y, Xu X, Yang Y, Zhang P, et al. Topoisomerase 2 Inhibitor Etoposide Promotes Interleukin-10 Production in LPS-Induced Macrophages via Upregulating Transcription Factor Maf and Activating PI3K/Akt Pathway. Int Immunophar (2021) 101(Pt A):108264. doi: 10.1016/j.intimp.2021.108264

CrossRef Full Text | Google Scholar

6. Cohen J, Blumenthal A, Cuellar-Partida G, Evans DM, Finfer S, Li Q, et al. The Relationship Between Adrenocortical Candidate Gene Expression and Clinical Response to Hydrocortisone in Patients With Septic Shock. Intensive Care Med (2021) 47(9):974–83. doi: 10.1007/s00134-021-06464-5

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Zhang Z, Pan Q, Ge H, Xing L, Hong Y, Chen P. Deep Learning-Based Clustering Robustly Identified Two Classes of Sepsis With Both Prognostic and Predictive Values. EBioMedicine (London) (2020) 62:103081. doi: 10.1016/j.ebiom.2020.103081

CrossRef Full Text | Google Scholar

8. Zhang Z, Zhang G, Goyal H, Mo L, Hong Y. Identification of Subclasses of Sepsis That Showed Different Clinical Outcomes and Responses to Amount of Fluid Resuscitation: A Latent Profile Analysis. Crit Care (2018) 22(1):347. doi: 10.1186/s13054-018-2279-3

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Seymour CW, Kennedy JN, Wang S, Chang CH, Elliott CF, Xu Z, et al. Derivation, Validation, and Potential Treatment Implications of Novel Clinical Phenotypes for Sepsis. JAMA (2019) 321(20):2003–17. doi: 10.1001/jama.2019.5791

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Liu D, Yu Z, Yin J, Chen Y, Zhang H, Xin F, et al. Effect of Ulinastatin Combined With Thymosin Alpha1 on Sepsis: A Systematic Review and Meta-Analysis of Chinese and Indian Patients. J Crit Care (2017) 39:259–66. doi: 10.1016/j.jcrc.2016.12.013

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Meng C, Qian Y, Zhang WH, Liu Y, Song XC, Liu H, et al. A Retrospective Study of Ulinastatin for the Treatment of Severe Sepsis. Med (Baltimore) (2020) 99(49):e23361. doi: 10.1097/MD.0000000000023361

CrossRef Full Text | Google Scholar

12. Stark R, Grzelak M, Hadfield J. RNA Sequencing: The Teenage Years. Nat Rev Genet (New York) (2019) 20(11):631–56. doi: 10.1038/s41576-019-0150-2

CrossRef Full Text | Google Scholar

13. Singer M, Deutschman CS, Seymour CW, Shankar-Hari M, Annane D, Bauer M, et al. The Third International Consensus Definitions for Sepsis and Septic Shock (Sepsis-3). JAMA (2016) 315(8):801–10. doi: 10.1001/jama.2016.0287

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Meng WT, Qing L, Li CZ, Zhang K, Yi HJ, Zhao XP, et al. Ulinastatin: A Potential Alternative to Glucocorticoid in the Treatment of Severe Decompression Sickness. Front Physiol (Lausanne) (2020) 11:273. doi: 10.3389/fphys.2020.00273

CrossRef Full Text | Google Scholar

15. Huang H, Hu PF, Sun LL, Guo YB, Wang Q, Liu ZM, et al. Treatment of Patients With Covid-19 With a High Dose of Ulinastatin. Exp Ther Med (2022) 23(2):121. doi: 10.3892/etm.2021.11044

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Kechin A, Boyarskikh U, Kel A, Filipenko M. Cutprimers: A New Tool for Accurate Cutting of Primers From Reads of Targeted Next Generation Sequencing. J Comput Biol (2017) 24(11):1138–43. doi: 10.1089/cmb.2017.0096

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Kim D, Langmead B, Salzberg SL. HISAT: A Fast Spliced Aligner With Low Memory Requirements. Nat Methods (2015) 12(4):357–60. doi: 10.1038/nmeth.3317

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL, et al. StringTie Enables Improved Reconstruction of a Transcriptome From RNA-Seq Reads. Nat Biotechnol (2015) 33(3):290–5. doi: 10.1038/nbt.3122

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Pertea G, Pertea M. GFF Utilities: GffRead and GffCompare [Version 1; Peer Review: 3 Approved] [Internet]. F1000Research (2020) 9(304). doi: 10.12688/f1000research.23297.1

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Love MI, Huber W, Anders S. Moderated Estimation of Fold Change and Dispersion for RNA-Seq Data With Deseq2. Genome Biol (2014) 15(12):550. doi: 10.1186/s13059-014-0550-8

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Langfelder P, Horvath S. WGCNA: An R Package for Weighted Correlation Network Analysis. BMC Bioinf (2008) 9:559. doi: 10.1186/1471-2105-9-559

CrossRef Full Text | Google Scholar

22. Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, et al. Clusterprofiler 4.0: A Universal Enrichment Tool for Interpreting Omics Data. Innovation (N Y) (2021) 2(3):100141. doi: 10.1016/j.xinn.2021.100141

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Aibar S, González-Blas CB, Moerman T, Huynh-Thu VA, Imrichova H, Hulselmans G, et al. SCENIC: Single-Cell Regulatory Network Inference and Clustering. Nat Methods (2017) 14(11):1083–6. doi: 10.1038/nmeth.4463

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Ru Y, Kechris KJ, Tabakoff B, Hoffman P, Radcliffe RA, Bowler R, et al. The Multimir R Package and Database: Integration of microRNA-Target Interactions Along With Their Disease and Drug Associations. Nucleic Acids Res (2014) 42(17):e133. doi: 10.1093/nar/gku631

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Liao Y, Smyth GK, Shi W. Featurecounts: An Efficient General Purpose Program for Assigning Sequence Reads to Genomic Features. Bioinformatics (2014) 30(7):923–30. doi: 10.1093/bioinformatics/btt656

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Butler A, Hoffman P, Smibert P, Papalexi E, Satija R. Integrating Single-Cell Transcriptomic Data Across Different Conditions, Technologies, and Species. Nat Biotechnol (2018) 36(5):411–20. doi: 10.1038/nbt.4096

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Satija R, Farrell JA, Gennert D, Schier AF, Regev A. Spatial Reconstruction of Single-Cell Gene Expression Data. Nat Biotechnol (2015) 33(5):495–502. doi: 10.1038/nbt.3192

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene Set Enrichment Analysis: A Knowledge-Based Approach for Interpreting Genome-Wide Expression Profiles. Proc Natl Acad Sci USA (2005) 102(43):15545–50. doi: 10.1073/pnas.0506580102

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust Enumeration of Cell Subsets From Tissue Expression Profiles. Nat Methods (2015) 12(5):453–7. doi: 10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Patil I. Visualizations With Statistical Details: The “Ggstatsplot” Approach. J Open Source Software (2021) 6(61):3167. doi: 10.21105/joss.03167

CrossRef Full Text | Google Scholar

31. Jin S, Guerrero-Juarez CF, Zhang L, Chang I, Ramos R, Kuan CH, et al. Inference and Analysis of Cell-Cell Communication Using CellChat. Nat Commun (2021) 12(1):1088. doi: 10.1038/s41467-021-21246-9

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Huang S, Guan X, Chen J, Ou Yang B. Clinical Study and Long-Term Evaluation of Immunomodulation Therapy on Trauma, Severe Sepsis and Multiple Organ Dysfunction Syndrome Patients. Zhongguo Wei Zhong Bing Ji Jiu Yi Xue (2006) 18(11):653–6.

PubMed Abstract | Google Scholar

33. Wang J, Zhou J, Bai S. Combination of Glutamine and Ulinastatin Treatments Greatly Improves Sepsis Outcomes. J Inflammation Res (2020) 13:109–15. doi: 10.2147/JIR.S234122

CrossRef Full Text | Google Scholar

34. Pan W, Wei N, Xu W, Wang G, Gong F, Li N, et al. MicroRNA-124 Alleviates the Lung Injury in Mice With Septic Shock Through Inhibiting the Activation of the MAPK Signaling Pathway by Downregulating Mapk14. Int Immunophar (2019) 76:105835. doi: 10.1016/j.intimp.2019.105835

CrossRef Full Text | Google Scholar

35. Jang JC, Li J, Gambini L, Batugedara HM, Sati S, Lazar MA, et al. Human Resistin Protects Against Endotoxic Shock by Blocking LPS-TLR4 Interaction. Proc Natl Acad Sci USA (2017) 114(48):E10399–408. doi: 10.1073/pnas.1716015114

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Wang Q-L, Xing W, Yu C, Gao M, Deng L-T. ROCK1 Regulates Sepsis-Induced Acute Kidney Injury via TLR2-Mediated Endoplasmic Reticulum Stress/Pyroptosis Axis. Mol Immunol (2021) 138:99–109. doi: 10.1016/j.molimm.2021.07.022

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Diao N, Zhang Y, Chen K, Yuan R, Lee C, Geng S, et al. Deficiency in Toll-Interacting Protein (Tollip) Skews Inflamed Yet Incompetent Innate Leukocytes In Vivo During DSS-Induced Septic Colitis. Sci Rep (2016) 6:34672. doi: 10.1038/srep34672

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Sim H, Jeong D, Kim HI, Pak S, Thapa B, Kwon HJ, et al. CD11b Deficiency Exacerbates Methicillin-Resistant Staphylococcus Aureus-Induced Sepsis by Upregulating Inflammatory Responses of Macrophages. Immune Netw (2021) 21(2):e13. doi: 10.4110/in.2021.21.e13

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Saiga H, Ueno M, Tanaka T, Kaisho T, Hoshino K. Transcription Factor MafB-Mediated Inhibition of Type I Interferons in Plasmacytoid Dendritic Cells. Int Immunol (2022) 34(3):159–72. doi: 10.1093/intimm/dxab103

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Jiang J, Zhang J, Fu K, Zhang T. Function and Mechanism Exploration of Zinc Finger Protein 64 in Lung Adenocarcinoma Cell Growth and Metastasis. J Recept Signal Transduct Res (2021) 41(5):457–65. doi: 10.1080/10799893.2020.1825490

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Steen CB, Liu CL, Alizadeh AA, Newman AM. Profiling Cell Type Abundance and Expression in Bulk Tissues With CIBERSORTx. Methods Mol Biol (2020) 2117:135–57. doi: 10.1007/978-1-0716-0301-7_7

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Janols H, Bergenfelz C, Allaoui R, Larsson AM, Rydén L, Björnsson S, et al. A High Frequency of MDSCs in Sepsis Patients, With the Granulocytic Subtype Dominating in Gram-Positive Cases. J Leukoc Biol (2014) 96(5):685–93. doi: 10.1189/jlb.5HI0214-074R

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Agrati C, Sacchi A, Bordoni V, Cimini E, Notari S, Grassi G, et al. Expansion of Myeloid-Derived Suppressor Cells in Patients With Severe Coronavirus Disease (COVID-19). Cell Death Diff (2020) 27(11):3196–207. doi: 10.1038/s41418-020-0572-6

CrossRef Full Text | Google Scholar

44. Marais C, Claude C, Semaan N, Charbel R, Barreault S, Travert B, et al. Myeloid Phenotypes in Severe COVID-19 Predict Secondary Infection and Mortality: A Pilot Study. Ann Intensive Care (2021) 11(1):111. doi: 10.1186/s13613-021-00896-4

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Trikha P, Carson WE. Signaling Pathways Involved in MDSC Regulation. Biochim Biophys Acta (2014) 1846(1):55–65. doi: 10.1016/j.bbcan.2014.04.003

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Aarts CEM, Kuijpers TW. Neutrophils as Myeloid-Derived Suppressor Cells. Eur J Clin Invest (2018) 48 Suppl 2:e12989. doi: 10.1111/eci.12989

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Flores R, Döhrmann S, Schaal C, Hakkim A, Nizet V, Corriden R, et al. The Selective Estrogen Receptor Modulator Raloxifene Inhibits Neutrophil Extracellular Trap Formation. Front Immunol (2016) 7:566. doi: 10.3389/fimmu.2016.00566

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Toufiq M, Roelands J, Alfaki M, Syed AhamedKabeer B, Saadaoui M, et al. Annexin A3 in Sepsis: Novel Perspectives From an Exploration of Public Transcriptome Data. Immunology (2020) 161(4):291–302. doi: 10.1111/imm.13239

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Mui L, Martin CM, Tschirhart BJ, Feng Q. Therapeutic Potential of Annexins in Sepsis and COVID-19. Front Pharmacol (2021) 12:735472. doi: 10.3389/fphar.2021.735472

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Bonaventura A, Carbone F, Vecchié A, Meessen J, Ferraris S, Beck E, et al. The Role of Resistin and Myeloperoxidase in Severe Sepsis and Septic Shock: Results From the ALBIOS Trial. Eur J Clin Invest (2020) 50(10):e13333. doi: 10.1111/eci.13333

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Kumar S, Payal N, Srivastava VK, Kaushik S, Saxena J, Jyoti A. Neutrophil Extracellular Traps and Organ Dysfunction in Sepsis. Clin Chim Acta (2021) 523:152–62. doi: 10.1016/j.cca.2021.09.012

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Brandes F, Borrmann M, Buschmann D, Meidert AS, Reithmair M, Langkamp M, et al. Progranulin Signaling in Sepsis, Community-Acquired Bacterial Pneumonia and COVID-19: A Comparative, Observational Study. Intensive Care Med Exp (2021) 9(1):43. doi: 10.1186/s40635-021-00406-7

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Tian G, Jin X, Wang Q, Ye T, Li G, Liu J. Recent Advances in the Study of Progranulin and Its Role in Sepsis. Int Immunophar (2020) 79:106090. doi: 10.1016/j.intimp.2019.106090

CrossRef Full Text | Google Scholar

54. Terryn J, Verfaillie CM, Van Damme P. Tweaking Progranulin Expression: Therapeutic Avenues and Opportunities. Front Mol Neurosci (2021) 14:713031. doi: 10.3389/fnmol.2021.713031

PubMed Abstract | CrossRef Full Text | Google Scholar