哈萨克羊骨骼肌差异表达基因的筛选及功能分析

李 雷, 王 悦, 李晓悦, 刘开平, 戴继鸿, 李村院, 胡圣伟

(石河子大学生命科学学院,新疆石河子 832003)

绵羊由于其具有能够为人们提供毛、皮、肉、奶等各类产品的优势,目前已成为我国乃至世界上具有重大经济价值的家畜动物之一[1]。同时,作为畜牧业大省,绵羊是新疆维吾尔自治区畜牧业发展的重要家畜之一,也是当地牧民的主要收入来源[2]。近些年,随着经济的快速发展,人们对含有较高水平必需氨基酸和较低水平胆固醇的羊肉的需求愈来愈大。虽然我国是世界上最大的羊肉生产和消费大国,但与肉羊养殖业发达的国家相比,我国绵羊有羊肉品质不高、产肉率低、生长速度慢等缺点,这均对我国绵羊养殖行业的可持续发展有着不可忽视的制约作用[3]。因此,选育出高产肉量和高肉质性状的肉羊品种是我国畜牧工作者不懈的追求。

骨骼肌是绵羊非常重要的器官,其良好的生长发育不仅保障了绵羊身体运动系统的正常运行,且是影响绵羊产肉量和产肉品质的直接靶器官。因此,绵羊羊肉产量和品质等经济性状培育和改良研究的重点在于深入了解绵羊肌肉的发育过程。在肌肉发育过程中,骨骼肌细胞的增殖与分化起着近乎决定性的作用[4-5]。众所周知,骨骼肌的发育和形成要经过极其复杂的生物学过程,绵羊出生后肌肉的生长并不表现于肌纤维数量的增加,而在于肌纤维的变大、变粗[6]。简单地理解,胎儿出生后,肌纤维的数量不会有净增加,仅能在原有的基础上使肌纤维变粗,所以胎儿期是骨骼肌发育的主要阶段[7]。Scicchitano等研究表明,胰岛素样生长因子1在调节肌原细胞增殖,肌小管和肌纤维细胞的合成代谢和肌肉的再生中发挥着核心作用[8]。Youngson等研究表明,Pax3是肌肉早期发育过程中必须的调控因子,其表达可抑制肌肉的分化[9]。Cesana等报道了肌肉特异性长非编码RNA linc-MD1在小鼠和人类成肌细胞中充当ceRNA来控制肌肉分化的时间,linc-MD1的下调会使肌肉分化程序延迟,而过表达与肌肉分化程序的预期相关[10],这些报道均为本次研究提供了重要的参考。

选取身处新疆高原牧场的哈萨克羊作为研究对象,该羊种具有耐低温、产肉多、增膘快、体型大、耐粗饲等优点,作为一种典型的肉脂两用型绵羊在新疆被大量养殖。深入研究调控该绵羊骨骼肌发育过程中的调控基因有助于提高该羊的产肉量、改善肉质。因此,本研究利用高通量测序结合生物信息学分析鉴定了哈萨克羊在胎儿期和成年期2个不同时期骨骼肌中差异表达的基因,旨在为深层次研究绵羊骨骼肌发育奠定基础。

1 材料与方法

1.1 试验动物

本项目于2021年10月在石河子大学生命科学学院完成,涉及动物的所有试验程序均经石河子大学动物卫生委员会批准。本研究在石河子市屠宰场采集了3头成年绵羊(母羊,2～3岁,体质量50～55 kg)的背最长肌,同时,采集了3头自然交配 110 d 后妊娠母羊的子宫,在无菌实验室中解剖子宫取出胎儿。后快速在胴体的最后一个胸椎两侧取下约12 g背最长肌。所有样品立即保存在不含RNA酶的超低温保存管中,于液氮中迅速冷冻,之后保存于-80 ℃冰箱直至抽提RNA。

1.2 文库构建及测序

将冷冻的组织从-80 ℃冰箱中取出,于液氮预冷后的研钵中研磨,后将其置于新的不含RNA酶的冻存管中,使用Trizol抽提总RNA。提取的RNA纯度和浓度使用RNA 6000 Nano kit进行检测。每个样品取等量的RNA制备材料,建库时使用3个同一时期的样本合并建库[11]。接着使用Epicentre Ribo-zeroTMrRNA试剂盒去除样本中的核糖体RNA。严格按照NEBNext UltraTM制备试剂盒说明书的步骤生成测序文库。利用AMPure XP系统纯化文库片段,选择长度为150～200 bp的cDNA片段。选择后的cDNA使用USER酶处理长度,接着利用高保真DNA聚合酶、随机引物进行PCR扩增。使用Agilent Bioanalyzer 2100系统进行文库鉴定。最后使用Illumina Hiseq 2500平台对制备好的文库进行测序产生125 bp配对末端Reads用于后续分析。

1.3 序列处理和分析流程

