We used GEOquery package of R software (version 4.0.0, http://rproject.org/) (Davis and Meltzer, 2007) to download the expression profile of reliable sample source for patients with sepsis-induced lung injury (SILI) from GSE28750 dataset from GEO database (https://www.ncbi.nlm.nih.gov/geo/) (Sutherland et al., 2011). The chip platform was based on the GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array, and the samples in the dataset were all derived from Homo sapiens. There was a total of 10 sepsis samples and 20 normal samples in GSE28750 dataset. We read the original data of GSE28750 by “affy” package (Gautier et al., 2004), obtained the gene expression matrix of the data set, and normalized the data by “limma” package. Setting species as H. sapiens, the single cell sequencing (scRNA-seq) of sepsis was downloaded from Single Cell Portal database (https://singlecell.broadinstitute.org/single_cell). The Seurat R package (Version 4.0) was used to process single-cell data, and 29 samples were included to create Seurat objects for our analysis (Butler et al., 2018; Stuart et al., 2019; Hao et al., 2021; Stuart et al., 2021). We used the doubletFinder_v3 function to remove the two-cell effect and the LogNormalize method to normalize the data. After controlling for the relationship between average expression and dispersion, highly variable genes were identified in individual cells. Whereafter, variable genes were identified as inputs, Principal Component Analysis (PCA) were carried out and significant principal components based on the ElbowPlot function were finally identified (Supplementary Figure S1). According, the “elbow criterion,” the first 15 principal components were selected as statistically significant inputs of the t-Distributed Neighbor Embedding (t-SNE).

We selected differentially expressed genes (DEGs) of GSE28750 through the limma package (Ritchie et al., 2015). The ggplot2 package was used to map the volcano of DEGs (Ito and Murphy, 2013) and the pheatmap package was used to map the heatmap of DEGs (Ito and Murphy, 2013) to display the distribution of DEGs. In addition, we used the removeBatchEffect function in the limma package to remove the batch effect and DEGs met the requirements of the adj.p.value < 0.05 and |log2FC|>1. We downloaded and collated lists of immune-related genes from Pubmed and MSigDB databases, and then manually searched the original literature and reviews containing these genes (Reyes et al., 2020). Finally, after excluding extraneous genes and adding other reported genes, DEGs of sepsis were screened out and intersected with immune-related genes to obtain immune-related differentially expressed genes (IRDEGs).

First, we used the FindClusters function to cluster the cells and identify the cell types of the clusters.In order to verify the annotation of single cell, HumanPrimaryCellAtlasData was used to annotate the cell types through SingleR (Version 1.8.1) (Aran et al., 2019). In order to annotate cells more accurately, we sorted out marker genes of all kinds of cells according to previously published literature to identify cell types, and further searched for differential marker genes among cell subsets.

Cell differentiation was inferred using the Monocle package of R (version 2.22.0) (Trapnell et al., 2014). An integrated gene expression matrix from each cell type was first exported from the Seurat object to Monocle to construct the cell data set. We used the variable genes defined by dispersionTable function and then sequenced the cells using the setOrderingFilter function. Finally, the Darter method was used for dimensionality reduction, and the orderCells function was used to estimate the arrangement of cells along the trajectory. Based on the clustering characteristics and marker gene analysis, the differentiation time locus of cell subsets in single cell data set was obtained.

Gene set variation analysis (GSVA) is a non-parametric, unsupervised algorithm (Hänzelmann et al., 2013). We analyzed the data based on the GSVA package (version 1.42.0) of R. GSVA algorithm transformed gene expression data from the expression matrix of a single gene as a feature to the expression matrix of a specific gene set as a feature (Hänzelmann et al., 2013). The gene set corresponding to each feature was calculated using rank statistics similar to K-S test, and the expression matrix was converted into Enrichment Score (ES) matrix for feature (Hänzelmann et al., 2013). GSVA enrichment score corresponding to each sample for each feature could be obtained, which would facilitate further statistical analysis (Hänzelmann et al., 2013).

In an animal experiment, 40 male C57BL6 mice aged 6–8 weeks were employed (Purchased from Guangdong Medical Laboratory Animal Center). After 7-day adaptive feeding, these rats were divided into four groups: sham group (n = 10), sham + anti-IL-7 group (n = 10), sepsis group and sepsis + anti-IL-7 group (n = 10). The rats in sepsis group were anesthetized by isoflurane, and a 1-cm midline incision was made along the abdominal line. The abdominal cavity was then opened layer by layer to find the ileum. A tight ligature was made at 1 cm from the distal end of ileum, and two punctures were made 0.5 cm away from the distal end with a 7-gauge needle. A small amount of feces was squeezed into the abdominal cavity. The ileum was placed back to its normal position in the abdominal cavity and the abdomen was then closed in two layers. After the operation, 1 mL of pre-heated physiological saline was injected subcutaneously to each mouse for fluid supplementation, and 10 mg/kg of tramadol hydrochloride was injected subcutaneously as an analgesic. It was injected once every 12 h within 48 h after the operation. Except for the blind intestinal ligation and perforation, the other operations were the same as those of the surgical group in the sham operation group. 48 h later, the septic rats showed symptoms of depression, decreased activity, and decreased appetite. After the modelling was completed, IL-7 (2.5 μg/day) was administered i.p. daily for 5 days, and death was performed on the sixth day. The spleen tissues were collected after the mice were sacrificed. Recombinant human IL-7 was purchased from Abcam Company and diluted with sterile distilled water. The animal experiments were performed following the National Institutes of Health Guide for the Care and Use of Laboratory Animals and approved by the Animals Care and Use Committee of Sun Yat-sen University (Guangzhou, China).

We extracted proteins from the spleen of mice. Tissue samples were lysed in RIPA buffer adding protease and phosphatase inhibitors. Protein concentrations of the supernatants were determined after centrifugation using BCA Protein Assay Kit (Beyotime). We performed Western blotting assays according to the following steps: 1) Proteins were separated by electrophoresis on NuPAGE^® 4%–12% Bis-Tris Gel (Invitrogen, Life Technologies) for 40 min at 200 V and transferred onto polyvinylidene difluoride (PVDF) membranes (Invitrogen, Thermo Fisher; 30 V, 1 h); 2) The membranes were blocked with Blocking Solution of the Western Breeze^® Chromogenic Western blot Immunodetection Kit (Invitrogen, Thermo Fisher) for 30 min on a rotary shaker at room temperature; 3) The membranes were then incubated overnight at 4°C with the following primary antibodies: GAPDH (ProteinTech, Wuhan, China, 1:5000), IL-7R (Abcam, 1:1000), and CD4 (Abcam, 1:1000). Goat anti-rabbit IgG or goat anti-mouse IgG (Invitrogen) were used as the secondary antibodies; 4)Target proteins were visualized by Molecular Imager ChemiDoc XRS System (Bio-Rad) with super Electro-Chemi-Luminescence (ECL) plus kits (Beyotime); 5) The protein bands were analyzed with Image-Pro Plus 6.0 software (Media Cybernetics); 6) Relative protein expression levels were expressed as the ratio of the band intensity of the target protein to that of GAPDH.

