梅花鹿鹿茸全長轉錄組測序及鹿茸產量相關基因的挖掘

摘" 要： 旨在獲得梅花鹿鹿茸組織全長轉錄本和結構信息，并進一步挖掘出影響梅花鹿鹿茸產量的功能基因。本研究采取飼養管理條件一致的6頭2鋸齡健康雄性東大梅花鹿，取其相同部位鹿茸組織。根據茸重將鹿茸分為高、低產兩組，每組3個重復。利用PacBio Sequel平臺以及Illumina平臺分別對高、低產鹿茸混合樣品進行全長轉錄組測序。結合高、低產兩組二代測序技術，使用SUPPA、KOBOS-I、Cytoscape等工具對可變剪接、功能富集及Hub基因等預測分析。高、低產組測序分別獲得有效插入片段36 074 385 bp和32 413 962 bp；轉錄因子注釋顯示高產組中bHLH家族轉錄本數量要高于TF_bZIP家族，低產組中則相反。高產鹿茸組共存在6 644個可變剪接事件，而低產鹿茸組則為7 626個可變剪接事件；經比對，COL1A1在多個融合基因中起主效作用；兩組樣品中共篩選出207個差異表達基因，KEGG富集顯示差異表達基因主要集中在甲狀旁腺激素的合成、分泌以及膽固醇代謝等通路；鑒定到11個Hub基因，分別是APOB、SH3GL2、SYT1、ANGPTL4、FGF21、CACNA1E、SERPINF2、CACNB4、KLK4、LRP2、KCNA1。本研究結果為進一步的功能基因鑒定工作和分子育種策略的發展提供了新的理論基礎和研究路徑。

關鍵詞： 梅花鹿鹿茸；全長轉錄組測序；產茸性狀；功能基因；差異分析

中圖分類號：S813.1

文獻標志碼：A

文章編號：0366-6964（2024）12-5549-18

doi： 10.11843/j.issn.0366-6964.2024.12.020

開放科學（資源服務）標識碼（OSID）：

收稿日期：2024-03-14

基金項目：吉林省農業關鍵核心技術示范推廣（產業技術體系）項目（JARS-2024-0801-04）；吉林省畜禽遺傳資源開發利用及科技提質增效項目（202404-4）；吉林省自然科學基金（YDZJ202201ZYTS463）

作者簡介：閔祥玉（1998-），女，山東濰坊人，碩士生，主要從事動物繁育原理與技術研究，E-mail：mxy10504@163.com

*通信作者：王桂武，主要從事鹿遺傳育種研究，E-mail： wangguiwu2005@163.com；鄭軍軍，主要從事鹿遺傳育種研究，E-mail：zhengjunjun@caas.cn

Full-length Transcriptome Sequencing of Sika Deer Antler and Mining of Antler Yield-related Genes

MIN" Xiangyu， WEI" Jiali， XU" Biao， LIU" Huitao， ZHENG" Junjun*， WANG" Guiwu*

（Institute of Special Animal and Plant Sciences of Chinese Academy of Agricultural Sciences， Changchun 130112，" China）

Abstract： The study aimed to obtain the full-length transcript and structure information of the antler tissue of Sika deer， and to further mine the functional genes that affect the antler production of Sika deer. In this study， six two-sained healthy male Dongda Sika deer with the same feeding and management conditions were used to collect antler tissues from the same parts of the deer. According to the antler weight， the antler was divided into 2 groups of high and low production with 3 replicates per group. The whole length transcriptome of antler with high and low production were sequenced using PacBio Sequel platform and Illumina platform， respectively. Combined with the next generation sequencing technology of high and low yield groups， the alternative splicing， functional enrichment and Hub genes were predicted and analyzed by SUPPA， KOBOS-I and Cytoscape tools. The effective insertion fragments of 36 074 385 bp and 32 413 962 bp were obtained in the high and low yield groups， respectively. Transcription factor annotation showed that the number of bHLH family transcripts in the high yield group was higher than that of TF_bZIP family， but the opposite was observed in the low yield group. There were 6 644 alternative splicing events in the high production group and 7 626 alternative splicing events in the low production group. COL1A1 played a major role in multiple fusion genes. A total of 207 differentially expressed genes were screened out from the two groups of samples. KEGG enrichment showed that the differentially expressed genes were mainly concentrated in the synthesis and secretion of parathyroid hormone and cholesterol metabolism pathways. Eleven Hub genes were identified， including APOB， SH3GL2， SYT1， ANGPTL4， FGF21， CACNA1E， SERPINF2， CACNB4， KLK4， LRP2， KCNA1. This study provides a new theoretical basis and research path for further functional gene identification and the development of molecular breeding strategies.

Key words： sika deer antler; full-length transcriptome sequencing; anlter producing characters; functional gene; difference analysis

*Corresponding authors： WANG Guiwu， E-mail： wangguiwu2005@163.com; ZHENG Junjun， E-mail：zhengjunjun@caas.cn

梅花鹿是一種重要的特種畜禽，而鹿茸自古以來都是梅花鹿最重要的經濟性狀。鹿茸作為一種重要的動物藥材，具有抗疲勞、增強免疫力、壯腎陽以及改善部分婦科疾病等功效［1-4］。目前市場中，鹿茸仍處于供不應求的局面，所以通過測序等方法深入挖掘控制鹿茸產量的相關基因具有重要意義。

