低温胁迫下民猪骨骼肌的转录调控分析

张冬杰，汪 亮，马 红，李忠秋，王文涛，刘 娣

(黑龙江省农业科学院畜牧研究所，哈尔滨 150086)

低温寒潮是自然界中最常见的一种自然灾害，对动物、植物、微生物等都会造成破坏性影响。家养畜禽虽有人为保护，能将这种伤害程度减轻，但目前的规模化养殖模式，以及引入品种为主的生产方式，都使得畜禽对环境温度的敏感性提高。猪是我国重要的家养畜种之一，与牛、羊等其他大型家养牲畜相比，其对环境温度的变化更为敏感。1998年Lossec等界定了不同生长阶段猪的临界温度，其中断奶仔猪为26～28 ℃，体重60 kg育肥猪为18～20 ℃，低于这些温度区间，均被视为低温环境。低温环境是导致新生仔猪死亡和发病的主要原因，会造成围产期母猪死亡率上升、肉品质下降。因此，猪与环境温度之间的关系依旧是养殖业不容忽视的问题之一。

在过去的20年里，测序技术的飞速进步使得人们能够研究各种内部和外部刺激下基因表达变化的全基因组模式，人们发现低温胁迫对每个物种的转录组均产生了广泛的影响。在对一种南极细菌的深度测序中发现，与初级代谢相关的基因在低温下受到抑制，而与转录调控蛋白和信号转导蛋白相关的基因受到促进，乙醇氧化途径对细菌在低温下生长起着重要作用。对葡萄柚低温耐受性的全基因组分析显示，光合作用、防御、细胞壁和次级代谢相关转录本下调，而膜蛋白、脂质代谢、植物激素和冷响应转录因子上调。短期冷暴露(3 d，4 ℃)改变了小鼠腹股沟部白色脂肪(iWAT)中参与脂肪酸伸长、三酰基甘油、鞘脂和甘油三酯合成的基因表达与生物学通路，此外，与代谢综合征、神经退行性疾病和衰老相关的基因及表达途径均发生了改善性变化。遭受重度低温的大鼠髂腰肌中有91种mRNAs 的表达量是正常、轻度以及中度低温组的2倍 以上，这些mRNAs参与了一系列的生物学过程，包括对应激和脂质的反应以及细胞对缺氧的反应，发现、、41和4是重度低温下特异表达的基因。对冷驯化下人的股外侧肌进行转录组测序发现，冷驯化改变了756个基因的表达量，表达量变化最高的基因参与了神经细胞和肌肉细胞之间的收缩和信号转导，这些上调基因与运动训练造成的基因上调存在显著重叠，尤其是与细胞外基质重塑相关的基因和通路。

民猪是生活在我国东北地区的一个古老的地方猪种，受环境气候影响，它具有显著的耐寒特性，能够在-30～30 ℃的广泛温度区间内保持正常的生理反应。本研究的目的是筛选和鉴定民猪骨骼肌内参与低温胁迫应答的基因及调控因子。为此，对经历了58 d低温胁迫的民猪骨骼肌组织进行了RNA-Seq分析。

1 材料与方法

1.1 试验材料

本试验在黑龙江省农业科学院畜牧研究所民猪养殖场进行。2020年11月9日，将6头6月龄体重相近的雌性民猪随机分成2组，每组3头个体。试验开始前，6头个体全部在室内有供暖的猪舍内饲养，温度控制在(18±2) ℃，试验开始后，一组置于室外的半敞篷舍内饲养，温度随外界环境温度变化而变化，第1天的环境温度为5 ℃/-5 ℃(最高温度/最低温度)，饲养至2021年1月6日结束，此时的环境温度为-15 ℃/-24 ℃(图1)，相同时间内，另一组始终留在温暖的舍内饲养。两组均保证自由采食和饮水，共处理58 d。试验结束后，屠宰全部个体，取背最长肌置于液氮中保存备用。

图1 试验期间环境温度变化情况Fig.1 Environment temperature changes during the experiment

1.2 民猪背最长肌组织总RNA的提取

采用常规的酚-氯仿法提取总RNA，使用NanoDrop 2000&8000 微量分光光度计对所提取的RNA进行纯度检测，要求OD/ OD≥ 1.8，使用安捷伦的2100 Bioanalyzer进行浓度与完整性检测，要求浓度不低于100 ng·μL，28 S / 18 S ≥ 1.5。

1.3 文库构建和上机测序

