口腔白斑病组织中环状RNA 的差异表达谱分析

徐思明， 宋雨翰， 邵彦雄， 陶岚， 周海文

1.上海交通大学医学院附属第九人民医院·口腔医学院口腔黏膜病科，国家口腔疾病临床研究中心，上海市口腔医学重点实验室，上海市口腔医学研究所，上海（200011）； 2. 上海交通大学医学院附属第九人民医院·口腔医学院口腔综合科，国家口腔疾病临床研究中心，上海市口腔医学重点实验室，上海市口腔医学研究所，上海（200011）

口腔白斑（oral leukoplakia，OLK）是一种发生在口腔黏膜上以白色为主、不能擦除的斑片或斑块，临床或病理上不能诊断为其他任何疾病的口腔黏膜潜在恶性疾病。流行病学显示OLK 的恶性转化率高达3%～18%，非均质的OLK 恶性转化率甚至高达15%～40%［1］。OLK 的病因较为复杂，发病机理仍不清楚，目前临床上没有特效治疗方法，因此研究OLK 发病机制并确定新的潜在分子标志物具有重要意义。环状RNA（circRNA）是具有共价闭环结构的非编码RNA［2］。有报道表明circRNA 在肿瘤、神经系统、心血管疾病、炎症及免疫疾病中发挥重要的作用［3⁃7］。研究表明circRNA可作为有效的生物标志物和治疗靶点［8⁃11］。但目前还没有OLK 相关circRNA 的报道。本研究通过高通量测序描绘OLK 和正常口腔黏膜（normal oral mucosal，NOM）组织中circRNA 的表达谱，并通过对筛选出的显著差异表达的circRNA 进行验证分析，同时构建circRNA⁃miRNA⁃mRNA 的竞争性内源RNA（competing endogenous RNA，ceRNA）网络，为circRNA 在OLK 发病机制中的机理研究提供新的研究方向。

1 资料和方法

1.1 研究对象

选取2018 年9 月至11 月于上海第九人民医院口腔黏膜科接受治疗的OLK 患者的白斑组织26例，男女比例1∶1。患者的组织样本都通过组织病理学检查确诊为OLK。26 例NOM 组织来自整复外科及口腔外科拔牙的手术患者。组织纳入标准为不伴有系统性疾病，且3 个月内未接受过免疫刺激治疗（包括类固醇治疗），无放疗、化疗史。所有的组织样本取样后都立即放入液氮保存。提供组织样本的患者均于术前签署知情同意书，本研究已通过上海第九人民医院伦理委员会审批［审批编号：（2012）21］。从收集的52 例样本中选取6 例OLK 组织及6 例NOM 组织进行高通量测序，样本信息见表1，并采用qRT⁃PCR 验证这52 例组织中筛选重点circRNA 的表达情况。

1.2 主要试剂和仪器

TRIzol 裂解液（Life Technologies，美国）；Ribo⁃Zero rRNA Removal Kits（Illumina，美国）；TruSeq Stranded Total RNA Library Prep Kit（Illumina，美国）；SuperScript Ⅲ逆转录试剂盒（Invitrogen，美国）；qRT⁃PCR 引物（上海云序生物工程技术有限公司）；qPCR SYBR Green master mix（上海云序生物科技有限公司）；RNaseR（epicentre，美国）；RNase inhibitor（Millipore，美国）。NanoDrop ND⁃1000 仪器（Thermo Fisher Scientific，美国）；BioAnalyzer 2100 仪器（Agilent Technologies，美国）；Illumina HiSeq 4000测序平台（上海云序生物科技有限公司）；实时荧光定量PCR 仪（Thermo Fisher，美国）；Gene Amp PCR Sys⁃tem 9700（应用生物系统公司，美国）。

表1 高通量测序口腔白斑组织的临床特点Table 1 The clinical characteristics of patients with oral leukoplakia for sequencing

1.3 RNA 提取及质控

使用TRIzol 试剂盒提取总RNA，使用NanoDrop ND⁃1000 仪测量每个样品的RNA 浓度。以吸光光度（OD）260/280值作为RNA 纯度指标。使用Ribo⁃Zero rRNA Removal Kits 移除总RNA 中的rRNA。采用TruSeq Stranded Total RNA Library Prep Kit 预处理RNA，构建测序文库。采用BioAnalyzer 2100 仪器进行文库质控和定量。根据Illumina 测序说明，将10 pM 文库变性为单链DNA 分子，在Illumina flowcell上捕获，原位扩增为簇，并在Illumina HiSeq 测序仪上采用双端模式（PE mode）进行150 cycle 测序。