为保证后续生物信息学分析的质量,本研究使用FastQC软件对下机的fastq格式的原始数据进行质量控制。首先,去除包≥10%无法确定碱基信息的读数、Phred评分<5%的读数及含有接头序列的读数,得到高质量的清洁读数用于后续分析。此外,利用FastQC软件估测序列质量的好坏(Q20值),对质量进一步控制。使用Hisat2(v2.0.5)将清洁读数比对到绵羊参考基因组,Cufflinks软件将每个样品比对上的读数进行转录本的组装来拼接转录物,使用Cuffmerge去除含有多条重复转录本的基因并选择最长的一条转录本作为单基因簇,由此得到每个样品的单基因簇集,并对其进行后期分析[12-14]。

1.4 差异表达基因的筛选

转录本数据收集到后,使用Cuffdiff软件进行定量并分析差异表达的基因,FPKM均一化处理测序深度和基因长度不同的碎片,对基因表达的信息进行量化[15]。使用DESeq2 R包(1.16.1)中的负二项分布模型计算差异表达的基因,且通过q值的阈值来调整P值,阈值控制在q<0.01和|log2(差异倍数)|>1。

1.5 差异表达基因的富集分析

筛选出差异表达的基因后,使用GOseq R软件包,利用超几何检验P值,对这些基因进行GO富集分析,将其与GO数据库中的各个类别进行映射,分析映射在生物学过程、细胞组分及分子功能的基因数目。P值经过Bonferroni校正,校正后GO术语的P值<0.05,则认为是显著富集[16]。Cluster Profiler R软件包用于对差异表达的基因进行KEGG通路富集分析,同样显著富集的通路也通过超几何检验,P<0.05被认为显著富集。

2 结果与分析

2.1 哈萨克羊骨骼肌中RNA-seq分析的基本数据

为深入了解调控哈萨克羊骨骼肌发育的潜在主效基因,本研究使用转录组技术对胎儿期和成年期骨骼肌中的基因表达谱进行分析。使用Illumina HiSeq 2500对胎儿期和成年期哈萨克羊骨骼肌中基因的表达模式进行配对末端RNA测序。总地来说,在本次转录组研究中,共获得了206 390 670条原始序列。为保证生物信息学分析的准确性,使用FastQC对这些原始序列进行质控,按照方法中所述的质控流程处理后,用于后续分析的有199 001 830条清洁读数,数据有效率在96.4%以上。进一步对这些数据进行质量评估,发现Phred数值>20的高质量数据占总数的95.89%以上,GC含量占47.08%,结果表明测序质量良好,可用于后续分析。

2.2 哈萨克羊胎儿期和成年期骨骼肌中转录组的比较分析

为计算每个基因的表达水平,本研究计算了参考序列的总读数中清洁读数映射的数量及其在清洁读数中所占的百分比,包括多重映射和唯一映射的读数。为消除基因长度和测序深度带来的偏差,使用FPKM均一化mRNA的表达量。然后对胎儿期和成年期这2个不同时期样本中基因的表达水平进行分析。由图1-A可知,共有2 692个差异表达基因,其中,包含2 251个上调基因和441个下调基因。为从总体水平检查2组样品中基因表达的差异模式,使用2组样本中差异表达的mRNA进行层次聚类分析,由图1-B可知,虽未聚集在较明显的大类上,但差异基因可粗略地分为8个小类。进一步分析2组中差异表达的基因,由表1可知,某些调控肌肉发育的基因在胎儿期高显著表达。这些结果表明,肌肉中基因表达模式在胎儿期和成年期具有较大的差异,且两者中差异最大的基因是调控肌肉发育的基因。

表1 哈萨克羊胎儿期和成年期骨骼肌中差异表达最显著的前10个基因

2.3 哈萨克羊胎儿期和成年期骨骼肌中差异表达基因的聚类分析

为进一步理解差异表达基因的生物学功能,本研究使用GO和KEGG富集分析这些基因。GO注释结果表明,差异表达的基因在360个GO术语中被显著富集(P<0.05),其中,包含了275个生物过程、56个细胞成分和29个分子功能基因。其中,在前20个被显著富集的术语中,有解剖结构发育、发育过程、解剖结构形态发生、细胞周期、微管细胞骨架组织等生物过程,也有绑定等分子功能;由图2-A可知,蛋白质细胞外基质等细胞组分目录被显著富集(P<0.05)。结果表明,在绵羊骨骼肌发育的过程中可能涉及大量的细胞发育和细胞形态改变。

同时,本研究使用KEGG数据库对这些差异基因进行富集分析其通过基因产物的功能及相互作用的网络图,以期进一步理解差异基因的生物学功能。由图2-B可知,差异基因被显著富集(P<0.05)的共有32个术语,其中,包含了ECM-受体相互作用、PI3K-Akt信号通路、细胞周期、p53信号通路、黏着斑、钙信号通路、糖胺聚糖生物合成-硫酸肝素/肝素、间隙连接、PI3K-Akt信号通路、cGMP-PKG 信号通路和类固醇生物合成等能够调控骨骼肌发育的信号途径。

3 讨论与结论

