雞肌肉發(fā)育相關(guān)miR-133a-5p分析和關(guān)鍵靶基因篩選

賈富民 賈羽晴 夏晶晶 稅 斐 曹程程 張?zhí)┛?丁元飛 郭立平 耿照玉 金四華

(安徽農(nóng)業(yè)大學(xué) 動(dòng)物科技學(xué)院,合肥 230036)

MicroRNAs(miRNAs)作為一類內(nèi)源性非編碼小RNA,長(zhǎng)度約為22 bp,廣泛表達(dá)于各類細(xì)胞和組織中,并在生物發(fā)育中發(fā)揮調(diào)節(jié)作用[1],其發(fā)揮作用通常是通過與3′非翻譯區(qū)(3′UTR)結(jié)合來(lái)抑制mRNA的翻譯,調(diào)節(jié)靶基因的表達(dá)水平[2]。肌肉是由骨骼肌的肌肉細(xì)胞(也稱為肌肉纖維),心肌和平滑肌組成。其中,骨骼肌被認(rèn)為是最豐富,占成年動(dòng)物體重的45%～60%[3-4]。肌肉發(fā)育一般分為2個(gè)階段,胚胎期和出生后[5]。在胚胎期,肌肉祖細(xì)胞經(jīng)過分化和增殖形成成肌細(xì)胞,然后融合形成多核肌管。最后,肌管發(fā)育成具有收縮特性的肌纖維[6]。幾種肌源性miRNA,例如miR-206[7]、miR-1和miR-133a-3p[8],以及miR-203[9]、miR-181[10]和miR-128[11-12],被證明通過抑制靶基因表達(dá)來(lái)調(diào)節(jié)肌肉發(fā)育。miR-206在骨骼肌中特異表達(dá),通過抑制與肌肉細(xì)胞融合相關(guān)的連接蛋白43(Cx43)基因的表達(dá),促進(jìn)肌肉細(xì)胞的分化[13]。miR-128-3p促進(jìn)了雞骨骼肌衛(wèi)星細(xì)胞(SMSCs)的分化[1]。miR-1通過靶向組蛋白去乙?；?(HDAC4)促進(jìn)肌肉生成,HDAC4是肌肉基因表達(dá)的轉(zhuǎn)錄抑制因子[8]。在小鼠中,miR-1和miR-133a-3p可以通過負(fù)反饋調(diào)控機(jī)制與肌細(xì)胞增強(qiáng)子2(MEF2)和血清反應(yīng)因子(SRF)一起調(diào)節(jié)骨骼肌的生長(zhǎng)[8]。

大量研究報(bào)道,miR-133a對(duì)骨骼肌發(fā)育有著重要影響,其中,miR-133a-3p是肌肉中特異表達(dá)的miRNA,對(duì)肌纖維增殖與分化具有重要的調(diào)控作用[14]。研究發(fā)現(xiàn),運(yùn)動(dòng)后大量外小體攜帶的miR-133a-3p和miR-133b等microRNA(miRNA)參與調(diào)節(jié)骨骼肌組織的增殖和分化[15]。Guo等[16]研究證實(shí)雞miR-133a-3p通過靶向PRRX1調(diào)節(jié)成肌細(xì)胞增殖和分化。miR-133a-3p和miR-1a-3p通過抑制靶基因表達(dá)來(lái)調(diào)節(jié)Akt/mTOR/S6K信號(hào)通路,從而促進(jìn)小鼠成肌細(xì)胞分化和骨骼肌生長(zhǎng)[17]。miR-133a-3p促進(jìn)了山羊肌肉干細(xì)胞細(xì)胞(MuSC)分化,且熒光素酶報(bào)告基因分析和功能分析確定miR-133a-3p直接靶向成纖維細(xì)胞生長(zhǎng)因子受體1(FGFR1)[18]。而關(guān)于miR-133a-5p,同樣被證明影響肌肉發(fā)育過程。Fu等[19]通過構(gòu)建不同時(shí)期miRNA-mRNA相互作用網(wǎng)絡(luò),發(fā)現(xiàn)miR-133a-5p可能在雞乳房肌肉發(fā)育中發(fā)揮重要作用。Chen等[20]研究證明miR-133a-5p在雞成肌細(xì)胞中高表達(dá),且抑制成肌細(xì)胞增殖和分化。通過已有研究發(fā)現(xiàn),miR-133a-5p對(duì)雞肌肉發(fā)育存在負(fù)面影響,但其具體機(jī)制尚不明晰。因此,本研究利用miRNA作用機(jī)制,結(jié)合生物信息學(xué)分析方法,篩選miR-133a-5p影響雞肌肉發(fā)育的關(guān)鍵靶基因。

