Single-cell landscape of bronchoalveolar immune cells in patients with immune checkpoint inhibitor-related pneumonitis

Abstract

The pathophysiology of immune checkpoint inhibitor-related pneumonitis remains incompletely understood. We conducted single-cell and T-cell receptor transcriptomic sequencing on the bronchoalveolar lavage fluid from five patients with grade ≥2 immune checkpoint inhibitor-related pneumonitis. Our analyses revealed a prominent enrichment of T cells in the bronchoalveolar lavage fluid of patients with immune checkpoint inhibitor-related pneumonitis. Within the CD4 + T cell subset, Tfh-like T cells were highly enriched and exhibited signatures associated with inflammation and clonal expansion. Regulatory T cells were also enriched and displayed enhanced inhibitory functions. Within the CD8 + T-cell subset, effector memory/tissue-resident memory T cells with an elevated cytotoxic phenotype were highly infiltrated. Among myeloid cells, alveolar macrophages were depleted, while pro-inflammatory intermediate monocytes were elevated. Dendritic cells demonstrated enhanced antigen presentation capabilities. Cytokines CXCR4, CXCL13, TNF-α, IFN-α, IFN-γ, and TWEAK were elevated. Through a comprehensive single-cell analysis, we depicted the landscape of immune checkpoint inhibitor-related pneumonitis.

Introduction

Immune checkpoint inhibitors (ICIs), such as PD-1 and PD-L1 monoclonal antibodies, have revolutionized the treatment paradigm of various malignancies. The Food and Drug Administration (FDA) has approved multiple ICIs products for the treatment of a variety of cancer types¹. However, as the use of ICIs continues to expand, oncologists increasingly face immune-related adverse events (irAEs), which constitute a group of autoimmune-like disorders triggered by ICIs^2,3. IrAEs can affect nearly any organ or system, with widely varying onset times^4,5,6. While most irAEs are mild, a subset of patients experiences severe irAEs, such as ICI-related pneumonitis (ICP) and ICI-related myocarditis (ICM), necessitating intensive immunosuppressive therapy⁷. These life-threatening toxicities not only compromise patients’ quality of life but also impair their long-term survival.

ICP is one of the most common and lethal irAEs, accounting for 35% of anti-PD-1 monotherapy related deaths⁷. Additionally, 10–20% of patients who develop ICP do not respond to conventional corticosteroid-based therapy, posing a thorny challenge for clinicians^8,9. To date, there remains a lack of effective treatments based on the pathophysiological mechanisms of ICP.

Previous reports have suggested the involvement of T lymphocytes in the pathogenesis of ICP. Suresh et al. reported an enrichment of CD4+central memory T cells in the bronchoalveolar fluid (BALF) of ICP patients, with a Th1-skewed immune response characterized by IFN-γ and TNF-α secretion¹⁰. Wang et al. profiled the cellular and cytokine landscape within the alveolar compartment of ICP patients and observed an increase in Th17 cells and IL-17A in BALF, with Th17 cells known for their secretion of the pro-inflammatory cytokine IL-17, which has been implicated in autoimmune disorders such as psoriasis^11,12,13. PD-1/PD-L1 axis blockade has been shown to induce a Th2-Th1 phenotype shift in cancer patients¹⁴. Overall, a predominance of Th1/Th17-mediated inflammation likely contribute to the development of ICP. Suresh’s study also suggested a potential role of regulatory T cells (Tregs) in the onset of ICP, as Tregs in ICP BALF expressed fewer co-inhibitory molecules such as CTLA-4 and PD-1, indicating reduced inhibiting function¹⁰.

Despite these fundamental studies, little is known about specific T-cell subsets and other immune cells that may play a role in the pathogenesis of ICP. In our study, we employed single-cell RNA sequencing (scRNA-seq) with high resolution to profile the immune cell landscape of ICP and elucidate the molecular pathways involved in its pathophysiology. This research aimed to deepen our understanding of ICP’s mechanisms and provide insights for the development of novel diagnostic and therapeutic approached to address this clinically significant adverse event.

Methods

Patients and data sources

We prospectively enrolled patients who had received ICIs and subsequently developed grade≥2 ICP before the initiation of immunosuppressive therapy from May 1st, 2022, to December 31st,2022. The inclusion criteria for enrolled patients were as follows: 1) age≥18 years; 2) histologically diagnosed with a single solid tumor (excluding primary lung cancer); 3) received at least one dose of anti-PD-1/PD-L1 treatment; 4) met the physiological requirements for bronchoscopy under general anesthesia; 5) ICP diagnosis was confirmed by two experienced physicians based on clinical, laboratory and radiographic findings and graded according to the Common Terminology Criteria for Adverse Events version 5.0 (CTCAE v5.0). Exclusive criteria included: 1) a diagnosis of other respiratory diseases within the past 3 months, such as cancer lung metastasis, asthma, chronic obstructive pulmonary disease, interstitial pulmonary disease, and pleural effusion, as determined by imaging examination, pulmonary function testing, allergen testing, etc.; 2) a history of lung radiotherapy or exposure to drugs with pulmonary toxicity (e.g., bleomycin, antibody-conjugated drugs) within the past 3 months; 3) diagnosis of hematologic or rheumatic diseases (e.g., thrombocytopenia, psoriasis, vasculitis) within the past 3 months, irrespective of medication requirements; 4) local or systemic application of glucocorticoids or immunosuppressants within the past 2 weeks; 5) suspected comorbid pulmonary infection. We also used single-cell gene expression data of BALF from four healthy subjects and nine patients diagnosed with coronavirus disease-2019 (Covid-19) pneumonia as control. These datasets were publicly accessible at https://www.ncbi.nlm.nih.gov/gds/?term=GSE145926.

