冷应激后民猪背脂中差异lncRNA的筛选与分析

马骏杰，汪 亮，刘 娣，张冬杰*，杨秀芹*

(1.东北农业大学动物科学技术学院，哈尔滨 150030； 2.黑龙江省农业科学院畜牧研究所，哈尔滨 150086；3.农业农村部种养结合重点实验室，哈尔滨 150086)

低温灾害是养殖生产中经常会遇到的一种自然灾害，防范不及时会造成巨大的经济损失。已有研究表明，哺乳动物如果突然暴露在低温环境下或长时间处于低温环境中均会增加冻伤或患低体温症的风险，损害认知能力，并引发心血管死亡[1-3]。机体在应对低温伤害时，首先会启用颤栗性产热维持体温恒定，但随着时间的延长，将转变为非颤栗性产热，即代谢产热，这是一种通过提高组织代谢率来增加产热的方式。非颤栗性产热作用最强的组织是褐色脂肪组织，但在机体内含量很少，机体主要以白色脂肪为主。随着人们对脂肪研究的深入，发现白色脂肪在机体遭受冷应激后也会发生褐色化，增加产热[4]，脂肪组织是研究哺乳动物冷应激相关能量代谢的最佳模型[5]。脂肪组织受交感神经系统控制，由酪氨酸羟化酶(TH)阳性儿茶酚胺能神经元介导，这些神经元从椎旁交感神经节向脂肪组织传导[6-7]。冷暴露刺激交感神经释放儿茶酚胺，儿茶酚胺激活脂肪细胞和基质细胞中的肾上腺素能受体，引发脂解、白色脂肪褐变和适应性产热[8]。

猪与其他哺乳动物间存在一个显著的不同，即它的解偶联蛋白1(uncoupling protein, UCP1)基因，一个决定褐色脂肪细胞生成的关键基因，在进化过程中因未知原因发生了突变，使其无法产生褐色脂肪[9]。有研究认为,UCP1同家族的UCP3可在冷环境下替代UCP1发挥功能，促进猪白色脂肪代谢产热[10]。也有人将鼠的UCP1和猪的过氧化物酶体增殖物激活受体γ共激活因子1α(peroxisome proliferator-activated receptor γ coactivator-1α, PGC-1α)共同过表达发现，可将猪的白色前体脂肪细胞褐色化[11]。民猪是生活在中国北方地区的一个地方猪种，脂肪含量高，具有明显的耐寒性，在遭受同样冷刺激的情况下，表型和生理变化都优于大白猪[12]，是研究冷应激下脂肪代谢的良好素材。

长链非编码RNA(long non-coding RNA, lncRNA) 属于非编码RNA的一种，在发育和基因表达中发挥着复杂而精确的调控功能。目前，人们对lncRNA的认识还处在初级阶段，起初被认为是基因组转录的“噪音”，是RNA聚合酶II转录的副产物，不具有生物学功能。后来陆续发现，lncRNA可参与X染色体沉默、基因组印记、染色质修饰、转录激活、转录干扰以及核内运输等多种重要的调控过程。哺乳动物基因组序列中4%～9%的序列产生的转录本是lncRNA(编码蛋白RNA的比例是1%)。虽然关于lncRNA的研究进展迅猛，但是绝大部分lncRNA的功能仍不清楚。

本研究拟以遭受24 h冷应激后的民猪为试验材料，利用转录组测序技术分析背脂内lncRNA的表达变化情况，对其靶基因进行预测并开展联合分析，探讨民猪耐寒性强的分子机理。

1 材料与方法

1.1 试验动物

选择6月龄雌性民猪6头，由黑龙江省农业科学院畜牧研究所民猪养殖场提供。试验动物随机分成两组，一组置于正常舍内饲养，设为对照组，样品编号为MC1、MC2和MC3；另一组置于室外半敞篷舍内饲养，设为试验组，样品编号为ML1、ML2和ML3。对照组环境温度控制在(18±2)℃，试验组环境温度控制在(-18±5)℃，两组饲喂相同的饲粮，均为自由采食、饮水。处理24 h后，屠宰，取背部皮下脂肪，置于液氮中带回实验室，-80 ℃保存备用。

