野生盘羊×巴什拜羊杂交二代尾部脂肪差异表达lncRNA的鉴定与分析

2022-05-31苏晓慧何海迎吐尔逊阿依牙力坤马海玉刘玲玲刘武军

中国畜牧兽医 2022年5期

苏晓慧，何海迎，吐尔逊阿依·牙力坤，王琼,2，马海玉，刘玲玲，刘武军

(1.新疆农业大学动物科学学院，乌鲁木齐 830052；2.吐鲁番职业技术学院，吐鲁番 838000)

目前，市场上大部分脂尾型绵羊都是由瘦尾型绵羊经过长期选择进化而来[1-2]，但是，由于当代人的饮食和社会结构的变化，绵羊尾部脂肪被视为一种能量浪费，现阶段小尾甚至无尾成为当下绵羊尾部育种的目标性状[3]，因此，了解绵羊尾部脂肪沉积的分子机制对绵羊育种至关重要。近年来常通过研究长链非编码RNA(long non-coding RNA，lncRNA)的作用来探讨动物脂肪沉积的过程，lncRNA是一类长度>200 nt的非编码RNA分子，在多种层面上调控基因的表达水平[4-5]。Huang等[6]研究发现，莱芜猪和长白猪背最长肌内脂肪中存在55个差异lncRNAs，其中XLOC_046142、XLOC_004398和XLOC_015408可通过靶向MAP活化蛋白激酶2(MAPK2)等调节猪肌内脂肪沉积；李富原等[7]研究了不同喂养期朗德鹅肝脏、腹部脂肪和胸肌中lncRNA的差异表达，结果发现，LOC106047490在3种组织中的表达量均显著下降，推测可能是高糖、胰岛素和不饱和脂肪酸诱导的，表明LOC106047490与脂代谢、糖代谢有关；Zhang等[8]对鸡腹前脂肪细胞lncRNA研究发现，在前脂肪细胞中有1 336个lncRNAs在不同分化阶段差异表达，且lncRNA是通过顺式作用于邻近的编码基因控制腹部前脂肪细胞的分化，说明lncRNA是参与调控动物脂肪沉积过程的重要因子，但报道主要集中在猪、鸡等家养动物。Ma等[9]通过对兰州肥尾羊、小尾寒羊和小尾藏羊进行高通量RNA测序发现，在尾部脂肪沉积过程中发挥重要作用的核心lncRNAs有：TCONS_00372767、TCONS_00171926、TCONS_00054953和TCONS_00373007；王鑫悦等[10]对野生盘羊与巴什拜羊杂交二代群体(F2代羊)的脂尾形状进行了K-means聚类分析，结果分成4个类型：小尾型、中尾型、大尾型和肥尾型，是理想的尾部脂肪变异和进化的研究模型。鉴于此，本试验利用全转录组测序技术对新疆2个不同脂尾型绵羊群体巴什拜羊(肥尾型)和F2代羊(小尾型)的尾部脂肪进行测序，筛选差异表达lncRNA，对其进行生物信息学分析并表达验证，以期为解析绵羊脂肪沉积调控机制提供一定的理论依据。

1 材料与方法

1.1 材料

试验动物来自于新疆塔城裕民县，随机挑选相同饲养条件下6月龄巴什拜羊及野生盘羊×巴什拜羊的杂交二代羊各6只，屠宰后立即采集尾部脂肪组织，并将新鲜样品放入2 mL无菌无酶的低温保存管中暂存于液氮中，样品带回实验室后-80 ℃保存备用。杂交二代是以野生盘羊为父本、巴什拜羊为母本进行杂交，获得F1代，再以巴什拜羊为父本、F1代为母本进行杂交，获得F2代。

1.2 方法

1.2.1 总RNA提取及检测采用Trizol法提取尾部脂肪组织总RNA，利用1.0%琼脂糖凝胶电泳检测RNA完整性，使用NanoDrop分光光度仪检测RNA浓度，巴什拜羊样品编号为CG1～CG6；F2代羊样品编号为TG1～TG6。

1.2.2 文库构建和测序总RNA去除rRNA后纯化回收，被随机打断成短片段，以片段化去除rRNA的RNA为模板，合成第一链cDNA，形成双链结构，利用T4 DNA聚合酶和Klenow DNA聚合酶将DNA的黏性末端修复为平末端，3′-末端加碱基A并加接头，经过PCR扩增获得最终测序文库，文库质检合格后采用Illumina Hiseq4000进行测序。

