蔣蘭， 張雪艷， 王俊茵， 李靜

(1. 四川大學(xué)生命科學(xué)學(xué)院，生物資源與生態(tài)環(huán)境教育部重點(diǎn)實(shí)驗(yàn)室，成都610065；2. 峨眉山景區(qū)管委會(huì)，峨眉山生物多樣性保護(hù)研究所，四川峨眉山614200)

豹貓Prionailurusbengalensis屬食肉目Carnivora貓科Felidae豹貓屬，地理分布極為廣泛，從亞洲南部和西南部的熱帶雨林到中國(guó)北方和俄羅斯遠(yuǎn)東地區(qū)的森林(Hughes，1998；Woolfetal.，2002)，在中國(guó)，根據(jù)地理分布豹貓被分為5個(gè)亞種(Rossetal.，2015)。盡管豹貓資源豐富，但人類活動(dòng)使其棲息地不斷縮小，同時(shí)對(duì)其皮毛的需求也導(dǎo)致豹貓的生存受到嚴(yán)重威脅(Hughes，1998)。目前中國(guó)的野生豹貓種群數(shù)量正在下降，分布范圍也逐漸縮小，一些原有的分布區(qū)域，如寧夏、青海、山東等地區(qū)的豹貓數(shù)量已非常稀少。《中國(guó)瀕危動(dòng)物紅皮書》(汪松，1998)將豹貓列為易危種。

目前國(guó)內(nèi)外有關(guān)豹貓的文獻(xiàn)報(bào)道集中在遺傳多樣性和系統(tǒng)發(fā)育等研究。研究人員分別利用12s RNA、cytb序列和cDNA序列探討西表島貓Prionailurusiriomotensis和豹貓系統(tǒng)發(fā)育關(guān)系，結(jié)果表明二者分歧時(shí)間較近，西表島貓應(yīng)為豹貓的一個(gè)亞種(Masudaetal.，1994；Suzukietal.，1994)。Masuda和Yoshida(1995)基于cytb序列認(rèn)為馬島貓Prionailurusbengalensiseuptilurus也是豹貓的一個(gè)亞種。Saka等(2018)發(fā)現(xiàn)西表島和馬島的豹貓種群的MHC多樣性遠(yuǎn)低于家貓F(tuán)eliscatus種群，暗示島嶼上的豹貓由于近親繁殖或地理隔離，其對(duì)病原體的抵抗力降低。Ko等(2018)基于微衛(wèi)星分析了韓國(guó)豹貓的遺傳多樣性，發(fā)現(xiàn)其平均等位基因數(shù)和雜合度分別為3.8和0.41，多樣性低于其他12種貓科動(dòng)物57個(gè)種群的平均值。白素英等(2004)利用隨機(jī)擴(kuò)增多態(tài)性DNA(RAPD)研究中國(guó)6個(gè)豹貓種群結(jié)構(gòu)發(fā)現(xiàn)，各種群的關(guān)系與種群間的地理距離及亞種劃分一致。目前國(guó)內(nèi)關(guān)于豹貓的研究集中在分類、分布、食性以及解剖等方面(張淑云等，1988；白素英等，2004；Xiongetal.，2016)，而關(guān)于其基因表達(dá)的研究還未見報(bào)道。

