Impact of chemotherapy and immunotherapy on the composition and function of immune cells in COVID-19 convalescent with gynecological tumors

Single cell landscape of PBMCs of patients with gynecologic tumor and COVID-19

By performing 10X Genomics scRNA-seq, we studied the transcriptomic profiles of each PBMCs from 4 gynecological tumor patients recovering from COVID-19 and undergoing cancer treatment (TCs), 4 COVID-19 convalescent (Covs) patients without cancer and 3 healthy donors (HDs) with no cancer history (Figure 1A and Supplementary Table 1). Using the unified single-cell transcriptomic analysis pipeline (see Methods), totally, 111,798 cells across all HDs, Covs and TCs, were integrated into an unbatched dataset. On average, there are 10,163 cells per sample.

Figure 1. Differences in cell compositions by single-cell transcriptomes of PBMCs. (A) Schematic of the study design. (B) t-SNE plot of cells colored by cell types identified from HDs (n = 3), Covs (n = 4), TCs before treatment (n = 4). (C) Expression of selected canonical markers using violin plots in each cell type. (D) Proportion of each cell type in each sample. Bars are colored by cell types. (E) Box plots of proportion of cell types in each group. Shown are Memory CD4⁺ T, Memory CD8⁺ T, NK and Memory B cells. Samples are shown in different colors. Horizontal lines represent median values, with a maximum of 1.5× interquartile range. t-tests were conducted between each group and p-values indicated. t test was performed using R (version 4.0). ^*p < 0.05, ^**p < 0.01, ^***p < 0.001.

Using unsupervised clustering of t-distributed stochastic neighbor embedding (t-SNE), 19 cell populations were identified by matching the expression of canonical gene markers for different cell types (Figure 1B, 1C and Supplementary Figure 1). The relative proportions of each immune cells were analyzed, revealing the discrepancies in cell composition between TCs and Covs and to compare them with that of HDs (Figure 1D, 1E). While significant differences were observed between COVID-19 patients (TCs and Covs) and HDs, indicating that SARS-CoV-2 immune memory was maintained in convalescent patients, there were limited differences between TCs and Covs. Compared to that in HDs, there is a significant increase COVID-19 patient in the proportion of memory CD8⁺ T cells (TCs and Covs), while memory CD4⁺ T cells were mildly elevated in TCs and Covs although there was no statistical difference (Figure 1E). SARS-Cov-2 specific memory T cells inducted after infection is a crucial participant in long-term protection [21]. Specific memory CD4⁺ T cells are necessary to stimulate potent B cell responses, resulting in subsequent antibody affinity maturation [22]. However, no significant difference in memory B cells was observed among three groups. This was not surprising as upon viral clearance, B cells will no longer be stimulated and proliferation is therefore terminated [23]. And a rapid decline in antibodies has been reported in mild cases [24]. Natural killer (NK) cell, which is a population of cytotoxic lymphocytes taking part in innate immunity against viral infection and tumor [25], were decreased in TCs, albeit not significantly (Figure 1E). Whether the reduction in NK cells is a risk factor for severe disease and mortality among patients with cancer and COVID-19 remains to be further determined.

Functional changes of T cells and B cells between convalescent Covs and TCs

Next, we performed hierarchical clustering to delineate the molecular differences of each cell type in Covs and TCs, according to the relative changes in gene expression compared to HDs. All T cell types and innate immune cell types, including NK, Monocyte, and DC cells clustered together based on disease group instead of by cell-types. with exceptions lied in B cells, hematopoietic stem cells (HSCs) and megakaryocytes (Figure 2A). Moreover, the difference between the two groups across cell types was limited (Figure 2A), which indicates the molecular characteristics of PBMCs in TCs are similar to that in Covs.

