Reversing Immunosenescence with Senolytics to Enhance Tumor Immunotherapy

Human specimens

This study includes three independent clinical trials for HNSCC: OOC-001 (NCT04718415), REDUCTION-I (NCT05582265), and COIS-01 (NCT05724329). Eligible HNSCC patients received prior approval from the Institutional Review Board (IRB) at Sun Yat-sen Memorial Hospital, Sun Yat-sen University, and were required to provide written informed consent to participate in the study. All patients were diagnosed with HNSCC by expert histopathologists, and all invasive procedures and systemic treatments were ordered by the patients’ attending physicians, individually determined based on the response after the administration of pharmacological drugs. The eight patients used for the scRNA-Seq study were part of the OOC-001 clinical trial and received the described treatment regimen. Specifically, all patients had their primary tumor tissue and peripheral blood collected through biopsy at the time of enrollment, and the peripheral blood was collected after sintilimab in combination with carboplatin and nab-paclitaxel treatment before surgery. Among the 7 patients (SH1-SH7) who underwent surgery, 3 patients are MPR, and 4 patients are noMPR. Paired biopsy samples of primary tumor tissue and peripheral blood were collected at baseline before treatment, and peripheral blood samples were collected 3 weeks after the first dosing cycle. Except for SH1, peripheral blood samples were collected from the other patients before and after the first medication cycle for matching. For patient SH1, blood samples were collected only before the first medication cycle. Peripheral blood samples were collected after the second medication cycle for SH1, SH2, SH3, and SH4, and peripheral blood samples were collected for SH3, SH4, and SH7 after the third medication cycle. Paraffin sections from two ongoing clinical study patient cohorts were used for immunofluorescence staining of senescent markers P16 and P21. This included 33 MPR patients and 17 noMPR patients from OOC-001, as well as 37 MPR patients and 24 noMPR patients from REDUCTION-I. Another set from COIS-01 collected paraffin sections from 7 patients for immunofluorescence staining before and after treatment.

Sample collection and processing

Peripheral blood mononuclear cells (PBMCs) were isolated according to the manufacturer’s protocol. Briefly, before the start of treatment, 5 mL of fresh peripheral blood was collected in EDTA anticoagulated tubes and subsequently extracted with Lymphoprep™ (Stem Cells) solution. After centrifugation, carefully transfer lymphocytes to a new 50 mL tube and wash with PBS. Next, incubate the lymphocytes for 10 minutes to lyse the blood cells. Finally, cells were resuspended in sorting buffer containing PBS plus 2% FBS.

For tumor biopsy samples, cut tissue with a diameter of about 5 mm at the junction of the tumor and normal tissue as a representative sample. All 7 tumor biopsy samples were collected as a discovery cohort through routine clinical practice and then stored in MACS Tissue Storage Solution (Miltenyi). Fresh biopsy samples of oral HNSCC were minced and isolated using a human tumor isolation kit (Miltenyi), followed by enzymatic digestion on a rotor at 37°C for 45–60 min. The detached cells were then passed through 70 μm and 40 μm cell filters (Biosharp) and centrifuged at 300 × g for 10 min. Pelleted cells were suspended in red blood cell lysis buffer (Thermo-Fisher) and incubated on ice for 2 min, then washed and resuspended in phosphate buffer supplemented with 2% fetal bovine serum (FBS, ExCell) PBS twice in saline (PBS, Gibco). The viability of all samples was confirmed to be >90% using the trypan blue (Thermo-Fisher) exclusion method.

Single-cell RNA-seq library and TCR-seq library preparation and sequencing

According to the manufacturer’s protocol, scRNA-seq libraries were constructed. In summary, the density of cells was determined after washing once with PBS containing 0.04% bovine serum albumin (BSA, Invitrogen). Next, 2×10⁵ cells were loaded on a 10x Genomics GemCode Single cell instrument (10x Genomics) that generates single-cell Gel Bead-In-EMulsion (GEMs). Libraries were generated and sequenced from the cDNAs with the Chromium Next GEM Automated Single Cell 3’ cDNA Kitv3.1 (10x Genomics). Upon dissolution of the Gel Bead in a GEM, primers were released and mixed with cell lysate and Master Mix. Barcoded, full length cDNAs were then reverse-transcribed from polyadenylated mRNA. Full-length, barcoded cDNAs were then amplified by quantitative real-time PCR (qRT PCR) for library construction. Finally, single-cell RNA libraries were sequenced by an Illumina HiSeq X Ten sequencer with 150 bp paired-end reads.

