從頭組裝豚鹿血液轉(zhuǎn)錄組揭示雌雄基因表達(dá)差異

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

（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）。豚鹿在我國僅分布于云南省西南部，曾一度被認(rèn)為已在國內(nèi)滅絕（胡婧等，2011）。近年來，野生豚鹿種群數(shù)量急劇減少，自2008 年起已被世界自然保護(hù)聯(lián)盟（IUCN）列為瀕危（EN）物種（Timminset al.，2015）。雌雄豚鹿在外形上區(qū)別較大：雄性更大且有鹿角。研究雌雄豚鹿基因表達(dá)的特征并揭示鹿角生長發(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é)和遺傳育種等研究提供了參考。但該基因組并未對豚鹿的基因進(jìn)行功能注釋，豚鹿基因轉(zhuǎn)錄表達(dá)方面的研究尚未見報道。

轉(zhuǎn)錄組可對不同組織或樣品間的基因表達(dá)量進(jìn)行定量分析、識別新的轉(zhuǎn)錄本和進(jìn)行差異表達(dá)基因（differentially expressed genes，DEGs）鑒定等（Jeukenset al.，2010），還可以進(jìn)行單核苷酸多態(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進(jìn)行了4個發(fā)育階段的多組織轉(zhuǎn)錄組分析，揭示了其生命周期中最復(fù)雜、最重要的階段是青少年過渡期；Yang 等（2016）對馬鹿Cervus canadensis鹿茸組織進(jìn)行從頭轉(zhuǎn)錄組組裝和分析，揭示了鹿茸的再生機(jī)制。為了解豚鹿基因表達(dá)的基本特征，本研究采用RNA-seq技術(shù)對成都動物園圈養(yǎng)的16 只豚鹿（6 雌10 雄）的血液進(jìn)行轉(zhuǎn)錄組測序、從頭組裝，并進(jìn)行基因功能注釋、代謝通路分析以及雌雄差異表達(dá)基因分析，以期豐富豚鹿的基因注釋信息并為其遺傳多樣性保護(hù)提供潛在的理論依據(jù)。

1 材料和方法

1.1 實驗材料

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

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

測序獲得的raw reads 進(jìn)一步進(jìn)行質(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 進(jìn)行從頭組裝，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ù)器進(jìn)行注釋（Moriyaet al.，2007）。

1.4 基因表達(dá)量分析及DEGs篩選

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

2 結(jié)果

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

16 個血液樣本提取RNA，利用Illumina HiSeq平臺進(jìn)行高通量測序，一共獲得了118.6 G raw reads，包括 51 068 480 條 raw reads，去除被接頭污染和低質(zhì)量的序列后，得到了48 472 722 條共112.6 G clean reads。使用Trinity 將過濾后的clean reads進(jìn)行從頭組裝，一共得到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ì)量良好，可用于進(jìn)一步的分析（表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 條）和細(xì)胞組分（5 855 條）。在生物過程中，細(xì)胞過程和代謝過程占主導(dǎo)；在分子功能中，轉(zhuǎn)錄本集中在連接和催化活性；在細(xì)胞組分中，轉(zhuǎn)錄本數(shù)量最豐富的類別是細(xì)胞、細(xì)胞部分和細(xì)胞器（圖1：B）。

COG 中共有39 077條轉(zhuǎn)錄本被分配到25個功能類別中，數(shù)量最多的類別是一般性功能預(yù)測，其次是信號轉(zhuǎn)導(dǎo)機(jī)制、氨基酸運輸與代謝，而分布最少的是染色質(zhì)的結(jié)構(gòu)和動力學(xué)、細(xì)胞外結(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條），其次是自然殺傷細(xì)胞介導(dǎo)的細(xì)胞毒性（ko04650，1 880條）和B細(xì)胞受體信號傳導(dǎo)途徑（ko04662，1 823條）。

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

2.3 血液表達(dá)量最豐富的基因

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

表2 豚鹿血液轉(zhuǎn)錄組中表達(dá)量前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é)干細(xì)胞多能性的信號通路（ko04550），參與的基因為賴氨酸乙酰轉(zhuǎn)移酶6A（KAT6A）。雄性顯著上調(diào)的DEGs 的GO 富集顯示，在分子功能類別中，DEGs主要富集在分子功能調(diào)節(jié)劑（GO：0098772）和酶調(diào)節(jié)活動（GO：0030234）；與細(xì)胞成分有關(guān)的DEGs 主要富集在與細(xì)胞骨架部分（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）、破骨細(xì)胞分化（ko04380）、NOD 樣受體信號通路（ko04621）和甘油磷脂代謝通路（ko00564）等通路上（圖2：B）。其中，NOD 樣受體信號通路、mTOR 信號通路、甘油磷脂代謝等通路含有部分免疫及性別相關(guān)基因，如參與免疫防御過程的TNFRSF25基因、與精子發(fā)育相關(guān)的精子鞭毛蛋白1（Spef1）基因，此外，與精子鞭毛形態(tài)形成有關(guān)的基因（Dnah17）以及睪丸表達(dá)基因（TEX264）也是雄性個體中上調(diào)的DEGs。這些雄性中高表達(dá)的DEGs可能與其生理特性相關(guān)。

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

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

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

3 討論

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

豚鹿血液轉(zhuǎn)錄組大量表達(dá)信號轉(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)的吞噬作用的大量基因。吞噬作用是重要的免疫機(jī)制之一，巨噬細(xì)胞在FcγR 受體作用下識別外源性感染物，觸發(fā)特有免疫應(yīng)答反應(yīng)（Park，2003）。在表達(dá)量前10 的基因中，表達(dá)水平最高的是MHC Ⅰ類抗原復(fù)合體基因，MHC 分子在免疫應(yīng)答過程中參與抗原識別，是動物免疫系統(tǒng)的重要組成部分（Kasahara，1999）。

