續苓健骨方對骨質疏松模型大鼠miRNA表達譜的作用研究

李生強 陳賽楠 謝冰穎 陳娟 謝麗華 葉云金 黃景文 葛繼榮

福建省中醫藥科學院骨質疏松證候基因組學研究室，福建 福州 350003

骨質疏松癥(osteoporosis，OP)可造成骨折及骨折并發癥引起的殘疾甚至死亡，骨質疏松的防治已成為老年人健康的重大課題。近年來，中藥治療骨質疏松癥得到廣泛認可。續苓健骨方是以中醫藥理論遣方用藥形成的治療OP的經驗方。前期的動物實驗表明，續苓健骨方可提高骨質疏松模型大鼠的骨密度[1]，改善骨生物力學性能[2]，提升血鈣含量，降低血磷含量[3]，并對其機制進行了一定探索[4-5]。miRNA參與了骨的生長、發育、代謝及藥物治療等各個過程[6-7]。本研究采用二代測序方法，對骨質疏松大鼠模型在續苓健骨方治療后進行腰椎miRNA測序，探討小分子非編碼RNA在治療過程中的變化，為臨床應用該方治療OP提供實驗室依據。

1 材料與方法

1.1 材料

1.1.1實驗動物：6月齡SPF級SD雌性大鼠18只，購自上海斯萊克實驗動物有限責任公司，實驗動物許可證號：SCXK(滬)2007-0005，合格證編號：2007000530078。飼養于福建省中醫藥科學院比較醫學實驗中心，實驗動物許可證：SYXK(閩)2016-0005。

1.1.2實驗藥物與試劑：續苓健骨方(國家發明專利號ZL201510784227.4)由川續斷、白術、陳皮、赤芍、紅花、甘草等12味中藥組成，共122 g/劑。中藥飲片購自福建省醫藥公司，按傳統中藥煎熬方法，每劑取汁 81.33 mL，生藥含量為1.5 g/mL，4 ℃冰箱保存備用。大鼠灌胃量計算方法：按成人50 kg：122/50=2.44 g/(kg·d)，大鼠灌胃量：2.44×6.25倍=15 g/kg。

