雞肌肉發育相關miR-133a-5p分析和關鍵靶基因篩選

賈富民 賈羽晴 夏晶晶 稅 斐 曹程程 張泰康 丁元飛 郭立平 耿照玉 金四華

(安徽農業大學 動物科技學院,合肥 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的作用機制研究提供更多新思路。

4 結 論

miR-133a-5p在雞肌肉中顯著表達,其在肌肉中發揮生物學功能可能是通過調節靶基因SOX5、PRTFDC1、THRB、THBS2和ARRDC3表達實現的,其中,SOX5和THRB尤為重要。