利用RNA-seq 技术筛选火毛撒坝猪和杜洛克猪背最长肌差异基因

李再磊，窦必成，赵玉香，邱丙姗，吴正明，翁亚烦，孙婷婷，赵桂英*

（1.云南农业大学动物科学技术学院，云南昆明 650201；2.云南楚雄南华县农业农村局，云南楚雄 675200）

我国地方猪种资源是养猪生产可持续发展的基础，近年来随着组学技术的不断发展，地方猪种资源的保护及开发利用越来越受到重视。大多数云南地方猪具有适应性强和肉质好等特性，李志勋等对火毛撒坝猪、高黎贡山猪和杜洛克猪生产性能的比较研究发现，火毛撒坝猪的熟肉率极显著高于杜洛克猪，贮存损失极显著低于杜洛克猪；火毛撤坝猪的背最长肌、股二头肌肌纤维直径显著高于杜洛克猪；火毛撒坝猪肌肉的胆固醇含量极显著低于高黎贡山猪和杜洛克猪。应用RNA-seq技术研究调控中国地方猪优质性状的分子机制对猪肉质性能测定的研究已有很多，但针对杜洛克猪和撒坝猪的背最长肌差异基因研究与筛选、火毛撒坝猪优良肉质性状遗传机制的研究未见报道。本研究通过RNA-seq技术对杜洛克猪与火毛撒坝猪之间的差异表达基因进行筛选与注释，分析2 个猪种之间肉质遗传机制的差异，寻找差异表达基因及信号通路，确定调控肉质性状的关键基因，为探究云南特有地方猪种与优良外来猪种的差异以及对云南地方猪的发展提供参考依据。

1 材料与方法

1.1 材料 选用10 头70 日龄胎次和产仔日龄相近的火毛撒坝猪和杜洛克猪，相同条件下饲养至260 日龄，随机屠宰各4 头。采集背最长肌样品，同一猪种4 头猪样品混成1 份进行总RNA 提取，采用天根生化科技（北京）有限公司RNA 提取试剂盒（DP419），用2.0%的琼脂糖凝胶电泳检测RNA，Agilent 2100 RNA 6000 Nano Kit 进行浓度精确定量，将检测合格的样品进行转录组测序。

1.2 实验方法 测序reads 的比对及基因表达量统计：利用STAR 软件（版本2.5.3a）将火毛撒坝猪和杜洛克猪获得的干净读数比对到参考基因组Sus scrofa 11.1（ftp://ftp.ensembl.org/pub/release-99/variation/gvf/sus_scrofa/）上，得到全面的转录本信息，并统计比对结果。再通过featureCounts（Subread-1.5.1，Bioconductor）对映射到每个基因外显子区域的读数进行计数，然后计算RPKM（Reads per Kilobase per Million Reads）值，即每一百万个比对上基因组的读段中比对到转录本的每一千个碱基上的读段个数，作为基因表达量的衡量指标。

基因差异表达分析：差异表达基因是指样品间上调表达基因和下调表达基因的汇总。本研究使用edgeR 软件（版本3.12.1）分析基因在各样本中的差异表达情况，差异基因通过筛选阈值|logFC|≥1 且错误发现率（False Discovery Rate，FDR）≤0.05 为表达量显著差异的基因。差异表达基因的GO 富集分析：利用GO 数据库（Gene Ontology,http://www.geneontology.org）将基因按照其参与的生物过程（Biological Process，BP）、分子功能（Molecular Function，MF）及细胞组分（Cellular Component，CC）3 个方面进行分类注释。通过GO term 注释的差异基因，计算每个term 的基因列表和基因数目，应用超几何检验，找出与整个基因组背景相比，差异表达基因中显著富集的GO 条目，从而找出差异表达基因显著相关的生物学功能，＜0.05 时，该条目作为显著条目。

差异表达基因的KEGG 富集分析：KEGG（Kyoto Encyclopedia of Genes and Genomes）数据库（http://www.genome.jp/kegg/），即京都基因和基因组百科全书，是通路分析的主要数据库。利用KOBAS 软件对差异表达的基因进行通路分析，找出与整个基因组背景相比FDR＜0.05 的通路为差异表达基因显著富集的通路，其中FDR 值是校正后的值。

2 结果与分析

2.1 总RNA 质量 提取2 个猪种背最长肌组织的RNA，检测浓度，总RNA 的浓度均大于40 ng/μL（杜洛克猪150 ng/μL、火毛撒坝猪139.5 ng/μL），OD值均处于1.8～2.2，结果说明实验提取的总RNA 质量高可以用来测序。

2.2 测序数据质控及比对情况 使用Illumina HiSeq 平台进行测序，结果显示火毛撒坝猪和杜洛克猪背最长肌总reads 读取数分别为130 567 646 条、85 043 169 条，clean reads 分别为95 362 436.5、82 899 945.6，有 效读数率为74.8%、97.4%；Q30 低于0.1% 的为99.1%和92.5%，GC 碱基含量分别为50.6% 和54.1%；约95.81%的过滤后序列被定位到杜洛克猪的参考基因组，低于撒坝猪的97.13%（表1）。

