CN104636638A - 不同品种猪背最长肌差异表达基因的筛选与注释方法 - Google Patents

Abstract

本发明公开了不同品种猪背最长肌差异表达基因的筛选与注释方法，包括以下步骤：(1)采样与测序；(2)原始数据处理；(3)差异基因的筛选；(4)功能显著性富集分析；(5)Pathway显著性富集分析。本发明采用高通量转录组测序技术，构建不同品种猪背最长肌组织的基因表达谱，获得差异表达基因，利用GO分析和pathway分析对差异基因进行筛选，确定调控肉质性状的关键基因，加深对肌肉生长发育和脂肪代谢调控的了解，为之后研究形成国内外品种在肉质和生长性状方面的差异的分子机理提供分子生物学证据。

Description

不同品种猪背最长肌差异表达基因的筛选与注释方法

技术领域

本发明涉及一种不同品种猪背最长肌差异表达基因的筛选与注释方法。

背景技术

猪肉是人类动物蛋白的主要来源，我国是目前世界上生猪生产和猪肉消费的第一大国。改善肉品质已经成为猪肉产业及满足消费者口感、健康和营养需求的关键因素。骨骼肌是动物体内最丰富的组织，约占动物体重的40％，猪肌肉生长发育的差异是引起不同品种猪肉质差异的重要原因。由于肉质性状形成的复杂性及不能(或不易)活体测定，传统的选育方法很难对肉质性状取得较大遗传进展。因此采用分子生物学手段，深入分析肉质性状遗传机制，为其改良提供给了机会，从而参考制定合理遗传改良方案是后基因组时代家畜动物育种研究的热点。

随着后基因组时代的到来，转录组学、蛋白质组学、代谢组学等各种组学技术相继出现，其中转录组学是率先发展起来以及应用最广泛的技术。转录组是在一个或多个细胞中表达的RNA转录物的集合。通过转录组分析可以帮助我们在整体水平上研究细胞中基因转录的情况及转录调控规律。过去的转录组研究通常采用Sanger测序和cDNA芯片方法，而现在的转录组研究越来越多采用高通量测序技术。新一代高通量转录组测序(RNA-seq)与传统的基因芯片相比具有高通量、低成本、高灵敏度、重复性好且无需已知参考序列等优点，已逐步取代基因芯片成为转录组研究的主要手段。该技术已经应用于小鼠、海鲈鱼和水稻等的RNA测序。此外，2013年杨亚岚利用高通量测序技术筛选出了不同猪种间的差异表达基因并对其进行了功能注释和分析。

中外猪种比较发现，国外猪种具有生长速度快、瘦肉率高和饲料转化率高等优点，但往往存在肌内脂肪少、肉色灰白和系水力低等缺点，我国地方猪种具有肉色鲜红、系水力强、肌内脂肪含量高和肌纤维直径小等优点。

发明内容

本发明要解决的技术问题是提供一种基于RNA-seq技术对不同品种猪背最长肌差异表达基因的筛选与注释方法。

为了解决上述技术问题，本发明采用的技术方案是：不同品种猪背最长肌差异表达基因的筛选与注释方法，包括以下步骤：

(1)采样与测序：按不同猪种，采取成年期的背最长肌组织，每个猪种采集三头个体，样品采集后-80℃贮存，提取各样品总RNA，采用Solexa测序技术对样品进行测序；

(2)原始数据处理：将测序仪产生的原始图像数据经Base Calling转化为序列数据，并称之为原始序列数据，其结果以FASTQ文件格式存储；由于原始序列数据可能包含低质量序列、接头序列等杂质序列数据，不能直接用于生物信息分析，所以原始序列数据必需经过数据处理转换为高质量序列数据，利用Tophat(version：2.0.6)软件的spliced mapping算法对高质量序列数据进行猪基因组比对，采用基因组版本为Sscrofa10.2；

(3)差异基因的筛选：基因表达量的计算使FPKM法；根据基因的表达量(FPKM值)计算该基因在不同样本间的差异表达倍数；差异表达基因定义为FDR≤0.01且倍数差异在2倍及以上的基因；