為深入挖掘梅花鹿產茸性狀相關基因，目前已進行了部分研究。梅花鹿和馬鹿鹿茸生長與血液IGF1濃度之間存在正相關性，GH及IGF1基因是促進鹿茸生長的主要調控因子［5］；COL1A1基因在快速生長期以及在間充質組織、前軟骨組織、軟骨組織中對鹿茸生長具有促進作用［6］；KGF基因在鹿茸生長及組織修復中扮演關鍵角色［7］；RXFP2基因作為雄激素受體的下游靶基因，參與角柄發生，是鹿茸再生生長的關鍵信號［8］。MTNR1A基因的多態性和突變位點在梅花鹿的生茸性狀中起著重要作用［9］；YWHAE基因及其編碼的14-3-3ε蛋白在梅花鹿生茸過程中可能扮演著重要角色［10］。Park等［11］鑒定到了130個蛋白，其中MAPK1、SMUBP2和ZFP28等可能與鹿茸生長密切相關。即使目前已經對梅花鹿鹿茸進行了相關研究，但并未完全闡明鹿茸生長機制，鹿茸產量調控機制仍需深入探究。

二代基因組測序、轉錄組測序等已經廣泛應用到鹿茸相關研究中［12-14］，但仍未對控制鹿茸產量的機制做出合理的解釋。可變剪接（alternative splicing）已經在牛、羊等畜禽中廣泛研究［15-18］，而在鹿科動物中鮮見報道。目前，三代轉錄組測序是檢測可變剪接事件最準確、有效的方法［19］，而主流的三代轉錄組測序平臺包括PacBio Sequel平臺和Nanopore平臺。在本研究中，利用PacBio Sequel及Illumina平臺進行混合測序，對梅花鹿高產茸以及低產茸的全長轉錄組進行對比，對兩組可變剪接事件進行統計分析。以期為深入理解梅花鹿產茸機制以及相關功能基因挖掘奠定基礎。

1" 材料與方法

1.1" 材料、試劑與儀器

梅花鹿二杠茸采購自長春市東大鹿業有限公司。

Trizol試劑、RNA 6000 Nano試劑盒、SMARTerTM PCR cDNA Synthesis試劑盒、TB Green Premix Ex TaqTM Ⅱ試劑盒、PrimeSTAR GXL DNA 聚合酶、SMRTbell Template Prep試劑盒、反轉錄試劑盒，其余試劑均為國產分析純。

電泳儀（Eppendorf公司）、熒光定量PCR儀（Thermo fisher scientific公司）、超微量分光光度計（天根生化有限公司）、凝膠成像儀（Bio-Rad公司）、Agilent 2100生物分析儀（Agilent公司）、Pacbio 測序平臺（Pacific Bio sciences 公司）、Illumina 測序平臺（Illumina公司）。

1.2" 方法

1.2.1" 樣本采集

試驗選取同一鹿圈、飼養管理條件一致、2鋸齡的6頭健康東大梅花鹿，根據茸重分為高產茸組（H）和低產茸組（L），樣品重量信息見表1。兩組樣品的重量在統計學上具有顯著差異（P≤0.01）（圖1）。在鹿茸生長50天時，取茸頭同一部位樣品切割后，PBS沖洗，并對樣品進行編號分裝，－80℃保存待測。

1.2.2" 總RNA提取

取高產組和低產組樣品分別置于1.5 mL無RNase離心管中，按照Trizol試劑說明書提取各樣本的總RNA，使用電泳儀、凝膠成像儀和生物分析儀評估RNA質量，合格者保留備用。

1.2.3" 三代全長轉錄組PacBio Sequel 測序

根據SMARTerTM PCR cDNA Synthesis試劑盒參照說明書對高產組及低產組RNA進行反轉錄合成cDNA，以cDNA為模板進行PCR擴增富集。隨后根據SMRTbell Tem-plate Prep試劑盒說明書對PCR擴增產物進行修復擴增純化，獲得測序文庫。測序工作主要由北京諾禾致源科技股份有限公司完成。

1.2.4" 二代高通量Illumina 測序

RNA樣品檢測合格后，用帶有Oligo（dT）的磁珠富集真核生物mRNA。隨后將mRNA打斷成短片段，以mRNA為模板，用六堿基隨機引物合成一鏈cDNA，然后合成二鏈cDNA，隨后利用AMPure XP beads純化雙鏈cDNA。純化的雙鏈cDNA先進行末端修復、加A尾并連接測序接頭，再用AMPure XP beads進行片段大小選擇。最后進行PCR擴增，并用AMPure XP beads純化PCR產物，得到最終的測序文庫。測序工作主要由北京諾禾致源科技股份有限公司完成。

1.2.5" 熒光定量PCR

從Hub基因中隨機找6個基因采用熒光定量PCR儀對兩組樣品進行熒光定量PCR驗證，每組3個重復，且每個樣品設置3個重復孔，根據NCBI中基因的序列設計引物（表2），以GAPDH作為內參基因進行矯正計算，結果采用2-ΔΔCT方法進行分析。

1.3" 數據處理

1.3.1" 全長轉錄組數據處理

采用PacBio官方軟件包SMRTlink處理PacBio測序原始下機數據，獲得高質量的consensus序列進行分析。利用獲得的Illumina高質量數據對三代reads進行測序錯誤校正，同時采用GMAP（genomic mapping and alignment program）軟件［20］將Polished consensus序列比對到梅花鹿參考基因組［21］上，并對未比對到參考基因組的轉錄本序列進行數據庫功能注釋；采用軟件SMRTlink v8.0對輸出數據進行過濾和處理；使用下機數據中subreads.bam文件經過CCS算法，即對單分子多次測序序列進行自我糾錯，得到CCS序列；通過檢測CCS是否包含5'-primer、3'-primer、poly-A對CCS進行分類找出FLNC序列和nFL序列；將同一轉錄本的FLNC序列使用hierarchical n*log（n）算法聚類，得到consensus序列；最后對得到的全長序列進行polish，終獲得Polished consensus序列進行后續分析。

1.3.2" 新基因的功能注釋

為了研究新基因的注釋信息，將新基因進行基因本體論（gene ontology，GO）、東京基因與基因組百科全書（kyoto encyclopedia of genes and genomes，KEGG）、非冗余數據庫（nonredundant protein database，NR）、蛋白質序列（swiss-port protein database，Swiss-Prot）、蛋白質真核同源（clusters of eukaryotic orthologous groups，KOG）、蛋白質家族數據庫（protein families database，PFAM）和核苷酸序列數據庫（nucleotide sequene database，NT）7個公共數據庫的功能注釋。

