梅花鹿鹿茸快速生長期環狀RNA鑒定及生物信息學1)

韓若冰 李和平

(東北林業大學，哈爾濱，150040)

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

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

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

1 材料與方法

1.1 樣品采集

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

1.2 RNA提取及文庫構建、測序

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

1.3 環狀RNA的鑒定與特征分析

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

1.4 表達量計算與差異分析

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

1.5 富集分析

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

1.6 circRNA-miRNA-mRNA調控網絡的構建

運用RNAhybrid (v2.1.2) 、svm_light (v6.01)，Miranda (v3.3)、TargetScan (v7.0)對最顯著的10個環狀RNA進行靶向miRNA的預測，結合本研究前期測得的miRNA數據[8]，從中篩選差異表達的miRNA。除此之外，利用mirTarBase預測得到與這些差異表達的miRNA具有靶向關系的mRNA，并從中篩選出差異表達的mRNA，進而得到miRNA-mRNA靶向作用關系，最終利用cytoscape[9]構建與鹿茸生長發育相關的circRNA-miRNA-mRNA調控網絡。

2 結果與分析

2.1 環狀RNA的鑒定

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

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

表1 環狀RNA的類型分布

表2 環狀RNA的長度分布

表3 環狀RNA的染色體分布

2.2 環狀RNA的差異表達分析

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

紅色表示表達量上調的環狀RNA；綠色表示表達量下調的環狀RNA；黑色表示表達量沒有差異。

表4 差異最顯著的前10個環狀RNA

2.3 差異表達基因的富集分析

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

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

環狀RNA來源基因的KEGG結果顯示共富集到81條信號通路上，包括粘著斑、PI3K-Akt信號通路、MAPK信號通路、TNF信號通路、FoxO信號通路(圖2)。其中，有13個差異顯著的環狀RNA的6個來源基因顯著富集到與細胞增殖、細胞周期等相關的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 差異表達環狀RNA的來源基因的KEGG富集分析

2.4 控網絡的構建

通過預測得到的環狀RNA、miRNA、mRNA之間的靶向互作關系構建最終的調控網絡，我們發現在10個差異最顯著的環狀RNA中，有2個環狀RNA預測得到的miRNA在鹿茸發育的前期和快速生長期的表達量差異不顯著。因此，最終的調控網絡由8個差異表達的環狀RNA和與其具有靶向互作關系的25個差異表達的miRNA所構成，通過進一步預測，得到這些差異表達的miRNA有120個存在靶向關系的差異表達的mRNA(圖3)。

圖3 差異表達的circRNA -miRNA-mRNA調控網絡

3 討論

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

本研究首次鑒定了鹿茸生長發育過程中的環狀RNA，得到鹿茸發育前期及快速生長期的間充質組織中共有7 790個環狀RNA，其中1 901個環狀RNA僅在鹿茸發育前期具有特異性表達，1 478個環狀RNA僅在鹿茸的快速生長期有特異性表達，有208個環狀RNA為顯著差異表達。我們猜測這些在不同生長發育時期有特異性表達的環狀RNA在鹿茸的某個生長發育時期發揮著獨特的調控作用。除此以外，我們發現一些環狀RNAs在鹿茸發育的2個時期均為顯著高表達(表達量高于其他環狀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)，這些高表達的環狀RNA通過與miRNA及mRNA之間的靶向作用關系對鹿茸的生長發育發揮著某些重要的調控作用，但具體的作用機制還需進一步的功能試驗來探究。

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