劉運，柳方圓，余建秋，牛李麗，劉選珍，李靜*

（1. 四川大學(xué)生命科學(xué)學(xué)院，生物資源與生態(tài)環(huán)境教育部重點實驗室，成都610065；2. 成都動物園，成都野生動物研究所，成都610081）

豚鹿Axis porcinus屬偶蹄目Artiodactyla 鹿科Cervidae 花鹿屬Axis，包括 3 個亞種：印度亞種A. porcinus annamiticus、南亞亞種A. porcinus kuhli和指名亞種A. porcinus porcinus，其中，印度亞種分布在巴基斯坦、尼泊爾、印度、孟加拉國和緬甸，南亞亞種分布在印度，指名亞種分布在中國、泰國、老撾、柬埔寨和越南（Tanushree & Mathur，2000）。豚鹿在我國僅分布于云南省西南部，曾一度被認為已在國內(nèi)滅絕（胡婧等，2011）。近年來，野生豚鹿種群數(shù)量急劇減少，自2008 年起已被世界自然保護聯(lián)盟（IUCN）列為瀕危（EN）物種（Timminset al.，2015）。雌雄豚鹿在外形上區(qū)別較大：雄性更大且有鹿角。研究雌雄豚鹿基因表達的特征并揭示鹿角生長發(fā)育相關(guān)基因?qū)τ诶斫馄渖飳W(xué)特性具有重要意義。

通過構(gòu)建系統(tǒng)發(fā)育樹，Abbas 等（2017）發(fā)現(xiàn)豚鹿具有遺傳多樣性較低的特點；Wang 等（2017）采用重測序技術(shù)，以單核苷酸多態(tài)性（SNPs）為標(biāo)記研究了成都動物園圈養(yǎng)豚鹿的遺傳多樣性也支持遺傳多樣性低的結(jié)論；Wang 等（2019）首次報道了豚鹿的全基因組序列，為鹿科動物比較基因組學(xué)和遺傳育種等研究提供了參考。但該基因組并未對豚鹿的基因進行功能注釋，豚鹿基因轉(zhuǎn)錄表達方面的研究尚未見報道。

轉(zhuǎn)錄組可對不同組織或樣品間的基因表達量進行定量分析、識別新的轉(zhuǎn)錄本和進行差異表達基因（differentially expressed genes，DEGs）鑒定等（Jeukenset al.，2010），還可以進行單核苷酸多態(tài)性的鑒定，開發(fā)大量簡單重復(fù)序列標(biāo)記和識別可變剪接（Garget al.，2011；Gaoet al.，2012）。Malenfant等（2013）通過構(gòu)建白尾鹿Odocoileus virginianus血液轉(zhuǎn)錄組識別單核苷酸多態(tài)性，共鑒定出66 596 個SNP，豐富了其遺傳資源信息；Jia 等（2016）對梅花鹿Cervus nippon進行了4個發(fā)育階段的多組織轉(zhuǎn)錄組分析，揭示了其生命周期中最復(fù)雜、最重要的階段是青少年過渡期；Yang 等（2016）對馬鹿Cervus canadensis鹿茸組織進行從頭轉(zhuǎn)錄組組裝和分析，揭示了鹿茸的再生機制。為了解豚鹿基因表達的基本特征，本研究采用RNA-seq技術(shù)對成都動物園圈養(yǎng)的16 只豚鹿（6 雌10 雄）的血液進行轉(zhuǎn)錄組測序、從頭組裝，并進行基因功能注釋、代謝通路分析以及雌雄差異表達基因分析，以期豐富豚鹿的基因注釋信息并為其遺傳多樣性保護提供潛在的理論依據(jù)。

1 材料和方法

1.1 實驗材料

16 份豚鹿血液樣品（10 份雄性和6 份雌性）于2017—2018年采自成都動物園圈養(yǎng)種群，?80 ℃保存，送諾禾致源公司（北京，中國）使用Illumina HiSeqTM 2000 測序儀進行150 bp 雙向測序，獲得raw reads。

1.2 數(shù)據(jù)過濾及從頭組裝