Apoptosis was quantified by flow cytometry using the TUNEL assay. Flow cytometric analysis (50,000–100,000 events/sample) was per_x0002_formed on FACScan. Tissues were then rinsed with PBS and treated with a TUNEL reaction mixture according to the kit instructions. After rinsing in PBS, the converter-POD was added to the tissue, covered with a glass slide or sealing film, and reacted in a dark wet box for 30 min.

The procedure Elisa (Invitrogen, Thermo Fisher) assays procedure were roughly as follows: dilute the corresponding protein with a coated solution, seal the enzymic label plate with plastic film, and place overnight in a 4°C refrigerator. The enzymic label plate was evenly wrapped at room temperature, washed with PBST lotion for 3 times, added BSA sealer, sealed with plastic film, and placed in a 37°C water bath for 1 h. They were incubated with primary and secondary antibodies successively, washed with PBST lotion for 3 times, and added TMB and 3% H₂O₂ at 37°C for 10 min of dark color development. Add H2SO4 to terminate the color rendering. Then, the optical density (OD) of each well was measured at the wavelength of 450 nm by enzyme-labeled instrument.

We first drew up the overall technical route of this study, which was shown in Figure 1.

The data of DEGs were summarized and sorted out according to GEO data platform, and the data of the analyzed sepsis samples and normal samples were compared. We removed the probes corresponding to multiple molecules from one probe, and when the probes corresponding to the same molecule were encountered, only the probe with the largest signal value was retained. After data preprocessing for GSE28750, a total of 21,655 molecules were filtered, of which 1198 molecules met |log2(FC) |>1 and p.adj<0.05. Under the above threshold (|log2(FC) |>1 and p.adj<0.05), compared with the control group, 599 genes were upregulated, and 599 genes were downregulated in the sepsis group. We used R software to extract the mutual differentially expressed genes from the gene expression matrix, as shown in the heat map (Figure 2A) and volcano map (Figure 2B). We manually searched the original literature and reviews containing immune-related genes, excluded irrelevant genes and added other reported genes, and intersected them with the list of differentially expressed genes with p-value less than 0.05, thus obtaining a total of 64 genes, as shown in the Venn diagram (Figure 2C). For the unsupervised analysis, we plotted the single-cell transcriptome on t-distributed Stochastic Neighbour Embedding (tSNE) diagram (Figures 2D, E). T cells, NK cells, Mono cells, Megakaryocyte cells, DC cells, and B formed six separate clusters on the tSNE diagram (Figure 2E). Single-cell differential expression analysis (SCDE) identified the genes that are differentially expressed between different cell types (Figure 2F).

