梅花鹿鹿茸快速生長期環(huán)狀RNA鑒定及生物信息學(xué)1)

韓若冰 李和平

(東北林業(yè)大學(xué)，哈爾濱，150040)

鹿茸的快速生長及其每年的周期性再生現(xiàn)象在哺乳動物中極為罕見。在梅花鹿鹿茸的快速生長期，其生長速度可達(dá)到12.5 mm·d-1[1]，不僅如此，鹿茸的生長速度堪比癌細(xì)胞，卻極少發(fā)生癌變，這些獨特的生物學(xué)特性使得鹿茸受到研究者的廣泛關(guān)注。了解鹿茸生長發(fā)育的分子調(diào)控機制對于鹿茸生物學(xué)的研究甚至是癌癥治療的進(jìn)一步研究均具有重要的科學(xué)意義。

迄今為止，已有眾多研究關(guān)注鹿茸生長發(fā)育的分子調(diào)控機制。隨著高通量測序技術(shù)的迅速發(fā)展，選用梅花鹿鹿茸作為試驗材料，從轉(zhuǎn)錄組水平上探究鹿茸生長發(fā)育的研究應(yīng)運而生。通過對不同生長時期的梅花鹿鹿茸進(jìn)行轉(zhuǎn)錄組測序，得到了幾種可能對鹿茸生長發(fā)揮重要調(diào)控作用的轉(zhuǎn)錄因子：PER1、EGFR1、GAS1[2]。通過對骨化期梅花鹿鹿茸的轉(zhuǎn)錄組測序得到了Sp7、SOX4、MMP9、MMP13等30種與鹿茸軟骨內(nèi)骨化相關(guān)的基因[3]。然而關(guān)于非編碼RNA在鹿茸生長發(fā)育過程中的研究還較為匱乏，主要集中于miRNA的研究，Yao et al.[4]通過高通量測序技術(shù)鑒定得到一些參與鹿茸軟骨發(fā)生和軟骨內(nèi)骨化的miRNA(miR-3072-5p、miR-1600、miR-34-5p、miR-6889-5p、miR-6729-5p)。

環(huán)狀RNA是一類在轉(zhuǎn)錄后水平具有重要基因表達(dá)調(diào)控作用的內(nèi)源性非編碼RNA[5-6]，不具有5’末端及3’末端。近些年來已經(jīng)成為1個新的研究熱點。然而，對于鹿科動物中環(huán)狀RNA的研究還較少。本研究運用高通量測序技術(shù)鑒定鹿茸2個生長發(fā)育時期(鹿茸生長前期、快速生長期)的環(huán)狀RNA并進(jìn)行生物信息學(xué)分析，試圖挖掘?qū)β谷咨L發(fā)育發(fā)揮調(diào)控作用的環(huán)狀RNA，旨在為梅花鹿的多組學(xué)研究及其生長發(fā)育機理的探究提供參考。

1 材料與方法

1.1 樣品采集

試驗材料采集自哈爾濱金地鹿業(yè)有限公司的3只健康成年雄性梅花鹿鹿茸的頂端間充質(zhì)組織。在其鹿茸生長至小鞍子茸型期，取左茸頂端組織的間充質(zhì)組織為鹿茸生長前期樣品；在其生長至二杠茸型時，采集右茸頂端組織的間充質(zhì)組織為快速生長期樣品。樣品采集后，將間充質(zhì)組織切成小塊，迅速放入液氮中保存?zhèn)溆谩?/p>

1.2 RNA提取及文庫構(gòu)建、測序

按照TRIzol試劑盒的操作說明書提取6個間充質(zhì)組織樣品的總RNA，使用分光光度計檢測提取總RNA的濃度、純度并使用Agilent Bioanalyzer檢測總RNA的完整性。提取得到的總RNA放置于冰箱-80 ℃保存以防止RNA降解。得到質(zhì)檢合格的RNA后，去除其中的核糖體RNA，合成cDNA的第一鏈和第二鏈并且進(jìn)行PCR擴增，利用Illumina HiseqTM 4000對構(gòu)建的文庫進(jìn)行測序。

