Single-cell transcriptome analysis reveals cellular reprogramming and changes of immune cell subsets following tetramethylpyrazine treatment in LPS-induced acute lung injury

The LPS-induced ALI model in mice

A total of 12 six-to-eight-week adult male C57BL/6J mice (18–22 g) were purchased from the Shanghai Model Organisms Center, Inc. and housed in a specific-pathogen-free facility under a light/dark cycle of 12 h/12 h, with access to food and water ad libitum. The mice were housed at a density of five mice per cage to ensure adequate space for activity and optimal housing conditions. No mice were excluded from the analysis. The experimental procedures were conducted in accordance with the guidelines and the Animal Welfare Act Regulations of Shanghai Pulmonary Hospital, Tongji University. The animal treatment protocol followed the methods described in a previous study (Jiang et al., 2021). Briefly, the mice were randomly divided into four groups of three mice: a control group (injected intraperitoneally with 0.2 ml of 0.9% saline) and a model group administered with intraperitoneal injection of LPS in saline (0.2 ml) at a dosage of 30 mg/kg to induce ALI. The LPS+Dex group and LPS+TMP group were treated with intraperitoneal injection of either 0.2 ml of Dex (5 mg/kg) or 0.2 ml of TMP (50 mg/kg) 1 h before LPS administration. After a 72-h treatment period, the mice were anesthetized with an intraperitoneal injection of ketamine (80 mg/kg) and xylazine (10 mg/kg), then lung tissues were collected for subsequent analysis. All the animals were bred under specific pathogen-free conditions Shanghai Pulmonary Hospital Laboratory Animal Center, in accordance with hospital guidelines (Ethics approval number: Q20-269Y). All animal studies were approved by the Institutional Animal Care and Use Committee of Tongji University.

Histological analysis

Lung tissues were isolated, fixed in 4% paraformaldehyde for 48 h and dehydrated using different concentrations of ethanol solution, from 70% to 100% for 40 min. The samples were then embedded in paraffin and cut into 4 μm sections. Hematoxylin-eosin (H & E) staining was performed on the sections. The slides were dyed with hematoxylin for 5 min and then they were rinsed with water, and then after returning to a blue coloration, the slides were stained with eosin for 1–3 min. The degree of lung damage was assessed by two independent technicians who were blinded to the experimental groups using a recently published criterion (Qiang et al., 2020). The histological images were observed using an optical microscope (Nikon, Tokyo, Japan).

Enzyme-linked immunosorbent assay

Blood samples were obtained from the orbital sinus and these were allowed to clot for 30 min at 37 °C. The samples were then centrifuged at 3,500 rpm for 10 min at 4 °C. The concentrations of IL-1β and TNF-α were determined using commercially available enzyme-linked immunosorbent assay (ELISA) kits (Abcam, Cambridge, United Kingdom), by following the manufacturer’s protocols. The absorbance at 450 nm was measured at the end of the procedures using a microplate reader.

Cell counting in BALF

The mice were anesthetized, and a surgical incision was made in the anterior neck to expose the trachea. Bronchoalveolar lavage fluid (BALF) was collected by instilling and withdrawing 0.5 mL of saline three times through an endotracheal tube. The collected BALF samples were combined and centrifuged at 3,500 rpm for 10 min at 4 °C. The resulting cell pellet was resuspended, and red blood cells were removed by using a lysis buffer. This process was repeated twice, with centrifugation at 3,500 rpm for 10 min at 4 °C each time. The final cell pellet was resuspended in 100 µL of saline. The cell types were identified and then they were quantified by using Wright-Giemsa staining.

Immunohistochemistry

The lung tissues were fixed in formalin, embedded in paraffin and then they were sectioned to 4 µm thickness. The sections were mounted on glass slides and heated at 65 °C overnight. After de-paraffinization and rehydration, the endogenous peroxidase activity was blocked by using 3% hydrogen peroxide for 20 min. To minimize non-specific binding, the slides were incubated with 5% normal goat serum for 1 h at room temperature. The sections were then incubated overnight at 4 °C with primary antibodies against CD45, CD8 and IL-6 (ab40763, ab217344 and ab9324, respectively; Abcam, Cambridge, UK) at a 1:250 dilution. After washing with TBST, the slides were incubated with HRP-conjugated goat anti-rabbit secondary antibody for 1 h at room temperature. The sections were then stained with DAB and counterstained with hematoxylin.

