鸡肌肉发育相关miR-133a-5p分析和关键靶基因筛选

2023-07-15贾富民贾羽晴夏晶晶曹程程张泰康丁元飞郭立平耿照玉金四华

中国农业大学学报 2023年8期

贾富民贾羽晴夏晶晶税斐曹程程张泰康丁元飞郭立平耿照玉金四华

(安徽农业大学动物科技学院,合肥 230036)

MicroRNAs(miRNAs)作为一类内源性非编码小RNA,长度约为22 bp,广泛表达于各类细胞和组织中,并在生物发育中发挥调节作用[1],其发挥作用通常是通过与3′非翻译区(3′UTR)结合来抑制mRNA的翻译,调节靶基因的表达水平[2]。肌肉是由骨骼肌的肌肉细胞(也称为肌肉纤维),心肌和平滑肌组成。其中,骨骼肌被认为是最丰富,占成年动物体重的45%～60%[3-4]。肌肉发育一般分为2个阶段,胚胎期和出生后[5]。在胚胎期,肌肉祖细胞经过分化和增殖形成成肌细胞,然后融合形成多核肌管。最后,肌管发育成具有收缩特性的肌纤维[6]。几种肌源性miRNA,例如miR-206[7]、miR-1和miR-133a-3p[8],以及miR-203[9]、miR-181[10]和miR-128[11-12],被证明通过抑制靶基因表达来调节肌肉发育。miR-206在骨骼肌中特异表达,通过抑制与肌肉细胞融合相关的连接蛋白43(Cx43)基因的表达,促进肌肉细胞的分化[13]。miR-128-3p促进了鸡骨骼肌卫星细胞(SMSCs)的分化[1]。miR-1通过靶向组蛋白去乙酰化酶4(HDAC4)促进肌肉生成,HDAC4是肌肉基因表达的转录抑制因子[8]。在小鼠中,miR-1和miR-133a-3p可以通过负反馈调控机制与肌细胞增强子2(MEF2)和血清反应因子(SRF)一起调节骨骼肌的生长[8]。

大量研究报道,miR-133a对骨骼肌发育有着重要影响,其中,miR-133a-3p是肌肉中特异表达的miRNA,对肌纤维增殖与分化具有重要的调控作用[14]。研究发现,运动后大量外小体携带的miR-133a-3p和miR-133b等microRNA(miRNA)参与调节骨骼肌组织的增殖和分化[15]。Guo等[16]研究证实鸡miR-133a-3p通过靶向PRRX1调节成肌细胞增殖和分化。miR-133a-3p和miR-1a-3p通过抑制靶基因表达来调节Akt/mTOR/S6K信号通路,从而促进小鼠成肌细胞分化和骨骼肌生长[17]。miR-133a-3p促进了山羊肌肉干细胞细胞(MuSC)分化,且荧光素酶报告基因分析和功能分析确定miR-133a-3p直接靶向成纤维细胞生长因子受体1(FGFR1)[18]。而关于miR-133a-5p,同样被证明影响肌肉发育过程。Fu等[19]通过构建不同时期miRNA-mRNA相互作用网络,发现miR-133a-5p可能在鸡乳房肌肉发育中发挥重要作用。Chen等[20]研究证明miR-133a-5p在鸡成肌细胞中高表达,且抑制成肌细胞增殖和分化。通过已有研究发现,miR-133a-5p对鸡肌肉发育存在负面影响,但其具体机制尚不明晰。因此,本研究利用miRNA作用机制,结合生物信息学分析方法,筛选miR-133a-5p影响鸡肌肉发育的关键靶基因。

1 材料与方法

1.1 材料

在RNAcentral(https:∥rnacentral.org/)数据库查询鸡(Gallusgallus)、人(Homosapiens)、斑马鱼(Daniorerio)、小鼠(Musmusculus)、猪(Susscrofa)、山羊(Caprahircus)等不同物种的miR-133a-5p的序列。

从美国国立生物技术信息中心(NCBI)的GEO数据库(https:∥www.ncbi.nlm.nih.gov/geo/),下载数据集(GSE124418[21]、GSE133401[22]和GSE93855[23]),获取数据集中鸡不同组织的miR-133a-5p测序读数(Counts)、转录组Counts和关键靶基因的相对表达量。

从NCBI数据库中下载关键靶基因的序列,用于后续分析。

1.2 方法

1.2.1miR-133a-5p靶基因预测

利用Targetscan数据库(https:∥www.targetscan.org/vert_80/)和miRDB数据库(http:∥mirdb.org/mirdb/index.html)分别预测鸡miR-133a-5p的靶基因,并对结果取交集,减少假阳性,提高预测结果可信度。

