绵羊组织特异性长链非编码RNA的鉴定与功能分析

赵 琛，刘晓懿，王喜宏

( 西北农林科技大学 动物科技学院，陕西 杨凌 712100)

长链非编码RNA(Long non-coding RNA, lncRNA)是指长度超过200个核苷酸，且不具备蛋白编码能力的RNA[1]。lncRNA在染色质调控、转录调控和转录后调控中具有重要作用[2]，主要通过招募染色质重塑复合物，介导染色体间的相互作用，作为转录因子“向导”或“诱饵”，充当miRNA的“海绵”等方式发挥抑制或激活作用[3]。近年来高通量测序方法的迅速发展使得对家畜lncRNA的研究增加，Bakhtiarizadeh等[4]分析绵羊8种组织的RNA-seq数据并报告了第一个基因间lncRNA(lincRNA)数据集(n=325)。Bush等[5]根据绵羊3.1版本参考基因组(Oar_v3.1)在11个绵羊组织中鉴定出11646个新的lncRNA。但是与其他哺乳动物(例如人，小鼠，牛和猪)相比，绵羊基因组中注释的lncRNA很少。NONCODE数据库是一个比较全面的ncRNA相关注释的数据库，包含人类，小鼠，大鼠，牛，鸡和猪等17个物种的ncRNA信息[6]，但是该数据库中缺少绵羊的lncRNA数据集。本研究以绵羊为研究对象，以高通量测序为手段，筛选并鉴定新型的lncRNA及组织性lncRNA，旨在发现对家畜重要经济性状有调控作用的lncRNA，丰富家畜lncRNA数据库，为lncRNA在育种工作中的应用提供借鉴和参考。

1.1 试验数据

绵羊(Ovisaries)RNA-seq测序数据包括特克赛尔羊♂×哈萨克羊♀杂交绵羊F1代个体、哈萨克羊、特克塞尔羊♂×苏格兰黑脸羊♀的杂交F1代个体的192个组织样本。由肌肉(半腱肌、半膜肌、背最长肌)、皮肤(腹侧皮肤、体侧皮肤)、脂肪(尾脂、皮下脂肪)、心脏、肺、肝脏、肾脏、小肠、脾脏构成，包含新生、出生后1周、8周、3个月、5个月、6个月与成年的数据。原始数据下载自NCBI SRA数据库，BioProject号为PRJNA485657[7]。

1.2 序列比对与转录本组装

使用Trimmomatic(版本0.36)[8]对每个样本的转录组测序原始数据进行质控，具体参数如下：LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:40。去除接头序列和过滤低质量的读段后，通过STAR软件(版本2.5.1)[9]将样本读段比对到绵羊参考基因组(Oar_v4.0)[10]。利用StringTie软件(版本1.3.4)[11]将比对结果进行样本的组装，并通过String Tie --merge将多个样本组装合并，得到绵羊的转录本集合。

1.3 lncRNA的鉴定

lncRNA的鉴定步骤：(1)去除长度小于200 bp，外显子数目小于2的转录本[12]；
(2)筛去蛋白编码基因；
(3)利用CPC2软件[13]，CNCI软件[14]预测候选lncRNA的编码能力，去除评分不合格的转录本；
(4)使用blastx与UniProt数据库(版本2016_04)[15]和Pfam数据库(版本29)[16]中的已知蛋白质进行比对，认为E值≤10-3，比对长度≥10个氨基酸和一致性≥95%的转录本为“潜在的编码基因”，将此类转录本去除[17]。保留筛选下来的转录本作为绵羊候选lncRNA转录本。统计lncRNA与编码基因的外显子数目、GC含量、表达量等特征。

1.4 表达量评估

采用StringTie软件(版本1.3.4)[11]计算lncRNA在每一个组织中的表达量，使用每千碱基外显子百万片段数(fragments per kilobase of exons per million fragments mapped，FPKM)来评估表达水平，并保留至少在任意一个组织中有20%以上的样本FPKM表达量>1的转录本[7]，作为绵羊lncRNA的组织表达谱。

1.5 组织特异性分析