1.3.3" 可變剪接（AS）及融合基因分析

利用SUPPA軟件［22］對AS進行分析，將鑒定的AS事件分為7類：外顯子跳躍（SE）、外顯子互斥（MX）、5'端可變剪接（A5）、3'端可變剪接（A3）、內含子滯留（RI）、起始外顯子可變剪接（AF）、終止外顯子可變剪接（AL）。融合基因是指兩個或多個基因的編碼區首尾相連，置于同一套調控序列（包括啟動子、增強子、核糖體結合序列、終止子等）控制之下構成的嵌合基因。本研究將獲得的融合基因序列在NCBI數據庫（https：//www.ncbi.nlm.nih.gov/）中進行序列比對。

1.3.4" 差異表達基因篩選

利用DESeq軟件的readcount數據分析組間基因表達差異。用DEGseq軟件進行組間差異表達分析。P ≤0.05和|Fold change|gt;1作為差異表達基因的篩選閾值。

1.3.5" GO和KEGG功能富集分析

本研究使用KOBAS-I （http：//bioinfo.org/kobas/）在線軟件對差異表達基因進行GO、KEGG功能富集分析，P≤0.05 時，則表明該通路顯著富集。

1.3.6" 蛋白互作網絡分析（PPI）

在STRING數據庫中對所有差異表達基因進行蛋白互作分析，結果導入到Cytoscape3.7.1軟件中進行可視化。然后用CytoNCA插件進行Hub基因分析，參數為默認值。

1.3.7" 數據統計

采用excel進行數據錄入，采用SPSS 22.0統計學軟件進行數據分析，多組比較采用單因素方差分析，以P≤0.05表示統計學上差異顯著，柱形圖采用GraphPad Prism 8.0繪制。

2" 結" 果

2.1" 全長轉錄組數據統計

采用PacBio Sequel測序平臺分別對高、低產組梅花鹿鹿茸進行全長轉錄組測序，高產組共獲得了148.46 GB的Subreads，有效插入片段數目36 074 385 bp，平均長度為4 116 bp，N50為4 548 bp；低產組共獲得了135.06 GB的Subreads，有效插入片段數目32 413 962 bp，平均長度4 167 bp，N50為4 534 bp。

2.2" 轉錄本數據庫注釋

高產組鹿茸共檢測到16 133條全長轉錄本，其中已知轉錄本數目1 620條，已知基因新轉錄本數目13 834條，新基因的新轉錄本數目679條，低產組鹿茸共檢測到19 742條全長轉錄本，其中已知轉錄本數目1 790條，已知基因新轉錄本數目17 001條，新基因的新轉錄本數目951條。本研究對對應物種信息的NR數據庫、蛋白數據庫KOG注釋信息、劃分基因功能的GO數據庫以及分析基因通路的KEGG數據庫進一步分析。

2.2.1" NR功能注釋

在本研究中，高產組中共有15 895條Unigenes被注釋到NR數據庫中，其中有12.61%的新基因被注釋到綿羊（Ovis aries）上、12.02%的新基因被注釋到黃牛（Bos taurus）上、11.01%的新基因被注釋到山羊（Capra hircus）上。低產組中共有19 380條Unigenes被注釋到NR數據庫中，其中有12.33%的新基因被注釋到黃牛（Bos taurus）上、12.00%的新基因被注釋到綿羊（Ovis aries）上、11.23%的新基因被注釋到山羊（Capra hircus）上。對比可知，高、低產組在物種相似度上差異不顯著（圖2）。

2.2.2" KOG功能注釋

高產組在KOG數據庫中共有12 561條Unigenes被注釋，分別注釋到26個小類中，其中被注釋到的Unigenes最多的是一般功能小類（2 084條），隨后是信號傳導機制（1 907條）以及細胞外結構（1 289條），而未命名蛋白被注釋到的Unigenes最少（11條），其次是次級代謝產物的生物合成、轉運和分解代謝（47條）；而低產組共有15 268條Unigenes被注釋，同樣注釋到26個小類中，被注釋到較多的小類同高產組相同分別為一般功能小類（2 642條）、信號傳導機制（2 421條）以及細胞外結構（1 421條），被注釋到較少的與高產組略有不同，分別是未命名蛋白（17條）以及細胞運動性功能（47條）（圖3）。

2.2.3" KEGG功能注釋

高、低產組在KEGG數據庫中分別注釋到15 757條Unigenes和19 188條Unigenes，其中均包含新陳代謝、基因信息處理、環境信息處理、細胞進程、生物系統以及人類疾病6個一級通路。在兩組數據中均為信號傳導通路占比最高分別為2 443條和2 866條，耐藥性、抗菌性占比最低，均為1條（圖4）。

2.2.4" GO功能注釋

在本研究中，高產組在GO數據庫中被注釋到的Unigenes為12 180條，而低產組為14 519條，均可分為生物過程、細胞組分、分子功能3個大類。在高產組生物過程的24個小類中以細胞過程（5 110條）定位的較多，以行為（1條）、定位的較少；在細胞組分的18個小類中以細胞（3 348條）定位的最多，以突觸（11條）定位的最少；在分子功能的11個小類中以結合活性（8 346條）定位的居多，以金屬合酶活性（1條）、定位的較少。低產組生物過程的23個小類中以細胞過程（6 156條）、定位的較多，以節律過程（6條）定位的較少；在分子功能的10個小類以結合活性（9 904條）定位的居多，以抗氧化活性（58條）定位的較少；在細胞組分的18個小類中以細胞組成（3 994條）定位的較多，以突觸（10條）定位的最少（圖5）。