1.3 環(huán)狀RNA的鑒定與特征分析

利用TopHat[7]將過濾后的數(shù)據(jù)與馬鹿參考基因組進(jìn)行比對，在比對結(jié)果中提取未比對上的數(shù)據(jù)，截取每一條未比對上的數(shù)據(jù)的兩端以得到短序列讀段，再次比對到參考基因組，并利用find_circ[5]軟件鑒定得到環(huán)狀RNA。鑒定條件為：只保留有且只有1個清晰斷點的環(huán)狀RNA；每條讀長的2個短序列讀段比對到基因組上的位置重疊不能超過2 bp且只允許2 bp錯配；獨特的讀長需大于2條；環(huán)狀RNA的長度小于100 kbp等。鑒定得到環(huán)狀RNA后，對環(huán)狀RNA類型、長度及環(huán)狀RNA在參考基因組上的分布進(jìn)行統(tǒng)計。

1.4 表達(dá)量計算與差異分析

采用每千個堿基的轉(zhuǎn)錄每百萬映射讀取的轉(zhuǎn)錄本(RPM)計算環(huán)狀RNA的表達(dá)量。根據(jù)得到的環(huán)狀RNA的表達(dá)量對鹿茸2個生長發(fā)育時期的環(huán)狀RNA進(jìn)行差異分析，從而得到顯著差異表達(dá)的環(huán)狀RNA，顯著差異表達(dá)的環(huán)狀RNA的篩選條件設(shè)置為P<0.05，|log2FC|>1。

1.5 富集分析

為探究環(huán)狀RNA在鹿茸生長發(fā)育過程中發(fā)揮的作用，對環(huán)狀RNA的來源基因進(jìn)行GO富集分析和KEGG功能注釋分類。在差異表達(dá)基因中顯著性富集的通路以及差異表達(dá)基因中顯著富集的GO條目篩選閾值均設(shè)置為q≤0.05。

1.6 circRNA-miRNA-mRNA調(diào)控網(wǎng)絡(luò)的構(gòu)建

運用RNAhybrid (v2.1.2) 、svm_light (v6.01)，Miranda (v3.3)、TargetScan (v7.0)對最顯著的10個環(huán)狀RNA進(jìn)行靶向miRNA的預(yù)測，結(jié)合本研究前期測得的miRNA數(shù)據(jù)[8]，從中篩選差異表達(dá)的miRNA。除此之外，利用mirTarBase預(yù)測得到與這些差異表達(dá)的miRNA具有靶向關(guān)系的mRNA，并從中篩選出差異表達(dá)的mRNA，進(jìn)而得到miRNA-mRNA靶向作用關(guān)系，最終利用cytoscape[9]構(gòu)建與鹿茸生長發(fā)育相關(guān)的circRNA-miRNA-mRNA調(diào)控網(wǎng)絡(luò)。

2 結(jié)果與分析

2.1 環(huán)狀RNA的鑒定

通過對測序得到的原始數(shù)據(jù)進(jìn)行過濾，得到的Q20為98.93%～99.05%，Q30為96.25%～96.58%，GC含量為51.35%～53.34%。從比對結(jié)果中提取未比對上的數(shù)據(jù)并截取其兩端，得到的短序列讀段為45865210～60425688，比對到參考基因組上的比對率為63.83%～67.32%。