測序獲得的raw reads 進一步進行質(zhì)量控制。當(dāng)reads中未知堿基（N）的含量超過該條reads堿基數(shù)的10%以及reads 中的低質(zhì)量堿基數(shù)超過該條reads 堿基數(shù)的50%時，去除此對reads，再使用Trimmomatic（Bolgeret al.，2014）去除含有接頭污染的reads，最后得到高質(zhì)量的clean reads。使用Trinity（Haaset al.，2013）將 clean reads 進行從頭組裝，Trinity 組裝獲得的轉(zhuǎn)錄組，一般含有相似的冗余序列，再使用Cd-hit-est（Li & Godzik，2006）去除冗余序列，得到的非冗余序列用于后續(xù)分析。

1.3 功能注釋

為了獲得轉(zhuǎn)錄本的注釋信息，使用Blastx 將所有的轉(zhuǎn)錄本比對到 NR、GO、COG、Swiss-Prot 和KEGG數(shù)據(jù)庫（E-value≤10?5）。為了解基因參與的生物學(xué)通路，所有的轉(zhuǎn)錄本上傳到在線的KEGG 自動注釋服務(wù)器進行注釋（Moriyaet al.，2007）。

1.4 基因表達量分析及DEGs篩選

使用RSEM 1.2.2（Li& Dewey，2011）對每個樣本基因表達量水平進行估計。表達量的結(jié)果使用FPKM 值表示。使用 DEseq2（Loveet al.，2014）進行不同樣本之間的DEGs分析，篩選標(biāo)準(zhǔn)為log2|FC|>1且P<0.05。獲得DEGs 集后，使用Bioconductor R語言包ClusterProfiler 3.8.1（Yuet al.，2012）對其進行GO 功能富集和KEGG 通路注釋，以進一步評估DEGs的功能。

2 結(jié)果

2.1 轉(zhuǎn)錄組測序和組裝

16 個血液樣本提取RNA，利用Illumina HiSeq平臺進行高通量測序，一共獲得了118.6 G raw reads，包括 51 068 480 條 raw reads，去除被接頭污染和低質(zhì)量的序列后，得到了48 472 722 條共112.6 G clean reads。使用Trinity 將過濾后的clean reads進行從頭組裝，一共得到615 548條unigenes，長度為200~18 935 bp，平均長度為600.48 bp，Con?tig N50 為803 bp，GC 含量平均值為46.89%。大多數(shù)轉(zhuǎn)錄本的長度為200~400 bp（260 743條），1 000~2 000 bp（50 248 條）的次之。由此可見，拼接組裝后的轉(zhuǎn)錄組質(zhì)量良好，可用于進一步的分析（表1）。

表1 豚鹿轉(zhuǎn)錄本組裝結(jié)果統(tǒng)計Table 1 Summary of transcriptome assembly of Axis porcinus

2.2 轉(zhuǎn)錄本注釋

與NR 的比對結(jié)果顯示，有118 326 條轉(zhuǎn)錄本（37.93%）獲得注釋信息；在Swiss-Prot 中注釋到的轉(zhuǎn)錄本為96 171 條（30.83%）；COG 中注釋到的轉(zhuǎn)錄本為39 077條（12.53%）；在KEGG中獲得注釋信息的轉(zhuǎn)錄本共42 428 條（6.94%）。在四大數(shù)據(jù)庫中均獲得注釋信息的轉(zhuǎn)錄本有6 307條（圖1：A）。

GO分類顯示15 971條轉(zhuǎn)錄本被分配到了生物過程（8 604 條）、分子功能（10 594 條）和細胞組分（5 855 條）。在生物過程中，細胞過程和代謝過程占主導(dǎo)；在分子功能中，轉(zhuǎn)錄本集中在連接和催化活性；在細胞組分中，轉(zhuǎn)錄本數(shù)量最豐富的類別是細胞、細胞部分和細胞器（圖1：B）。

COG 中共有39 077條轉(zhuǎn)錄本被分配到25個功能類別中，數(shù)量最多的類別是一般性功能預(yù)測，其次是信號轉(zhuǎn)導(dǎo)機制、氨基酸運輸與代謝，而分布最少的是染色質(zhì)的結(jié)構(gòu)和動力學(xué)、細胞外結(jié)構(gòu)和核結(jié)構(gòu)（圖1：C）。