1.2 方法

1.2.1 背部脂肪总RNA的提取与检测 采用常规的Trizol法提取总RNA，每个个体提取1份，共6份样品。1%的琼脂糖凝胶电泳检测RNA样品是否有降解及杂质，凯奥K5500分光光度计检测样品纯度，安捷伦2100 RNA Nano 6000 试剂盒检测RNA样品的完整性和浓度。

1.2.2 lncRNA文库的构建及测序 每个样品取3 μg总RNA作为起始量构建lncRNA文库。使用Ribo-ZeroTMMagnetic Gold试剂盒去除rRNA，向反应体系中加入碎片化缓冲液(fragmentation buffer)使RNA成为短片段，再以片段后的RNA为模板，用六碱基随机引物合成cDNA第一链，并加入缓冲液、dNTPs、RNase H和DNA Polymerase I 合成cDNA第二链。使用QiaQuick PCR试剂盒纯化合成的cDNA，经末端修复、加碱基A，加测序接头后，再次回收目的片段。加入UNG酶消化cDNA 第二链，并进行PCR扩增和回收纯化。最后用构建好的文库进行Illumina测序。

1.2.3 民猪背脂内受冷应激诱导的lncRNA的筛选与分析 本研究中lncRNA的基本筛选条件：1)转录本长度大于等于200 bp，外显子个数大于等于2；2)计算每条转录本的reads覆盖度，筛去所有样本中均小于5的转录本；3)利用gffcompare同猪的基因组注释文件进行比较，筛除已知的mRNA及其他非编码RNA(rRNA、tRNA、snoRNA、snRNA等)；4)根据比较结果中的class_code信息(如果class_code="u",则表明是lincRNA；如果class_code="i",则表明是intronic lncRNA；如果class_code="x"，则表明是anti-sense lncRNA)筛选潜在的各类lncRNAs。筛选完成后，同时使用4个编码潜能分析软件CNCI、CPC、PFAM和CPAT判断所筛选出的lncRNA是否具有编码潜能。

采用DEseq软件包进行差异表达分析，选取|log2Ratio|≥1和q<0.05的基因作为显著差异表达筛选条件，得到上、下调基因个数。利用R软件(版本号：V3.1.1)对不同样本中差异表达的mRNA进行层次聚类分析，利用NCBI、Uniprot、GO和KEGG等数据库对差异表达基因进行注释，获得差异表达基因详细描述信息。利用Ensembl、NCBI、GENCODE和HGCN等数据库进行lncRNA的注释。

1.2.5 民猪背脂内受冷应激诱导的差异表达mRNA 的GO功能分析及KEGG通路分析 利用猪的GO注释数据库直接进行GO分析。统计GO数据库中第三层条目中差异表达基因(区分上调表达和下调表达)在该条目里的个数，并计算百分比。对KEGG中每个Pathway应用超几何检验进行富集分析，找出差异表达基因显著性富集的Pathway。

1.2.6 民猪背脂内受冷应激诱导的lncRNA与mRNA的联合分析 对lncRNA和mRNA进行结构、表达量和组织特异性的比较分析。对差异表达的lncRNA的靶基因开展cis和trans两种方式的靶标预测，通过靶基因间接预测其功能。将lncRNA 相邻位置(上、下游50 kb)的蛋白编码基因筛选出来作为以cis方式作用的靶基因，根据lncRNA 同mRNA表达量的相关系数(相关系数corr≥0.9)筛选以trans方式作用的靶基因。对差异lncRNA的靶基因进行功能富集分析，分析方法同差异mRNA的分析。根据识别出的靶基因，利用Cytoscape软件绘制lncRNA-mRNA的调控网络图。