轉(zhuǎn)錄組測(cè)序技術(shù)(RNA-seq)是一種廣泛應(yīng)用于研究基因表達(dá)調(diào)控的分析方法，通過有效識(shí)別和量化RNA轉(zhuǎn)錄本，從而獲得其特定生理狀況的轉(zhuǎn)錄組信息(Fengetal.，2012；Spradlingetal.，2013)。RNA-seq使得從整體水平分析動(dòng)物基因結(jié)構(gòu)及基因功能成為可能，更有利于發(fā)現(xiàn)生物學(xué)過程、免疫反應(yīng)基因、疾病產(chǎn)生過程中的基因表達(dá)譜變化以及免疫受體多樣化機(jī)制(Wang & Kirkness，2005)。由于RNA-seq技術(shù)對(duì)樣本新鮮程度要求較高，一定程度上阻礙了其在野生動(dòng)物相關(guān)研究中的應(yīng)用。本研究基于1只意外死亡的野生豹貓個(gè)體，采用RNA-seq對(duì)其大腦、心臟、腎臟、肝臟、肺和骨骼肌 6個(gè)組織進(jìn)行轉(zhuǎn)錄組測(cè)序，通過從頭組裝、注釋，報(bào)道了一個(gè)高質(zhì)量的豹貓轉(zhuǎn)錄組。這些數(shù)據(jù)為理解豹貓的遺傳背景及基因表達(dá)譜提供了寶貴的資源，也為這一小型貓科動(dòng)物的基因組注釋和深入開展保護(hù)遺傳學(xué)研究奠定了基礎(chǔ)。

1 材料與方法

1.1 樣品采集與測(cè)序

豹貓各器官的組織樣品均來自于1只在峨眉山被汽車意外撞死的雌性成年個(gè)體，被管理人員發(fā)現(xiàn)后提供給本實(shí)驗(yàn)室。共采集了6個(gè)不同器官的組織樣品，分別為大腦、心臟、腎臟、肝臟、肺和骨骼肌，新鮮的組織樣品迅速置于液氮中保存，其余樣品-80 ℃儲(chǔ)存。樣品的總RNA提取和轉(zhuǎn)錄組建庫測(cè)序工作由北京諾禾致源生物信息科技有限公司完成，采用Illumina HiSeq 4000測(cè)序平臺(tái)對(duì)6個(gè)樣品的cDNA文庫分別進(jìn)行150 bp的雙端測(cè)序。

1.2 質(zhì)控與過濾

轉(zhuǎn)錄組組裝前需對(duì)原始序列進(jìn)行質(zhì)量控制，過濾標(biāo)準(zhǔn)為去除含N比例大于10%、低質(zhì)量(質(zhì)量值<5)堿基數(shù)比例大于50%、含有接頭污染的序列，最終得到高質(zhì)量的clean reads用于后續(xù)分析。

1.3 轉(zhuǎn)錄組從頭組裝

使用Trinity(Haasetal.，2013)對(duì)過濾后的clean reads進(jìn)行denovo組裝，僅輸出片段長(zhǎng)度大于300 bp的contig序列。由于Trinity組裝獲得的轉(zhuǎn)錄組中通常含有大量相似的冗余序列，故最后使用Cd-hit-est (Li & Godzik，2006)去除冗余序列，得到非冗余轉(zhuǎn)錄本用于后續(xù)分析。

1.4 功能注釋

為了獲得較為全面的基因功能信息，利用blastx(Camachoetal.，2009)將Unigene與NCBI非冗余蛋白質(zhì)(NR)數(shù)據(jù)庫、euKaryotic Ortholog Groups(KOG)數(shù)據(jù)庫、Swiss-Prot數(shù)據(jù)庫進(jìn)行比對(duì)，E值設(shè)為1E-5。再根據(jù)NR數(shù)據(jù)庫的比對(duì)結(jié)果進(jìn)行Gene Ontology(GO)和Kyoto Encyclopedia of Genes and Genomes(KEGG)注釋的分類統(tǒng)計(jì)，分別通過在線工具Web Gene Ontology Annotation Plot(WEGO)(Yeetal.，2006)和KEGG Automatic Annotation Server(KAAS)完成(Moriyaetal.，2007)。在使用WEGO進(jìn)行分類統(tǒng)計(jì)之前，需先將轉(zhuǎn)錄本與NR數(shù)據(jù)庫比對(duì)的結(jié)果導(dǎo)入Blast2GO(Conesaetal.，2005)搜索相關(guān)的GO功能條目。

1.5 轉(zhuǎn)錄本的表達(dá)量和組織特異性分析

2 結(jié)果

2.1 測(cè)序和質(zhì)控