After reverse transcription and cell barcoding, the cDNA was de-emulsified and purified using silane magnetic beads, followed by PCR amplification. The amplified cDNA was then used for 5’ gene expression library construction and TCR V(D)J targeted enrichment amplification using the Illumina® bridge. Libraries prepared according to the manufacturer’s user guide were then purified and analyzed for quality assessment. Single-cell RNA and TCR V(D)J libraries were sequenced by an Illumina HiSeq X Ten sequencer using 150 bp paired-end reads.

Bulk TCR-seq library preparation and sequencing

Whole Blood Samples Blood was collected in EDTA tubes and sent to isolated PBMCs, and DNA was extracted for TCRβ analysis by immunoSEQ ⁶⁵. Briefly, extracted genomic DNA was amplified using a bias-controlled multiplex PCR system, followed by high-throughput sequencing. Raw data processing and analysis were performed by immunoSEQ. Subsequently, we used the International ImMunoGeneTics Information System (IMGT; http://imgt.cines.fr) database for comparison. After the alignment, the relative abundance of each TCRβ CDR3 sequence was clarified and calculated. Furthermore, we use batch correction to eliminate batch effects across different datasets. Subsequently, multiple TCR data statistics were performed.

Bulk ATAC and RNA isolation and sequencing

Isolation, library construction and sequencing of bulk ATAC and RNA were performed at Gene Denovo Biotechnology Co. (Guangzhou, China). For ATAC sequencing, CD4 naïve cells from C57 mouse tongue tumors were extracted. Tn5 transposase was added to the nuclear suspension. After the reaction was completed, the DNA fragment was purified; the amplified product was then used for PCR amplification. The fragments were purified using AMPure XP magnetic beads (Beckman Coulter, Brea, CA, USA) to construct a sequencing library. After the library construction was completed, Agilent 2100 (Agilent, Santa Clara, CA) was used to detect the quality of the library. Libraries that pass the quality inspection will be used for on-machine sequencing (Novaseq 6000) to obtain sequence information of the open chromatin region fragments to be tested.

Extraction of CD4 naïve cells from C57 mouse tongue tumors. Total RNA was extracted using Trizol reagent kit (Invitrogen) according to the manufacturer’s protocol. RNA quality was assessed on an Agilent 2100 Bioanalyzer (Agilent Technologies). After total RNA was extracted, eukaryotic mRNA was enriched by Oligo(dT) beads. Then the enriched mRNA was reversely transcribed into cDNA by using NEBNext Ultra RNA Library Prep Kit for Illumina (NEB #7530, New England Biolabs). The resulting cDNA library was sequenced using Illumina Novaseq6000 by Gene Denovo Biotechnology Co. (Guangzhou, China).

scRNA-Seq data processing

The Cell Ranger toolkit (version 3.1.0) provided by 10x Genomics was applied to aggregate raw data, filter low-quality reads, align reads to the human reference genome (GRCh38), assign cell barcodes and generate a UMI matrix. The toolkit Scanpy (version 1.9.3) ⁶⁶ was used to analyze single-cell RNA sequencing (scRNA-Seq) data. Specifically, the original UMI matrix was processed to filter out genes expressed in less than 10 cells and cells with less than 200 genes. We further quantified the number of genes and UMI counts per cell and maintained thresholds for high-quality cells of 1000-50,000 UMIs, 400-6,000 genes and less than 10% mitochondrial gene counts, ensuring that most heterogeneous cell types were Incorporate into downstream analyses. Scrublet (version 0.2.3) ⁶⁷ was then applied to each sequencing library to remove potential doublets, using default parameters. The normalized expression matrix was calculated from raw UMI counts normalized to total counts per cell (library size) and then scaled by 1e4 and log-transformed.

Dimensionality reduction and unsupervised clustering of scRNA-seq data

Dimensionality reduction and unsupervised clustering were performed according to standard workflows in Scanpy ⁶⁶. Briefly, use highly variable genes (HVG) based on the highly_variable_genes default parameters. Unwanted sources of variation, including total counts, and mitochondrial gene counts were further regressed from the normalized expression matrix. Primary component analysis (PCA) was performed on the variable gene matrix to reduce noise and the top 50 components were used for downstream analysis. To correct for batch effects across different patients, we applied BBKNN to generate a batch-balanced k-nearest neighbor (KNN) map, which identifies the top neighbors of each cell in each batch individually, rather than the entire pool of cells ⁶⁷. The Leiden algorithm is then applied to such nearest neighbor graphs to detect communities and find cell clusters ⁶⁸. It is worth noting that the same principal components are also used for nonlinear dimensionality reduction to generate uniform distributions for visualization in Uniform Manifold Approximation and Projection (UMAP). After the first round of unsupervised clustering, we annotated each cell cluster according to canonical immune cell markers and identified the main immune cell types, including T cells, NK cells, B cells, myeloid cells, fibroblasts, endothelial cells and tumor cells. Use the rank_genes_groups function with parameter method Wilcoxon to detect marker gene.

TCR sequences assembly

Applying the TCR sequence assembly Cell Ranger toolkit (version 3.1.0) provided by 10x Genomics will perform FASTQ sequence quality filtering, sequence alignment, V(D)J assembly and TCRαβ pairing. Only TCRs containing paired TRAV-CDR3-TRAJ and TRBV-CDR3-TRBJ-TRBC chains were considered valid and retained for downstream analysis. Each cell is assigned a pair of α and β chains with the highest UMI count. Cells with identical TCR pairs are defined as clonal and are thought to originate from a common ancestor.

Calculation of Diversity Index

Characterize the immune repertoire by examining diversity and clonality. Immune repertoire analysis involves a variety of diversity indices, such as Chao1 index, Clone Diversity Index, Richness index and Shannon index ²⁴.

Bulk RNA-sequencing and ATAC-sequencing data processing

For bulk RNA-Seq data, reads were filtered using FASTQ (0.23.3) default parameters, and then clean reads were aligned to the mouse genome GRCm39 and gencode. vM30 annotations were using STAR (2.7.10b) and using RSEM (1 .3.1) rsem-calculate-expression is in quantitative analysis of gene expression. Difference analysis was done using DESeq2 (1.42.0).

For bulk ATAC data, use FASTQ (0.23.3) for quality control, then use bwa (0.7.17) to align the reads to the reference genome (GRCm39), and peak using Genrich (0.6.1, https://github.com/jsh58/Genrich). Peak is annotated to gencode. vM30 through annotatePeaks.pl of HOMER (4.11). The difference peak is obtained by BEDTools (v2.31.1) subtract. (https://github.com/jsh58/ATAC-seq2)

Generation of Immunosenescence gene set

Our own senescence-related gene set was generated by combining genes reported to be enriched in senescent immune cells in previous studies and experimentally validated in at least human or mouse cells. We screened 1745 studies, but after removing studies reporting duplicates, case reports, non-human genes, and non-high-throughput sequencing, a list of 6 studies was developed (Senmayo ⁴⁶, GSE157007 ²⁴, HRA000624 ⁶⁹, HRA000395 ⁷⁰, HRA000615 ⁷¹ and HSA000203 ⁷² and 1 database (https://genomics.senescence.info/cells/)). We screen immune cells and compare samples from young and old people to obtain differential genes that are highly expressed in the elderly, remove duplicate genes, and screen genes that appear in at least two data sets. Based on the enrichment of relevant gene pathways and machines Learning method, 154 genes constituting the immunosenescence gene set were identified (Supplementary Table 6).

GO enrichment analysis and Gene set enrichment analysis (GSEA)

Using the R package ClusterProfiler (version 4.8.3), gene function annotation and GO were performed on genes that were down-regulated and up-regulated in T cells and B cells in the MPR and noMPR groups, as well as differential genes between macrophages and monocytes and other cells. Analysis, a program that supports statistical analysis and visualization of functional profiles of genes and gene clusters.

GSEA was performed by the GSEA software (http://software.broadinstitute.org/gsea/index.jsp) ⁷³. Gene sets used in this article were c2.cp.kegg.v6.2.symbols.gmt downloaded from the Molecular Signatures Database (MSigDB, http://software.broadinstitute.org/gsea/msigdb/index.jsp).

Gene Set Score Analysis

The scoregenes function in Scanpy is used to calculate the module score of a gene expression program in a single cell. First, all analyzed genes were binned according to their average expression, and control genes were randomly selected from each bin. Then, the average expression value of the gene set at the single-cell level minus the aggregate expression of the control gene set was calculated. The gene set was obtained from the immunosenescence gene set constructed by ourselves. The genes in each gene set are listed in Supplementary Table 6.

Cell-cell interaction

We used the CellChat package (1.6.1), which allows analysis of scRNA-Seq data. We studied the interaction between immune cells (i.e. T cells and B cells) and myeloid cells (i.e. macrophages and monocytes), one of which contributes ligands or receptors during the interaction, using CellChat Signaling pathway networks are analyzed and visualized.

Trajectory inference of immune cell subsets across tissues

The status of immune cell subpopulations in the periphery and in tumors is dynamic, and they may differentiate into different cell states and exert different biological functions in different patients. We use the pseudo-temporal inference algorithm Monocle 2 for trajectory analysis to reconstruct the cell differentiation trajectory of immune cells across tissues, reveal the progression of cells and reconstruct the trajectory of cells progressing through biological processes under study.

Survival analysis

Survival analysis uses TCGA and HNSCC data to evaluate the prognostic performance of different immune cell populations. Gene expression data and clinical information of HNSCC patients were downloaded from cBioPortal (https://www.cbioportal.org/study/clinicalData). Use Cibersort and Xcell to extract T cells and B cells subtypes and group them based on the expression levels of characteristic genes. Specifically, HNSCC patients were divided into four groups based on the median proportion of immune cell subsets and the median mean expression of a given signature gene. In the analysis of myeloid cells, a subset of cells within a single cell is characterized using a set of signature genes, defined based on a combination of typical marker genes for major immune cell types and cluster-specific marker genes. Specifically, for different macrophage clusters, we used the typical IGF1⁺ Mφ marker genes “CD276”, “APOE”, “TREM2”, “HTRA1”, “PTGR1”, “RBP1”, “IGF1”, “CDKN1A”, “CD81” define its signature. HNSCC patients were divided into two groups based on the median of the mean expression of the given signature. Survival analysis was performed using the COX proportional hazards model implemented in the R package “survival”.

Hierarchical Clustering

In the end, cell subpopulations and characteristic genes were used to construct a neoadjuvant treatment model for patients before surgery, and a hierarchical clustering method was used to create groups. Cophenetic distance was used to evaluate the objects in the group so that they were like each other and different from objects in other groups. Verify different cell subpopulations and gene signatures to evaluate the efficacy of neoadjuvant therapy in patients before surgery.

TNBC single cell RNA sequencing data processing

TNBC single-cell RNA sequencing data comes from GSE169246. Patients who received combined neoadjuvant chemotherapy and immune therapy before surgery were screened, and the patients were divided into PR (response) and SD (non-response) groups based on the corresponding imaging conditions. The rest of the analysis process as before, similar cell subpopulations were evaluated after dimensionality reduction clustering, and senescent immune cells were evaluated on the cell subpopulations. Hierarchical clustering was used to verify the effect of different cell subpopulations and gene characteristics on patients receiving neoadjuvant therapy before surgery to evaluate.

Animal experiments

The animal experiments in this study were conducted with the approval of the Animal Ethics Committee of Sun Yat-sen University and the Institutional Animal Care and Use Committee at the Shenzhen Institutes of Advanced Technology, Chinese Academy of Sciences (CAS). For all animal experiments, mice were housed under controlled conditions of ∼50% humidity, ∼25 °C temperature and a 12h light/dark cycle. Isoflurane was used for anaesthesia. At the end of all experiments, animals were euthanized under CO₂ anaesthesia.

4-NQO-induced mouse HNSCC model

Female C57 mice (12 months old) were housed and drinking water containing 50 mg/L 4-NQO was given for 12 weeks. All diet and water were made available ad libitum. Subsequently, the oral cavities of the mice were examined. 24 mice with lesions were randomly divided into 4 groups. Normal drinking water was given. Two days later, αPD-1 was administered intraperitoneally at a dosage of 10 mg/kg. Then, mice were euthanized. Tongues were collected for pathological examination and IF.

Establishment of cell lines and tissue culture

Female C57 mice (6-8 weeks old, weighing 18-20 g) were obtained from Gempharmatech (Guangdong, China). 4-nitroquinoline 1-oxide (4-NQO) was diluted in drinking water to a final concentration of 50 mg/L for the animals, with weekly replacements. After 12 weeks, all animals were returned to normal water until week 28. At week 28, animals were euthanized for tissue extraction. Lesions were dissected, and cells were isolated to establish the 4MOC cell line.

SA-β-gal staining

Fresh tissues were fixed in 10% neutral buffered formalin (NBF) for 3–4 h and then transferred to 30% sucrose overnight. Tissues were then embedded in cryo-embedding medium (OCT) and cryosectioned at 5 μm for staining of SA-β-gal (pH 6.0) at 37 °C for 16–24 h in SA-β-gal staining solution (40 mM citric acid in sodium phosphate buffer, 5 mM K₄[Fe(CN)₆]3H₂O, 5 mM K₃[Fe(CN)₆], 150 mM sodium chloride, 2 mM magnesium chloride and 1 mg ml⁻¹ X-gal dissolved in N,N-dimethylformamide). Slides were imaged at 20× with a Panoptiq slide scanner (ViewSiq).

RNA isolation and qRT–PCR

Tissues were collected from euthanized mice and flash-frozen in liquid nitrogen. Tissues were homogenized using FastPrep-24 homogenizer (MP Biomedicals) and total RNA was isolated by Trizol extraction according to manufacturer’s specifications (Thermo Fisher). Total RNA was quantified using a Nanodrop spectrophotometer (Thermo Fisher) and 1 μg of total RNA used to generate cDNA via the Transcriptor First Strand cDNA synthesis kit (Roche) according to the manufacturer’s specifications. Gene expression changes in (Cdkn2a) p16^Ink4a, (Cdkn1a) p21^Cip1 were quantified by qRT–PCR reactions using 20 μl reaction volumes and a StepOne thermocycler (Thermo Fisher) with input of 50 ng total RNA per reaction. For each sample, reactions were performed in duplicate. Data were analysed by the ΔΔCt method and expression was normalized to Gapdh.

Quantitation of 8-OHdG DNA lesions

Tissues from mice were analysed for 8-OHdG levels using the ELISA kit (Abcam) according to manufacturer’s specifications.

Immunoblotting

Snap-frozen livers from mice were incubated in RIPA buffer (Thermo-Fisher) on ice for 30 min after being homogenized with a FastPrep-24 homogenizer. Samples were centrifuged at 17,000g for 15 min at 4 °C. Supernatant was resuspended in 2× SDS loading buffer and 50 μg of total protein run on a 4–15% SDS–PAGE gel (Bio-Rad) before being transferred to nitrocellulose membrane. Membranes were blocked for 1 h in 10% milk TBS-T solution at room temperature before incubation in anti-γH2AX (Novus Biologicals, NB100-384, 1:2,000) and anti-GAPDH (Abcam, ab8425, 1:5,000) antibody at 4 °C overnight. After washing, samples were incubated in either horse anti-mouse HRP (Cell Signaling Technology, 7076S) or goat anti-rabbit HRP secondary antibody (Thermo-Fisher, 656120, 1:2,000) in 5% milk TBS-T solution for 3 h before washing and visualization with ECL (Thermo-Fisher).

TIL isolation and flow cytometry

Tumors were dissected, minced, and resuspended in complete media (DMEM with 10% FBS and 1% antibiotics) supplemented with Collagenase-D (1 mg/mL; Roche) and incubated at 37 °C for 30 min with shaking to form a single-cell suspension. Tissue suspensions were washed with fresh media and passed through a 100-µm strainer. Samples were washed with PBS and immediately processed for live/dead cell discrimination using BD HorizonTM Fixable Viability Stain 510. Cell surface staining was done for 30 min at 4 °C with the following antibodies: CD45 (30-F11) (1:100), CD3 (145-2C11) (1:400), CD8a (53-6.7) (1:100), CD4 (RM4-4) (1:400), NK1.1 (PK136) (1:400), CD62L (MEL-14) (1:200) and CD44 (IM7) (1:100). All flow cytometry data acquisition was done using Beckman Cytexpert and analyzed using FlowJo software (v10.8.1).

Immunofluorescence and image quantification

Briefly, tissues were harvested, fixed, and paraffin embedded. Slides were stained for P16, P21, CD4, CD8a, CCR7, PAN-CK, LCK, CD62L and S100A11 antibodies (Supplementary Materials). Quantification of immune infiltration was done using QuPath, an open-source software for digital pathology image analysis. For the quantification, at least three regions of interest (ROI) were selected for each condition and the percentage of positive cells was calculated.

DNA sequencing and genomic signature analysis

Raw sequencing data were aligned to the mm10 reference genome using BWA. Somatic mutations were identified by comparing the sequencing data from each cancer sample to the sequencing data from a normal tissue derived from the tail of one of the mice (all mice were genetically identical). To ensure robustness of the results, a consensus variant calling strategy was leveraged in which somatic mutations were identified using three independent bioinformatics tools: Strelka2, Varscan2, and GATK4 Mutect2. Somatic mutational profiles were derived using the immediate sequencing context by evaluating the base 5′ and the base 3′ to each single point mutation.