N6-甲基腺嘌呤调控因子对宫颈癌预后及免疫浸润的作用①

郭依琳 白洋洋 王 璐 徐 臻 王习亮 王武亮 （郑州大学第二附属医院妇产科，郑州 450014）

宫颈癌是女性生殖系统最常见的恶性肿瘤，其 中90%发生于发展中国家［1］。2020 年，全球宫颈癌新发病例604 127 例，死亡病例341 831 例［2］。随着宫颈癌筛查系统逐步开展及宫颈癌疫苗普及，早期宫颈癌预后较好，但仍有部分宫颈癌患者首诊时已为晚期，失去了最佳治疗时机。因此阐明宫颈癌发生发展的分子机制对改善宫颈癌患者预后具有重要临床意义。

m6A 是真核生物最丰富的编码和非编码RNA甲基化修饰［3］。m6A 是位于腺嘌呤碱基第6 位氮原子的甲基化修饰，主要分布于终止密码子附近和3′-非翻译区，定位于保守的RRACH 序列（R=G 或A；H=A，C或U），在RNA加工成熟、剪接和转录调控中发挥重要作用［4-5］。m6A 受甲基转移酶复合物（writers）、去甲基化酶（erasers）和m6A 读取蛋白（readers）调控［6］。最近研究提出m6A 调控因子具有免疫调节功能，参与肿瘤免疫微环境（tumor immune microenvironment，TIME）的免疫调控，但其表达对宫颈癌免疫浸润和预后的研究未见报道。因此，本研究利用TCGA 和GEO 数据库分析m6A 调控因子在宫颈癌中的表达，建立LASSO Cox回归预后模型，利用ssGSEA 计算免疫细胞浸润程度，为宫颈癌患者临床诊疗和预后改善提供依据。

1 资料与方法

1.1 TCGA 和GEO 数 据 获 取 从TCGA 数 据 库（http：//cancergenome.nih.gov/）下载306例宫颈癌组织和3 例正常组织RNA 测序转录组数据（FPKM 格式）及相应临床信息。从GEO 数据库（https：//www.ncbi.nlm.nih.gov/geo/）下载55 例宫颈癌组织和17 例正常组织（GSE52903）芯片转录组数据及相应临床信息。采用R 软件（4.1.2）“limma”包将TCGA数据集中FPKM 格式的RNA 测序转录组数据转换为TPM 格式，并与GSE52903 数据集合并，对合并数据进行矫正。

1.2 m6A 调控因子选择及相关性分析 选择25个公认的m6A 调控因子，包括8 个writers（METTL3、METTL5、METTL14、METTL16、WTAP、ZC3H13、RBM15、RBM15B）、14 个readers（YTHDC1、YTHDC2、YTHDF1、YTHDF2、YTHDF3、HNRNPC、FMR1、LRPPRC、HNRNPA2B1、IGF2BP1、IGF2BP2、IGF2BP3、ELAVL1、RBMX）和3 个erasers（FTO、ALKBH5、ALKBH3）。采用R 软件“limma”包对25 个m6A 调控因子进行Wilcox 秩和检验比较宫颈癌组织和正常组织表达差异。“pheatmap”包绘制m6A 调控因子表达和临床特征热图。“corrplot”包比较m6A 调控因子间的相关性。

1.3 LASSO 预后模型构建及验证 采用R 软件“Survival”包对25 个m6A 调控因子进行单因素Cox回归分析，“forestplot”包绘制森林图。选取5个预后相关调控因子（P<0.05），使用LASSO Cox 回归算法建立宫颈癌预后特征和风险评分，风险评分=Codfi×xi，N 为m6A 调控因子数量，Codfi 为系数，xi 为各基因Z-score 转换的相对表达。该公式可计算每例宫颈癌患者风险评分，并将其分为高风险组和低风险组，绘制相应Kaplan-Meier 曲线。绘制ROC 曲线和计算AUC 检验预后模型的预测效能。利用单变量和多变量Cox回归分析评价预后模型和临床特征对宫颈癌预后的价值。