We also established a retrospective patient cohort from our institute for cytokine profiling. Stored BALF samples from ICP patients, COVID-19 pneumonia patients, and healthy individuals were retrieved. Eligible patients were those diagnosed with any solid tumors, who received anti-PD-1/PD-L1 therapies, and who developed any grade of ICP.

Ethics approval and informed consent

The study was conducted in accordance with the Declaration of Helsinki and received approval from the ethics committee of Peking University Cancer Hospital (Approval number: 2022KT125). Patients signed written informed consent for the bronchoscopy performed and for the anonymous use of their clinical and molecular data for research purposes. Written informed consent was obtained from patients for the publication of any data included in this article.

BALF sample preparation

All enrolled patients in the prospective cohort underwent standard bronchoscopy under general anesthesia to collect BALF samples. Generally, the bronchoscope was guided to the ICP-affected lung lobe or segment, and preheated sterile saline was injected into the airway multiple times. A total of 30 ml of lavage fluid was retrieved from each patient, with 5–10 ml reserved for routine clinical tests and the remainder used for scRNA-seq. The collected BALF samples were encapsulated using 20 ml centrifuge tubes and immediately placed on ice.

The cooled BALF samples were mixed with RPMI 1640 complete medium and filtered through a 100 mm nylon cell strainer to obtain single-cell suspensions. To minimize the potential for blood cell contamination, red blood cells were removed using Red Blood Cell Lysis Buffer. After two washes with 1x PBS (Invitrogen), the cell pellets were resuspended in sorting buffer (PBS supplemented with 1% FBS). Based on fluorescence-activated cell sorter (FACS) analysis, single CD45+cells from each patient’s samples were sorted using anti-human CD45 antibodies (Biolegend) into 1.5 ml low binding tubes (Eppendorf) with 50 ml of sorting buffer for scRNA-seq.

Single-cell library preparation and sequencing

Cell suspensions were barcoded through the 10x Chromium Single Cell platform using Chromium Next GEM Single Cell 5’ Kit v2. Single cells were washed once with PBS containing 0.04% bovine serum albumin (BSA) and resuspended in PBS containing 0.04% BSA to a final concentration of 500–1200 cells/ml as determined by Rigel S2 cell counter (Countstar). The loaded cell numbers ranged from 300–500000 aiming for 300–14000 single cells per reaction. 5000 cells were captured in droplets to generate nanoliter-scale Gel bead in EMulsion (GEMs). GEMs were then reverse transcribed in a C1000 Touch Thermal Cycler (BioRad) programmed at 53 °C for 45 min, 85 °C for 5 min, and held at 4 °C. After reverse transcription and cell barcoding, emulsions were broken and cDNA was isolated and purified with Cleanup Mix containing DynaBeads and SPRIselect reagent (ThermoFisher Scientific), followed by PCR amplification. Amplified cDNA was then used for 5’ gene expression library construction. For RNA-seq library construction, amplified cDNA was fragmented and end-repaired, double-sided size-selected, PCR-amplified with sample indexing primers, and repeated double-sided size-selection. For T-cell receptor (TCR) library construction, TCR transcripts were enriched from amplified cDNA by PCR. Subsequently, enriched PCR product was fragmented and end-repaired, size-selected, PCR-amplified with sample-indexing primers, and size-selected. Libraries prepared according to the manufacturer’s instructions were then purified and profiled for quality assessment. Single-cell RNA and TCR V(D)J libraries were sequenced by an Illumina NovaSeq 6000 sequencer with 150 bp paired-end reads.

Cytokine measurement by cytometric bead array

Cytokines in the BALF, including CXCR4, CXCL13, IL-17A, TWEAK, IFN-γ, IL-6, IFN-α, and TNF-α, were detected using Luminex200 (Luminex Corporation). Briefly, BALF samples were thawed, centrifuged at 1000 g at 4 °C for 15 min, and the supernatant was collected. Antibodies were coated onto the surface of magnetic microbeads to form a solid-phase carrier for cytokine binding. Blended microbeads (10 μL), cytokine standards (100 μL), and supernatant samples (100 μL) were added to the microwells. Biotinylated antibodies were then added, and after washing away the unbound biotinylated antibodies, PE-labeled streptavidin was introduced. Fluorescence excitation allowed the median fluorescence intensity (MFI) to quantify the concentration of different cytokines using a Luminex analyzer. Each sample was measured in duplicate, and the average value was taken as the final concentration.

Bioinformatic analysis

