基于RNA-Seq 技术鉴定撒坝猪与杜洛克猪肝脏组织差异表达基因

吴正明，邱丙姗，李志勋，程 娜，翁亚烦，李再磊，赵桂英*

（云南农业大学动物科学技术学院，云南昆明 650201）

我国很多地区都分布有独具特色的地方猪种资源，它为我国养猪生产可持续发展奠定了基础，也是满足人类需求的一份重要基因库。近年来，随着组学技术的不断发展，地方猪种资源的保护及开发利用越来越受到重视。有相关研究表明，大多数中国地方猪具有肉质好、适应性强和繁殖力高等特性[1]。研究调控中国地方猪种优质性状的分子机制，发掘和鉴定中国地方猪种优质性状的主效基因对中国地方猪种资源的保护及开发利用具有重要意义。转录组（Transcriptome）指某一物种在某一特定生理状态下的单个细胞或细胞群体内几乎全部转录本，即转录后所有mRNA 的总称[2]，它不仅是研究基因结构、功能和不同基因间表达差异的基础，也是探索未知功能基因的方法之一[3]。转录组测序即RNASeq（Ribonucleic Acid Sequencing），是一种与高通量测序技术相结合从而对转录组进行分析的技术[4]，其具有高灵敏度、成本低和不受物种限制等特点，是一种被运用在转录组学上比较精确的测定分析方法[5]。现今RNA-Seq 普遍运用在各种畜禽重要分子机制和鉴定相关功能基因的研究领域。基于RNA-Seq 技术，研究者采用了不同的策略开展了不同品种猪各组织间存在的差异表达基因（Differentially Expressed Genes，DEGs）的鉴定和筛选，并对不同猪种间的种质差异进行相关分析。Liu 等[6]利用RNA-Seq 技术对民猪和大白猪背最长肌建立了转录组文库，通过基因表达图谱鉴定了1 371 个DEGs，为研究不同品种猪肉质差异的形成机制提供了参考。Ropka 等[7]通过对5 个品种944 个样本的RNA-Seq 数据分析，找到了2 个与肌纤维直径和密度极显著相关的单核苷酸多态性位点（Single Nucleotide Polymorphisms，SNPs），说明基于RNA-Seq数据的SNPs 分析可用于寻找和定位基因。Liu 等[8]采用RNA-Seq 技术分别对妊娠期不同日龄（40、55、63、70、90 d）通城猪与约克夏猪的骨骼肌转录图谱进行了综合比较，分别检测到了1 317 个和691 个DEGs，结果显示55 d 是导致通城猪和约克夏猪发育差异的关键阶段，并阐明了2 个品种间发育差异的6 个关键候选基因。Huang 等[9]比较了金华猪和长白猪肝脏中环状RNA 的差异表达情况，挖掘到336 个DEGs，结合GO（Gene Ontology）分析发现其宿主基因大部分与代谢途径相关。上述研究为转录组层次阐明猪种间的差异表达奠定了分子基础。尽管通过RNA-Seq 技术挖掘了大批与生产和肉质性状相关的候选基因，但只在国外猪种和我国少数地方猪种的研究中有相关报道，而云南地方猪种的相关研究更是鲜有报道，大量优良地方猪种的遗传基础仍很薄弱，仍需深入研究。本研究拟选用生产性能和肉质性能存在较大差异的云南地方猪——撒坝猪（肉脂兼用型品质）和国外引进的杜洛克猪（瘦肉型品种）为材料，采用RNA-Seq 技术筛选撒坝猪和杜洛克猪肝脏组织的DEGs，并对筛选的DEGs 进行GO 功能富集分析、KEGG 通路富集分析和可变剪接分析，对调控撒坝猪生长性状和肉质性状的相关候选基因进行挖掘，为进一步阐明中国地方猪种生产及肉质性能遗传机制及探索肝脏组织的转录谱提供理论依据。

1 材料和方法

1.1 实验动物和样品制备 选择胎次、产仔日龄相近的撒坝猪和杜洛克猪肥猪各4 头，饲养于云南农业大学动物科学技术学院实习猪场，在相同饲养条件下饲养至260日龄（撒坝猪体重近90 kg，杜洛克猪体重近110 kg）后屠宰，屠宰后用无RNA 酶的冻存管采集肝脏样品，样品在低温条件下运回实验室，并存于-80℃的冰箱中。