1 材料與方法

1.1 材料

在RNAcentral(https:∥rnacentral.org/)數(shù)據(jù)庫(kù)查詢雞(Gallusgallus)、人(Homosapiens)、斑馬魚(Daniorerio)、小鼠(Musmusculus)、豬(Susscrofa)、山羊(Caprahircus)等不同物種的miR-133a-5p的序列。

從美國(guó)國(guó)立生物技術(shù)信息中心(NCBI)的GEO數(shù)據(jù)庫(kù)(https:∥www.ncbi.nlm.nih.gov/geo/),下載數(shù)據(jù)集(GSE124418[21]、GSE133401[22]和GSE93855[23]),獲取數(shù)據(jù)集中雞不同組織的miR-133a-5p測(cè)序讀數(shù)(Counts)、轉(zhuǎn)錄組Counts和關(guān)鍵靶基因的相對(duì)表達(dá)量。

從NCBI數(shù)據(jù)庫(kù)中下載關(guān)鍵靶基因的序列,用于后續(xù)分析。

1.2 方法

1.2.1miR-133a-5p靶基因預(yù)測(cè)

利用Targetscan數(shù)據(jù)庫(kù)(https:∥www.targetscan.org/vert_80/)和miRDB數(shù)據(jù)庫(kù)(http:∥mirdb.org/mirdb/index.html)分別預(yù)測(cè)雞miR-133a-5p的靶基因,并對(duì)結(jié)果取交集,減少假陽(yáng)性,提高預(yù)測(cè)結(jié)果可信度。

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

利用基因本體數(shù)據(jù)庫(kù)(Gene Ontology,GO)(http:∥geneontology.org/),對(duì)靶基因進(jìn)行GO功能分析;利用KOBAS在線數(shù)據(jù)平臺(tái)(http:∥bioinfo.org/kobas/),對(duì)靶基因進(jìn)行京都基因與基因組百科全書(Kyoto encyclopedia of genes and genome,KEGG)通路分析。其中,GO功能分析包括生物學(xué)過程(Biological process,BP)、分子功能(Molecular function,MF)和細(xì)胞組成成分(Cell component,CC)。最后,利用Chiplot(https:∥www.chiplot.online/)將富集結(jié)果可視化。

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

利用STRING在線分析工具,分析靶基因編碼蛋白的相互作用關(guān)系,并利用MCODE(Cytoscape的插件)鑒定網(wǎng)絡(luò)圖中有意義的模塊。

1.2.4miR-133a-5p組織表達(dá)分析

整理GSE124418數(shù)據(jù)集中成熟清遠(yuǎn)麻雞(Qingyuan partridge chicken)不同組織的miR-133a-5p的測(cè)序讀數(shù),利用以下公式將讀數(shù)進(jìn)行TPM(Transcripts per million)歸一化。再將歸一化的表達(dá)數(shù)據(jù)(TPM)進(jìn)行單因素方差分析,判斷肌肉相對(duì)各組織的表達(dá)量是否存在顯著差異,并對(duì)結(jié)果進(jìn)行可視化。

(1)

式中:T代表miRNA讀數(shù);N代表樣本總miRNA讀數(shù)。

1.2.5miR-133a-5p關(guān)鍵靶基因篩選

