Identification and characterization of anthracycline-induced chemoimmunomodulation in triple-negative breast cancer

Triple-negative breast cancer (TNBC) is the deadliest breast cancer subtype with a median survival < 24 months in advanced cases. While TNBC treatment has advanced, chemotherapy remains a cornerstone of curative treatment. Despite its central role in TNBC treatment, the molecular drivers underpinning the immunomodulatory effects of chemotherapy (chemoimmunomodulation; CIM), which enable long-term efficacy and synergy with other therapeutic modalities, remain understudied, thereby limiting efforts to optimize chemotherapeutic regimens. This study aims to identify and characterize anthracycline-induced CIM states to address this challenge. To achieve this, we applied our Chemoimmunomodulation Induction Classifier (CIMIC) pipeline to delta gene expression values (ΔGE; Δlog 2 (TPM+1)) derived from TNBC cell lines (N = 6; 3 biological replicates each) profiled by bulk RNA-sequencing pre- and post-48-hour IC 30 doxorubicin exposure. CIMIC is an iterative unsupervised clustering pipeline that classifies samples into distinct groups based on their induction of 3,100 genes spanning 19 CIM pathways. Using CIMIC, we identified two discrete CIM trajectories in our samples, a functional one (Fun-CIM; N = 3 cell lines) and a dysfunctional one (Dys-CIM; N = 3 cell lines) that were highly distinct (silhouette = 0.95) and stable as determined via standard clustering metrics across multiple bootstraps; proportion of ambiguous clustering = 0.00; cluster consensus score = 1). The CIM states were characterized by significantly differential induction of antitumoral inflammatory markers, including chemoattractants ( CXCL14 ) and immune mediators ( IFNB1 and IL12A/B ) (Fun-CIM vs Dys-CIM all fold change (FC) ≥ 1.3; p < 0.05; FDR < 0.15), and protumoral markers, including tumor surveillance inhibitors ( THBS1 and TGFB2 ), immunosuppressive cell population mediators ( NNMT ), and tumor-intrinsic stress adaptation signals ( EIF2A, YARS1, and XPOT ) (Dys-CIM vs Fun-CIM all FC ≥ 1.3; p < 0.05; FDR < 0.15). Overrepresentation analysis of induced genes identified an enrichment of proteostasis and mitochondrial homeostasis programs within the Dys-CIM Group and metabolic rewiring within the Fun-CIM group (FDR < 0.01). Altogether, CIMIC enabled classification of distinct CIM induction states from paired pre- and post-treatment transcriptomic analysis, allowing interrogation of CIM in a manner not captured via traditional chemoresistance studies alone, the dissection of CIM heterogeneity, and the discovery of underappreciated molecular programs that may underlie dysfunctional CIM in TNBC. To our knowledge, CIMIC is the first unsupervised pipeline specifically designed to classify CIM trajectories. Future work will focus on characterizing baseline molecular features that may influence these CIM induction states, toward the goal of informing personalized TNBC chemotherapeutic regimen.