Single-cell capture, cDNA library preparation and sequencing processing

Twelve mice were randomly assigned to four experimental groups: Control, LPS, LPS+Dex, and LPS+TMP for single-cell transcriptomic analysis. Single-cell sequencing achieved an average depth of 31,987 reads per cell. The statistical power of this experimental design, calculated using RNASeqPower (https://rodrigo-arcoverde.shinyapps.io/rnaseq_power_calc/), was determined to be 0.91. Fresh lung tissues were resected from mice and stored in tissue preservation solution (21903-10; Shanghai Biotechnology Corporation, Shanghai, China). The tissues were placed into a sterile culture dish and rinsed with 5 mL of 1 X DPBS precooled at 4 °C. This was repeated 1–2 times to remove any residual tissue protective solution. The moistened tissues were placed in a new 1.5 mL EP tube with 500 μL of digestive solution and then these were cut into small pieces. The shredded tissues were transferred into a 50 mL centrifuge tube containing 5 mL of digestive solution, and this was gently shaken at 37 °C in a water bath. The first sample was collected after 15 min of digestion, and fresh digestive solution was added and the tissues was digested for a further 20 min. The digested cell suspensions were pooled and these were passed twice through pre-wetted 40 μm filter membranes. The filtrates were collected into new 50 mL centrifuge tubes and centrifuged at 2,000 rpm at room temperature for 10 min. The supernatants were discarded and 2 mL of precooled red blood cell lysis solution was added to the pellets. The cell suspensions obtained were gently dispersed and they were incubated at room temperature for 3 min. After the red blood cells were lysed, the remaining cells were centrifuged at 1,500 rpm for 5 min at room temperature. The supernatants were discarded and the cell pellets were re-suspended in PBS containing 2% FBS.

The dissociated cells were then stained with 0.4% trypan blue (Cat. no. 14190144; Thermo Fisher Scientific, Waltham, MA, USA) to check for viability on a Countess® II Automated Cell Counter (Thermo Fisher Scientific, Waltham, MA, USA) and then the number of cells obtained were counted. A total of 10× library preparation and sequencing bead solution with unique molecular identifier (UMI) and cell barcodes was loaded into the cell suspension close to the saturation point. This ensured that each cell was paired with a bead in a Gel Beads-in emulsion (GEM). After addition of cell lysis buffer, the poly-adenylated RNA molecules hybridized to the beads. The beads were retrieved into a single tube for reverse transcription. After cDNA synthesis, each cDNA molecule was tagged at the 5′-end (that is, the 3′-end of a messenger RNA transcript) with an UMI and a cell label indicating its cell of origin. The 10× beads that were obtained were then subjected to second-strand cDNA synthesis, adaptor ligation and universal amplification. Sequencing libraries were prepared using randomly interrupted whole-transcriptome amplification products to enrich the 3′-ends of the transcripts linked with the cell barcode and UMI. All the remaining procedures, including the library construction, were performed according to the standard manufacturer’s protocol (CG000206 RevD). The sequencing libraries were quantified using a High Sensitivity DNA Chip (Agilent, Santa Clara, CA, USA) on a Bioanalyzer 2100 and the Qubit High Sensitivity DNA Assay (Thermo Fisher Scientific, Waltham, MA, USA). The libraries were sequenced on a NovaSeq 6000 (Illumina, San Diego, CA, USA) system by using 2x150 chemistry.

Single-cell RNA-seq data processing and filtering

The single-cell (sc) RNA-seq data were processed using the Cell Ranger 3.0.1 pipeline. FASTQs generated from the Illumina sequencing output were aligned to the mouse genome, version mm10, using the STAR algorithm. Next, the gene-barcode matrices for each sample was generated by counting UMIs and filtering the non-cell associated barcodes. Finally, we generate a gene-barcode matrix containing the barcoded cells and gene expression counts. This output was then imported into the Seurat (v3.0.2) R toolkit for quality control and downstream analysis of our scRNA-seq data. All functions were run with default parameters, unless otherwise specified. We excluded cells with fewer than 200 or more than 5,000 detected genes, and selected each gene had to have at least one UMI aligned in at least three cells. The expression of mitochondria genes was calculated using the PercentageFeatureSet function in the Seurat software package. To remove low activity cells, those with more than 10 percent of mitochondrial expressed genes were excluded. The normalized data (obtained using the NormalizeData function in the Seurat package) was performed to extract subsets of variable genes. The variable genes were identified by controlling for the strong relationship between variability and average expression. Next, we integrated the data from different samples after identifying ‘anchors’ between our datasets by using the FindIntegrationAnchors and IntegrateData in the Seurat software package.

Clustering, visualization and identification of cell types

We performed principal component analysis (PCA) and reduced the data to the top 30 PCA components after scaling the data. We visualized the clusters on a 2D map produced by using uniform manifold approximation and projection (UMAP). Identification of the cell types and subtypes was by nonlinear dimensional reduction. Cells were clustered using graph-based clustering of the PCA reduced data by using the Louvain method after computing a shared nearest neighbor graph. For sub-clustering, we applied the same procedure of scaling, dimensionality reduction and clustering to the specific set of data and this process was usually restricted to one type of cell). For each cluster, we used the Wilcoxon rank-sum test to find the significant deferentially expressed genes when comparing the remaining clusters. SingleR and known marker genes were used to identify the cell types in our samples.