1.2.7 民猪背脂内受冷应激诱导的差异lncRNA的qPCR验证 利用Roche实时荧光定量PCR仪LightCycler480Ⅱ对试验组显著高表达的10个lncRNAs进行SYBR Green实时定量检测，引物信息见表1。反应体系：cDNA样品0.5 μL，2×SYBR Green PCR Mixture 10 μL，特异性上、下游引物各0.5 μL，灭菌水补充至20 μL。反应程序：95 ℃ 10 min； 95 ℃ 15 s，60 ℃ 30 s，72 ℃ 30 s，40个循环。根据2-ΔΔCt法计算各模板中目的基因相对于内参基因β-actin的表达量。

表1 qPCR检测用引物序列

2 结 果

2.1 冷应激下民猪背脂的转录组测序

为了全面分析民猪在遭受冷应激后背脂中表达发生变化的lncRNA，对试验组和对照组共计6个个体的RNA开展了全转录组去rRNA测序。每个样品平均获得87.5百万个raw reads，质量控制后获得约83.0百万个clean reads。与猪的参考基因组序列(Sus scrofa 11.1)比对，平均94.4%的clean reads可以比对到基因组上，大约有75.7百万个reads可比对到唯一位置。通过转录本的富集和重建总计获得了24 305个转录本(FPKM≥0.01)，包括18 772个已知的转录本。

2.2 民猪背脂内 lncRNA的筛选与鉴定

是否具有编码潜能为判断新转录本是否为lncRNA 的关键条件，将CNCI、CPC、PFAM和CPAT软件筛选出的lncRNA的交集作为后续分析的数据集，共计7 856个(图1)。将本研究所获得的lncRNA与mRNA进行转录本长度、外显子个数、组织特异性以及表达量的对比分析。结果发现，lncRNA和mRNA的长度都集中在2 900 bp左右，整体趋势一致(图2A)。主要的lncRNA有2个外显子(最多可达31个外显子，平均是2.99个)，显著低于编码蛋白基因(最高可达118个外显子，平均为9.64个)[14](图2B)。与mRNA相比，lncRNA表现出更强烈的组织特异性(图2C)。同等条件下，mRNA的表达水平要相对高于lncRNA(图2D)。 这与以往在猪、鼠和人上的研究结果一致[15-17]。

2.3 受冷应激诱导的lncRNA和mRNA的筛选与分析

利用FPKM值定量估计基因表达值，以对照组为参照，筛选出受冷应激诱导的166个lncRNAs，其中上调表达的78个，下调表达的88个。但仅有2个下调表达的为已知lncRNAs，其余全部为新发现的lncRNAs。269个mRNAs发生了显著变化，其中上调表达的141个，下调表达的128个。显著上调表达基因中，双孔K+通道蛋白(KCNK3)变化倍数最大，为7.49倍，其次是激酶非催化C-叶结构域蛋白(KNDC1)和减数分裂特异性端粒结合蛋白(TERB1)，分别达到了5.92倍和5.44倍。显著下调表达基因中，果糖二磷酸醛缩酶B(ALDOB)变化倍数最大，为6.24倍，其次是Slit引导配体1(SLIT1)和微管相关蛋白6(MAP6)，分别下调了5.98倍和5.20倍。

聚类分析往往用于反映不同试验条件下样本差异表达RNA的变化模式，可以很直观反映出不同条件下RNA的表达量变化情况。利用R软件(v3.1.1)对试验组和对照组间差异表达的RNA进行层次聚类分析，结果见图3，表达量的变化用颜色的变化表示，蓝色表示表达量较低，黄色表示表达量较高。由聚类图可知，组内个体间表达模式基本一致，说明重复性较好，但组间差异较大，尤其是lncRNA 变化明显，说明冷应激对lncRNA的影响较大。

2.4 差异表达基因的GO注释与KEGG通路分析