Figure 2. Analysis of functional changes between TCs and Covs. (A) Hierarchical clustering using the Pearson Correlation Coefficient (PCC) of a normalized transcriptome change between disease groups and HDs at cell type resolution. The color intensity indicates the PCC and the color bars above the heatmap indicate the cell type and disease group. (B) Differentially expressed genes in TCs compared to Covs in T and B cells. Red dots represent genes upregulated in TCs (logFC ≥ 0.25 and adjusted p < 0.05), while blue dots represent genes downregulated in TCs (logFC ≤ −0.25 and adjusted p < 0.05). Genes with |logFC| > 0.5 were labeled by gene symbols. (C, D) Enriched GO pathways of COVID-19 groups in T (C) and B (D) cells (left 4 columns: DEGs for TCs and Covs compared to HDs, right 4 columns: DEGs between TCs and Covs). The color intensity indicates the enrichment p-values and the point size indicates the ratio of gene enrich in each pathway.

Further investigation of the transcriptomic changes of adaptive immune cells during recovery from SARS-CoV-2 infection was conducted by comparing the expression patterns of T and B cells in Covs or TCs with that of HDs. Consistent with previous reports [26, 27], we found T cell and lymphocyte activation, lymphocyte and myeloid cell differentiation, interferon (IFN) responses and response to virus pathways in T cell types were significantly upregulated in Covs and TCs compared to HDs (Figure 2B, up and 2C). Interestingly, in TCs IL-6 and TNF production were suppressed compared to Covs and HDs potentially due an interaction between the immune system and the tumor (Figure 2C). IL-6 can stimulate B cell proliferation, differentiation, and antibody production [28], increase antitumor activity of CD4⁺ and CD8⁺ T cells [29], and tumor necrosis factor (TNF) contributes to activation of antigen presenting cell (APC) functions, the recruitment and activation of effector CD4⁺ T cells, effector CD8⁺ T cells, and the recruitment of NK cells as well [30]. Given above, our results indicate that TCs have a unique tumor-related molecular background, which is characterized by IL6&TNF associated immune suppression.

For B cells, similar to T cells, elevated and down-regulated pathways in Covs and TCs were basically the same compared to HDs, which suggested a consistent response of adaptive immunity to SARS-CoV-2 infection between Covs and TCs (Figure 2B, bottom and 2D). Cytokine production and nuclear factor (NF)-κB signaling pathways were significantly enriched in COVID-19 patients (Covs and TCs), while defense response to virus, B cell activation and differentiation, viral genome replication, and response to IFN-γ pathways were unexpectedly reduced in COVID-19 patients (Figure 2D). In accordance with the alteration in the proportion of memory B cells in Figure 1E, the immune response of B cells to SARS-CoV-2 had likely declined in convalescent patients due to virus clearance. During the outbreak in 2003, it was reported that those infected with SARS-CoV generated durable T cell responses lasting six years. In contrast, long-term memory B cells were absent [31]. These results suggest that compared with humoral immunity, cellular immunity may play a more important role in prevention of reinfection.

Change of immune cell composition and function after anti-tumor treatment

Long-term interruption of anti-tumor treatment risks tumor progression. In the meanwhile, patients may experience anxiety and stress during this time period. Conversely, immunosuppressive chemotherapy in convalescent COVID-19 patients risks inducing COVID-19 reactivation. To investigate the effect of recent cancer treatment, including cytotoxic chemotherapy and immunotherapy on immunity in convalescent COVID-19 patients, we performed 10X Genomics of 7 PBMCs from 3 COVID-19 convalescent patients with active gynecologic cancers pre- and post- treatment. Two of them received chemotherapy, of which one patient received a paclitaxel and cisplatin regimen (TP), another received a cisplatin, doxorubicin, and cyclophosphamide (PAC) regimen, and the other received immunotherapy (Supplementary Table 2). Cancer treatment of all three patients had been delayed due to pandemic (Supplementary Figure 2).