(4)功能显著性富集分析：把所有差异表达基因向Gene Ontology数据库的各条目映射，计算每个条目的基因数目，然后应用超几何检验，找出与整个基因组背景相比，在差异表达基因中显著富集的GO条目，其计算公式为：

$P=1-\underset{i=0}{\overset{m-1}{\mathrm{\Σ}}}\frac{\left(\begin{array}{c}M\\ i\end{array}\right)\left(\begin{array}{c}N-M\\ n-i\end{array}\right)}{\left(\begin{array}{c}N\\ n\end{array}\right)}---\left(1\right)$

公式(1)中，N为基因组中具有GO注释的基因数目；n为N中差异表达基因的数目；M为基因组中注释为某特定GO条目的基因数目；m为注释为某特定GO条目的差异表达基因数目；计算得到的P值通过Bonferroni校正之后，以校正后的P值＜0.01为阈值，满足此条件的GO条目定义为在差异表达基因中显著富集的GO条目；通过GO功能显著性富集分析能确定差异表达基因行使的主要生物学功能；

(5)Pathway显著性富集分析：Pathway显著性富集分析以KEGG Pathway为单位，应用超几何检验，找出与整个基因组背景相比，在差异表达基因中显著性富集的Pathway；该分析的计算公式同GO功能显著性富集分析的公式(1)，其中N为芯片中具有Pathway注释的基因数目；n为N中差异表达基因的数目；M为芯片中注释为某特定Pathway的基因数目；m为注释为某特定Pathway的差异表达基因数目；P＜0.01的Pathway定义为在差异表达基因中显著富集的Pathway；通过Pathway显著性富集能确定差异表达基因参与的最主要生化代谢途径和信号转导途径。

本发明的有益效果是：

采用高通量转录组测序技术，构建不同品种猪背最长肌组织的基因表达谱，获得差异表达基因，利用GO分析和pathway分析对差异基因进行筛选，确定调控肉质性状的关键基因，加深对肌肉生长发育和脂肪代谢调控的了解，为之后研究形成国内外品种在肉质和生长性状方面的差异的分子机理提供分子生物学证据。

附图说明

下面结合附图和具体实施方式对本发明作进一步详细的说明。

图1是本发明实施例的皖南花猪和大约克猪差异表达基因的GO功能富集分布图。

具体实施方式

1材料与方法

1.1试验材料与测序

实验所需组织样品来自两个不同猪种分别为皖南花猪和大约克猪，采取成年期的背最长肌组织，每个猪种采集三头个体，样品采集后-80℃贮存，提取各样品总RNA，送交上海伯豪生物技术有限公司采用Solexa测序技术对样品进行测序。

1.2数据分析

1.2.1原始数据处理。

测序仪产生的原始图像数据经Base Calling转化为序列数据，我们称之为原始序列数据，结果以FASTQ文件格式存储。由于原始序列数据可能包含低质量序列、接头序列等杂质序列数据，不能直接用于生物信息分析，所以原始序列数据必需经过数据处理转换为高质量序列数据，利用Tophat(version：2.0.6)软件的spliced mapping算法对高质量序列数据进行猪基因组比对，采用基因组版本为Sscrofa10.2。

1.2.2差异基因的筛选。

基因表达量的计算使FPKM法。根据基因的表达量(FPKM值)计算该基因在不同样本间的差异表达倍数。差异表达基因定义为FDR≤0.01且倍数差异在2倍及以上的基因。

1.2.3GO功能显著性富集分析

把所有差异表达基因向Gene Ontology数据库的各条目映射，计算每个条目的基因数目，然后应用超几何检验，找出与整个基因组背景相比，在差异表达基因中显著富集的GO条目，其计算公式为：

$P=1-\underset{i=0}{\overset{m-1}{\mathrm{\Σ}}}\frac{\left(\begin{array}{c}M\\ i\end{array}\right)\left(\begin{array}{c}N-M\\ n-i\end{array}\right)}{\left(\begin{array}{c}N\\ n\end{array}\right)}$