Raw sequencing data were aligned and quantified with the Cell Ranger Software v6.0 (10×Genomics) using the human reference genome GRCh38 v5.0.0. Seurat v4.1.0 software was used for initial data integration and analysis. All filtered data were required to meet the following criteria: 1) a median of greater than 1000 genes per cell; 2) sequencing saturation falling within the range of 60% to 90%; 3) a minimum fraction of reads in cells exceeding 70%.

The original expression matrix and the single-cell expression matrix from the GSE145926 dataset, including patients with Covid-19 pneumonia and healthy controls, were integrated into Scanpy python toolkit v1.7.2. Genes with less than 3 expressed cells and cells with less than 100 expressed genes were filtered out. The calculate_qc_metrics function was used to calculate the proportion of mitochondrial genes in each cell:

$${{Mt}}_{i}=\frac{{mitochondria\; Gene\; Counts}}{{Total\; Gene\; Counts}}$$

(1)

Cells were further screened based on the following criteria:

• 1. the proportion of mitochondrial genes should be less than 25%.

• 2. ${g}_{i}$ represents the number of genes detected in the i^th cell, and N represents any cell in the sample.

$${G}_{i}={log}_{10}{g}_{i}$$

(2)

$${G}_{i}\, > \,{\rm{Median}}\left(\left\{{G}_{i}|i\in N\right\}\right)-3\times {\rm{MAD}}\left(\left\{{G}_{i}|i\in N\right\}\right)$$

(3)

The scrublet package was used to process the single-cell expression matrix and calculate the Doublet score for each cell. To minimize the potential bias, we over-clustered the cells and calculated the average double cell score for each cell cluster. Cell clusters with a double cell score greater than 0.6 or cell clusters expressing more than 1 canonical marker genes were filtered out.

Dimension reduction and clustering were performed using the Scanpy python toolkit v1.7.2. Data normalization, logarithmic transformation, variable genes calling, scaling and principal component analysis (PCA) were performed. To mitigate batch effects, data were processed using the batch-balanced k-nearest neighbors (BBKNN) method. Uniform manifold approximation and projection for dimension reduction (UMAP) embedding plots, t-distributed stochastic neighbor embedding (t-SNE), and diffusion maps were used for data visualization. Main cell types were annotated based on established canonical marker genes. Unsupervised clustering of different cell subtypes was manipulated based on logistic regression and machine learning, utilizing the CellTypist package, with reference to the Human lung atlas, Immune allow, and Immune all high cell sets¹⁵. The annotation results were further proofread by canonical marker genes.

After data normalization, the relative proportions of each cell subtype across different populations were calculated. To account for variations in cell numbers between different samples, the normalized relative cell number of the j^th sample was calculated first:

$${C}_{k,j}^{{norm}}=\frac{{C}_{k,j}}{\sum _{k\in R}{C}_{k,j}}$$

(4)

$R$ represents the set of all cell subtypes, and ${C}_{k,j}$ represents the cell number in the k^th cell subtype of the j^th sample.

The relative numbers of cell subtypes in the BALF of patients with ICP and healthy controls, ${C}_{k,P}^{{norm}}$ and ${C}_{k,L}^{{norm}}$, were calculated:

$${C}_{k,P}^{{norm}}=\sum _{j\in P}{C}_{k,j}^{{norm}}$$

(5)

$${C}_{k,L}^{{norm}}=\sum _{j\in L}{C}_{k,j}^{{norm}}$$

(6)

$P$ represents patients with ICP, and $L$ represents health controls.

Finally, the relative proportions of different cell subtypes were obtained:

$${P}_{k,P}=\frac{{C}_{k,P}^{{norm}}}{{C}_{k,P}^{{norm}}+{C}_{k,L}^{{norm}}}\times 100 \%$$

(7)

$${P}_{k,L}=\frac{{C}_{k,L}^{{norm}}}{{C}_{k,P}^{{norm}}+{C}_{k,L}^{{norm}}}\times 100 \%$$

(8)

${P}_{k,P}$ represents the relative proportion of cell subtype $k$ in the BALF of patients with ICP, and ${P}_{k,L}$ represents the relative proportion of cell subtype $k$ in the BALF of healthy controls.

Genes with an adjust p-value (P_adj) < 0.05 and |log₂ (fold change) | > 0.5 were selected as differential expressed genes (DEGs). Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) and gene set enrichment analysis (GSEA) were performed using the clusterProfiler v3.8.1 and the gsea v3.0 in the R package.

The Cellchat package was used to manipulate cell communication analysis. Representative cell types were selected by randomly sampling 500 cells for each cell subtype. The sampled cells were analyzed for their cell communication intensity based on the expression of specific ligands and receptors. Cell communications that existed in more than five cells were retained to minimize the interference effect of low-abundance cell communications.

The clipped and unclipped expression matrices were obtained using the velcyto tool. Cell velocity was calculated through the scvelo tool, and the PAGA function of the Scanpy package was employed to generate directed developmental pathway diagrams for different cell subtypes.

The cytokine activation states of each cell subtype were estimated using the CytoSig software, a penalized ridge regression-based cytokine signaling analyzer¹⁶.

The STARTRAC algorithm, based on Shannon entropy, was used to characterize CD4 + T cell developmental state transitions, considering transcriptional similarity between cells and cell velocity¹⁷.

Statistics

