m6A调节因子在肺腺癌和肺鳞癌中的生物学差异

莘 华,刘 智,张洁莉

在我国,肺癌位居恶性肿瘤的第一位,且发病率及病死率占全球的三分之一[1],其中以非小细胞肺癌(non-small cell lung cancer, NSCLC)最常见,占肺癌的80%～85%。NSCLC最主要的两种病理组织学亚型是肺腺癌(Lung adenocarcinoma,LUAD)和肺鳞癌(Lung squamous cell carcinoma,LUSC),二者在临床表现上相似,但在发病机制、治疗方案及患者预后等方面都存在着明显的差异[2]。N6-甲基腺苷(N6-methyladenosine,m6A)修饰是在RNA分子上的一种表观遗传修饰,是一种动态可逆的修饰模式,参与调控信使RNA(mRNA)的稳定性、转运、翻译和降解等[3]。m6A甲基化修饰过程主要依靠三种类型酶参与完成,即甲基化转移酶(Writer)、去甲基化酶(Eraser)和甲基化识别酶(Reader)。既往研究表明,m6A甲基化修饰与NSCLC的发生发展具有密切的关系。Li等[4]发现,m6A去甲基化酶FTO可以通过调节USP7 mRNA的m6A修饰水平进而促进肺癌细胞的生长。Jin等[5]研究表明,m6A去甲基化酶ALKBH5通过减少YTHDFs介导的YAP表达和抑制miR-107/LATS2介导的YAP在NSCLC中的活性来抑制肿瘤的生长和转移。Du等[6]对METTL3的SUMO化的研究发现,其可与METTL14和WTAP发生相互作用,进而抑制甲基化转移酶活性,从而使mRNA中m6A甲基化的修饰水平下降。同时发现,在NSCLC细胞中将METTL3敲除时,可以直接影响肿瘤细胞的增殖与转移。METTL3在NSCLC中还可以与eIF3h发生相互作用,进而增强翻译,促使癌基因转化形成密集的多核糖体,该发现可以作为NSCLC潜在的治疗靶点[7]。除此之外,一些研究还表明m6A调节因子的表达与NSCLC预后显著相关[8-10]。然而,对于m6A调节因子在非小细胞肺癌不同病理亚型中发挥的生物学功能的差异性的研究尚未有报道。本研究利用生物信息学的方法,基于LUAD和LUSC的表达谱数据,通过比较分析26个m6A调节因子在LUAD和LUSC肿瘤与正常组织中的表达、通路活性及生存预后,明确m6A调节因子在LUAD和LUSC发挥调控功能时的差异性,为实现针对NSCLC具体病理组织亚型的精准治疗和寻找潜在靶点提供理论基础。

1 材料与方法

1.1 材料 利用癌症基因组图谱(the cancer genome atlas,TCGA)数据库,检索含有人源的LUAD和LUSC的mRNA表达谱数据集,分别选取含有LUAD和LUSC肿瘤组织和正常组织对照研究的表达谱数据作为研究对象。纳入研究的LUAD表达谱数据共包含576例样本,其中肿瘤组织样本518例,正常组织样本58例,58例有配对样本,515例具有临床信息;LUSC表达谱数据共包含554例样本,其中肿瘤组织样本503例,正常组织样本51例,51例有配对样本,501例具有临床信息。

1.2 m6A调节因子的差异表达分析 下载的LUAD和LUSC表达谱数据经RNA测序期望最大化(RNA-seq by expectation maximization,RSEM)方法分析,同时进行归一化和批次差异校正处理后,获得研究对象的表达矩阵,筛选出26个m6A调节因子在LUAD和LUSC的表达谱数据,随后利用ggplot2包对LUAD和LUSC的肿瘤组和正常组织进行差异表达分析26个m6A调节因子的基因信息(表1)。

表1 26个m6A调节因子基因信息

1.3 m6A调节因子参与的通路活性分析 10个肿瘤相关经典通路包括结节性硬化症蛋白复合物/哺乳动物雷帕霉素靶点(tuberous sclerosis complex/mammalian target of rapamycin,TSC/mTOR)、受体酪氨酸激酶(receptor tyrosine kinase, RTK)、Ras蛋白/丝裂原活化蛋白激酶(ras/mitogen-activated protein kinase,RAS/MAPK)、磷脂酰肌醇-3-激酶/丝氨酸/苏氨酸激酶(PI3K/AKT)、激素雌激素受体(estrogen receptor,ER)、雄激素受体(androgen receptor,AR)、上皮细胞-间充质转化(epithelial-mesenchymal transition,EMT)、DNA损失应答(DNA damage response,DDR)、细胞周期(cell cycle)和凋亡(apoptosis)被用来进行通路活性分析。来自于癌症蛋白质组图谱(the cancer proteome atlas,TCPA)的反相蛋白微阵列(reverse phase protein array,RPPA)数据被用来计算10个肿瘤相关通路的通路活性评分(pathway activity score,PAS)[11]。以基因表达量的中位数作为阈值将样本分为高和低两组,当PAS(目标基因高表达)大于PAS(目标基因低表达)时,认为目标基因可能对某个通路有激活作用,否则对某个通路有抑制作用[12]。

