黃甜甜,段碧晗,葉迎春,王碩,高凌
甲狀腺癌是一種常見的內分泌系統惡性腫瘤,也是人類最常見的腫瘤之一。在世界許多地區,特別是發達國家,甲狀腺癌的發病率呈上升趨勢。分化型甲狀腺癌在甲狀腺癌中占90%以上,包括甲狀腺乳頭狀癌和濾泡型甲狀腺癌。然而,甲狀腺癌通常是無癥狀的,芬蘭一項經典研究發現,隱匿性甲狀腺癌的患病率高達30%[1]。因此,了解甲狀腺癌發病的分子機制對早期診斷、準確治療及改善預后都至關重要。
N6-甲基腺苷(N6-methyladenosine,m6A)是一種甲基化修飾,可發生在RNA腺嘌呤上,是大多數真核生物最普遍的mRNA修飾形式[2]。在人類及其他許多物種中,m6A在調節基因表達、蛋白質翻譯、細胞行為和生理狀態等方面發揮著重要作用[3]。m6A甲基化是一個動態可逆的過程,根據13個m6A甲基化調控因子的功能不同,可以將其分為三類[4-5]:“書寫者”介導RNA的甲基化修飾過程,包括METTL3、METTL14、KIAA1429、WTAP、RBM15和ZC3H13;“擦除者”介導RNA的去甲基化過程,包括FTO和ALKBH5;“閱讀者”負責讀取RNA甲基化信息,參與下游RNA的翻譯和降解,包括YTHDC1、YTHDC2、YTHDF1、YTHDF2和HNRNPC。“擦除者”之一的FTO與人類肥胖和心理發展有關,而ALKBH5已被證明以去甲基化依賴的方式影響小鼠的精子發生[6-7],這表明在各種生理過程中m6A扮演著重要的角色。
越來越多的研究表明,m6A在許多復雜的人類疾病發生發展中發揮著重要作用[8-9]。關于m6A與胃癌、肝癌、乳腺癌、肺鱗狀細胞癌等多種腫瘤疾病的關系已有較多報道[10-13],而m6A在甲狀腺癌發生和進展過程中的作用仍不十分清楚。本研究利用TCGA數據庫中的RNA測序數據,系統分析了13個廣泛報道的m6A RNA甲基化調控因子在甲狀腺癌中的表達。此外,使用LASSO Cox回歸分析,篩選出3個m6A RNA甲基化調控因子構建LASSO風險回歸模型,預測甲狀腺癌患者的生存和預后,報道如下。
1.1 臨床資料 從TCGA (https://cancergenome.nih.gov/)數據庫中獲取了甲狀腺癌組織510例和正常甲狀腺組織58例的RNA-seq轉錄組數據及相應甲狀腺癌患者臨床病理數據。在分析甲狀腺癌亞組的總體生存率和臨床病理因素時,排除了缺少年齡、性別、病理分級和隨訪數據的樣本,最終保留了501例甲狀腺癌樣本,其中甲狀腺癌患者中男314例,女187例;年齡≤60歲171例,>60歲330例;生存狀態:死亡182例,存活319例;臨床分級:1級36例,2級184例,3級281例;臨床分期:Ⅰ期79例,Ⅱ期152例,Ⅲ期205例,Ⅳ期65例;原發腫瘤:T139例,T2111例,T3215例,T4136例;遠處轉移:M0451例,M150例;淋巴結轉移:N0153例,N1138例,N2104例,N3106例。
1.2 生物信息學分析方法 主要使用Active-Perl 5.24.3 Build 2404(https://www.perl.org)、Rv3.6.0 (https://www.r-project.org/)軟件對數據進行轉換、統計分析和計算。使用R軟件包中的“pheatmap” 軟件包比較甲狀腺癌組織和正常甲狀腺組織中13個m6A RNA甲基化調控因子的表達。使用“Consensus Cluster Plus” 軟件包通過一致性聚類分析方法將甲狀腺癌組織分為2個亞組,探究m6A RNA甲基化調控因子在甲狀腺癌中的作用。主成分分析“PCA”包檢測基于m6A相關基因表達的腫瘤分組的可行性。為了探討m6A RNA甲基化調控因子的預后價值,使用Lasso Cox回歸分析方法來確定與總生存率顯著相關的3個基因,并得到最終風險評分公式。根據每位患者的風險評分,將其分為高危組和低危組,并使用R包中的“survival”“survminer”“pheatmap”軟件包比較2組總體生存率和臨床病理特征。使用2組間的ROC曲線來檢驗生存模型的預測效率,并使用單因素和多因素Cox回歸分析方法確定風險評分的預后價值。
1.3 統計學方法 主要采用R軟件對數據進行統計分析。以Wilcoxon秩和檢驗比較m6A RNA甲基化調控因子在腫瘤組織和正常組織中的表達情況。根據m6A RNA甲基化調控因子的表達情況將患者分為2個亞組后,采用卡方檢驗比較聚類 1和聚類2的性別、年齡、WHO分期和病理分級。在總體生存分析中,采用Kaplan-Meier法檢測甲狀腺癌亞組間、高危組和低危組間的總體生存率。使用ROC曲線分析評估風險評分對5年生存率的預測效率。通過單因素和多因素Cox回歸分析確定風險評分的預后價值。P<0.05為差異有統計學意義。
2.1 m6A RNA甲基化調控因子在甲狀腺癌中的表達情況 分析TCGA數據庫中501例甲狀腺癌樣本和58例正常甲狀腺組織樣本中13個m6A RNA甲基化調控因子的表達情況。除YTHDF2外,12個調控因子均表現出顯著的差異表達(圖1A、B)。與正常組織比較,甲狀腺癌組織的METTL3、YTHDC2、WTAP、ALKBH5、YTHDF1、METTL14、YTHDC1、FTO、KIAA1429、ZC3H13、RBM15顯著下調(P<0.001),HNRNPC顯著上調(P<0.001)。為進一步了解13個m6A RNA甲基化調控因子之間的關系,構建了PPI網絡圖來結果顯示,WTAP、HNRNPC和YTHDC1似乎處于中心位置,并與其他m6A RNA甲基化調控因子相互作用 (圖1C)。此外,大多數m6A RNA甲基化調控因子都顯示出正相關關系。其中,ZC3H13與KIAA1429、KIAA1429與METTL14、METTL14與YTHDC1等3對調控因子呈較強的正相關關系(ρ=0.78),只有HNRNPC和ALKBH5呈負相關關系(ρ=-0.18)(圖1D)。

注:A.甲狀腺癌與正常甲狀腺組織中13個m6A RNA甲基化調控因子的基因表達水平(藍色代表正常組織,粉紅色代表甲狀腺癌組織)。紅色越深表示基因表達上調越明顯,綠色越深代表基因表達下調越明顯;B.小提琴圖可視化m6A RNA甲基化調控因子在甲狀腺癌中的表達差異;C.PPI網絡圖評估m6A RNA甲基化調控因子之間的相互作用;D.Spearman相關分析每個m6A RNA甲基化調控因子之間的關系
2.2 m6A RNA甲基化調控因子的一致性聚類分析 采用一致性聚類分析對501例甲狀腺癌組織進行分組,分組標準如下:(1)累積分布函數增長率較慢;(2)聚類后各組樣本量足夠大;(3)組內相關性高而組間相關性低。在TCGA數據集中,根據m6A RNA甲基化調控因子的表達相似性,最終將甲狀腺癌組織分為2個亞組(圖2A~C)。對甲狀腺癌2個亞組的主成分分析驗證表明,總RNA在亞組中的表達具有特異性,即聚類1似乎有向中心聚集的趨勢,而聚類2則有向四周擴散的趨勢(圖2D),這說明根據m6A RNA甲基化調控因子的一致性聚類分析可以將甲狀腺癌很好地區分為2個亞組。但2個亞組之間存在很多重疊的區域,這表明甲狀腺癌2個亞組之間具有一些共同點。

注:A.k = 2~9的一致性聚類累積分布函數;B.k =2~9時累積分布函數曲線下面積的相對變化;C.當k = 2時,一致性聚類分析將TCGA甲狀腺癌數據分成2個不同的亞組;D.TCGA數據集總RNA表達譜的主成分分析,紅色代表聚類1,藍色代表聚類2
2.3 聚類結果與臨床轉歸及臨床病理特征的關系 進一步分析甲狀腺癌患者的總體生存率及臨床病理特征發現,聚類1的總體生存率明顯高于聚類2(圖3A)。聚類1的4年生存率約為97%,聚類2約為87%。此外,在聚類1中HNRNPC明顯高表達,而ALKBH5的表達明顯降低,進一步說明這兩種調控因子之間的負相關關系。與聚類 1比較,聚類 2與男性、高臨床分期和高T狀態顯著相關(圖3B)。總而言之,聚類結果與甲狀腺癌患者的總體生存率及臨床腫瘤惡性程度密切相關。

注:A.聚類1總體生存率明顯高于聚類2(P<0.05);B.聚類 1和聚類 2臨床病理特征比較差異有統計學意義(P<0.05);與聚類1相比,aP<0.05,bP<0.01


注:A.TCGA數據集中與總體生存率相關的13個m6A RNA甲基化調控因子的單變量Cox回歸分析,根據風險比、95%置信區間確定了3個基因(P<0.05,HR>1);B、C.在TCGA數據集中選擇靶基因的過程;D.根據風險評分將甲狀腺癌患者分為高危組和低危組的總體生存曲線
2.5 甲狀腺癌危險評分與預后及臨床病理特征的關系 ROC曲線結果表明,風險評分對甲狀腺癌患者的5年生存率有較高的預測準確率(AUC=0.717)(圖5A)。為進一步了解高危組與低危組甲狀腺癌患者3個m6A RNA甲基化調控因子的表達情況及其臨床病理特征的關系,再對不同分組與臨床病理特征(性別、年齡、生存狀態、T狀態)進行相關性分析并繪制其m6A RNA甲基化調控因子的表達熱圖(圖5B)。結果顯示,與低危組比較,高危組RBM15、KIAA1429、FTO表達明顯增高,且高危組和低危組甲狀腺癌患者在T狀態之間存在顯著差異。

注:A.ROC曲線顯示風險評分的5年生存率預測價值;B.熱圖顯示3種m6A RNA甲基化調控因子在甲狀腺癌患者高危組和低危組的表達水平以及臨床病理特征差異;與甲狀腺癌患者高危組相比,aP<0.05,bP<0.01
2.6 甲狀腺癌的獨立預后指標分析 單因素Cox回歸分析顯示,年齡、風險評分、分期、T狀態均與總體生存率相關。多因素Cox回歸分析結果顯示,年齡和風險評分可以作為甲狀腺癌的獨立預后指標,即甲狀腺癌的風險隨著年齡和風險評分的增加而增加,見圖6。

注:A.TCGA數據集中臨床病理因素和風險評分與患者總體生存之間的單變量Cox回歸分析;B.TCGA數據集中臨床病理因素和風險評分與患者總體生存之間的多變量Cox回歸分析
甲狀腺癌是最常見的內分泌系統惡性腫瘤,近40年來在世界范圍內發病率呈上升趨勢。超聲檢查的進步可能是造成這種趨勢的原因[14],然而,導致這一增長的其他原因尚不清楚。雖然甲狀腺癌相對于其他惡性腫瘤具有較好的惰性和良好的長期生存率,但甲狀腺癌可發生在任何年齡,通常對年輕人的影響更大[15]。m6A RNA甲基化調控因子在腫瘤的發生發展中起重要作用,且已被證實與多種癌癥密切相關,靶向m6A RNA甲基化的抗癌藥物在臨床的應用前景非常廣闊[16]。因此,有必要探討m6A RNA甲基化調控因子與甲狀腺癌的關系或對甲狀腺癌預后的預測價值。
在本研究中,首先分析了m6A RNA甲基化調控因子在甲狀腺癌和正常甲狀腺組織中的表達情況及其與不同臨床病理的關系。甲狀腺癌組織和正常組織中多達12個m6A RNA甲基化調控因子的表達比較差異有統計學意義。
急性髓系白血病患者中METTL3表達升高,并發揮致癌作用。METTL3通過增加cMYC、BCL-2和PTEN的翻譯,抑制細胞分化和凋亡,促進細胞增殖。MELLT3還激活PI3K / AKT通路,控制細胞分化和自我更新[17]。肝細胞癌中METTL3表達升高,METTL3在體外可促進肝癌細胞的生長、遷移和集落形成,在體內可影響肝癌的成瘤性、生長和肺轉移,且METTL3水平升高提示預后不良[16]。METTL3在卵巢癌中也高表達,且METTL3高表達與腫瘤分級、FIGO分期及總生存率顯著相關[18]。近期,有學者進行了METTL3促進甲狀腺癌進展的相關研究[19],研究表明,METTL3誘導TCF1 mRNA m6A甲基化,上調TCF1蛋白水平,TCF1通過激活Wnt通路進一步刺激腫瘤細胞轉移,從而加重甲狀腺癌的惡性進展,提示METTL3是一種潛在的癌基因。然而,本研究顯示,METTL3在甲狀腺癌患者中的表達遠低于正常組織,提示METTL3可能阻礙了甲狀腺癌的進展。這表明METTL3在甲狀腺癌發生發展中的作用可能還需要更多的研究來證明。
在探索m6A RNA甲基化調控因子是否對甲狀腺癌有預后價值的過程中,進行了Cox單變量分析和LASSO回歸分析,鑒定篩選出RBM15、KIAA1429和FTO等3個基因。研究發現,RBM15、KIAA1429和甲狀腺疾病的關系研究甚少,而FTO與甲狀腺癌密切相關。
RBM15對基因敲除模型小鼠中多種組織的發育至關重要,尤其對維持造血干細胞長期穩態和巨核細胞、B細胞的分化具有重要作用[20]。此外,RBM15參與了染色體易位t (1;22),產生與急性巨核細胞白血病相關的RBM15-MKL1融合蛋白[21]。除參與急性巨核細胞白血病的發生外,RBM15還是X染色體沉默的重要因素,并已證實在乳腺組織中表達[22]。研究表明,雌激素可獨立增加RBM15的表達,并與乳腺疾病的發生有關。對于甲狀腺癌,女性的發病率高于男性,這是否與雌激素信號傳導有關還有待進一步研究。
KIAA1429已被證實通過不同途徑影響多種癌癥,包括通過依賴于GATA3的m6A轉錄后修飾促進肝細胞癌的進展,以及通過直接靶向c-jun mRNA調節胃癌細胞增殖[23]。此外,一項關于乳腺癌的研究表明,KIAA1429可通過調節CDK1表達而成為乳腺癌的致癌因子[24]。各種CKD抑制劑被廣泛應用于臨床治療惡性腫瘤,包括甲狀腺癌[25]。因此,KIAA1429與甲狀腺癌的關聯機制值得進一步探討。
FTO在乳腺癌、急性髓系白血病和膠質母細胞瘤中高表達,并促進這些癌癥的進展,現階段已開發出相關藥物[10,26-27]。大黃酸是第一種有效的FTO抑制劑,通過競爭性結合催化區域和單鏈RNA底物抑制FTO。大黃酸還能有效抑制m6A去甲基化,提高體外m6A水平[28]。然而,本研究發現FTO在甲狀腺癌中表達水平較低。此外,它在膀胱癌中也呈低表達。m6A調控基因的致癌作用似乎存在爭議。這可能是由于癌癥的異質性,相同的基因在不同的癌癥中扮演不同的角色,導致m6A調節基因在不同癌癥中的表達差異。FTO更多的是作為一個與肥胖、糖尿病和代謝綜合征相關的基因進行研究,在一些大型的研究中發現,BMI增加已成為甲狀腺癌一個新的危險因素。
綜上,m6A RNA甲基化調控因子在甲狀腺癌中異常表達并具有預后價值,為進一步探討m6A甲基化在甲狀腺癌中的作用和尋找相關治療靶點提供了重要證據。
利益沖突:所有作者聲明無利益沖突
作者貢獻聲明
黃甜甜:設計研究方案,實施研究過程,論文撰寫;高凌:提出研究思路,分析試驗數據,論文審核;段碧晗:實施研究過程,資料搜集整理,論文修改;葉迎春:進行統計學分析;王碩:課題設計,論文撰寫