其中，N为基因组中具有GO注释的基因数目；n为N中差异表达基因的数目；M为基因组中注释为某特定GO条目的基因数目；m为注释为某特定GO条目的差异表达基因数目。计算得到的P值通过Bonferroni校正之后，以校正后的P值＜0.01为阈值，满足此条件的GO条目定义为在差异表达基因中显著富集的GO条目。通过GO功能显著性富集分析能确定差异表达基因行使的主要生物学功能。

1.2.4Pathway显著性富集分析

Pathway显著性富集分析以KEGG Pathway为单位，应用超几何检验，找出与整个基因组背景相比，在差异表达基因中显著性富集的Pathway。该分析的计算公式同GO功能显著性富集分析，在这里N为芯片中具有Pathway注释的基因数目；n为N中差异表达基因的数目；M为芯片中注释为某特定Pathway的基因数目；m为注释为某特定Pathway的差异表达基因数目。P＜0.01的Pathway定义为在差异表达基因中显著富集的Pathway。通过Pathway显著性富集能确定差异表达基因参与的最主要生化代谢途径和信号转导途径。

2结果与分析

2.1皖南花猪和大约克猪背最长肌的肉质性状

表1皖南花猪和大约克猪背最长肌的肉质性状

从表1可以看出，皖南花猪背最长肌的肉色、PH值、大理石纹和肌内脂肪显著高于大约克猪(P＜0.05)；皖南花猪背最长肌的滴水损失和肌纤维直径显著低于大约克猪(P＜0.05)。

2.2测序质量评估

在表1中，经高通量测序，皖南花猪(WH)每个样本得到的原始序列数据数分别为65，584，126、64，327，738和59，244，502，其中高质量序列数据有效率达到94.4％、94.3％和94.2％。大约克猪(YY)每个样本得到的原始序列数据数分别为57，969，134、68，254，168和72，298，142。其中皖南花猪和大约克猪每个样本的高质量序列数据有效率分别达到94.4％、94.3％、94.2％、93.8％、93.7％和93.8％。可见所测样品中高质量序列数据的比例都占到原始序列数据Rawreads的90％以上，表明测序质量较高，所得高质量序列数据数目能够满足试验需要。

表2皖南花猪和大约克猪每个样品的测序结果

2.3差异表达基因及其功能注释

2.3.1皖南花猪和大约克猪差异基因筛选

基因表达量以RPKM法计算。以FDR小于或等于0.01和差异倍数(Foldchange)绝对值大于或等于2为进行筛选，差异表达基因共有347个，其中以皖南花猪为对照组，大约克猪与之相比，上调的基因有94个，下调的基因有253个。

2.3.2皖南花猪和大约克猪差异表达基因的GO功能富集分析

GO是一个采用动态更新的标准词汇表来描述基因和其产物功能的数据库，目前被广泛应用于生物的转录组数据分析研究中，GO总共分为3大功能类，分别描述基因的分子功能、所处的细胞位置和参与的生物过程。背最长肌差异表达基因中已注释的基因有6937个，GO功能富集得到分类条目有487个。其中涉及细胞组分、分子功能以及生物学过程的基因所占比例分别为12.7％、23.9％和63.4％，可以推测，绝大部分差异表达基因在一些生物学过程扮演重要角色。差异表达基因参与的显著性生物过程(P＜0.01)主要包括糖酵解、肌肉纤维的发育、心肌收缩、横纹肌收缩和骨骼肌收缩等；细胞所处位置主要表现在肌原纤维相邻两肌节交接处的盘状(Z盘)结构、肌钙蛋白复合、肌浆网和线粒体基质等；其分子功能主要表现在肌肉的结构组成、钙离子通道调节活性、腺苷酸激酶活性和AMP活性等(图1)。

2.3.3皖南花猪和大约克猪差异表达基因的KEGG富集分析