獲取GSE133401數(shù)據(jù)集中成熟來(lái)航雞(Leghorn chicken)不同組織的轉(zhuǎn)錄組數(shù)據(jù),利用R包(DESeq2)分析肌肉相對(duì)其它組織的差異表達(dá)基因,并將肌肉相對(duì)于其他組織(肺部、肝臟和腎臟)的下調(diào)基因取交集,得到肌肉相對(duì)下調(diào)的公共基因集。所有Venn圖均使用Hiplot(https:∥hiplot-academic.com/)繪制。

將上述基因集與miR-133a-5p靶基因集再一次取交集,得到關(guān)鍵靶基因。利用以下公式將關(guān)鍵靶基因的讀數(shù)數(shù)據(jù)進(jìn)行FPKM(Fragments per kilobase million)歸一化,得到關(guān)鍵靶基因的表達(dá)譜,并用Chiplot進(jìn)行聚類熱圖可視化。

(2)

式中:Exon Mapped Fragments為基因的讀數(shù);Total Mapped Fragments為樣本的總讀數(shù);Exon Length為基因的有效長(zhǎng)度。

1.2.6miR-133a-5p關(guān)鍵靶基因組織表達(dá)驗(yàn)證與結(jié)合位點(diǎn)評(píng)估

獲取數(shù)據(jù)集GSE93855中低海拔地區(qū)清遠(yuǎn)麻雞的關(guān)鍵靶基因的各組織相對(duì)表達(dá)量,觀察肌肉相對(duì)各組織的表達(dá)差異。最后,使用RNA22 v2(https:∥cm.jefferson.edu/rna22/Interactive/)在線軟件評(píng)估m(xù)iR-133a-5p與關(guān)鍵靶基因3′UTR序列匹配的可靠性和能量穩(wěn)定性。

1.3 數(shù)據(jù)分析

miR-133a-5p與其關(guān)鍵靶基因的組織表達(dá)數(shù)據(jù),均使用Excel整理完成,隨后使用GraphPad Prism 9進(jìn)行單因素方差分析。

2 結(jié)果與分析

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

經(jīng)查詢,雞miR-133a-5p位于2號(hào)染色體上,對(duì)比雞、人、斑馬魚、小鼠、豬和山羊等物種的miR-133a-5p的序列,發(fā)現(xiàn)miR-133a-5p在進(jìn)化中相對(duì)保守,成熟體序列基本一致(表1)。

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

2.2 miR-133a-5p靶基因預(yù)測(cè)

經(jīng)Targetscan數(shù)據(jù)庫(kù)預(yù)測(cè),獲得雞miR-133a-5p靶基因839個(gè),經(jīng)miRDB數(shù)據(jù)庫(kù)預(yù)測(cè)得到雞miR-133a-5p靶基因394個(gè),為降低假陽(yáng)性,將二者取交集,共得到157個(gè)交集靶基因(圖1)。

圖1 不同軟件預(yù)測(cè)miR-133a-5p靶基因交集Venn圖

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

為探討miR-133a-5p靶基因功能作用,使用GO數(shù)據(jù)庫(kù)對(duì)其進(jìn)行GO富集分析(圖2)。對(duì)于BPs,靶基因主要參與調(diào)控或負(fù)調(diào)控生物合成、代謝等相關(guān)的條目中,包括細(xì)胞生物合成過程的調(diào)控、大分子生物合成過程的負(fù)調(diào)控和細(xì)胞代謝過程的調(diào)節(jié)等。在CCs富集中,靶基因顯著富集于核、細(xì)胞內(nèi)解剖結(jié)構(gòu)和細(xì)胞內(nèi)細(xì)胞器。而就MFs類別而言,靶基因富集結(jié)果最豐富的術(shù)語(yǔ)就是結(jié)合,包括酶結(jié)合、核酸結(jié)合、DNA結(jié)合及染色體結(jié)合等。

對(duì)靶基因進(jìn)行KEGG途徑通路富集(圖3),顯著富集包括ErbB信號(hào)通路(gga04012,P<0.01)、MAPK信號(hào)通路等信號(hào)通路(gga04010,P<0.01),也存在磷脂酰肌醇信號(hào)系統(tǒng)(gga04070,P<0.05)、肌動(dòng)蛋白細(xì)胞骨架的調(diào)節(jié)(gga04810,P<0.05)等通路。