1.2 总RNA 的提取和文库的构建 用Trizol 试剂分别提取每头撒坝猪和杜洛克猪肝脏组织中的总RNA；Nanodrop、琼脂糖凝胶电泳、Qubit、Aglient 2100 等软件检测总RNA 的浓度、纯度和完整性，提取完整性和纯度检测符合要求后的总RNA 样本。向Oligo（dT）的磁珠富集过的mRNA 中加入fragmentationbuffer，把mRNA 打断片段。以mRNA 为模板、六碱基为随机引物合成第1 链cDNA；加入dNTPs、缓冲液、DNA polymerase I 和RNase H 合成第2 链cDNA。先 对cDNA 进行纯化并洗脱、做末端修复和连接测序接头，然后用琼脂糖凝胶电泳进行片段大小选择，最后进行PCR 扩增富集cDNA 即完成文库制备。

1.3 数据质控 为了知晓测序数据质量的优劣，对在Illumina 测序平台上完成测序的转录组数据下机后进行质控，质控包含3 个方面：碱基质量的分布、碱基平衡性和重复序列水平。对原始序列进行质量评估的目的是通过评估数据质量的情况来引导是否有必要进行下游分析。碱基平衡性分析可以清晰检测到数据中是否具有AT、GC 分散情况，正常情况下4 种碱基出现频率是一致的，并且位置上也没有差异，当个别位置碱基的比例发生偏斜现象时，往往提示有超出比例的核酸污染。测序深度越高，越容易产生一定程度的重复，这属于正常现象；但如果重复的程度很高，可能有大量核酸受到了污染。

1.4 基因差异表达与富集分析 以每百万序列中来自于某基因每千碱基长度的序列数（Reads Per Kilobase per Million reads，RPKM）作为基因表达量的衡量水平，计算公式如下：

为了对假阳性率（False Positive rate）的控制，需结合假发现率（False Discovery Rate，FDR）即P-value和RPKM 值倍数变化（Fold Change，FC）对DEGs 进行筛选。以|logFC|>1 且FDR<0.05 为条件筛选出撒坝猪与杜洛克猪肝脏组织表达量显著差异的基因。为了比较2 个猪种间基因表达量的差异，以lgRPKM 值进行聚类分析。使用edgeR[7]对撒坝猪与杜洛克猪肝脏组织转录组数据进行差异分析，并利用DESeq[8]软件对DEGs进行GO 富集分析和KEGG 通路显著富集分析。

1.5 差异剪接事件检测 rMATS 软件具有对转录组数据进行自动检测和差异可变剪接分析的优点，因此可用该软件对成对样品进行可变剪接事件分析，主要通过rMATS 统计模型的方式检测不同样品的可变剪切类别及相应的表达量（ψ），并设定|Δψ|=|ψ1 − ψ2|<0.05 为差异表达阈值。

2 结果

2.1 转录组测序数据 利用RNA-Seq 进行撒坝猪和杜洛克肝脏组织转录组测序，每个文库产生1 亿多个原始序列，过滤掉低质量、带接头和重复的原始序列后，8 个RNA-Seq 文库中均获得8 900 万个以上的过滤序列，经过滤处理后，约有97.13%的序列被定位到杜洛克猪的参考基因组（Sus scrofa 11.1），低于撒坝猪的97.34%。杜洛克猪的外显子、内含子和基因组间区域分别占85.83%、3.62%和10.28%。而约克夏猪区域的外显子、内含子和基因组区间分别占89.56%、4.23% 和0.588%。2 种猪所比对到参考基因组上的序列比例错误率小于0.001（Q30）的均大于99%。结果表明：RNASeq 得到的数据利用率高，所选择的参考基因组能够满足本实验后续分析的需求。

2.2 筛选DEGs 从差异基因火山图（图1）中可以非常直观并且合理地筛选出撒坝猪和杜洛克猪肝脏组织的DEGs。其中x 轴为lgFC 值，y 轴为lgFDR 值。结果显示从撒坝猪和杜洛克猪肝脏组织间共鉴定出2 326 个差异表达基因，在撒坝猪上调显著表达基因有1 049 个，在杜洛克猪显著上调表达的基因有1 277 个。

本研究同时还通过热图显示了撒坝猪和杜洛克猪差异表达基因的表达水平（图2），可以清楚地观察到这2 个品种猪肝脏组织的基因表达量存在差异。