1.4 m6A调节因子的预后生存分析 根据病例编号(Case_submitter_id)分别合并从TCGA数据库中下载的LUAD和LUSC的mRNA表达谱数据和临床生存数据,筛选出26个m6A调节因子并将其表达值的中位值作为阈值,将肿瘤样本分为高表达组和低表达组。使用R软件中的forestplot、survival和survminer包,应用单因素Cox分析(Univariate Cox analysis)和Kaplan-Meier 生存分析方法,逐个基因分析其对患者总生存期(Overall survival,OS)的影响,计算风险比( Hazard ratio,HR),分析m6A调节因子的表达对LUAD和LUSC患者OS的影响差异。

1.5 统计学处理 应用R软件(v 4.0.5)进行差异表达分析和通路活性分析,并采用配对样本秩和检验(wilcoxon test)估算P-value值,预后生存分析采用时序检验(logrank test)进行统计学分析,P<0.05为差异有统计学意义。

2 结 果

2.1 m6A调节因子在LUAD和LUSC的表达差异 通过配对样本wilcoxon检验比较分析26个m6A调节因子在58对LUAD和51对LUSC的肿瘤和正常组织中的表达差异情况,结果显示,METT5、VIRMA和RBM15、YTHDF1、EIF3A、IGF2BP1、IGF2BP3、HNRNPA2B1、HNRNPC和RBMX在LUAD和LUSC的肿瘤组织中相比于正常组织均显著上调表达,METTL14、METTL16和ZC3H13和FTO在LUAD和LUSC的肿瘤组织中相比于正常组织均显著下调表达。除此之外,ZCCHC4、ALKBH1和YTHDF2仅在LUAD的肿瘤组织和正常组织中的表达水平差异有统计学意义,在LUSC的肿瘤组织和正常组织中,表达水平差异无统计学意义。类似地,YTHDC2和IGF2BP2仅在LUSC的肿瘤组织和正常组织中表达水平存在统计学的意义(图1)。

图1 m6A调节因子在LUAD和LUSC肿瘤和正常组织的表达分析

2.2 m6A调节因子在LUAD和LUSC的通路活性差异 26个m6A调节因子在10个经典的肿瘤相关通路活性分析发现,有18个和16个m6A调节因子对应地在LUAD和LUSC中参与了至少一个肿瘤相关通路的激活或抑制。其中,METTL3、METTL5、METTL14、ZC3H13、RBMX和FTO仅在LUAD中的Apoptosis、Cell cycle、EMT和RAS/MAPK等肿瘤相关通路中发挥激活或抑制功能;EIF3A、YTHDC2和RBM15B仅在LUSC中的PI3K/AKT和TSC/mTOR信号通路中发挥抑制或激活作用。除此之外,相同的基因在LUAD和LUSC中参与调控不同的通路活性,例如:YTHDF1在LUAD中参与Cell cycle的激活,在LUSC中则参与TSC/mTOR通路的激活;YTHDF2在LUAD和LUSC中均可参与TSC/mTOR通路的激活,在LUAD中该基因同时发挥对Apoptosis通路的抑制作用(图2)。

2.3 m6A调节因子在LUAD和LUSC的生存预后价值 利用TCGA数据库中515例LUAD和501例LUSC患者的临床信息,结合单因素COX回归分析和Kaplan-Meier生存分析方法预测m6A调节因子与患者OS的关系。分析结果显示,HNRNPC、IGF2BP3、HNRNPA2B1、YTHDF2、METTL5、IGF2BP2、VIRMA和IGFBP1调节因子与LUAD的预后显著相关,仅FTO、IGF2BP1和IGF2BP3与LUSC的预后显著相关(表2)。Kaplan-Meier生存分析结果提示,HNRNPC、IGF2BP3、HNRNPA2B1、METTL5、IGF2BP2、VIRMA和IGFBP1高表达与LUAD患者较长的OS呈负相关,而YTHDF2高表达与LUAD患者较长的OS呈正相关(表3)。相反,在LUSC中,IGFBP1和IGF2BP3高表达与患者具有更长的OS显著正相关,同时,FTO低表达与LUSC患者较长的OS显著正相关(表3)。

3 讨 论