Statistical analysis

Statistical analysis was conducted using the R 4.0.3 software package. The two-sided Wilcoxon rank-sum test and the one-way ANOVA followed by Dunnett’s test were used to evaluate the statistical significance between groups. A p-value of 0.05 or less was considered statistically significant.

LPS-induced ALI is mitigated by treatment with TMP

We first established a model of ALI by utilizing LPS and administered either Dex or TMP to determine their protective effects. The traditional anti-inflammatory drug Dex was used as a positive control. Subsequently, we collected lungs and dissociated tissues for scRNA-seq analysis. Integrated analysis and expression profiling was conducted to investigate the potential mechanism of TMP protection against ALI. Histological examination with H & E staining demonstrated that LPS-induced ALI, characterized by pulmonary alveolar wall thickening and inflammatory cell infiltration, was attenuated significantly by either Dex or TMP treatment (Fig. 1A). We also determined the protective effects of TMP on the immune inflammatory response in the lung tissues following LPS treatment, by staining for immune cell markers and measurement of cytokine production. The results showed that the LPS induced the expression of CD45, CD8 and IL-6, which was significantly mitigated by treatment with either Dex or TMP (Fig. 1B). Consistent results were obtained when determining the levels of inflammation-related cytokines in the mouse blood samples. We found that the concentrations of TNF-α and IL-1β in the blood, which were significantly raised by LPS, were subsequently decreased by treatment with TMP (Figs. 1C and 1D). In addition, LPS led to a significant increase of total cell number in the BALF, which was attenuated by treatment with either Dex or TMP (Fig. 1E). These results indicated that TMP could efficiently protects against LPS-induced ALI.

scRNA-seq analysis of LPS-damaged lung tissues after TMP treatment

After removing low-quality scRNA-seq cells and reducing the possible batch effects by using the Harmony algorithm, we obtained 116,213 high-quality scRNA-seq profiles of various cell types from all our samples. As shown in Fig. 2A, uniform manifold approximation and projection (UMAP) analysis of cell types revealed seven major cell clusters: myeloid cells (48,650), B cells (30,458), T cells (12,072), NK cells (10,605), fibroblasts (6,252), endothelial cells (4,101) and epithelial cells (4,065) (Fig. 2A). The canonical marker genes used for dividing the major cell clusters are presented in Figs. 2B and 2C. We found that there were no unique cell subpopulations and the seven major cell clusters were present in different groups (Figs. 2D and 3A). Quantification of the cell clusters revealed that LPS significantly increased the percentage of myeloid and NK cells, while reducing the percentage of B and T cells. However, either Dex or TMP treatment did not show any protective effects against LPS with respect to the cell cluster percentages. Similar results were observed for the percentage of endothelial and epithelial cells (Figs. 3B and 3C). These results suggest that TMP may exert its protective effects by influencing specific cell subpopulations.

