單細胞轉錄組測序揭示化療增強乳腺癌細胞亞精胺代謝激活腫瘤免疫

劉子菲，黃鵬翰，江詠雪，曾文鋒

乳腺癌發病率日益增加，世界衛生組織國際癌癥研究機構2020 年癌癥報告顯示，乳腺癌超越了肺癌成為世界上最普遍的惡性腫瘤，并位于癌癥相關死亡人數的前5 位［1］。乳腺癌是全球包括中國在內的多個國家和地區女性發病率最高的惡性腫瘤。在中國，乳腺癌發病率約占中國女性惡性腫瘤總發病率的15%，而且發病率和死亡人數呈逐年上升趨勢［2，3］。2020 年世界衛生組織公布的我國癌癥數據顯示，乳腺癌是我國第四大女性癌癥死亡原因，僅次于肺癌、結直腸癌和胃癌［3］。而且乳腺癌年輕化的問題在我國尤為突出，在小于45 歲的中國女性人群中乳腺癌已經成為首要的癌癥死亡原因［2］。目前乳腺癌的治療主要依靠外科手術、放療、化療、內分泌治療和靶向治療。其中新輔助化療是乳腺癌患者的重要根治性手段之一，研究其抑制腫瘤進展的作用和機制有助于我們加深腫瘤化療的理解和臨床上合理地選擇化療藥物。近年來的研究表明，化療除了對腫瘤本身的殺傷作用外，還可能還可以調節腫瘤微環境代謝改變［4-6］及激活抗腫瘤免疫［7］。