每个样品取1 μg总RNA作为起始量构建lncRNA 文库。首先使用Ribo-off rRNA Depletion试剂盒去除rRNA，然后向反应体系中加入碎片缓冲液(fragmentation buffer)使RNA片断化成为短片段，再以片断后的RNA为模板，利用随机引物进行逆转录，实现cDNA第一链合成；随后加入第二链标记缓冲液(2nd Strand Marking Buffer)、第二链/末端修复混合酶(2nd strand/end repair enzyme mix) 合成cDNA第二链，经末端修复、加碱基A、加测序接头后，通过磁珠筛选回收目的片段，加UNG酶消化cDNA二链，并进行PCR扩增，最后通过磁珠纯化回收目的片段，从而完成整个文库制备工作。使用NovaSeq 6000 S4对构建好的文库进行双末端测序。

1.4 民猪背最长肌转录组测序后原始数据的处理

测序得到的原始下机序列含有测序接头序列以及低质量序列，需要对这些序列进行过滤，包括去除接头污染的Reads，低质量的Reads，含N比例大于5%的Reads，以及与rRNA匹配的Reads，得到高质量的Clean Reads后再进行后续分析。

1.5 民猪背最长肌中lncRNA的鉴定

首先对转录本进行筛选，保留那些长度≥200 bp，外显子个数≥2且reads覆盖度≥5的转录本；然后利用在线软件gffcompare同猪的注释文件进行比较，去掉与mRNA及其他非编码RNA(rRNA、tRNA、 snoRNA、snRNA等)相匹配的转录本；根据比较结果中的class_code信息(“u”、“i”、“x”)筛选潜在的lincRNA、 intronic lncRNA和anti-sense lncRNA； 最后利用CNCI、CPC、PFAM和CPAT这4个分析软件对筛选出的lncRNA是否具有编码潜能进行预测，利用这4个软件均判别为非编码(non-coding)的转录本定义为最终的lncRNA数据集。

1.6 民猪背最长肌中基因和lncRNA的表达量分析

使用FPKM定量估计基因表达值(附文献)，采用DEseq2进行试验组和对照组间差异表达基因分析，满足|logfold change|≥1和q<0.05的差异基因被认为达到显著水平。根据两组间上、下调基因情况绘制差异表达基因的火山图。

1.7 民猪背最长肌中受低温胁迫诱导的基因功能分析

GO富集分析：针对GO数据库中的二级条目，计算每个条目的基因数目，应用超几何检验找出与整个基因组背景相比，差异表达基因显著富集的GO条目。通路富集分析：对KEGG中每个通路应用超几何检验进行富集分析，找出差异表达基因显著性富集的通路。

1.8 民猪背最长肌中受低温胁迫诱导的lncRNA与mRNA联合分析

寻找差异表达lncRNA的靶基因，进而通过靶基因间接预测其生物学功能。目前已知lncRNA和靶基因间主要通过顺式()和反式()两种方式进行调控，通过lncRNA的功能与其在基因组座位上临近的蛋白编码基因的相关性进行调控靶基因的预测，将lncRNA相邻位置(上、下游50 kb)的蛋白编码基因筛选出来作为靶基因。通过lncRNA 的功能不依赖于和编码基因的位置关系，而是与共表达基因的相关性进行靶基因的预测。根据lncRNA同mRNA表达量的相关性系数(斯皮尔曼相关系数，corr ≥0.9)进行筛选。筛选后获得的靶基因开展功能注释、GO富集和通路富集分析。

1.9 差异表达基因的real-time PCR验证

利用Roche实时荧光定量PCR仪LightCycler480Ⅱ对试验组显著高表达的9个基因和2个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法计算各模板中目的基因相对于内参基因-的表达量。

表1 qRT-PCR检测用引物信息

2 结 果

2.1 低温环境下民猪背最长肌的转录组测序结果

原始下机数据经过滤处理后，每头个体平均获得98 M的clean reads，与基因组的比对效率为96.12%，其中，38.94%的序列比对到外显子区域，30.50%的序列比对到内含子区域，30.56%的序列比对到基因间区,筛选出的lncRNA经CNCI、CPC、PFAM蛋白结构域和CPAT分析后，取四者间的交集，共获得10 220个lncRNAs。

2.2 低温环境下民猪背最长肌中差异表达基因的筛选与分析