KEGG数据库记录细胞中基因产物的功能以及基因产物的相互作用网络，基于KEGG通路的分析有助于我们进一步了解基因的生物学功能。对皖南花猪和大约克猪背最长肌间差异表达基因进行KEGG pathway信号通路分析，结果表明KEGG数据库注释的差异表达基因有6620个，Pathway显著性富集分析发现这些基因参与了146条生物学代谢通路，其中有22条存在显著性富集(P＜0.01)。其中参与脂肪代谢途径的差异基因数目最多有37个，主要包括脂肪酸代谢、PPAR信号通路、丙酸代谢、脂肪酸生物合成、脂肪细胞因子信号转导通路、脂肪酸伸长率和不饱和脂肪酸生物合成，其次是代谢途径、碳代谢、糖酵解和丙酮酸代谢等。说明脂肪代谢相关基因的差异表达可能是皖南花猪和大约克猪肉质差异的重要原因。(表3)

表3皖南花猪和大约克猪背最长肌差异基因显著富集pathway

3、讨论

本试验选用在生长速度和肉质性状等方面存在较大差异的皖南花猪和大约克猪的背最长肌进行高通量转录组测序，筛选与肉质性状相关的基因。结合GO功能富集分析和pathway通路富集分析对高通量测序技术获得的皖南花猪与大约克猪背最长肌间差异表达基因数据进行了生物信息学分析。结果发现能量代谢中的糖酵解代谢、脂肪沉积代谢和肌肉发育等显著通路和生物过程，其可能与皖南花猪和大约克猪背最长肌肉质性状的差异有关。

糖酵解代谢是影响肉质性状形成的重要生理过程。糖酵解潜能是指动物在屠宰前肌肉含有的发生糖酵解产生乳酸的底物(如肌糖原、葡萄糖-6-磷酸和葡萄糖)以及乳酸的总量。乳酸是影响肉质性状的重要生化指标，对pH和肉色都有显著影响。pH是重要的肉质性状，肉色、保水力、熟肉率等诸多性状都直接或间接地决定于pH的高低。终端pH越低、糖酵解潜力越高，则肉色越苍白，保水力越差。本试验发现。GO功能分析的最显著生物过程均与糖酵解代谢有关，参与这些的差异基因包括GPI、GAPDH、LDHA、PGK1、ENO3、PFKM、TPI1、PGAM2，且这8个基因在大约克猪背最长肌中的表达量显著高于皖南花猪，说明皖南花猪的PH值显著高于大约克猪，这与表1中皖南花猪的PH值显著高于大约克猪相一致。此外，段艳宇等研究发现猪背最长肌GP与pH24h值呈显著的负相关(r＝-0.39，P＜0.001)、与b值和L值呈显著正相关(b值：r＝0.36，P＜0.001、L值：r＝0.13，P＜0.001)，结果表明背最长肌GP越高，则肉的终pH值趋向更低，肉色也更加苍白。PGK是真核细胞糖酵解过程中一种主要的组成酶，它在将葡萄糖转变丙酮酸的过程中将1，3-二磷酸甘油酸分子C1上的高能磷酸键转移给ADP生成ATP。徐德全研究表明PGK1在IMF含量低的大白猪中高表达，梅山猪、通城猪中低表达，可能是影响IMF沉积的下调基因，这与本试验中PGK1在大约克猪中高表达，在皖南花猪中低表达的结果相一致。TPI1基因则主要参与糖异生、糖酵解和脂肪酸生物合成。Laville等通过对猪背最长肌肌浆蛋白质的双向凝胶电泳发现嫩度较高肉样中的TPI1含量显著高于嫩度较低的肉样(P＜0.05)。

