基于RNA-Seq对心肌肌钙蛋白IR193H突变限制性心肌病小鼠转录组的分析

干意，封叶，赵唯安，刘振国，张蕾，田杰

限制性心肌病(restrictive cardiomyopathy，RCM)是儿童心肌病的一种类型，占小儿心肌病的2%～5%[1]，是一种由于心肌僵硬度升高导致以舒张功能严重受损为主要特征的心肌病，表现为心室舒张末容积正常或缩小，心室壁厚度正常或轻度增加，而收缩功能大多正常或仅有轻度受损[2-3]。儿童RCM多为特发性，病因不明，约30%病例有家族发病倾向，多为常染色体显性遗传。已有研究发现心肌肌小节蛋白、结蛋白、心肌肌钙蛋白I(cardiac troponin I，cTnI)、β-肌球蛋白重链等基因突变与RCM有关[4]，其中，cTnI作为重要的心肌结构蛋白，其突变与原发性心肌病关系尤为密切[5-6]。研究发现，6种类型的cTnI突变可导致RCM(L144Q、R145W、A171T、K178E、D190G、R192H)，其中两种类型的突变(K178E、R192H)临床症状最为严重[7]。本课题组拥有黄旭培教授成功构建的来源于人cTnI基因突变的RCM小鼠动物模型(在小鼠系193位点，cTnI R193H)，该模型小鼠与携带该突变的RCM患者(在人192位点，cTnI R192H)有相似的病理改变、症状和表型[8-9]，其在6～8周时，超声即可检测到舒张功能的改变[10]，但具体致病机制尚不明确。本研究利用高通量测序的方法，分析RCM小鼠与正常小鼠差异表达的基因，筛选出可能参与致病过程的分子及信号通路，为RCM发病机制研究提供新的思路和方向。

1 材料与方法

1.1 主要试剂和设备 R193H(美国佛罗里达大西洋大学黄旭培教授赠送)C57和野生型(WT)C57 3月龄小鼠(购自重庆医科大学实验动物中心，饲养于重庆医科大学附属儿童医院SPF动物中心)，RNA提取试剂盒(BioTeke，北京)，反转录试剂盒(TaKaRa，日本)，SYBR Green试剂盒(KAPA，美国)。实验室常用设备由重庆医科大学儿研所提供。动物实验符合伦理要求，并经重庆医科大学医学研究伦理委员会批准。

1.2 小鼠心脏组织总RNA提取 准备3月龄R193H和WT小鼠，利用CO2麻醉处死后，剪开小鼠胸腔，灌洗心脏，取心脏组织，利用高纯度总RNA快速提取试剂盒提取心脏组织总RNA。总RNA样本质量由上海裂冰生物医药有限公司采用Agilent2200仪器进行检测。RNA分子完整数(RNA integrity number，RIN)＞8.0的RNA用于cDNA文库构建。样本和RNA均放入–80℃冰箱保存。

1.3 基于RNA-Seq技术的差异表达基因筛选基因的表达量计算主要采用RPKM方法(Reads per kb per million reads)，计算公式为：RPKM=106C/(NL/103)。其中RPKM为基因表达量，C为唯一比对到该基因上读段(reads)数，N为唯一比对到参考基因上的总读段数，L为该基因编码区的碱基数。RPKM法消除了基因长度和测序量差异对计算基因表达的影响，计算得到的基因表达量可直接用于比较不同样品间的基因表达差异[11]。采用国际公认算法EBSeq软件包对表达量(Counts)进行差异筛选。以差异倍数(fold change，FC)≥1.2或≤0.833、P＜0.05为标准筛选差异表达基因，其中错误检出率(false discovery rate，FDR)＜0.05。利用Cluster软件，以欧氏距离为距离距阵计算公式，对差异表达基因和样本同时进行分层聚类分析，聚类结果用Java Treeview 1.1.6 R2软件显示。

