两种蒙古绵羊25日龄胚胎转录组学差异分析

陈 陆,苏 红,王佳业,王大清,4,曹贵方*

(1.内蒙古农业大学兽医学院,内蒙古呼和浩特 010018;2.内蒙古自治区基础兽医学重点实验室,内蒙古呼和浩特 010018;3.动物胚胎与发育工程自治区高等学校重点实验室,内蒙古呼和浩特 010018;4.内蒙古自治区农牧业科学院,内蒙古呼和浩特 010010)

中国内蒙古地区是肉羊的主要养殖地,大约有17个类群,统称为蒙古羊(Mongolia sheep)。其中呼伦贝尔短尾羊主要分布在呼伦贝尔草原,因其尾短、产肉性能好、瘦肉率高等特点成为牧民养殖的良好选择,具有很高的经济价值和育种价值。而乌珠穆沁羊产于内蒙古锡林郭勒草原,属肉脂兼用型地方品种。其尾脂肥厚,具有耐严寒耐干旱等特点。然而,在当今现代化的生产系统中,绵羊被饲养在集约化或半精养的饲喂系统中,尾部脂肪沉积比积累等量的瘦肉组织需要更大的能量成本。因此,鉴定调节蒙古绵羊尾脂肪沉积的候选基因和分子途径在经济上具有重要意义。

先前的研究表明,不同品种的绵羊尾脂肪代谢与mRNA[1](信使核糖核酸)、miRNA(微小核糖核酸)和lncRNA(长链非编码RNA)有关。但是目前大量的工作都是围绕发育成熟的尾部脂肪组织展开。例如,Zhang W[2]等应用RNA-Seq技术测定了2个中国绵羊品种(肥臀阿勒泰羊和细尾新疆细毛羊)尾部脂肪组织的转录组谱,发现198个DEGS与脂肪代谢、脂肪细胞因子、PPAR信号传导通路和多种代谢途径有关。Wang X[3]等通过比较肥尾绵羊(哈萨克绵羊)和短尾绵羊(藏羊)脂肪组织转录组数据,筛选出与脂肪代谢相关的基因(NELL1和FMO3)。通路分析发现与脂肪酸代谢和细胞连接相关的途径在肥尾羊脂肪组织样本中上调。很少有研究通过比较不同尾型的蒙古绵羊胚胎发育期间的差异来确定与尾脂肪形成相关的重要功能基因。事实上,在绵羊胚胎发育期尾巴组织中已囤积了大量的脂肪组织,且25日龄时绵羊胚胎尾部生长加快,尾部体节已不清楚。

复杂且精密的细胞间信号传导和内在的转录程序决定着胚胎发育过程中尾组织的发生与形成。通过RNA测序(RNA-Seq)进行基因表达谱分析,为了解脂肪沉积的潜在机制提供了机会。因此,本研究从生物信息学角度揭示呼伦贝尔短尾羊和乌珠穆沁羊25日龄胚胎基因表达差异,研究两种胚胎基因、通路的变化等,以确定短尾的决定性基因以及控制蒙古绵羊尾脂肪沉积的潜在分子机制。同时有利于完善蒙古绵羊胚胎的发育情况,为研究早期胚胎发育的生命活动提供补充性材料。

1 材料与方法

1.1 材料

1.1.1 试验用动物 呼伦贝尔短尾羊选自内蒙古自治区呼伦贝尔市鄂温克族自治县,乌珠穆沁羊选自内蒙古自治区锡林浩特市东乌珠穆沁旗,为了保证获取胚胎的时间一致,将两种绵羊运至同一个牧场(呼和浩特市武川县种羊牧场),1个月后,母羊和公羊已适应当地气候环境,开始进行试验。

