基于转录组挖掘不同性别贵州白山羊肌肉生长发育的关键基因及PI3K/ATK信号通路分析

安清明，王 星，吴震洋，孟金柱，赵园园，宋兴超

(贵州省梵净山地区生物多样性保护与利用重点实验室，铜仁学院 农林工程与规划学院，贵州 铜仁 554300)

当下社会，人们饮食习惯不断发生改变，对健康饮食和养生保健的关注度越来越高，意味着对肉质食品的质量要求变得更高，这也为我国畜禽肉质性状的遗传改良提出了新的要求。贵州白山羊是我国云贵高原特有的地方山羊(Caprahircus)品种，是发展当地畜牧业的重要种质资源。研究贵州白山羊优良胴体肌肉性状的分子调控机理，挖掘、鉴定控制白山羊胴体肌肉生长发育的主效基因，并应用现代分子标记育种与传统育种法相结合的新途径，对改良贵州白山羊胴体肌肉生长性状具有重要意义。转录组测序技术(Transcriptome Sequencing，RNA-Seq)是基于Illumina Hiseq测序平台能够快速检测特定物种在全基因组水平上基因表达和调控的重要方法，诸如新基因预测、基因结构和功能优化、可变剪接等，目前广泛应用于生物科学领域，尤其在畜禽重要性状分子机制研究中已被广泛应用[1-2]。

动物肌肉组织的生长发育状况直接决定了动物的经济效益。肌肉组织的生长发育指标又可以从生长速度和肉品质两方面评价，肌肉的生长速度直接决定动物胴体质量和产肉量，而肉品质通常用嫩度、吸水力、风味、营养成分和肌内脂肪含量等指标进行评定，直接影响动物产出肉的评级和价格。而大多数决定肌肉生长发育性状的因素都具有遗传效应，受相应基因调控，如骨骼肌生长发育、肌内脂肪含量、嫩度等指标[3-6]。已有研究报道，孟祥忍等[7]对1周龄不同肉牛品种的背最长肌进行差异基因筛选，最终筛选到2个可能与肌肉嫩度相关的功能基因。赵伟明等[8]对不同肉牛背膘脂肪进行转录组测序分析，最终筛选出1 296个差异表达基因。孟宪然等[9]对不同性别、年龄的绒山羊背最长肌进行差异基因筛选和验证，最终筛选出3个可能与肉品质相关的候选基因。王位等[6]对绵羊背最长肌进行转录组分析，最终筛选出6个与肉质性状相关的候选基因。刘远等[10]对不同山羊背最长肌进行比较转录组分析，挖掘出707个可能与骨骼肌生长发育和肌内脂肪含量相关的新基因，这些研究均为家畜的胴体肌肉生长发育遗传机理提供了一定研究基础。

目前，虽然对家畜动物胴体肌肉生长相关基因进行了一定筛选，但是有关地方山羊品种胴体肌肉生长发育的研究仍比较匮乏，尤其对不同性别山羊胴体肌肉生长发育的研究更为稀少。贵州白山羊作为云贵高原特有品种，是发展贵州乃至西南地区畜牧业的重要种质资源，值得对其优良胴体肌肉生长发育的遗传特性进行深入挖掘研究。

本试验选择同等饲养水平下、不同性别的贵州白山羊作为研究对象，采用RNA-Seq技术对二者背最长肌的表达基因进行筛选分析，鉴定关键基因和筛选重要信号通路，以期为贵州白山羊种质资源的开发利用提供理论依据。

1 材料和方法

1.1 试验动物及样品采集

试验所需贵州白山羊由贵州省铜仁市沿河土家族自治县中寨镇永珍牧业养殖合作社提供。所选山羊在同等水平下采用半封闭式饲养，选取2周龄健康的公、母羊各3只，按照国家实验动物屠宰标准进行屠宰并测定屠宰性能指标(宰前活质量、胴体质量、屠宰率、净肉重和净肉率)，采集6只白山羊(公、母各3只)第10～12根胸椎位置背最长肌装入冻存管并移至液氮中进行冷冻保存，待运回实验室后置于-80 ℃冰箱中保存，用于转录组测序和RT-qPCR。

1.2 总RNA提取、文库构建及测序

