肝细胞癌模型小鼠中外泌体来源RNA剪接基因的功能及差异表达分析

阮文慧，唐玉莲，倪安妮，李根亮，李曙波，陆瑞群，黄海玲

（右江民族医学院生物化学与分子生物学教研室，广西 百色533000）

肿瘤微环境（tumor microenvironment，TM）是肿瘤的内部生存环境，肿瘤前微环境（pretumor microenvironment，PTM）是健康生理微环境（physiological microenvironment，PM）到TM之间的过渡阶段［1］。与健康组织相比，TM中出现了大量的异常表达基因［1-2］。外泌体能通过其携带的基因参与多种生理、病理过程，TM中的外泌体与癌症的发生发展密切相关［3-4］。然而，外泌体来源的RNA剪接基因（exosome-derived RNA-splicing genes，ERGs）在模型小鼠的肝细胞癌微环境中发生发展的具体生理途径仍少有报道。

RNA剪接调控着细胞周期的进展，并在癌症和多种疾病中起着重要作用［5］。mRNA前体的选择性剪接是一种基本的调控机制。经过选择性剪接可产生多个不同的mRNA转录本，生成不同甚至相反功能的蛋白质异构体［6］。真核细胞中最丰富的核糖体RNA与mRNA剪接位点上的序列互补，经RNA剪接酶的识别，切割和剪接RNA。剪接体的磷酸化可能通过调节RNA剪接参与HCC（hepatocellular carcinoma，HCC）的转移［7］。RNA聚合酶相关蛋白则能参与mRNA前体的选择性剪接并诱导人肝癌细胞凋亡［8］。可见RNA剪接与癌症息息相关。

为探究这一问题，本研究以PM作为对照，通过构建肝癌模型小鼠，分别提取PM、PTM和TM三组样本中的总RNA进行RNA测序（RNA-seq）和生物信息学分析，并构建外泌体来源RNA剪接基因作为比对数据库与测序数据比对，筛选差异表达的外泌体来源RNA剪接基因（ERDs）和分析ERDs的功能及异常表达的可能原因。研究结果从分子机制的角度，为HCC发生发展过程中外泌体来源RNA剪接基因参与HCC的病理机制提供基础资料。

1 材料和方法

1.1 细胞培养H22肝癌细胞系（购自南京凯基生物技术有限公司）经体外细胞培养至指数增长期，3 000 r·min-1离心、洗涤、稀释至终浓度105个·mL-1时用于构建肝癌模型小鼠。

1.2 肝癌模型小鼠的构建及样本获取100只健康昆明小鼠（鼠龄7周左右，体重22～25 g）随机分成两组，即实验组小鼠和阴性对照组小鼠（PM组）。实验组小鼠（n=60）给予每只双前腋下注射指数生长期时的H22肝癌细胞系0.2 mL（隔天注射一次，重复3次），其余的小鼠（n=40）注射不含H22的生理盐水同法注射作为对照，分别正常饲养4周后分三组样品进行收集，实验组小鼠中成瘤小鼠分别取瘤体和肝未癌变部分为肝癌模型小鼠的TM和PTM，对照组小鼠中取健康肝组织细胞作为健康小鼠的肝组织（PM）［1-2］。实验动物的使用符合国家及右江民族医学院实验动物伦理有关章程，在实验实施时严格执行实验动物安全制度及相关规章。

1.3 RNA提取和测序各样本内提取总RNA（n=20），同组内等量取各样本总RNA混合成混合总RNA样本［1-2］，提取方法按照本实验室之前的实验方法进行，并送公司测序。RNA质量检测、逆转录、文库构建、RNA-seq等按常规程序完成［1］。

1.4 测序数据差异表达基因分析以小鼠RNA基因组为背景对测序结果进行比对、注释和较正，统计分析各基因RPKF（Reads Per Kilo bases per Million reads）值的差异，同一基因不同样本间的RPKF比值≥2或≤0.5，且错误发现率（false discovery rate，FDR）≤0.05的 基 因 作 为 差 异 表 达 基 因（differentially expressed genes，DEGs）。

1.5 ERDs的筛选使用NCBI数据库（https：//www.ncbi.nlm.nih.gov/）以（"exosome"or"exosomal"）and（"splicing"or"spliceosome"or"spliceosomal"）为检索词，在NCBI中的Nucleotide数据库中进行检索。将检索结果通过DAVID在线转换成Ensembl gene ID和Official gene symbol后，将转换结果与测序数据中的DEGs进行比较，筛选出样本间差异表达的外泌体源RNA剪切基因（exosome-derived RNA-splicing related DEGs，ERDs）。