民猪在经历58 d的低温胁迫后，背最长肌中共有102个基因发生了显著变化，其中86个基因发生了显著上调，16个基因发生了显著下调(|logfoldchange|≥1，q<0.05)。186个lncRNAs发生了显著变化，其中112个lncRNAs发生了显著上调，74个 lncRNAs发生了显著下调(图2)。由于猪的lncRNA数据库还不完善，因此这186个lncRNAs均显示为novel lncRNA。发生显著变化的、且有注释的蛋白编码基因共计86个。其中8个基因的变化倍数在3倍以上，27个基因的变化倍数在2～3倍之间，51个基因的变化倍数在1～2倍之间。上调变化倍数最大的是神经降压素受体2(neurotensin receptor 2，2)基因，上调了7.18倍；下调变化倍数最大的是RAS样蛋白家族11A(RAS like family 11 member A，11A)基因，下调倍数为2.75倍(表2)。

图2 低温胁迫下民猪背最长肌内差异表达基因的火山图Fig.2 Volcanic plot of differentially expressed genes in the longissimus dorsi of Min pigs under low temperature stress

表2 低温胁迫下民猪背最长肌内发生显著变化的基因

2.3 民猪低温环境下背最长肌中差异表达基因的功能分析

对低温环境下民猪背最长肌中差异表达基因进行GO富集分析，发现94个差异基因显著富集到94个生物过程条目、12个分子功能条目和2个细胞组分条目(图3，仅列出前10个条目)。在生物过程(biological process, BP)这个本体中，富集程度最高的是骨骼肌细胞分化过程(skeletal muscle cell differentiation)、学习记忆(learning of memory)和MAPKKK活性的激活(activation of MAPKKK activity)；在分子功能(molecular function, MF)这个本体中，富集程度最高的是RNA聚合酶II近端启动子序列特异性DNA结合(RNA polymerase II proximal promoter sequence-specific DNA binding)、近端启动子序列特异性DNA结合(proximal promoter sequence-specific DNA binding)和RNA聚合酶II调节区DNA结合(RNA polymerase II regulatory region DNA binding)；在细胞组分(cellular component, CC)这个本体中，仅有2个条目被显著富集，分别是CHOP-ATF3复合体(CHOP-ATF3 complex)和CHOP-C/EBP复合体(CHOP-C/EBP complex)。

图3 差异表达基因的GO富集分析Fig.3 Go enrichment analysis of differentially expressed genes

对差异表达基因进行通路富集分析，发现仅3条通路被显著富集，分别是细胞凋亡通路(apoptosis, q=0.001 98)、直肠癌通路(colorectal cancer, q=0.001 98) 和癌症中的转录误调节通路(transcriptional misregulation in cancer, q=0.036 44)。在细胞凋亡通路中，7个基因、、45B、45G、45A、3、2L11发生了显著上调。

2.4 民猪低温环境下背最长肌中差异表达lncRNA的靶基因预测

lncRNA不编码蛋白，主要通过或方式作用于编码蛋白基因。本研究发现，民猪在经历58 d的低温环境后，背最长肌内186个差异表达的lncRNAs顺式调控251个靶基因，反式调控1 377个靶基因。通过对lncRNA与mRNA之间的互作分析发现两者间的调控关系错综复杂。比如MSTRG.13828同时调控67个基因，其中3个为调控，64个为调控，它所调控的基因215在两组间达到了显著差异；基因2AIP1则同时受到14个lncRNAs的反式调控。表3列出了变化水平前10位的上调表达lncRNAs的靶基因情况及调控方式。

表3 低温胁迫下民猪背最长肌内变化倍数前10位的lncRNAs和其靶基因情况

2.5 民猪低温环境下背最长肌中差异表达lncRNA的靶基因功能分析

对差异表达lncRNA的靶基因进行GO功能注释，发现在生物过程本体中，有81个条目被显著富集，其中富集程度最高的是核酸代谢过程(nucleic acid metabolic process)，在分子功能本体中有28个条目被显著富集，其中核酸结合(nucleic acid binding)富集程度最高；在细胞组分本体中，有19个条目被显著富集，其中富集程度最高的在细胞内部分(intracellular part)，每个本体中富集程度前10的条目情况见图4。

A. 生物过程显著富集的前10个条目结果; B.分子功能显著富集的前10个条目结果; C.细胞组分显著富集的前10个条目结果。横坐标表示富集程度，纵坐标表示GO条目。每个点表示该GO条目的富集程度，颜色越趋近于红色表示富集程度越高。每个点的大小表示富集到该GO条目的基因的个数，点越大表示富集到该GO条目的基因越多，反之则越少。Rich ratio的计算公式为：(该通路的差异基因/所有的差异基因)/(注释到该通路的基因/所有能被注释到的基因)A. Results of the top 10 terms with significant enrichment in biological processes; B. Results of the top 10 terms with significant enrichment in molecular function; C. Results of the top 10 terms with significant enrichment in cellular components. The abscissa represents the degree of enrichment and the ordinate represents the GO term. Each point indicates the enrichment degree of the GO term， the closer the color to red, the higher the enrichment degree. The size of each point indicates the number of genes enriched into the GO term, the larger the point, the more genes enriched into the GO term, and vice versa. Rich ratio: (differential genes in this pathway / all differential genes) / (genes annotated to this pathway / all genes that can be annotated)图4 GO条目FDR值富集图Fig.4 FDR value enrichment map of GO terms

