基于高通量測序分析頭頸腫瘤中CircRNAs的表達譜差異

王榴倩，沈志森，顧姍姍，鄧紅霞，李贊

頭頸腫瘤（HNC）是來源口腔上皮、咽、喉、食道、鼻腔、唾液腺和甲狀腺的惡性腫瘤[1]，5年生存率僅為40%～50%[2]，主要危險因素是吸煙、飲酒、咀嚼檳榔及病毒感染等[3]。近年來高通量測序技術以其高靈敏度、高特異性、高完整性和低成本逐漸取代基因芯片成為基因研究的主要技術[4]。CircRNAs是一類內源性非編碼RNA（ncRNA）[5]，已有研究發現CircRNAs的表達失調在HNC的發生與進展中發揮關鍵作用[6-8]，但是CircRNAs在HNC中的功能機制少見報道。因此本研究旨在通過高通量測序探索HNC中的新型CircRNAs，結合實驗驗證找到有意義的分子標記物，為HNC的臨床早期診斷和精準治療提供科學依據?，F報道如下。

1 資料與方法

1.1 一般資料 收集2012年1月至2019年12月在寧波大學附屬李惠利醫院耳鼻咽喉頭頸外科及湖南省腫瘤醫院接受手術的155例HNC患者的腫瘤組織及癌旁組織，病理均診斷為HNC。其中18例HNC患者的腫瘤組織及癌旁組織標本用于高通量測序分析，其他137例HNC患者的腫瘤組織和癌旁組織樣本用于qRT-PCR驗證部分CircRNAs的差異表達。本研究經寧波大學倫理審查委員會批準。

1.2 方法

1.2.1 總RNA提取和高通量測序分析將所取得的新鮮標本立即浸泡在含甲醛的低溫冷凍管中，并保存在－80℃的冰箱中以備后續提取總RNA。NanoPhotometer?分光光度計用于檢測樣品純度（IMPLEN，CA，USA），Qubit?3.0熒光儀（Life Technologies，CA，USA）用于檢測RNA樣本濃度，Agilent 2100 RNA Nano 6000檢測試劑盒（Agilent Technologies，CA，USA）用于檢測RNA樣本的完整性和濃度。每個樣本以3 g總RNA為起始量構建CircRNAs文庫。使用Ribo-zero?Gold試劑盒從樣品中去除核糖體RNA（rRNA），根據Illumina（NEB，Ispawich，USA）的說明選擇不同的索引標簽建立文庫。測序策略為PE150。通過bcl2fastq2軟件將Base調用轉換為Illumina高通量測序原始圖像數據的序列讀取，即原始數據。使用Perl腳本處理原始數據，以確保后續分析中使用數據的質量。DEG seq或DE seq用于有或沒有重復的兩個樣本的差異表達分析。q≤0.05和|log2_ratio|≥1的基因被鑒定為差異表達基因。

1.2.2 基因本體（GO）和KEGG通路富集分析 對差異表達基因進行GO（http://geneontology.org/）富集分析，根據每個GO Term所包含的差異表達基因個數，應用超幾何檢驗方法，找出與整個基因組背景相比顯著富集的GO Term。計算得到的P值通過多重檢驗校正之后，選取q＜0.05的GO term作為在差異表達基因中顯著富集的GO Term。KEGG（http://www.kegg.jp/）是全基因組及代謝途徑方面的數據庫。對KEGG中每個通路應用超幾何檢驗進行富集分析，對得到的P值通過多重檢驗進行校正，以q＜0.05為閾值，滿足此條件的通路定義為在差異表達基因中顯著富集的通路。

1.2.3 qRT-PCR驗證 用TRIzol試劑（Life Technologies，USA）從腫瘤和癌旁組織樣品中提取總RNA。根據Promega公司的GoScript TM反轉錄系統反轉錄試劑盒的說明，將總RNA逆轉錄為cDNA。在MX3005P熒光定量PCR儀上按照說明書進行qRT-PCR反應，反應體系為：5 l cDNA，5.5 l游離水，12.5 l 2xqPCR主混合液，1 l正向引物，1 l反向引物。對照基因GAPDHmRNA作為管家內參基因。qRT-PCR實驗過程如下：95℃預變性10 min，95℃變性15 s，59℃退火30 s，72℃延伸30 s，共40個循環。qRT-PCR實驗引物如下：hsa_circ_0012212：forward，5’-CCAGCACTGGGCATGGAG-3’and reverse，5’-GGCGTTCTAGGTCCTGCTC-3’；hsa_circ_0016148：forward，5’-TCTTCAAGGGAATCCTCCGC-3’and reverse，5’-CCTTAACCAGCAGACTGGGG-3’；The primer’s sequences of glyceraldehyde-3-phosphate dehydrogenase(GAPDH)，as a control：forward，5’-CATGAGAAGTATGACAACAGCCT-3’and reverse，5’-AGTCCTTCCACGATACCAAAGT-3’。利用Ct法分析靶基因的表達水平，△Ct=Ct（靶基因）－C（t內參基因）。

1.3 統計方法 采用GraphPad Prism 7.0、Microsoft Excel 2010及SPSS 18.0軟件進行分析，計量資料以均數±標準差表示，采用配對t檢驗；基因表達水平與臨床資料的相關性采用單因素方差分析和非參數檢驗；建立受試者工作特征（ROC）曲線評估hsa_circ_0012212和hsa_circ_0016148的診斷價值。P＜0.05為差異有統計學意義。

2 結果