表1 测序数据在参考基因组上的比对情况

2.3 差异表达基因筛选 本研究只对2 个猪种间差异显著的基因进行分析。在火毛撒坝猪中上调的基因，在杜洛克猪中是下调的，反之亦然。2 个猪种间共筛选出3 683 个差异显著表达基因，火毛撒坝猪与杜洛克猪相比，有1 516 个上调基因，2 167 个下调基因，其中火毛撒坝猪主要的上调基因有和等，下调基因有等（表2、图1）。

表2 富集于GO 和KEGG 通路内与脂肪相关的5 个相关基因

图1 火毛撒坝猪背最长肌组织差异表达基因筛选情况

2.4 GO 功能富集分析 火毛撒坝猪与杜洛克猪背最长肌上调差异表达基因中显著富集的GO 条目共有297 个，最显著的前20 个GO 条目如表3 所示。生物过程主要显著性富集在电子传递链、呼吸电子传递链、mRNA加工等条目；分子功能（MF）主要显著富集在RNA 结合等条目；细胞组分（CC）主要显著性富集在呼吸链、线粒体、核糖核蛋白复合物、细胞内细胞器部分等条目。

表3 火毛撒坝猪背最长肌上调差异表达基因显著富集的条目

对火毛撒坝猪与杜洛克猪的差异表达基因富集的条目分析，背最长肌下调差异表达基因中显著富集的GO条目共有217 个，最显著的前20 个GO 条目如表4 所示。生物过程主要显著性富集在细胞对生物刺激的反应、响应DNA 损伤的信号转导、细胞分解代谢过程的负调节、蛋白质分解代谢过程的负调节等条目；分子功能（MF）主要显著富集在核糖体的结构成分、结构分子活性、酶结合、GTP 酶结合等条目；细胞组分（CC）主要显著性富集在胞质部分、胞质核糖体、核糖体亚基、胞质大核糖体亚基、胞质小核糖体亚基等条目。

表4 火毛撒坝猪背最长肌下调差异表达基因显著富集的条目

2.5 KEGG 富集分析 如表5、6 所示，上、下调差异表达基因显著富集的通路分别为23 条、25 条，共48 条。在火毛撒坝猪背最长肌组织中，上调基因参与的23 条代谢通路中主要的代谢通路有mRNA 监测途径、脂肪酸代谢等。而下调差异表达基因参与的25 条生物学代谢通路中，主要的代谢通路有6 条与脂类代谢相关的通路，分别是脂肪酸代谢、PPAR 信号通路、不饱和脂肪酸的生物合成、脂肪酸降解、脂肪酸延长、甘油磷脂代谢。与疾病相关的通路，分别是金黄色葡萄球菌感染、利什曼病、核糖体、病毒性心肌炎、B 细胞受体信号通路、哮喘、用于IgA 生成的肠道免疫网络等。

表5 火毛撒坝猪背最长肌上调差异表达基因显著富集的通路

3 讨 论

本研究以GO 和KEGG 2 大数据库的数据资源为基础，对其进行了通路分析，有48 条显著富集通路，其中有6 条与脂类代谢相关。大部分通路内富集有关脂肪酸代谢、沉积的基因，其中和ACAD 这5 个与脂肪沉积、代谢和脂肪细胞分化相关的基因均上调。

表6 火毛撒坝猪背最长肌下调差异表达基因显著富集的通路

对火毛撒坝猪与杜洛克猪背最长肌的差异表达基因进行GO 功能富集分析，上、下调差异表达基因中，生物过程富集到的GO 条目最多。由此可以得出，生物过程是火毛撒坝猪、杜洛克猪在背最长肌中差异表达基因行使的主要生物学功能。通过GO 生物学功能和KEGG 通路富集分析表明，基因产物主要涉及IRS1/Akt/FoxO1 信号通路、腺苷50-单磷酸激活蛋白激酶（AMPK）的级联作用、脂质代谢和氨基酸代谢途径，这些途径有助于民猪的脂肪酸代谢、肌内脂肪沉积和骨骼肌生长。通过对葡萄牙阿连特雅诺（Alentejano，AL）和比萨罗（Bísaro，BI）本地猪种的脂肪组织进行转录组分析，结果表明共鉴定出458 个差异表达的基因，AL 中较高脂肪沉积的决定因素与脂肪酸从头合成的增加有关，这正是由于和等关键基因的上调所导致的；此外通路分析表明通过和，使得AL 中胆固醇合成提高。使用RNA-seq 对猪的肝脏组织进行转录组测序，结果鉴定出等55 个差异表达基因，这些基因属于与脂质和脂肪酸代谢有关的经典途径和基因网络。本研究在火毛撒坝猪和杜洛克猪中均未富集到胆固醇合成代谢直接相关的功能条目，但富集到了脂肪酰乙酰辅酶A 生物合成过程、酰基辅酶A 脱氢酶活性、基辅酶A 结合、脂肪酸-氧化、脂质氧化、细胞脂质分解代谢过程等与脂代谢相关的GO 条目，这些结果与前人研究结果相似。脂肪酸能够对胆固醇的代谢产生影响，单不饱和脂肪酸可诱导胆固醇酯化进而维持机体胆固醇在较低水平，其中棕榈油酸（C16:1）和油酸（C18:1）是胆固醇酯中含量最高的脂肪酸，对肌肉脂质沉积等有着重要作用。而乙酰辅酶A 是含有NADPH 的非脂质前体，是胆固醇从头生物合成的底物，与胆固醇关系密切。由此可见，以上这些GO 条目可能在火毛撒坝猪低胆固醇形成的分子机制中具有重要的功能。