Continuous variables were compared using the Wilcoxon rank sum test. The Benjamini & Hochberg method was used for p-value adjustment when making multiple comparisons. A two-sided alternative hypothesis at the 5% significance level was adopted. All data analyses, statistical analyses and graphing were performed using Python v3.6.10, R v3.6.3, and GraphPad PRISM v9.

Results

Characterization of the immune landscape in ICP patients

To comprehensively characterize the immune landscape of the BALF in patients with ICP, we recruited five patients who received ICIs and developed grade≥2 ICP before receiving any immunosuppressive therapy. Table 1 presents the demographic and clinical characteristics of these patients. After collecting fresh BALF, CD45+immune cells were isolated through FACS. In total, we obtained 58022 high-quality single cells, which were sequenced by 10x Genomics 5’ mRNA and TCR sequencing kits (Fig. 1A). To contrast the differences between the BALF in ICP patients and that in healthy controls, we integrated our data with that of a published study containing three healthy subjects and nine patients diagnosed with Covid-19 pneumonia¹⁸. We used the BBKNN integration method to merge cells generated from different patients or sources, and subsequently annotated cells with CellTypist.

Table 1 Demographic and clinical characteristics of patients with ICP

Fig. 1: The immune landscape of the BALF of patients with ICP.

A Diagram of the workflow for single-cell sequencing and UMAP embedding plot displaying immune cells generated from the BALF of patients with ICP, Covid-19 pneumonia, and healthy controls. B UMAP embedding plots of immune cells colored according to the expression of canonical marker genes. C, D UMAP embedding plots of immune cells grouped into 8 major types and 23 cell subtypes. E Bar plots indicating the differences in the main cell types among different patient groups. ICP, n = 5; Covid-19, n = 8; Healthy, n = 3. Statistical testing was performed by a Kruskal-Walis test. Data were presented as mean values (central line) +/- interquartile range (border). (Created with BioRender.com).

Subsequently, we classified all isolated immune cells into eight major classes: B cells (CD79A), CD4 + T cells (CD3D, CD4), CD8 + T cells (CD3D, CD8A), NK cells (GNLY), macrophages (LYZ), monocytes (NAMPT, IFITM2), and dendritic cells (CD74) (Fig. 1B, C). These major classes were further subdivided into 23 subclasses (Fig. 1D). The detailed information on the annotation genes for all cell types and subtypes were shown in Supplementary data 1-3. To investigate the abundance of different immune cells in the BALF from different patient sources, we calculated the relative proportion of patients within each major cell class. Notably, CD4 + T cells and CD8 + T cells were significantly enriched in the BALF of patients with ICP compared to that in healthy controls or patients with Covid-19 pneumonia (Fig. 1E). Meanwhile, a relative reduction of macrophages was shown in the BALF of patients with ICP and patients with Covid-19 pneumonia (Fig. 1E). The quantities of other major cell classes did not significantly differ across the different populations.

Enrichment of Tfh-like T cells in the BALF of patients with ICP

To gain insight into the role of CD4 + T cells in the pathogenesis of ICP, we subclustered CD4 + T cells and annotated them into five subtypes base on their specific gene expression features (Fig. 2A–C). Follicular helper T (Tfh) cells are a subset of CD4 + T cells characterized by high expression of CXCR5, CXCL13, and immune modulatory molecules such as PD-1 and ICOS. Within the CD4 + T cell population, we identified a subset that expressed both Th1 signature genes (IFNG, CXCR3) and Tfh signature genes (CXCL13), which we defined as Tfh-like T cells. These Tfh-like T cells constituted a significantly higher percentage of CD4 + T cells in the BALF of patients with ICP compared to other groups of patients (Fig. 2D). Additionally, T effector memory (Tem)/effector helper T cells, characterized by pro-inflammatory signatures (RPL41, RPS27, ANXA1, CXCR4), were enriched in the BALF of both patients with ICP and patients with Covid-19 pneumonia (Fig. 2D). Tregs and exhausted T (Tex) helper cells showed a higher percentage in the BALF of patients with ICP than in other group of patients (Fig. 2D).

Fig. 2: Phenotypes, quantitative and functional analyses of CD4 + T cells.

A, B UMAP embedding plots of CD4 + T cells colored according to patient groups and cell labels. C Dot plots displaying canonical marker genes of CD4 + T cell subclusters. The sizes of the dots are proportional to the fraction of the cells, and the color represents the normalized gene expression level. D Bar plots indicating the differences in the cell subtypes within the total cell population and within the CD4 + T cell populations among different patient groups. ICP, n = 5; Covid-19, n = 8; Healthy, n = 3. Statistical testing was performed by a Kruskal-Walis test. Data were presented as mean values (central line) +/- interquartile range (border). E UMAP embedding plots displaying the RNA velocities of CD4 + T cells. Cells were colored according to their cell labels and the clonal expansion ability. F Histogram showing the relative abundance of CD4 + T cells with different clonal expansion ability. Cells were grouped according to patient groups and cell subtypes. G Box plots indicating clonal expansion (upper) and developmental transition (lower) of CD4 + T cells quantified by STARTRAC for each sample divided by cell subtypes or patient groups. Two-sided Wilcoxon tests were performed to make comparisons. The center lines correspond to the median values, and the boxes correspond to the interquartile range. H Heatmap of cytokine signaling activities in CD4 + T cells among patients with ICP, patients with Covid-19 pneumonia, and healthy controls.