1.2.3 转录组测序数据分析筛选得到原始数据，与NCBI中参考基因组进行序列比对，获取比对数据，使用FastQC软件对序列质量(QC)进行评估，过滤含有适配器、高N含量或大量低质量碱基的reads，得到clean reads。利用Hisat2工具将reads映射到参考基因组[11]，使用StringTie[12]对Hisat2比对结果进行拼接，并通过ASprofile[13]软件获得备选拼接的类型和每个样本的相应表达量。转录水平的计算采用FPKM(Fragments Per Kilobase of transcript per Million fragments mapped)片段法，采用DESeq软件[14]分析样本组间差异表达，将log2|FoldChange|≥2且FDR<0.05作为筛选标准。分析巴什拜羊和F2代羊尾部脂肪差异表达lncRNA，以及这些差异表达lncRNA功能的注释分析。

1.2.4 lncRNA靶基因预测利用顺式(cis)或者反式(trans)方法对差异表达的lncRNA进行靶基因预测。选取与lncRNA距离10 kb的上游或下游基因作为顺式靶基因，采用Miranda方法对反式靶基因进行预测。

1.2.5 GO功能和KEGG通路富集分析 GO数据库用于预测和阐明基因产物在分子功能、生物过程及细胞组分方面的功能[15]。这些基因首先被映射到数据库中的GO项，然后以校正后的P<0.05为阈值，计算每个GO项的基因数，以确定显著富集的GO项。利用KEGG进行通路富集分析[16]，确定基因参与的主要生化和信号通路，以校正后的P<0.05为阈值，确定显著富集的KEGG通路。

1.2.6 引物设计挑选显著差异的7个与脂肪沉积相关的lncRNAs，根据测序所得序列设计引物；从NCBI数据库中查找β-actin内参基因序列(登录号：HM067830.1)设计引物，引物序列见表1。

表1 引物信息

1.2.7 实时荧光定量PCR 使用浓度和纯度合格的总RNA，根据PrimeScript RT Reagent Kit with gDNA Eraser说明书反转录为cDNA。以cDNA为模板进行目的lncRNA的PCR扩增。PCR反应体系20 μL：2×PCR MasterMix 10 μL，上、下游引物各0.4 μL，cDNA 1.5 μL，RNase-free ddH2O 7.7 μL。取4 μL PCR产物与2 μL 6×Loading Buffer混匀，用1.5%琼脂糖凝胶电泳检测，180 V 25 min，使用凝胶成像系统采集图像。

以β-actin为内参基因，采用CFX96型荧光定量PCR仪进行实时荧光定量PCR检测。PCR反应体系20 μL：TB Green Premix ExTaqⅡ 10 μL，上、下游引物各0.5 μL，cDNA 1 μL，RNase-free ddH2O 8 μL。PCR反应程序：95 ℃ 30 s；95 ℃ 5 s，54 ℃ 30 s，共39个循环。每个样品重复3次，采用2-ΔΔCt法计算目的lncRNA的相对表达量。

1.2.8 统计分析运用SPSS 25.0软件的单因素方差分析(One-way ANOVA)，组间采用独立样本t检验。结果以平均值±标准误表示，P<0.05表示差异显著；P<0.01表示差异极显著。

2 结果

2.1 总RNA检测及测序数据统计

提取的总RNA具有较好的完整性，28S和18S条带明显，5S条带暗淡，未发生降解，可用后续试验。

对巴什拜羊和F2代羊的尾部脂肪分别进行RNA-Seq高通量测序后，共获得1 438 500 450条原始序列数据(raw reads)；对原始数据进行过滤后，共获得1 338 966 230条高质量序列数据(clean reads)；Q20在所有样品中均是99%以上，说明单个碱基的错误率非常低；GC含量在47%～49%之间(表2)，更加验证了测序结果的可靠性。

表2 测序数据统计表

2.2 序列比对

由表3可知，各样品的原始序列数据与参考基因组的比对率在85.28%～91.06%之间，单一映读取率在64.11%～71.57%之间，多个映读取率在19.12%～21.43%之间，说明巴什拜羊和F2代羊相似性较高，测序结果质量可靠，可用于后续分析。

表3 基因组比对情况

2.3 差异表达lncRNA的筛选

使用DESeq2软件筛选巴什拜羊和F2代羊中尾部脂肪的差异表达lncRNA，并绘制火山图，结果见图1。由图1可知，共筛选到728个差异lncRNAs，与巴什拜羊相比，F2代羊中270个lncRNAs表达下调，458个lncRNAs表达上调。