肌纤维类型及组成被认为是调控肉质的关键。肌纤维类型也显著影响肌肉的生长。猪出生后骨骼肌中主要含4种肌纤维，根据肌球蛋白重链的不同亚型将其分为I型、IIA型、IIB型和IIX型。本试验发现GO功能分析的显著性生物过程均与肌肉纤维发育有关，参与此过程的差异基因包括MLC2V、CHRNB1，且这些基因在皖南花猪背最长肌中的表达量显著高于大约克猪。这说明皖南花猪与大约克猪背最长肌的肌纤维组成存在显著差异，可能与两品种猪肉质差异和生长速度差异有关，这与表1中的大约克猪的肌纤维直径显著大于皖南花猪的结果相一致。同时也与以前的报道基本一致，杨晓静等通过比较肉质欠佳的瘦肉型猪大白猪和肉质鲜美的肥胖型猪二花脸猪背最长肌中I、IIA、IIX和IIB四种肌纤维类型组成差异及其发育性变化。结果表明，二花脸猪背最长肌含有I和IIA型肌纤维显著高于大白猪，而大白猪含有IIB型肌纤维显著高于二花脸猪，这与二花脸猪优良肉质有关，此外酵解型肌纤维比例的增加必将导致生长速度加快；郭佳等分析了金华猪和长白猪背最长肌肌纤维组成差异，结果表明金华猪优良的肉质可能与其背最长肌中丰富I和IIA型肌纤维含量有关，而长白猪背最长肌中更多的IIB型肌纤维含量可能是其生长快而肉质欠佳的主要原因；以及中国地方猪如莱芜猪、梅山猪与国外猪种肌纤维组成差异，也得到类似结果。在显著性生物过程中还有一些与肌肉组织发育有关的基因，包括CHRNB1、ACTN2、TCAP、CSRP3、TNNT1、PGAM2、ARG2。PGAM2基因编码的蛋白参与糖酵解过程，与横纹肌收缩有关，Tang等研究表明PGAM2与骨骼肌的生长发育确实相关。TCAP通过与筒箭毒碱的相互作用，影响生肌细胞前体分泌筒箭毒碱，控制骨骼肌的生长。CSRP3又叫肌肉特异性LIM蛋白(muscle LIM protein，MLP)是骨骼肌发育的一种正性调节因子，邱海芳等研究表明CSRP3与肌肉发育有关。

脂肪性状影响哺乳动物的肉质，由于脂肪水平受品种和年龄的影响，因此哺乳动物体脂的调控是一个复杂的过程。猪脂肪组织的沉积还具有品种差异，不同品种的脂肪沉积能力相差较大，表现为脂肪型猪种的脂肪沉积能力强于瘦肉型猪种。在pathway差异显著的通路中，发现了一些与脂肪代谢显著相关的信号通路，分别是脂肪酸代谢、丙酸代谢、PPAR信号通路、脂肪细胞因子信号通路、脂肪酸伸长率、脂肪酸生物合成和不饱和脂肪酸生物合成，参与这些通路的差异基因包括UBC、RXRG、LDHA、HADHA、ACADM、ACADL、CPT1B、CPT2、ACSL1、ACSL3、ACAA2、ECI1、FABP3、CD36、ACSS1、ACACB、ACSS3、PRKAG1、TECR。说明皖南花猪和大约克猪的背最长肌在脂类代谢方面存在显著差异。此外，李明洲等研究发现脂肪型猪种太湖猪在不同月龄时脂肪细胞体积、IMF含量等指标均高于瘦肉型猪种长白猪，且在3月龄后差距逐渐明显。杨帆通过对猪生长、胴体及肉质与候选基因集的关联分析发现丙酸代谢通路与脂肪酸代谢通路内基因间的互作效应显著影响猪生长、胴体及肉质性状。在筛选到的差异表达基因中，FABP3、ACSL1和ACACB是参与脂类代谢的基因且在皖南花猪背最长肌中的表达量都显著高于大约克猪。FABP3又称H-FABP主要功能是参与细胞内的脂肪酸运输，可将脂肪酸从细胞膜上运到脂肪酸氧化和甘油三酯及磷酸的合成位置。有研究表明，FABP3与畜禽肌内脂肪含量呈正相关，可影响肉质的嫩度、风味和多汁性，故将FABP3基因作为改善肉质的候选基因之一，这与表1中皖南花猪背最长肌的肌内脂肪显著高于大约克猪相一致。ACSL1和ACACB基因参与胰岛素抵抗和肥胖等代谢过程，在牛的肌肉组织中，ACSL1基因影响多不饱和脂肪酸的沉积，抑制ACACB的表达可以促进肌肉组织的脂肪酸氧化，改善全身葡萄糖动态平衡。在这些差异基因中目前已有研究证明其对于猪的肉质性状有重要影响的基因主要有FABP3、CPT1B、CD36、CSRP3、ACSS1、PGAM2、PGK1和PFKM等。还有一些基因尚未见相关报道如TECR、GPI和MLC2V，但是其已知功能均与调控细胞生长和脂肪酸生物合成及氧化相关，可作为肉质性状候选基因进行深入研究和验证。