肺鳞癌和肺腺癌是非小细胞肺癌的常见病理组织亚型,其中周围型腺癌居多,鳞癌多为中央型[13]。由于两者之间的基因组、转录组和免疫等特征不同,其对应治疗方案、预后也均存在较大差异[14-17],因此,针对肺癌患者不同病理亚型进行精准治疗具有重要意义。同时,肿瘤与正常组织中基因的异常表达与调控的差异,对肿瘤的发生、发展及生存预后也具有重要的影响。目前,m6A调节因子通过调节mRNA的甲基化水平影响肿瘤相关基因的表达影响肿瘤患者的预后已被广泛报道[18-20]。临床对NSCLC中关于m6A调节因子的调控功能也有很多研究,但大多没有进行细致的病理分型[21-23]。另外,由于各自研究中纳入分析的数据集不同,数据的处理方法不同,所得到的结果存在一定差异。本研究选择了TCGA数据库中含有LUAD和LUSC肿瘤和正常组织的表达谱数据和临床患者信息作为研究对象,从表达、通路活性及预后3个方面比较了26个m6A调节因子在LUAD和LUSC中发挥的生物学作用的差异性。

对m6A调节因子在NSCLC中发挥的生物学功能的研究主要目的是明确m6A调节因子的功能机制及对NSCLC预后的影响,阐释其对相关治疗方式耐药的潜在机制,探索其成为潜在靶点的可能性。本研究对26个m6A调节因子的差异表达分析发现,ZCCHC4和ALKBH1仅在LUAD的肿瘤组织中显著上调表达,而在LUSC的肿瘤组织和正常组织中,表达水平无明显差异。既往研究提示,ZCCHC4在多个人类癌症组织中异常高表达,并与肝细胞癌(HCC)、胰腺癌和结肠癌患者的不良预后和化疗耐药相关[24]。研究[25]发现,ALKBH1的高表达与肿瘤的发生发展密切相关,ALKBH1可以通过调节m6A RNA的去甲基化进而促进肺癌的发展。ZCCHC4、ALKBH1及其他m6A调节因子在LUAD和LUSC中的差异,将可能成为LUAD和LUSC患者区别治疗或相同治疗产生不同应答结果潜在的靶点或潜在的分子机制。

m6A调节因子在LUAD和LUSC中的差异表达模式将对参与调控的肿瘤相关通路活性产生不同的影响。在我们的研究中,YTHDF2在LUAD的肿瘤组织中显著上调表达,而在LUSC的肿瘤组织和正常组织中,表达则无明显差异。同时,YTHDF2在LUAD和LUSC中均可参与TSC/mTOR通路的激活,但在LUAD中该基因同时发挥对凋亡通路的抑制作用,可能是导致LUAD在治疗过程中肿瘤细胞耐药的诱因之一。类似地,IGF2BP2仅在LUSC的肿瘤组织中显著上调表达,在LUAD的肿瘤和正常组织中则无明显差异,而IGF2BP2在LUSC中仅参与激素雌激素受体信号通路的抑制,在LUAD中则同时参与凋亡、细胞周期和EMT信号通路的激活及PI3K/AKT和RAS/MAPK信号通路的抑制,该结果进一步强调了不同病理亚型中m6A调节因子发挥功能的差异性。

对于NSCLC不同的病理亚型,相同的m6A调节因子对患者预后的影响也不相同。Kazuo等[26]发现,YTHDF2高表达与NSCLC更好的生存预后相关,而本研究发现,YTHDF2在LUAD中高表达与较好的预后生存相关,在LUSC中则表现为与预后无显著相关性。除此之外,Han等[27]发现,IGF2BP2高表达与较短的OS和较差的预后相关,本研究发现,IGF2BP2虽然在LUSC的肿瘤组织中显著上调表达,但却与预后没有显著的相关性。相反IGF2BP2在LUAD的肿瘤和正常组织中的表达水平虽然无明显差异,但其IGFBP2高表达却与更差的预后显著相关。本研究发现,FTO基因高表达仅与LUSC的不良预后相关,与LUAD的预后无显著相关性。与Wang等[28]发现,FTO基因通过诱导NELL2表达来抑制E2F1的m6A甲基化修饰从而加速NSCLC的进展有一定的差异。本研究中存在一些m6A调节因子在表达水平和通路活性上无明显差异,且与LUAD和LUSC患者的预后无相关性,如ALKBH5这类m6A调节因子可能通过其他形式,在 LUAD和LUSC的发生、发展中起到重要调控作用,通过影响PD-L1的表达和免疫微环境进而对肿瘤的免疫治疗应答疗效产生不同的治疗结局[29],相关的结果尚待进一步的研究证实。

本研究中存在一些不足之处,主要包括配对样本数量相对较少,且样本数据来源于TCGA数据库,其中包含关于中国人群的数据有限,临床信息不够完善,无法针对不同病理亚型临床特征进行更细致的分析。

综上所述,本研究应用生物信息学的相关分析方法对26个m6A调节因子在LUAD和LUSC的表达谱数据进行挖掘。有助于对m6A调节因子在NSCLC的发生、发展和转移的分子过程进一步认识,实现针对LUAD和LUSC的精确分层,同时为研究LUAD和LUSC潜在的生物标志物及靶点提供参考依据,随着后续数据库的不断完善及我们对m6A调节因子的深入研究,会针对性的开展相关内容的分析。