2.3" 轉錄因子注釋分析

根據注釋結果，將注釋到轉錄本數目最多的前30個轉錄因子家族進行記錄。其中高產組注釋較多的家族有zf-C2H2（262個）、ZBTB（110個），而zf-GATA、HSF、IRF、CG-1家族注釋較少，均只有3個。低產組注釋較多的家族有zf-C2H2（358個）、ZBTB（114個），而CG-1、NF-YA、ESR-like家族注釋較少，均為4個。值得一提的是高產組中bHLH家族注釋數量（39個）要高于TF_bZIP家族（38個），而低產組中則相反，TF_bZIP家族注釋數量（46個）高于bHLH家族（40個）（圖6）。

2.4" AS分析

在發生可變剪接的基因比例中，高產組SE、AF占比較高，分別為20.48% 和15.13%，MX占比最低，為1.85%。低產組占比最高的同樣為SE和AF占比較高20.47% 和15.26%，最低為MX，占比為1.63%（表3）。

2.5" 融合基因分析

在融合基因檢測中，高產組共檢測到26個融合基因，其中7個融合基因來自同一個染色體，占比約為26.92%，19個來自不同染色體，占比約為73.08%；低產組共檢測到28個融合基因，其中9個融合基因來自同一個染色體，占比約為32.14%， 19個來自不同染色體，占比約為67.86%。通過對兩組融合基因的對比，發現即使是不同染色體來源的融合基因中兩組仍有相似，但是高產組中有5個融合基因為其獨有，而低產組中有7個融合基因為其獨有，通過進一步對這12個融合基因進行數據庫比對發現，COL1A1基因在3個融合基因中均起到主要作用。

2.6" 差異表達基因分析

將每個樣品測序數據經過DESeq標準化后進行差異分析，以Plt;0.05和|Fold change|gt;1為篩選閾值，兩組樣品中共篩選到207個差異表達基因，其中59個上調基因，148個下調基因（圖7）。

為了深入了解高、低產鹿茸中差異基因的表達模式，本研究對篩選出的差異表達基因進行了聚類分析。KEGG富集結果顯示：差異表達基因主要與甲狀旁腺激素的合成、分泌和作用相關通路（parathyroid hormone synthesis， secretion and action）、鈣信號通路（calcium signaling pathway）、膽固醇代謝通路（cholesterol metabolism）、精、脯氨酸代謝（arginine and proline metabolism）等通路相關（圖8）。GO富集結果顯示：差異表達基因主要與生物生長的負調控（negative regulation of multicellular organism growth）、細胞外空間（extracellular space）、鈣離子調節的胞吐作用（calcium-ion regulated exocytosis）等條目相關（圖9）。

對篩選出的207個差異表達基因進行PPI蛋白互作分析（圖10）。以BC（betweenness centrality）值為主要篩選條件，篩選出11個Hub基因，依次是APOB、SH3GL2、SYT1、ANGPTL4、FGF21、CACNA1E、SERPINF2、CACNB4、KLK4、LRP2、KCNA1。隨機選取APOB、SYT1、ANGPTL4、FGF21、SERPINF2、CACNB4六個基因對其進行熒光定量PCR驗證。結果與測序結果相吻合（圖11）。

3" 討" 論

長期以來，鹿茸一直是梅花鹿生產過程中的首要經濟性狀［23］。為了獲得更加準確的梅花鹿鹿茸轉錄組數據，以及對高、低產鹿茸轉錄組進行對比，深入挖掘鹿茸生長相關基因。在本研究中，利用PacBio Sequel測序以及Illumina測序相結合的方法，共獲得了高、低產鹿茸約283.52 Gb的數據。同

內圈為11個Hub基因

The 11 Hub genes in the inner circle

時，高、低產組分別共檢測到16 133、19 742條全長轉錄本，在數據庫注釋中，兩組轉錄組在KOG數據庫存在明顯差異。測序數據經DESeq標準化后進行差異分析共得到207個差異表達基因，這些基因可能與鹿茸產量有關。

Zhang等［24］將生長階段以及骨化階段的鹿茸分別進行測序，以期深入了解鹿茸快速生長和骨化的相關機制。其中，該研究對7大數據庫注釋顯示NT數據庫占比最高，其次為NR數據庫；在NR數據庫和KEGG數據庫中兩組差異并不明顯。但在KOG數據庫中，在次級代謝產物的生物合成、轉運和分解代謝小類中的注釋略有差別，可能是因為某些代謝產物影響著鹿茸生長，但仍需進一步試驗驗證；GO數據庫注釋顯示，高、低產兩組在結合活性注釋到的基因均最多，這與在小尾寒羊中的結果一致［25］，這說明結合活性在影響鹿茸產量的過程中可能發揮了重要作用。

本研究中，高、低鹿茸產組分別注釋到718、907個轉錄因子。其中兩組中轉錄因子注釋最多的家族均為zf-C2H2，該家族是人類最大的一個轉錄因子家族，也是哺乳動物最大的家族之一［26］。其次均為ZBTB家族，該家族具有基因轉錄功能，對于細胞的增殖、分化、凋亡和功能調節等方面發揮著重要的作用［27］。而兩組在bHLH與TF_bZIP兩個家族的注釋中略有差別，高產組中bHLH家族數量要高于TF_bZIP家族，低產組中則相反。在bHLH轉錄因子家族中，有16個成員僅在高產組中被發現，而另外17個成員則僅存在于低產組中。這些轉錄因子很有可能與鹿茸生長相關。目前，在鹿科動物中，并沒有對該類轉錄因子進行深入研究，但在牦牛中［28］，bHLH轉錄因子參與轉錄調節、DNA結合、神經元分化以及肌肉器官發育等多個生物學過程。而在山羊中［29］，該家族廣泛參與肌肉細胞生成與器官發生、細胞增殖與分化以及晝夜節律等相關生長發育過程。由此可知，該家族或許與發育和增殖有關，因此可為研究梅花鹿鹿茸生長提供重要依據。