1.1.2 主要试剂 Trizol,逆转录酶,Invitrogen公司产品;琼脂糖,Biowest公司产品;TAE,Solarbio公司产品;镁离子打断试剂盒,E.coliDNA polymerase I,UDG酶(NEB公司产品;dUTP Solution,TaqDNA 聚合酶,Thermo Fisher公司产品;RNAfast200RNA快速提取试剂盒,上海飞捷公司产品,PrimeScript RT reagent Kit with gDNA Eraser,Takara公司产品;ChamQ Universal SYBR qPCR master mix试剂盒,南京诺唯赞公司产品;组织保存液,Miltenyi公司产品。

1.1.3 主要仪器 NanoDrop ND-1000,Thermo Fisher公司产品;Bioanalyzer 2100,Agilent公司产品;Centrifuge 5425R,Eppendorf 公司产品;VeritiProTM Thermal Cycler A48141,ABI公司产品;NovaseqTM6000,Illumina公司产品; Real-Time System CFX96,Bio-Rad公司产品;Oligo(dT)磁珠61002,Thermo Fisher公司产品。

1.2 方法

1.2.1 胚胎样本采集 试验前,用黄体酮阴道海绵栓塞处理20只乌珠穆沁和40只呼伦贝尔3岁健康母羊。在去掉栓塞同时肌肉注射PMSG(330 IU/只),发情母羊在12 h后交配。妊娠25 d分别采集乌珠穆沁羊和呼伦贝尔短尾羊的胚胎(图1A,图1B)。使用PBS洗涤胚胎3次,立即在室温下转移到组织保存液中。

A.呼伦贝尔短尾羊25日龄胚胎;B.乌珠穆沁羊25日胚胎。比例尺1∶1 000

1.2.2 RNA提取与建库测序 用Trizol,根据厂商提供的操作方案对总样品的RNA进行分离和纯化。然后用NanoDrop ND-1000通过测量RNA溶液在260 nm处的吸光值对总RNA的量与纯度进行质控。再通过Bioanalyzer 2 100对RNA的完整性进行检测。Oligo(dT)磁珠通过两轮的纯化对其中带有PolyA(多聚腺苷酸)的mRNA进行特异性捕获。再利用镁离子打断试剂盒进行片段化。将片段化的RNA在逆转录酶的作用下合成cDNA。然后使用E.coliDNA polymerase I,与RNase H进行双链合成,同时在二链中掺入dUTP Solution,将双链DNA的末端补齐为平末端。再在其两端各加上1个A碱基,再次利用磁珠对其片段大小进行筛选和纯化。以UDG酶消化二链,再通过PCR,使其形成片段大小为300 bp±50 bp的文库。最后,使用illumina NovaseqTM6 000按照标准操作对其进行双端测序,测序模式为PE150。

1.2.3 生物信息学分析 下机原始数据格式为fastq,使用fastp(https://github.com/OpenGene/fastp)软件对下机原始数据进行质控,包括去除接头、重复序列和低质量序列,参数为默认参数。使用HISAT2(https://ccb.jhu.edu/software/hisat2)将测序数据比对到基因组上(Ovis aries ARS-UI_Ramb_v2.0),得到文件格式为bam。使用StringTie软件(https://ccb.jhu.edu/software/hisat2) 对基因或转录本进行组装并用FPKM定量(FPKM=total-exon-fragments/mapped-reads(millions)×exon-length(kB)]),使用R包edgeR(https://bioconductor.org/packages/release/bioc/html/edgeR.html) 对样本之间进行显著差异分析,将差异倍数>2倍或<0.5倍且P value<0.05定义为差异基因。最后使用DAVID软件(https://david.ncifcrf.gov/) 对基因进行GO和KEGG富集分析。

1.2.4 RT-qPCR 使用RNAfast200从E25从乌珠穆沁和呼伦贝尔短尾羊中分离RNA。使用 ChamQ Universal SYBR qPCR master mix 进行qRT-PCR反应。在CFX96 PCR系统(Bio-rad)上使用ChamQ Universal SYBR qPCR进行实时PCR。使用管家基因GAPDH作为内部对照,基于ΔΔCT方法计算基因表达。随机挑选10个差异表达基因进行RT-qPCR反应,引物序列如表1所示。使用单因素方差分析和t检验进行统计分析(*P<0.05,**P<0.01,***P<0.001)。

表1 RT-qPCR引物信息

2 结果

2.1 RNA测序数据

从胎体样本的整个转录组中(n=3个胚胎25日龄的乌珠穆沁羊和n=3个胚胎25日龄的呼伦贝尔短尾羊),平均产生4 870万条读数/样本。经过质量控制,平均保留4 699万条读数/样本,平均有效测序量为96.48%,测序错误率小于0.01的平均概率为99.91%。此外,绝大部分的读数被映射到绵羊参考基因组(OvisariesARS-UI-Ramb-v2.0)。序列数据以Bioproject编号为GSE186602保存在NCBI数据库的序列读取档案中。(https://www.ncbi.nlm.nih.gov/geo/info/linking.html,2023年10月22日发布)

2.2 U(E25)和H(E25)胚胎之间的差异表达基因

在H( E25)与U( E25)的DEGs的比较中,具有显著性差异的基因有2721个,在呼伦贝尔短尾羊胚胎中有2320个上调基因,在乌珠穆沁羊胚胎中有401个上调基因(图2)。

蓝色(左)表示在乌珠穆沁羊中上调,红色(右)表示在呼伦贝尔短尾羊中上调

2.3 功能分析

通过对H(E25)中上调基因进行GO富集分析发现,这些基因可以被富集到132个GO功能注释中(P<0.01),包括86个生物学过程,占65.2%;29个分子功能,占21.9%;17个细胞组件,占12.9%。丰富的GO术语突出了发育生物学的几个关键生物过程,如轴突引导、突触发育及成熟、神经元投射、神经元分化等。通过对U(E25)中上调基因进行GO富集分析发现,这些基因可以被富集到81个GO功能注释中(P<0.01),包括40个生物学过程,占49.4%;23个分子功能,占28.4%;18个细胞组件,占22.2%。突出的生物过程有白细胞介素反应、皮肤分化、表皮发育、皮肤屏障的建立和胶原蛋白的合成等(见表2,表3)。将H(E25)和U(E25)中前20个的GO术语绘制成GO富集气泡图(图3)。

表2 H(E25)和U(E25)中显著性最高的20个GO功能注释

表3 H(E25)和U(E25)中显著性最高的20个上调基因

A.H(E25)中前20个GO Term:1.氧气输送; 2.载氧活性; 3.血红蛋白复合体; 4.氧气结合; 5.学习; 6.轴突触; 7.细胞骨架; 8.血红素生物合成法; 9.血影蛋白相关细胞骨架; 10.血液凝固; 11.蛋白质磷酸化; 12.细胞内信号传导; 13.核膜; 14.金属离子结合; 15.轴突导向; 16.纤维蛋白溶解; 17.跨膜转运蛋白活性; 18.含卟啉化合物生物合成方法; 19.跨膜运输;20.转录调控,DNA模板。B.U(E25)中前20个GO Term:1.细胞因子活性; 2.免疫反应; 3.基底膜; 4.核小体; 5.胶原蛋白结合; 6.趋化因子活性; 7.细胞趋化作用; 8.胞内 9.细胞外基质结构成分; 10.胞外; 11.大核糖体亚单位; 12.信号转导受体活性的调节; 13.含胶原蛋白的细胞外基质; 14.细胞外间隙; 15.胶原蛋白三聚体; 16.细胞外基质; 17.胶原纤维组织; 18.转录; 19.核糖体; 20.核糖体结构成分。点的颜色代表富集分析的P值,即富集的显著性,pvalue小于等于0.05表示显著富集。

对H(E25)上调基因进行KEGG富集分析发现,这些基因共被富集到56条信号通路(P<0.05),发现了Axon guidance通路和Hedgehog信号通路等参与胚胎早期发育的重要途径。对U(E25)上调基因进行KEGG富集分析发现,这些基因共被富集到40条信号通路(P<0.05),发现了一些重要的信号通路,如MAPK信号通路、蛋白质消化吸收和脂肪细胞分化及脂肪酸代谢相关的通路等。对差异表达的基因进行KEGG富集分析,选取富集分析结果中前20条信号通路绘制散点图(图4)。

通路名称:1.补体和凝血级联; 2.铁死亡; 3.糖尿病并发症中的AGE-RAGE信号通路; 4.癌症信号通路; 5.AMPK信号通路; 6.调节干细胞多能性的信号通路; 7.癌症中的胆碱代谢; 8.胰岛素抵抗; 9.非小细胞肺癌; 10.胆固醇代谢; 11.慢性粒细胞白血病; 12.苯丙氨酸代谢; 13.胆汁分泌; 14.卟啉和叶绿素代谢; 15.肝细胞癌; 16.磷脂酰肌醇信号传导系统; 17.急性髓系白血病; 18.PI3K-Akt信号通路; 19.MAPK信号通路; 20.p53信号通路。圆点大小表示富集到某个KEGG通路中差异基因的数量多少;圆点颜色表示不同的P值;Rich factor表示属于某个KEGG通路的差异基因个数/属于此KEGG通路的总基因数,Rich factor越大,KEGG通路富集程度越高。点的颜色代表富集分析的P值,即富集的显著性,P value小于等于0.05表示显著富集。

2.4 RT-qPCR分析

在mRNA测序数据中随机选取5个上调和5个下调的基因进行RT-qPCR反应,用来验证mRNA测序数据的准确性。10个差异表达mRNA的RT-qPCR结果显示,其表达水平差异显著(P<0.05),且与测序结果中表达的趋势一致,表明mRNA测序结果具有可靠性(图5)。

*P<0.05,**P<0.01,***P<0.001

3 讨 论

探究蒙古绵羊尾部脂肪堆积的分子机制,为进一步发展新品种育种提供理论基础。功能富集分析表明,大多数DEGs直接或间接参与脂质代谢相关途径。这意味着蒙古绵羊尾巴中脂肪沉积的分子机制由复杂的基因网络控制着。在上调DEGs中丰富的GO项和KEGG通路中,脂肪细胞分化与脂质代谢直接相关。根据相关的术语筛选到的DEG包括ZFP36和JUNB。ZFP36积极参与抑制脂解因子(例如IL-6)的表达,从而防止脂肪细胞降解并有利于其积累[3]。JUNB是一种转录因子(TF),其作用已在脂质代谢和脂肪细胞分化中被记录。以上研究表明ZFP36和JUNB的上调可能导致乌珠穆沁羊胚胎尾部的脂肪沉积过多。

GO分析表明,与白细胞介素反应相关的DEG在乌珠穆沁羊胚胎样本中上调。白细胞介素除了参与炎症反应外,最近的研究表明白细胞介素在脂质代谢中具有重要作用[4]。在本数据中,GO分析关于白细胞介素-1,2,6和8的反应显著富集,这与先前的发现一致,即IL-6[5-6]和IL-8[7],可通过对脂肪组织的直接作用来调节脂质代谢。与此反应相关的DEG包括JUNB和金属蛋白酶组织抑制剂2(TIMP2)等基因,TIMP2已被证明是人类肥胖的相关基因[8-9],可以抑制金属蛋白酶在脂肪细胞分化中的作用。与细尾(西藏)绵羊品种相比,发现在肥尾绵羊(哈萨克斯坦)的脂肪组织中TIMP2的表达量更高[10]。因此,这些基因可能与白细胞介素相互作用,以调节肥尾绵羊品种的脂肪细胞分化和脂肪沉积。

在乌珠穆沁羊胚胎样本中,观察到与MAPK信号通路相关的基因表达增加,包括生长停滞和DNA损伤诱导型β(GADD45B),丝裂原活化蛋白激酶8(MAP3K8)和FOS基因等。MAPK途径直接参与脂肪生成[11]或间接的诱导这些基因的表达。据黄[5]等人的报道,与荷斯坦品种的牛相比,MAPK信号通路在一种具有高脂肪量的牛品种中显著上调。与此同时,在肥尾与细尾藏羊品种的miRNA相关mRNA靶标中也富集了这一途径[12]。因此,目前研究的结果提供了进一步的证据,证明MAPK信号通路是有助于调节尾组织中脂肪沉积基因网络的重要组成部分。

3.1 呼伦贝短尾羊胚胎发育相关基因分析

随着神经系统的发育,新分化的神经元需要将其轴突向突触目标延伸,以形成功能性神经回路[13]。在胚胎发育中存在着各种导向信息影响着轴突的伸展方向。解读这些导向信息的是位于轴突末端高度能动和敏感的结构—生长锥,轴突在其延伸的顶端(生长锥)表达导向受体,并通过整合其环境中存在的吸引和排斥导向信息来导航至其目标[14]。在该研究中,H(E25)上调基因被富集到Axon guidance通路以及可以影响轴突和树突发育的基因FEZ1,CRMP1和DCC。

Hedgehog信号通路是一个高度进化保守的信号级联,可以协调胚胎发育过程中的多个关键步骤,包括胚胎模式、器官形态发生和调节细胞增殖、分化和迁移[15-16]。Hedgehog/Gli途径在刺激未分化祖细胞的细胞分裂和将细胞交付特定命运中起着关键作用。在中枢神经系统中,Hedgehog途径不仅是神经管模式的主要调节者,而且在中后脑的形成和小脑的生长模式中发挥一定的作用。Gli2A和Gli1可以控制前期中期后期脑内腹侧神经元的规格,以及晚期小脑内的小脑叶状结构和颗粒细胞群的扩张。除了神经管,Hedgehog途径对于正常的骨骼发育至关重要,Hedgehog通过3种Gli蛋白调节骨形成和修复。Gli2可以诱导软骨血管形成,这是未成熟软骨转化为成熟矿化组织的关键步骤[17]。有研究表明,Gli2突变的小鼠在骨骼发育中,颅面骨骼、椎骨、四肢长骨和胸骨处会表现出各种缺陷[18]。

3.2 乌珠穆沁羊胚胎发育相关基因分析

在U(E25)中上调基因较多富集到与胶原生物合成与共价交联有关的过程,富集到的差异基因包括COL1A1、COL1A2、COL3A1、COL4A5、COL4A6、COL5A2、COL5A1、COL6A3、COL11A2、COL12A1、COL22A1和COL24A1。胶原蛋白对组织的完整性至关重要。纤维状胶原是脊椎动物中最丰富的胶原,尽管胶原在组织中的含量非常少,但它们有助于组织分子结构的完整、塑造形状以及机械性能的实施,如皮肤的拉伸强度和韧带的抗拉性[19-20]。

还有部分上调基因被显著富集到与蛋白质消化吸收有关的过程,一些常见的基因,比如DPP4,也称为CD26,是一种Ⅱ型跨膜蛋白,以非经典分泌机制从细胞膜释放。这种外肽酶选择性地降解多种底物,包括肠促胰岛素激素、生长因子和细胞因子等[21]。XPNPEP2是氨肽酶蛋白家族成员,具有氨基酸切除功能,可将蛋白质水解为氨基酸[22]。ACE2作为重要的氨基酸转运因子[23],可以调节包括营养物质和药物在内的多种重要内源和外源物质的跨膜运输。

在这项研究中,通过RNA-Seq数据比较了2个具有不同尾型的蒙古绵羊品种,鉴定出一些含有脂质沉积相关基因的候选途径。结果表明,除了影响脂质代谢的途径外,与白细胞介素的反应和MAPK信号通路相关的一系列功能术语可能与脂质代谢相关的途径相互作用,有助于脂肪在乌珠穆沁羊尾巴中的沉积。仅通过基因表达谱很难了解脂肪沉积的具体调控途径。因此,需要基于研究结果进一步揭示羊尾脂肪沉积的确切遗传机制。同时,比较了2种蒙古绵羊胚胎的发育情况,发现H(E25)突出了神经和骨组织的发育,而U(E25)侧重于与胶原和蛋白质相关的生命过程。