19 cell populations were identified based on the expression of canonical cell-type gene markers (Figure 3A) and the relative percentage of naïve CD4⁺ and CD8⁺ T cells decreased after treatment respectively in all three patients (Figure 3B, 3C). Of note, the relative abundance of CD14⁺ and CD16⁺ monocyte decreased sharply at 7 days after TP treatment (Figure 3B, 3C, up) and were restored at Day 27 after treatment. A similar change was observed in the patient treated with PAC. Due to severe bone marrow suppression 7 days after PAC treatment, this patient’s sample failed to meet sequencing requirements (Figure 3B, 3C, medium). The two chemotherapy regimens had little effect on the relative proportion of T and B cells (Figure 3B, 3C). In conformity to a previous report which demonstrated that PD-L1 blockade induced proliferation of CD8⁺ T and NK cell subsets [32, 33], NK cells, cytotoxic T cells, and NKT cells were massively increased in the third patient at Day 27 after treatment with anti-PD-L1 (Figure 3B, 3C, bottom).

Figure 3. Compositional and functional changes of cells pre- and post-treatment in COVID-19 patients with tumor. (A) t-SNE plot of cells colored by cell types identified from COVID-19 patients with tumor pre- and post-treatment. (B) Proportion of each cell type in each person before and after treatment. Bars are colored by cell types. (C) t-SNE plot of cells in each person at different time point before and after treatment. (D) Hierarchical clustering using the Pearson Correlation Coefficient (PCC) of a normalized transcriptome change between post- and pre-treatment in lymphocytes at cell type resolution. The color intensity indicates the PCC and the color bars above the heatmap indicate the cell type and disease group. (E, F) Enriched GO pathways of COVID-19 patients with tumor after treatment in T (E) and B (F) cells (top 4 columns: downregulated GO pathways after treatment, bottom 4 columns: upregulated GO pathways after treatment). The color intensity indicates the enrichment p-values and the point size indicates the ratio of gene enrich in each pathway. (G) Box plot of cytotoxic scores for each patient in cytotoxic T, NKT and NK cells. ^*p < 0.05, ^**p < 0.01, ^***p < 0.001.

Next, to evaluate molecular changes in each cell type after treatment, we performed hierarchical clustering according to relative gene expression changes compared with the pre-treatment sample. Interestingly, for all cell types of PBMC, cells clustered based on treatment groups instead of by cell-types (Figure 3D and Supplementary Figure 3A). Obviously, molecular characteristics of PBMCS varied significantly among treatment groups, which is regardless of cell-type. Hence, we attempted to identify variations in T and B cells in relevant biological functions after treatment through gene ontology (GO) analyses. Most importantly, the enriched pathways of up-regulation and down-regulation after anti-PD-L1 treatment were exactly the opposite of that after chemotherapy in T, B and monocyte cells (Figure 3E, Supplementary Figure 3B). We observed that genes upregulated after anti-PD-L1 treatment in T cells were involved in processes including T and lymphocyte cell activation, T cell receptor signaling pathway, and neutrophil mediated immunity (Figure 3E), and in B cells were enriched in neutrophil mediated immunity, B cell activation, and lymphocyte proliferation (Figure 3F). Neutrophil mediated immunity pathway was also enriched in monocyte and dendritic cells (DC) after anti-PD-L1 treatment (Supplementary Figure 3B, 3C). Significant activation of neutrophil mediated immunity in T and B cells after anti-PD-L1 treatment suggests that neutrophils may act as APCs via direct interaction with T and B cells or exert regulatory effects on adaptive immunity [34, 35].

In accordance with the DEGs enrichment results, we found that CD8⁺ T cells, NKT cells, and NK cells all showed significantly higher cytotoxic scores after anti-PD-L1 treatment in TC3 patient, while all three subsets had lower cytotoxic scores after TP treatment in TC1 patient (Figure 3G). No obvious changes were found in TC2 patient, except that the cytotoxic score in NK cells was significantly decreased after PAC treatment (Figure 3G).

Change of clone status of T cells after treatment in TCs

CD8⁺ T cells confront virus-infected cells with the viral antigen recognition function of T cell, which is meditated by T cell receptor (TCR). The immune response towards specific antigen is mainly determined by the diversity as well as the size of T cell receptor (TCR) repertoire [36, 37]. To gain insight into the change of clonal status of individual T cells across three conditions, TCR sequences were reconstructed based on original TCR sequencing. Briefly, most cells in all T cell subsets have matched TCR information. The exception lied in NKT and γδT cells (Figure 4A). Among T cell subsets, there are varying degrees of clonal expansion observed. Large clonal expansions with a clonal size ≥30 were mostly concentrated in cytotoxic T cells (Figure 4A), which indicated that effector T cells underwent dynamic and clonal responses in time of SARS-CoV-2 infection [26, 38].