KEGG 中注釋的42 428 條轉(zhuǎn)錄本被劃分到502 個代謝通路中，參與信號轉(zhuǎn)導(dǎo)通路（ko09132）的最多，共有7 144 條；其次是病毒感染（ko09172，4 857 條）和與免疫系統(tǒng)相關(guān)的通路（ko09151，4 064 條）（圖1：D）。根據(jù)免疫系統(tǒng)分類，免疫相關(guān)基因被映射到21 個免疫通路：涉及FcγR介導(dǎo)的吞噬作用（ko04666）的最多（1 894條），其次是自然殺傷細胞介導(dǎo)的細胞毒性（ko04650，1 880條）和B細胞受體信號傳導(dǎo)途徑（ko04662，1 823條）。

圖1 轉(zhuǎn)錄本與數(shù)據(jù)庫同源性比對注釋結(jié)果Fig. 1 Results of homology search of transcripts against database

2.3 血液表達量最豐富的基因

根據(jù)NR數(shù)據(jù)庫的比對并使用FPKM 值統(tǒng)計轉(zhuǎn)錄本表達水平，表達量最高的前10 個轉(zhuǎn)錄本主要包括主要組織相容性復(fù)合體（MHC）Ⅰ類抗原、幾種核糖體蛋白和細胞色素氧化酶等基因，其中，HBA基因也是血液高表達的基因，編碼了血紅蛋白的α亞基（表2）。

表2 豚鹿血液轉(zhuǎn)錄組中表達量前10個轉(zhuǎn)錄本Table 2 The top 10 most abundant transcripts in the blood transcriptome of Axis porcinus

2.4 雌雄DEGs與富集分析

共鑒定到1 793 個性別DEGs。雄性顯著上調(diào)的有1 781 個，顯著下調(diào)的有14 個，顯著下調(diào)的僅富集到1 條KEGG 信號通路：調(diào)節(jié)干細胞多能性的信號通路（ko04550），參與的基因為賴氨酸乙酰轉(zhuǎn)移酶6A（KAT6A）。雄性顯著上調(diào)的DEGs 的GO 富集顯示，在分子功能類別中，DEGs主要富集在分子功能調(diào)節(jié)劑（GO：0098772）和酶調(diào)節(jié)活動（GO：0030234）；與細胞成分有關(guān)的DEGs 主要富集在與細胞骨架部分（GO：0044430）和微管組織中心（GO：0005815）有關(guān)的條目；在生物過程中，DEGs主要富集在免疫系統(tǒng)過程（GO：0002376）和囊泡介導(dǎo)轉(zhuǎn)運（GO：0016192）（圖2：A）。KEGG 富集顯示，顯著上調(diào)的DEGs 主要富集在血管內(nèi)皮生長因子（VEGF）信號通路（ko04370）、趨化因子信號通路（ko04062）、溶酶體（ko04142）、mTOR 信號通路（ko04150）、破骨細胞分化（ko04380）、NOD 樣受體信號通路（ko04621）和甘油磷脂代謝通路（ko00564）等通路上（圖2：B）。其中，NOD 樣受體信號通路、mTOR 信號通路、甘油磷脂代謝等通路含有部分免疫及性別相關(guān)基因，如參與免疫防御過程的TNFRSF25基因、與精子發(fā)育相關(guān)的精子鞭毛蛋白1（Spef1）基因，此外，與精子鞭毛形態(tài)形成有關(guān)的基因（Dnah17）以及睪丸表達基因（TEX264）也是雄性個體中上調(diào)的DEGs。這些雄性中高表達的DEGs可能與其生理特性相關(guān)。

圖2 雄性豚鹿上調(diào)差異表達基因富集Fig. 2 Enrichment of up-regulated differentially expressed genes in male Axis porcinus

在雄性中顯著高表達的DEGs 中，還發(fā)現(xiàn)了一些過去被報道參與鹿角生長發(fā)育相關(guān)的基因，包括與癌細胞生長有關(guān)的基因Rela、ADAMTS2，與神經(jīng)嵴生長分化有關(guān)的基因OLIG1、Snai1，與軟骨組織分化有關(guān)的基因PTH1R、TGFBI、MMP9，這些基因的表達量均顯著高于雌性（圖3）。