本研究兩組樣本中檢測到的可變剪接事件均為SE占比最高，這與山羊［30］、伊拉兔［31］等的分析結果一致，這或許可以說明SE事件在動物體中最為普遍，發揮的功能更加廣泛。車隆等［32］的研究發現，SE可能與骨骼肌細胞的分化、橫紋肌組織發育調節等密切相關，這些數據說明SE可能與分化、發育調節等生物過程相關，所以可能在鹿茸生長中發揮作用。此外，兩組發生可變剪接時間的比例雖然相似，但仍存在較多差異基因，這些發生相同可變剪接事件的不同基因很可能與鹿茸生長相關。

在融合基因檢測中，低產組中不同染色體上融合的基因相較于高產組來說占比更高，更容易引起基因表達異常以及蛋白質結構的改變。同時在不同染色體來源的融合基因中，高產組中有5個融合基因為其獨有，而低產組中有7個融合基因為其獨有。通過對兩組特有的融合基因序列進行數據庫比對發現，COL1A1基因在3個融合基因均發揮主效作用。COL1A1基因編碼的I型膠原蛋白α1鏈是細胞外基質（ECM）的主要結構蛋白，在多個組織中發揮作用。該基因不僅作為結構蛋白起到物理支撐作用，還作為關鍵因子參與多種細胞信號傳導。COL1A1通過與細胞表面的整合素受體結合，可以激活多種信號傳導途徑，如MAPK/ERK、PI3K/Akt等，這些信號途徑參與調節細胞的增殖、分化和生存［33］。在腫瘤細胞相關的研究中，過表達COL1A1基因，能夠促進腫瘤細胞的增殖，并影響其侵襲［34］。COL1A1還可能通過影響細胞凋亡相關信號途徑，參與調節細胞的生存和死亡。綜上所述，COL1A1基因可能通過影響細胞增殖、分化和凋亡等過程，來對鹿茸生長發揮作用，該基因可以作為影響茸重性狀的候選基因。融合基因可以加速遺傳物質的重組和傳遞，從而縮短育種周期，為梅花鹿育種工作提供新思路。

GO和KEGG富集分析顯示，與鹿茸產量相關的差異基因主要與鈣離子調節的胞吐作用、甲狀旁腺激素的合成、分泌和作用相關通路、膽固醇代謝通路以及精、賴氨酸通路等生物過程相關。目前，甲狀腺激素、甲狀旁腺激素已被證實可以與其他激素共同調控骨骼的生長［35］。由膽固醇轉化而成的1，25-二羥基維生素D3能夠促進破骨細胞分化，調節鈣的重吸收［36］，研究證實，破骨細胞對于鹿茸生長有著密切關系。鹿茸生長過程中鈣離子調節和膽固醇代謝發揮著不可或缺的作用，這與Van Der Eems等［37］的研究結果一致。張然然等［38］提出精氨酸在鹿茸生長發育過程中會有顯著變化。以上結果表明，甲狀旁腺激素的合成、分泌和作用相關通路、膽固醇代謝通路和精賴氨酸代謝通路等可以作為影響鹿茸產量的關鍵通路。

2018年，Hu等［39］采用PacBio公司的RS系統對梅花鹿三杈茸高產、低產兩組混合RNA進行全長轉錄組測序，同時高產、低產兩組共6個樣品進行Illumina測序。獲得了ACTN1、BAG6、COL1A2等26個與鹿茸生長相關的基因，與本研究并沒有重合基因，這可能是由于基因表達在鹿茸不同生長期會存在差異而二杠茸與三杈茸生長期不同，發揮作用的基因不同所導致的。本研究對差異基因進行蛋白互作分析獲得了11個Hub基因。其中，LRP2作為一種多配體內吞受體，在脂質轉運中起著至關重要的作用，它通過介導多種脂質成分的細胞攝取和轉運，維持脂質穩態［40］。APOB是血液中脂質運輸的關鍵蛋白，同時可以參與脂蛋白的代謝過程。吉琳等［41］研究結果同樣表明APOB在脂類代謝中發揮重要作用且鹿茸蠟片中脂質代謝產物含量最高［42］。ANGPTL4基因可以通過影響甘油三酯和高密度脂蛋白的水平，參與調節脂質代謝［43］。FGF21可以通過促進脂肪酸氧化、抑制脂肪酸合成、調節脂蛋白代謝和抗動脈粥樣硬化作用，在脂質代謝中發揮著重要的調節功能［44］。LRP2、APOB、ANGPTL4和FGF21通過不同的機制協同工作，維持脂質穩態，并影響脂質的攝取、轉運、合成和代謝過程。在本研究中，這些基因共同富集到膽固醇代謝通路，膽固醇代謝可以調節IGF1，IGF1刺激膽固醇合成，故這些基因可能通過調節脂質代謝來影響鹿茸生長。KLK4可以通過釋放細胞外基質中的生長因子和細胞因子，從而影響細胞生長、分化和遷移［45］。CACNB4基因編碼的CaMKIIβ亞基是細胞內鈣離子信號通路的關鍵調節因子，通過感應和響應鈣離子濃度變化，激活并參與調控細胞內多種生理過程和信號傳導［46］。KLK4與CACNB4可能是通過調節鈣離子相關功能在鹿茸生長中發揮作用。鹿茸的再生是包括神經、血管等組織的完全再生［47］，SYT1基因編碼的突觸結合蛋白1作為一種鈣離子感受器，在神經遞質釋放和突觸小泡停泊過程中發揮關鍵調控作用［48］。SH3GL2在神經系統中廣泛存在［49］，因此，SH3GL2和SYT1可能通過影響神經組織調節鹿茸的再生。通過進一步比對發現，KEGG富集結果與11個Hub基因具有一致性。主要富集在甲狀旁腺激素的合成、分泌和作用相關通路、鈣離子相關通路以及膽固醇代謝通路等相關通路。