对差异表达lncRNA的靶基因进行通路富集分析，发现单纯疱疹病毒1感染通路(herpes simplex virus 1 infection)、MAPK信号通路(MAPK signaling pathway)和范可尼贫血通路(fanconi anemia pathway)被显著富集。

2.6 测序结果的实时定量PCR验证

为了验证RNA测序结果的准确性，选择了8个 显著上调表达的基因、、5、1、1、1、1、1和6个显著上调表达的lncRNAs，2个显著下调表达的基因146和3以及1个显著下调表达的lncRNA73004进行了qRT-PCR试验验证。结果显示，qPCR定量检测值与测序值相比，变化趋势基本一致(图5)，表明测序结果可信。

图5 qRT-PCR验证结果Fig.5 qRT-PCR validation results

3 讨 论

自然界的动植物都会受到来自环境的低温胁迫影响，植物对低温胁迫的响应包括膜脂类型和去饱和水平的变化。动物因其生理结构和遗传基础比植物复杂的多，因此其低温耐受机理还未完全解析清楚。目前研究主要集中在两个方面，一个是低温环境下白色脂肪的褐色化机理，另一个是低温环境下骨骼肌的产热机理。骨骼肌是哺乳动物成年个体机体内最丰富的组织，产热速率的微小变化就会引起全身产热的重大变化。在对小鼠、鱼类和鸟类的研究中发现，肌浆网Ca循环是一种不会干扰ATP合成的产热机制。骨骼肌内还存在丰富的神经组织，能够对循环的肾上腺素能因子作出反应；富含线粒体，具有较高的脂质氧化能力。缺乏褐色脂肪组织的动物(鸟类和猪)或褐色脂肪含量较少的动物(成人和大型哺乳动物)，都以骨骼肌内依赖于肌质蛋白的非战栗性产热为主。

本研究以经历58 d低温胁迫的民猪骨骼肌为试验材料，采用RNA-seq技术对其进行差异表达基因的筛选与分析，发现4个与神经系统相关的基因发生了显著上调，如上调倍数最高的2基因，它表达水平的增加可使小鼠对热痛觉和化学痛觉产生抵抗力。紧随其后的基因，是一种神经元基因，对于哺乳动物大脑中持久的信息储存至关重要，并且与各种神经系统疾病有关。1作为一种转录因子，在神经元细胞损伤以及脑损伤中发挥作用，它的高表达会抑制AMPK信号通路的活性，进而抑制神经元细胞活性、促进炎症发生、加速细胞凋亡。1是一种内源性钙调磷酸酶抑制剂，在唐氏综合征(人类神经系统的最常见染色体疾病)中的表达增加了3倍，它通过抑制神经生长因子受体的内吞作用损害交感神经元的营养支持。

本研究中，溶质转运蛋白超家族(solute carrier，SLC)的7个基因受到低温诱导，发生了显著上调，包括2A4、2A5、4A10、19A2、20A1、28A1和38A2。该家族介导细胞与外界或细胞内部各类溶质的跨膜运输，调节基本化学物质(包括水分、营养素、激素和许多药物)的跨膜运输，维持细胞内稳态，保持细胞膜的完整性，以保护细胞的代谢活动和遗传信息免受环境影响。SLC转运蛋白包含52个基因家族和近400个不同的转运蛋白基因。不同家族的SLC负责转运不同的基质，比如SLC2家族主要负责转运葡萄糖，SLC4家族主要负责碱基的转运，SLC19家族主要介导叶酸和硫铵两个重要水溶性维生素的运输，SLC20负责转运钠依赖性磷酸盐，以维持体内的磷酸盐平衡，SLC28家族专门负责核苷酸转运，SLC38家族的转运体存在于身体的所有细胞类型中，介导依赖Na的小中性氨基酸的净摄取和流出。由此可见，长时间的低温胁迫使机体内多种物质的转运发生了显著变化。