圖3 部分差異表達基因在豚鹿中的FPKM 表達量Fig. 3 FPKM values of some differentially expressed genes of Axis porcinus

3 討論

目前關(guān)于豚鹿的基因組雖已有報道（Wanget al.，2019），但該基因組缺乏基因注釋信息。轉(zhuǎn)錄組的從頭組裝能夠為基因表達提供重要信息，對基因組的補充注釋具有重要意義。本研究基于16只豚鹿個體進行血液RNA-seq 測序，從頭組裝了一個血液轉(zhuǎn)錄組，N50 為803 bp，組裝質(zhì)量較好?；诓煌臄?shù)據(jù)庫對轉(zhuǎn)錄本進行注釋，結(jié)果顯示，僅有37.93%的序列獲得了注釋信息，低于水牛Bubalus bubalis（77.91%）（Denget al.，2016）、山 羊Capra hircus（45.41%）（Liuet al.，2013）等的注釋率。這可能是由于參與組裝的豚鹿個體多，測序的轉(zhuǎn)錄本中含有較多短序列；其次，注釋信息不完善或可能存在豚鹿特有的新基因等。

豚鹿血液轉(zhuǎn)錄組大量表達信號轉(zhuǎn)導(dǎo)、病毒感染和免疫系統(tǒng)相關(guān)基因，這與血液的重要生理功能密切相關(guān)。KEGG 注釋顯示許多轉(zhuǎn)錄本注釋到與免疫相關(guān)的通路，這與大熊貓Ailuropoda melano?leuca（Duet al.，2015）、水牛（Denget al.，2016）等的結(jié)果一致，包括參與FcγR 介導(dǎo)的吞噬作用的大量基因。吞噬作用是重要的免疫機制之一，巨噬細胞在FcγR 受體作用下識別外源性感染物，觸發(fā)特有免疫應(yīng)答反應(yīng)（Park，2003）。在表達量前10 的基因中，表達水平最高的是MHC Ⅰ類抗原復(fù)合體基因，MHC 分子在免疫應(yīng)答過程中參與抗原識別，是動物免疫系統(tǒng)的重要組成部分（Kasahara，1999）。

通過比較雌雄性的血液轉(zhuǎn)錄組，雄性中顯著上調(diào)的DEGs 有1 781 個，遠多于顯著下調(diào)的基因數(shù)量（14 個）。其中，雄性顯著高表達的基因包括一些與精子形成、睪丸功能相關(guān)基因，如Dnah17與精子鞭毛發(fā)育有關(guān)（Shaet al.，2020）。此外，在KEGG 通路富集結(jié)果中，雄性上調(diào)的DEGs 富集到mTOR 信號通路，該通路是哺乳動物體內(nèi)與細胞凋亡、細胞生長、免疫調(diào)節(jié)以及疾病發(fā)生等過程相關(guān)的重要通路（El Hianiet al.，2019），Mok等（2013）證實mTOR 信號通路對睪丸的發(fā)育調(diào)控具有重要作用，并且參與了精原細胞的分裂增殖。參與該通路的基因高表達可能與雄性個體生殖發(fā)育有關(guān)。雄性顯著高表達的基因中還涉及許多免疫相關(guān)基因，如ZC3H12A基因的缺失與自身免疫性疾病有關(guān)（Cifuenteset al.，2010）。TNFRSF25參與免疫防御和炎癥過程（Jaeckleet al.，2020），PSMB8編碼免疫蛋白酶體的催化亞基，在人類白細胞抗原Ⅰ類呈遞系統(tǒng)中起關(guān)鍵作用，與自身炎癥性疾病相關(guān)（Kitamuraet al.，2011），腫瘤壞死因子（TNF）是免疫的中樞調(diào)節(jié)劑，通過介導(dǎo)細胞存活和細胞死亡誘導(dǎo)信號，參與免疫系統(tǒng)的發(fā)育和正常運作（Web?ster & Vucic，2020）。而雌性只鑒定到1 個高表達的基因（KAT6A）與免疫相關(guān)，該基因的上調(diào)對于單核細胞向巨噬細胞分化和促炎細胞因子的產(chǎn)生很重要（Wiesel-Motiuk & Assaraf，2020）。這些結(jié)果表明，豚鹿在自身免疫方面存在顯著性別差異，這可由哺乳動物的免疫系統(tǒng)具有廣泛的性別二態(tài)性來解釋（Gal-Ozet al.，2019），不同的因素包括遺傳因素和激素介質(zhì)都可以解釋免疫反應(yīng)中基于性別的差異（Oertelt-Prigione，2012）。