1.4 circRNA 测序数据分析

通过Illumina HiSeq 4000 测序仪测序，并使用circBase 数据库和circ2Traits 对所鉴定的环状RNA进行注释［12⁃13］。参考倍数变化（fold change，FC）和P 值进行筛选，选择FC ≥2.0，P ＜0.05 作为筛选差异circRNA 的阈值。circRNA 高通量测序由上海云序生物有限公司完成检测及数据初步分析。

1.5 qRT⁃PCR、酶耐受及sanger 测序验证

根据高通量测序结果，选取表达差异程度较显著的circRNA（FC 越大，P 值越小越佳；原始信号表达值越平均越佳），并结合功能分析结果，共选取10 个circRNA，采用测序的6 对样本进行qRT⁃PCR 验证。对提取的RNA 使用SuperScript Ⅲ逆转录试剂盒逆转录为cDNA 文库。选取GAPDH 为内参，使用ViiA 7 实时荧光定量PCR 仪，采用qPCR SYBR Green master mix 试剂盒进行qRT⁃PCR。本研究所用引物序列由上海生工提供，引物序列见表2。采用2⁃ΔΔCt法对目的基因的相对表达水平进行计算分析，得到各目的基因差异表达情况。根据PCR 选取circRNA 进行酶耐受实验，使用RNase处理提取的RNA，37 ℃水浴20 min，结果合格者进行sanger 测序验证。另选OLK 与NOM 的样本组织各20 例，对筛选合格circRNA circHLA⁃C 进行qRT⁃PCR 进行验证。

表2 qRT⁃PCR 引物信息Table 2 The primers used for qRT⁃PCR experiments

1.6 ceRNA 网络的构建及功能富集分析

通过TargetScan 和miRanda 预测目的circRNA下游的miRNA、mRNA，利用生物信息学软件Cyto⁃scape（v2.8.0）构建ceRNA 网络图。

对差异表达基因进行GO 富集分析和KEGG 功能通路分析，预测差异circRNA 的功能和参与的功能通路。GO 分析主要归类为三部分，包括：生物学过程、细胞成分和分子功能。GO 分析是根据差异表达circRNA 的来源基因和GO 注释列表对GO条目进行筛选。

1.7 统计学方法

数据分析采用SPSS 19.0 和Graphpad Prism 8.0完成。计量资料符合正态分布，用±s 表示，两组间比较采用独立t 检验。使用ROC 曲线分析circH⁃LA⁃C 诊断的特异性与灵敏性，并用spearman 等级相关分析circHLA⁃C 与OLK 异常增生程度相关性。P＜0.05 为差异具有统计学意义。

2 结 果

2.1 OLK 与NOM 间差异circRNA 表达谱及表达特点

聚类分析图显示了OLK 与NOM 间的差异cir⁃cRNA 表达谱（图1a）。根据统计标准，绘制散点图和火山图（图1b、1c）。在OLK 中，共鉴定出366 个显著差异的circRNA，其中65 个上调，301 个为下调。在这366 个circRNA 中，发现28 个新的cir⁃cRNA，其余的circRNA 在circbase 中已有过报道，且差异表达circRNA 大多数为外显子circRNA（图2a、2b），提示差异circRNA 具有重要的生物学功能。circRNA 的长度集中在小于500 bp 和500～2 000 bp 两组内（图2c）。这些差异表达circRNA 几乎分布在所有染色体上，包括22 条常染色体和X染色体（图2d），提示差异circRNA 分布较为广泛。

Figure 1 Expression profiles of circRNAs in OLK versus NOM tissues图1 OLK 与NOM 组织中circRNA 的表达谱

未经处理和分析的测序数据经标准化后，已上传至GEO 网站，GEO 序列号为GSE131568。

2.2 进一步验证目的circRNA