To further understand the expansion and transition of CD4 + T cells, we conducted single-cell TCR sequencing. We observed stronger clonal expansion of CD4 + T cells in the BALF of patients with ICP compared to patients with Covid-19 pneumonia (Fig. 2F). Among CD4 + T cell subtypes, Tfh-like T cells exhibited the strongest clonality and expansion ability (Fig. 2G). Additionally, the STARTRAC results indicated enhanced transition and expansion ability in most CD4 + T cells in the BALF of ICP patients, with Tfh-like T cells showing the most robust transition ability (Fig. 2G).

We also examined the cytokine signaling pattern in our integrated datasets. Anti-virus responses-related cytokines, including IL-6, IFN-1, and IFN-G, were upregulated in the BALF of patients with Covid-19 pneumonia (Fig. 2H). In contrast, anti-inflammatory cytokines such as BMP6 and TGFB1 were upregulated in the BALF of healthy controls (Fig. 2H). Interestingly, anti-inflammatory cytokine IL-10 and Activin A, a member of the TGF-β family, were upregulated in the BALF of patients with ICP (Fig. 2H). Notably, tumor necrosis factor (TNF)-α, IL-17, and TNF-like weak inducer of apoptosis (TWEAK), a member of the TNF family, were significantly elevated in the BALF of patients with ICP (Fig. 2H).

Activation of conventional CD4 + T cells and Tregs in the BALF of patients with ICP

Further investigation of Tregs in ICP revealed that they were more activated in the BALF of ICP patients compared to healthy controls (Fig. 3A–D). In contrast, Tregs in the BALF of patients with Covid-19 pneumonia highly expressed IFN-related genes, such as IFIT1 and IFI27 (Fig. 3A, B).

Fig. 3: Phenotypes, quantitative and functional analyses of regulatory and conventional T cell.

A Volcano plots displaying the DEGs of Tregs between patients with ICP, patients with Covid-19 pneumonia, and healthy controls. Statistical tests were performed by a two-sided Wilcoxon test, with adjusted P value < 0.05 and log₂ (fold change) ≥ 0.5. B Lollipop plots showing the KEGG and GO biological process enrichment results for DEGs in Tregs. C GSEA analysis indicated that genes ranked in Tregs from patients with ICP and healthy controls were enriched in the pathway “cytokine receptor interaction”. Statistical analysis was performed by permutation test. D Heatmap of core enrichment genes in the pathway “cytokine-cytokine receptor interaction” from the GSEA analysis of Tregs among patients with ICP, patients with Covid-19 pneumonia, and healthy controls. E Volcano plots displaying the DEGs of conventional T cells between patients with ICP, patients with Covid-19 pneumonia, and healthy controls. Statistical tests were performed by a two-sided Wilcoxon test, with adjusted P value < 0.05 and log₂ (fold change) ≥ 0.5. F Lollipop plots showing the KEGG and GO biological process enrichment results for DEGs in conventional T cells. G GSEA analysis indicated that genes ranked in conventional T cells from patients with ICP and healthy controls were enriched in the pathway “leukocyte trans-endothelial migration”. Statistical analysis was performed by permutation test. H Heatmap of core enrichment genes in the pathway “leukocyte trans-endothelial migration” from the GSEA analysis of conventional T cells among patients with ICP, patients with Covid-19 pneumonia, and healthy controls.

DEGs and enrichment analyses of conventional CD4 + T cells in the BALF of patients with ICP indicated a feature related to Th17 cell differentiation and stronger trans-endothelial migration ability compared to healthy controls (Fig. 3E–H). Conventional CD4 + T cells in the BALF of patients with Covid-19 pneumonia indicated a higher cytosolic DNA-sensing ability instead, highlighting their specific roles in anti-virus responses (Fig. 3F).

Overactivated CD8 + T cells in the BALF of patients with ICP

We analyzed the gene expression signatures of CD8 + T cells, including T effective memory (Tem)/CD45RA re-expressed terminally differentiated effective memory (Temra) cytotoxic T cells, and Tem/tissue-resident memory (Trm) cytotoxic T cells (Fig. 4A, B). The primary compartment of CD8 + T cells in the BALF of patients with ICP was Tem/Trm cytotoxic T cells, which constituted a higher percentage in the BALF of patients with ICP than in other patient groups (Fig. 4C). The expansion ability of Tem/Trm cytotoxic T cells was robust in the BALF of patients with ICP (Fig. 4D, E).

Fig. 4: Phenotypes, quantitative and functional analyses of CD8 + T cells.