1.4 m6A 调控因子在组织中的表达验证 人类蛋白图谱数据库（Human Protein Atlas，HPA，https：//www.proteinatlas.org/）分析预后相关m6A 调控因子（HNRNPC、LRPPRC、ELAVL1 和FMR1）在宫颈癌组织及配对的正常宫颈组织的表达。收集郑州大学第二附属医院就诊的15 例宫颈癌患者肿瘤组织和正常宫颈组织，qRT-PCR 检测METTL3表达，本研究经郑州大学第二附属医院伦理委员会批准（2021040），采用Trizol 法提取RNA，测定RNA 浓度和纯度，逆转录为cDNA，得到的cDNA 以ACTIN基因为内参进行PCR 反应。PCR 反应体系为SYBR®Premix 10 μl、PCR正向引物 （10 μmol/L） 0.5 μl、PCR反向引物（10 μmol/L）0.5 μl、DNA 模板1 μl、ddH2O 8 μl；PCR 反应条件为95 ℃ 1 min；95 ℃ 15 s，57 ℃30 s，共40个循环。引物序列为：METTL3 F：5'-CCAGCACAGCTTCAGCAGTTCC-3'，R：5'-GCGTGGAGATGGCAAGACAGATG-3'。

1.5 m6A 调控因子聚类分型 采用R 软件“ConsensusClusterPlus”包对25 个m6A 调控因子进行样本聚类分析。参数设定如下：iterations=50，resample rate=80%，分析方法为Pearson 相关性分析。选择k=2 对样本进行分型，将样本分为m6Acluster A组和m6Acluster B 组。主成分分析算法（principal component analysis，PCA）比较两种聚类分型m6A 调控因子表达谱。

1.6 GSEA 分析 对聚类分型的m6Acluster A 和m6Acluster B 进行GSEA 分析。从MSigDB 数据库（http：//www.broadinstitute.org/gsea/msigdb/index.jsp）下载“c2.cp.kegg.v7.4.symbols.gmt”基因集。R 软件“GSVA”包进行GSEA 分析，调整后的P<0.05 为差异有统计学意义。采用“clusterProfiler”包对m6A调控因子进行功能注释，明确模型中潜在的生物学功能和参与的信号通路。

1.7 免疫浸润分析 从Pornpimol Charoentong 获得TIME 浸润免疫细胞基因集［7］。该基因集包括多种人类免疫细胞亚型，如活化的CD8+T细胞、活化的树突状细胞、巨噬细胞、自然杀伤T 细胞、调节性T细胞等。ssGSEA计算聚类分型的两组样本TMIE中免疫细胞浸润相对丰度。

2 结果

2.1 宫颈癌中m6A 调控因子表达和相关性 剔除重复及不完整数据后，共收集宫颈癌样本356例，正常样本20 例。在376 例样本中对25 个m6A 调控因子进行分析，结果显示，12 个m6A 调控因子表达在宫颈癌样本和正常样本中存在差异，其中RBM15、IGF2BP2、IGF2BP3、FMR1、YTHDF1、METTL3 和YTHDF2 表达上调，ZC3H13、METTL14、YTHDC1、FTO、ALKBH3 表达下调（P<0.05，图1A、B）。计算25个m6A调控因子表达的相关性，发现m6A调控因子表达均呈正相关，其中METTL14和YTHDC1正相关性最为显著（r=0.76，图1C、D）。

2.2 m6A 调控因子对宫颈癌预后的评估价值 为研究m6A 调控因子对宫颈癌预后的影响，对25 个m6A 调控因子进行单变量Cox 回归，结果显示FMR1（P=0.017，HR：0.680，95%CI：0.496～0.933）、HNRNPC（P=0.016，HR：2.065，95%CI：1.145～3.725）、LRPPRC（P=0.011，HR：1.711，95%CI：1.133～2.584）、ELAVL1（P=0.028，HR：0.569，95%CI：0.343～0.942）和METTL3（P=0.038，HR：1.551，95%CI：1.026～2.346）是独立的预后因素（图2A）。对5 个m6A 调控因子进行Kaplan-Meier 生存分析，结果显示HNRNPC、LRPPRC 和METTL3 均为危险调控因子，其高表达与较差的总生存期（overall survival，OS）显著相关；而ELAVL1 和FMR1 为保护调节因子，其高表达与较好的OS显著相关（图2B～F）。

2.3 宫颈癌风险预后模型构建及组织样本验证

选 择HNRNPC、LRPPRC、METTL3、ELAVL1 和FMR1 建立风险预后模型，通过LASSO Cox 回归模型计算宫颈癌患者风险评分（图3A、B）。HNRNPC、LRPPRC、METTL3、ELAVL1 和FMR1 的系数分别为0.686 4、0.383 6、0.428 5、-0.880 4 和-0.434 4。根据中位风险评分将患者分为高风险组和低风险组。高风险组OS 与低风险组OS 相比明显降低（P=2.554e-08，图3C）。ROC 曲线评估预后特征的特异度和敏感度，AUC=0.734，表明预测性能良好（图3D）。为验证筛选的5个m6A 调控因子，从HPA数据库分析HNRNPC、LRPPRC、ELAVL1 和FMR1在正常组织和宫颈癌组织中的免疫组化表达，发现LRPPRC和HNRNPC在宫颈癌组织中的表达高于正常组织，ELAVL1和FMR1在正常组织中的表达高于宫颈癌组织（图4A～H）。qRT-PCR 检测收集的15对宫颈癌和正常组织中METTL3 表达，发现宫颈癌组织中METTL3表达更高（P<0.000 1，图4I）。