第一圈:顯著富集路徑,圈外代表基因數(shù)量的標(biāo)尺。第二圈:背景基因通路的數(shù)量和P值。第三圈:miRNA靶點(diǎn)的路徑數(shù)基因。第四圈:帶有背景網(wǎng)格線的每條路徑的豐度因子。每個(gè)網(wǎng)格表示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靶基因的蛋白相互作用

將靶基因?qū)胫罶TRING在線軟件,保持默認(rèn)參數(shù),分析其編碼蛋白之間的相互作用,隨后利用Cytoscape軟件的MCODE插件,篩選網(wǎng)絡(luò)圖的中樞基因。共包括4個(gè)核心模塊,17個(gè)基因,如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)。

中心區(qū)域相同顏色與形狀代表同一模塊。 The center area with the same color and shape represents the same module.

2.5 miR-133a-5p在雞各組織中的表達(dá)

從GSE124418數(shù)據(jù)集中獲取miR-133a-5p的Counts數(shù)據(jù),TPM歸一化后,獲得肌肉、心臟、肺部、肝臟、脾臟和腎臟的miR-133a-5p表達(dá)譜。由分析數(shù)據(jù)可知,與肺、脾臟、腎臟和肝臟相比,肌肉中的miR-133a-5p表達(dá)量顯著上調(diào)(P<0.001)(圖5)。

數(shù)據(jù)表示為3份樣品的平均值±標(biāo)準(zhǔn)差。***P<0.001。 Data are expressed as the mean ± SD of triplicate samples. ***P<0.001.

2.6 肌肉相對(duì)各組織的差異表達(dá)基因

對(duì)GSE133401中獲得雞各組織的轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行差異分析,以|log2FC|>2,矯正P<0.05,為篩選條件,得到肌肉相對(duì)于肺部的下調(diào)基因?yàn)? 974個(gè),肌肉相對(duì)于肝臟的下調(diào)基因?yàn)? 290個(gè),肌肉相對(duì)于腎臟的下調(diào)基因?yàn)? 780個(gè)。由于缺少脾臟數(shù)據(jù),僅取三者交集,得到肌肉相對(duì)下調(diào)的公共基因?yàn)? 790個(gè)(圖6)。

圖6 肌肉相對(duì)其他組織下調(diào)基因Venn圖

2.7 miR-133a-5p關(guān)鍵靶基因篩選

依據(jù)miRNA與靶基因負(fù)相關(guān)的作用機(jī)制,將miR-133a-5p靶基因與肌肉相對(duì)下調(diào)的公共基因取交集,篩選出miR-133a-5p影響肌肉發(fā)育的關(guān)鍵靶基因(圖7)共5個(gè),SOX5、PRTFDC1、THRB、THBS2和ARRDC3。

圖7 miR-133a-5p關(guān)鍵靶基因篩選Venn圖

2.8 GSE133401數(shù)據(jù)集中miR-133a-5p關(guān)鍵靶基因表達(dá)

用熱圖展示miR-133a-5p關(guān)鍵靶基因表達(dá)譜(圖8),各組織的不同樣本間重復(fù)性好,肌肉與各組織組間差異顯著(P<0.05)。聚類顯示,ARRDC3和PRTFDC1首先聚成一束,隨后分別與THBS2、SOX5和THRB依次聚類成束。

圖8 miR-133a-5p關(guān)鍵靶基因不同組織相對(duì)表達(dá)量熱圖

2.9 GSE93855數(shù)據(jù)集中miR-133a-5p關(guān)鍵靶基因的表達(dá)

從GSE93855數(shù)據(jù)集中獲取關(guān)鍵靶基因的相對(duì)表達(dá)量,ARRDC3和THRB的表達(dá)量均是在肌肉中最低;PRTFDC1的表達(dá)量在肝臟中最低,其次為肌肉,其原因可能是在肝臟中受其他生物學(xué)功能影響所致;而THBS2的表達(dá)在各組織中均處于一個(gè)較低的水平;同時(shí),該數(shù)據(jù)集中未能發(fā)現(xiàn)SOX5(圖9)。