To identify differences in biological processes between sepsis samples and normal samples, based on gene expression profile data from the GSE28750 dataset, we carried out gene set enrichment analysis (GSEA) using R packages (Table 1; Figures 3A–F). GESA results showed that sepsis samples and normal samples were related to the biological phenomena as below: Neutrophil Degranulation (Figure 3A), Modulators of Tcr Signaling and T Cell Activation (Figure 3B), IL 17 Pathway (Figure 3C), T Cell Receptor Signaling Pathway (Figure 3D), Ctl Pathway (Figure 3E), Immunoregulatory Interactions Between A Lymphoid and A Non-Lymphoid Cell (Figure 3F).

We performed GO analysis of immune-related genes, the results showed that intersection genes were mainly associated with T cell receptor binding, immune receptor activity, SH2 domain binding, signaling receptor complex adaptor activity, chemokine activity and other biological phenomena (Figure 4A). The results showed that intersection genes were mainly associated with T cell receptor binding, immune receptor activity, SH2 domain binding, signaling receptor complex adaptor activity, chemokine activity and other biological phenomena. The results of KEGG analysis showed that intersection genes were enriched in Primary immunodeficiency, Hematopoietic cell lineage, PD-L1 expression and PD-1 checkpoint pathway in cancer, T cell receptor signaling pathway, PD-L1 expression and PD-1 checkpoint pathway in cancer, T cell receptor signaling pathway, and Cytokine-cytokine receptor interaction pathways (Figures 4B, C). Single-cell GSVA analysis showed that the function of immunorelated cells in sepsis was mainly enriched in ABC transporters, Acute myeloid leukemia, Adherens Junction, Adipocytokine Signaling Pathway, Alanine Aspartate and Glutamate Metabolism (Figure 4D). Finally, we identified the key genes to be mainly concentrated in T cells after Aucell scoring, Ucell scoring, ssgsea scoring and singscoring (Supplementary Figure S2).

