基于RNA-Seq技术挖掘贵州白山羊公羊背最长肌生长发育下调基因

安清明，王 星，吴震洋，孟金柱，赵园园 (．铜仁学院 农林工程与规划学院 贵州省梵净山地区生物多样性保护与利用重点实验室，贵州 铜仁 554300；2．甘肃农业大学 草食动物生物技术重点实验室，甘肃 兰州 730070)

目前虽然已有研究对家畜动物的肌肉生长进行了筛选，获得了部分相关的候选基因，但关于地方山羊品种肌肉生长发育的研究较少，尤其对不同性别的山羊肌肉生长发育还有待深入研究。因此通过研究不同性别贵州白山羊背最长肌组织中差异表达基因，挖掘和鉴定控制白山羊肌肉组织生长发育的相关基因，对改良贵州白山羊胴体肌肉生长发育具有重要意义。本研究采用RNA-Seq技术对不同性别贵州白山羊背最长肌进行分析，筛选与公羊肌肉生长相关的下调基因，并对筛选的下调基因进行GO、KEGG等分析，并通过RT-PCR进行验证，挖掘出可调节贵州白山羊肌肉生长发育的关键候选基因，以期为贵州白山羊在肉质和肌肉生长性状方面的分子机理研究提供理论依据。

1 材料与方法

1.1 样品采集样品采自贵州省铜仁市沿河土家族自治县华珍牧业有限责任公司，选取同等饲养条件下2周岁龄、健康的白山羊6只(公、母羊各3只)，屠宰后采集背最长肌置于液氮中迅速冷冻，带回实验室后置于-80℃保存备用。

1.2 背最长肌RNA提取、文库构建及转录组测序按照TRIzol法提取每只白山羊总RNA，用RNAprep pure Tissue Kit(TianGen，天根生化科技(北京)有限公司)纯化，采用Agilent Bioanalyzer 2100检测完整性，Qubit 2.0 Flurometer测量浓度，待样品合格后送北京诺和致源生物信息科技有限公司构建文库，应用Illumina Hiseq 2500平台进行转录组测序。