2.4 lncRNA靶基因预测

lncRNA可调控临近mRNA表达，利用lncRNA与蛋白编码基因的共定位预测lncRNA的靶基因，以便进一步对获得的lncRNA进行功能研究，其中上、下游100 kb作为共定位的阈值，结果显示，728个差异表达的lncRNAs靶向17 850个靶基因，其中MSTRG.37980、MSTRG.38164、MSTRG.36912、MSTRG.36913、MSTRG.8169、MSTRG.24995和MSTRG.31389可以靶向与脂肪沉积和脂肪代谢相关的基因，如SCD、GPAM、THRSP和FASN等。lncRNA靶向与脂肪沉积和脂肪代谢相关基因的部分结果见表4。

Y轴>1.3且X轴≥1区域代表上调；Y轴>1.3且X轴≤-1区域代表下调The area of Y>1.3 and X≥1 means up-regulated；The area of Y>1.3 and X ≤-1 means down-regulated图1 差异表达lncRNA火山图Fig.1 Volcano diagram of differentially expressed lncRNA

表4 与脂肪沉积相关的部分lncRNAs及其靶基因

2.5 GO功能和KEGG通路富集分析

通过GO功能分析，差异显著lncRNA的靶基因参与到634个功能分类，注释到生物过程、细胞组分和分子功能三大类，主要富集于T细胞受体信号通路(T cell receptor signaling pathway)、T细胞受体复合体(T cell receptor complex)、SH3/SH2适配器的活动(SH3/SH2 adaptor activity)、正向调节T细胞增殖(positive regulation of T cell proliferation)、GTPase活性的正调控(positive regulation of GTPase activity)、阴性胸腺T细胞选择(negative thymic T cell selection)、肽-酪氨酸磷酸化的负调控(negative regulation of peptidyl-tyrosine phosphorylation)、MHCⅡ类蛋白复合物(MHC class Ⅱ protein complex)、膜的积分分量(integral component of membrane)、先天性免疫应答(innate immune response)、炎症应答(inflammatory response)、免疫突触(immunological synapse)、免疫反应(immune response)、质膜的外侧(external side of plasma membrane)、对病毒的防御反应(defense response to virus)、染色质沉默(chromatin silencing)、趋化性(chemotaxis)、chemokine-mediated信号通路(chemokine-mediated signaling pathway)、细胞表面(cell surface)、B细胞受体信号通路(B cell receptor signaling pathway)。其中显著富集到膜的积分分量、GTPase活性的正调控、炎症应答、免疫反应、质膜外侧、对病毒的防御反应等，最显著富集到了膜的积分分量(图2)。

KEGG通路富集结果显示，巴什拜羊和F2代羊尾脂组织差异表达lncRNA的靶基因共富集到76条通路，主要涉及病毒致癌作用(viral carcinogenesis)、Ⅰ型糖尿病(type Ⅰ diabetes mellitus)、T细胞受体信号通路(T cell receptor signaling pathway)、全身性红斑狼疮(systemic lupus erythematosus)、金黄色葡萄球菌感染(Staphylococcusaureusinfection)、原发性免疫缺陷(primary immunodeficiency)、破骨细胞分化(osteoclast differentiation)、自然杀伤细胞介导的细胞毒性(natural killer cell mediated cytotoxicty)、麻疹(measles)、利什曼病(Leishmaniasis)、IgA生产的肠道免疫网络(intestinal immune network for IgA production)、造血玻璃纸血统(hematopoietic cel lineage)、移植物抗宿主病(graft-versus-host disease)、Fcγ-R介导的吞噬作用(Fc gamma R-mediated phagocytosis)、细胞因子受体相互作用(cytokine-cytokine receptor interaction)、趋化因子信号通路(chemokine signaling pathway)、B细胞受体信号通路(B cell receptor signaling pathway)、自身免疫性甲状腺疾病(autoimmune thyroid disease)、同种异体移植物排斥(allograft rejection)、醇中毒(alcoholism)，主要通路见图3。

2.6 差异表达lncRNA的实时荧光定量PCR验证

