基于高通量测序分析头颈肿瘤中CircRNAs的表达谱差异

王榴倩，沈志森，顾姗姗，邓红霞，李赞

头颈肿瘤（HNC）是来源口腔上皮、咽、喉、食道、鼻腔、唾液腺和甲状腺的恶性肿瘤[1]，5年生存率仅为40%～50%[2]，主要危险因素是吸烟、饮酒、咀嚼槟榔及病毒感染等[3]。近年来高通量测序技术以其高灵敏度、高特异性、高完整性和低成本逐渐取代基因芯片成为基因研究的主要技术[4]。CircRNAs是一类内源性非编码RNA（ncRNA）[5]，已有研究发现CircRNAs的表达失调在HNC的发生与进展中发挥关键作用[6-8]，但是CircRNAs在HNC中的功能机制少见报道。因此本研究旨在通过高通量测序探索HNC中的新型CircRNAs，结合实验验证找到有意义的分子标记物，为HNC的临床早期诊断和精准治疗提供科学依据。现报道如下。

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患者临床资料的相关性