共獲得171 370 884條150 bp的raw reads，共計(jì)51.4 Gb。質(zhì)控過濾后得到154 489 126條(90.15%)clean reads，共計(jì)45.42 Gb，每個(gè)組織的數(shù)據(jù)量為5.54～13.46 Gb(表1)，表明測(cè)序數(shù)據(jù)質(zhì)量符合要求。

表1 測(cè)序結(jié)果和質(zhì)量控制Table 1 Summary of sequencing and quality control

2.2 轉(zhuǎn)錄組從頭組裝

使用Trinity將clean reads進(jìn)行denovo組裝，最終得到369 246條轉(zhuǎn)錄本，轉(zhuǎn)錄本長(zhǎng)度301～31 696 bp，平均長(zhǎng)度1 465 bp，Contig N50為2 660 bp(表2)。拼接的序列中，155 985條(42.24%)序列長(zhǎng)度在1 000 bp以上，組裝質(zhì)量較高。大多數(shù)轉(zhuǎn)錄本長(zhǎng)度為300～400 bp和1 000～2 000 bp(圖1)。使用Cd-hit-est對(duì)Trinity組裝的轉(zhuǎn)錄本去冗余，最終得到269 412條非冗余轉(zhuǎn)錄本用于注釋分析。

表2 豹貓轉(zhuǎn)錄組組裝結(jié)果Table 2 Summary of transcriptome assembly of Prionailurus bengalensis

2.3 轉(zhuǎn)錄本的功能注釋

使用blastx將去冗余后的轉(zhuǎn)錄本與NR、Swiss-Prot、KOG和KEGG 4個(gè)公共數(shù)據(jù)進(jìn)行同源性比對(duì)注釋。NR蛋白數(shù)據(jù)庫注釋到轉(zhuǎn)錄本110 997條(41.13%)；Swiss-Prot蛋白數(shù)據(jù)庫注釋到轉(zhuǎn)錄本 98 192條(36.39%)；KOG數(shù)據(jù)庫注釋到轉(zhuǎn)錄本97 616條(36.17%)；KEGG數(shù)據(jù)庫獲得注釋信息轉(zhuǎn)錄本65 895條(24.42%)。4個(gè)數(shù)據(jù)庫注釋結(jié)果的韋恩圖顯示，共注釋到轉(zhuǎn)錄本114 517條(42.44%)，63 860 條轉(zhuǎn)錄本在所有數(shù)據(jù)庫中都有注釋(圖2：A)。

NR數(shù)據(jù)庫比對(duì)表明，73%的序列有較強(qiáng)的同源性(E值≤1.0E-30)，其中，25%的序列比對(duì)E值為0，39%的為0～1.0E-45，9%的為1.0E-45～1.0E-30(圖2：B)。相似度分布顯示，83%的序列比對(duì)相似度≥80%，僅有17%的為30%～80%(圖2：C)。物種間比對(duì)結(jié)果顯示，豹貓的轉(zhuǎn)錄本與貓科動(dòng)物的相似性高，約49%的序列與貓科動(dòng)物基因組(家貓：23%，金錢豹Pantherapardus：13%，獵豹Acinoyxjubatus：8%，虎Pantheratigris：5%)己經(jīng)注釋的基因同源(圖2：D)。

轉(zhuǎn)錄本的GO富集顯示，33 262條轉(zhuǎn)錄本注釋到247 376個(gè)GO功能條目(圖3：A)，包括生物過程(113 399條，45.84%)、分子功能(37 773條，15.27%)和細(xì)胞組分(96 204條，38.89%)。生物過程中，細(xì)胞過程(20 431條)、代謝過程(16 659條)和生物調(diào)節(jié)(12 227條)最多；細(xì)胞組分中，細(xì)胞(18 428條)、細(xì)胞部分(18 261條)和細(xì)胞器(13 859條)最豐富；分子功能中，最豐富的為連接(18 061條)和催化活性(11 917條)。