Figure 4. Clonal dynamics of T cells pre- and post-treatment in COVID19 patients with tumor. (A) t-SNE plot of T cells colored by cell types (top) and clone status (bottom). (B) The proportion of each clone status (top) and clone type (bottom). Bars are colored by the top 20 most abundant clones in each patient at different time points. (C) Scatterplots comparing TCR clone frequencies pre- and post-treatment. Shown are amplified clones (n > 1). (D) t-SNE plot of single cells colored by selected TCR clones. Shown are the top 5 most abundant clones before treatment in each patient at different time point. (E) Phenotypes of single cells belonging to the same TCR clone. Shown are the top 20 most abundant clones before treatment in each patient at different time point.

Next, the global clone abundance was investigated based on the comparison of pre- and post- treatment frequencies of each clone in each patient (Figure 4B, 4C and Supplementary Figure 4A). We found that TP and anti-PD-L1 treatment enhanced the clonal expansion of cytotoxic T cells, while a slightly reduction was observed after PAC treatment (Figure 4B). Interestingly, most clones that were highly expanded before therapy were preserved after treatment (Figure 4C and Supplementary Figure 4A). Integration of scRNA-seq data revealed that more than 90% of highly expanded clones (n > 4) remained after treatment and were even maintained in a similar frequency after treatment (Figure 4D, 4E and Supplementary Figure 4B). Namely, across all patients, TOP 5 or even TOP 20 clones were still enriched in the same cell subsets (Figure 4D, 4E). Taken together, TP and anti-PD-L1 therapy had limited effects on the clonal status of T cells.

Clinical information and data access

Approval of this study was given by the Institutional Review Board of Tongji Hospital, Tongji Medical College, Huazhong University of Science and Technology (TJ-IRB20200405). Patients were informed of necessary information in advance and an inform consent for each patient is signed. All blood samples in this study was acquired by collecting residual samples in clinical diagnostic tests, which ensured no extra burden to patients. The medical history is gained by accessing electronic medical record system which is under informed consent of patients. All data and information were used and analyzed according to ethical guidelines.

As June 7, 2020, whole blood samples of 4 COVID-19 convalescents with gynecological cancer, 4 COVID-19 convalescents without cancer and 3 healthy donors as normal controls were enrolled and 10X scRNA-seq was performed with PBMCs isolated from their blood samples. Routine laboratory tests were acquired through medical record and included in the analysis.

10X genomics scRNA library construction sequencing

Blood samples were collected using heparin tubes (Becton, Dickinson and Co.) for later management. For PBMC isolation and extraction, blood samples were overlaid with Ficoll-Paque Plus medium (GE Healthcare). Each heparin tube was centrifugated by density gradient, followed by Ca/Mg-free PBS resuspension and wash. The whole sample collection and processing procedure was mostly within 1h and was strictly limited to 4 h at most. The viability of PBMC cells were identified using 0.4% Trypan blue coloring. After microscope examining, samples with ≥80% viability was qualified to be incorporate into single cell library construction. In brief, in order to generate single-cell gel beads in emulsion (GEMs), the single cell suspension was loaded together with regents and Gel Beads by implementation of Chromium™ Single Cell 5’ Reagent Version 2 Kit (10x Genomics, Pleasanton, CA). Each Gel Bead contains barcoded oligonucleotides. The above operations are performed in Chromium™ Controller following manufacturer’s instructions (GemCode Technology). Inside each GEM, cells are lysed with a lysis buffer followed by barcoded reverse transcription. For each single cell, the reverse transcription of polyadenylated mRNA was distinguished by the unique barcode. Fragmented cDNA was end repaired by adding A-tailing to the 5’ end. The adaptor ligated DNA was then subjected to double sided SPRI selection. A final SPRI purification was conducted after sample index PCR. The resulting single-cell RNA-seq library was sequenced on the MGISEQ-2000 sequencer (BGI, Shenzhen, China) to insure qualification.