數(shù)據(jù)表示為3份樣品的平均值±標(biāo)準(zhǔn)差。*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關(guān)鍵靶基因結(jié)合位點(diǎn)驗(yàn)證

為進(jìn)一步評(píng)估雞miR-133a-5p與篩選關(guān)鍵靶基因的調(diào)控關(guān)系,將miR-133a-5p的序列和SOX5、PRTFDC1、THRB、THBS2和ARRDC3的3′UTR序列提交至RNA22 v2分析,結(jié)果發(fā)現(xiàn)僅有SOX5和THRB與雞miR-133a-5p存在結(jié)合位點(diǎn),且同時(shí)滿足P<0.05(表2)。

表2 miR-133a-5p與靶基因結(jié)合位點(diǎn)

3 討 論

本研究對(duì)比不同物種miR-133a-5p的序列,發(fā)現(xiàn)其21 bp的堿基在不同物種的成熟序列均一致。說(shuō)明miR-133a-5p在不同物種之間較為保守,可能存在相似功能。本研究重點(diǎn)關(guān)注miR-133a-5p對(duì)肌肉發(fā)育的影響,研究表明,miR-133a-5p是人橫紋肌特異性miRNA,與骨骼肌功能障礙有關(guān)[24],也是大鼠肌肉減少癥潛在的診斷靶點(diǎn)[25],在雞[20]和豬[26]的肌肉中均高表達(dá)。同時(shí),miR-133a-5p被證明抑制雞肌肉發(fā)育[20],但具體作用機(jī)制尚不明晰。因此,進(jìn)一步探究其影響肌肉發(fā)育的分子機(jī)制勢(shì)在必行。

通過2種生物軟件預(yù)測(cè)靶基因,取交集后[27],得到包含157個(gè)靶基因的數(shù)據(jù)集。GO富集分析顯示,miR-133a-5p靶基因主要參與生物合成的調(diào)控和細(xì)胞代謝等相關(guān)條目,并且與酶、核酸、DNA及染色體結(jié)合相關(guān)。KEGG通路富集包含MAPK信號(hào)通路等多個(gè)信號(hào)通路,而顯著富集的肌動(dòng)蛋白細(xì)胞骨架等通路被證明與肌肉發(fā)育直接相關(guān)[28]。這也從側(cè)面佐證了miR-133a-5p影響雞肌肉發(fā)育。

本研究構(gòu)建靶基因的蛋白相互作用,并且篩選出4個(gè)核心模塊,涉及17個(gè)基因,結(jié)合后續(xù)的轉(zhuǎn)錄組測(cè)序數(shù)據(jù),發(fā)現(xiàn)JARID2(Jumonji and AT-rich interaction domain containing 2)、PDIK1L(PDLIM1 interacting kinase 1 like)和ZCCHC8(Zinc finger CCHC-type containing 8)雖然不屬于肌肉相對(duì)下調(diào)基因,但在肌肉中也有下調(diào)趨勢(shì),說(shuō)明三者具有成為miR-133a-5p影響肌肉發(fā)育關(guān)鍵靶基因的潛力。目前也有學(xué)者探究了部分基因?qū)∪獍l(fā)育的影響,如Adhikari等[29]發(fā)現(xiàn)JARID2通過直接抑制Wnt拮抗劑調(diào)節(jié)Wnt信號(hào)通路來(lái)調(diào)節(jié)肌肉分化。Bordini等[30]在研究與產(chǎn)生肌肉疾病相關(guān)的共表達(dá)基因簇(即模塊)和關(guān)鍵基因座時(shí),篩選到了ZCCHC8。

