《MEDIATORS OF INFLAMMATION》：Single-Cell Transcriptomics Reveals Biomarkers for NK Cell Dysfunction in Endometriosis-Associated Immune Dysregulation
1. 引言
子宫内膜异位症（EM）是一种慢性炎症性疾病，特征为子宫内膜样组织在子宫腔外异位生长，全球约10%的育龄妇女受其影响。该疾病常导致严重盆腔痛和不孕，显著损害生活质量并给医疗系统带来沉重负担。诊断常被延迟，症状出现后平均滞后7-12年。尽管已有广泛研究，EM的发病机制仍不清楚。越来越多证据表明免疫失调在疾病起始和进展中起关键作用。
在EM中已记录的免疫异常中，自然杀伤（NK）细胞功能障碍已成为主要病理生理机制。NK细胞是先天免疫系统的主要效应细胞，执行免疫监视并清除病毒感染的、转化的、错位或转移的细胞，从而维持组织稳态。在生理条件下，这些细胞识别并清除通过经血逆流进入腹腔的子宫内膜细胞，防止异位植入和后续增殖。然而在EM中，NK细胞损伤构成了一种复合表型，其中细胞毒性活性减弱、受体表达谱重编程、趋化反应减弱，使得患病妇女的腹腔NK细胞表现出杀伤和迁移能力显著降低，且趋化缺陷在整个月经周期持续存在。其结果是形成一个允许的免疫环境，使得本应被清除的异位子宫内膜细胞实现免疫逃逸，从而在子宫外位点存活、植入和扩张。
现有证据指向局部微环境中复杂的细胞和分子相互作用是NK细胞功能障碍的基础。雌激素介导的机制可能通过下调CD16和颗粒酶B表达减弱NK活性；此外，子宫内膜基质细胞与其他免疫细胞之间的相互作用可进一步损害NK细胞功能。然而，NK细胞失调的精确分子通路和关键基因尚未完全阐明。
本研究应用整合多组学框架研究NK细胞功能障碍的分子机制及其对EM进展的贡献。通过结合单细胞RNA测序（scRNA-seq）与批量转录组数据集，我们构建了子宫内膜异位病变的细胞图谱，并鉴定了与免疫失调相关的分子特征。通过体外功能实验，我们进一步阐明了ENTPD1在促进子宫内膜基质细胞迁移和侵袭中的作用，从而为EM发病机制提供机制见解，并提出潜在治疗靶点以实现早期干预和个性化治疗策略。
2. 方法
2.1. 数据获取与预处理
本研究获取了一个自建的scRNA-seq数据集，包含三个EM患者的病变组织和三个同源正常子宫内膜样本，这些样本来自经伦理批准且获得知情同意的患者组织。样本经酶解制备单细胞悬液，细胞活性大于85%，使用10x Genomics Chromium Single Cell 3′ v3试剂盒构建文库；在Illumina NovaSeq 6000平台上测序，平均每个细胞约50,000条读数。原始数据使用CellRanger（v6.1.1）进行比对（GRCh38/hg38）、条形码解复用、UMI去重复并生成表达矩阵。
批量RNA-seq数据从GEO数据库下载，包括GSE105765（8个正常子宫内膜和9个EM）、GSE7305（10个正常和10个EM）和GSE6364（16个正常和21个EM），分别用作训练集、验证集和测试集。数据使用“affy”R包（v1.70.0）和RMA算法进行标准化和预处理。
2.2. 单细胞数据分析
质量控制标准包括：每个细胞检测到200-5000个基因；在少于3个细胞中检测到的基因被排除；线粒体基因占比小于30%；UMI大于100。后续分析使用Seurat R包（v4.0.5）进行：数据对数标准化，通过“FindVariableFeatures”筛选2000个高变特征；通过“Find Integration Anchors”和“IntegrateData”整合数据以校正批次效应；标准化后进行主成分分析（PCA），前40个主成分用于“FindNeighbors”、“FindClusters”（分辨率0.2）聚类，并通过“RunUMAP”可视化。
细胞类型基于CellMarker 2.0数据库的标记基因进行注释。通过“FindAllMarkers”识别簇特异性差异表达基因（DEGs），阈值设为：log2FC大于0.5，min.pct等于0.5，校正后p值小于0.05。提取免疫细胞亚群进行集中分析。使用clusterProfiler（v4.0.5）和org.Hs.eg.db（v3.13.0）进行基因本体（GO）和京都基因与基因组百科全书（KEGG）功能富集分析。
使用Monocle3（v1.0.0）推断拟时序轨迹：将Seurat对象转换为CellDataSet，使用“learn_graph”推断轨迹结构，应用“order_cells”按拟时序排序细胞。然后使用“graph_test”检测动态基因（q值小于0.05）。
2.3. 批量RNA分析与生物标志物筛选
使用limma（v3.48.3）在批量数据集中筛选DEGs（log2FC大于1，校正后p值小于0.05），并与scRNA-seq差异基因取交集。使用支持向量机递归特征消除（SVM-RFE；e1071 v1.7–9）和最小绝对收缩和选择算子（LASSO；glmnet v4.1–2）进行特征选择，结合10折交叉验证。通过单变量和多变量逻辑回归验证枢纽基因。基于rms（v2.0）构建列线图，pROC（v1.18.0）绘制受试者工作特征（ROC）曲线，决策曲线分析（DCA）（v2.0）进行DCA。
使用CIBERSORT（v1.03）估计免疫浸润；通过Spearman相关分析评估基因与免疫细胞的关系。使用GSVA（v1.40.1）对Hallmark和KEGG通路进行评分，使用clusterProfiler实施GSEA，标准为标准化富集分数（NES）大于1且校正后p值小于0.05。基于STRING（v11.5）构建蛋白质相互作用网络（PPI），在Cytoscape（v3.9.1）中可视化。通过比较毒理基因组学数据库（CTD）筛选潜在药物。
2.4. 细胞培养与转染
人子宫基质细胞（HUSCs）购自尚恩生物技术（武汉，中国；货号SNP-H371）。永生化人子宫内膜基质细胞（ihESCs，模拟EM样行为）购自伊默生物技术（厦门，中国；货号IM-H540）。两种细胞系均经供应商通过短串联重复谱分析认证，并提供确认身份和无支原体状态的分析证书。细胞在接收后10代内使用。细胞在完全子宫内膜基质细胞培养基（由各自供应商提供）中培养，添加10%胎牛血清、100 U/mL青霉素和100 μg/mL链霉素，在37°C、5% CO2下孵育。
过表达（OE）实验使用Lipofectamine 3000（Invitrogen）转染pcDNA3.1-ENTPD1或空载体（OE-NC）；敲低实验使用靶向ENTPD1的siRNA（si-ENTPD1）或阴性对照（si-NC）（终浓度50 nM；序列：si-ENTPD1）。ihESCs的转染效率通常大于80%，HUSCs为70%–85%。转染后48小时通过qPCR验证转染效率。
2.5. 实时定量PCR（qRT-PCR）
使用TRIzol（Invitrogen）提取总RNA，使用PrimeScript RT试剂盒（Takara）进行反转录；在QuantStudio 5系统（Applied Biosystems）上使用SYBR Green Master Mix进行定量。引物如下：ENTPD1：正向5′-AGGTGTGGATATCAGCCTGTA-3′，反向5′-CTTCTCTCCGAGATCCCTTCC-3′；GNLY：正向5′-ATATCGTCTCCCAGATGCACT-3′，反向5′-AGCTTTCCTCTCCAAGTTGAT-3′；PRF1：正向5′-CCACTCACAGGCAGCCAA-3′，反向5′-GGAGATGAGCCTGTGGTAAG-3′；GAPDH：正向5′-GAAGGTGAAGGTCGGGAGTC-3′，反向5′-GAAGATGGTGATGGGATTTC-3′。
2.6. Transwell侵袭实验
将转染后的ihESCs以5×104个细胞接种于Matrigel包被的Transwell（孔径8 μm，Corning）上室，加入无血清培养基；下室加入含20% FBS的培养基。孵育24小时后，固定浸润细胞，用0.1%结晶紫染色；在显微镜（200×）下随机拍摄五个视野并计数。实验独立重复三次。
2.7. 划痕愈合实验
将转染后的ihESCs培养于6孔板至90%融合，用无菌枪头划痕形成“伤口”，继续在无血清培养基中培养；在0小时和24小时使用倒置显微镜成像。使用ImageJ（v1.53）测量划痕宽度，计算伤口愈合率：愈合率（%）=（初始宽度-24小时宽度）/初始宽度×100。实验独立重复三次。
2.8. 统计分析
使用R v4.3.2和GraphPad Prism 9.0进行统计分析。使用Wilcoxon秩和检验与Benjamini–Hochberg错误发现率（FDR）校正进行差异表达分析；调整后p值小于0.05且|log2FC|大于0.5视为显著。对于涉及多组的体外实验，使用单因素方差分析（ANOVA）后接Tukey事后检验；对于两组比较，使用未配对双尾Student t检验或Mann–Whitney U检验。功能实验中的多重比较使用Bonferroni校正进行校正。所有报告的p值在进行了多重检验的地方均经过调整，调整后p值小于0.05被认为具有统计学意义。
3. 结果
3.1. EM中单细胞转录组景观的构建与表征
为描绘EM的细胞景观，我们对三名患者的异位病变组织和三名对照个体的配对在位子宫内膜进行了scRNA-seq。经过严格的质量控制过滤，保留了线粒体基因含量与总转录本计数呈负相关、同时检测到的基因数与UMI计数呈正相关的细胞。在鉴定出的2000个高变基因中，观察到细胞毒性和促炎介质（GNLY、IL1B、CXCL8、S100A9和LYZ）的显著表达。
PCA显示子宫内膜异位症样本和对照样本之间清晰分离，肘部图指导选择前20个主成分进行下游分析。在分辨率为0.8的无监督聚类中，识别出11个转录不同的簇（簇0-11），在UMAP上可视化，并基于已建立的典型标记和SingleR参考映射注释为八种主要细胞类型：成纤维细胞（簇0、1和3）、上皮细胞（簇2）、内皮细胞（簇4）、T细胞（簇5）、NK细胞（簇6）、巨噬细胞（簇7）、平滑肌细胞（簇8）和组织干细胞（簇9），次要簇10-11代表低质量或过渡群体。
EM和对照样本之间的细胞组成存在明显差异。跨细胞类型的差异表达分析突出显示免疫亚群中炎症和细胞毒性分子（SLPI、WFDC2、IL1B、CXCL8、S100A9、LYZ和GNLY）上调，而成纤维细胞和平滑肌区室中细胞外基质和Wnt通路调节因子（COL1A1、COL1A2、DCN、SFRP2和SFRP4）下调。这些发现共同建立了一个全面的单细胞图谱，揭示了子宫内膜异位病变中深刻的免疫激活和基质重塑。
3.2. 已鉴定细胞群的分子特征和功能富集
通过标记基因的空间分布验证了细胞类型注释的准确性。上皮细胞特异性表达EPCAM，而内皮细胞富含血管标记。平滑肌细胞以TAGLN和ACTA2表达为特征；成纤维细胞高表达COL1A1、COL1A2和DCN；NK/T细胞表达细胞毒性标记NKG7、GNLY、GZMB和PRF1。单核细胞显示独特的LYZ表达。GO和通路分析显示在免疫反应、粒细胞趋化性、细胞因子-受体相互作用、白细胞分化和ECM组织方面显著富集。
3.3. 轨迹推断揭示疾病相关的细胞动力学
细胞的拟时序排序揭示了疾病微环境中的连续状态转换。沿推断轨迹，免疫效应分子逐渐增加，而ECM/纤维化相关基因（COL1A1、COL1A2、DCN和SPARCL1）逐渐减少，表明病变形成过程中存在“免疫激活-基质重塑”连续体。单细胞差异基因与轨迹相关标记的整合鉴定出361个关键基因，其细胞类型特异性表达模式如图所示。
3.4. 跨平台验证鉴定稳健的疾病标志物
为验证我们的单细胞发现，我们分析了两个独立的批量RNA-seq队列（GSE105765和GSE7305），鉴定了大量DEGs。与单细胞衍生的关键基因取交集，得到20个一致失调的候选基因，在数据集中显示可重复的表达模式。蛋白质-蛋白质相互作用网络分析将这些基因定位为免疫调节和细胞因子信号传导中的中枢枢纽。
3.5. 细胞毒性相关基因显示出强大的判别能力
单变量分析鉴定出NKG7、GZMB、SCGB2A1、EPCAM、ENTPD1、GNLY、PRF1、PTGIS、GATA6和STAR与疾病状态显著相关。这些基因在子宫内膜异位组织和对照组织之间表现出明显的表达差异，单基因ROC曲线显示出优异的判别能力（大多数基因AUC大于等于0.80）。
3.6. 基于机器学习的特征选择建立三基因诊断特征
SVM-RFE鉴定出一个最优的8基因子集，达到约90%的交叉验证准确率。LASSO回归进一步细化特征集，揭示了稳定的系数路径，并鉴定出权重最高的基因。基于这些分析，我们开发了一个包含GNLY、PRF1和ENTPD1的列线图，用于EM风险评估。
3.7. 构建并验证用于EM的诊断列线图
为提供EM风险预测的定量工具，我们使用训练队列（GSE105765）开发了一个包含三个枢纽基因（GNLY、PRF1和ENTPD1）的列线图。在列线图中，GNLY和PRF1（反映受损的细胞毒性活性）表达较低以及ENTPD1表达较高导致总分更高和预测的EM概率更大。
DCA证明，与全治疗或无治疗策略相比，列线图在广泛的阈值概率范围（10%–90%）内赋予更优的临床净收益，证实了其临床效用。校准曲线显示，在训练队列（GSE105765）、内部验证队列（GSE7305）和独立外部测试队列（GSE6364）中，列线图预测概率与实际EM发生之间具有极好的一致性，偏差校正线与理想参考线紧密对齐。
ROC分析进一步证实了强大的判别性能，在训练队列中AUC值为0.837（95% CI 0.783–0.892），在验证队列（GSE7305）中为0.671（95% CI 0.595–0.746），在测试队列（GSE6364）中为0.777（95% CI 0.705–0.835）。这些结果共同验证了三基因列线图作为EM可靠且临床适用的诊断工具。
3.8. 免疫浸润景观与特征基因相关
对整合的批量RNA-seq队列应用CIBERSORT反卷积。该分析揭示了免疫微环境的深刻重塑，EM病变与正常子宫内膜之间多个免疫细胞亚群的比例存在显著差异。值得注意的是，EM样本显示静息NK细胞（p小于0.001）和活化的CD8+T细胞显著减少，同时M2巨噬细胞（p小于0.01）、静息肥大细胞和中性粒细胞比例增加，这与有利于病变持续存在的免疫抑制和促纤维化环境一致。对从scRNA-seq鉴定的主要细胞类型进行免疫相关通路特征的GSVA进一步突出了EM中细胞毒性活性受抑制。与对照组相比，与EM相关的簇中细胞毒性效应程序，包括“活化的NK细胞”、“CD8+T细胞记忆/活化”和“γδ T细胞”的通路评分显著降低，而与替代性巨噬细胞极化和调节性T细胞活性相关的特征相对富集。
Spearman相关分析揭示了三个枢纽基因的不同免疫浸润模式。GNLY和PRF1与活化的NK细胞（r约0.70–0.78）和活化/记忆CD8+T细胞（r=0.45–0.65）呈强正相关，但与M2巨噬细胞和调节性T细胞呈负相关。在批量样本中，两个基因均与NK细胞浸润强相关（GNLY：r=0.78，p=5.8×10?42；PRF1：r=0.73，p=1.5×10?47）。相比之下，ENTPD1与静息NK细胞（r=0.781，p小于0.001）、单核细胞、嗜酸性粒细胞和M2巨噬细胞呈最强正相关，同时与活化的细胞毒性淋巴细胞呈较弱/负相关，并且与CD8+T细胞浸润相关性最高（r=0.81，p=6.7×10?48）。因此，GNLY和PRF1主要标记EM中受损的细胞毒性活性，而ENTPD1尽管主要在基质中表达，但与CD8+T细胞浸润和抑制性髓系群体相关，可能通过腺苷介导的免疫调节促进NK细胞功能障碍。
3.9. 单细胞分辨率映射分析免疫细胞特异性表达
为验证从批量RNA-seq分析中鉴定的枢纽基因的免疫细胞特异性表达模式，我们使用公开可用的EM及相关组织（GSE213216）数据集进行了独立的scRNA-seq分析，该数据集包括子宫内膜异位囊肿、异位和在位子宫内膜、未受累卵巢和无EM腹膜。如图所示，UMAP降维显示18个不同的细胞簇，成功注释为疾病微环境内的主要细胞群体，包括上皮细胞、成纤维细胞、内皮细胞、巨噬细胞、单核细胞、CD4+T细胞、CD8+T细胞、B细胞和NK细胞。细胞组成分析表明成纤维细胞是最丰富的细胞类型（26%），其次是B细胞（15%）、单核细胞（13%）、NK细胞（13%）、内皮细胞（12%）、CD8+T细胞（10%）、巨噬细胞（8%）和CD4+T细胞（3%），突出了NK细胞在EM相关微环境中相对较高的浸润和潜在功能重要性。
小提琴图清晰显示了枢纽基因的细胞类型特异性表达。ENTPD1（CD39）主要在髓系亚群（巨噬细胞和单核细胞）和CD8+T细胞中表达，在NK细胞和B细胞中表达极少。相比之下，GNLY和PRF1在NK细胞和CD8+T细胞中表达最高，与其在细胞毒性功能中的既定作用一致。点图进一步量化了这些基因在免疫亚群中的平均表达和表达细胞的百分比。值得注意的是，NK细胞显示GNLY和PRF1的平均表达最高，并且在细胞毒性效应细胞中表达细胞分布最广，尽管其丰度中等。ENTPD1主要在巨噬细胞、单核细胞和内皮细胞中显示强表达，支持其在免疫抑制性腺苷信号传导中的作用。这些单细胞分辨率发现独立证实了枢纽基因的免疫细胞特异性表达模式，强调了NK细胞群在EM微环境中的主导细胞毒性潜力以及髓系和内皮区室的免疫抑制贡献。
3.10. 通路富集和药物重定位机会
Hallmark和KEGG富集分析，结合GSEA，显示与我们的特征共表达的基因在干扰素反应、IL2-STAT5信号传导、TNFα-NFκB信号传导、炎症反应、缺氧、补体级联和Th1/Th2分化中显著富集，而氧化磷酸化被抑制。药物-基因相互作用网络分析鉴定出多种可能靶向ENTPD1/GNLY/PRF1的临床可用化合物，包括白藜芦醇、布洛芬、吲哚美辛、罗非考昔、达那唑、米非司酮和孕酮，提示治疗重定位机会。
3.11. ENTPD1促进子宫内膜基质细胞迁移
为功能验证我们的计算发现，我们首先比较了原代HUSCs和ihESCs中三个特征基因的转录水平。虽然GNLY和PRF1在细胞类型之间无显著差异，但ENTPD1在ihESCs中显著升高（2.0倍，p小于0.001），与我们的批量和单细胞分析一致。
细胞功能实验证明，ENTPD1过表达（OE-ENTPD1）显著增强伤口愈合和细胞迁移，而ENTPD1敲低（si-ENTPD1）产生相反效果。具体而言，与对照（OE-NC，p小于0.01）相比，OE-ENTPD1使划痕闭合增加30%–40%，迁移细胞数增加1.7倍（p小于0.001）。相反，si-ENTPD1使划痕闭合相对于si-NC减少30%–40%（p小于0.001）并减少迁移细胞数（p小于0.01）。Transwell迁移/侵袭实验证实了这些发现：ENTPD1过表达显著增加跨膜细胞数，而沉默显著抑制迁移（OE-ENTPD1 p小于0.01，si-ENTPD1 p小于0.001）。
总之，ENTPD1在子宫内膜基质细胞中显示出明确的促迁移和伤口愈合功能，而GNLY和PRF1在基质细胞中未显示显著变化，支持其在细胞毒性免疫细胞中的细胞类型特异性表达。这些实验验证，结合我们的多队列、多尺度分析，确立了ENTPD1作为EM发病机制中的关键功能分子和有前景的治疗靶点。
4. 讨论
EM是一种多因素妇科疾病，特征为子宫内膜样组织异位生长，伴有慢性盆腔痛和不孕。除了经血逆流，汇聚的证据指出免疫功能障碍——特别是NK细胞毒性减弱——是病变持续和进展的驱动因素。多项研究报告患者外周血和腹腔液中NK细胞毒性降低，促进了免疫逃逸和异位子宫内膜细胞在腹腔内的持续增殖。从机制上讲，功能受损涉及激活性和抑制性NK受体表达失衡以及免疫抑制性细胞因子环境。
本研究应用多组学框架描绘EM中NK功能障碍及其与基质生物学的相互作用。scRNA-seq生成了子宫内膜异位病变的高分辨率细胞图谱，与批量转录组整合后产生了三基因诊断特征（GNLY、PRF1和ENTPD1）。功能扰动确立了ENTPD1在子宫内膜基质细胞迁移和侵袭中的效应作用，将免疫-基质相互作用定位为疾病进展的核心过程。这些结果勾勒出免疫逃逸和局部微环境重塑作为相互关联的轴，对早期干预和个体化治疗具有转化意义。
4.1. 细胞毒性免疫功能障碍的分子基础：GNLY和PRF1
GNLY和PRF1是NK细胞和细胞毒性T淋巴细胞的核心细胞毒性效应物，通过孔形成和颗粒酶递送介导靶细胞死亡。在单细胞分辨率下，GNLY和PRF1主要定位于NK和CD8+T细胞区室，它们的表达与病变中活化的细胞毒性群体相关，强调了这些细胞在EM免疫监视中的贡献。GNLY和PRF1在基质细胞中与对照相比无显著差异，表明基质病理不是由这些细胞毒性效应物的直接调节驱动的。免疫细胞内细胞毒性基因程序减少与病变细胞的免疫逃逸和持续的腹腔内扩张一致。鉴于它们作为NK活性指标的既定作用，病变中GNLY和PRF1减少提供了NK杀伤功能受损的分子证据。这种损伤与免疫抑制介质升高以及抑制性NK受体表达增加同时发生，与NK耗竭和功能下降一致。
GNLY和PRF1的细胞分辨率改变与先前观察到的EM患者外周血和腹腔液中NK细胞毒性减弱一致。受体水平的变化进一步支持这一图景：EM中已记录到NKG2D减少和NKG2A增加。NKG2D识别应激诱导的配体，受体表达减少和可溶性配体升高损害靶标识别和杀伤。通过NKG2A–HLA-E和LILRB1–HLA-G的抑制信号在病变中增强，抑制NK效应功能。总之，这些失衡在病变微环境内建立了一个功能失调的NK状态。NK细胞上检查点分子增加的报道为EM中细胞毒性功能障碍提供了额外背景。
4.2. ENTPD1作为双功能介质
ENTPD1（CD39）是一种外切核苷酸酶，将细胞外ATP/ADP水解为AMP，后者是腺苷的前体，从而有助于免疫抑制信号传导。在肿瘤中，CD39/CD73-腺苷轴通过腺苷受体接合和下游细胞毒性反应抑制来减弱抗肿瘤免疫。关于EM中外切核苷酸酶异常表达的报道表明了嘌呤能代谢通路的参与，为ENTPD1介导的免疫抑制提供了背景。在本研究中，ENTPD1在基质细胞中表达升高， coupled with其功能上促进基质细胞侵袭和迁移，表明它既驱动病变行为又塑造免疫环境。腺苷信号通过腺苷受体减弱NK脱颗粒和细胞因子产生，降低细胞毒性输出。这些数据支持一个双功能模型，其中ENTPD1调节基质侵袭性并在动态、相互作用的微环境内强制执行局部免疫抑制。
细胞因子介导的抑制，特别是涉及TGF-β和IL-10，构成了EM免疫病理的一个充分描述的维度。腹腔液中IL-6和IL-10升高抑制NK细胞毒性并使NK表型偏向免疫调节。在此背景下，腺苷能通路提供了免疫抑制的互补途径；ENTPD1轴和细胞因子信号可以协同运作，维持一个有利于异位细胞存活和增殖的免疫许可生态位。在治疗上，CD39/CD73-腺苷通路在肿瘤学中正在积极研究，为EM的机制引导干预提供了概念模板。
在诊断界面，三基因特征在独立队列中表现出稳健性能。当代护理仍常依赖腹腔镜检查且诊断时间漫长