图2 差异表达基因热图

2.3 DEGs 的功能分析 将2 猪种筛选到的2 326 个表达显著差异的基因与GO 数据库和KEGG 数据库进行对比。在2 个数据库中共有926 个DEGs 被注释，其中在GO 数据库中共检测到521 个差异表达基因，在KEGG数据库中共检测到651 个差异表达基因。

2.3.1 GO 富集分析 将521 个DEGs 注释到GO 数据库进行功能分类，分类内容包括：分子功能（Molecular Function，MF）、生物过程（Biological Process，BP）及细胞组分（Cellular Component，CC）3 个方面。其中MF 分成7 个小类，BP 分成33 个小类，CC 分成12 个小类。在BP 分类中参与代谢过程的DEGs 最多；在CC分类中参与调控细胞内环境和细胞外环境的DEGs 最多；在MF 分类中调控催化活性的DEGs 最多（图3）。

图3 差异表达基因GO 注释

首先，通过与整个基因组背景相比并应用超几何检验，挑选出撒坝猪与杜洛克猪肝脏组织的DEGs 显著富集的GO 条目。图4 为按P-value 值排序选择的20 个最显著GO 条目作出的GO 富集分析散点图。用-log10（P-value）表示富集的显著性（横轴），该值越大富集得越显著，纵轴表示富集到的GO 条目数。圆点大小及颜色的深浅表示该GO 条目包含的差异基因数目和富集程度。

图4 差异表达基因GO 富集分析散点图

2.3.2 KEGG 富集分析 首先，通过与整个基因组背景相比并应用超几何检验，筛选出撒坝猪和杜洛克猪肝脏组织DEGs 中显著性富集的通路。共有47 条显著富集的通路，其中有24 条通路发生上调，有23 条通路发生下调。对DEGs 做KEGG 通路富集后，统计出它们参与的KEGG 代谢通路，如图5 所示：横轴表示富集的显著性，P值越高表示富集得越显著，纵轴表示富集的KEGG 通路数。圆点大小表示KEGG 通路包含的差异基因数目，颜色深浅表示KEGG 富集因子富集的程度。这里按P值进行排序，挑选显著性最高的前20 个KEGG 条目作出相应的KEGG 通路富集散点图；20 条有关生产和肉质性能的代谢通路具体情况如表1 所示。

表1 撒坝猪与杜洛克猪之间表达差异显著的有关生产性能及肉质性能的相关通路

图5 差异表达基因KEGG 富集分析散点图

2.4 可变剪接分析 本研究通过对撒坝猪和杜洛克猪肝脏组织DEGs 可变剪接事件的检测发现，共有5 种可变剪接事件类型：外显子互斥（Mutually Exclusive Exons，MEX）、外显子跳跃（Skipped exon，SE）、内含子保留（Retained intron，RI）、5'端可变剪接（Alter native 5’splice site，A5SS）和3'端可变剪接（Alternative 3'splice site，A3SS）。交叉点读数（Junction reads）中检测到总可变剪接事件数分别为21 487、3 489、438、746 和1 319；交叉点读数检测到的撒坝猪和杜洛克猪差异可变剪接事件比例分比为369:266、820:498、21:20、7:24 和14:22；目标读数（Reads on target）和交叉点读数检测到的总可变剪接事件数分别为21 491、3 489、439、747 和1 329；目标读数和交叉点读数检测到的撒坝猪和杜洛克猪存在差异可变剪接事件比例为389:282、817:495、20:20、9:29 和13:20（表2）。

表2 可变剪接事件分析

3 讨 论

本研究通过对GO 数据库和KEGG 数据库显著富集通路的分析，初步筛选出PPAR 信号通路等20 条有关生产和肉质性能的代谢通路，撒坝猪和杜洛克猪肝脏组织间共鉴定出2 326 个差异表达基因，其中大部分通路内富集的DEGs 与脂类代谢调控有关，包含脂蛋白酯酶（Lipoprotein Lipase,LPL）、长链脂酰辅酶A 合成酶3（Long-chain Acyl-CoA Synthetase3，ACSL3）、6-磷酸果糖激酶（6-phesphofructokinase,PFKL）、脂肪酸结合蛋白4（Fatty Acid Binding Protein 4,FABP4）、瘦素（Leptin，LEP）、脂肪转运蛋白4（Latty Acid Trans port Proteins 4，FATPs4）、脂联素（Adiponectin，ADIPOQ）、脂肪酸结合蛋白1（Fatty Acid Binding Protein1，FABP1）、长链辅酶A 脱氢酶（Long-chain Acyl-CoA dehydrogenase，ACADL）等。