GO富集分析结果表明，153个基因被富集到生物过程(BP)，7个基因被富集到分子功能(MF)，细胞组分(CC)中无显著富集的基因。在生物过程里，57个条目被显著富集，很多基因被富集到防御反应、对生物刺激的反应以及病毒的反应等条目下(图4)。在分子功能里只有“双链RNA结合”条目被富集。

KEGG通路富集分析表明，269个基因富集到156个信号通路中，但均未达到显著富集的程度。按照q值大小，排在前5位的信号通路分别是PPAR信号通路、视黄醇代谢通路、细胞外基质受体相互作用通路、果糖和甘露糖的新陈代谢通路和糖酵解/糖异生通路，这5个通路基本都与脂类和糖类代谢相关。

2.5 受冷应激诱导lncRNA的靶基因预测与功能分析

通过相关系数分析发现，208个差异表达的mRNAs与144个差异表达的lncRNAs共表达，其中193个靶基因预测为受lncRNA的反式调控(trans)，仅有15个靶基因预测为顺式调控(cis)。对这208个靶基因开展GO功能分析和KEGG通路分析，结果发现，MF分类下富集6个条目，分别是poly(A)RNA结合、RNA结合、SMAD结合、杂环化合物结合、有机环状化合物结合和核酸结合；CC分类下富集27个条目，排在前5位的是细胞内膜结合细胞器、细胞内部分、膜结合细胞器、核斑点、核质部分；BP分类下富集52个条目，排在前5位的是RNA加工、mRNA加工、RNA剪接、细胞代谢过程和初级代谢过程。KEGG通路分析中，显著富集的通路有2个，分别是真核生物的核糖体生物发生(ribosome biogenesis in eukaryotes)和基底细胞癌(basal cell carcinoma)。

2.6 受冷应激诱导的lncRNA与关键靶基因的信号通路

利用Cytoscape软件绘制的lncRNA-mRNA共表达网络表明，47个lncRNAs处于网络的中心，调控114个基因(图5)，1个lncRNA可同时调控多个靶基因，如MSTRG.90687可同时调控UCP3、TTC39A、SLC37A4、KLF11、ADAMTS18和PLIN5；多个lncRNAs也可同时调控多个靶基因，这些lncRNA和基因可能在机体应对冷应激方面起到关键作用。

2.7 qRT-PCR验证

为了验证RNA测序结果的准确度，随机选择10个差异表达的lncRNAs进行 qRT-PCR试验验证。结果显示，qPCR定量检测值与测序值相比，变化趋势基本一致(图6)，表明测序结果可信。

3 讨 论

lncRNA一般指不具有蛋白编码能力，转录本长度超过200个核苷酸的RNA分子，与细胞周期和分化、发育、生殖、性别调控、衰老以及多种人类疾病密切相关。其在白色脂肪沉积及褐色脂肪代谢产热方面均起着重要的调控作用[18-19]。研究发现，冷应激下lncRNA uc.417的表达显著下调，可解除对p38-MAPK信号通路磷酸化水平的抑制，促进线粒体产热基因的表达和氧化呼吸[20]，AK079912和GM13133的表达水平显著上调，促进白色脂肪棕色化[21]。由此可见，冷应激下，不同lncRNA行使着不同的生物学功能。

本研究以冷应激下民猪的背部脂肪为试验材料，采用转录组测序技术对lncRNAs的表达变化情况进行了分析，发现冷应激后上调表达的lncRNA数量(78个)略低于下调表达的lncRNA数量(88个)，但mRNA上调表达的数量(141个)略高于下调表达的数量(128个)，两者的变化趋势并不完全一致。对差异表达基因进一步分析发现，多个干扰素刺激基因(ISGs)，如RSAD2、ISG20、ISG15、MEGF6、DDX60、MX2、MX1和IFI6等均发生了显著上调，它们在宿主抗病毒防御过程中发挥着关键作用[22]。此外，还有与病毒感染、血液中免疫细胞的数量和分布相关的DHX58[23]，与先天免疫应答相关的PTX3[24]和抗猪瘟病毒的OASL基因[25]的转录水平也出现了显著上升。前人在急性冷暴露的大鼠白色脂肪细胞内也发现了多个与免疫相关的基因受到诱导[26]，表明冷刺激诱发了机体的天然免疫系统。GO注释结果也提示，民猪在遭受冷应激24 h后，并没有马上启动非颤栗性产热过程，而是先启动了机体的防御和免疫应答反应，抵抗冷刺激的同时，提升机体免疫力。