1.3 数据处理与分析将获得的原始Reads应用FastQC软件(http://www.bioinformatics.babraha- m.ac.uk/pr-ojects/fastqc/)进行质量评估，去除数据中带有接头、含N和低质量的原始片段，保证获得的数据具有品质较高的Clean Reads。应用Trinity(http://trinityrnaseq.sourceforge.net/)对Clean Reads进行组装，获得unigenes。利用HISAT2软件将获得数据与已知山羊参考基因组进行比对，获得序列在参考基因组上的相关信息。应用DESeq2软件对获得mRNA进行统计分析，筛选差异表达的下调基因mRNA。利用GO seq软件对筛选的下调表达基因进行GO分析，应用KOBAS软件进行KEGG通路分析。

1.4 cDNA及引物合成将转录组测序剩余总RNA反转录为cDNA：按照北京全式金生物技术有限公司反转录试剂盒说明进行，42℃孵化15 min，再85℃加热5 s。从转录组测序筛选出的差异下调表达基因中随机挑选4个基因进行验证(KLF11、TNS1、SIK2和IGF1R)，以GAPDH基因为参考基因，引物由生工生物工程(上海)股份有限公司合成，引物序列见表1。

表1 荧光定量引物信息

1.5 荧光定量PCR验证通过荧光定量验证转录组测序中贵州白山羊公羊(MF)与母羊(FL)背最长肌中差异表达下调mRNA的相对表达量。每个样品设3个重复，应用北京全式金生物技术有限公司试剂盒进行相对定量检测分析，反应体系：cDNA模板1 μL，2×qPCR super mix 10 μL，上、下游引物各0.5 μL，RNA-free H2O补齐至20 μL。反应条件：94℃ 1 min，94℃ 10 s，59℃ 30 s，72℃ 10 s，43个循环。使用2-△△Ct法计算各基因间的相对表达水平。

2 结果

2.1 RNA-seq数据比对分析对采集6个背最长肌样品的转录组数据进行分析，对原始获得的Raw data将带接头、含N和低质量的Reads片段过滤筛选后，最终在MF中获得95739024条Clean reads，结果中碱基质量在Q30以上的占92.86%；在FL中获得89 806 996条Clean reads，结果中碱基质量在Q30以上的占89.86%(图1)。将获得的Clean reads与山羊参考基因组序列进行比对，结果表明MF中比对效率为94.18%，FL比对效率为94.41%，其中在参考序列中能够唯一比对到位置的Reads数量分别占86.24%和86.07%，能够多点定位的Reads分别占7.94%和8.34%，表明测序质量高，数据可靠，可用于后续分析(表2)。

表2 测序数据与参考山羊基因组的序列比对结果

图1 MF和FL背最长肌中原始数据的分类

2.2 差异表达基因的筛选、功能注释和富集分析

2.2.1下调差异表达基因的筛选 以FPKM值作为衡量转录基因下调差异表达水平的指标，分析6个样本共得到33 896个转录本，对MF和FL的FPKM进行标准化后，筛选出FPKM≥1的转录本有25 089个。应用DESeq2软件对获得的25 089个mRNA进行差异表达筛选，设定参数FPKM≥1，且P<0.05时，共获得514差异表达下调基因，其中100个为新发现的基因(转录本)，表3列出MF背最长肌中FPKM数值最高的前20个基因，表4列出MF相较于FL背最长肌中差异表达最高的10个下调基因及其功能。

表3 MF背最长肌中FPKM数值最高的20个基因

表4 可能与贵州白山羊背最长肌发育相关的候选基因

2.2.2差异表达基因的GO功能富集 应用Goseq软件对筛选出的514个差异表达下调基因进行GO功能富集分析，结果共分为2大类21个功能亚类，主要参与生物学过程(biological process，BP)和细胞组分(cellular component，CC)，而分子功能(molecular function，MF)的功能富集在本研究中未检测到，结果见表5。

表5 贵州白山羊MF和FL中上调表达基因GO分析

2.2.3下调差异表达基因KEGG信号通路分析 为了筛选出与贵州白山羊背最长肌生长发育相关的信号通路，应用Kobas软件对下调差异表达基因进行KEGG信号通路分析，FDR≤0.05的通路为差异表达基因中显著富集的信号通路，结果表明被注释的差异基因共参与145个信号通路，其中参与核糖体信号通路的基因最为富集，共有27个被注释基因，图2为下调表达基因中KEGG信号通路结果中显著性最高的20个通路。

图2 贵州白山羊MF和FL背最长肌中下调表达基因KEGG通路富集散点图

富集在核糖体信号通路中的差异蛋白位点发现共有25个(包括27个注释基因)，核糖体信号通路图见图3。25种差异蛋白分别分布在核糖体的60S大亚基和40S小亚基中，其中60S大亚基核糖体蛋白15种(包括17个基因)，40S小亚基核糖体蛋白9种(包括10个基因)，其中S38e小亚基核糖体蛋白未注释到对应基因(表6)。

图3 核糖体信号通路

表6 核糖体信号通路中的差异蛋白位点及基因

2.2.4下调差异表达基因RT-PCR验证分析 通过功能分析，在差异表达的下调基因中，随机选择4个可能与贵州白山羊公羊背最长肌生长发育呈下调趋势的调控基因KLF11、TNS1、SIK2和IGF1R进行验证，具体信息见表7。采用RT-PCR进行验证这4个基因的相对表达量(图4)。

**表示P<0.05图4 候选基因在不同性别白山羊MF和FL背最长肌中的相对表达量

表7 贵州白山羊背最长肌负调控候选基因

3 讨论

肌肉的生长发育受到多种基因和信号通路的调节，同时还受品种、性别、饲养环境、营养等诸多因素的影响，同时肌肉的生长发育直接影响动物的体质量和肉品质。因此，本研究以同等饲养条件下、同年龄段的贵州白山羊公、母羊背最长肌为研究材料，基于RNA-Seq测序技术筛选与肌肉生长发育相关的基因，结合相关生物信息学功能分析，筛选出公羊背最长肌对于母羊背最长肌存在差异的514个下调基因，其中100个为新发现的关联基因，这些差异表达基因可能对不同性别贵州白山羊骨骼肌的生长发育具有一定影响，值得进一步深入研究。

本研究通过功能富集分析发现，差异最为富集的蛋白均属于核糖体信号通路。核糖体是机体细胞中重要的细胞器，是合成蛋白质的重要场所，核糖体中最主要的成分为核糖体蛋白，其广泛分布于各种组织中，在核糖体的翻译和蛋白质的合成中发挥着重要作用。研究还发现，核糖体蛋白同时还具有调控基因转录、调控mRNA的翻译、同时还能影响细胞分化、增殖和凋亡等作用[11]。本研究在核糖体通路中共发现了25种差异表达的蛋白，分布于核糖体的两种不同亚基中，这些差异蛋白是否会影响背最长肌中核糖体的蛋白质生物合成过程，进而影响肌肉的生长发育尚未可知，有待进一步研究。

在筛选出的514个差异下调基因中，随机选择4个可能与贵州白山羊公羊背最长肌肌肉生长发育呈下调趋势的调控基因KLF11、TNS1、SIK2和IGF1R进行RT-PCR验证，结果与RNA-Seq分析结果一致，表明RNA-Seq分析结果可靠，因此筛选出的514个关联基因可用于进一步分析。本研究随机选择的4个下调基因中，KLF11基因编码的转录因子KLF11属于Kruppel转录因子家族(Kruppel like factorsfamily，KLF)成员，该家族转录因子成员众多，目前共发现18个成员，且功能复杂多样，在不同组织和不同生理条件下参与细胞分化、凋亡、衰老、细胞间黏附、糖脂代谢和脂肪细胞分化等生命活动[12-15]。KLF11基因有4个外显子和3个内含子组成，共编码512个氨基酸，编码蛋白相对分子质量为74 kDa。已有研究表明，KLF11基因编码的转录因子可以与不同作用因子结合，发挥不同的转录调控作用。有研究表明，KLF11可以调节TGH-β信号通路的下游信号，从而参与由TGH-β信号通路调控的生理功能，如细胞分化凋亡、抑制细胞增殖等[16-17]。LOMBERK等[18]研究发现KLF11还能参与调控葡萄糖摄取、糖酵解和脂质代谢。王江林等[19]研究发现，通过干扰KLF11基因可促进山羊肌内前体脂肪细胞分化。综上，KLF11基因可能是与白山羊肌肉生长发育的一个重要候选基因。

TNS1基因(tensin 1)编码蛋白称张力蛋白1，属于黏着蛋白家族，该家族蛋白是一种多结构域蛋白，其主要参与细胞膜交互、构建肌动蛋白细胞骨架，因此又被称作支架蛋白，同时还能参与细胞黏着、运动和生长等生物学活动[20]。已有的研究报道主要集中在TNS1基因与癌症的关联分析中，如ZHAN等[21]研究发现TNS1基因高表达状态与乳腺癌细胞的转移呈负相关。曹一鑫等[22]研究发现TNS1基因高表达状态与胃癌肿瘤细胞转移呈负相关。任琎等[23]研究发现TNS1基因可以抑制非小细胞肺癌细胞增殖和侵袭。在家畜动物中，仅靳聪飞等[24]在牛肌肉组织中克隆获得了TNS1基因全部序列，通过生物信息学分析预测TNS1基因编码蛋白不属于分泌蛋白，磷酸化位点均在丝氨酸、苏氨酸和酪氨酸残基上。目前尚未见TNS1基因影响肌肉生长发育的相关报道，因此，TNS1基因是否可以作为影响动物肌肉组织的生长发育的候选基因，仍需进一步深入研究。

SIK2基因(salt-inducible kinase 2)编码蛋白称为盐诱导激酶2，是一种丝氨酸/苏氨酸蛋白激酶，属于AMP活化蛋白激酶(adenosine monophosphate-activated protein kinase， AMPK)亚家族。研究发现SIK2基因在脂肪组织中的表达最高，且在脂肪细胞中过表达SIK2时，可以显著降低细胞核类固醇调节元件结合蛋白1(sterol regulatory element binding protein-1，SREBP-1)水平，进而抑制多个参与脂肪合成代谢调节基因的表达，如乙酰辅酶A羧化酶2基因(ACC2)和脂肪酸合成酶基因(FAS)[25-26]。也有研究表明，SIK2可以通过磷酸化过程参与胰岛素和胰高血糖素对血糖浓度的调节过程中[27]。综上，SIK2基因在脂质代谢和糖代谢中发挥重要的作用，该基因可能是通过参与到调节肌肉中脂肪的生长，进而影响肌肉的嫩度，可以作为调节白山羊肌肉生长发育的候选基因。

胰岛素样生长因子1受体(insulin-like growth factor 1 receptor，IGF1R)是一种表面受体蛋白，研究表明，IGF1R主要通过介导IGF1和IGF2信号参与调节多种组织及器官的生长发育，如IGF1R基因缺失会导致小鼠胚胎死亡或肌肉发育迟缓，而过表达则会导致肌肉组织中肌肉发育肥大[28]。有研究通过抑制IGF1R表达可以减弱肌肉的再生能力及细胞的增殖[29-30]。FÜRSTENBERGER等[31]研究发现IGF1基因具有可以调节细胞的增殖分化和细胞凋亡的作用，而IGF1R是IGF1发挥具体功能的关键受体。而KLOTING等[32]通过降低小鼠脂肪组织中IGF1R的活性可以增加脂肪质量。YANG等[33]研究发现，IGF1R可以通过AKT信号通路促进肌肉的分化。在畜牧生产中，占思远等[34]研究发现IGF1R基因在山羊肌肉组织和肌肉细胞中均为高表达基因，且IGF1R与IGF2R可以协同促进肌肉组织发育和肌细胞增殖。马懿磊等[35]研究表明IGF1R基因在胎牛和成年牛时期的肌肉中的表达存在差异，且干扰IGF1R基因能够抑制牛肌细胞的分化。综上，IGF1R基因具有调节肌肉生长发育的功能，是研究肌肉生长发育的候选基因。

综上，本研究基于RNA-Seq技术筛选出贵州白山羊公羊背最长肌的差异下调基因514个，其中100个为新发掘基因或转录本，发现参与核糖体信号通路的基因最为富集，为27个基因。通过功能分析及RT-PCR验证，初步认为RNA-Seq技术筛选出的相关下调基因数据真实可靠，可作为贵州白山羊公羊背最长肌生长发育相关的重要候选基因，对进一步研究白山羊肌肉生长发育的遗传机制提供了依据。