选取每种组织的表达量的平均值代表该组织，对候选lncRNA转录本计算τ组织特异性指数。τ指数是基因特异或广泛表达的指标，取值范围为0-1，其中1表示仅在一个组织中特异表达，而0表示所有组织的广泛表达[18]。给定基因的τ指数可以用以下公式计算:

式中n是被研究组织的数量，xi是给定组织的表达谱。基于τ计算组织特异性分数，认为τ≥0.8且在所有组织中表达量排名前两位的为具有组织特异性的lncRNA，分别计算肌肉、脂肪、皮肤特异性的lncRNA。

1.6 靶基因功能预测

lncRNA的顺式作用主要通过影响邻近蛋白编码基因的表达发挥调控功能[19]。lncRNA的反式作用主要通过识别与候选lncRNA共表达的mRNA，调控远距离mRNA的转录激活与表达[20]。(1)顺式作用预测：皮尔森相关系数是用来反映两个变量线性相关程度的统计量。提取位于lncRNA上游10 kb，下游100 kb内的编码基因作为靶基因[21]。(2)反式作用预测：采用皮尔森相关系数法分析样本间lncRNA与蛋白编码基因的相关性，取相关性较高的蛋白编码基因作为靶基因(r≥0.90或r≤-0.90, Pvalue<0.01)[22]。

1.7 功能富集分析

使用R软件中的ClusterProfiler软件包[23]对基因进行KEGG通路富集分析，设置adjust.p<0.05作为显著性阈值。使用Metascape[24]对候选靶基因进行基因本体论(gene ontology, GO)富集分析，设置Pvalue<0.05作为显著性阈值。

1.8 lncRNA-mRNA互作网络的构建

利用STRING数据库(版本11.0)[22]分别对肌肉特异性lncRNA、皮肤特异性lncRNA、脂肪特异性lncRNA反式作用的靶基因数据集进行蛋白质相互作用(Protein-Protein Interaction, PPI)网络预测。整合PPI网络与共表达的lncRNA-mRNA来构建最终的互作关系。利用Cytoscape软件(版本3.7)[25]中的ClusterONE插件(版本1.0)[26]检测网络互作中的显著模块。设置Pvalue≤0.05，表达模块中最少的基因数量为3。使用Cytoscape软件[25]对lncRNA-mRNA互作网络进行可视化。

1.9 数量性状位点(QTL)分析

绵羊QTL数据(Oar_v4.0)来自AnimalQTLdb[27]。通过bedtools intersect[28]比较QTL和lncRNA在全基因组中的位置，预测lncRNA的功能。起始位置和结束位置都在QTL区域内的lncRNA被认为成功注释在QTL中[29]。

2.1 lncRNA的预测结果与组织特异性分析

利用软件CNCI、CPC2与蛋白质数据库Pfam和UniProt分别预测出94 721条、75 580条、58 136条和94 987条lncRNA，四者的交集为36 340条lncRNA，认为其为高可信度的绵羊lncRNA(图1A)。对lncRNA和蛋白编码基因的基本特征进行比较分析，lncRNA的GC含量为47.17%，低于编码基因的GC含量48.79%(图1B)；
lncRNA外显子平均数目为3.45个，低于蛋白编码基因的12.8个；
但外显子个数在2～10区间中lncRNA为35 423个，多于编码基因的23 466个(图1C)；
lncRNA和蛋白编码基因的平均表达值分别为0.38和4.49，lncRNA的平均表达量低于蛋白编码基因(图1D)。

图1 绵羊lncRNA与蛋白编码基因的基本特征比较

基于RNA-Seq数据的主成分分析结果表明(图2)，相似类型的组织样本更倾向于聚集在一起，如半腱肌、半膜肌、背最长肌三种肌肉组织聚集，体侧皮肤、腹侧皮肤两种皮肤组织聚集，尾脂和皮下脂肪两种脂肪组织聚集，呈现出明显的组织特异性。试验共获得肌肉特异性lncRNA 213个，皮肤特异性lncRNA 467个，脂肪特异性lncRNA 295个。

图2 组织样本的主成分分析

2.2 肌肉特异性lncRNA