Identification of changes in immune cell subsets after TMP treatment in LPS-induced ALI

To investigate whether any immune cell sub-clusters were involved in TMP protection against the immune response induced by LPS, we classified the three immune cell types (T, myeloid and NK cells) into sub-clusters (Ren et al., 2021). As shown in Fig. 4A, the T cells were divided into nine sub-clusters based on the marker genes, and a heatmap analysis of the top three genes for each T cell cluster is presented in Figs. 4B. We found that LPS significantly increased the percentage of CD8+ naive T (CD8+Tn) cells, while decreasing the percentage of follicular helper T (Tfh) cells. Treatment with either Dex or TMP significantly decreased the number of CD8+Tn cells, with both drugs having similar inhibitory effects. TMP treatment also increased the number of Tfh cells when compared to the LPS group, while Dex did not have this effect (Fig. 4C). In other cell sub-clusters, LPS also increased the cell number of Double Negative T cells (DNT), Effector T cells (Teff) and Regulatory T cells (Treg), whereas there were no significantly inhibitory effects of TMP on LPS-induced cell number change found in these cell subsets (Fig. S1).

Additionally, the myeloid cells were divided into 11 sub-clusters (Fig. 5A). The top three genes for each cell clusters of myeloid cells are presented in a heatmap (Fig. 5B). We found that LPS-increased a subset of cells within the C3 sub-cluster, named the Vcan+ monocyte, and this was reduced by Dex treatment. However, TMP had an only slight inhibitory effect on the percentage of the Vcan+ monocytes (Fig. 5C). In the other cell subsets of myeloid cells, LPS led to increased percentages of inflammatory Macrophages (Inflam-Mφ) and resident Tissue Macrophages (RTM-Mφ) cells, whereas no significantly inhibitory effects of TMP were found in these two cell subsets (Fig. S2). Furthermore, the NK cells were divided into 10 sub-clusters (Fig. 6A), and heatmap analysis of the top three genes for each of these was performed (Fig. 6B). LPS significantly induced the Siva+ NK cells and this is a subset of cells within the C5 sub-cluster. This effect was attenuated by either Dex or TMP treatment, particularly by Dex treatment (Fig. 6C). In the other cell subsets of NK cells, LPS induced the increased number of Khdc1a-positive, Klra1-positive and Mki67-positive NK cells, while treatment with either Dex or TMP only slightly inhibited the effects of LPS in MKI67-positive NK cells (Fig. S3). These results demonstrated the regulatory effects of TMP on the immune cell subsets in LPS-induced ALI.

Identification of molecular changes caused by TMP in epithelial cell subsets

Epithelial cells are considered as the main cell type damaged in ALI (Long, Mallampalli & Horowitz, 2022; Peters et al., 2014). Therefore, we investigated which epithelial cell subsets were affected by TMP in its protective effects against cell damage (Wang et al., 2021; Sauler et al., 2022). As shown in Fig. 7A, the epithelial cells were divided into eight sub-clusters based on marker genes. The top three genes for each cell clusters of epithelial cells are shown in a heatmap (Fig. 7B). We found that these cell subsets were generally present in different samples and groups, suggesting clustering of epithelial cells based on cell types from different mice (Figs. S5A and S5B). LPS significantly decreased the number of Krt14 high basal cells, which is a subset of cells within the C1 sub-cluster, and this decrease was reversed by TMP treatment. However, no significant differences were observed between the LPS and LPS+Dex groups (Fig. 7C). In other cell subsets, no significant cytotoxic effects were observed in the LPS group, except for the mixed lineage cell subset. Treatment with Dex led to further reduction in cell numbers, and no significant protective effects were found in the mixed lineage cell subset following the treatment with TMP when compared with the LPS group (Fig. S4). These results indicated that the treatment with TMP protected against LPS-induced cell damage in Krt14 high basal epithelial cell subset.