Single-cell RNA-seq data processing

The single cell RNA-seq sequencing data were processed by implementing Cell Ranger (version 3.0.1, 10x Genomics) with the GRCh38 human as reference genome and the gene expression matrices were analyzed with Seurat (version 3.2.2). Cells were considered as disqualified and were removed from data processing if they met one of the following criteria: 1) UMIs below 500 or above 30,000; 2) less than 200 genes expressed; 3) more than 15% of UMIs from mitochondrial genes on the UMI counts per cell. After removal of disqualified cells, the gene expression matrix of each sample was normalized through the “NormalizeData” and “FindVariableFeatures” function with default parameters. The normalized data were then integrated using “FindIntegrationAnchors” and “IntegrateData” function with parameters set to “nfeatures = 300, dims = 1:30”. Ribosomal reads and mitochondria protein were eliminated when the integrated dataset was scaled. For PCA analysis, cells were clustered by running “FindNeighbors” and “FindClusters”, with parameter set to “k.param = 10, resolution = 1”. To visualize the clusters, data was plotted on a tSNE space, which incorporated the top 30 principle components.

Cell type annotation

Unbiased cell type annotation was performed with the “FindAllMarker” function in Seurat. Clusters were annotated according to expressions of identified markers as well as some canonical markers of particular cell types. For better insight of the data, HSC cells and clusters which express more than 2 canonical cell-type markers were excluded from further analysis.

Identification of differential expressed genes (DEGs) analysis and GO enrichment

Differential gene expression was carried out using “FindMarkers” function with MAST algorithm in Seurat. The false discovery rate was estimated by employing Benjamini and Hochberg procedure. Differential expressed genes (DEGs) were filtered by adopting a minimum logFC of 0.25 and a maximum FDR value of 0.05. GO enrichment of DEGs was conducted using function “enrichGO” in ClusterProfier. The parameters were set to “OrgDb = org.Hs.eg.db, ont = ’BP’, pAdjustMethod = ’BH’”.

Hierarchical clustering of gene expression profiles at cell type resolution

Unsupervised hierarchical clustering of gene expression changes was carried out following method described previously [45]. In brief, the normalization of UMI count was conducted by multiply total UMI count of each cell type by 100,000. For each disease group, the gene expression was divided by the values in the HDs for TCs and Covs or those in the pre-treatment group respectively. Highly variable genes were determined by calculating the top 3000 standard deviation and then the data was log2-transformed. Hierarchical clustering of genes was performed according to the Pearson’s correlation coefficient (PCC) of the preciously determined highly variable genes.

Defining cell cytotoxic scores

The cytotoxic scores were calculated according to the average expression of 12 cytotoxicity-associated genes (GZMB, GZMA, GZMH, PRF1, IFNG, KLRK1, KLRB1, KLRD1, GNLY, NKG7, CTSW, CST7,). The calculation was based on the function “AddModuleScore” in Seurat with default parameters.

TCR V(D)J sequencing and data analysis

The TCR V(D)J gene segments were enriched from amplified cDNA libraries for subsequent analysis using Chromium Single-Cell V(D)J Enrichment kit (10X Genomics) following the manufacturer’s protocols. The resulting libraries were sequenced using the MGISEQ-2000 sequencer (BGI, Shenzhen, China). Demultiplex, quantification and clonotype assignment were carried out by Cell Ranger (v2.2.0) vdj pipeline with default parameters and with GRCh38 genome as reference. The frequency of each barcode was calculated and a matrix of barcode information and corresponding clonal type information was generated. Further analysis was performed on cells with following prerequisites: 1) at least one productive TCR alpha chain (TRA); 2) at least one productive TCR beta chain (TRB). A clonal type wad defined for each unique TRA(s)-TRB(s) pair.

Statistical analysis

Statistical analysis was performed with Wilcoxon rank-sum test and t test using R (version 4.0) in this study, with p-values indicated as ^*p < 0.05, ^**p < 0.01, ^***p < 0.001.