目前,随着高通量测序技术的快速发展和成本的降低,该技术在家畜功能基因挖掘和优良性状遗传基础的研究中发挥着越来越多的作用。研究表明,利用转录组学手段可更加全面解析不同状态或不同处理下差异性状的调控基因,从而更容易揭示不同性状背后的遗传基础[17]。Jäger等利用RNA-seq技术研究了骨折延迟愈合的羊模型,发现了调控羊骨愈合过程的关键基因,进一步减少了羊分子模型中基因序列限制的缺失[18]。Fan等通过研究与绵羊毛色相关的皮肤转录组特征,扩展本研究对羊皮肤生理学和黑色素生成的复杂分子机制的理解[19]。Kogelman等适应转录组学在绵羊中检测到许多与肌肉和肠道寄生虫抵抗相关的通路[20]。这些研究者使用转录组测序技术对绵羊不同性状调控机制的解析丰富了本研究对这些性状背后遗传机制的认识,对绵羊的遗传育种具有显著的促进作用。然而,目前对于具有良好环境耐受性的哈萨克羊肌肉发育的分子调控机制仍不清楚。

作为哺乳动物最重要的器官,骨骼肌对动物的生理活动本身具有重要意义,同时这也是家畜动物优质蛋白质的重要来源。骨骼肌的发育直接影响着家畜的产肉量和肉质,并且骨骼肌的发育主要是在胚胎期和胎儿期,因此,深入研究胎儿骨骼肌发育过程中的主效基因对利用遗传育种的手段、选育优质高产的肉羊品种是必不可少的。本研究使用高通量测序技术在胎儿期和成年期哈萨克羊的骨骼肌中鉴定出了2 692个差异表达基因(P<0.05),结果也表明在胎儿期和成年期哈萨克羊骨骼肌中的基因表达模式有较大差异。进一步分析发现,在这些差异表达的基因中骨膜蛋白、肌球蛋白轻链4、Rho GTPase激活蛋白36、支链氨基酸转氨酶1、原纤蛋白2等在2组中的表达具有显著差异,且均在差异最显著的前10个基因中。进一步分析发现,骨膜蛋白能够支持肌生成和骨骼肌再生[21];肌球蛋白轻链4可结合钙离子,促进肌肉发育并参与横纹肌的收缩[22];Rho GTPase激活蛋白36是原性转录因子MyoG53的主要靶基因,在调控骨骼肌发育过程中发挥着重要的作用,并且该基因的表达受时间限制[23];支链氨基酸转氨酶1能够促进细胞增殖,并调控肌肉的发育和修复[24];原纤蛋白2缺失的小鼠中肌纤维的横截面会降低,胶原蛋白的交联水平也会降低[25]。此外,这5个基因均在胎儿期骨骼肌中高表达,表明这些基因在胎儿期绵羊骨骼肌中发挥了重要的调节作用。本研究结合这些基因的表达量和已知的功能推断,这5个基因在调控哈萨克羊骨骼肌发育过程中发挥着主要的调节作用,尽管具体的机制还需要未来利用分子生物学手段进行进一步研究。

此外,利用GO注释和KEGG富集分析的结果表明,大多数的基因涉及发育过程、微管发育、细胞周期等在发育过程中发挥重要调节作用的GO类别。同时,KEGG通路分析还显示,差异表达的mRNA参与了与骨骼肌发育直接相关的信号通路,包括ECM-受体相互作用、PI3K-Akt信号通路、细胞周期、p53信号通路、黏着斑、钙信号通路、糖胺聚糖生物合成-硫酸肝素/肝素、间隙连接、信号通路、cGMP-PKG信号通路和类固醇生物合成。ECM-受体相互作用信号通路参与了骨骼肌从胚胎阶段到衰老的骨骼肌的发育,并且参与了骨骼肌祖细胞的成肌过程[26]。p53信号通路对于成肌细胞增殖和凋亡是必不可少的[27]。钙信号通路控制肌肉发育,调控卫星细胞对骨骼肌的贡献[28]。骨骼肌中的cGMP-PKG信号通路和PI3K-Akt信号通路受营养调节,有助于肌细胞的生长[29]。类固醇生物合成信号通路通过调节类固醇激素的生物合成,从而参与骨骼肌的生长、发育和成熟[30]。这些结论表明,绵羊骨骼肌的发育涉及多种生物途径的共同调控。

总体来讲,本研究利用RNA-Seq技术对哈萨克羊胎儿期和成年期骨骼肌进行了转录组测序,构建了绵羊骨骼肌的表达谱,分析了胎儿期和成年期绵羊骨骼肌中差异表达的基因,其中,筛选到POSTN、MYL4、ARHGAP36、BCAT1、FBN25个基因潜在的调控哈萨克羊骨骼肌发育的主效基因。此外,这些差异基因主要集中在诸如ECM-受体相互作用、PI3K-Akt、钙信号通路和类固醇生物合成等直接参与骨骼肌发育的信号通路中,这表明这些差异是调控绵羊骨骼肌发育的关键因素。此外,本研究仍存在局限性,未来可利用分子生物学方法揭示绵羊骨骼肌发育的分子机制、利用荧光定量PCR试验验证骨骼肌不同生长发育状态下基因的表达量。