机体在遭受冷刺激后，会发生局部血液循环不畅，组织缺血缺氧，甚至发生坏死。本团队在前期研究中也发现，大白猪在遭受重度冷应激后，其耳组织和体表会发生红肿与坏死。本研究同样筛选到多个与呼吸、氧化应激和血液循环相关的基因。如差异倍数最大的KCNK3基因(log2fold change=7.49)，该基因功能丧失是肺动脉高压发病的关键所在[27]。CCDC38基因在支气管上皮细胞中表达，与一秒率(FEV1/FVC，第一秒用力呼气量占所有呼气量的比例)相关[28]。GAS2L2的变异会导致纤毛不动综合征，明显特征为呼吸道反复感染[29]。与氧化应激相关的KNDC1和LIPG[30-31]，与血管生成相关的C1QL3和与造血相关的ETV7均发生了显著上调[32-33]。但是，与神经系统相关的多个基因则发生了显著下调，如诱导神经元突起生长的重要导向分子Slit1[34]，诱导神经元微管卷曲的腔内蛋白MAP6[35]，与阿尔兹海默症相关的TREM2[36]，与突触可塑性相关的Irfn2[37]等。促进细胞增殖的DUSP4、HK2和S100A4也发生了显著下调。

与糖代谢相关的基因，如调节糖异生过程的限速酶PCK1[38]和对葡萄糖消耗和脂肪酸转换起重要作用的PDK4[39]发生上调，但介导果糖代谢的ALDOB则发生下调。与脂类代谢相关的基因，仅发现了ANGPTL4、KLF11、NUPR1、HERC5等发生了上调，差异倍数在2～3倍之间，但脂解和白色脂肪褐色化的一些标志基因并未检测到显著变化。促进脂肪代谢的UCP3上调了1.26倍，但调节脂肪储存、加快新陈代谢的Leptin却下降了3.2倍。可见，冷刺激的确改变了民猪皮下脂肪的代谢，但具体作用机制还需进一步研究。

本研究通过对lncRNA和mRNA的联合分析发现，仅有7%的靶基因受lncRNA的顺式调控，更多的靶基因是受lncRNA反式调控，这与Zhan等[40]在对山羊发育过程中骨骼肌lncRNA靶基因的预测结果基本一致(6.8%)。据此推测，lncRNA与靶基因间更多的是反式作用调控。lncRNA与mRNA构建的共表达网络中，既存在1个mRNA可以结合多个lncRNAs，也存在1个lncRNA可以与多个mRNAs发生作用，导致它们彼此之间的作用关系非常复杂。以lncRNA MSTRG.90687为例，它可以同时与6个基因发生作用，包括与脂肪代谢相关的UCP3、KLF11和PLIN5，与繁殖相关的TTC39A、ADAMTS18，与糖原代谢病相关的SLC37A4。而与脂肪代谢相关的Leptin基因也可同时与6个不同的lncRNAs (MSTRG.79145、MSTRG.109887、MSTRG.157809、MSTRG.27578、MSTRG.29107和MSTRG.67063)结合。本研究所筛选的多为新发现的lncRNAs，哪些lncRNAs 在冷应激下发挥脂肪代谢调控作用，还需要进一步的分子生物学试验验证。

4 结 论

冷应激会造成民猪皮下脂肪内lncRNA与mRNA表达水平的变化，对机体的先天免疫系统、呼吸、氧化应激、血液循环和神经系统的敏感性等生理过程造成影响，使脂类代谢发生改变。受冷应激诱导的lncRNA与靶基因之间存在复杂的互作关系。