研究表明，惡性腫瘤中多胺生物合成酶的表達增加，并且在癌細胞內出現精胺和亞精胺蓄積?；熌艽龠M癌細胞多胺代謝通路激活，但該癌細胞代謝通路對微環境的影響還未被闡明。傳統觀點認為，精胺和亞精胺的富集與惡性腫瘤免疫抑制相關［8-10］，但隨著研究的深入，近年來有研究者發現多胺代謝的亢進是引起自身免疫性疾病的原因之一，例如銀屑?。?1］。此外，補充亞精胺可以導致T 淋巴細胞線粒體活性增強［12］并增強抗腫瘤免疫，而且亞精胺可以作為疫苗的佐劑刺激免疫反應［13］。除了適應性免疫外，亞精胺對固有免疫細胞例如自然殺傷細胞的作用也被逐漸揭示［14］。亞精胺代謝促進腫瘤免疫存在爭議，乳腺癌亞精胺代謝與腫瘤免疫之間的關系及其調控機制尚不明確?；熥鳛槟[瘤多胺代謝通路的激活手段之一，其是否有參與調控乳腺癌亞精胺代謝從而影響腫瘤免疫也有待進一步探討。

1 資料與方法

1.1 數據資料收集

于GEO 基因表達數據庫（https：//www.ncbi.nlm.nih.gov/geo/）檢索新輔助化療乳腺癌單細胞測序相關數據信息。經檢索，符合的數據集為GSE189302［15］，該數據來源于自發成瘤鼠腫瘤單細胞測序。為擴大分析范圍，我們收入了人乳腺癌單細胞測序數據集GSE176078［16］進行數據驗證。

1.2 單細胞測序數據分析

利用R 包Seurat［17，18］進行單細胞數據預處理，兩個數據集的細胞注釋內容都依據原文或在線數據提供的注釋文獻，并從測序數據中篩選出上皮細胞進行進一步聚類。利用R 包Harmony［17］對數據進行去批次處理。

1.3 單細胞數據基因表達可視化

利用R 包Nebulosa［19］以及Seurat 中的VlnPlot、FeaturePlot 函數對單細胞數據中亞精胺代謝關鍵基因的表達水平進行可視化處理。

1.4 單細胞樣本分組與差異基因分析

在單細胞測序數據分析的基礎上，利用R 包Seurat 中的AddMouduleScore 函數將GOBP_SPERMIDINE_METABOLIC_PROCESS（https：//www.gseamsigdb.org/gsea/msigdb/human/geneset/GOBP_SPER -MIDINE_METABOLIC_PROCESS.html）數據集中的基因合并添加到GSE176078 單細胞測序數據中，根據該數據的加權平均數將患者分為癌細胞亞精胺代謝高（high：CID45171；CID4530N；CID4535；CID4523；CID4066；CID3921；CID3948；CID44991；CID3941；CID4513）與亞精胺代謝低（low：CID4471；CID4290A；CID3963；CID4067；CID4465；CID44971；CID4515；CID4461；CID4495；CID4463）兩組。利用R 包Seurat 中的FindMarkers 函數分析癌細胞亞精胺代謝高低樣本中各類淋巴細胞的差異表達基因。利用R包ggplot2（H.Wickham.ggplot2：Elegant Graphics for Data Analysis.Springer-Verlag New York，2016。）可視化差異表達基因并繪制火山圖。并運用ggrepel 包顯示出P<0.05 且log2FC>0.5 的差異表達基因的基因名稱。

1.5 功能富集分析

利用R 包ClusterProfiler［20］以獲得差異基因在細胞內的生物學過程（Biological Process，BP）情況，α 設為0.05，P 值校正方法選擇Benjamini &Hochberg 法，即對每個P-value 做校正，轉換為qvalue。q=P*n/rank，其中rank 是指P-value 從小到大排序后的次序，獲得基因本體論（GO）分析結果。

1.6 統計學方法

采用R4.1.0軟件和Rggpubr包進行分析。兩兩比較采用t檢驗，以P<0.05 為具有統計學差異。

2 結 果

2.1 化療導致乳腺癌細胞亞群分布改變

為了探究化療是否能夠影響乳腺癌細胞的亞精胺代謝，對GSE189302 數據集進行分析。GSE189302 將MMTV-Wnt1 腫瘤細胞注射到野生型C57BL6/J 小鼠的乳房脂肪墊使其成瘤，待腫瘤體積達到~200 mm3，化療組每日腹腔注射多柔比星4 mg/kg 共5 天，對照組不做處理，末次給藥48 小時后分別獲取腫瘤組織進行單細胞轉錄組測序［15］。該數據集包括1 例化療組（GSM5698657）及1 例對照組（樣本GSM5698658）的單細胞轉錄組測序結果，將兩者合并進行UMAP 聚類分析，得到包括癌細胞及免疫細胞的共10個細胞群，其中3群細胞根據上皮細胞標記基因注釋為癌細胞（圖1A）。提取癌細胞進行再聚類，獲得9 個細胞亞群（圖1B）。如圖1C、D 所示，化療組與對照組的癌細胞分布存在明顯差異，化療后癌細胞0、1、7 群減少，2-6 群增加，其中第2、5 群細胞呈現“從無到有”的顯著改變。

圖1 GSE189302 化療前后乳腺癌細胞聚類分析 A：腫瘤組織所有細胞UMAP 聚類散點圖；B：癌細胞的UMAP 聚類散點圖；C：化療組及對照組癌細胞的UMAP 聚類散點圖；D：兩組癌細胞的亞群堆疊圖

2.2 化療誘導乳腺癌亞精胺代謝旺盛細胞亞群形成

為了進一步探究化療前后癌細胞的改變，對化療組和對照組的基因表達做差異分析，從而獲得兩組差異基因，以差異倍數>0.5，校正P值<0.05作為篩選閾值。詳見圖2A 火山圖，其中化療組上調基因有342 個，下調基因有406 個。分析差異基因發現，亞精胺代謝重要限速酶鳥氨酸脫羧酶1（OCD1）［21］表達水平上調，繪制基因表達水平小提琴圖如圖2B 所示；亞精胺代謝的另一重要分子精胺氧化酶（SMOX）雖然未達到篩選閾值，但其基因表達水平在兩組之間仍存在顯著差異（圖2C），這些差異提示化療后癌細胞亞精胺代謝水平升高。為了探究此改變是否與前述癌細胞亞群分布變化有關，在癌細胞聚類散點圖中分別顯示OCD1 及SMOX 基因表達水平，以及根據亞精胺代謝數據集對各樣本的癌細胞進行綜合評分（圖2D~G），對照圖1 中細胞亞群的位置，可以看到第5 群癌細胞具有OCD1、SMOX 高表達，并且亞精胺代謝水平顯著升高。以上數據表明多柔比星化療誘導乳腺癌細胞中亞精胺代謝旺盛亞群形成。

2.3 乳腺癌患者腫瘤細胞亞精胺代謝水平存在差異

為了探究人類乳腺癌細胞亞精胺代謝情況，對GSE176078 數據集進行相關分析。該數據集包含26 例未經化療的原發性乳腺癌組織的單細胞測序結果，包括11 例ER 陽性、5 例HER2 陽性以及10例三陰型乳腺癌，并提供細胞注釋信息［16］。我們將其中6 例未注釋癌細胞群的樣本予以剔除，其余20 例樣本的癌細胞UMAP 聚類結果如圖3A、B所示。根據亞精胺代謝數據集對各樣本的癌細胞進行綜合評分，箱式圖3D 顯示各樣本癌細胞亞精胺代謝水平存在差異。進而以評分中位值作為cut-off 值，將樣本分為亞精胺代謝high、亞精胺代謝low 兩組（圖3C、E）。20 例樣本的臨床特征及分組信息如表1、表2 所示。

表2 GSE176078 亞精胺代謝high、low 兩組的臨床特征

圖3 GSE176078 乳腺癌細胞亞精胺代謝水平 A：各樣本癌細胞的UMAP 聚類散點圖；B：各分子亞型癌細胞的UMAP 聚類散點圖；C：亞精胺代謝high、low 兩組癌細胞的UMAP 聚類散點圖；D：各樣本亞精胺代謝水平；E：兩組的亞精胺代謝水平

2.4 乳腺癌亞精胺代謝與腫瘤浸潤T 細胞的關系

近期有研究報道，亞精胺能夠通過增加脂肪氧化磷酸化增強CD8+T 細胞的線粒體活性和細胞毒性功能［12］。我們猜測癌細胞代謝產生的亞精胺能夠作用于腫瘤浸潤T 細胞，影響抗腫瘤免疫。為了驗證假說，根據GSE176078 所提供的細胞注釋信息，分析亞精胺代謝high、low 兩組的免疫細胞浸潤情況（表3）。如圖4A、B，亞精胺代謝high的腫瘤具有相近T 細胞浸潤比例，并且T 細胞各亞群的比例差異不大。GSE176078數據集的原作者根據基因表達特點對CD4+、CD8+T 淋巴細胞進行了再聚類，并將其注釋為具有不同功能特點的4 個CD4亞群和5個CD8亞群［16］。比較發現亞精胺代謝high、low 兩組之間亞群分布存在差異（圖4C、D）。尤其值得注意的是，在亞精胺代謝high 的腫瘤中，高表達T 細胞激活及效應相關基因ZFP36［22］的CD8+亞群（c4）明顯增多，以抑制性免疫檢查點LAG3為特點的CD8+亞群（c8）、以衰老標志物GZMK［23］為特點的CD8+亞群（c5）明顯減少（圖4C），提示乳腺癌細胞亞精胺代謝水平對CD8+T 細胞的基因表達和功能存在影響。

表3 GSE176078 亞精胺代謝high、low 兩組的免疫浸潤情況

圖4 GSE176078 亞精胺代謝high、low 兩組的免疫浸潤情況 A：兩組腫瘤組織細胞分類堆疊圖；B：T 細胞亞群堆疊圖；C：CD8+T 細胞亞群堆疊圖；D：CD4+T 細胞亞群堆疊圖

2.5 乳腺癌亞精胺代謝對CD8+T 細胞功能的影響

為了進一步驗證上述發現，我們對CD8+T 細胞做進一步分析。首先對亞精胺代謝high、low 兩組CD8+T 細胞的基因表達做差異分析，從而獲得兩組差異基因，以差異倍數>0.5，校正P值<0.05作為篩選閾值。如圖5A 火山圖所示，亞精胺代謝high 組有58 個基因上調，59 個基因下調，我們將上調基因予以標注，其中包括與T 細胞激活或效應功能相關的基因CD69、DUSP1、FOS、CXCR4［22］，基因表達水平如圖5C-F 所示。LAG3、PDCD1 是重要的免疫檢查點，能夠抑制免疫微環境［12，24］，其在亞精胺代謝high 的腫瘤浸潤CD8+T 細胞中表達量顯著降低（圖5G、H）。上述基因表達差異提示癌細胞的亞精胺代謝水平與腫瘤浸潤CD8+T 細胞的激活和功能呈正相關。利用基因本體（gene ontology，GO）數據庫富集亞精胺代謝high 組上調基因，發現與抗原處理與呈遞（antigen process and presentation）、抗原受體介導的信號通路（antigen receptor-mediated signaling pathway）、T 細胞受體信號通路（T cell recepter signaling pathway）、IFN-γ 介導的信號通路（interferon-gamma-mediated signaling pathway）等適應性免疫功能相關通路有關，進一步證明癌細胞的亞精胺代謝增強腫瘤浸潤CD8+T 細胞功能的作用。

圖5 GSE176078 亞精胺代謝high、low 兩組CD8+T 細胞的基因表達情況 A：差異基因火山圖；B：基因本體（GO）富集分析；C~H：相關基因表達水平差異

3 討 論

乳腺癌已經成為全球發病率最高的惡性腫瘤之一。新輔助化療為乳腺癌患者的治療帶來了更多的選擇和提高了患者預后。該術前化療手段的臨床獲益也引發了研究人員思考除了直接殺傷腫瘤細胞外，化療是否能重塑腫瘤微環境以增強抗腫瘤免疫。一直以來多胺，包括腐胺、精胺、亞精胺都被認為與免疫抑制形成相關［21］。但近幾年，多胺對免疫更復雜的作用逐漸被揭示。例如對于抗腫瘤免疫中起關鍵作用的CD8+T 細胞，諾貝爾獎獲得者Tasuku Honjo 團隊［12］發現在衰老的小鼠體內補充亞精胺增強CD8+T 細胞的代謝能力產生更多ATP，外源補充亞精胺能實現抑制免疫細胞衰老并促進抗腫瘤免疫。然而Hibino 等［9］人的研究顯示，腫瘤細胞產生的內源性亞精胺能夠通過抑制TCR 聚集從而阻止其下游信號傳導，多胺合成抑制劑能夠增強CD8+T 細胞的功能。上述多胺對CD8+T 細胞功能影響的矛盾揭示了多胺代謝通過多種機制從不同方向上影響抗腫瘤免疫的復雜性，無論是外源補充多胺或是體內多胺代謝對于抗腫瘤免疫的調控都需要更多探索。而我們的研究揭示了，亞精胺代謝的激活與乳腺癌腫瘤免疫激活存在正相關，豐富了我們對多胺代謝與腫瘤免疫的理解。

既往分析腫瘤代謝與免疫微環境的串擾，因技術的局限只能從整個腫物中提取核酸或代謝物進行檢測，隨著單細胞測序技術的普及，我們可以分析癌細胞代謝通路變化對腫瘤浸潤免疫細胞的調控。本研究發現化療可以增強乳腺癌細胞的多胺代謝，并且與腫瘤微環境中浸潤的淋巴細胞激活相關。我們猜想化療誘導腫瘤死亡的過程中，可能引起了癌細胞多胺代謝增強并釋放到環境引起淋巴細胞激活。本研究為開展后續體外實驗提供生物信息學依據。受限于標本數量和腫瘤類型，我們的發現尚未做到多癌種分析。后續的研究我們將收集更多的人乳腺癌及其他癌種化療前后的單細胞測序數據，并進一步開展功能實驗?？傊?，本研究證明了化療能促進乳腺癌細胞亞精胺代謝增強，同時癌細胞亞精胺代謝強與效應性淋巴細胞的浸潤和激活呈現正相關。本研究的發現為乳腺癌新輔助治療引起的癌細胞免疫原性死亡原理提供了理論支持。