在鹿茸發(fā)育的前期和快速生長期共鑒定得到7 790個環(huán)狀RNA。環(huán)狀RNA主要來源于蛋白質(zhì)編碼基因的已注釋的外顯子、單一外顯子、基因反義鏈、內(nèi)含子、外顯子和內(nèi)含子和基因間區(qū)[5，10]。對鑒定得到的環(huán)狀RNA類型進(jìn)行了統(tǒng)計(表1)，其中，大多數(shù)環(huán)狀RNA來自于蛋白質(zhì)編碼基因的基因間區(qū)及外顯子、內(nèi)含子區(qū)，來自于內(nèi)含子區(qū)的環(huán)狀RNA數(shù)目最少，而環(huán)狀RNA的長度主要分布于1～1 000之間(表2)。對環(huán)狀RNA在染色體上的分布進(jìn)行統(tǒng)計(表3)，結(jié)果表明7 790個環(huán)狀RNA分布在參考基因組的35條染色體上，其中分布在CM008042.1上的環(huán)狀RNA數(shù)量最少，而分布在CM008012.1上的環(huán)狀RNA數(shù)量最多。

表1 環(huán)狀RNA的類型分布

表2 環(huán)狀RNA的長度分布

表3 環(huán)狀RNA的染色體分布

2.2 環(huán)狀RNA的差異表達(dá)分析

通過對鑒定得到的2個發(fā)育時期的環(huán)狀RNA表達(dá)量進(jìn)行差異基因的篩選，共得到208個顯著差異表達(dá)的環(huán)狀RNA，包括顯著上調(diào)的39個差異表達(dá)環(huán)狀RNA以及顯著下調(diào)的169個差異表達(dá)環(huán)狀RNA(圖1)。差異最顯著的10個環(huán)狀RNA的詳細(xì)信息見表4。

紅色表示表達(dá)量上調(diào)的環(huán)狀RNA；綠色表示表達(dá)量下調(diào)的環(huán)狀RNA；黑色表示表達(dá)量沒有差異。

表4 差異最顯著的前10個環(huán)狀RNA

2.3 差異表達(dá)基因的富集分析

對差異環(huán)狀RNA的來源基因進(jìn)行富集分析可用于環(huán)狀RNA的功能預(yù)測。GO富集分析結(jié)果表明，環(huán)狀RNA的來源基因富集到166條與分子功能相關(guān)的GO條目上，顯著富集到生物調(diào)控，細(xì)胞過程、代謝過程等GO條目上；富集到156條與細(xì)胞組分相關(guān)的GO條目上，顯著富集到細(xì)胞、細(xì)胞部分等GO條目上；富集到980條與生物過程相關(guān)的GO條目上，包括結(jié)合與催化活性等相關(guān)的GO條目(表5)。

表5 差異表達(dá)的circRNA的來源基因的GO富集分析

環(huán)狀RNA來源基因的KEGG結(jié)果顯示共富集到81條信號通路上，包括粘著斑、PI3K-Akt信號通路、MAPK信號通路、TNF信號通路、FoxO信號通路(圖2)。其中，有13個差異顯著的環(huán)狀RNA的6個來源基因顯著富集到與細(xì)胞增殖、細(xì)胞周期等相關(guān)的PI3K-Akt信號通路上，分別是novel_circ_004050、novel_circ_004052、novel_circ_004059、novel_circ_004060、novel_circ_004062、novel_circ_004064、novel_circ_004069、novel_circ_004074、novel_circ_004863、novel_circ_002034、novel_circ_004646、novel_circ_003420、novel_circ_006430。

圖2 差異表達(dá)環(huán)狀RNA的來源基因的KEGG富集分析

2.4 控網(wǎng)絡(luò)的構(gòu)建

通過預(yù)測得到的環(huán)狀RNA、miRNA、mRNA之間的靶向互作關(guān)系構(gòu)建最終的調(diào)控網(wǎng)絡(luò)，我們發(fā)現(xiàn)在10個差異最顯著的環(huán)狀RNA中，有2個環(huán)狀RNA預(yù)測得到的miRNA在鹿茸發(fā)育的前期和快速生長期的表達(dá)量差異不顯著。因此，最終的調(diào)控網(wǎng)絡(luò)由8個差異表達(dá)的環(huán)狀RNA和與其具有靶向互作關(guān)系的25個差異表達(dá)的miRNA所構(gòu)成，通過進(jìn)一步預(yù)測，得到這些差異表達(dá)的miRNA有120個存在靶向關(guān)系的差異表達(dá)的mRNA(圖3)。