将-80 ℃冰箱中的背最长肌置于冰盒中解冻，利用RNAprep pure Tissue Kit (TianGen，北京天根生化科技有限公司)提取6个样品背最长肌组织的总RNA，公羊样本标注为R1、R2和R3，母羊样本标注为E1、E2和E3。为保证数据完整和准确性，采用Agilent Bioanalyzer 2100检测完整性，Qubit 2.0 Flurometer测量RNA浓度，RNA样品检测合格后送北京诺禾致源生物信息科技有限公司进行文库构建及测序，测序主要通过Illumina Hiseq 2500平台完成。

1.3 数据处理及分析

利用FastQC软件(http：//www.bioinformatics.babraha m.ac.uk/projects/fastqc/)对测序得到的原始数据进行质量评估，去除原始数据中含有带接头和低质量的片段，保证获得的数据具有品质较高的有效测序数据(Clean reads)。利用TopHat 2软件将获得的高质量序列数据与已知山羊参考基因组进行比对，获取测序拼装序列在参考基因组上的相关位置信息以及样品的序列特征信息。利用Cufflinks软件，对转录本数据进行组装、注释和合并，用DESeq 2软件分析mRNA的差异表达数据，利用基因本体(Gene ontology，GO)数据库对筛选出的表达差异基因进行GO分析，利用KEGG数据库进行DEGs的信号通路分析。

1.4 荧光定量PCR验证

为验证测序结果的可靠性，对差异表达的mRNA进行验证。挑选6个差异表达的候选基因，包括3个上调基因和3个下调基因，进行RT-qPCR验证。利用反转录试剂盒(购自北京全式金生物技术有限公司)将高通量测序使用剩余样品的总RNA进行反转录产生所需cDNA。应用Primer 5.0软件在线设计所挑选基因的引物，使用甘油醛-3-磷酸脱氢酶(Glyceraldehyde-3-phosphate dehydrogenase，GAPDH)基因作参考基因，引物由上海生工合成(表1)。

表1 荧光定量引物基因Tab.1 The list of gene primers of RT-qPCR

RT-qPCR验证不同性别贵州白山羊背最长肌中差异表达mRNA的相对表达水平时采用3个样本重复，采用TransStart®Tip Green qPCR SuperMix(北京全式金生物技术有限公司)进行相对定量扩增分析，扩增体系20 μL：cDNA模板1 μL，2×Transstart®Tip Green qPCR super mix 10 μL，上、下游引物各0.5 μL，RNA-free H2O 8 μL。RT-qPCR反应程序为：94 ℃预变性1 min；94 ℃ 10 s，59 ℃ 30 s，72 ℃ 10 s，46个循环。使用2-ΔΔCt法计算各基因间的相对表达水平。

2 结果与分析

2.1 贵州白山羊屠宰性能差异分析

测定不同性别白山羊2周龄屠宰性能，公羊群体的宰前活质量、胴体质量和净肉质量显著高于母羊，屠宰率和净肉率差异虽然不显著，但公羊仍要高于母羊，该结果表明，公羊的屠宰性能要优于母羊(表2)。

表2 不同性别贵州白山羊屠宰性能差异分析Tab.2 Slaughter performance differences between ram and ewe of Guizhou wihte goat

2.2 RNA-Seq数据比对分析

测序文库经过质量检测，各样本均符合要求，最终6个样本经高通量测序共获得78.99 Gb Clean bases Data，单个样本Clean bases Data均达到12.45 Gb以上，获得Clean reads Date在83 030 104～95 739 024条。结果中碱基质量Q30以上的百分比为89.86%～92.86%，碱基测序错误率均为2%，表明测序质量高、测序数据准确可靠，可用于后续分析试验。

将获得的Clean reads与参考山羊基因组序列进行比对分析，结果表明，各样品与参考基因组的比对效率为93.75%～94.79%，在参考序列上有唯一比对位置的reads数量占Clean reads的84.14%～87.21%，有多个比对位置reads占6.55%～10.65%(表3)。

表3 RNA-Seq与参考山羊基因组的序列比对Tab.3 RNA-Seq and mapping to the reference genome

2.3 差异表达基因的筛选

以FPKM值作为衡量转录本和基因表达差异水平的指标，分析得到6个样本共25 089个表达基因的FPKM值，表4列出在公羊中表达量较高的前20个基因。同时应用DESeq 2软件设定参数≥FPKM≥1，MF-FPKM/FL-FPKM > 1，P<0.05，对获得25 089个mRNA进行差异表达分析，共获得1 077个差异表达的mRNA，其中有194个为新发现基因(转录本)；其中公羊背最长肌中表达显著上调的基因为563个，显著下调的为514个。聚类分析热图显示，公羊和母羊样本间差异表达基因重复性较好，差异基因分别聚类在一起(图1)。

表4 公羊和母羊背最长肌中表达量最高的20个基因Tab.4 Top 20 highly expressed genes in ram and ewe longissimus doris tissue

A图中红色代表显著上调基因，绿色代表显著下调基因；B图中红色代表表达量上调基因，蓝色代表表达量下调基因。Red represents significantly up-regulated gene and green represents significantly down-regulated gene in Fig.A；Red signal represents genes with up-regulated expression，blue signal represents genes with down-regulated expression gene in Fig.B.

2.4 GO功能富集分析

GO注释系统可分为三大类，分别为：生物学过程(Biological process，BP)、细胞组分(Cellular component，CC)和分子功能(Molecular function，MF)。通过Goseq软件对筛选出的1 077个差异表达基因进行GO功能富集分析，结果显示，其中587个基因为显著富集(Q-value≥0.05)，这些基因在3个分支中均有分布，共有35个功能亚类包括14个生物学过程(BP)、20个细胞组分(CC)、1个分子功能(MF)(图2)。

1.RNA加工；2.ncRNA加工；3.核糖核蛋白复合体生物合成；4.mRNA加工；5.通过酯交换反应的RNA剪接；6.以膨大的腺苷为亲核试剂，通过酯交换反应的RNA剪接；7.通过剪接体的mRNA剪接；8.ncRNA代谢过程；9.RNA剪接；10.信使RNA代谢；11.核糖核蛋白复合体亚基组织；12.核糖核蛋白复合物组装；13.核酸代谢过程；14.核糖体生物合成；15.细胞核部分；16.核腔；17.细胞核；18.细胞内细胞器内腔；19.细胞器内腔；20.细胞内膜结合细胞器；21.细胞膜内腔；22.膜结合细胞器；23.细胞内细胞器部分；24.细胞内的细胞器；25.细胞器；26.核仁；27.细胞器部分；28.细胞内；29.细胞内的部分；30.核糖核蛋白；31.大分子复合体；32.无膜细胞器；33.胞内无膜细胞器；34.核质；35.RNA结合。

2.5 差异表达基因的KEGG信号通路分析

通过Kobas软件对差异表达基因进行分析，FDR≤0.05的通路为表达基因中显著富集的信号通路，结果显示，可被注释的表达基因共参与243个通路，其中参与磷脂酰肌醇3激酶(PI3K-AKT)信号通路相关的基因最为富集，共注释到32个基因，图3为差异表达基因的KEGG通路富集结果中显著性Q值最小的20个通路。

图3 背最长肌中差异表达基因KEGG通路富集散点图Fig.3 Scatter plot of enriched KEGG pathway in longissimus doris tissue

在被注释到PI3K/AKT信号通路中的32个基因中，ITGAV、LAMA4、LAMB1等17个基因为显著上调基因，CBL基因为显著下调基因(表5)。对通路中的32个基因进行共表达分析，发现COL4A1和COL4A2基因共表达估值最高，且在人、普通牛和小鼠中的共表达估值最高，ITGAV基因共表达最为丰富，与12个蛋白共表达。对通路中注释到的32个基因进行蛋白互作网络分析，结果表明，CSF3基因编码蛋白互作效应最强，ITGAV基因编码蛋白网络互作数最多，LPAR6基因编码蛋白网络互作数最少(图4)。

表5 参与磷脂酰肌醇3/激酶信号通路的基因Tab.5 Genes involved in PI3K/AKT signaling pathway

图4 PI3K/AKT通路注释基因共表达分析(A)蛋白网络互作分析(B)Fig.4 Coexpression analysis of annotated genes (A)and protein-protein network interaction(B) in PI3K/AKT signaling pathway

2.6 RT-qPCR验证分析

为了验证转录组测序筛选差异表达基因的准确性，本研究从差异表达基因中筛选出6个候选基因进行验证(表6)，其中FHL3、WFIKKN2、SOX6基因在公羊中为上调基因，QSOX2、MYH2、LAP基因为下调基因，RT-qPCR结果显示，这6个候选基因的mRNA表达趋势与转录组测序结果一致，且FHL3、WFIKKN2、SOX6基因在公羊背最长肌中的表达量均显著高于母羊(P<0.05)，QSOX2、MYH2、LAP基因在公羊背最长肌中的表达量均显著低于母羊(图5)(P<0.05)。

表6 可能与不同性别山羊背最长肌发育相关的候选基因Tab.6 Candidate genes may associate with longissimus doris tissue development in ram and ewe

不同小写字母表示差异显著(P <0.05)。Different lowecase letters show significant difference(P <0.05).

3 结论与讨论

本试验选用同等饲养条件下、同年龄段、不同体质量的贵州白山羊公、母羊背最长肌为研究对象，对其进行转录组测序，筛选与胴体肌肉生长发育相关的基因，结合相关生物学功能分析，共发现了1 077个差异表达基因，其中有194个为新发现基因或转录本；其中563个基因在公羊背最长肌中表达显著上调，514个表达显著下调，这些差异表达基因可能与贵州白山羊的背最长肌生长发育相关，值得进一步深入挖掘研究。

对差异基因进行GO功能注释和KEEG通路富集分析，发现PI3K/AKT信号通路是一条重要的信号通路，该通路在细胞增殖、生长、代谢和生存中起着重要调控作用[11]，通路中磷脂酰肌醇3激酶(PI3K)主要分为三类，其中Ⅰ类最为重要，其主要与生长因子介导的PI3K信号通路密切相关，其中p85α亚基可以介导胰岛素的调控进而影响葡萄糖转运[12]，还可以通过与p110α亚基结合影响p110α亚基的活性[13]。研究还发现，PTEN基因编码蛋白与PI3K的p85亚基之间的相互作用对PI3K/AKT信号通路的稳态起着重要作用[14]，本研究在PI3K/AKT信号通路中筛选到PTEN基因，也筛选到了编码2种亚基的基因PIK3CA和PIK3R1。在该通路中蛋白激酶B(AKT)是另一个关键因子，其重要的亚型AKT2主要在骨骼肌、棕色脂肪中表达，而且研究发现，活化后的AKT可以激活上百种底物，达到调控细胞生长、增殖、代谢、运动等功能[15-16]。吕欣等[17]分析也认为，通过影响PI3K/AKT信号通路中的蛋白表达、诱导通路的磷酸化，可以促进肌卫星细胞增殖分化、抑制凋亡，从而影响骨骼肌的再生修复。而本研究对不同山羊背最长肌生长的转录组分析筛选出PI3K/AKT信号通路是最为富集的信号通路，这也证明该通路是调控胴体肌肉生长发育的重要信号通路，值得进一步深入研究。

在PI3K/AKT信号通路筛选出的32个注释基因中，ITGAV基因编码蛋白(Integrin subunit alpha V，整合素αV亚基)在蛋白互作网络中互作蛋白数最为丰富，该基因属于整合素α链家族成员，研究表明，编码蛋白可通过精-甘-天冬序列(Arg-Gly-Asp，R-G-D)与多种细胞外基质配体相结合，从而介导细胞与胞外基质的信号传导，在调控肿瘤血管生成、肿瘤的迁移和侵袭等过程中发挥了重要作用[18]。马春辉[19]通过转录组分析发现，ITGAV基因可能在大鼠副侧韧带损伤恢复过程中发挥重要作用。本研究在PI3K/AKT信号通路筛选出的注释基因中，包括了COL家族的4个基因，分别为COL4A1、COL4A2、COL5A3和COL6A6，其中COL4A1和COL4A2基因编码蛋白共表达估值最高。已有研究表明，COL4A1、COL4A2基因分别编码Ⅳ型胶原蛋白的α1链和α2链，主要在细胞基底膜上表达，与细胞网状支架形成相关[20]，综上，本研究在通路中筛选出的注释基因可能都与细胞的生长发育相关，可进一步作为山羊胴体肌肉生长发育的候选基因。

通过转录组测序筛选出的1 077个差异表达基因中，筛选出6个可能与贵州白山羊背最长肌生长发育相关差异表达候选基因(上调基因3个，下调基因3个)，分别为FHL3、WFIKKN2、SOX6、QSOX2、MYH2和LAP基因，进行实时荧光定量PCR验证，验证结果均与RNA-Seq分析结果一致。在上调基因中，FHL3基因的cDNA全长843 bp，编码281个氨基酸，转录蛋白属于LIM蛋白家族[21]。FHL3基因主要在骨骼肌中高表达，FHL3蛋白可以通过和肌动蛋白相结合对肌动蛋白成束起到负调控的作用[22]，同时发现，FHL3基因在成肌细胞中能够诱导纤维解聚，说明FHL3是激动蛋白细胞骨架动态中的重要调控因子[22]。在畜牧生产中，有研究发现，FHL3基因的突变与猪肉的眼肌面积、肉质性状、肌内脂肪、肌肉嫩度等肉质生产性状相关[23-24]。WFIKKN2基因编码蛋白是一种分泌蛋白，可抑制蛋白酶活性。研究发现，WFIKKN2基因编码蛋白可以通过结合生长分化因子GDF11而参与调控骨骼肌的生长发育过程[25]。在畜牧业生产中，孙玲[26]发现，WFIKKN2基因与MSTN基因随骨骼肌的发育变化而表现出相反的表达规律。Wang等[27]研究发现，WFIKKN2基因变异与绵羊的肩部瘦肉量相关。SOX6基因最早从成年小鼠睾丸中分离得到[28]，其编码蛋白属于一种转录因子，后续发现其在脊椎动物发育中起着重要调控作用[29]，如SOX6基因编码蛋白参与维持肌肉组织的正常生理功能[30]，参与软骨和骨骼肌的生成等[31]。同时还发现，SOX6基因在细胞增殖和分化中起重要作用，如Jackson等[32]发现，Sox6是调节哺乳动物骨骼肌纤维类型分化的关键因素之一，Quiat等[33]发现，在小鼠骼肌中Sox6蛋白的缺失会导致慢肌纤维显著增加，肌肉颜色变鲜红。

在下调基因中，QSOX2基因相关研究仍较少，目前还未发现该基因与哺乳动物肌肉生长或个体发育相关的研究报道，仅Okada等[34]研究发现，QSOX2基因rs12338076突变可能与日本人群个体身高相关。在畜牧业生产中，骨架(体高)比较大的动物具有更为良好的育肥前景，且畜禽的日增质量直接受肌肉发育的影响。MYH2基因属于MYH基因家族，该家族基因能够编码一系列与肌肉生长发育相关的MyHC蛋白亚型，其中MYH2基因编码的亚型蛋白主要存在于快速氧化型肌纤维中，且在骨骼肌中也有发现，研究表明，MYH2基因是一种肌节调控基因，主要与动物的肌肉发育相关[35-37]，而肌纤维是影响肉品质的一个重要因素；李娜[38]通过GWAS研究发现，MYH2基因可作为影响猪肉中肌纤维性状和肌内脂肪的重要候选基因。LAP基因编码蛋白属于一类对亮氨酸残基具有最高切割效率的氨肽酶，这类蛋白在调节细胞极性和细胞增殖中具有重要作用[39]，研究发现，LAP基因可以通过影响泌乳细胞的功能调节奶牛的产奶量、牛奶的乳脂率和乳蛋白含量[40]，亦有研究通过全基因组关联分析发现，LAP基因与绵羊的体质量显著相关[41]。

综上，本研究共筛选出影响不同性别贵州白山羊背最长肌生长发育和肉品质的差异表达基因1 077个，其中194个为新发掘基因(转录本)；其中在公羊背最长肌中563个基因表达显著上调，514个基因显著下调。通过功能分析及RT-qPCR验证及PI3K/AKT信号通路分析表明，筛选出的基因可作为影响不同性别贵州白山羊背最长肌生长发育和肉品质的重要参考基因，对进一步深入探究白山羊胴体肌肉生长发育的遗传调控机制具有重要意义。