To examine the molecular changes in basal cells after TMP treatment, we also conducted differential gene expression analysis and performed GO, KEGG and GSEA analyses by using the differentially expressed genes to derive meaningful biological insights from our large-scale genomic dataset. GO analysis revealed that the down-regulated genes were enriched in DNA repair, intrinsic apoptotic signaling pathways in response to DNA damage and intrinsic apoptotic signaling pathways. The up-regulated genes were enriched in mitotic spindle organization, mitotic cell cycle checkpoint and cell cycle phase transition. This indicated that the treatment with TMP might regulate the cellular physiological processes related to cell survival and cell cycle progression in basal epithelial cells (Fig. S5C). KEGG analysis showed that the down-regulated genes were enriched in purine metabolism and the P53 signaling pathway, while the up-regulated genes were enriched in the ErbB, TGF-β and MAPK signaling pathways in basal cells (Fig. S5D). GSEA analysis indicated that the FOXO, ErbB and TGF-beta signaling pathways were positively associated with TMP treatment in basal cells (Fig. S5E). These results showed that the ErbB and TGF-β signaling pathways are likely to play important roles in mediating the TMP-protective effect against cell damage induced by LPS in basal epithelial cells. To identify the critical transcription factors affected by TMP in epithelial cells, we performed scenic analysis to determine the changes in their activities, especially in Krt14 high basal cells. We found that the activity of FOSL1, an essential transcription factor involved in regulating several cellular physiological processes, was increased significantly in Krt14 high basal cells (Fig. 7D). In addition, TMP treatment attenuated the LPS-induced decrease in FOSL1 expression in Krt14 high basal cells (Fig. 7E). These results suggest that FOSL1 is likely to act as a crucial downstream target in TMP protection of basal epithelial cells against LPS-induced cell death.

Identification of molecular changes caused by TMP in endothelial cell subsets

Endothelial cell dysfunction is also a key pathological change in ALI (Matthay et al., 2019; He et al., 2021). Therefore, we examined which subsets of endothelial cells could be affected by TMP (Geldhof et al., 2022; Lou et al., 2022). The endothelial cells were divided into six sub-clusters (Fig. 8A). A heatmap of the top three genes for each cell clusters of endothelial cells are shown (Fig. 8B). We found that these six sub-clusters were generally present in different samples and groups (Figs. S7A and S7B). As shown in Fig. 8C, LPS treatment decreased the number of LECs, which was a subset of cells within the C6 sub-cluster. TMP treatment significantly repressed the effects of LPS and increased the percentage of LECs, while treatment with Dex led to a decreased number of these cells. In the other cell subsets, LPS also significantly decreased the cell percentage of Atypical Capillary (aCap) and arterial endothelial cells, whereas no significant protective effects were found after treatment with TMP (Fig. S6). These results indicated that TMP could protect the LECs from damage, while Dex treatment may have negative effects on them.

Differential gene expression analysis was performed on the LECs, and GO, KEGG and GSEA analyses were conducted by using the differentially expressed genes. GO analysis showed that the down-regulated genes were enriched in the proteasomal protein catabolic process and Ras protein signal transduction, while the up-regulated genes were enriched in mitochondrial organization and response to oxidative stress in the LECs. This implied that TMP might affect the mitochondrial functions and the status of oxidative stress (Fig. S7C). KEGG analysis revealed that the down-regulated genes were enriched in protein processing in the endoplasmic reticulum and ubiquitin-mediated proteolysis, while the up-regulated genes were enriched in DNA replication, proteasome and RNA degradation (Fig. S7D). GSEA analysis indicated that biosynthesis of the amino acids as well as cysteine and methionine metabolism were positively associated with TMP treatment in the LECs (Fig. S7E). These results indicated that TMP treatment might influence protein synthesis and modification as well as amino acid metabolism. Scenic analysis revealed that JunB, a transcription factor involved in regulating oxidative stress, was activated significantly in LECs (Fig. 8D). TMP treatment significantly increased the expression of JunB in LECs, while Dex treatment did not have the same effect (Fig. 8E). These results suggested that TMP was likely to provide better protection for the survival of the endothelium in LPS-induced ALI.