A, B UMAP embedding plots of CD8 + T cells colored according to patient groups and cell labels. C Bar plots indicating the differences in the cell subtypes within the total cell population and within the CD8 + T cell populations among different patient groups. ICP, n = 5; Covid-19, n = 8; Healthy, n = 3. Statistical testing was performed by a Kruskal-Walis test. Data were presented as mean values (central line) +/- interquartile range (border). D Histogram showing the relative abundance of CD8 + T cells with different clonal expansion ability. Cells were grouped according to patient groups. E Clonal expansion (left) and developmental transition (right) of CD8 + T cells quantified by STARTRAC for each sample divided by patient groups. Two-sided Wilcoxon tests were performed to make comparisons. The center lines correspond to the median values, and the boxes correspond to the interquartile range. F Volcano plots displaying the DEGs of CD8 + T cells between patients with ICP, patients with Covid-19 pneumonia, and healthy controls. Statistical tests were performed by a two-sided Wilcoxon test, with adjusted P value < 0.05 and log₂ (fold change) ≥ 0.5. G Lollipop plots showing the KEGG and GO biological process enrichment results for DEGs in CD8 + T cells. H GSEA analysis indicated that genes ranked in CD8 + T cells from patients with ICP and healthy controls were enriched in the pathway “regulation of leukocyte migration”. Statistical analysis was performed by permutation test. I Heatmap of core enrichment genes in the pathway “regulation of leukocyte migration” from the GSEA analysis of CD8 + T cells among patients with ICP, patients with Covid-19 pneumonia, and healthy controls. J GSEA analysis indicated that genes ranked in CD8 + T cells from patients with ICP and healthy controls were enriched in the pathway “response to ketone”. Statistical analysis was performed by permutation test. K Heatmap of core enrichment genes in the pathway “response to ketone” from the GSEA analysis of CD8 + T cells among patients with ICP, patients with Covid-19 pneumonia, and healthy controls. L Heatmap of cytokine signaling activities in CD8 + T cells among patients with ICP, patients with Covid-19 pneumonia, and healthy controls.

DEGs analysis of CD8 + T cells reflected distinct pathogenic mechanisms between cytotoxic T cells in ICP and Covid-19 pneumonia. In the BALF of patients with ICP, cytotoxic T cells exhibited a more vigorous migration ability and upregulated critical metabolic pathways associated with autoimmune disease activity (Fig. 4F–J). Interestingly, genes stimulated by ketones were also upregulated in CD8 + T cells of patients with ICP (Fig. 4J). In contrast, anti-virus signatures were profoundly expressed in CD8 + T cells of patients with Covid-19 pneumonia (Fig. 4F–G).

We also profiled the cytokine signaling pattern of CD8 + T cells. Anti-inflammatory cytokines, including BMP2, BMP6, and TGFB1, were upregulated in CD8 + T cells in the BALF of healthy controls (Fig. 4L). Certain pro-inflammatory cytokines, including TWEAK, were upregulated in CD8 + T cells in the BALF of patients with ICP (Fig. 4L). Interestingly, the anti-inflammatory cytokine IL10 was also upregulated in CD8 + T cells in patients with ICP (Fig. 4L).

Elevation of antigen presentation ability of dendritic cells in the BALF of patients with ICP

Myeloid cells were annotated and subdivided into seven subtypes based on their canonical marker gene expressions. Macrophages and monocytes were classified into three subtypes: classical monocytes (IFITM2, IFITM3), intermediate macrophages (CCL2), and alveolar macrophages (C1QA, APOC1) (Fig. 5A–D). Dendritic cells (DCs) were divided into four subtypes: plasma DC (GZMB, TCF4), DC1 (LGALS2), DC2 (CDC1), and migratory DCs (LAMP3) (Fig. 5A–D). Intermediate macrophages and classical monocytes were scarce in the BALF of healthy populations, and almost all macrophages in the BALF of healthy individuals were tissue-resident alveolar macrophages (Fig. 5E). The enrichment analysis of DEGs indicated an enhanced migration capability of macrophages in the BALF of patients with ICP (Fig. 5F). The cytokine signaling pattern of macrophages in patients with ICP was similar to that of CD8 + T cells (Fig. 5G).

Fig. 5: Phenotypes, quantitative and functional analyses of myeloid cells.

A–C UMAP embedding plots of myeloid cells colored according to patient groups, main labels, and cell labels. D Dot plots displaying canonical marker genes of myeloid cell subclusters. The sizes of the dots are proportional to the fraction of the cells, and the color represents the normalized gene expression level. E Bar plots indicating the differences in the cell subtypes of myeloid cells within the total cell population among different patient groups. ICP, n = 5; Covid-19, n = 8; Healthy, n = 3. Statistical testing was performed by a Kruskal-Walis test. Data were presented as mean values (central line) +/- interquartile range (border). F Lollipop plots showing the GO biological process enrichment results for DEGs in macrophages. G Heatmap of cytokine signaling activities in myeloid cells among patients with ICP, patients with Covid-19 pneumonia, and healthy controls.

Notably, the percentage of DC2 was significantly higher in the BALF of patients with ICP compared to patients with Covid-19 pneumonia, indicating its potential role in the pathogenesis of ICP (Fig. 6A). In the BALF of patients with ICP, both DC1 and DC2 highly expressed major histocompatibility complex (MHC) II genes, showing an augmented antigen presentation ability in ICP (Fig. 6B, C). Cell communication analysis suggested an enhancement in the MHC-II signaling between conventional dendritic cells (cDC1 and cDC2) and CD4 + T cells in the BALF of patients with ICP, illustrating an augmented interaction between antigen-presenting cells and helper T cells, which may contribute to the autoimmune response in ICP (Fig. 6D, E).