4" 結" 論

本研究基于三代全長轉錄組以及二代高通量混合測序，更加準確地挖掘出與梅花鹿茸重相關的bHLH轉錄因子家族、SE可變剪接、融合基因中的主效基因COL1A1以及207個候選差異基因。11個樞紐基因均富集在甲狀旁腺激素、鈣離子調節以及膽固醇代謝等通路，可能通過與IGF1互作調控鹿茸的生長。但與IGF1的互作機制還需進一步探究。本研究為探索茸重基因提供了新的見解，為后續的功能基因驗證和分子育種指明了方向。

參考文獻（References）：

［1］" LI C Y，YANG F H，SHEPPARD A.Adult stem cells and mammalian epimorphic regeneration-insights from studying annual renewal of deer antlers［J］.Curr Stem Cell Res Ther，2009，4（3）：237-251.

［2］" LI C Y，LI Y，WANG W Y，et al.Deer antlers：the fastest growing tissue with least cancer occurrence［J］.Cell Death Differ，2023，30（12）：2452-2461.

［3］" ZHANG P，GUO Z H，MA L，et al.Investigation of anti-fatigue effect and simultaneous determination of eight nucleosides in different parts of velvet antler in red deer and sika deer［J］.Chem Biodivers，2020，17（2）：e1900512.

［4］" 胡艷紅，顏" 鑫，雷" 燕，等.鹿茸的化學成分、藥理作用與臨床應用研究進展［J］.遼寧中醫藥大學學報，2021，23（9）：47-52.

HU Y H，YAN X，LEI Y，et al.Research progress on chemical constituents，pharmacological effects and clinical application of Lurong （Cervi Cornu Pantotrichum）［J］.Journal of Liaoning University of TCM，2021，23（9）：47-52.（in Chinese）

［5］" 李光玉.鹿茸生長與血液IGF-1濃度年周期變化規律［D］.長春：中國農業科學院特產研究所，2003.

LI G Y.Changing law of annual cycle of antler growth and blood IGF-1 concentration［D］.Changchun：Institute of Special Animal and Plant Sciences of CAAS，2003.（in Chinese）

［6］" 張芙蕊，韓若冰，郭夢雅，等.梅花鹿茸生長過程中頂端不同組織COL1A1基因啟動子DNA甲基化模式及差異分析［J］.中國農業大學學報，2022，27（10）：153-163.

ZHANG F R，HAN R B，GUO M Y，et al.DNA methylation and differential analysis of COL1A1 gene promoter in the different tissues of antler top during the antler growing in sika deer （Cervus nippon）［J］.Journal of China Agricultural University，2022，27（10）：153-163.（in Chinese）

［7］" 郭夢雅，韓若冰，邢海華，等.梅花鹿KGF的基因克隆及在不同時期頂端茸皮組織的表達分析［J］.中國農業大學學報，2022，27（6）：154-161.

GUO M Y，HAN R B，XING H H，et al.Cloning of KGF gene from sika deer and expression analysis of antler tip skin tissues at different stages［J］.Journal of China Agricultural University，2022，27（6）：154-161.（in Chinese）

［8］" BA H X，WANG X，WANG D T，et al.Single-cell transcriptome reveals core cell populations and androgen-RXFP2 axis involved in deer antler full regeneration［J］.Cell Regen，2022，11（1）：43.

［9］" 熊麗容.梅花鹿MTNR1A基因多態性及其與產茸量的相關性研究［D］.武漢：華中農業大學，2011.

XIONG L R.The study of MTNR1A gene in the influence of sika deer and its polymorphism［D］.Wuhan：Huazhong Agricultural University，2011.（in Chinese）

［10］" 劉" 慧，王桂武，王權威，等.梅花鹿YWHAE基因編碼區全長cDNA的克隆及序列分析［J］.農業生物技術學報，2017，25（5）：833-841.

LIU H，WANG G W，WANG Q W，et al.cDNA coding region cloning and sequence analysing of sika deer （Cervus nippon） YWHAE Gene［J］.Journal of Agricultural Biotechnology，2017，25（5）：833-841.（in Chinese）

［11］" PARK H J，LEE D H，PARK S G，et al.Proteome analysis of red deer antlers［J］.Proteomics，2004，4（11）：3642-3653.

［12］" 鞠" 妍，劉華淼，何金明，等.基于高通量測序技術的敖魯古雅馴鹿鹿茸轉錄組分析［J］.中國畜牧獸醫，2020，47（7）：1981-1989.

JU Y，LIU H M，HE J M，et al.Transcriptome analysis of antler in Aoluguya reindeer based on high-throughput sequencing technology［J］.China Animal Husbandry amp; Veterinary Medicine，2020，47（7）：1981-1989.（in Chinese）

［13］" 楊曉光.馬鹿（Cervus elaphus）鹿茸角頂端茸皮與軟骨組織轉錄組研究［D］.哈爾濱：東北林業大學，2015.

YANG X G.Velvet skin and cartilage transcriptomes at the antler tips of Red deer Cervus elaphus［D］.Harbin：Northeast Forestry University，2015.（in Chinese）

［14］" 幺寶金.梅花鹿鹿茸頂端組織轉錄組分析及不同生長期差異基因篩選［D］.長春：吉林大學，2012.

YAO B J.Transcriptome analysis and differentially expressed genes screening of antler tips at different growth stages［D］.Changchun：Jilin University，2012.（in Chinese）

［15］" 崔迎迎，芒" 來，李" 蓓，等.MEI1基因可變剪切事件對蒙古馬精子生成的調控作用［J］.畜牧獸醫學報，2022，53（4）：1096-1108.