1.2.2miR-133a-5p靶基因功能富集分析

利用基因本体数据库(Gene Ontology,GO)(http:∥geneontology.org/),对靶基因进行GO功能分析;利用KOBAS在线数据平台(http:∥bioinfo.org/kobas/),对靶基因进行京都基因与基因组百科全书(Kyoto encyclopedia of genes and genome,KEGG)通路分析。其中,GO功能分析包括生物学过程(Biological process,BP)、分子功能(Molecular function,MF)和细胞组成成分(Cell component,CC)。最后,利用Chiplot(https:∥www.chiplot.online/)将富集结果可视化。

1.2.3miR-133a-5p蛋白相互作用分析

利用STRING在线分析工具,分析靶基因编码蛋白的相互作用关系,并利用MCODE(Cytoscape的插件)鉴定网络图中有意义的模块。

1.2.4miR-133a-5p组织表达分析

整理GSE124418数据集中成熟清远麻鸡(Qingyuan partridge chicken)不同组织的miR-133a-5p的测序读数,利用以下公式将读数进行TPM(Transcripts per million)归一化。再将归一化的表达数据(TPM)进行单因素方差分析,判断肌肉相对各组织的表达量是否存在显著差异,并对结果进行可视化。

(1)

式中:T代表miRNA读数;N代表样本总miRNA读数。

1.2.5miR-133a-5p关键靶基因筛选

获取GSE133401数据集中成熟来航鸡(Leghorn chicken)不同组织的转录组数据,利用R包(DESeq2)分析肌肉相对其它组织的差异表达基因,并将肌肉相对于其他组织(肺部、肝脏和肾脏)的下调基因取交集,得到肌肉相对下调的公共基因集。所有Venn图均使用Hiplot(https:∥hiplot-academic.com/)绘制。

将上述基因集与miR-133a-5p靶基因集再一次取交集,得到关键靶基因。利用以下公式将关键靶基因的读数数据进行FPKM(Fragments per kilobase million)归一化,得到关键靶基因的表达谱,并用Chiplot进行聚类热图可视化。

(2)

式中:Exon Mapped Fragments为基因的读数;Total Mapped Fragments为样本的总读数;Exon Length为基因的有效长度。

1.2.6miR-133a-5p关键靶基因组织表达验证与结合位点评估

获取数据集GSE93855中低海拔地区清远麻鸡的关键靶基因的各组织相对表达量,观察肌肉相对各组织的表达差异。最后,使用RNA22 v2(https:∥cm.jefferson.edu/rna22/Interactive/)在线软件评估miR-133a-5p与关键靶基因3′UTR序列匹配的可靠性和能量稳定性。

1.3 数据分析

miR-133a-5p与其关键靶基因的组织表达数据,均使用Excel整理完成,随后使用GraphPad Prism 9进行单因素方差分析。

2 结果与分析

2.1 miR-133a-5p保守性分析

经查询,鸡miR-133a-5p位于2号染色体上,对比鸡、人、斑马鱼、小鼠、猪和山羊等物种的miR-133a-5p的序列,发现miR-133a-5p在进化中相对保守,成熟体序列基本一致(表1)。

表1 不同物种的miR-133a-5p序列

2.2 miR-133a-5p靶基因预测

经Targetscan数据库预测,获得鸡miR-133a-5p靶基因839个,经miRDB数据库预测得到鸡miR-133a-5p靶基因394个,为降低假阳性,将二者取交集,共得到157个交集靶基因(图1)。

图1 不同软件预测miR-133a-5p靶基因交集Venn图

2.3 miR-133a-5p靶基因富集分析

为探讨miR-133a-5p靶基因功能作用,使用GO数据库对其进行GO富集分析(图2)。对于BPs,靶基因主要参与调控或负调控生物合成、代谢等相关的条目中,包括细胞生物合成过程的调控、大分子生物合成过程的负调控和细胞代谢过程的调节等。在CCs富集中,靶基因显著富集于核、细胞内解剖结构和细胞内细胞器。而就MFs类别而言,靶基因富集结果最丰富的术语就是结合,包括酶结合、核酸结合、DNA结合及染色体结合等。

对靶基因进行KEGG途径通路富集(图3),显著富集包括ErbB信号通路(gga04012,P<0.01)、MAPK信号通路等信号通路(gga04010,P<0.01),也存在磷脂酰肌醇信号系统(gga04070,P<0.05)、肌动蛋白细胞骨架的调节(gga04810,P<0.05)等通路。