为了验证RNA-Seq数据的准确性，选择巴什拜羊和F2代羊尾部脂肪之间差异表达的7个lncRNAs进行了实时荧光定量PCR验证，结果见图4。由图4可知，MSTRG.37980、MSTRG.38164、MSTRG.36912、MSTRG.36913、MSTRG.24995、MSTRG.8169和MSTRG.31389的实时荧光定量PCR与RNA-Seq测序结果趋势相同，说明RNA-Seq结果可靠。实时荧光定量PCR检测结果显示，7个差异表达lncRNAs中有5个(MSTRG.37980、MSTRG.38164、MSTRG.36912、MSTRG.36913和MSTRG.24995)在F2代羊中表达下调，其中MSTRG.37980、MSTRG.38164和MSTRG.36912表达量极显著低于巴什拜羊(P<0.01)，MSTRG.36913和MSTRG.24995表达量显著低于巴什拜羊(P<0.05)；MSTRG.8169和MSTRG.31389在F2代羊中表达上调，且显著高于巴什拜羊(P<0.05)。

图2 差异表达lncRNA靶基因的GO功能富集散点图Fig.2 Rich distribution map of GO function with differentially expressed lncRNA target genes

图3 差异表达lncRNA靶基因的KEGG通路富集散点图Fig.3 Rich hub diagram of KEGG pathway with differentially expressed lncRNA target genes

3 讨论

脂肪不仅是一种重要的储能物质，而且在调控动物体能量平衡过程中发挥着重要作用，并与动物产肉量、肉品质、肉风味等产肉性状密切相关。随着人们生活水平的不断提高，饮食结构发生改变，对食物质量的要求也越来越高，也逐渐发现食用过多的脂肪对人体健康有害，对现代消费者而言尾部脂肪并不理想，并降低了绵羊肉的价值，因此小尾型绵羊备受青睐[17-18]。但是，脂肪沉积是一个复杂的过程，涉及脂肪细胞的增殖与分化、发生与凋亡、甘油三酯的合成与水解等生化过程，受到多种分子和途径的调控。此外，各种生物大分子之间也存在相互作用，如RNA-RNA、DNA-蛋白质、RNA-蛋白质、蛋白质-蛋白质等，因此关于脂肪沉积相关机制的探索十分重要。

目前已对动物脂肪沉积进行过转录组测序研究，Miao等[19]研究了lncRNA在金华猪和长白猪肌内脂肪组织中的表达水平，发现119个差异表达lncRNAs，其中有6个共表达lncRNAs参与脂肪沉积和脂质代谢通路；郑竹清等[20]通过CPC(coding potential calculator)和CPAT(coding potential assessment tool)预测发现1 472个lncRNAs在山羊脂肪细胞中表达，其中29个在内蒙古绒山羊成熟前后肌内脂肪细胞中显著表达，推测这29个lncRNAs影响细胞的分泌和生长调节，从而在脂肪细胞的发育中发挥重要作用；晁天乐[21]以杜泊羊和小尾寒羊为材料，检测到47个lncRNAs参与调节十二指肠和大肠的转录后表达，发现这些lncRNAs可能在碳水化合物代谢和脂质代谢中发挥主要的作用；Ma等[9]对兰州肥尾羊、小尾寒羊和藏羊的尾部脂肪进行高通量测序，在差异表达的6个lncRNAs中发现4个lncRNAs(TCONS-00372767、TCONS-00171926、TCONS-00054953和TCONS-00373007)可能作为核心lncRNA在尾部脂肪沉积过程中发挥关键作用。本研究通过转录组测序共筛选得到728个差异lncRNAs，共鉴定到17 850个靶基因，并筛选到至少7个与绵羊尾脂沉积相关的lncRNAs，说明lncRNA在动物脂肪沉积中发挥着至关重要的调控作用。本研究采集的是不同尾型的巴什拜羊和F2代羊的尾部脂肪组织，这2种羊与前人研究所用的试验材料有较近的亲缘关系，相对降低了种间差异性，保证了结果的准确性和可靠性。

本研究鉴定到的靶基因可富集到76条信号通路，多条与脂代谢相关(如趋化因子信号通路)，且一些差异靶基因对绵羊尾部脂肪沉积的调控作用得到了部分验证，如，MSTRG.37980的靶基因SCD编码硬脂酰辅酶A去饱和酶，此酶在PPAR信号通路中发挥关键作用，以激活脂肪细胞的成脂分化[22-23]，且在动物肝脏和脂肪组织中高度表达[24]；MSTRG.38164的靶基因GPAM可以催化磷脂生物合成过程中甘油三酯合成的第一反应，调节磷脂和甘油三酯的水平[25-26]，且GPAM基因在所有脂肪组织中均有表达，在肝脏和脂肪组织中表达最高[27]。本试验结果为从lncRNA的角度分析绵羊尾部脂肪在转录组中的差异及绵羊脂肪沉积的进一步研究提供理论依据。