文庫構建、測序相關試劑：NEB Next?Poly(A) mRNA 磁性分離模塊、NEB Multiplex Small RNA Library Prep Set(美國New England Biolabs公司)；RiboZero Magnetic Gold Kit (Human/Mouse/Rat) (美國Epicentre公司)；TruSeq Rapid SR Cluster Kit (#GD-402-4001，美國Illumina公司)。

1.1.3儀器：Discovery W 雙能 X 線骨密度儀(變異系數1.0 CV%，精度0.25%，美國Hologic 公司)；NanoDrop ND-1000(美國Thermal Scientific公司)；生物芯片分析系統Agilent 2100 Bioanalyzer(美國Agilent公司)；Illumina NextSeq 500 測序儀(美國Illumina公司)。

1.2 方法

1.2.1動物分組、造模及給藥：大鼠按隨機數字表分成3組：假手術組、模型組及續苓健骨方組，每組6只。根據參考文獻[1]造模。術后 30 d 開始灌胃給藥，續苓健骨方 15 g/kg，上午灌胃；假手術組和模型組給予生理鹽水 15 g/kg，連續給藥12周。各組動物均普通飼料飼養，自由飲水和活動。

1.2.2取材：最后一次灌胃完成2 h后，以10%水合氯醛麻醉后，取左側脛骨和第一、二腰椎。脛骨用于骨密度檢測；每組隨機選3只大鼠腰椎用于提取RNA及測序。

1.2.3骨密度檢測：采用雙能X線骨密度儀測定各組大鼠左側脛骨骨密度，操作過程嚴格按儀器使用說明進行。

1.2.4樣品RNA提取及質檢：從大鼠第一腰椎提取總 RNA，Nanodrop ND-1000測定RNA 260/280的濃度及蛋白質污染情況，采用變性瓊脂糖凝膠電泳檢測RNA純度及完整性。

1.2.5文庫構建及測序：構建好文庫，用Agilent 2100 Bioanalyzer進行文庫質量測定，混合好的不同樣品的測序文庫，通過0.1 M NaOH變性生成單鏈DNA，然后在Illumina NextSeq 500測序儀上進行51循環測序。文庫構建及測序均委托上海康成生物工程有限公司完成。

1.2.6測序數據分析：采用 miRDeep2 0.0.8 reads 軟件把trimmed miRNA 進行已知 miRNA定量和新miRNA預測。使用R軟件edgeR 3.18.1進行差異表達計算，篩選差異miRNAs，在線VENNY2.1分析共同差異miRNAs。利用miRNA靶基因數據庫統計top10差異miRNA的靶基因，并進行靶基因GO功能顯著性富集分析、靶基因Pathway顯著性富集分析。

1.3 統計學方法

使用edgeR篩選差異表達miRNA時，設定閾值為1.5倍差異，P<0.05且組內CPM均值≥1。采用 SPSS 17.0 軟件進行數據分析，計量資料用均數±標準差表示，根據數據是否正態分布，選擇非參數檢驗或t檢驗比較組間差異，P<0.05 表示差異有統計學意義。

2 結果

2.1 各組大鼠脛骨骨密度比較

與假手術組相比，模型組大鼠骨密度顯著下降(P<0.01)；續苓健骨方組骨密度較模型組顯著提高(P<0.01)；續苓健骨方組與假手術組相比無顯著性差異(P>0.05)，見表1。

表1 各組大鼠脛骨骨密度值比較

2.2 miRNA表達譜組間差異分析

各組之間差異miRNAs數見表2。VENNY 2.1在線進一步分析模型組與假手術組的差異表達基因在續苓健骨方組中表達情況。模型組表達上調，在續苓健骨方組表達下調的miRNAs有50條；模型組表達下調，在續苓健骨方組表達上調的miRNAs有60條，即模型組共有110條異常表達的miRNAs在續苓健骨方組得到糾正。表3列出了這110條差異miRNAs中，續苓健骨方治療之后上調或下調的top 10 miRNAs，它們在假手術組、模型組、續苓健骨方組的表達變化及其已知功能。

表2 組間差異miRNA數比較

表3 前10位差異表達 miRNAsTable 3 Top 10 of differential expressed miRNAs

2.3 實時熒光定量PCR驗證miRNAs表達

根據miRNA測序的結果，從表3中隨機選取了4個miRNA進行定量PCR驗證，內參選擇 snRNA U6。結果表明，rno-miR-136-3p、rno-miR-26b-5p在模型組中顯著提高(P<0.01)，在續苓健骨方組中顯著下降(P<0.01)；rno-miR-485-5p、rno-miR-138-5p在模型組中顯著降低(P<0.01)，在續苓健骨方組中顯著上升(P<0.01)。定量PCR的結果與測序結果一致(圖1)。

圖1 定量PCR驗證部分測序結果Fig.1 Partial sequencing results verified by real-time PCR注：與假手術組比較，*P<0.01；與模型組比較，?P<0.01

2.4 差異miRNAs靶基因預測及GO富集分析

基于miRNA數據庫，從110條miRNA中選擇經過續苓健骨方治療后上調或及下調的top 10 miRNAs，進行靶基因預測，然后對靶基因進行Gene Ontology(GO)分析，包括分子功能(molecular function，MF)、細胞組成(cellular component，CC)及參與的生物學過程(biological process，BP)。圖2分別列出了差異最顯著的前10個GO條目。這些靶基因參與了蛋白結合、多肽結合、轉錄因子調節、RNA聚合酶轉錄因子活性等過程。

圖2 GO分析柱形圖 A：表達上調miRNAs靶基因參與的GO分析；B：表達下調miRNAs靶基因參與的GO分析Fig.2 Column map of enrichment analysis. A: Go analysis of target genes involved in up regulated miRNAs; B: Go analysis of target genes involved in down regulated miRNAs.

2.5 差異miRNAs靶基因信號通路分析

根據預測的靶基因，用其在KEGG數據庫進行信號通路分析。分別選取前10個信號通路，根據P值繪制成柱狀圖(圖3)。表達上調的miRNA的靶基因參與的信號通路包括甘油磷脂代謝通路、磷脂酰肌醇信號通路、HIF-1信號通路、甲狀腺激素信號通路、Hedgehog信號系統、細胞內吞信號通路等；表達下調的miRNA的靶基因參與的信號通路主要包括內質網蛋白質加工信號通路、谷氨酸突觸信號通路、GnRH信號通路、自噬調節、CAMP及MAPK等信號通路等。這些信號通路很可能參與了續苓健骨方治療骨質疏松模型大鼠的治療過程。

2.6 預測新的miRNAs

通過使用Dicer對miRNA前體處理的簡單模型，miRDeep2預測新的miRNAs。表4列出了miRDeep2得分最高的10個新miRNAs。

圖3 KEGG信號通路分析 A：表達上調miRNAs靶基因參與的信號通路；B：表達下調miRNAs靶基因參與的信號通路Fig.3 Pathway analysis on KEGG. A: signaling pathways involved in up regulated miRNAs; B: signaling pathways involved in down regulated miRNAs

表4 新miRNA表達數據展示Table 4 Display of expression data of new miRNAs

3 討論

中醫藥治療骨質疏松的miRNA表達譜的研究較少。趙清等[20]應用Illumina高通量測序技術檢測中藥骨康治療絕經后骨質疏松癥(postmenopausal osteoporosis，PMOP)SD大鼠血漿miRNAs，認為差異miRNAs可能在PMOP發生發展以及中藥治療中起重要的調控作用。嚴芳娜[21]用miRNA芯片檢測黃精多糖對大鼠骨髓間充質干細胞誘導分化為成骨細胞和骨髓源性單核巨噬細胞誘導分化為破骨細胞過程中miRNAs表達譜的變化，認為miRNA的表達對誘導成骨細胞和破骨細胞均有重要作用。陳哲[22]運用高通量測序技術檢測并qPCR驗證滋腎降糖丸可能通過miRNAs對Wnt通路進行調控，進而影響成骨分化、實現防治糖尿病骨質疏松的作用。

本動物實驗中，續苓健骨方組大鼠骨密度則顯著提升。從差異miRNAs來看，續苓健骨方對在模型組出現異常表達的miRNAs起了較好的恢復作用。經12周的治療后，共有110條miRNAs的異常表達得到糾正。與假手術組相比，續苓健骨方組僅有8個異常表達的miRNAs，這表明無論是骨密度還是miRNA層面，續苓健骨方對骨質疏松癥均具有良好的治療作用。

在篩選到的差異miRNAs中，許多miRNAs均與骨質疏松癥或骨代謝密切相關，如Chen等[11]發現miRNA-136-3p在酒精誘導的小鼠骨量減少模型中的表達顯著降低，而且它能靶向PTEN調節血管生成和骨形成并改善酒精誘導的骨量減少。Panach等[23]從臨床患者血液中篩選到miR-122-5p，定量PCR證明其與骨質疏松性骨折相關；Mandourah等[13]采用芯片從臨床患者外周血篩選到差異表達miR-122-5p，定量PCR驗證它與低骨密度及脆性骨折相關，認為是一個潛在的診斷標記物；Liao等[24]也發現miR-122-5p負性調控SPRY2，升高受體酪氨酸激酶(RTK)活性，從而促進成骨細胞的增殖和分化。Zhang等[14]發現在骨質疏松癥患者中miR-485-5p表達水平上調，而在成骨分化過程中，miR-485-5p受到抑制，后又通過數據庫預測及熒光素酶報告基因證實miR-485-5p靶向Ostex促進骨質疏松。此外，本文在測序的基礎上，通過軟件預測得到了高信度的新miRNAs，不過其具體功能還有待進一步深入研究。

miRNA是機體重要的調控體系，對相關信號通路的調控是miRNA發揮作用的重要方式。本研究對top10差異miRNA的靶基因進行信號通路分析，發現靶基因參與HIF-1、甲狀腺激素、Hedgehog、GnRH、自噬、cAMP及MAPK等信號通路。這些信號通路均參與了骨代謝的過程。

綜上所述，本研究采用Illumina測序結合生物信息學分析，建立了卵巢切除大鼠模型骨組織miRNA表達譜，分析了續苓健骨方對骨質疏松模型大鼠腰椎miRNAs作用，篩選出差異表達miRNAs。top10 miRNAs靶基因參與了蛋白結合、多肽結合、轉錄調控、RNA聚合酶轉錄因子活性等生物學過程，及HIF-1、甲狀腺激素、Hedgehog、GnRH、自噬、cAMP及MAPK等信號通路，對于從miRNA角度探討骨質疏松癥的形成及中藥治療骨質疏松癥的分子機制具有重要參考價值。