Fig. 6: Phenotypes, quantitative and functional analyses of DCs.

A Bar plots indicating the differences in the cell subtypes of DCs within the total cell population among different patient groups. ICP, n = 5; Covid-19, n = 8; Healthy, n = 3. Statistical testing was performed by a Kruskal-Walis test. Data were presented as mean values (central line) +/- interquartile range (border). B Heatmap of core enrichment genes in the pathway “autoimmune thyroid disease from the GSEA analysis of DC1 among patients with ICP, patients with Covid-19 pneumonia, and healthy controls. C Heatmap of core enrichment genes in the pathway “antigen processing and presentation” from the GSEA analysis of DC2 among patients with ICP, patients with Covid-19 pneumonia, and healthy controls. D Cell communication analysis suggested an augmented antigen presentation ability of DCs in patients with ICP. E Cell communication analysis suggested a strong signaling interaction between DCs and cytotoxic T cells or Tfh-like T cells.

Distinct cytokine profile characterizes the BALF of ICP

In the validation cohort, a total of 28 individuals with ICP, 10 patients with COVID-19 pneumonia, and 6 healthy individuals were enrolled. The clinical characteristics of these patients and healthy subjects are detailed in Table 2.

Table 2 Clinical characteristics of participants in the validation cohort

The cytometric bead array results indicated a marked elevation in pro-inflammatory cytokines within the BALF of ICP patients. Specifically, the levels of CXCR4, CXCL13, TNF-α, IFN-α, IFN-γ, and TWEAK were significantly elevated, corroborating the presence of an overactive immune response in ICP. These findings indirectly support the notion that Tfh-like T cells, hyperactivated CD8 + T cells, and effector memory T cells are enriched in the BALF of patients with ICP (Fig. 7).

Fig. 7: Distinct cytokine profile characterizes the BALF of ICP.

Cytokines CXCR4, CXCL13, TNF-α, IFN-α, IFN-γ, and TWEAK were significantly higher in the BALF of patients with ICP than in healthy controls. Statistical testing was performed by unpaired t-test. Data were presented as mean values (central line) +/- interquartile range (quartile line) in the violin plots. ^**P value < 0.01; ^****P value < 0.0001; COV, covid-19 pneumonia; HC, healthy controls.

Discussion

In this study, we conducted comprehensive single-cell transcriptomic analyses to profile the immune cell atlas of BALF from patients with ICP. Our findings are of paramount importance as they illuminate the underlying mechanisms of irAEs and guide the identification of effective immunomodulatory treatments for ICP. Our results reveal the immune microenvironment in ICP characterized by CD8+Tem/Trm cytotoxic T cells and CD4+ Tfh-like T cells with unique pro-inflammatory gene signatures. Myeloid cells also contribute significantly to autoimmune responses by enhancing antigen presentation to T cells.

We meticulously delineated the immune cell landscape of BALF in ICP and compared it among three populations: patients with ICP, patients with Covid-19 pneumonia, and healthy individuals. This comparison was crucial because ICP and Covid-19 pneumonia share certain similarities, such as acute pneumonitis, interstitial lung involvement, and overactivation of the adaptive immune system. However, our results indicate significant differences in gene expression patterns among various immune cells between ICP and Covid-19 pneumonia.

Our analysis revealed that CD8+ and CD4 + T lymphocytes accounted for 54% of all CD45+immune cells in the BALF of patients with ICP, compared to only 10% in healthy controls. This suggests that T cell abnormalities play a dominant role in the pathogenesis of ICP. Our finding s build upon earlier studies that also reported prominent T cell infiltration in ICP, while they did not characterize the precise phenotype and transition state of these T cells^19,20.

Building on this background, our study identified specific cell subtypes associated with ICP, such as CD8+Tem/Trm cytotoxic T cells and CD4+Tfh-like T cells. In the BALF of patients with ICP, CD8+cytotoxic T cells highly expressed cytotoxic genes such as GZMA, GZMB, and GZMK, while conventional CD4 + T cells displayed a pro-Th1 inflammatory phenotype, marked by high expression of TNFA, IFNG, IL1A, and IL17A, indicating a type 1 inflammatory response. This is distinct from the antiviral pathway observed in the BALF of patients with Covid-19 pneumonia, marked by high expression of interferon-related genes. Tfh-like T cells, which are highly enriched in the BALF of patients with ICP, exhibited both an inflammatory phenotype and strong antigen-driven clonal proliferation. These cells are characterized by high expression of inhibitory molecules such as PD-1 and play a crucial role in the maturation of both CD8 + T cells and B cells. Consequently, it is possible that ICIs could promote a dysregulated Tfh-like cell response, contributing to a proinflammatory milieu. This type of cells has also been recognized to be associated with synovitis in anti-PD-1 treatment-related arthritis²¹.

A cross-sectional study of soluble cytokine levels in the BALF of patients with different lung diseases corroborated our findings by showing significantly elevated IFN-γ in ICP patients, distinct from other pulmonary diseases¹⁹.