通過比較雌雄性的血液轉(zhuǎn)錄組，雄性中顯著上調(diào)的DEGs 有1 781 個，遠(yuǎn)多于顯著下調(diào)的基因數(shù)量（14 個）。其中，雄性顯著高表達(dá)的基因包括一些與精子形成、睪丸功能相關(guān)基因，如Dnah17與精子鞭毛發(fā)育有關(guān)（Shaet al.，2020）。此外，在KEGG 通路富集結(jié)果中，雄性上調(diào)的DEGs 富集到mTOR 信號通路，該通路是哺乳動物體內(nèi)與細(xì)胞凋亡、細(xì)胞生長、免疫調(diào)節(jié)以及疾病發(fā)生等過程相關(guān)的重要通路（El Hianiet al.，2019），Mok等（2013）證實mTOR 信號通路對睪丸的發(fā)育調(diào)控具有重要作用，并且參與了精原細(xì)胞的分裂增殖。參與該通路的基因高表達(dá)可能與雄性個體生殖發(fā)育有關(guān)。雄性顯著高表達(dá)的基因中還涉及許多免疫相關(guān)基因，如ZC3H12A基因的缺失與自身免疫性疾病有關(guān)（Cifuenteset al.，2010）。TNFRSF25參與免疫防御和炎癥過程（Jaeckleet al.，2020），PSMB8編碼免疫蛋白酶體的催化亞基，在人類白細(xì)胞抗原Ⅰ類呈遞系統(tǒng)中起關(guān)鍵作用，與自身炎癥性疾病相關(guān)（Kitamuraet al.，2011），腫瘤壞死因子（TNF）是免疫的中樞調(diào)節(jié)劑，通過介導(dǎo)細(xì)胞存活和細(xì)胞死亡誘導(dǎo)信號，參與免疫系統(tǒng)的發(fā)育和正常運作（Web?ster & Vucic，2020）。而雌性只鑒定到1 個高表達(dá)的基因（KAT6A）與免疫相關(guān)，該基因的上調(diào)對于單核細(xì)胞向巨噬細(xì)胞分化和促炎細(xì)胞因子的產(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)皮生長因子信號通路和破骨細(xì)胞分化等通路上。破骨細(xì)胞分化的調(diào)節(jié)對正常的骨重塑至關(guān)重要（Sharmaet al.，2006）。血管內(nèi)皮生長因子（VEGF）具有促進(jìn)血管形成和骨形成的作用，刺激內(nèi)皮細(xì)胞的分裂增殖，能促進(jìn)新生血管的形成并增加血管通透性（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）的研究認(rèn)為鹿角生長和腫瘤細(xì)胞發(fā)生具有相似的發(fā)育模式，鹿角組織的快速生長特性與腫瘤發(fā)生途徑有關(guān)。因此該原癌基因在雄性中的高表達(dá)可能與鹿角的快速生長有關(guān)，以滿足鹿角快速生長的代謝需求。ADAMTS2基因是ADAMTS家族成員，Wang 等（2019）指出ADAMTS家族成員中的ADAMTS2、ADAMTS4、ADAMTS12、ADAMTS14、ADAMTS17和ADAMTS18是鹿角的特異性基因，在鹿角中的表達(dá)量很高。它們通過控制細(xì)胞外基質(zhì)的結(jié)構(gòu)和功能以及調(diào)節(jié)腫瘤微環(huán)境來抑制癌細(xì)胞的生長（Jinet al.，2007），該基因在雄鹿中的高表達(dá)可能是為了防止與鹿角生長相關(guān)的細(xì)胞因快速增殖而發(fā)生癌變。OLIG1基因與神經(jīng)嵴分化途徑有關(guān)（Betancuret al.，2010），神經(jīng)嵴細(xì)胞具有迅速增殖和分化為軟骨和神經(jīng)細(xì)胞的潛力，Wang 等（2019）發(fā)現(xiàn)鹿角基因都是從神經(jīng)嵴干細(xì)胞發(fā)育而來，且Carlson 等（2016）的研究表明，OLIG1基因突變導(dǎo)致了牛的無角性。因此猜測OLIG1基因可能也在豚鹿鹿角的進(jìn)化和發(fā)育過程中起著至關(guān)重要的作用。Snai1基因是與神經(jīng)嵴細(xì)胞遷移相關(guān)的基因，被證實在胎兒羊角中高表達(dá)（Wanget al.，2019），PTH1R和TGFBI基因在鹿角軟骨組織中高表達(dá)，并且能調(diào)節(jié)軟骨細(xì)胞和破骨細(xì)胞的增殖和分化（Faucheuxet al.，2002，2004），MMP9基因在鹿角軟骨細(xì)胞中高度表達(dá)，參與基質(zhì)降解（Faucheuxet al.，2001）。這些與鹿角發(fā)育相關(guān)DEGs 的表達(dá)量均在雄鹿中上調(diào)，表明其可能在雄性豚鹿鹿角形成和發(fā)育中發(fā)揮著重要作用，但它們在鹿角生長過程中的具體調(diào)控機(jī)理還有待進(jìn)一步挖掘探討與功能驗證。本研究結(jié)果不僅為豐富豚鹿轉(zhuǎn)錄表達(dá)提供了基礎(chǔ)數(shù)據(jù)，也為揭示偶蹄目和鹿科動物性別差異和鹿角形成機(jī)制奠定了基礎(chǔ)。