分析miR-133a-5p在雞各組織中的表達(dá)譜,發(fā)現(xiàn)相較于肺、腎臟和肝臟等組織,miR-133a-5p在肌肉的表達(dá)顯著上調(diào),這符合miR-133a-5p參與肌肉發(fā)育的生物學(xué)過程的預(yù)期。Chen等[26]在研究豬肌肉時(shí),同樣發(fā)現(xiàn)了miR-133a-5p高表達(dá)。而依據(jù)miRNA與其靶基因負(fù)相關(guān)的作用機(jī)制,本研究結(jié)合雞轉(zhuǎn)錄組數(shù)據(jù)分析,篩選到miR-133a-5p影響肌肉發(fā)育的5個(gè)關(guān)鍵靶基因,即SOX5、PRTFDC1、THRB、THBS2和ARRDC3。這些基因都是miR-133a-5p影響肌肉發(fā)育的潛在橋梁。Gordon等[31]研究表明ARRDC3表達(dá)的改變可能通過改變自噬激活來(lái)調(diào)節(jié)合成代謝和分解代謝刺激后的骨骼肌質(zhì)量。Havis等[32]研究提到肌腱的進(jìn)一步分化需要肌肉收縮來(lái)維持雞肢體中THBS2的表達(dá)。PRTFDC1與肌肉相關(guān)的研究較少,值得后續(xù)進(jìn)一步挖掘。

后續(xù)觀察關(guān)鍵靶基因在GSE93855數(shù)據(jù)集中的表達(dá)譜,ARRDC3、THRB和PRTFDC1的表達(dá)量均在肌肉中顯著降低。最后,使用RNA22 v2對(duì)miR-133a-5p與關(guān)鍵靶基因3′UTR序列匹配的可靠性和能量的穩(wěn)定性進(jìn)行評(píng)估,結(jié)果發(fā)現(xiàn)僅有SOX5、THRB存在P<0.05的結(jié)合位點(diǎn)。有學(xué)者研究發(fā)現(xiàn),在小鼠發(fā)育過程中,SOX5在軟骨細(xì)胞中表達(dá)強(qiáng)烈[33],但在魚類的肌原細(xì)胞[34]、爪鼠[35-36]以及成人骨骼肌[37]中也檢測(cè)到SOX5的表達(dá),這表明SOX5可能在肌肉細(xì)胞中發(fā)揮作用。Della等[38]研究結(jié)果也表明,SOX5作為轉(zhuǎn)錄抑制因子間接增強(qiáng)爪蟾的肌源性轉(zhuǎn)錄。人體中TH受體(THRA或THRB)的表達(dá)對(duì)甲狀腺激素(TH)影響骨骼肌收縮功能、肌肉生成和生物能量代謝十分重要[39]。Mitchell等[40]研究發(fā)現(xiàn)甲狀腺激素受體β(THRB)突變與甲狀腺激素抵抗(RTH)疾病相關(guān),而在RTH患者的肌肉中,基礎(chǔ)線粒體底物氧化增加,以ATP合成形式產(chǎn)生的能量減少。Johnsen等[41]發(fā)現(xiàn)綿羊的晚期妊娠營(yíng)養(yǎng)不良(LG-UN)增加了心肌和背最長(zhǎng)肌的TH受體(THRA和THRB)的基因表達(dá)。上述研究表明,SOX5與THRB在不同物種中皆對(duì)肌肉發(fā)育存在積極影響,再結(jié)合GSE93855的表達(dá)譜數(shù)據(jù)、RNA22 v2的評(píng)估結(jié)果以及miR-133a-5p對(duì)于雞肌肉發(fā)育的抑制作用可推測(cè),SOX5與THRB可能是miR-133a-5p影響雞肌肉發(fā)育最為重要的中間橋梁。通過上述分析,為后續(xù)探索miR-133a-5p抑制肌肉發(fā)育的分子機(jī)制和驗(yàn)證試驗(yàn)提供了理論指導(dǎo),還可為今后探索microRNAs的作用機(jī)制研究提供更多新思路。

4 結(jié) 論

miR-133a-5p在雞肌肉中顯著表達(dá),其在肌肉中發(fā)揮生物學(xué)功能可能是通過調(diào)節(jié)靶基因SOX5、PRTFDC1、THRB、THBS2和ARRDC3表達(dá)實(shí)現(xiàn)的,其中,SOX5和THRB尤為重要。