第一圈:显著富集路径,圈外代表基因数量的标尺。第二圈:背景基因通路的数量和P值。第三圈:miRNA靶点的路径数基因。第四圈:带有背景网格线的每条路径的丰度因子。每个网格表示0.1。 The first circle: indicates significant enrichment path, outside the circle represents the number of genes. The second circle: indicates number of background gene pathways and P values. The third circle: indicates path number genes of miRNA targets. The fourth circle: indicates abundance factor of each path with background grid lines. Each grid represents 0.1.

2.4 miR-133a-5p靶基因的蛋白相互作用

将靶基因导入至STRING在线软件,保持默认参数,分析其编码蛋白之间的相互作用,随后利用Cytoscape软件的MCODE插件,筛选网络图的中枢基因。共包括4个核心模块,17个基因,如PPARGC1A(PPARG coactivator 1 alpha)、HIF1A(Hypoxia inducible factor 1 alpha subunit)、CREB1(cAMP responsive element binding protein 1)、SIRT1(Sirtuin 1)和RPS6KB1(Ribosomal protein S6 kinase B1)等(图4)。

中心区域相同颜色与形状代表同一模块。 The center area with the same color and shape represents the same module.

2.5 miR-133a-5p在鸡各组织中的表达

从GSE124418数据集中获取miR-133a-5p的Counts数据,TPM归一化后,获得肌肉、心脏、肺部、肝脏、脾脏和肾脏的miR-133a-5p表达谱。由分析数据可知,与肺、脾脏、肾脏和肝脏相比,肌肉中的miR-133a-5p表达量显著上调(P<0.001)(图5)。

数据表示为3份样品的平均值±标准差。***P<0.001。 Data are expressed as the mean ± SD of triplicate samples. ***P<0.001.

2.6 肌肉相对各组织的差异表达基因

对GSE133401中获得鸡各组织的转录组数据进行差异分析,以|log2FC|>2,矫正P<0.05,为筛选条件,得到肌肉相对于肺部的下调基因为4 974个,肌肉相对于肝脏的下调基因为4 290个,肌肉相对于肾脏的下调基因为4 780个。由于缺少脾脏数据,仅取三者交集,得到肌肉相对下调的公共基因为1 790个(图6)。

图6 肌肉相对其他组织下调基因Venn图

2.7 miR-133a-5p关键靶基因筛选

依据miRNA与靶基因负相关的作用机制,将miR-133a-5p靶基因与肌肉相对下调的公共基因取交集,筛选出miR-133a-5p影响肌肉发育的关键靶基因(图7)共5个,SOX5、PRTFDC1、THRB、THBS2和ARRDC3。

图7 miR-133a-5p关键靶基因筛选Venn图

2.8 GSE133401数据集中miR-133a-5p关键靶基因表达

用热图展示miR-133a-5p关键靶基因表达谱(图8),各组织的不同样本间重复性好,肌肉与各组织组间差异显著(P<0.05)。聚类显示,ARRDC3和PRTFDC1首先聚成一束,随后分别与THBS2、SOX5和THRB依次聚类成束。

图8 miR-133a-5p关键靶基因不同组织相对表达量热图

2.9 GSE93855数据集中miR-133a-5p关键靶基因的表达

从GSE93855数据集中获取关键靶基因的相对表达量,ARRDC3和THRB的表达量均是在肌肉中最低;PRTFDC1的表达量在肝脏中最低,其次为肌肉,其原因可能是在肝脏中受其他生物学功能影响所致;而THBS2的表达在各组织中均处于一个较低的水平;同时,该数据集中未能发现SOX5(图9)。

数据表示为3份样品的平均值±标准差。*P<0.05,**P<0.01,***P<0.001,****P<0.000 1。 Data are expressed as the mean ±SD of triplicate samples. *P<0.05, **P<0.01, ***P<0.001, ****P<0.000 1.

2.10 miR-133a-5p关键靶基因结合位点验证

为进一步评估鸡miR-133a-5p与筛选关键靶基因的调控关系,将miR-133a-5p的序列和SOX5、PRTFDC1、THRB、THBS2和ARRDC3的3′UTR序列提交至RNA22 v2分析,结果发现仅有SOX5和THRB与鸡miR-133a-5p存在结合位点,且同时满足P<0.05(表2)。

表2 miR-133a-5p与靶基因结合位点

3 讨论