KEGG 数据库中丰富的通路信息有助于去了解基因的生物学功能，如代谢途径、遗传信息传递以及细胞过程等一些复杂的生物功能。差异表达基因进行KEGG通路富集分析，在火毛撒坝猪的上调差异表达基因中均富集到了与脂类代谢相关的PPAR 信号通路，说明火毛撒坝猪与杜洛克猪在脂类代谢方面的差异可能是由于PPAR 信号通路的一些相关基因表达差异引起的，这与杨建敏研究杜洛克猪与槐猪的背最长肌转录组分析结果相似。

对PPAR 信号通路进行分析，和这5 个与脂肪的沉积、代谢和脂肪细胞分化的相关基因均上调，且这些基因的表达量均是火毛撒坝猪高于杜洛克猪。李优磊等对调控猪皮下脂肪与肌内脂肪差异沉积中PPAR 信号通路的作用进行研究，结果显示富集指数是PPAR 信号通路最高。长链脂酰辅酶A 合成酶（）是生物合成途径中不可或缺的。对脂肪细胞分化关键转录因子的活性起到调节的长链脂酰辅酶A 合成酶3（）参与脂肪酸的吸收。Bu 等使用RNA 干扰技术抑制大鼠原代肝脏细胞基因的表达，结果显示基因显著降低了等脂肪生成相关的转录因子的活性。在单不饱和脂肪酸（MUFA）的形成中，硬脂酰辅酶A 去饱和酶（）会限制其速度。现有研究显示有2 个亚型，其中主要在脂肪组织中表达，主要在脑组织中表达，但在脑外组织的功能尚未阐明。陶璇等采用荧光定量PCR（qRT-PCR）方法检测在藏猪和杜洛克猪里脊中关于IMF 沉积相关基因的表达差异，结果表明：180日龄的藏猪mRNA 表达量极显著高于杜洛克猪。火毛撒坝猪背最长肌中基因表达量也高于杜洛克猪，在脂肪沉积能力较强的猪种上基因均显著富集，说明基因可能是调控猪脂肪沉积的重要候选基因。脂肪酸结合蛋白（）家族的过表达会增加脂肪酸的转运。为肝脏型脂肪酸结合蛋白，Bamester 等对34 头猪的肝脏组织做共表达网络和表达基因组全关联分析，鉴定猪肉品质性状功能调控机制中涉及的DNA 变异和分子途径，通过DNA 变异与基因多重相互作用发现基因可能决定背膘中C18:2（n-6）脂肪酸与C18:3（n-3）的比例；FABP4为脂肪性脂肪酸结合蛋白，陈滇黔等在研究山羊肉质相关主效基因时得到FABP4 为其影响肉质的主效基因。李宝钧等研究发现在绵羊脂肪组织中的mRNA 表达最高，与其他组织差异显著。脂肪酸结合蛋白7（）位于PPAR 通路中，在细胞生长分化、脂质代谢过程中发挥作用。在本实验中显著富集的通路中大多筛选出脂肪酸结合蛋白（）家族的基因，且均在火毛撒坝猪中有较高表达，与以往的研究一致。乙酰辅酶A 脱氢酶（）家族催化脂肪酸氧化，长链辅酶A 脱氢酶（）在长链脂肪酸-氧化中起着表达调控和活性调节的作用。乔木等对长白猪、大约克夏猪和梅山猪进行基因分型和脂肪沉积性状关联分析，得出基因与6～7 腰椎间膘厚、肩臀膘厚和背最长肌肌内脂肪含量存在显著相关（＜0.05）的结论。戢爽研究显示延黄牛和延边黄牛的脂肪组织和背最长肌中基因表达差异不显著。以上结果提示PPAR 信号通路可能是火毛撒坝猪脂肪沉积能力的关键通路。

火毛撒坝猪下调差异表达基因主要富集在以下几个通路：利什曼病、病毒性心肌炎、哮喘、用于IgA 生成的肠道免疫网络等，这些通路与疾病相关，说明火毛撒坝猪具有较好的抗病力，该结果为地方猪抗病育种研究提供了方向。

4 结 论

本次研究利用RNA-seq 技术对火毛撒坝猪和杜洛克猪背最长肌的转录组进行测序分析，揭示了撒坝猪和杜洛克猪背最长肌的整体表达特征，共筛选出3 683 个差异基因。火毛撒坝猪和杜洛克猪相比，有1 516 个上调基因，2 167 个下调基因。本研究找出了与肉质性状显著相关的脂类代谢等代谢通路，同时发现与金黄色葡萄球菌感染、利什曼病有关的生物学通路其差异表达基因下调，说明火毛撒坝猪具有更强的抗病力。