根據(jù)同源性比對(duì)KOG數(shù)據(jù)庫，97 616條(36.17%)轉(zhuǎn)錄本分別在25種類別中被注釋到(E值≤1E-5)(圖3：B)，信號(hào)轉(zhuǎn)導(dǎo)機(jī)制和一般功能預(yù)測(cè)的轉(zhuǎn)錄本最多，其次是翻譯后修飾、蛋白轉(zhuǎn)換、分子伴侶和轉(zhuǎn)錄。

KEGG注釋結(jié)果顯示有65 895條(24.42%)轉(zhuǎn)錄本被分配到386個(gè)KEGG通路上，其中與人類疾病相關(guān)的轉(zhuǎn)錄本(26 761條)最多，其次是生物系統(tǒng)(22 997條)(圖3：C)。在二級(jí)分類中，參與信號(hào)轉(zhuǎn)導(dǎo)的轉(zhuǎn)錄本(10 039條)和傳染?。翰《?10 003條)的轉(zhuǎn)錄本數(shù)量最多。

2.4 轉(zhuǎn)錄本的組織特異性表達(dá)

轉(zhuǎn)錄本TSI值的頻數(shù)分布圖顯示(圖4：A)，僅17條轉(zhuǎn)錄本TSI值<0.15，它們?cè)谒?個(gè)組織中表達(dá)差異很小，為低組織特異性基因，其中可能包括了一些重要的管家基因；而39.65%(66 763條)的轉(zhuǎn)錄本TSI值在0.15～0.85之間，為中等組織特異性基因，它們的表達(dá)量在各組織中存在差異，TSI值越大差異越大；60.34%的轉(zhuǎn)錄本TSI值>0.85，為高組織特異性基因，其在部分組織中的表達(dá)與其他組織存在明顯差異，其中24.5%(41 354條)的轉(zhuǎn)錄本TSI值達(dá)到1.0，它們嚴(yán)格地在特定的某一種組織中表達(dá)。由表達(dá)量熱圖(圖4：B)可見，這3類轉(zhuǎn)錄本在各組織中的表達(dá)模式不同。

2.5 各組織高表達(dá)轉(zhuǎn)錄本

為了解這些組織特異表達(dá)基因是否為各組織高表達(dá)基因，統(tǒng)計(jì)了6個(gè)組織中表達(dá)量最高的10條轉(zhuǎn)錄本(表3)，共包括39條轉(zhuǎn)錄本，其TSI值為0.48～1.0(平均0.83)。其中骨骼肌中10條轉(zhuǎn)錄本FPKM值均超過10 000，而腦組織中9條轉(zhuǎn)錄本FPKM值為1 213.09～3 074.01，僅1條>10 000，顯示不同組織中基因表達(dá)量差異較大。這些轉(zhuǎn)錄本中，6條轉(zhuǎn)錄本(T_90962_c0_g2_i7、T_92170_c4_g8_i1、T_92170_c4_g7_i1、T_90443_c2_g7_i2、T_92619_c3_g1_i2和T_87893_c2_g3_i5)在骨骼肌和心臟中高表達(dá)，4條轉(zhuǎn)錄本(T_90962_c0_g2_i7、T_84299_c4_g5_i2、T_92588_c6_g2_i6和T_92288_c9_g4_i1)在腎臟和肝臟中高表達(dá)。39條轉(zhuǎn)錄本中26條(61.5%)為高組織特異性表達(dá)，包括8條在骨骼肌中高表達(dá)，7條在肺中高表達(dá)，6條在肝臟中高表達(dá)，腎臟、心臟、大腦中高表達(dá)的分別有5條、4條、3條；其他13個(gè)高表達(dá)轉(zhuǎn)錄本的組織特異性較低，如COX1在除肺以外的5個(gè)組織中都高表達(dá)。

表3 豹貓6個(gè)組織中表達(dá)量最高的10條轉(zhuǎn)錄本的NR注釋結(jié)果及其組織特異性Table 3 NR annotation of the top 10 transcripts with the highest FPKM value in the 6 tissues of Prionailurus bengalensis and TSI analysis