雌雄豚鹿的DEGs 還揭示了與鹿角生長發(fā)育相關(guān)的基因。KEGG 富集發(fā)現(xiàn)，雄鹿中上調(diào)的DEGs顯著富集到血管內(nèi)皮生長因子信號通路和破骨細胞分化等通路上。破骨細胞分化的調(diào)節(jié)對正常的骨重塑至關(guān)重要（Sharmaet al.，2006）。血管內(nèi)皮生長因子（VEGF）具有促進血管形成和骨形成的作用，刺激內(nèi)皮細胞的分裂增殖，能促進新生血管的形成并增加血管通透性（Apteet al.，2019）。鹿角的生長需要血管大量供血，因此富集到VEGF信號通路的DEGs 可能與鹿角的生長形成有關(guān)。通過基因的數(shù)據(jù)庫比對信息，在雄性上調(diào)的DEGs中，鑒定了可能與鹿角生長相關(guān)的基因：Rela、ADAMTS2、OLIG1、Snai1、PTH1R、TGFBI、MMP9。Rela是一種原癌基因（Lu & Yarbrough，2015），Wang等（2019）的研究認為鹿角生長和腫瘤細胞發(fā)生具有相似的發(fā)育模式，鹿角組織的快速生長特性與腫瘤發(fā)生途徑有關(guān)。因此該原癌基因在雄性中的高表達可能與鹿角的快速生長有關(guān)，以滿足鹿角快速生長的代謝需求。ADAMTS2基因是ADAMTS家族成員，Wang 等（2019）指出ADAMTS家族成員中的ADAMTS2、ADAMTS4、ADAMTS12、ADAMTS14、ADAMTS17和ADAMTS18是鹿角的特異性基因，在鹿角中的表達量很高。它們通過控制細胞外基質(zhì)的結(jié)構(gòu)和功能以及調(diào)節(jié)腫瘤微環(huán)境來抑制癌細胞的生長（Jinet al.，2007），該基因在雄鹿中的高表達可能是為了防止與鹿角生長相關(guān)的細胞因快速增殖而發(fā)生癌變。OLIG1基因與神經(jīng)嵴分化途徑有關(guān)（Betancuret al.，2010），神經(jīng)嵴細胞具有迅速增殖和分化為軟骨和神經(jīng)細胞的潛力，Wang 等（2019）發(fā)現(xiàn)鹿角基因都是從神經(jīng)嵴干細胞發(fā)育而來，且Carlson 等（2016）的研究表明，OLIG1基因突變導(dǎo)致了牛的無角性。因此猜測OLIG1基因可能也在豚鹿鹿角的進化和發(fā)育過程中起著至關(guān)重要的作用。Snai1基因是與神經(jīng)嵴細胞遷移相關(guān)的基因，被證實在胎兒羊角中高表達（Wanget al.，2019），PTH1R和TGFBI基因在鹿角軟骨組織中高表達，并且能調(diào)節(jié)軟骨細胞和破骨細胞的增殖和分化（Faucheuxet al.，2002，2004），MMP9基因在鹿角軟骨細胞中高度表達，參與基質(zhì)降解（Faucheuxet al.，2001）。這些與鹿角發(fā)育相關(guān)DEGs 的表達量均在雄鹿中上調(diào)，表明其可能在雄性豚鹿鹿角形成和發(fā)育中發(fā)揮著重要作用，但它們在鹿角生長過程中的具體調(diào)控機理還有待進一步挖掘探討與功能驗證。本研究結(jié)果不僅為豐富豚鹿轉(zhuǎn)錄表達提供了基礎(chǔ)數(shù)據(jù)，也為揭示偶蹄目和鹿科動物性別差異和鹿角形成機制奠定了基礎(chǔ)。