图4 组织中5个m6A调控因子表达验证Fig.4 Validation of five m6A regulators expressions in tissues

2.4 风险预后模型在宫颈癌中的临床应用价值

绘制5 个m6A 调控因子在高、低风险组中的表达及临床特征热图，结果发现其表达在国际妇产科联盟（Federation International of Gynecology and Obstetrics，FIGO）分期（P<0.05）和生存状态（P<0.001）中差异有统计学意义（图5A）。对风险预后模型分别进行单变量和多变量Cox 回归模型分析，发现 FIGO 分期、风险预后评分与OS 相关（P<0.001），可独立预测宫颈癌患者预后（图5B、C）。

图5 风险预后模型在宫颈癌中的临床应用价值Fig.5 Clinical value of prognostic risk signature in cervical cancer

2.5 宫颈癌中存在的两种m6A 修饰模式 一致性聚类分析显示，Delta 面积在k>3 时增大不明显（图6A）。分型矩阵显示k=2 时无样品特别小的分组，分组内相关性高，分组间相关性低，聚类稳定性最佳（图6B）。因此，基于m6A调控因子的表达确定了两种m6A 修饰模式，命名为m6Acluster A（n=224）和m6Acluster B（n=152，图6C）。生存分析表明，两种m6A 修饰模式在OS 上差异无统计学意义（P=0.299，图6D）。热图显示25 个m6A 调控因子在m6Acluster A中的表达高于m6Acluster B（图6E）。

图6 25个m6A调节剂介导的m6A甲基化修饰模式Fig.6 m6A methylation modification patterns mediated by 25 m6A regulators

2.6 两种m6A 修饰模式生物学功能差异和免疫浸润分析 GSVA 富集分析提示两种m6A 修饰模式存在功能通路差异。m6Acluster A 与RNA 降解、细胞周期、核苷酸修复（切除修复和错配修复）、剪接体等KEGG 通路相关，而m6Acluster B 则与PPAR 通路、细胞色素P450代谢、花生四烯酸代谢、钙信号通路等KEGG 通路相关（图7A）。为研究两种m6A 修饰模式对宫颈癌TIME 的影响，评估了m6A 调控因子表达上调的m6Acluster A 和表达下调的m6Acluster B 免疫细胞的浸润水平，结果发现，m6Acluster B在免疫细胞浸润中显著富集，包括激活的B细胞、激活的CD8+T 细胞、嗜酸性粒细胞、未成熟的B 细胞等，两种修饰模式显示在TIME 中存在不同免疫细胞浸润特征（图7B）。

图7 两种m6A修饰模式生物学功能差异和免疫浸润分析Fig.7 Biological functional differences and immune infiltration of two m6A modification patterns

3 讨论

宫颈癌已成为女性第四大常见癌症死亡原因，严重威胁女性健康［8］。既往研究证实人乳头瘤病毒（human papillomavirus，HPV）感染与宫颈癌密切相关［9］。HPV 筛查和疫苗接种是预防宫颈癌的有效策略［10］。手术及放化疗是早期无转移宫颈癌的主要治疗手段。但传统治疗方法在晚期或复发性宫颈癌的治疗中未取得令人满意的效果，存在预后差、生活质量低等问题［11］。m6A 甲基化是最常见的RNA 修饰，在肿瘤增殖、侵袭和转移中起重要作用，但作用机制尚不明确［12］。因此本研究通过TCGA 和GEO 公共数据库下载宫颈癌转录组表达数据和相应临床数据，通过生物信息学方法和组织验证实验探究m6A 甲基化修饰与宫颈癌发生发展及TIME 免疫细胞浸润的关系，为宫颈癌患者临床诊疗和预后改善提供依据。