猪是人类较为理想的医学模型，通过建立与人类疾病相似的猪医学模型，研究疾病发生的基因调控或作用机制，可为人类疾病防治提供科学依据。本试验中发现了7个在肥厚型心肌病和6个在扩张型心肌病通路中差异表达的基因，其中包括MYH7、MYBPC3、TNNT2和TNNI3等已知的与肥厚型心肌病相关的基因以及已知的LMNA、SCN5A、CTLA-4等基因与扩张型心肌病有关，这也预示着我们可以考虑通过构建肥厚型心肌病病猪和扩张型心肌病病猪模型，来研究这两种病的致病机理以及开发临床用药等。

4、结论

本试验利用Solexa测序技术对皖南花猪和大约克猪背最长肌进行转录组测序。通过比对分析，筛选出了347个差异表达基因，其中上调的基因有94个，下调的基因有253个。对表达差异基因进行GO功能分析、KEGG通路分析筛选出了部分与肉质性状相关的候选基因，如TECR、GPI和MLC2V，在GO和pathway中发现了与肉质显著相关的糖酵解代谢、脂类代谢和肌肉发育相关生物过程和代谢通路，表明肉质性状不仅受到遗传因素的影响还受到宰后糖酵解代谢的影响，初步揭示了皖南花猪和大约克猪之间肉质性状差异的原因，为挖掘新基因和深入研究猪肉质性状遗传因素等奠定了基础。

以上所述的本发明实施方式，并不构成对本发明保护范围的限定。任何在本发明的精神和原则之内所作的修改、等同替换和改进等，均应包含在本发明的权利要求保护范围之内。

Claims (1)

1.不同品种猪背最长肌差异表达基因的筛选与注释方法，包括以下步骤：

(1)采样与测序：按不同猪种，采取成年期的背最长肌组织，每个猪种采集三头个体，样品采集后-80℃贮存，提取各样品总RNA，采用Solexa测序技术对样品进行测序；

(2)原始数据处理：将测序仪产生的原始图像数据经Base Calling转化为序列数据，并称之为原始序列数据，其结果以FASTQ文件格式存储；将原始序列数据经过数据处理转换为高质量序列数据，利用Tophat软件的splicedmapping算法对高质量序列数据进行猪基因组比对，采用基因组版本为Sscrofa10.2；

(3)差异基因的筛选：基因表达量的计算使FPKM法；根据基因的表达量(FPKM值)计算该基因在不同样本间的差异表达倍数；差异表达基因定义为FDR≤0.01且倍数差异在2倍及以上的基因；

(4)功能显著性富集分析：把所有差异表达基因向Gene Onto logy数据库的各条目映射，计算每个条目的基因数目，然后应用超几何检验，找出与整个基因组背景相比，在差异表达基因中显著富集的GO条目，其计算公式为：

$P=1-\underset{i=0}{\overset{m-1}{\mathrm{\Σ}}}\frac{\left(\underset{i}{M}\right)\left(\underset{n-i}{N-M}\right)}{\left(\underset{n}{N}\right)}---\left(1\right)$

公式(1)中，N为基因组中具有GO注释的基因数目；n为N中差异表达基因的数目；M为基因组中注释为某特定GO条目的基因数目；m为注释为某特定GO条目的差异表达基因数目；计算得到的P值通过Bonferroni校正之后，以校正后的P值＜0.01为阈值，满足此条件的GO条目定义为在差异表达基因中显著富集的GO条目；

(5)Pathway显著性富集分析：Pathway显著性富集分析以KEGG Pathway为单位，应用超几何检验，找出与整个基因组背景相比，在差异表达基因中显著性富集的Pathway；该分析的计算公式同GO功能显著性富集分析所用公式(1)，其中N为芯片中具有Pathway注释的基因数目；n为N中差异表达基因的数目；M为芯片中注释为某特定Pathway的基因数目；m为注释为某特定Pathway的差异表达基因数目；P＜0.01的Pathway定义为在差异表达基因中显著富集的Pathway。