本研究对比不同物种miR-133a-5p的序列,发现其21 bp的碱基在不同物种的成熟序列均一致。说明miR-133a-5p在不同物种之间较为保守,可能存在相似功能。本研究重点关注miR-133a-5p对肌肉发育的影响,研究表明,miR-133a-5p是人横纹肌特异性miRNA,与骨骼肌功能障碍有关[24],也是大鼠肌肉减少症潜在的诊断靶点[25],在鸡[20]和猪[26]的肌肉中均高表达。同时,miR-133a-5p被证明抑制鸡肌肉发育[20],但具体作用机制尚不明晰。因此,进一步探究其影响肌肉发育的分子机制势在必行。

通过2种生物软件预测靶基因,取交集后[27],得到包含157个靶基因的数据集。GO富集分析显示,miR-133a-5p靶基因主要参与生物合成的调控和细胞代谢等相关条目,并且与酶、核酸、DNA及染色体结合相关。KEGG通路富集包含MAPK信号通路等多个信号通路,而显著富集的肌动蛋白细胞骨架等通路被证明与肌肉发育直接相关[28]。这也从侧面佐证了miR-133a-5p影响鸡肌肉发育。

本研究构建靶基因的蛋白相互作用,并且筛选出4个核心模块,涉及17个基因,结合后续的转录组测序数据,发现JARID2(Jumonji and AT-rich interaction domain containing 2)、PDIK1L(PDLIM1 interacting kinase 1 like)和ZCCHC8(Zinc finger CCHC-type containing 8)虽然不属于肌肉相对下调基因,但在肌肉中也有下调趋势,说明三者具有成为miR-133a-5p影响肌肉发育关键靶基因的潜力。目前也有学者探究了部分基因对肌肉发育的影响,如Adhikari等[29]发现JARID2通过直接抑制Wnt拮抗剂调节Wnt信号通路来调节肌肉分化。Bordini等[30]在研究与产生肌肉疾病相关的共表达基因簇(即模块)和关键基因座时,筛选到了ZCCHC8。

分析miR-133a-5p在鸡各组织中的表达谱,发现相较于肺、肾脏和肝脏等组织,miR-133a-5p在肌肉的表达显著上调,这符合miR-133a-5p参与肌肉发育的生物学过程的预期。Chen等[26]在研究猪肌肉时,同样发现了miR-133a-5p高表达。而依据miRNA与其靶基因负相关的作用机制,本研究结合鸡转录组数据分析,筛选到miR-133a-5p影响肌肉发育的5个关键靶基因,即SOX5、PRTFDC1、THRB、THBS2和ARRDC3。这些基因都是miR-133a-5p影响肌肉发育的潜在桥梁。Gordon等[31]研究表明ARRDC3表达的改变可能通过改变自噬激活来调节合成代谢和分解代谢刺激后的骨骼肌质量。Havis等[32]研究提到肌腱的进一步分化需要肌肉收缩来维持鸡肢体中THBS2的表达。PRTFDC1与肌肉相关的研究较少,值得后续进一步挖掘。

后续观察关键靶基因在GSE93855数据集中的表达谱,ARRDC3、THRB和PRTFDC1的表达量均在肌肉中显著降低。最后,使用RNA22 v2对miR-133a-5p与关键靶基因3′UTR序列匹配的可靠性和能量的稳定性进行评估,结果发现仅有SOX5、THRB存在P<0.05的结合位点。有学者研究发现,在小鼠发育过程中,SOX5在软骨细胞中表达强烈[33],但在鱼类的肌原细胞[34]、爪鼠[35-36]以及成人骨骼肌[37]中也检测到SOX5的表达,这表明SOX5可能在肌肉细胞中发挥作用。Della等[38]研究结果也表明,SOX5作为转录抑制因子间接增强爪蟾的肌源性转录。人体中TH受体(THRA或THRB)的表达对甲状腺激素(TH)影响骨骼肌收缩功能、肌肉生成和生物能量代谢十分重要[39]。Mitchell等[40]研究发现甲状腺激素受体β(THRB)突变与甲状腺激素抵抗(RTH)疾病相关,而在RTH患者的肌肉中,基础线粒体底物氧化增加,以ATP合成形式产生的能量减少。Johnsen等[41]发现绵羊的晚期妊娠营养不良(LG-UN)增加了心肌和背最长肌的TH受体(THRA和THRB)的基因表达。上述研究表明,SOX5与THRB在不同物种中皆对肌肉发育存在积极影响,再结合GSE93855的表达谱数据、RNA22 v2的评估结果以及miR-133a-5p对于鸡肌肉发育的抑制作用可推测,SOX5与THRB可能是miR-133a-5p影响鸡肌肉发育最为重要的中间桥梁。通过上述分析,为后续探索miR-133a-5p抑制肌肉发育的分子机制和验证试验提供了理论指导,还可为今后探索microRNAs的作用机制研究提供更多新思路。