肝脏是研究脂质代谢的重要组织，也是脂肪酸合成最主要的场所，因此肝脏是研究代谢过程、免疫、生长和肉质的首选组织[10]。李建平等[11]通过对育成猪饲喂不同脂肪源饲料饲养屠宰后,采肝脏组织样品做转录组差异分析，结果表明对肝脏组织差异基因的GO 富集分析和KEGG 富集分析中均有与脂类代谢相关基因被注释；Ge 等[12]选取低肌内脂肪酸和高肌内脂肪酸的巢湖鸭肝脏组织为实验样本进行转录组测序，筛选出9 个与脂肪沉积的相关基因；本研究结果和以往的研究结果大致相同，表明肝脏确实是研究脂质代谢的重要组织。

撒坝猪LPL基因表达量较高，是杜洛克猪的28.3倍，提示它和脂肪沉积能力较高的性状有关。脂蛋白脂酶（Lipoprotein Lipase，LPL）位于猪14 号染色体上，属于水解酶家族，LPL对脂质在机体内的代谢及相关代谢物的转运发挥重大作用，它能将血乳糜微粒和极低密度脂蛋白中的三酰甘油分解为脂肪酸和甘油，并将脂肪酸和甘油转运至机体相应组织中储存[13]。Wang 等[14]对藏猪和滇南小耳猪两个中国地方猪与大白和长白猪两个外来猪背最长肌进行同位素标记的相对和绝对定量，找到LPL等10 个与脂肪沉积能力相关的基因，其中藏猪-滇南小耳猪组LPL基因表达量是大白猪-长白猪组的2.84 倍，从而LPL被认为与藏族和滇南小耳猪较高的脂肪沉积能力有关；本研究结果与以往的研究结果具有一致性。

ADIPOQ基因是被人们所发现唯独与肥胖呈负相关的激素，ADIPOQ基因通过受体1（AdipoR1）、受体2（AdipoR2）以及T-黏钙蛋白相结合来调控机体内脂肪的含量，当生物体内脂肪含量过高时，就会减少ADIPOQ基因的表达量和ADIPOQ在血液中的含量[15]。Cirera 等[16]分别采集哥本哈根迷你猪和杜洛克猪背最长肌、背部脂肪、肝脏组织，测定ADIPOQ基因在2猪种间3 个组织中的表达量，结果显示：ADIPOQ基因哥本哈根迷你猪和杜洛克猪背最长肌、背部脂肪、肝脏组织中的表达量并没有明显差异，但ADIPOQ基因在杜洛克猪中表达量稍高于哥本哈根迷你猪；而本研究表明ADIPOQ基因在杜洛克猪中的表达量极显著高于撒坝猪，可能是撒坝猪脂肪沉积能力高于哥本哈根迷你猪所导致的。

PPAR 信号通路可能是影响撒坝猪脂肪沉积能力的关键通路。过氧化物酶体增值激活受体是甾体激素受体超家族成员，通过与类色素x 受体相结合而转变为异二聚体，它是调控脂肪细胞的分化、脂肪酸转运、炎症以及调节糖代谢等的相关基因[17-18]。PPARα、PPARβ和PPARγ是PPAR 的3 种不同亚型；其中PPARα和PPARβ主要在肝脏及脂肪酸氧化水平较高的组织中发挥作用，PPARγ主要作用于脂肪细胞，影响其分化和生长，也是I 型糖尿病药物治疗的靶点[19]。

4 结 论

本研究对撒坝猪和杜洛克猪肝脏组织的转录组分析，共鉴定出2 326 个DEGs，结合GO 分析和KEGG通路富集分析，本研究初步认为LPL、ACSL3、DIPOQ等9 个基因和PPAR 信号通路等16 条代谢通路在撒坝猪的肉质形成中发挥了重要作用，研究结果为揭示撒坝猪肉质和生长性状的形成机制提供了重要的参考数据，并提高了对脂质沉积相关分子机制的认识。