We used STRING tool to conduct protein-protein interaction (PPI) analysis for immune-related genes and visualized the number of interactions between each immune-related gene. The larger the genes involved, the larger the degree of nodes, and the thicker the lines, the larger the number of betweenness (Figure 5A). Subsequently, we utilized the Maximal Clique Centrality (MCC), Maximum neighborhood component (MNC), and Density of Maximum Neighborhood Component (DMNC) algorithms based on cytoHubba plug-in of Cytoscape software to screened seven hub genes (Figures 5B–E). Intersection of hub genes screened by the above three algorithms was selected and the screening results described the seven hub genes as CD28, CD3D, CD2, CD4, IL7R, LCK and CD3E (Figures 5B–E).

We visualized the expression values of the other six hub genes except CD2 in the dataset GSE28750 and showed the differences (Figures 6A–F). The results showed the expression level of CD3D (Figure 6A), CD3E (Figure 6B), CD4 (Figure 6C), CD28 (Figure 6D), IL7R (Figure 6E), and LCK (Figure 6F) was statistically different between sepsis and normal groups. The lower expression of the hub genes was observed in sepsis samples, which was consistent in all the six hub genes (Figures 6A–F).

Based on CIBERSORT algorithm, the immune cell infiltration analysis was conducted on the sepsis samples and control samples in the GSE28750 data, and the proportion of immune cells in each sample was analyzed to obtain the proportion of 22 kinds of immune cells in the sepsis samples and control samples (Figure 7A). Further, we carried out an intergroup comparison of 22 types of immune cells between sepsis samples and control samples (Figure 7B). We observed the significant difference of B cells naïve, plasma cells, T cells CD8, T cells CD4 naïve, T cells CD4 memory resting, T cells CD4 memory activated, T cells gamma delta, NK cells resting, Monocytes, Macrophages M0, Eosinophils, Neutrophils cells between sepsis samples and control samples (Figure 7B). Further, the correlation between hub genes and immune cells were analyzed according to gene expression and the abundance of corresponding immune cells, and the interaction relationship was visualized. It could be seen that IL7R was negatively correlated with T cells CD4 memory activated (Figure 7C), while positively correlated with T cells CD4 memory resting (Figure 7D) and T cells CD8 (Figure 7E). LCK was negatively correlated with T cells CD4 memory activated (Figure 7F), while positively correlated with T cells CD4 memory resting (Figure 7G) and T cells CD8 (Figure 7H).

According to 3.4 analysis results showed that the differential gene score in T cells was higher. Next, we were able to capture the differentiation of T cells into Memory T cells, NK T cells, CD4⁺ T cells, CD8⁺ T cells, Th T cells and Treg T cells (Figure 8A) in the Trajectory analysis for T cells. The results of pseudotime analysis were displayed in Figures 8B–E.

We found that the cell counts of CD4 T cells were highest in sham group while lowest in sepsis group (Figure 9A). IL-7 antibody treatment could increase the counts of CD4 T cells (Figure 9A), which may be the cause of increased apoptosis (Figure 9B). The concentration of IFN-γ in peripheral blood were lowest in sepsis group and IL-7 antibody treatment could increase the concentration of IFN-γ (Figure 9C). The concentration of IL-6 in peripheral blood were highest in sepsis group and IL-7 antibody treatment could decrease the concentration of IL-6 (Figure 9D). IL-7 antibody treatment could increase the concentration of IL-10 (Figure 9E). The concentration of TNF-γ in peripheral blood were highest in sepsis group and IL-7 antibody treatment could decrease the concentration of TNF-α (Figure 9F). The concentration of IL-7R in peripheral blood were highest in sepsis group and IL-7 antibody treatment could decrease the concentration of IL-7R (Figure 9G). IL-7 antibody treatment could also decrease the expression of IL-7R and CD4 (Figures 9H–J).