选取的10 个circRNA 包括4 个表达上调和6 个表达下调的circRNA（表3）。经验证，在这10 个circRNA 中有7 个qRT⁃PCR 结果与测序结果保持一致，且具有统计学意义（P ＜0.05，图3）。由于cir⁃cRNA 对核酸酶不敏感，通过qRT⁃PCR 验证的7 个circRNA 中选择了5 个疾病相关性高的circRNA 进行酶耐受实验。结果表示，circHLA⁃C、circPLIN4（perilip3⁃4）和circRNF13 可抵抗核糖核酸酶的水解（图3c），对这3 个circRNA 进行Sanger 测序。结果显示，差异表达最显著的circHLA⁃C 是真正的circRNA，且具有特定反式剪接位点。进一步对circHLA⁃C 进行了可视化分析，circHLA⁃C 由8 个外显子和7 个内含子组成，且反剪接位点为CTACCT⁃GG（图4）。

Figure 2 Distribution of the characteristics of significantly dysregulated circRNAs图2 显著表达差异circRNA 分布特点

表3 OLK 组织10 个显著表达差异circRNATable 3 Ten circRNAs differentially expressed in OLK

Figure 3 Ten circRNAs differentially expressed in OLK图3 OLK 组织10 个显著差异表达circRNA

2.3 OLK 中circRNA 的 功 能 富 集

在OLK 表达上调的circRNA 中，circHLA⁃C 的来源基因HLA⁃C 富集到多达95 个GO 功能条目，分别包括44 条生物学过程、45 条细胞组分和6 条分子功能通路（见本文OSID 码）。HLA⁃C 在生物学过程中可参与调控T 细胞介导的细胞毒性和免疫功能。在细胞组分功能富集分析中，发现HLA⁃C可调控Ⅰ类人类主要组织相容复合体（major histo⁃compatibility complex，MHC）。

通过来源基因的相关通路，进行KEGG 通路富集分析，预测circRNA 功能通路。与366 个差异基因相关的KEGG 通路共有139 条，其中富集到上调的circRNA 通路有60 条。OLK 作为一种癌前病变，在上调通路中共预测到3 条KEGG 通路参与癌变过程。肿瘤相关蛋白聚糖、肿瘤相关miRNA 和肿瘤相关通路都可参与调节胞外信号蛋白激酶/促分裂原活化蛋白激酶（extracellular⁃signal⁃regulated protein kinase/mitogen⁃activated protein kinase，ERK/MAPK1），二者已被证明是OLK 发病的关键因素［14⁃15］。301 个下调circRNA 参与的79 条KEGG 通路，其中显著调节的有泛素介导的蛋白水解、调节干细胞多能性的信号通路和mTOR信号通路（图5）。

2.4 ROC曲线分析circHLA⁃C对OLK的诊断意义

通过ROC 曲线分析评估circHLA⁃C 对OLK 的诊断意义。circHLA⁃C 曲线下面积（Area Under Curve，AUC）为0.955（95% CI：0.892～1.00，P＜0.001）（图6）。提示circHLA⁃C 具有成为OLK 生物标志物的可能性。

2.5 circHLA⁃C 与OLK 异常增生的相关性

Figure 4 Ring formation mechanism and back spliced site of circHLA⁃C图4 circHLA⁃C 的成环机制及反剪接位点

分析20例OLK组织circHLA⁃C的表达量与异常增生程度的相关性，结果表明circHLA⁃C 表达量与OLK轻、中度异常增生呈正相关（r＝0.478，P＝0.033）。

2.6 ceRNA 网络分析

本研究共筛选出7 个circRNA 建立ceRNA 网络（图7）。根据结合能力，列出了每个circRNA 调控的前5 个miRNA。生物信息学继续分析下游的miRNA⁃mRNA 网络，图中呈现了前5 个miRNA 能在下游抑制的mRNA。这表明circRNA 具有十分广泛的作用空间，可通过多种不同的通路对下游分子甚至是蛋白进行调控。

Figure 6 ROC curve analysis of circHLA⁃C in OLK图6 circHLA⁃C 在OLK 中的ROC 曲线分析

3 讨 论