1.4 基因功能分析 将分析得到的差异基因基于数据库NCBI(http://www.ncbi.nlm.nih.gov/)、UniProt(http://www.uniprot.org/)、GOC(http://www.geneontology.org/)，分别从生物学过程(biological process，BP)、细胞组分(celluar components，CC)、分子功能(molecular function，MF)这三个层面进行基因本体(gene ontology，GO)注释，得到基因参与的所有GO，采用Fisher检验计算每个GO的显著性水平(P值)，从而筛选出差异基因富集的显著性GO。其中，P＜0.05的值用灰白相间的表格标出。

1.5 功能调控网络构建 基于GO的层次结构，利用NovelBio云计算平台将所有GO之间的相互调控及从属关系整理成数据库，构建功能关系网络，总结实验影响的功能群体，以及显著性功能的内在从属关系。以差异基因所做GO分析中的显著性GO条目(P＜0.05)为研究对象进行功能调控分析，构建功能调控网络。

1.6 信号通路分析(pathway-analysis) 将筛选出的差异基因于KEGG(Kyoto encyclopedia of genes and genomes)数据库进行通路注释，得到差异基因参与的所有通路条目，采用Fisher检验计算通路的显著性水平(P值)，从而筛选出差异基因富集的显著性通路条目。P＜0.05的值用灰白相间的表格标出。

1.7 信号通路相互作用关系分析 选用显著性通路，根据KEGG数据库中通路之间的关联关系整合成显著性通路之间的信号转导网络。在此以差异基因所做的通路分析中的显著性通路条目(P＜0.05)进行信号通路相互关系网络构建。

1.8 引物的设计和合成 在NCBI的GenBank数据库(https://www.ncbi.nlm.nih.gov/genbank/)中检索Tifab、Ifi202b、Ctss、Itgb2、Pik3ap1、Arhgap21、Myocd、MT-ND5基因的mRNA序列。采用Primer Premier 5.0软件针对这些基因CDS核心编码区设计特异性引物，引物的碱基序列由上海生工生物工程股份有限公司合成。引物序列为：Tifab，正向5'-GCTGCCACAGTTGTCTCGCTAC-3'，反向5'-ATG CCAGAGAAGGAGATTCTGTTGATG-3'；Ifi202b，正向5'-ACTGACCAGTGTCAACTGTGAGATTG-3'，反向5'-TGACCTCCATGTAACTGTACCTCGTAG-3'；Ctss，正向5'-GCCAGCCAヰCCTCCヰCヰCヰC-3'，反向5'-TCTTGCCATCCGAATGTATCCTTGATC-3'；Itgb2，正向5'-GCAGCAGAAGGACGGAAGGAA C-3'，反向5'-ATGACCAGGAGGAGGACACCAATC-3'；Pik3ap1，正向5'-GCAGCCAACCCAGTACAGTT-3'，反向5'-GGGACAAATCCAGCCATAGA-3'；Arhgap21，正向5'-GAGAGCAGACAACGCCATCAGAAG-3'，反向5'-TCCGCCTCTTGGTCTCACTGTAC-3'；Myocd，正向5'-CTGCCAGACACCヰCACCGATG-3'，反向5'-CCACCAGCATCTTGTCCTTCTCAG-3'；MT-ND5，正向5'-CCTGGCAGACGAACAAGACATC C-3'，反向5'-CGAGGCヰCCGAヰACTAGGCATG-3'。

1.9 Q-PCR检测基因表达量 为验证RNA-Seq结果的可靠性，随机选择8个差异表达基因进行Q-PCR检测，分别从R193H和WT小鼠中提取总RNA，并利用反转录试剂盒(TaKaRa，日本)将其反转录成cDNA，利用SYBR Green试剂盒(Qiagen，德国)进行Real-time PCR，以β-actin为内参基因，利用实时荧光定量PCR仪(Bio-Rad，美国)测量吸光度，mRNA相对表达量分析采用2–ΔΔCt法。实验数据采用GraphPad Prism 6软件进行独立样本t检验并制图。P＜0.05为差异有统计学意义。

2 结 果

2.1 差异表达基因筛选 对R193H和WT两组共6个样本(生物学重复3个)进行分析，得出差异表达的基因共52个(表1)，其中28个基因上调，24个基因下调。对差异表达基因和样本同时进行分层聚类分析(图1)，可见表达模式相似的基因具有功能相关性。根据差异表达的基因构建火山图，对样本之间和基因之间的差异进行直观展示(图2)。

表1 RNA-Seq分析野生型和突变型R193H小鼠差异表达基因Tab.1 Analysis of the differentially expressed genes of wild type and mutant R193H mice by RNA-Seq

图1 R193H与WT基因差异表达模式聚类图Fig.1 Clustering diagram of R193H vs. WT gene expression pattern

图2 R193H与WT基因差异表达火山图Fig.2 R193H vs. WT genetic differentially expressed volcanic maps

2.2 利用生物信息学方法进行的基因功能分析从BP、CC、MF这3个层面进行GO注释(图3)，基于通路的分析进一步了解基因的生物学功能。将差异基因进行GO显著性富集(图4)，差异表达基因在生物过程分组中主要富集在翻译延伸(translational elongation)、固有免疫应答(innate immune response)、ATP合成耦合电子传输(ATP synthesis coupled electron transport)这三个生物学过程中，涉及转录翻译、免疫以及能量代谢相关基因。层次树状图显示出显著性功能的内在从属关系(图5)。

图3 GO功能分类Fig.3 GO function classification

图4 KEGG富集散点图Fig.4 KEGG enriched scatter plot

图5 显著性功能间的层次树形关系图Fig.5 The hierarchical tree diagram for significant functions

2.3 显著性差异基因的KEGG分析 根据差异表达的KEGG功能注释结果进行通路显著性富集分析(图6)，共得到16条显著富集的通路，涉及了氨酰基-tRNA生物合成通路(aminoacyl-tRNA biosynthesis)、NF-κB信号通路(NF-kappa B signaling pathway)、Toll样受体信号通路(Toll-like receptor signaling pathway)等。对显著性差异的通路条目进行调控网络构建，可以直观显示显著性功能之间的关系(图7)。

图6 显著通路富集度柱状图Fig.6 Significant pathway enrichment histogram

图7 显著性功能间的网络构建关系图Fig.7 Network construction diagram of significant functions

2.4 差异表达基因验证 随机选取8个差异表达的基因，每组取4个样本进行实时荧光定量PCR验证，分别为Tifab、Ifi202b、Ctss、Itgb2、Pik3ap1、Arhgap21、Myocd、MT-ND5。结果显示，突变组Tifab、Ifi202b、Ctss、Itgb2、Pik3ap1基因mRNA表达量与WT组比较明显升高，Arhgap21、Myocd、MT-ND5基因mRNA表达量与WT组比较明显降低(P＜0.05)，该结果与测序结果一致(图8)。

图8 野生型和突变型相关基因的mRNA水平Fig.8 mRNA levels of wild-type and mutant

3 讨 论

目前研究认为，RCM可能是特发性的、家族性的或继发于其他条件发生的疾病[12]。最新分子遗传技术的进展明显加深了对该疾病的理解[5]，在编码肌节蛋白的相关基因中已经鉴定了几个RCM特异性突变，其中cTnI作为重要的心肌结构蛋白，其突变可以导致RCM和家族性肥厚性心肌病[13]。本课题组拥有来源于人cTnI基因192位点突变的RCM小鼠动物模型(在小鼠系193位点，cTnI R193H)，模型小鼠与携带该突变的RCM患者有相似的病理改变、症状和表型[10,14]，从而验证了cTnI突变与RCM的关系。目前，RCM相关机制的研究颇多。有研究表明，cTnI突变可能引起心肌细胞钙敏性改变而导致RCM[15]，也有学者认为cTnI突变可能导致与其相互作用的基因或蛋白发生改变，从而导致心脏功能的改变[16]，但cTnI突变引起RCM的机制尚不明确。

转录组是特定组织或细胞在某一发育阶段或功能状态下转录出来的所有RNA的集合。RNA-Seq利用高通量测序技术对组织或细胞中所有RNA反转录而成的cDNA文库进行测序，通过统计相关reads数计算出不同RNA的表达量[17]，寻找出两组样本中差异表达的基因[18]。这种方法改变了转录序列的研究方式，主要应用于转录组学的研究。本课题组利用RNA-Seq技术探究cTnI在RCM发病机制中的生物学作用[19]。通过检测cTnI突变导致的RCM模型鼠和同月龄WT小鼠基因表达谱的差异，找出52个差异表达基因，随机挑选出8个差异表达基因利用Q-PCR进行验证，其结果初步验证了测序结果的准确性。对差异表达的基因进行GO分析发现，上调和下调的基因在细胞组分、分子功能、生物过程三个方面差异有统计学意义。其中上调的基因及产物主要分布于细胞表面和质膜外侧面，分子水平主要是具有结合活性，多参与到固有免疫应答的过程，涉及的基因有Cxcr3、Fcgr1、Icam1、Clec12a、Ctss、Lbp等。而下调的基因及产物主要分布于细胞骨架和线粒体中，分子功能主要是催化活性，参与的生物学过程主要是翻译延长和ATP合成耦联电子运输，涉及的基因有D2hgdh、ND5、Myom1、Ubr4、Spta1、Lrpprc、Arhgap21。这一结果表明，cTnI突变可能伴随免疫相关分子的上调和能量代谢的抑制。通过通路显著性富集分析，共得到16条显著富集的通路，涉及氨酰基-tRNA生物合成、NF-κB信号通路、Toll样受体信号通路。但是由于差异基因的数量并不大，通路分析的结果无统计学意义。因此，本研究后续重点分析了差异基因参与的生物学进程及参与的信号通路。通过分析该测序结果发现，差异表达的基因主要参与免疫应答、线粒体能量代谢、心肌细胞分化等过程，如Ifi202b、ND5、Itgb2、Myocd等基因。既往研究表明，cTnI某些序列是免疫调节序列，能通过免疫炎症反应导致扩张性心肌疾病[20]，而本研究结果显示，在cTnI突变所致的RCM动物模型中参与免疫应答的一系列基因表达量发生改变，提示免疫功能的改变可能在RCM的发病中具有重要作用。另有研究表明，通过转运蛋白和结合蛋白之间复杂的相互作用，调节Ca2+信号通路来调节线粒体的代谢活性，从而导致能量代谢受损可能是RCM发生的潜在机制[21]。本研究发现，线粒体能量代谢相关基因差异表达显著，提示在RCM发病中cTnI突变可能伴随着能量代谢改变。NADH泛醌氧化还原酶亚基5(ND5或MT-ND5，EC 1.6.5.3)由位于线粒体DNA的鸟嘌呤富集重链中的最大基因编码(mtDNA坐标：12337-14148)，参与的相关途径包括呼吸电子传递，通过化学渗透耦合进行ATP合成，以及通过解耦联蛋白产生热量，该基因的表达降低可能与线粒体能量代谢受损相关[22]，但其具体内在机制尚需深入研究。差异表达的基因中也发现有PI3KAP1基因，结合课题组前期研究证实cTnI与多种PI3K-Akt信号通路中的蛋白相结合[23]，表明PI3K-Akt信号通路可能通过cTnI参与RCM的发生。

综上所述，利用测序技术和特有的RCM模型，本研究初步筛选出了可能参与RCM发病机制的基因及其相关通路，可以通过寻找对差异表达基因有调控作用的药物，提高对该疾病的治疗效果。cTnI作为一种重要的心肌蛋白，作用机制十分复杂，其功能仍需进一步深入研究。