CUI Y Y，MANG L，LI B，et al.Regulatory role of MEI1 with alternative splicing event on spermatogenesis in mongolian horse［J］.Acta Veterinaria et Zootechnica Sinica，2022，53（4）：1096-1108.（in Chinese）

［16］" 羅朋娜.雞GLUT4基因可變剪切體鑒定與主要轉錄本的功能分析［D］.鄭州：河南農業大學，2023.

LUO P N.Identification of alternative splicing variants and functional analysis of major transcripts of GUT4 gene in chicken［D］.Zhengzhou：Henan Agricultural University，2023.（in Chinese）

［17］" 王進娥，麻俊淵，陳" 舟，等.牦牛發情期與妊娠期轉錄本及剪切事件研究［J］.甘肅農業大學學報，2022，57（6）：1-9.

WANG J E，MA J Y，CHEN Z，et al.Transcription and splicing events during estrus and pregnancy in the yak［J］.Journal of Gansu Agricultural University，2022，57（6）：1-9.（in Chinese）

［18］" 陳" 媛，蔡" 禾，李" 利，等.山羊TNNT3基因可變剪切及其對骨骼肌細胞分化的作用［J］.中國農業科學，2021，54（20）：4466-4477.

CHEN Y，CAI H，LI L，et al.Alternative splicing of TNNT3 and its effect on the differentiation of Mu SCs in goat［J］.Scientia Agricultura Sinica，2021，54（20）：4466-4477.（in Chinese）

［19］" HESTAND M S，AMEUR A.The versatility of SMRT sequencing［J］.Genes （Basel），2019，10（1）：24.

［20］" GASPERINI M，HILL A J，MCFALINE-FIGUEROA J L，et al.A genome-wide framework for mapping gene regulation via cellular genetic screens［J］.Cell，2019，176（1-2）：377-390.e19.

［21］" XING X M，AI C，WANG T J，et al.The first high-quality reference genome of sika deer provides insights into high-tannin adaptation［J］.Genom Proteom Bioinform，2023，21（1）：203-215.

［22］" ALAMANCOS G P，PAGS A，TRINCADO J L，et al.Leveraging transcript quantification for fast computation of alternative splicing profiles［J］.Proc Natl Acad Sci U S A，2015，21（9）：1521-1531.

［23］" 李浩東，閔祥玉，周" 雅，等.基于GBLUP等模型對梅花鹿（Cervus Nippon）生長相關性狀基因組選擇的預測準確性比較［J］.畜牧獸醫學報，2023，54（2）：608-616.

LI H D，MIN X Y，ZHOU Y，et al.Comparison of prediction accuracy of genomic selection for growth-related traits in sika deer （Cervus Nippon） based on GBLUP and other models［J］.Acta Veterinaria et Zootechnica Sinica，2023，54（2）：608-616.（in Chinese）

［24］" ZHANG R R，DONG Y M，XING X M.Comprehensive transcriptome analysis of sika deer antler using PacBio and Illumina sequencing［J］.Sci Rep，2022，12（1）：16161.

［25］" 張春蘭.小尾寒羊和杜泊羊臂二頭肌轉錄組及肌球蛋白輕鏈基因家族結構特征分析［D］.泰安：山東農業大學，2014.

ZHANG C L.Transcriptome analysis of small-tailed Han sheep and Dorper's biceps Brachii and structure characteristics of myosin light chain gene families［D］.Taian：Shandong Agricultural University，2014.（in Chinese）

［26］" MYRZABEKOVA M O，LABEIT S B，NIYAZOVA R Y.Features of miRNAs binding sites within the C2H2 ZNF family：a Bos taurus，Eqqus caballus，and Ovies aries comparative approach［J］.Int J Biol Chem，2020，13（1）：33-46.

［27］" 林" 潔.肝臟ZBTB20對膽固醇代謝的調節作用及其機制研究［D］.上海：中國人民解放軍海軍軍醫大學，2023.

LIN J.Study on the regulation of cholesterol metabolism by liver ZBTB20［D］.Shanghai：Naval Medical University，2023.（in Chinese）

［28］" 李楓梅.肉用牛bHLH轉錄因子家族成員的鑒定與分析［D］.阜陽：阜陽師范學院，2017.

LI F M.Identification and functional analysis of the basic helix-loop-helix （bHLH） transcription factors in the meat cattle［D］.Fuyang：Fuyang Normal University，2017.（in Chinese）

［29］" 魏振宇.山羊bHLH和PLAG基因家族分析及其與生長相關遺傳標記挖掘［D］.楊凌：西北農林科技大學，2022.

WEI Z Y.Analysis of goat bHLH and PLAG gene families and mining of genetic markers associated with growth［D］.Yangling：Northwest Aamp;F University，2022.（in Chinese）

［30］" 鮑晶晶，浦亞斌，馬月輝，等.綿羊不同發育階段背最長肌組織中可變剪接的鑒定與分析［J］.生物技術通報，2019，35（7）：32-38.

BAO J J，PU Y B，MA Y H，et al.Identification and analysis of alternative splicing in Longissimus dorsi of sheep at different development stages［J］.Biotechnology Bulletin，2019，35（7）：32-38.（in Chinese）

［31］" 吳倩雯，王" 露，王文佳，等.不同生長階段伊拉兔骨骼肌轉錄組數據的挖掘分析［J］.基因組學與應用生物學，2024，43（6）：1028-1038.

WU Q W，WANG L，WANG W J，et al.Mining analysis of skeletal muscle transcriptome data of Ira rabbits at different growth stages［J］.Genomics and Applied Biology，2024，43（6）：1028-1038.（in Chinese）

［32］" 車" 隆，葉" 琳，吳鵬飛，等.京海黃雞生長早期骨骼肌可變剪切分析［J］.中國家禽，2022，44（12）：21-27.

CHE L，YE L，WU P F，et al.Alternative splicing analysis of skeletal muscle for Jinghai Yellow chicken during early growth［J］.China Poultry，2022，44（12）：21-27.（in Chinese）

