Tumor-associated neutrophil precursors impair homologous DNA repair and promote sensitivity to PARP-inhibition

Abstract
Tumor evolution is one of the major mechanisms responsible for acquiring therapy-resistant and more aggressive cancer clones. Whether the tumor microenvironment through immune-mediated mechanisms might promote the development of more aggressive cancer types is crucial for the identification of additional therapeutical opportunities. Here, we identified a novel subset of tumor-associated neutrophils, defined as tumor-associated neutrophil precursors (PreNeu). These PreNeu are enriched in female highly proliferative hormone-dependent breast cancers and impair DNA repair capacity.  Mechanistically, succinate secreted by tumor-associated PreNeu inhibit homologous recombination, promoting error-prone DNA repair through non-homologous end-joining regulated by PARP-1. Consequently, breast cancer cells acquire genomic instability promoting tumor editing and progression. Selective inhibition of these pathways induces increased tumor cell killing in vitro and in vivo. Tumor-associated PreNeu score correlates with copy number alterations in highly proliferative hormone-dependent tumors from breast cancer patients. Treatment with PARP-1 inhibitors counteract the pro-tumoral effect of these neutrophils.

DESIGN:

Enclosed are raw count gene expression matrices from 3 single-cell sequencing experiments.
- Experiment 1: SMART-Seq2 RNA-sequencing of human mature neutrophils (MatNeu) and neutrophil precursors (PreNeu) derived from luminal B breast tumor biopsy (raw counts).
- Experiment 2: SMART-Seq2 RNA-sequencing of human mature neutrophils (MatNeu) and neutrophil precursors (PreNeu) sorted from luminal B breast tumor biopsies or derived from human cord blood mononuclear cells (vst-normalized gene expression values).
- Experiment 3: Single-cell RNA sequencing using BD Rhapsody on a highly proliferative ER+ breast biopsy (.RDS object).

METHODS:

Rhapsody single-cell data generation, processing and analysis

Sample processing: fresh biopsy single-cell suspension was enriched for CD45+ cells using the CD45 MicroBeads, human (Cat. 130-045-801; Miltenyi) following the manufacturer’s protocol. Cells were single cell captured for sequencing using the Rhapsody HT Single-Cell Analysis system and library was generated using the BD Rhapsody WTA Amplification Kit (Cat. 633801, BD Biosciences) with 8 cycles of amplifications, and sequenced using an Illumina NextSeq2000 instrument with a P2 flow cells and chemistry XLEAP obtaining around 8000 reads/cell.

Computational analysis: Tumor sample FASTQ files were aligned, and feature-barcode matrices were generated using the BD Rhapsody™ Sequence Analysis Pipeline on the Seven Bridges Genomics platform, with the GRCh38 genome assembly as the reference. The resulting data was analyzed using Seurat (v.5.1.0) [10.1038/nbt.4096, 10.1016/j.cell.2019.05.031, 10.1016/j.cell.2021.04.048] in R (v.4.4.2) [R Core Team (2020). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL: https://www.R-project.org/]. Quality control measures were applied to each dataset to remove low-quality cells and multiplets. Metrics such as cell counts, UMI counts per cell, genes detected per cell, mitochondrial gene count ratio, and ribosomal gene count ratio were inspected. Data normalization and scaling were performed, regressing out potential confounding factors (UMI counts per cell, genes detected per cell, mitochondrial gene count ratio, and ribosomal gene count ratio). A total of 30 principal components (PCs) was selected for Uniform Manifold Approximation and Projection (UMAP), which was used for dimensionality reduction and cell clustering. Clustering at a resolution of 0.6 yielded 20 clusters. Marker gene identification utilized a hurdle model designed for scRNA-seq data, implemented in the MAST statistical framework [10.1186/s13059-015-0844-5], with Bonferroni-adjusted p-values to correct for multiple testing. Markers were considered significant if they were expressed in at least 70% of cells in a cluster, had an adjusted p-value (pval_adj < 0.05), and displayed a log2 fold change (log2FC > 1). Major cellular populations were annotated based on marker genes, and the dataset was subsetted to retain only bona fide neutrophils. The subsetted cells underwent further graph-based clustering at a resolution of 0.6, resulting in 3 clusters. Marker-based annotation was applied for in-depth characterization. Transcriptomic data underwent zero-preserving imputation using the ALRA method [10.1038/s41467-021-27729-z]. The Seurat AddModuleScore function was employed to compute and evaluate the enrichment of neutrophil-related signatures. Pathway enrichment analysis was conducted using Metascape [10.1038/s41467-019-09234-6], considering pathways with a Benjamini-Hochberg-adjusted p-value < 0.05, involving at least three differentially expressed genes (DEGs), and a minimum enrichment score of 1.5. Finally, trajectory analysis was performed using a custom script based on Monocle3 [10.1038/s41586-019-0969-x] (v.1.2.7). Cell spatial coordinates were imported from the Seurat object to create the CellDataSet object required by Monocle. Trajectories were constructed using the learn_graph function, and cells were ordered along pseudotime.

Single-cell SMART sequencing (SMART-Seq2)

Sample processing: fresh biopsy single-cell suspension or cord blood derived-IMCs were sorted as described above into 96-well PCR plates containing cell lysis buffer. Samples collected in cell lysis buffer were used for RNA-seq library preparation with the NEBNext Single Cell / Low InputRNA Library Prep Kit for Illumina (NEB, E6420S), following the manufacturer’s protocol for single cells, with the following parameters adjusted: 17 cycles for cDNA amplification PCR, 11cycles for library enrichment PCR. Libraries were dual-indexed (NEBNext Dual Index Primers Set 1, NEBE7600S), and sequenced on an Illumina NextSeq 500 instrument with 75 cycles reagents. Quality controls and read mapping to the reference genome were performed using the same criteria described for bulk RNA-seq, except that reads having same start/end coordinates and identical nucleotide sequence were marked and deduplicated to avoid excessive bias due to PCR amplification. Differential expression was performed in R statistical environment using DESeq2 pipeline (v1.28.1). Cellswith less than 500K mapped reads were removed from the analysis. Library-size normalized data were transformed using the variance stabilizing transformation and batch effect between tumor- sorted PreNeu or LOX-1+ Neutrophils and cord blood derived-PreNeu or LOX-1+ Neutrophils was corrected using Combat-Seq. Cell-types of origin (PreNeu/PMN-MDSC), which were independently identified in both conditions, were set as covariates to preserve biological signal.

Computational analysis: Marker genes specific for PreNeu were identified from the single-cell data. We selected all genesbeing differentially expressed in PreNeu cells vs LOX-1⁺ neutrophils (FDR < 0.05). To identify genes being robustly expressed in this setting and reduce the possibility of selecting significant genes expressed at low levels, we restricted theanalysis to features according to their mean expression levels (basemean > 50) and then focused on genes showing selectiveupregulation in PreNeu (log2FoldChange>1). Filtered elements were used to generate a protein-protein interaction network through String Database. We determined marker genes by identifying a main subnetwork showing higher degreeof connectivity between nodes. We then selected 11 genes within this cluster based on their biological function related tothe regulation of immune system processes.