本研究首先分析了25 个m6A 甲基化调控因子在宫颈癌组织和正常宫颈组织中的表达和相关性，结果发现在宫颈癌组织中有7个显著上调的m6A调控因子（RBM15、IGF2BP2、IGF2BP3、FMR1、YTHDF1、METTL3和YTHDF2）和5个显著下调的m6A调控因子（ZC3H13、METTL14、YTHDC1、FTO、ALKBH3）。研究发现METTL14 可作为癌基因或抑癌基因在肿瘤中发挥不同作用，为m6A 修饰甲基转移酶相关肿瘤药物研发提供了思路，但目前未见其在宫颈癌中作用的研究［13］。IGF2BP 包括IGF2BP1、IGF2BP2 和IGF2BP3，是m6A阅读蛋白的一种。研究发现IGF2BP在正常和应激条件下均可增强mRNA 稳定性和翻译，导致Myc等致癌产物积累，在宫颈癌和肝癌细胞中敲除IGF2BP 可抑制肿瘤细胞增殖和侵袭［14］。FTO 是第一个被发现的m6A 去甲基化酶，通过调控m6A 修饰E2F1 和Myc 转录本，在宫颈癌细胞增殖和迁移中发挥重要致癌作用［15］。YTH 家族蛋白包括YTHDF1、YTHDF2、YTHDC1与YTHDC2。研究发现YTHDF1和YTHDF2在宫颈癌中表达增加，可促进宫颈癌细胞增殖和糖酵解［16-17］。

近年基于m6A 调控因子的预后模型被用于评估和预测各种肿瘤患者预后，指导后续个体化治疗方 案 制 定［18-20］。本 研 究 选 择5 个m6A 调 控 因 子HNRNPC、LRPPRC、METTL3、ELAVL1 和FMR1，通过LASSO Cox 回归建立了预后风险模型，利用ROC分析预测模型的特异度和敏感度，结果显示AUC=0.734，表明预测性能良好。Kaplan-Meier 生存曲线显示高风险组OS 明显低于低风险组。单变量及多变量Cox回归分析证明该模型可作为独立的预后因子。利用HPA 数据库和qRT-PCR 检测5 个m6A 调控因子在临床组织中的表达，与生物信息学预测结果一致，表明该模型可作为预测宫颈癌患者生存结局的工具。

本研究中，HNRNPC、LRPPRC 和METTL3 均为m6A 的危险调控因子，能增加宫颈癌患者死亡风险。多项研究已证实METTL3 在宫颈癌中高表达，可促进PDK4和HK2基因m6A 修饰，进而增加其mRNA 稳定性，促进宫颈癌发生发展［16，21］。LRPPRC是一种RNA 结合蛋白，在肺腺癌、乳腺癌、胃癌等肿瘤中表达上调，抑制肿瘤细胞凋亡，促进其侵袭和转移［22］。HNRNPC 能够促进乳腺癌细胞和前列腺癌细胞增殖和肿瘤生长［23-24］。ELAVL1 和FMR1 为m6A 调控因子的保护因子改善宫颈癌患者预后。ELAVL1 又称HuR基因，与卵巢癌、尿路上皮癌、食管癌、结直肠癌等多种肿瘤侵袭和不良预后相关［25］。FMR1 通常与YTHDF2 共同参与RNA 剪接、核输出和翻译等，同时被认为与脆性X综合征、卵巢早衰和自然流产等有关［26-27］。

为进一步研究宫颈癌中m6A 甲基化修饰情况，对25个m6A甲基化调节因子进行一致性聚类，将其分为m6Acluster A 和m6Acluster B。生存分析发现两组OS 差异无统计学意义，但25 个m6A 调控因子在m6Acluster A 中的表达明显高于m6Acluster B。GSEA 富集分析探讨m6Acluster A 和m6Acluster B两组间生物学功能和信号通路差异，m6Acluster A与RNA 降解、细胞周期、核苷酸修复（切除修复和错配修复）、剪接体等KEGG通路相关，而m6Acluster B则与PPAR 通路、细胞色素P450 代谢、花生四烯酸代谢、钙信号通路等KEGG 通路相关。ssGSEA 分析两组间肿瘤微环境中免疫细胞浸润水平，发现m6Acluster B 在免疫细胞浸润中显著富集，明显强于m6Acluster A，提示m6Acluster B 存在免疫激活和免疫细胞浸润，可能对应免疫炎症型的“热肿瘤”，对免疫治疗更敏感［28］。

综上，本研究系统证明了m6A 调控因子在宫颈癌中的表达改变和预后价值，构建了宫颈癌风险预后模型并验证其应用价值。同时，基于25个m6A调控因子确定了两种不同的m6A 修饰模式，并与TIME免疫细胞浸润特征相关，为指导免疫治疗策略提供了依据。