圖3 差異表達(dá)的circRNA -miRNA-mRNA調(diào)控網(wǎng)絡(luò)

3 討論

自從環(huán)狀RNA在植物類RNA病毒中首次被發(fā)現(xiàn)后[11]，在多種生物體內(nèi)均發(fā)現(xiàn)有環(huán)狀RNA的存在[12]，其表達(dá)具有時空特異性。根據(jù)競爭性內(nèi)源RNA機制(ceRNA)，環(huán)狀RNA通過與miRNA競爭性結(jié)合進(jìn)而調(diào)控下游靶基因的表達(dá)[13]，在多種腫瘤發(fā)生及組織器官的生長發(fā)育過程中均發(fā)揮著重要的調(diào)控作用，但在動物體中關(guān)于環(huán)狀RNA的研究還較少，且主要集中于環(huán)狀RNA的鑒定。

本研究首次鑒定了鹿茸生長發(fā)育過程中的環(huán)狀RNA，得到鹿茸發(fā)育前期及快速生長期的間充質(zhì)組織中共有7 790個環(huán)狀RNA，其中1 901個環(huán)狀RNA僅在鹿茸發(fā)育前期具有特異性表達(dá)，1 478個環(huán)狀RNA僅在鹿茸的快速生長期有特異性表達(dá)，有208個環(huán)狀RNA為顯著差異表達(dá)。我們猜測這些在不同生長發(fā)育時期有特異性表達(dá)的環(huán)狀RNA在鹿茸的某個生長發(fā)育時期發(fā)揮著獨特的調(diào)控作用。除此以外，我們發(fā)現(xiàn)一些環(huán)狀RNAs在鹿茸發(fā)育的2個時期均為顯著高表達(dá)(表達(dá)量高于其他環(huán)狀RNA 100倍以上，例如novel_circ_002441、novel_circ_007423、novel_circ_001328、novel_circ_007778、novel_circ_001447、novel_circ_005051、novel_circ_000884、novel_circ_005048、novel_circ_005975)，這些高表達(dá)的環(huán)狀RNA通過與miRNA及mRNA之間的靶向作用關(guān)系對鹿茸的生長發(fā)育發(fā)揮著某些重要的調(diào)控作用，但具體的作用機制還需進(jìn)一步的功能試驗來探究。

環(huán)狀RNA的功能分析主要通過其來源基因富集分析及與其具有靶向作用關(guān)系的miRNA和mRNA來進(jìn)行，顯著差異表達(dá)的環(huán)狀RNA來源基因的GO富集分析主要集中于生物調(diào)控、細(xì)胞過程、代謝過程、細(xì)胞代謝過程，而KEGG富集到代謝通路、粘著斑、PI3K-Akt信號通路、MAPK信號通路、FoxO信號通路上，說明差異顯著的環(huán)狀RNA主要涉及鹿茸發(fā)育的生物調(diào)控以及代謝等過程。在構(gòu)建的調(diào)控網(wǎng)絡(luò)中，1個環(huán)狀RNA可能會同時靶向多個miRNA，并且1個miRNA也可能會與多個環(huán)狀RNA具有靶向互作關(guān)系。其中，novel_circ_004646的靶向能力最強，靶向了13個差異表達(dá)的miRNA，而novel_circ_005407、novel_circ_000593只靶向了1個差異表達(dá)的miRNA。在與環(huán)狀RNA具有靶向關(guān)系的miRNA中，novel_circ_004646的其中1個靶向miR-457-x的靶向能力最強，與29個差異表達(dá)的mRNA均具有互作關(guān)系。但環(huán)狀RNA、miRNA、mRNA三者之間是否具有真實的靶向關(guān)系，其對鹿茸發(fā)育發(fā)揮著怎樣的調(diào)控作用以及如何發(fā)揮這些調(diào)控作用都需要進(jìn)一步的研究。