OLK 是最常见的口腔黏膜潜在恶性病变，研究表明OLK 发生癌变平均需要4 至5 年［16⁃17］，因此，OLK 的早期诊断和治疗非常重要。circRNA 曾被认为是分子内错剪接的副产物［18］。近年来，cir⁃cRNA 得到越来越多科研界研究者的关注，目前circRNA 已在真核细胞内确认是一类新的广泛存在、丰度高、保守性高且较为稳定的内源性非编码RNA［19⁃21］。研究表明circRNA_100290，circDOCK1和circ_0007059 可通过ceRNA 机制调控口腔鳞状细胞癌（oral squamous carcinoma cell，OSCC）的生长和 转 移［22⁃24］。hsa_circ_0005320，hsa_circ_0067531和hsa_circ_0000869 可调控口腔黑色素瘤的形成和转移［25］。

circRNA 可作为顺式作用元件调控其亲本基因的表达，促进基因转录［26］。本研究通过circRNA来源基因的功能预测circRNA 的功能。HLA⁃C 表达的变化可影响免疫系统对感染性和自身免疫性疾病的应答效果，且HLA⁃C 的进化使其与杀伤细胞免疫球蛋白样受体的相互作用更加强烈，从而影响自然杀伤细胞反应的调节［27］。研究表明KIR2DL1（+）⁃HLA⁃C2（+）基因型可能是OSCC 早期易感人群的危险因素［28］。

GO 功能富集分析显示差异表达circRNA 具有丰富的功能，多数差异表达circRNA 都可富集到细胞周期过程。推测调控细胞周期的基因表达的变化是OLK 发生的主要原因之一。HLA⁃C 可富集到功能T 细胞介导的细胞毒性和免疫、MHCⅠ类蛋白复合体。研究表明，MHCⅠ类分子与特异性T 细胞受体和CD8+细胞表面分子的识别有关。MHCⅠ类分子是OLK 中的关键分子，在OLK 中呈低表达，且与异常增生程度相关［29⁃30］。此外，OSCC 和伴异常增生OLK 中CD8+细胞数量明显多于无异常增生的OLK［31］。

KEGG 通路富集分析显示在60 种上调的cir⁃cRNA 通路中，富集分数排名前10 通路包括自然杀伤细胞介导的细胞毒性、B 细胞受体信号通路、EB病毒感染、蛋白多糖在癌症中的作用、T 细胞受体信号通路和移植物排斥等。在这些途径中，circH⁃LA⁃C 的靶基因参与了自然杀伤细胞介导的细胞毒性、EB 病毒感染和移植物排斥。这一结果提示OLK 可能与感染和免疫抑制因子有关，而circHLA⁃C 可能是其中的关键因素。与NOM 组织样本相比OLK 中circHLA⁃C 显著升高。ROC 曲线表明circH⁃LA⁃C 对OLK 的诊断具有较高的特异性和灵敏性，且circHLA⁃C 的表达量与轻、中度异常增生程度呈正相关。揭示了circHLA⁃C 可能成为新的生物标志物。

研究表明circRNA 可作为ceRNA 作用于细胞调控网络［32］。本研究构建的ceRNA 网络显示circHLA⁃C 可通过海绵机制与miR⁃9、miR⁃339⁃5p、miR⁃29a⁃3p、miR⁃26a⁃5p、miR⁃222⁃3p 进行互补配对，通过TargetScan、MiRanda 软件可预测circHLA⁃C 与以上miRNA 之间具有结合位点。有研究发现OLK 和OSCC 患者血清miR⁃9 表达较NOM 标本明显降低［33］。在OLK 和一些原发癌组织中，miR⁃339⁃5p 表达下调，作为肿瘤抑制因子参与癌症进展［34⁃35］。miR⁃29a⁃3p 和miR⁃26a⁃5p 在OLK 和癌组织中显著下调［36］。此外，与正常和OSCC 患者相比，OLK 中miR⁃222⁃3p 显著下调［37］。由此可见，circHLA⁃C 可能具有同一miRNA 的多个结合位点或多个不同miRNA 的结合位点。circHLA⁃C 可能就是通过海绵吸附这些miRNAs 来调控OLK 的发生发展。对于circHLA⁃C 与miRNAs 及OLK 下游靶基因的相互作用，今后还需要进行更深入的研究。

综上所述，本研究探讨了circRNA 在OLK 病因病理及发生发展中的作用，并揭示了circHLA⁃C 具有成为OLK 潜在分子标志物及治疗靶点的可能。

【Author contributions】Xu SM wrote and revised the article. Song YH, Shao YH, Tao L directed the data collection and analysis. Zhou HW designed the study. All authors read and approved the final manu⁃script as submitted.