［33］" 周明帥，溫曉艷，張" 艷，等.COL1A1基因在貴州黑山羊性腺軸中的表達及其對產羔相關基因的影響［J］.農業生物技術學報，2022，30（11）：2152-2162.

ZHOU M S，WEN X Y，ZHANG Y，et al.Expression of COL1A1 gene in the gonad axis of Guizhou black goat （Capra hircus） and its effect on lambing related genes［J］.Journal of Agricultural Biotechnology，2022，30（11）：2152-2162.（in Chinese）

［34］" LI S J，MEI L Y，HE C F，et al.Identification of a family with van der Hoeve's syndrome harboring a novel COL1A1 mutation and generation of patient-derived iPSC lines and CRISPR/Cas9-corrected isogenic iPSCs［J］.Human Cell，2024，37（3）：817-831.

［35］" KRONENBERG H M.Developmental regulation of the growth plate［J］.Nature，2003，423（6937）：332-336.

［36］" 周" 琳，崔" 磊，鄧" 敦.母豬肢蹄病的病因、危害與預防措施［J］.豬業科學，2023，40（9）：61-62.

ZHOU L，CUI L，DENG D.Causes，hazards and preventive measures of hoof and limb disease in sows［J］.Swine Industry Science，2023，40（9）：61-62.（in Chinese）

［37］" VAN DER EEMS K L，BROWN R D，GUNDBERG C M.Circulating levels of 1，25 dihydroxyvitamin D，alkaline phosphatase，hydroxyproline，and osteocalcin associated with antler growth in white-tailed deer［J］.Acta Endocrinol （Copenh），1988，118（3）：407-414.

［38］" 張然然，榮" 敏，董依萌，等.不同生長時期梅花鹿鹿茸代謝組分析［J］.畜牧獸醫學報，2022，53（12）：4518-4526.

ZHANG R R，RONG M，DONG Y M，et al.Metabolomic analysis of sika deer antler in different growth stages［J］.Acta Veterinaria et Zootechnica Sinica，2022，53（12）：4518-4526.（in Chinese）

［39］" HU P F，WANG T J，LIU H M，et al.Full-length transcriptome and microRNA sequencing reveal the specific gene-regulation network of velvet antler in sika deer with extremely different velvet antler weight［J］.Mol Genet Genomics，2019， 294（2）：431-443.

［40］" ZHAO Y Y，LIU W M，ZHAO X M，et al.Low-density lipoprotein receptor-related protein 2 （LRP2） is required for lipid export in the midgut of the migratory locust，Locusta migratoria［J］.J Integr Agric，2024，23（5）：1618-1633.

［41］" 吉" 琳，楊秋月，方斐旻，等.雞Apob基因組織表達與CRISPR/Cas9敲除系統的構建［J］.畜牧獸醫學報，2021，52（3）：630-640.

JI L，YANG Q Y，FANG F M，et al.Tissue expression and construction of CRISPR/Cas9 knockout system of Apob in chicken［J］.Acta Veterinaria et Zootechnica Sinica，2021，52（3）：630-640.（in Chinese）

［42］" 劉文媛.基于蛋白質組學和代謝組學的鹿茸區段劃分與生物活性物質的分析與鑒別［D］.北京：中國農業科學院，2020.

LIU W Y.Analysis and identification of bioactive compounds from four portions of velvet antlers in sika deer （Cervus nippon） based on proteomics and metabolomics［D］.Beijing：Chinese Academy of Agricultural Sciences，2020.（in Chinese）

［43］" CHEN Y Q，POTTANAT T G，SIEGEL R W，et al.Angiopoietin-like protein 4 （ANGPTL4） is an inhibitor of endothelial lipase （EL） while the ANGPTL4/8 complex has reduced EL-inhibitory activity［J］.Heliyon，2021，7（9）：e07898.

［44］" ZHANG Z S，XU L，QIU X，et al.Fibroblast growth factor 21 （FGF21） attenuates tacrolimus-induced hepatic lipid accumulation through transcription factor EB （TFEB）-regulated lipophagy［J］.J Zhejiang Univ-Sci B，2023，24（6）：485-495.

［45］" 謝" 珂，鄭衛紅，譚" 瀟.KLKs促進腫瘤增殖、遷移與侵襲的分子機制［J］.中國生物化學與分子生物學報，2018，34（9）：935-941.

XIE K，ZHENG W H，TAN X.Molecular mechanism of Kallikrein-related peptidases in the promotion of tumor proliferation，migration and invasion［J］.Chinese Journal of Biochemistry and Molecular Biology，2018，34（9）：935-941.（in Chinese）

［46］" DOLPHIN A C，LEE A.Presynaptic calcium channels：specialized control of synaptic neurotransmitter release［J］.Nat Rev Neurosci，2020，21（4）：213-229.

［47］" 周" 雅，張禾垟，劉琳玲，等.梅花鹿成纖維細胞因子受體2基因多態性及其與茸重性狀的關聯分析［J］.中國畜牧獸醫，2022，49（8）：3006-3014.

ZHOU Y，ZHANG H Y，LIU L L，et al.Analysis of FGFR2 gene polymorphism and their correlation with antler weight traits in sika deer［J］.China Animal Husbandry amp; Veterinary Medicine，2022，49（8）：3006-3014.（in Chinese）

［48］" COURTNEY N A，BAO H，BRIGUGLIO J S，et al.Synaptotagmin 1 clamps synaptic vesicle fusion in mammalian neurons independent of complexin［J］.Nat Commun，2019，10（1）：4076.

［49］" ZHU Y F，ZHANG X，WANG L，et al.Loss of SH3GL2 promotes the migration and invasion behaviours of glioblastoma cells through activating the STAT3/MMP2 signalling［J］.J Cell Mol Med，2017，21（11）：2685-2694.

（編輯" 郭云雁）