1.6 ERDs的功能分析利用DAVID 6.8在线软件（https：//david.ncifcrf.gov）对21个高表达的ERDs进行GO（Gene Ontology）和KEGG（Kyoto Encyclopedia of Genes and Genomes）分析，并对富集到的生物学程序（biological process，BP）、细胞组分（cellular component，CC）、分 子 功 能（molecular function，MF）和信号通路（KEGG pathway）等进行功能注释聚类（functional annotation clustering），其中以FDR≤0.05为差异有统计学意义。

1.7 ERDs编码蛋白的互作分析通过STRING 11.0版本数据库（https：//string-db.org/）在线对ERDs编码蛋白进行Co-expression互作分析，以评价其在表达上的相关性及其中发挥功能的关键分子。

1.8 ERDs的AS和SNV/INDEl与基因表达的相关性分析用Excel 2010和SPSS统计分析，对外泌体相关的ERDs中出现的AS（alternative splicing，AS）和单核苷酸位点变异（single nucleotide variants，SNV）及插入缺失标记（insertion/deletion，INDEL）进行卡方检验，其中P＜0.05且两样本间发生频率的倍数不小于2的为差异有统计学意义上。

1.9 RT-qPCR验证ERDs随机选择3个在三组样本都表达的ERDs进行RT-qPCR分析，以验证基因表达量和RNA-seq测序结果准确性，同一样本间和同一组不同样本间的实验重复3次。PCR条件参照文献，引物采用在线引物设计软件（https：//www.ncbi.nlm.nih.gov/tools/primer-blast/）进行设计，并由上海生工进行合成［1-2］。RT-qPCR实验内参基因为18sRNA。

2 结果

2.1 测序数据TM组、PTM组及PM组通过RNAseq分别获得了5.01 G base（35.28 M clean reads）、5.06 G base（35.36 M clean reads）和5.02 G base（35.32 M clean reads）的数据。三组样本所得clean reads对基因组的匹配率分别为92.65%、90.10%和89.20%，且三样本测序数据的Q20皆大于99%。

2.2 ERDs NCBI的nucleotide数据库中共检索到28 275个与小鼠外泌体RNA以及RNA剪接相关的核苷酸序列。经DAVID 6.8软件共成功注释到12 273个ERGs作为研究背景。在TM和PM之间总共有4 533个测序数据组差异表达基因（DEGs），与DEGs匹配的ERDs有151个，其中10个下调、141个上调。在TM和PTM之间总共有4 042个DEGs，筛选出的ERDs有169个，其中11下调、158个上调。两组对比样本中共有21个相同的ERDs都在TM中高表达。

2.3 ERDs的功能以上筛选出的21个相同的ERDs功能富集分析结果显示，它们主要富集到4个BP、5个CC、2个MF和1个KEGG pathway（图1）。它们主要参与剪接体的功能执行，尤其是在剪接体执行剪接过程中ploy（A）RNA结合、核来源的第二步剪接小体的催化作用以及mRNA前体的剪接和剪接体snRNP组装。结果表明，它们可能通过剪接体调控的信号通路（图2）促进HCC发生发展。

图1 ERDs的功能

图2 ERDs参与剪接体的过程（红星表示该基因在TM中表达上调）

2.4 ERDs编码蛋白的相互作用21个上调的ERDs编码蛋白的相互作用见图3。

图3 ERDs编码蛋白的互作关系

2.5 ERDs的AS和SNV/INDEL 21个在TM中高表达的ERDs的AS分析结果显示，Hnrnpa1和Lsm3两个基因发生了AS，且其与这两个基因一样，在TM中都上调表达，而在PTM和PM中都下调表达。可见，AS的发生可能与相应基因的表达呈正相关。

SNV和INDEL分析显示，TM与PM对比，差异基因发生了110个SNV和3个INDEL位点的改变；TM与PTM对比，差异基因发生了112个SNV和4个INDEL位点的改变。且两组对比样本的基因位点的改变都主要发生在基因的6个部位，即外显子区（exonic）、内含子区（intronic）、基因间区（intergenic）、基因下游区（downstream）、基因上游区（upstream）和剪接区（splicing）（图4）。在这些核苷酸位点发生改变的基因中，TM与PM对比中有16个RNA剪接相关的差异基因的SNV改变具有统计学意义，其中有12个ERDs的SNV改变与其基因的表达呈正相关（表1）。TM与PTM中有16个差异基因SNV改变有统计学意义，其中11个SNV改变与其基因的表达也呈正相关（表2）。

在TM与PM对比中，11个ERDs在TM中出现了50个SNV现象，但在PM无此位点改变（表3）。在TM与PTM对比中，12个ERDs在TM中出现了40个SNV现象，但PTM中没有出现这些位点改变（表4）。

图4两组对比样本中ERDs的SNV和INDEL在基因不同部位上发生的频率

表1 TM/PM中ERDs的SNV发生频率与基因表达的相关性

表2 TM/PTM中ERDs的SNV发生频率与基因表达的相关性