續(xù)表3

3 討論

作為一種分布廣泛的小型貓科動(dòng)物，過去關(guān)于豹貓的研究大多關(guān)注其種群分布、遺傳多樣性水平和系統(tǒng)發(fā)育關(guān)系等(Mukherjeeetal.，2010；Pateletal.，2017)。Ito等(2020)基于基因組SNP分析了豹貓遺傳多樣性，Bredemeyer等(2021)報(bào)道了豹貓的全基因組，但迄今為止關(guān)于豹貓基因表達(dá)調(diào)控相關(guān)研究仍非常缺乏。本研究采用RNA-seq對(duì)豹貓6個(gè)組織的轉(zhuǎn)錄組進(jìn)行了測(cè)序、從頭組裝、功能注釋及不同組織的基因表達(dá)譜分析，從而為豐富豹貓的遺傳資源、輔助基因組注釋、開發(fā)分子標(biāo)記等提供了重要數(shù)據(jù)，也對(duì)于理解豹貓的基因表達(dá)調(diào)控，開展保護(hù)遺傳學(xué)研究奠定了基礎(chǔ)。

本研究基于豹貓6個(gè)組織的大量轉(zhuǎn)錄本從頭組裝了豹貓轉(zhuǎn)錄組，能更全面地體現(xiàn)豹貓基因表達(dá)情況，組裝質(zhì)量良好：N50長(zhǎng)度為2 660 bp，42.24%的轉(zhuǎn)錄本長(zhǎng)度超過1 000 bp。研究顯示，越長(zhǎng)的轉(zhuǎn)錄本序列越完整，組裝質(zhì)量越高，更容易獲得較多關(guān)于基因的信息(Mengetal.，2015)。與發(fā)表的東北虎Pantheratigrisaltaica肺轉(zhuǎn)錄組(319 bp)和6個(gè)組織合并轉(zhuǎn)錄本(332 bp)(Luetal.，2016)、小尾寒羊和杜泊羊骨骼肌轉(zhuǎn)錄組(735 bp和706 bp)(張春蘭，2016)、梭子蟹Portunustrituberculatus雌性(954 bp)和雄性(823 bp)(Wangetal.，2018)轉(zhuǎn)錄組序列平均長(zhǎng)度相比，豹貓轉(zhuǎn)錄組質(zhì)量更高。GO、KOG和KEGG注釋結(jié)果顯示出豹貓各組織表達(dá)不同生理功能相關(guān)的多種基因。僅有42.44%(114 517條)的豹貓轉(zhuǎn)錄本具有注釋信息，超過50%的轉(zhuǎn)錄本無法比對(duì)到現(xiàn)有數(shù)據(jù)庫，這可能是由于豹貓的基因組和蛋白質(zhì)相關(guān)信息仍不完善，大量轉(zhuǎn)錄本可能是豹貓?zhí)赜械霓D(zhuǎn)錄本，深入研究這些新轉(zhuǎn)錄本可揭示豹貓獨(dú)特的生理、遺傳特征。

豹貓6個(gè)組織的TSI分布顯示，TSI<0.15的管家基因最少，僅17條(0.01%)，60.34%的轉(zhuǎn)錄本TSI>0.85，其表達(dá)具有高度組織特異性。這顯示盡管共有同一個(gè)基因組，但豹貓?jiān)诓煌M織中的基因轉(zhuǎn)錄和表達(dá)卻存在很大差異，以保障不同組織執(zhí)行不同的生理功能。Su等(2002)統(tǒng)計(jì)老鼠45個(gè)正常組織、Yanni等(2005)統(tǒng)計(jì)人類Homosapiens12種組織的TSI也顯示，TSI值為0.9～1.0的高組織特異性表達(dá)基因最多，這與本研究結(jié)果一致。但不同的是，人類中57%、老鼠中56%的基因都屬于中等組織特異性表達(dá)，即0.15≤TSI≤0.85。這可能是由于轉(zhuǎn)錄本注釋質(zhì)量的差異。人和老鼠具有高質(zhì)量基因組，轉(zhuǎn)錄本注釋完善，而豹貓缺乏相關(guān)的基因組和蛋白質(zhì)信息，本研究從頭組裝的轉(zhuǎn)錄本高達(dá)26萬余條，許多轉(zhuǎn)錄本并非全長(zhǎng)基因，存在大量難以注釋的短序列。

