环状RNA是一类新发现的非编码RNA,主要通过miRNA“海绵”机制
[1]发挥转录后基因表达调节作用
[2, 3]。其表达失调与恶性肿瘤的发生、发展密切相关
[4, 5]。既往研究表明环状RNA可作为部分肿瘤诊断和预后的标志物
[6, 7],是一个具有前景的肿瘤治疗靶点
[8, 9]。随着测序技术的发展加深了对环状RNA的理解,并促进其生物发生和作用机制的深入研究
[10-12]。 然而,在肾母细胞瘤(WT)中环状RNA表达谱及其进展的作用机制尚未阐明
[13, 14],同时由于WT缺少用于临床早期筛查、鉴别诊断和预后评估的分子标志物
[15-17]。因此,为了系统描绘环状RNA在WT中的表达模式和潜在作用。本研究基于临床WT样本进行高通量环状RNA测序,鉴定在WT中差异表达的环状RNA分子,首次展示WT中环状RNA表达谱。随后,基于临床来源的肿瘤样本队列进行表达分析,同时搜集临床队列的病理资料和随访结局,分析环状RNA分子与临床病理特征的相关性和对预后的影响,提示环状RNA在WT中的潜在作用,为WT寻找新的治疗靶点提供依据。
1 资料和方法
1.1 临床样本搜集
肾母细胞瘤肿瘤组织和邻近对照正常组织取自重庆医科大学附属儿童医院泌尿外科。所有标本未经过术前化疗且均由儿院病理科确诊。肿瘤的临床分期和病理分型依据COG方案。选择4对配对的肾母细胞瘤组织进行高通量环状RNA测序。其余搜集的肿瘤样本储存在液氮和多聚甲醛中用于后续实验分析。总结纳入本研究患儿的临床资料包括:性别、年龄、发病侧别、肿瘤体积、瘤体是否破裂、COG病理分型、术后分期和淋巴结是否转移等,同时进行常规随访记录患儿是否有肿瘤复发或死亡事件。本研究方案通过重庆医科大学附属儿童医院伦理委员会的批准(伦理批号:2022年伦审第50号)。
1.2 环状RNA测序
环状RNA测序由上海鲸舟基因科技有限公司完成。主要流程如下:首先使用RiboZero rRNA去除试剂盒(Epicenter)处理4对配对的肾母细胞瘤组织的总RNA样品。随后,将rRNA耗尽和RNase R酶切后的RNA样本进行片段化,用随机引物合成cDNA。纯化cDNA的PCR扩增产物,随后建库并进行质量控制。使用Novaseq6000仪器(Illumina,USA)进行双端测序,因此每个示例包含两个FASTQ文件。使用FASTQC软件评估每个FASTQ文件的测序质量。使用FASTP软件
[18]对序列进行过滤,获得可用于定量分析的有效数据“clean reads”。参考基因组(GRCH38/hg38)从UCSC基因组浏览器下载(
http://genome.ucsc.edu/)。使用ANNOVAR数据库将每个环状RNA的剪接端映射到基因组区域,然后将结果与参考序列和UCSC数据库进行比较,以注释环状RNA的功能元件。随后,采用CIRI软件
[19]对环状RNA进行预测和鉴定。在以上预测的基础上,根据环状RNA在染色体上的位置信息,合并全部样本中环状RNA的鉴定结果,对合并后结果进行重新编码ID。而后将其与circBase数据库
[20]进行匹配,从而识别出哪些是已知以及哪些是新发现的环状RNA分子。Circos软件用于构建环状RNA测序结果全貌的Circos图
[21]。对测序结果进行质控分析,测序数据质控结果如下
表1所示。数据量约15G/样本,每项碱基质量大于20(Q20)的比例不小于90%。所有的样品的测序量都满足标准,可以进行数据分析。
1.3 生物信息学分析
1.3.1 环状RNA的表达定量和差异表达分析
无论是已被数据库注释的还是新预测的环状RNA分子,目前大部分环状RNA都无法获得其完整序列。只能使用环状RNA序列上的“back-spicing”即反向剪切位点部分的“junction reads”数量来进行环状RNA定量
[22]。通常使用每十亿映射的剪接读数(SRPBM)值来定量环状RNA的表达水平
[23],其计算公式如下:
。环状RNA在肿瘤组织和配对正常组织间的差异表达趋势使用差异倍数(FC)表示,使用“edgeR”包进行差异表达分析;然后基于log2FC和
P值进行“火山图”展示,差异表达筛选标准定义为|log2FC|>2和
P值<0.05。
1.3.2 环状RNA泛癌表达分析
CIRI-HUB是最新开发的环状RNA可视化分析平台(
http://circatlas.biols.ac.cn/)。该平台整合了33种不同类型的肿瘤组织包括2187个肿瘤样本和680个正常样本的环状RNA表达谱。同时,该平台可用于对肿瘤和正常组织中环状RNA表达的个性化分析,并自动生成图形结果。
1.3.3 临床相关性分析
为了验证测序结果的有效性,首先对环状RNA表达谱进行“z值”均一化,进一步使用“prcomp”函数进行降维分析以获得降维后的矩阵。然后使用“stats”主成分分析(PCA)数据并使用“ggplot2”包绘图,以可视化环状RNA表达矩阵对肿瘤和正常组织的区分能力。此外,使用Kaplan-Meier生存曲线评估单个环状RNA分子对肾母细胞瘤预后的预测能力,使用受试者工作特性曲线(ROC)来评估单个环状RNA分子对肿瘤组织和正常组织的鉴别效能。为了分析单个环状RNA分子与临床病理特征的相关性,我们根据目的环状RNA的表达中位数将患儿分为高表达和低表达两组,随后采用卡方检验或t检验比较组间临床变量差异。
1.4 q-PCR实验
细胞RNA提取使用核酸快速提取试剂盒(BioFlux)。使用反转录试剂盒(TaKaRa)进行总RNA逆转录合成cDNA。使用实时荧光定量PCR试剂盒(TaKaRa)进行q-PCR定量分析,反应体系和反应程序根据试剂盒说明进行。每组设置3个复孔,GAPDH内参基因用于待检测基因的标准化,根据公式2
-ΔΔCT计算待检测基因在样本中的相对表达水平。实验过程中所需的引物由北京擎科生物公司设计合成,序列如下:(
表2)。
1.5 环状RNA分子环状结构的验证
1.5.1 RNase R酶消化实验
将RNA样本浓度控制在1500 ng/μL左右,安装试剂盒说明在冰上配制反应体系。使用逆转录PCR仪设置消化程序:37 ℃反应10~30 min→70 ℃灭活10 min;另取等量RNA样品(3 μg)加入DEPC水补充体系至20 μL作为RNase R阴性处理组。
1.5.2 琼脂糖凝胶电泳定性检测
对RNase R阳性处理组和阴性处理组同时进行逆转录反应,逆转录过程使用反转录试剂盒(TaKaRa)合成cDNA,方法与前述一致(逆转录体系设置为40 μL)。本部分实验使用可用于直接电泳的逆转录试剂盒(TaKaRa),安装试剂盒说明配置反应体系。使用1×TAE电泳缓冲液制作1%的琼脂糖凝胶,在微波炉中加热沸腾至澄清状态;冷却后缓慢倒入已装好梳子的磨具中,冷却至定性;每个样品上样量为5 μL;电泳完毕后,采用凝胶成像系统拍照保存。
1.5.3 q-PCR定量检测
反应体系和反应程序与前述实验方法一致。上机检测RNase R阳性处理组和阴性处理组内环状RNA、亲本基因和内参基因GAPDH的绝对表达水平。由于两组间RNA总量一致,因此将RNase R阴性处理组的GAPDH表达作为公共内参基因进行表达均一化,最后比较RNase R阳性处理组和阴性处理组内相关基因的表达差异。
1.6 统计学分析
分类变量以n(%)表示,连续变量使用均数±标准差表示,分别使用卡方检验和t检验进行统计学分析。使用Graph Pad Prism 9.5.0软件分析实验数据。所有的生物信息学分析使用R软件(版本4.0.3)。P<0.05表示差异具有统计学意义。
2 结果
2.1 基于高通量测序鉴定和注释肾母细胞瘤中的环状RNA
搜集4个配对的临床肿瘤组织,进行高通量环状RNA测序(
图1A)。将用于环状RNA鉴定的“clean reads”映射到参考基因组,共检测到23,978个环状RNA分子(
图1B)。染色体分布分析显示环状RNA广泛分布在1号染色体和10~19号染色体上。与环状RNA数据库进行比对,其中10,884个环状RNA已获得注释(
图1C)。
2.2 环状RNA的分类和差异表达模式
基于亲本基因序列来源,大多数(84.18%)环状RNA由外显子构成,11.72%由内含子构成,剩余4.10%映射到基因间区域(
图2A)。验证测序结果的有效性,针对环状RNA表达谱进行主成分分析,测序结果清晰地将肿瘤和正常组织区分开来(
图2B)。根据筛选参数(|log2FC|>2和P值<0.05),同时人工去除表达丰度为0的样本,初步鉴定出614个差异表达的环状RNA分子,其中包括269个上调和345个下调(
图2C)。
2.3 临床样本验证提示hsa_circ_0001900在WT中呈现出肿瘤特异性表达模式
根据环状RNA测序结果进行表达丰度排序(
图3A)。选择表达丰度最高的前6个环状RNA(环状RNA注释见
表3)在12对配对临床样本中进行验证,其中hsa_circ_0001900表达丰度最高且在临床来源的肿瘤组织中与测序结果表达趋势一致(
图3B)。诊断ROC曲线显示该分子具有良好的鉴别诊断能力,曲线下面积为0.72(
图3C)。
2.4 临床相关性分析提示hsa_circ_0001900的高表达与较大的肿瘤体积和较差的预后有关
hsa_circ_0001900在大部分恶性肿瘤组织中呈现出高表达(
图4A),在一个较大的临床样本队列中分析该分子的表达水平。差异表达分析显示该分子在肿瘤组织中特异性高表达(
图4B)。根据表达量中位数将该分子分为高低表达两组,临床相关性分析(
表4)显示高表达组患儿肿瘤体积明显更大(
P=0.025)。相关性分析(
图4C)显示hsa_circ_0001900的表达水平与肿瘤大小呈正相关(R=0.35,
P=0.040)。结合随访数据的生存分析(
图4D)显示高表达亚组具有较低的生存期(4年无复发生存率分别为81.25%
vs 85.11%,
P=0.496),但是差异无统计学意义,此结果考虑与临床样本量和随访时间有关,图中统计结果为初步发现, hsa_circ_0001900可能是WT进展的关键驱动基因。
2.5 环状RNA分子hsa_circ_0001900环状结构的注释和验证
根据环状RNA测序结果,hsa_circ_0001900是由亲本基因CAMSAP1的外显子2和外显子3通过“反向剪接”形成,包含425个碱基的共价闭环结构(序列如下):
>hsa_circ_0001900(chr9:135881633-135883078)
→ATAACATCCCTGAGGACCTCAGAGACCCTTTCTACG TTGACCAGTATGAGCAGGAGCACATTAAGCCGCCTGTTATCAAGCTTCTCCTGTCCAGCGAGCTGTACTGCCGTGTCTGCAGCCTCATCCTGAAAGGGGACCAGGTGGCCGCCTTACAGGGACACCAGTCTGTCATCCAGGCCCTGTCCCGGAAAGGGATCTATGTGATGGAGAGTGATGACACCCCCGTGACAGAGTCCGACCTCAGTCGCGCACCCATAAAAATGAGTGCCCACATGGCAATGGTGGATGCCCTGATGATGGCCTACACTGTGGAGATGATCAGCATCGAGAAGGTGGTGGCCAGTGTCAAGCGCTTCTCAACGTTCAGTGCCTCGAAAGAACTTCCGTACGACCTCGAGGATGCCATGGTGTTCTGGATCAACAAG→
该分子的成环剪接位点序列得到Sanger测序的验证,Sanger测序序列与环状RNA测序序列一致(
图5A)。环状RNA具有共价闭环结构,相对于线性基因更加稳定,更耐受RNase R酶的消化作用。采用Rnase R酶消化实验,对消化产物进行琼脂糖凝胶电泳定性分析。相对于线性亲本基因CAMSAP1和内参基因GAPDH,hsa_circ_0001900更加稳定(
图5B)。随后,对RNase R酶消化产物进行PCR定量分析,进一步验证hsa_circ_0001900分子的稳定性能(
图5C)。
3 讨论
儿童肾母细胞瘤的综合治疗措施已取得优异的治疗效果,然而,该肿瘤发生的驱动机制和肿瘤进展的分子机制仍然未完全阐明
[24-26]。先前的研究强调了环状RNA在肿瘤进展机制以及临床治疗中的潜在作用
[27-29]。环状RNA是一类新的内源性非编码RNA,该分子具有共价闭环结构,没有5’端帽子或3’端Poly A尾巴。现有研究表明,环状RNA通常由1~5个外显子通过首尾剪接构成,长度在几百到数千个核苷酸之间
[30]。环状RNA在真核生物中大量表达,并表现出组织或细胞特异性。研究表明,cricRNA含有大量的miRNA结合位点,具有miRNA海绵作用,参与调控STIM1基因的表达并参与肾间质纤维化
[31]。例如:cricWAC通过靶向miR-142,上调WWP1并激活PI3K/AKT途径诱导TNBC的化疗耐药性
[32]。这提示我们hsa_circ_0000190可能通过特定条件海绵机制参与肿瘤的进展。虽然环状RNA的产生和功能机制尚不完全清楚
[33, 34],但最近的研究表明,环状RNA可作为疾病诊断和治疗的潜在分子标记物。
迄今为止,还没有研究报道环状RNA在肾母细胞瘤中的详细全面表达。在本研究中,我们通过高通量测序技术系统地分析了肾母细胞瘤中环状RNA的表达谱,其中绝大部分由亲本基因的外显子环化形成。与已知的环状RNA数据库circBase(
http://www.circbase.org/)相比,在肾母细胞瘤中新鉴定超过一半的环状RNA分子(54.61%)。环状RNA的染色体分布与其他类型肿瘤相似
[35]。根据方法中描述的筛选标准,共鉴定出614个差异表达的环状RNA分子,其中包括269上调和345个下调。为了确定关键的环状RNA分子,基于表达水平和表达差异进行临床样本验证。结果提示新的环状RNA分子hsa_circ_0001900在肾母细胞瘤中高表达,且表达水平最高。泛癌分析显示该分子在多种肿瘤中呈现显著的差异表达,暗示该分子在肿瘤中的重要作用。为了证实hsa_circ_0001900测序结果的有效性。使用RNase R酶消化实验和PCR产物的Sanger测序进一步验证了该分子的“反向剪接”序列和环状结构。
环状RNA的稳定性、保守性和特异性使其成为一种潜在的肿瘤预后和诊断生物标志物。近年来,许多研究表明环状RNA在人体体液中稳定表达
[36],如唾液、血浆和外泌体,这也使环状RNA成为肿瘤液体活检的理想生物标志物。例如,相比非肿瘤组织和健康对照,Circ-ZEB1.33在肝细胞癌组织和血清样本中高表达,同时该分子的表达水平与TNM分期和肝细胞癌患者的总生存率有关,表明Circ-ZEB1.33可作为肝细胞癌患者有价值的预后生物标志物
[37]。hsa_circ_0000190在胃癌组织和血浆样本中表达下调,其表达水平与肿瘤大小、远处转移、淋巴结转移、TNM分期和CA19-9水平显著相关
[38]。在另一项研究中,hsa_circ_0000745在胃癌组织和血浆样本中低表达,临床分析显示胃癌组织中该分子的表达水平与肿瘤分化相关,而血浆中的表达水平与肿瘤淋巴结转移相关
[39]。通过将环状RNA与其他临床生物标志物结合,可以建立新的诊断模型提高诊断的准确性和特异性。血浆中hsa_circ_0000745联合癌胚抗原(CEA)水平的诊断曲线下面积为0.775,提示血浆hsa_circ_0000745联合CEA水平对胃癌有良好的诊断价值。基于以上背景,本研究也分析了hsa_circ_0001900在肾母细胞瘤组织中的表达与临床指标和预后的关联。临床相关性分析显示该分子与肾母细胞瘤体积呈正相关,同时hsa_circ_0001900高表达组患儿相比低表达组患儿的无复发生存率更低。通过文献检索发现,hsa_circ_0001900参与促进鼻咽癌、结肠癌、骨肉瘤、肝细胞癌和小细胞肺癌的进展和转移
[40-44]。以上结果充分提示hsa_circ_0001900可能是肾母细胞瘤进展的关键驱动基因。
综上所述,通过临床样本的测序分析和表达验证,本研究在肾母细胞瘤中鉴定出新的环状RNA分子hsa_circ_0001900。临床相关性分析提示该分子的表达水平与肾母细胞瘤大小和患儿预后有关,提示其对肾母细胞瘤生物学功能的潜在影响。