2.1 高通量測序識別CircRNAs差異表達譜 18例HNC患者的腫瘤組織及癌旁組織測序后共檢測到40577個CircRNAs，共714個CircRNAs在腫瘤組織與癌旁組織中的表達有差異（P＜0.05），其中有382個上調，332個下調。對714個表達顯著差異的CircRNAs進行分層聚類分析，結果見封二彩圖7a所示；火山圖譜進一步闡明了CircRNAs在腫瘤組織組和癌旁組織組中的表達差異（封二彩圖7b）。

圖7 差異表達CircRNAs的聚類圖（a）和火山圖（b）

2.2 CircRNAs的GO和KEGG通路分析 差異表達CircRNAs的生物過程變化主要發生在5-甲基胞嘧啶分解代謝過程和5-甲基胞嘧啶代謝過程中。分子功能的變化主要集中在（R）-檸檬烯6-單加氧酶活性、（S）-檸檬烯6-單加氧酶活性、（S）-檸檬烯7-單加氧酶活性和ATP結合。細胞組成分析揭示了與高爾基體有關的CircRNAs的差異表達。KEGG通路分析顯示，差異表達的CircRNAs在ErbB信號通路和腫瘤Micro RNAs信號通路顯著富集。此外，差異表達的CircRNAs也在胰島素維持、破骨細胞分化、神經營養因子信號通路和孕激素介導的卵母細胞成熟通路富集，排序前11的KEGG通路富集圖也顯示了相似的結果。

2.3 qRT-PCR驗證部分差異表達CircRNAs對高通量測序結果中差異表達的CircRNAs進行ceRNA預測分析，篩選出部分能參與構成CircRNA-miRNAmRNA互作關系的CircRNAs，從中選出表達明顯下調的hsa_circ_0012212、hsa_circ_0016148 二 條 CircRNAs。hsa_circ_0012212和hsa_circ_0016148在腫瘤組織中的表達水平均低于癌旁組織（均P＜0.05）（封二彩圖8）。

圖8 qRT-PCR驗證，a為hsa_circ_0012212，b為hsa_circ_0016148

2.4 差異表達CircRNAs與臨床資料的相關性 hsa_circ_0012212的表達與年齡有關（P＜0.05），而與腫瘤類別、吸煙、飲酒、分化程度、TNM分期、侵襲能力及淋巴結轉移均無關（均P＞0.05）。hsa_circ_0016148與腫瘤分類和淋巴結轉移有關（P＜0.05），與其他因素均無關（均P＞0.05）。見表1～2。

表1 hsa_circ_0012212表達水平（Ct）與HNC患者臨床資料的相關性

2.5 ROC曲線分析 hsa_circ_0012212的曲線下面積（AUC）為0.628，hsa_circ_0016148的AUC為0.596，hsa_circ_0012212作為HNC診斷分子標志物的敏感性更高。

3 討論

高通量測序和計算分析顯示環狀轉錄本在真核生物轉錄組中廣泛存在[9-10]，在各種組織和不同物種中大量穩定表達[11-14]，通過反向剪接外顯子或內含子的pre-mRNA新發現了大量的circRNAs[15-16]。CircRNAs的高度保守和穩定性（半衰期＞48 h）提示circRNAs在癌癥診斷方面可能優于miRNAs和長鏈非編碼RNA（lncRNAs）。

Ke等[17]發現circHIPK3在鼻咽癌組織和細胞系中高表達，circHIPK3通過保護ELF3免受miRNA-4288介導的沉默作用從而促進了鼻咽癌的進展，同時circHIPK3的耗竭顯著抑制了體內腫瘤的生長和轉移。Wu等[18]發現CircRNA hg19_circ_0005033可促進 CD133+CD44+喉癌干細胞的增殖、遷移和侵襲，提示CircRNA與癌癥免疫調節相關。Feng等[19]通過CircRNA基因測序發現，下咽癌腫瘤組織與鄰近正常組織中有173個CircRNA差異表達，其中上調的CircRNA 71個，下調的CircRNA 102個。還發現了CircRNA（hsa_circ_0008287和hsa_circ_0005027）和miRNA（has-miRNA-548c-3p）組成的ceRNA網絡，該網絡對ErbB和Hippo信號通路均有顯著影響。以上這些研究闡明了CircRNA在HNC中的作用，啟示研究HNC癌變和發展的分子機制有助于發現早期診斷的生物標志物和有效的治療靶點。

本研究結果顯示，差異表達的CircRNAs有714個，其中上調的CircRNAs有382個，下調的CircRNAs有332個。此外，還發現了一些circBase中沒有包含的CircRNAs，這些新發現的CircRNAs豐富了非編碼RNA家族。然后還發現CircRNAs的主要功能和通路與腫瘤相關，如ErbB信號通路和腫瘤MicroRNAs信號通路，為進一步研究HNC的發病機制提供了重要線索。在本研究中，通過選擇2個差異表達的CircRNAs進行qRT-PCR驗證，顯示在腫瘤組織中hsa_circ_0012212和hsa_circ_0016148表達下調。本研究結果發現在HNC腫瘤組織中低表達的hsa_circ_0012212與年齡相關，hsa_circ_0016148與腫瘤分類和淋巴結轉移相關（均P＜0.05），這提示不同的CircRNA有不同的診斷價值。ROC曲線分析顯示hsa_circ_0012212的AUC為0.628，比hsa_circ_0016148的0.596要高，這說明hsa_circ_0012212在HNC的分子診斷中具有更高的敏感性。本研究HNC中各腫瘤類型的樣本量不相同，在將來需要對各個腫瘤類型作擴大樣本量的研究。盡管hsa_circ_0012212和hsa_circ_0016148在HNC中的作用才剛被揭示，但是其作為診斷分子標志物和治療的潛在價值將為HNC的生物治療帶來更廣的應用前景。

表2 hsa_circ_0016148表達水平（Ct）與HNC患者臨床資料的相關性