為進(jìn)一步了解這些組織特異性表達(dá)的基因是否也是高表達(dá)的基因，統(tǒng)計(jì)了豹貓每個(gè)組織中表達(dá)量最高的前10條轉(zhuǎn)錄本。共計(jì)39個(gè)高表達(dá)的基因中，26個(gè)都屬于在特定組織高表達(dá)的基因，這些基因都是與該組織特定功能密切相關(guān)的重要基因。如肝臟中特異高表達(dá)的RBP4與肝臟的脂肪降解功能相關(guān)，人的多組織基因表達(dá)圖譜也顯示該基因在肝臟和脂肪中特異表達(dá)(Fagerbergetal.，2014)。肺中特異高表達(dá)的SFTPC編碼肺表面活性物質(zhì)蛋白C，它通過降低覆蓋肺部的液體的表面張力來維持肺組織的穩(wěn)定性，SFTPC突變能引發(fā)嬰幼兒間質(zhì)性肺病，其過表達(dá)可抑制肺癌細(xì)胞的增殖(Hayasakaetal.，2018；Lietal.，2019)。大腦中特異高表達(dá)的PCP4與神經(jīng)元的功能相關(guān)，小鼠胚胎發(fā)生過程中PCP4的過表達(dá)誘導(dǎo)神經(jīng)元提前分化成熟，而在成年期過表達(dá)會(huì)導(dǎo)致小鼠學(xué)習(xí)障礙(Mouton-Ligeretal.，2014；Reneltetal.，2014)。此外我們還鑒定了不同組織共同高表達(dá)的基因。如COX1編碼細(xì)胞色素c氧化酶Ⅰ，是線粒體呼吸鏈的末端組分，該酶可將電子從還原的細(xì)胞色素c轉(zhuǎn)移到分子氧，并有助于維持線粒體內(nèi)膜的電化學(xué)梯度(Baklouti-Gargourietal.，2013)，為ATP合成酶合成ATP提供電化學(xué)梯度(Dennerlein & Rehling，2015)，是細(xì)胞呼吸過程中重要的功能基因。COX1在心臟、肝臟、腎臟、大腦和骨骼肌中都是表達(dá)量前10的基因，在肺中的表達(dá)量也相對(duì)較高，證明該基因在各個(gè)組織中均發(fā)揮重要生理功能。TMSB4是一種廣泛分布于真核細(xì)胞中的多功能肽，與細(xì)胞的移動(dòng)性、分化能力及細(xì)胞分裂密切相關(guān)(Kuzan，2016)。TMSB4是肺、腎臟和大腦組織中共同高表達(dá)的基因，Du等(2015)發(fā)現(xiàn)其在大熊貓血液中也是高表達(dá)的基因之一，它可能也是維持細(xì)胞基本生理功能的重要基因。心臟和骨骼肌都屬于肌肉組織，具有收縮性，本研究還鑒定到6條僅在心臟和骨骼肌中高表達(dá)轉(zhuǎn)錄本，這些轉(zhuǎn)錄本可能和肌肉組織的收縮功能密切相關(guān)。之前人類不同組織的基因表達(dá)譜研究也顯示心臟和骨骼肌的表達(dá)譜相似性高(Shmuelietal.，2003)。此外腎臟和肝臟組織共同高表達(dá)的轉(zhuǎn)錄本有4條，這可能源于肝臟和腎臟之間有密切相關(guān)的生理功能。