2.6 RT-qPCR结果RT-qPCR结果显示（n=3），随机选择的3个ERDs的表达结果与RNA-seq结果一致（图5）。引物序列见表5。

表3 TM/PM中ERDs的SNV发生情况

表4 TM/PTM中ERDs的SNV发生情况

图5 ERDs在TM/PTM以及TM/PM中的相对表达情况（与对照组相比）

表5 RT-qPCR引物

3 讨论

在本研究中，我们通过从NCBI上找到RNA剪接有关基因和小鼠外泌体RNA基因构建背景基因，利用背景基因与高通量测序结果进行比对，筛选得到了21个在TM中高表达的外泌体源RNA剪接相关基因（ERGs）。这些差异基因包括：Fancd2、Ctnna1、Tnfrsf1b、DNAJB1、dis3、KIF23等，它们主要参与了催化作用的步骤2剪接小体、剪接体snRNP组装等生物学过程，这些生物学过程可能通过poly（A）RNA结合的方式以剪接体通路来促进HCC的发生发展。

真核细胞中剪接体执行的过程是基因表达的重要环节，一些疾病的起因可由剪接体蛋白的突变而影响转录本的剪接，甚至一些癌症也与剪接因子的错误调控有关［9］。例如：妇科恶性疾病如宫颈癌，剪接因子对常见Fancd2基因结构异常起重要作用，引起了复发患者中Fancd2基因突变率高，Fancd2基因突变形式通常为基因组的缺失，最终导致了Fancd2基因的低表达或功能缺失，导致宫颈癌预后不良易复发［10］；Hdm2剪接体与结肠癌细胞的增殖相关，其通过下调抑癌基因表达，促进结直肠癌的发生和发展［11］；在t（8，21）急性髓系白血病中，选择性剪切体AML1-ETO9a与细胞的恶性程度息息相关［12］；LSM剪接体蛋白过度表达导致肺癌细胞选择性剪接增多［13］。在本研究结果中，LSM剪接体蛋白基因中的Lsm2和Lsm3，在HCC发生发展过程中出现上调表达。Lsm2和Lsm3编码蛋白位于细胞核中起催化作用的第二步剪接小体中，通过poly（A）RNA结合，导致mRNA前体剪接异常，不能成为成熟的mRNA。并且Lsm2和Lsm3编码蛋白又与细胞内的核糖核蛋白复合体密切相关，核糖体的异常导致mRNA无法通过核糖体的作用转变成蛋白质行使生命活动的各种功能。

另外，RNA-结构元件在mRNA前体的剪接催化中起着关键作用，然而RNA催化结构的形成需要剪接体蛋白的协助。Rbm22作为一种剪接体蛋白在促进剪接体催化RNA元件活性构象中的具有潜在作用［14］。而Snrpd1参与催化mRNA前体剪接、调节多功能特异性剪接体组装及获得和维持多潜能剪接，与许多免疫性疾病有关，如：系统性红斑狼疮［15］。本研究中，催化作用的第二步剪接小体中的Rbm22和Snrpd1上调表达，使剪接体snRNP组装异常，从而影响poly（A）RNA结合，并导致mRNA在细胞质中的稳定性降低，mRNA被过早降解。这些基因都来源细胞核，且相互关联而影响剪接体通路。我们的研究中的通路图和蛋白互作图表明了这一点。

已知mRNA表达过程中的AS事件会产生不同的蛋白质，并且这种情况也发生在肝癌的发展过程中［16］。我们的AS分析结果发现，TM样本中Hnrnpa1和Lsm3内AS和基因表达均上调表达，也同样证明了这一点。可见，引起基因多态性的遗传因素AS与基因表达呈正相关。另外，SNV分析显示SNV改变与RNA剪接基因表达量具有很强的相关性，TM中上调表达RNA剪接基因多具有SNV现象，但在PTM和PM中无此SNV现象；基因位点突变和选择性剪接在健康组织和癌旁健康组织中并无发生，但在癌症发生后，由于RNA剪接基因上的这些基因位点的突变导致剪接体出现变异。因此，RNA剪接基因中的AS以及SNV突变可能导致异常剪接和poly（A）RNA异常结合，在HCC发生发展中起着重要作用。

综上所述，在肝癌细胞的发生发展过程中，在TM样本中上调表达21个ERGs影响了参与poly（A）RNA结合的Lsm2、Magoh、Eftud2、Rbm22等因子，导致剪接体通路出现异常，进而导致mRNA转录本的稳定性和mRNA的翻译效率受到影响。其中Pabpc1、Ppil1、Lsm2等还出现异常上调表达，使核中步骤2剪接小体异常催化；Snrpd1、Smn1、Gemin5和Rbm22等异常上调表达导致剪接体snRNP组装出现异常，这造成附加因子和剪接装置的错误剪接而参与HCC的发生发展，且这种异常剪接与ERDs自身的基因多态性改变有关。