多个与炎症和自身免疫反应相关的基因也发生了显著上调，如基因，它是表皮生长因子受体(EGFR)的配体，炎症脂质、细胞因子和激素等可刺激它的表达与释放，表达的慢性升高与不同病理条件相关，主要是炎症和肿瘤。基因属于细胞因子信号转导抑制因子(SOCS)家族，会抑制IL-15介导的JAK-STAT信号传导，它的上调表达会使先天免疫力下降。1具有免疫调节功能，在炎症状态下上调表达，特异性抑制IRF3的K63型多泛素化修饰，最终缓解体内过强的炎症反应、降低人体发生自身免疫性疾病的风险。1可促进肿瘤细胞增殖，抑制凋亡，还可以通过影响炎症因子的分泌来调节急慢性炎症。此外，上调倍数在4倍以上的33、都是与免疫相关的基因。由此可知，低温胁迫导致骨骼肌出现炎症反应，同时还影响了机体的免疫力。

此外，低温胁迫还可诱发组织的氧化应激反应，氧化应激被认为会促进肌肉萎缩。本研究发现，与氧化应激反应相关的基因也受到了诱导，如1、1、3、和等。

民猪在遭受低温胁迫后，肌肉组织内发生显著下调的基因共计16个，在数量上远远少于发生显著上调的基因。在功能上则较为分散，有抑制细胞增殖的11A基因，参与细胞膜结构组成的跨膜蛋白139(139)基因，具有调节TGF-β生物活性功能的潜在转化生长因子结合蛋白2(2)基因，参与糖代谢和神经节苷脂分解代谢的唾液酸酶(3)基因以及与昼夜节律相关的芳香烃受体核转运体样蛋白1()基因等。Dopico等在对16 000多人的血液和脂肪样本分析中发现，人类基因中有将近1/4的基因活性(5 136/22 822)会随着季节的变化而不同，一些基因在冬季更活跃，而另一些则在夏季更活跃。本研究发现,显著下调的基因表现出夏季活跃，而冬季不活跃的现象。

本研究中对差异基因的GO富集分析发现，骨骼肌细胞分化过程的富集程度最高，该生物过程中的5个基因发生了显著上调，分别为转录激活因子3(3)、MAF bZIP转录因子F()、Kelch家族成员40(40)、心肌锚蛋白重复域1(1)基因和生肌因子6(6)。低温环境会抑制猪只生长，增加机体能量代谢。Shima和Matsuda研究发现，小鼠和人类的成肌细胞在低温(30 ℃)下培养时不能表达肌生成素，肌细胞分化受到抑制。而且，来源于不同部位的肌细胞对低温环境的应答反应也存在差异。但本课题组在前期研究中发现，民猪在低温环境下并未出现明显的生长抑制现象，推测可能与该生物过程的活化相关。

对差异基因的通路富集分析发现，仅有3条通路被显著富集，分别为细胞凋亡、直肠癌、癌症中的转录误调节通路。癌症的本质是细胞分化和增殖异常、生长失去控制，除了2条富集到癌症相关通路上，还有1条是富集到细胞凋亡通路，这说明低温胁迫影响了骨骼肌细胞的正常生长。细胞凋亡通路中有3个45基因被富集，GADD45蛋白可参与DNA修复、细胞周期调控、衰老和基因毒性应激等多种细胞功能的调控，它们的上调代表着机体启动了对细胞的保护和修复功能。

在对差异表达lncRNA的分析发现，由于猪的lncRNA数据库还不完善，所以无法对它们进行功能注释，但在本研究中发现，lncRNA对靶基因的调控多以反式调控为主，且与靶基因间存在复杂的调控关系。对这些lncRNA预测的靶基因进行功能注释时发现，这些靶基因显著富集在3条生物学通路，其中的MAPK通路是信号从细胞表面传导到细胞核内部的重要传递者，包括MKKK、MKK和MAPK三种激酶，这三种激酶能够依次激活，共同调节着细胞的生长、分化、对环境的应激适应、炎症反应等多种重要的细胞生理/病理过程。本研究中有36个基因富集到该通路上，其中24个基因上调，12个基因下调，发生变化的基因主要集中在MAPK通路的JNK和p38/MAPK分支路线上，其中p38/MAPK是促进骨骼肌运动代谢适应的信号通路之一。据此推测，低温胁迫下骨骼肌内的MAPK通路受到lncRNA的转录调控。

4 结 论

低温胁迫对民猪的神经系统和免疫系统造成了影响，使细胞内的物质转运发生了变化，细胞的生长、分化也受到了影响，大量的lncRNA参与了骨骼肌的低温胁迫应答反应。但发生显著变化的基因差异倍数多集中在1～2倍之间，而且只有细胞凋亡通路、直肠癌通路和癌症中的转录误调节通路受到影响，推测低温胁迫并未对民猪的骨骼肌造成严重损伤。