基于绵羊的213个肌肉特异性lncRNA预测到了255个上下游基因。为了评估上下游基因的调控能力，从AnimalTFDB3.0数据库[30]中提取了绵羊转录因子集合，将其与顺式作用靶基因进行比较，发现有14个靶基因为绵羊转录因子。对肌肉特异性lncRNA上下游靶基因进行GO分析，发现显著富集在肌肉系统过程(19，GO:0003012)(P<0.05)和肌肉结构发育(21，GO:0007517)(P<0.05)。Six1为绵羊转录因子，是MSTRG.71513.3(7:69708802-69711134)的邻近基因，对骨骼肌发育和肌纤维类型转化具有调控作用[31]。肌肉特异性lncRNA共预测出63个反式作用的靶基因，GO富集注释结果显示，可注释到7个功能条目，显著富集是肌肉系统发育(16，GO:0003012)(P<0.05)，肌肉器官发育(11，GO:0007517)(P<0.05)，钙介导信号(7，GO:0019722)(P<0.05)，肌动蛋白丝运动的调节(3，GO:1903115)(P<0.05)。括号内包含注释到该条目的基因数目与GO编号(图3A)。

图3 肌肉与皮肤特异性lncRNA反式作用靶基因的GO富集分析

2.3 皮肤特异性lncRNA

基于绵羊的467个皮肤特异性lncRNA预测到了1006个上下游基因，有44个为绵羊转录因子。对靶基因进行GO分析，共富集到11个功能条目，其中最显著富集的为皮肤发育(28，GO:0043588)(P<0.05)，包括表皮生长因子家族基因AREG，骨形态发生蛋白受体1A型(BMPR1A)，角蛋白基因(KRT)，细胞间粘附分子(ICM)等。AREG位于MSTRG.68599.1(6:88881341-88903389)邻近，KRT(KRT9，KRT15，KRT17)位于MSTRG.12764.8(11:41165824-41274496)邻近，ICM(ICAM1，ICAM4，ICAM5)位于MSTRG.62954.1(5:12366881-12385388)邻近。皮肤特异性lncRNA共预测出334个反式作用的靶基因。GO富集注释结果显示(图3B)，lncRNA反式作用的靶基因可注释到32个功能条目，显著富集的前7位分别是角质化(28，GO:0070268)(P<0.05)，蜕皮周期(14，GO:0042303)(P<0.05)，牙齿发生(11，GO:0042476)(P<0.05)，表皮发育调节(10，GO:0045682)(P<0.05)，长链脂肪酰基辅酶A代谢(7，GO:0035336)(P<0.05)，蛋白交联(7，GO:0018149)(P<0.05)，皮肤屏障的建立(6，GO:0061436)(P<0.05)。进一步选取更显著的lncRNA-mRNA对(r≥0.95或r≤-0.95)进行共表达分析，发现可分为4个表达模块，黄色圆圈表示ClusterONE[26]识别的显著表达模块，红色圆圈表示组织特异性lncRNA，蓝色圆圈表示lncRNA反式作用靶基因。若仅关注包含lncRNA的表达模块，最显著的为MSTRG.39053.2(2:228846018-228866160)所在模块(P<0.05)(图4)，共10个基因，包含TGM1、IL36RN、PLA2G4E、THEM5、CERS3、PNPLA1、ALOX12B、KRT28、DSG3编码基因，对该模块富集分析发现这些基因显著富集在角质化(4，GO:0070268)(P<0.05)。

图4 皮肤lncRNA-mRNA共表达网络图

2.4 脂肪特异性lncRNA