Notably, both CD8+ and CD4 + T cells in ICP highly expressed TWEAK, a soluble cytokine secreted by various cell types including structural and immune cells. It pays a critical role in regulating multiple cellular activities such as proliferation, migration, differentiation, and apoptosis. The interaction between TWEAK and its ligand Fn14 can activate the downstream transcription factor NF-κB, leading to tissue inflammation and damage^22,23,24. Overactivation of the TWEAK/Fn14 pathway has been reported in various autoimmune diseases, such as psoriasis²³. Gupta, et al. demonstrated that TWEAK acts synergistically with TNF and IL-17A, promoting feedback inflammation in psoriasis²². Our cytokine analysis of BALF also confirmed an elevation of TWEAK in the BALF of ICP. We speculate that the TWEAK/Fn14 signaling pathway has the potential to serve as a hallmark of tissue damage and may be a promising therapeutic target for ICP. However, further investigations are required to determine whether blocking this pathway could prevent or alleviate ICP.

Additionally, we observed an increase in the quantities and inhibitory function of FOXP3+Treg cells (with highly expressed FOXP3, CTLA4, and TGFB) in the BALF of patients with ICP, which was not observed in patients with Covid-19 pneumonia. Activin A, a member of the TGF-β cytokine superfamily, is associated with the generation of FOXP3+Tregs and airway protective remodeling^25,26. Our CytoSig analysis indicated that this cytokine, along with IL10, was highly secreted by CD4 + T cells in ICP, aligning with the enhanced immunosuppressive function of Tregs in ICP. These findings suggest that Tregs play a role in neutralizing and restraining inflammatory responses triggered by cytotoxic and helper T cells in ICP. Given the high abundance of Tregs in the immune microenvironment of ICP, it’s reasonable to assume that preclinical models of irAEs based on Tregs depletion, such as conditional FOXP3 gene knockout mice, may not accurately simulate the real pathophysiology of ICP. Previous literature indicates that Tregs may play different roles in the context of various types of irAEs. While reductions in tissue-resident Tregs are observed in endocrine system-affected irAEs, Tregs typically undergo reactive expansion to counteract the existing immune response in ICIs-related enteritis, ICP, and arthritis^21,27,28,29.

Our study also yielded insights into myeloid cells in the BALF of ICP. Alveolar-resident macrophages (highly express C1QA, APOC1, APOE) are relatively depleted while pro-inflammatory monocyte-derived macrophages (highly express CCL2, IFITM3, ISG15) were elevated in the BALF of patients with ICP compared to healthy controls. Another group has also elucidated a profound interaction between intermediate monocytes and helper T cells via GM-CSF signaling in ICP, although our data did not shed light on this communication signal³⁰. These findings collectively suggest that monocytes might migrate into the BALF and differentiate into intermediate macrophages with more pro-inflammatory signatures. Additionally, for the first time, we observed an enhanced MHC II signaling in DCs in the BALF of patients with ICP, indicating that DCs might participate in self-antigen presentation to pathogenic T cells and mediate their activations. However, in vitro assays are needed to verify these causal hypotheses.

ICP is one of the most frequent fatal irAEs, accounting for one third of all ICIs-related deaths⁷. While most cases of ICP can be well managed with steroid-based therapies, these treatments come with long-term side effects and can adversely affect tumor control. Moreover, 20–30% of ICP cases are resistant to steroid-based therapies, and there is still a lack of precise treatment approaches based on clear scientific rationale³¹. Our study offers a potential mechanistic understanding of ICP and supports early-phase clinical trials investigating novel treatments targeting Tem/Trm cytotoxic T cells, Tfh-like T cells, and inflammatory molecules such as TWEAK. Our profiling data are insufficient to determine whether targeting these cells or molecules could eliminate irAEs while preserving antitumor efficacy. Therefore, controlled clinical trials with endpoints addressing both efficacy and toxicity are needed.

The strength of this study lies in its stringent eligibility criteria, as patients with any other lung diseases were excluded to avoid potential confounding effects on the immune microenvironment of the BALF. Consequently, all enrolled patients had a definite diagnosis of high-grade ICP mediated by anti-PD-1 agents. However, this study has several limitations. First, the small sample size and the absence of in-house data for healthy controls are two inherent defects. Despite batch correction and the application of algorithms to optimize the signal-to-noise ratio, the data still have limitations in their representativeness. Third, it is possible that certain pathophysiological processes are unique to patients who receive combination immunotherapy. Therefore, future investigations should take account into different grades of ICP, different treatment backgrounds, and integrating simultaneous changes in PBMCs during ICP development to achieve a more comprehensive understanding of the condition.

In summary, we have shed light on the mechanisms underlying ICP through comprehensive single-cell and TCR transcriptional sequencing. Our findings highlight the pivotal roles played by various immune cells, including Tfh-like T cells, Tem/Trm cytotoxic T cells, and DCs, in mediating the autoimmune responses that characterize ICP. Notably, we have identified the significance of inflammatory molecules such as TWEAK in the pathogenesis of ICP, which presents a potential therapeutic target for this condition. These insights offer valuable rationales for the development of innovative therapeutic approaches aimed at effectively managing and treating ICP.