基于绵羊的295个脂肪特异性lncRNA预测到了787个上下游基因，有43个为绵羊转录因子。对上下游基因进行GO分析，显著富集到脂质定位(22，GO:0010876)(P<0.05)。转录因子HOXC基因HOXC11、HOXC12、HOXC13位于LNC12062(LOC105611004，3：132242058-132246141)邻近。脂肪特异性lncRNA共预测出79个反式作用的靶基因。KEGG数据库注释结果显示，注释到8条通路(P<0.05)，注释基因数最多的前5位分别为甘油脂代谢(5)(P<0.05)、脂肪酸伸长率(3)(P<0.05)、不饱和脂肪酸的生物合成(3)(P<0.05)、缬氨酸，亮氨酸和异亮氨酸的降解(3)(P<0.05)、脂肪的消化吸收(3)(P<0.05)。GO富集注释结果显示，注释到11个功能条目，显著富集的是脂质代谢(14，R-HSA-556833)(P<0.05)，DAG和TAG的酰基链重塑(14，R-HSA-556833)(P<0.05)，跨化学突触的传递(7，R-HSA-112315)(P<0.05)，胆汁盐和有机酸，金属离子和胺化合物的转运(4，R-HSA-425366)(P<0.05)。lncRNA-mRNA共表达分析发现了1个显著表达模块(图5)，包含8个基因，其中包含2个预测的lncRNA——MSTRG.57079.1(3:166886406-167025566)、MSTRG.68891.1(6:100780317-100820136)。对该模块进行GO富集分析发现主要富集在酰基辅酶A代谢过程(4，GO:0006637)(P<0.05)，KEGG通路最显著富集于甘油三酯的生物合成(P<0.05)，这两个过程都与脂肪合成相关。括号内的数字表示注释到该条目的基因数。

图5 脂肪lncRNA-mRNA共表达网络图

2.5 QTL分析

绵羊QTL数据库[27]包含肌肉相关QTL 681个，脂肪相关QTL 231个。对组织特异性lncRNA与绵羊QTL数据库的重合关系进行检索，共发现14个lncRNA可能作为绵羊肌肉与脂肪相关性状的候选lncRNA。存在6个肌肉特异性lncRNA与肌肉相关QTL区域重合，其中有3个位于胴体肌肉重量QTL区域，1个位于肌肉密度QTL区域，2个位于体重QTL区域(表1)。11个脂肪特异性lncRNA与脂肪相关QTL区域重合，其中有5个位于尾脂沉积区域，6个位于全部脂肪区域(表1)。MSTRG.77930.1(9:93629571-93634360)与尾脂沉积QTL重合，其顺式作用基因之一为与脂质合成相关的DEPTOR，可通过抑制胰岛素信号mTORC1介导的反馈抑制作用来促脂肪合成[32]。

表1 位于QTL区域的组织特异性lncRNA

本研究共鉴定36 340个绵羊lncRNA，发现肌肉特异性lncRNA 213个，皮肤特异性lncRNA 467个，脂肪特异性lncRNA 295个，扩展了绵羊基因组lncRNA集合，为进一步分析提供了良好的基础。与编码基因相比，lncRNA的GC含量、表达丰度较低，但平均长度更长，并且组织特异性高。这些基本特征与前人报道的哺乳动物lncRNA相似[33]。

肌肉质量是决定畜禽经济价值的重要因素之一，肌肉发育和生长过程受遗传、环境和营养的影响。lncRNA已被证实通过顺式或反式作用调节肌肉生长发育[34]，很可能是重要的潜在的肌肉生长调节因子。Cai等[35]发现位于鸡Six1基因上游432 bp的lncRNA-Six1通过编码一种微肽激活Six1基因，发挥顺式调节作用，促进细胞增殖并参与肌肉生长。Mueller等[36]发现小鼠成肌转录因子MyoD基因的转录起始位点上游5 kb处编码的lncRNA MUNC，可以促进MyoD在骨骼肌生成。MSTRG.46395.2的靶基因是肌球蛋白1(MYOM1)，位于MYOM1转录下游区域。MYOM1在横纹肌的细胞分化中编码一种高度特异性的蛋白质标记，并可能在肌原纤维的组装或维持中起作用[37]。本研究发现MSTRG.46395.2在肌肉组织中特异表达，可能表明MSTRG.46395.2与MYOM之间可能存在调控关系，但是需要进一步研究以阐明确切的作用机理。Six1是参与骨骼肌发育的重要转录因子[38]，肌肉特异性lncRNA MSTRG.71513.3位于其下游7 kb。有研究表明Six1不仅可以调节滩羊的肌肉纤维类型，而且还可以调节肉的嫩度，可以作为选择肉嫩度的候选基因[39]。MSTRG.71513.3在肌肉组织中特异表达，推测MSTRG.71513.3与Six1之间可能存在调控关系，调控肌肉发育。但是确切的作用机理需要进一步研究。

毛囊发育和周期性生长过程中受多种信号通路调控，主要有WNT通路、转化生长因子(TGF)家族、丝裂原活化蛋白激酶(MAPK)家族、Notch信号通路、骨形成蛋白(BMP)信号通路等[40]。有研究发现绵羊胎儿皮肤初级毛囊中下调的lncRNA与角蛋白(KRT14和KRT15)、WNT信号通路(WNT16和SFRP1)等通路基因相互作用，进而影响表皮和羊毛发育[41]。本研究发现，皮肤特异性lncRNA的上下游基因与反式作用共表达基因主要与角质化、毛囊发育相关。Notch信号通路和毛囊的形态发生具有直接的关系，可能参与了毛囊发育的早期环节，并在毛干分化过程中起重要作用[40]。DLL4(Notch配体)位于MSTRG.70607.8(7:33364804-33478104)邻近，猜测该皮肤特异性lncRNA可能对毛囊发育具有调控作用。角蛋白是一种中间纤维形成蛋白，是复杂上皮细胞的主要成分[42]。有研究表明具有无效的KRT17基因的小鼠在出生后的第一周内无法形成完整的皮毛[43]，MSTRG.12764.8邻近基因包括多个KRT基因KRT9、KRT15、KRT17，猜测其对表皮角质化及毛囊发育有重要作用，但具体作用机制待进一步分析。

lncRNA在脂质代谢中的调节功能已经在不同生物体中得到了证实，包括人类[44]，小鼠[45]，猪[46]，牛[47]等。有研究表明，皮下脂肪组织中HOXC簇的表达上调[48]。还有研究发现，HOXC11，HOXC12，HOXC13等发育基因与绵羊尾部脂肪沉积有关[49]。本研究发现LOC105611004分别位于HOXC11，HOXC12，HOXC13附近，猜测该lncRNA可能对绵羊尾脂沉积有调控作用。还发现MSTRG.77930.1顺式调控的靶基因是DEPTOR。DEPTOR通过激活脂肪细胞中的过氧化物酶体增殖物激活受体γ自主促进脂肪生成，是脂肪形成的新调节剂[50]。lncRNA-mRNA共表达分析发现了2个lncRNA(MSTRG.57079.1和MSTRG.68891.1)出现在与脂肪相关的表达模块中，与酰基辅酶A代谢和甘油三酯的生物合成功能相关，有研究表明骨骼肌脂肪储存的增加通常伴随脂质代谢产物(例如甘油二酯和长链酰基辅酶A)的逐步积累[51]，推测这两个lncRNA可能对脂肪代谢具有调控作用。

本研究通过分析绵羊不同转录组的RNA-seq数据，获得了新的绵羊lncRNA数据集，并筛选出了具有肌肉、皮肤、脂肪组织特异性的lncRNA。基于顺式作用和反式作用两个方面，对lncRNA靶基因的功能进行分析，发现其可能在组织代谢、发育等多个功能中发挥调控作用，推测肌肉特异性lncRNA可能与肌肉生长功能相关，脂肪特异性lncRNA可能与脂肪合成相关，皮肤特异性lncRNA可能与表皮形成和毛囊发育相关。同时通过QTL筛选，获得了与绵羊肌肉、皮肤与脂肪性状相关的候选长链非编码RNA。本研究为探究lncRNA在绵羊重要经济性状的功能提供了新见解。

猜你喜欢反式绵羊特异性揭开反式脂肪酸的真面目中国生殖健康(2020年5期)2021-01-18数绵羊意林(儿童绘本)(2019年6期)2019-08-05数绵羊小学生必读(低年级版)(2019年12期)2019-04-20揭开反式脂肪酸的真面目中国生殖健康(2018年5期)2018-11-06奔跑的绵羊幼儿画刊(2018年7期)2018-07-24世卫组织建议减少饱和脂肪和反式脂肪摄入食品与机械(2018年5期)2018-07-16精确制导 特异性溶栓中国医疗保险(2017年5期)2017-05-17巧计得绵羊数学大王·低年级(2016年8期)2016-05-14全反式虾青素对环磷酰胺诱发小鼠睾丸氧化损伤的保护作用中国民族医药杂志(2016年5期)2016-05-09BOPIM-dma作为BSA Site Ⅰ特异性探针的研究及其应用三峡大学学报(自然科学版)(2016年6期)2016-04-16