胥文才，藺杼陽，肖雅丹，唐瑞祥，范振鑫，2，孟楊，2*

(1.四川大學(xué)生命科學(xué)學(xué)院，生物資源與生態(tài)環(huán)境教育部重點(diǎn)實(shí)驗(yàn)室，成都610065;2.四川大學(xué)生命科學(xué)學(xué)院，四川省瀕危動(dòng)物保護(hù)生物學(xué)重點(diǎn)實(shí)驗(yàn)室，成都610065)

中華蟾蜍Bufo gargarizans隸屬于無尾目Anura蟾蜍科Bufonidae蟾蜍屬Bufo，是中藥材蟾酥和蟾衣的藥源。蟾酥是蟾蜍耳后腺分泌物，蟾衣為蟾蜍去除內(nèi)臟后的表皮(國家藥典委員會(huì)，2010)，蟾蜍毒液是由蟾蜍耳后腺和皮膚腺不同類型腺體的混合產(chǎn)物組成(吳文英等，2011)。蟾酥含有豐富的化學(xué)物質(zhì)，研究表明，從不同種類蟾蜍中分離得到100多種化學(xué)成分，包括多肽、類固醇、吲哚生物堿等，其中蟾蜍毒素、蟾毒配基、甾醇類等為主要成分(趙大洲等，2006;陳麗萍等，2011;辛秀蘭等，2012)。在傳統(tǒng)中藥中，蟾蜍毒素被用于治療各種疾病已有數(shù)百年的歷史，現(xiàn)代研究(包括實(shí)驗(yàn)和臨床試驗(yàn))也揭示和支持了其中的分子機(jī)制。目前蟾酥已被廣泛用作治療心力衰歇、口腔炎、咽喉炎、咽喉腫痛、皮膚癌等(羅展雄等，2009;彭貝等，2011;楊宏梅，陳濤，2014;孫崇峰等，2018)。

轉(zhuǎn)錄組是特定組織或細(xì)胞在某一發(fā)育階段或功能狀態(tài)下轉(zhuǎn)錄的所有RNA的集合(祁云霞等，2011)。隨著高通量測(cè)序技術(shù)的不斷成熟，轉(zhuǎn)錄組高通量測(cè)序技術(shù)(RNA-Seq)在基因組學(xué)的研究中發(fā)揮了非常重要的作用，如大熊貓Ailuropoda melanoleuca(杜聯(lián)明等，2014)、美洲大蠊Periplaneta americana(晉家正等，2018)、彭波半細(xì)毛羊Ovis aries(張立等，2020)等。這項(xiàng)技術(shù)近期已應(yīng)用于中華蟾蜍研究中，如不同海拔蟾蜍基因表達(dá)的變化(Yang et al.，2017)、氟化物對(duì)蟾蜍骨骼發(fā)育的影響(Chao et al.，2018)。但多數(shù)研究都是針對(duì)不同個(gè)體的表達(dá)差異，而少有針對(duì)同一個(gè)體的不同組織，同時(shí)為了比較組織間的差異，要排除個(gè)體間的差異。本研究利用RNA-Seq對(duì)中華蟾蜍耳后腺與其他6個(gè)組織進(jìn)行了轉(zhuǎn)錄組測(cè)序和分析，并利用多種方法對(duì)差異表達(dá)基因進(jìn)行分析，探究相關(guān)基因富集的代謝通路，研究數(shù)據(jù)能支持后續(xù)對(duì)中華蟾蜍基因表達(dá)研究和藥用價(jià)值的開發(fā)。

1 研究方法

1.1 數(shù)據(jù)來源

中華蟾蜍個(gè)體于2018年9月捕捉于四川大學(xué)江安校區(qū)內(nèi)，取其中1只雌性個(gè)體，解剖取其肝臟、脾臟、心臟、卵泡、肌肉、腹部皮膚、右耳后腺7個(gè)組織樣本各2份，并迅速放于液氮中保存。樣品通過冷鏈運(yùn)輸?shù)街Z禾致源公司(北京，中國)進(jìn)行了總RNA的提取和RNA-Seq測(cè)序。使用高通量測(cè)序平臺(tái)(Illumina NovaSeq 6000)對(duì)這7個(gè)組織的cDNA文庫進(jìn)行轉(zhuǎn)錄組測(cè)序。質(zhì)控標(biāo)準(zhǔn)主要有2個(gè)限制條件:一是原始測(cè)序數(shù)據(jù)中每一條序列中的N含量不能超過這一條序列的總堿基數(shù)量的10%，如果超過則去掉這一條序列及其對(duì)應(yīng)的雙端序列;二是原始測(cè)序數(shù)據(jù)中每一條序列的低質(zhì)量、堿基質(zhì)量值小于5的堿基數(shù)量不能超過這條序列的總堿基數(shù)量的50%，如果超過則去掉這一條序列及其對(duì)應(yīng)的雙端序列。

1.2 轉(zhuǎn)錄組組裝

采用de novo對(duì)轉(zhuǎn)錄組進(jìn)行從頭組裝，使用Trinity v2.8.3(Grabherr et al.，2011)對(duì)高質(zhì)量的序列從頭拼接成轉(zhuǎn)錄組。在使用Trinity組裝時(shí)，設(shè)置只保留片段長(zhǎng)度在300 bp以上的contig序列(Haas et al.，2013)。針對(duì)組裝的轉(zhuǎn)錄本，在 BUSCO v3.0.1(Simao et al.，2015)中以脊索動(dòng)物門的保守序列為背景進(jìn)行轉(zhuǎn)錄組的質(zhì)量評(píng)估。再用hisat2 v2.1.0(Kim et al.，2019)將質(zhì)控?cái)?shù)據(jù)(clean reads)重新比對(duì)到組裝好的轉(zhuǎn)錄本。最后使用corset v1.06(Davidson＆Oshlack，2014)采用層次聚類的方法，重新篩選Unigene。

1.3 功能注釋

使用 BLAST v2.7.0(Altschul et al.，1990)將最終得到的所有轉(zhuǎn)錄本同源比對(duì)到NCBI NR(RefSeq non-redundant proteins)(Pruitt et al.，2005)、Swiss-Prot、COG數(shù)據(jù)庫，E值設(shè)為1e-05。用與NR數(shù)據(jù)庫比對(duì)結(jié)果分析轉(zhuǎn)錄本的相似性以及與哪些物種共享最多的轉(zhuǎn)錄本，在使用BLAST進(jìn)行比對(duì)的結(jié)果中，E值表明在隨機(jī)的情況下，其他序列與目標(biāo)序列相似度大于這條顯示的序列的可能性，所以E值越小，結(jié)果越接近真實(shí)情況。再將結(jié)果導(dǎo)入到Blast2GO v5.2(Conesa et al.，2005)進(jìn)行 GO 注釋，根據(jù)BLAST的比對(duì)結(jié)果，匹配相應(yīng)的GO功能條目，最后將結(jié)果提交到WEGO(Ye et al.，2006)網(wǎng)站上，對(duì)轉(zhuǎn)錄本按照細(xì)胞組分(cellular component)、分子功能(molecular function)和生物過程(biological process)進(jìn)行功能分類統(tǒng)計(jì)。使用 KAAS(KEGG automatic annotation server)(Moriya et al.，2007)對(duì)轉(zhuǎn)錄本進(jìn)行KEGG通路分析。

1.4 基因差異表達(dá)分析

使用Bioconductor R語言包edgeR v3.24.3(Robinson et al.，2010)和TMM對(duì)原始測(cè)序數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化，得到的結(jié)果進(jìn)行差異基因表達(dá)分析，并使用Benjamini-Hochberg對(duì)所有統(tǒng)計(jì)檢驗(yàn)結(jié)果的P值進(jìn)行多重假設(shè)檢驗(yàn)校正，通過控制FDR來決定P值的閾值。差異基因篩選標(biāo)準(zhǔn)為FDR＜0.05和|log2(fold change)|≥1。

使用R和python編寫統(tǒng)計(jì)分類代碼，對(duì)不同組織間的差異表達(dá)結(jié)果進(jìn)行篩選得到共同上調(diào)或下調(diào)的差異表達(dá)基因。使用Bioconductor R語言包c(diǎn)lusterProfiler v3.8.1(Yu et al.，2012)進(jìn)行差異表達(dá)基因的GO功能富集分析。在進(jìn)行KEGG差異表達(dá)基因功能富集時(shí)，以所鑒定出的差異表達(dá)基因序列為輸入文件，使用KOBAS(Xie et al.，2011)，選取和中華蟾蜍相似的同源物種為背景基因集進(jìn)行KEGG通路富集分析。富集分析的P值設(shè)為0.05。

2 結(jié)果

2.1 測(cè)序結(jié)果

使用RNA-Seq對(duì)中華蟾蜍7個(gè)組織進(jìn)行了轉(zhuǎn)錄組測(cè)序，獲得共計(jì)161 661 926條雙端測(cè)序的原始數(shù)據(jù)，進(jìn)一步質(zhì)控后，共獲得160 236 971條高質(zhì)量的質(zhì)控?cái)?shù)據(jù)，占原始數(shù)據(jù)的99.12%，以這部分高質(zhì)量的質(zhì)控?cái)?shù)據(jù)作為后面實(shí)際使用的數(shù)據(jù)(表1)。測(cè)序數(shù)據(jù)已經(jīng)提交到NCBI(GenBank登錄號(hào):PRJNA638036)。

2.2 轉(zhuǎn)錄本組裝

使用Trinity將7個(gè)組織樣本數(shù)據(jù)全部組裝，得到270 017條轉(zhuǎn)錄本，其中最長(zhǎng)的36 170 bp，最短的 279 bp，總長(zhǎng)度 294 515 519 bp，平均長(zhǎng)度1 090.73 bp，整體N50的長(zhǎng)度為1 845 bp，GC堿基數(shù)占總堿基數(shù)的45.26%，組裝較為完整。

BUSCO結(jié)果顯示，組裝的轉(zhuǎn)錄本與保守序列完全比對(duì)成功的比例為87.7%，其中成功比對(duì)一次的比例為82.9%，成功比對(duì)多次的比例為4.8%，同時(shí)還有6.8%的序列為部分轉(zhuǎn)錄本序列比對(duì)成功，而僅有5.5%的序列沒有比對(duì)上任何保守序列。對(duì)于從頭組裝的轉(zhuǎn)錄組，比對(duì)結(jié)果表明本實(shí)驗(yàn)組裝的轉(zhuǎn)錄組是較為完整且準(zhǔn)確的，可以用于后續(xù)的分析。

表1 1只中華蟾蜍的7個(gè)組織的雙端測(cè)序樣本Table 1 Paired read samples of 7 tissues of a Bufo gargarizans

2.3 轉(zhuǎn)錄組比對(duì)

使用hisat2 v2.1.0將原始測(cè)序數(shù)據(jù)與組裝轉(zhuǎn)錄本進(jìn)行比對(duì)，每個(gè)樣本的比對(duì)率均高于85.00%，最高的是心臟(89.75%)，最低的是右耳后腺(87.60%)(表2)。比對(duì)率側(cè)面反映了組裝結(jié)果的正確性，在從頭組裝中，越高序列比對(duì)率代表了越多的序列參與了轉(zhuǎn)錄本的組裝。最后使用corest，基于對(duì)比結(jié)果，使用層次聚類的方法將Trinity組裝出來的轉(zhuǎn)錄本重新聚類，挑選每個(gè)類中最長(zhǎng)的序列作為該類的代表，獲得Unigene轉(zhuǎn)錄組，共計(jì)132 117條，用作后續(xù)的分析。

表2 高質(zhì)量序列比對(duì)到組裝的中華蟾蜍轉(zhuǎn)錄組上的結(jié)果Table 2 Results of high-quality sequence mapping to the assembled transcriptome of Bufo gargarizans

2.4 轉(zhuǎn)錄組注釋

2.4.1 轉(zhuǎn)錄組注釋結(jié)果 將重新聚類篩選Unigene之后獲得的轉(zhuǎn)錄組(后文稱轉(zhuǎn)錄組)與公共數(shù)據(jù)庫(包括 NR、Swiss-Prot、KEGG、COG 和 GO)進(jìn)行同源性比對(duì)注釋，注釋到的轉(zhuǎn)錄本數(shù)分別為46 744(35.38%)、34 988(26.48%)、11 871(8.99%)、11 292(8.54%)和12 320(9.33%)(表3)，其中被各數(shù)據(jù)庫均注釋到的轉(zhuǎn)錄本數(shù)為5 281(4.00%)，共注釋47 533(35.98%)條轉(zhuǎn)錄本。

表3 Unigene功能注釋Table 3 Functional annotation of Unigene

2.4.2 同源性比對(duì)結(jié)果 與非冗余蛋白質(zhì)數(shù)據(jù)庫(NR數(shù)據(jù)庫)的比對(duì)結(jié)果統(tǒng)計(jì)發(fā)現(xiàn)，大量轉(zhuǎn)錄本注釋為電子預(yù)測(cè)的蛋白質(zhì):整個(gè)轉(zhuǎn)錄本注釋結(jié)果的E值分布反映了整體注釋的結(jié)果，26.8%的E值分布在 0～1e-150，8.9%在1e-150～1e-100，16.6%在1e-100～1e-50，47.7%在1e-50～1e-05。

從比對(duì)到的物種分布得出(圖1)，最顯著的為:熱帶爪蟾 Xenopus tropicalis、非洲爪蟾 Xenopus laevis和高山倭蛙Nanorana parkeri，反映了近緣物種在分子層面的物種相似性，而非洲爪蟾和熱帶爪蟾為模式物種，高比對(duì)率結(jié)果支持了后續(xù)分析的基因富集步驟中將其選為模式物種。

2.4.3 GO注釋 根據(jù)NR數(shù)據(jù)庫的同源性比對(duì)結(jié)果，對(duì)轉(zhuǎn)錄本進(jìn)行GO注釋，12 320條轉(zhuǎn)錄本序列被注釋，共得到24 124個(gè)GO功能條目，平均每條序列被分配到2個(gè)GO功能條目。GO功能條目主要分布在生物過程(8 272，34.28%)、分子功能(10 549，43.72%)和細(xì)胞組分(5 303，21.98%)。使用WEGO對(duì)其功能條目進(jìn)行了分類(圖2)。

在生物過程中，細(xì)胞過程(cellular process)和代謝過程(metabolic process)功能是主要被注釋到的條目;在分子功能中，結(jié)合(binding)是最有代表性的條目;而在細(xì)胞組分中，細(xì)胞(cell)和細(xì)胞組分(cell part)最多。

圖1 與NR數(shù)據(jù)庫比對(duì)結(jié)果的物種分布Fig.1 Species distribution compared with the NR database

圖2 轉(zhuǎn)錄本分配到的GO注釋分類Fig.2 Classification of GO annotations to which transcripts are assigned

2.5 差異基因分析

2.5.1 差異基因的篩選 使用edgeR包依據(jù)差異表達(dá)倍數(shù)|log2(fold change)|≥1以及數(shù)據(jù)的檢驗(yàn)值FDR＜0.05作為篩選條件，篩選得到了對(duì)應(yīng)符合要求的6組差異表達(dá)基因，分別來自腹部皮膚(圖 3:A)、肝臟(圖 3:B)、肌肉(圖 3:C)、卵泡(圖3:D)、脾臟(圖3:E)、心臟(圖3:F)與右耳后腺的比較。

根據(jù)右耳后腺與其他各個(gè)組織轉(zhuǎn)錄組差異分析的結(jié)果，得出在右耳后腺中差異表達(dá)上調(diào)和下調(diào)的基因以及數(shù)目(表4)。

因?yàn)榻M織特異性的存在，不同組織之間必然存

圖3 差異表達(dá)基因log2倍數(shù)變化分布(紅點(diǎn)為最顯著的1 000個(gè)基因)Fig.3 Distribution of log2(fold change)of differentially expressed genes(red dots are the most significant 1 000 genes)FC.倍數(shù)變化 fold change，CPM．每百萬的數(shù)目counts per million

表4 右耳后腺組織與其他各組織差異表達(dá)基因個(gè)數(shù)Table 4 Number of differentially expressed genes in the retroauricular gland and other tissues

在差異表達(dá)。為了驗(yàn)證耳后腺組織的特異性表達(dá)，提取各個(gè)組織與右耳后腺之間的表達(dá)差異上調(diào)和下調(diào)的基因并獲得彼此之前的交集，得到了對(duì)應(yīng)的各個(gè)組織的差異表達(dá)基因Venn圖。與其他組織相比，在右耳后腺中共同上調(diào)的基因?yàn)? 668個(gè)(圖4:左)、共同下調(diào)的基因?yàn)?60個(gè)(圖4:右)。

2.5.2 差異基因的富集 為深入了解差異表達(dá)基因所起的生物學(xué)作用，針對(duì)這3 668個(gè)上調(diào)基因和260個(gè)下調(diào)基因，共計(jì)3 928個(gè)顯著表達(dá)的差異基因，分別做了GO通路和KEGG通路的富集分析。

對(duì)于GO功能富集分析，以中華蟾蜍所有基因的GO注釋結(jié)果為背景基因集，分別以上調(diào)差異基因和下調(diào)差異基因?yàn)槟繕?biāo)基因集進(jìn)行富集，在目標(biāo)基因集中上調(diào)和下調(diào)分別有347個(gè)和42個(gè)基因有相應(yīng)的GO注釋，差異基因在注釋到的背景基因集中覆蓋率只有10%，這在非模式生物中情況較普遍。GO的富集結(jié)果以P-adjust＜0.05為篩選閾值(以BH法校正P值得到P-adjust)，P-adjust＜0.05為顯著表達(dá)，對(duì)于這些差異表達(dá)基因，上調(diào)基因顯著富集到了78個(gè)GO條目(圖5)。

在分子功能中，共富集到28個(gè)條目，最顯著的是肽酶抑制劑活性(GO:0030414)和肽酶調(diào)節(jié)物活性(GO:0061134)，P-adjust＜1e-15，其次為酶的抑制物活性(GO:0004857，P-adjust=1.11e-13)、角質(zhì)層的結(jié)構(gòu)組成(GO:0042302，P-adjust=5.10e-11)、結(jié)構(gòu)分子活性(GO:0005198，P-adjust=1.25e-09)等(圖6)。

在細(xì)胞組分中，共富集到9個(gè)條目(圖7)，最顯著的條目為胞外區(qū)域(GO:0005576，P-adjust=4.08e-06)、中間纖維 (GO:0005882，P-adjust=4.08e-06)、中間纖維細(xì)胞骨架(GO:0045111，P-adjust=4.08e-06)、角蛋白纖維(GO:0045095，P-adjust=2.15e-05)、膠原蛋白三聚物(GO:0005581，P-adjust=2.15e-05)。

圖4 組織差異表達(dá)上調(diào)基因(左)和下調(diào)基因(右)的Venn圖Fig.4 Venn diagrams of tissue differential expression up-regulated genes(left)and down-regulated genes(right)

圖5 差異表達(dá)基因富集到GO功能條目及分類Fig.5 Enrichment of differentially expressed genes into GO function entries and classification

在生物過程中，共富集到40個(gè)條目(圖8)，最顯著的條目中肽酶活性的負(fù)調(diào)控(GO:0010466)、蛋白水解的負(fù)調(diào)節(jié)(GO:0045861)、肽酶活性的調(diào)節(jié)(GO:0052547)P-adjust都小于1e-14，其次是水解酶活性的負(fù)調(diào)節(jié)(GO:0051346，P-adjust=3.49e-15)，調(diào)節(jié)蛋白水解作用(GO:0030162，P-adjust=7.32e-15)、分子功能的負(fù)調(diào)節(jié)(GO:0044092，P-adjust=9.37e-14)、催化活性的負(fù)調(diào)節(jié)(GO:0044092，P-adjust=9.37e-14)、細(xì)胞蛋白代謝過程的負(fù)調(diào)控(GO:0043086，P-adjust=1.76e-13)、蛋白質(zhì)代謝過程的負(fù)調(diào)控(GO:0051248，P-adjust=2.01e-13)。

下調(diào)基因富集到3個(gè)顯著的條目(圖9)，均在細(xì)胞組分中，包括網(wǎng)格蛋白接頭復(fù)合物(GO:0030131，P-adjust=0.018 770 6)、AP 型膜外套接頭復(fù)合物(GO:0030119，P-adjust=0.018 770 6)和網(wǎng)格蛋白外套(GO:0030118，P-adjust=0.020 008 4)。

圖6 上調(diào)差異表達(dá)基因富集到分子功能GO功能條目Fig.6 Up-regulated differentially expressed genes enriched to molecular function GO function entry

圖7 上調(diào)差異表達(dá)基因富集到細(xì)胞組分GO功能條目Fig.7 Up-regulated differentially expressed genes enriched into cell components GO function entry

為進(jìn)一步了解差異表達(dá)基因所參與的通路，將中華蟾蜍差異表達(dá)基因進(jìn)行了KEGG通路富集分析，用篩選出的差異表達(dá)基因序列，以中華蟾蜍的近緣物種非洲爪蟾和熱帶爪蟾的基因作為背景基因集，檢測(cè)這些差異表達(dá)基因所富集的 KEGG通路。

圖8 上調(diào)差異表達(dá)基因富集到生物過程GO功能條目Fig.8 Up-regulated differentially expressed genes enriched into biological process GO function entries

KEGG富集分析結(jié)果發(fā)現(xiàn)，差異表達(dá)基因顯著富集到40個(gè)通路上，分為6個(gè)大類，分別是新陳代謝相關(guān)(metabolism)、人類疾病(human diseases)、生物系統(tǒng)(organismal systems)、環(huán)境信息過程(environment information processing)、細(xì)胞過程(celluar processes)和基因信息過程(genetic information processing)。

有大量差異表達(dá)基因富集到的通路主要為新陳代謝過程中的全局和概覽通路(ko01100)、脂類代謝(ko00600)、糖的生物合成和代謝通路(ko00010);人類疾病相關(guān)的癌癥通路(ko05200)、傳染病通路(ko05130);生物系統(tǒng)中的內(nèi)分泌系統(tǒng)相關(guān)通路(ko04915)、免疫系統(tǒng)相關(guān)通路(ko05321);環(huán)境信息過程中的信號(hào)轉(zhuǎn)導(dǎo)通路(ko04012)、信號(hào)分子與相互作用(ko04060);細(xì)胞過程中的運(yùn)輸和分解代謝通路(ko04142);基因信息過程中的翻譯相關(guān)通路(ko03010)。

圖9 下調(diào)差異表達(dá)基因富集到細(xì)胞組分GO功能條目Fig.9 Down-regulated differentially expressed genes enriched into cell components GO function entry

圖10 差異表達(dá)基因富集的最顯著20個(gè)KEGG通路氣泡圖Fig.10 The most significant 20 KEGG pathway'bubble diagrams of differentially expressed gene enrichment

針對(duì)這些基因富集到的通路，按照富集度和顯著程度，挑選了其中最顯著的前20個(gè)通路進(jìn)行研究(圖10)?？梢钥吹接写罅坎町惐磉_(dá)基因富集到了核糖體通路上，富集程度較高的通路有肝細(xì)胞癌、腫瘤中MicroRNAs等與腫瘤相關(guān)的通路，而顯著性較高的包括藥物代謝、硝基甲苯降解、咖啡因代謝、鞘脂類的生物合成通路、糖磷脂生物合成以及黏蛋白型O-聚糖生物合成的相關(guān)通路。在KOBAS 3.0中，提交中華蟾蜍的差異表達(dá)基因序列，并分別以熱帶爪蟾和非洲爪蟾的基因作為背景基因集進(jìn)行了KEGG富集(圖11、圖12)，按照顯著性P值排列，并篩選了P-adjust＜0.05的結(jié)果。其中顯著性最高為核糖體的相關(guān)通路(P-adjust=6.56e-10)和鞘糖脂的生物合成通路(P-adjust=4.07e-06)，其他富集到的包括酪氨酸、亞油酸、煙酸和煙酰胺、花生四烯酸、鞘脂類等化合物代謝合成相關(guān)通路，以及苯丙氨酸、酪氨酸和色氨酸的生物合成，纈氨酸、亮氨酸和異亮氨酸的降解通路，以及蟾酥的主要成分甾族化合物(辛秀蘭等，2012)的生物合成(P-adjust=0.022 477 894)代謝通路。

圖11 差異表達(dá)基因富集到熱帶爪蟾的結(jié)果Fig.11 Enrichment of differentially expressed genes into Xenopus tropicalis

圖12 差異表達(dá)基因富集到非洲爪蟾的結(jié)果Fig.12 Enrichment of differentially expressed genes into Xenopus laevis

3 討論

以往研究表明了中華蟾蜍耳后腺特異分泌物——蟾酥的主要成分，并從臨床醫(yī)學(xué)、藥學(xué)等角度肯定了其具有重要作用。本研究中使用RNASeq對(duì)蟾蜍多組織的轉(zhuǎn)錄組進(jìn)行測(cè)序分析，得到了第一個(gè)由7個(gè)組織組成的較完整的中華蟾蜍轉(zhuǎn)錄組數(shù)據(jù)?；诖藬?shù)據(jù)進(jìn)行篩選獲得了代表其單獨(dú)基因的轉(zhuǎn)錄本，并對(duì)此轉(zhuǎn)錄組進(jìn)行功能注釋，完成了無參考基因組的中華蟾蜍轉(zhuǎn)錄組注釋工作。對(duì)與NR數(shù)據(jù)庫的比對(duì)結(jié)果進(jìn)行了詳細(xì)統(tǒng)計(jì)，發(fā)現(xiàn)大量被注釋到的轉(zhuǎn)錄本比對(duì)到了電子注釋預(yù)測(cè)的蛋白質(zhì)上，這是由于目前兩棲動(dòng)物的研究尚未深入，有大量蛋白質(zhì)的功能尚未經(jīng)試驗(yàn)驗(yàn)證，目前僅能依據(jù)同源的蛋白質(zhì)具有相似的功能來推斷新蛋白質(zhì)的潛在功能。

GO功能富集中，分子功能大致可分為2類，第一類為酶的調(diào)節(jié)相關(guān)通路，包括肽酶抑制物活性、肽酶調(diào)控分子活性、酶的抑制物活性、酶的調(diào)控分子活性、肽鏈內(nèi)切酶調(diào)控分子活性、肽鏈內(nèi)切酶抑制物活性、蛋白質(zhì)二硫鍵異構(gòu)酶活性等，從生物信息學(xué)的角度上驗(yàn)證了耳后腺組織存在大量調(diào)控和轉(zhuǎn)錄分泌蛋白活性及其調(diào)節(jié)物的基因的表達(dá);第二類為其他一些組織特異性相關(guān)的通路，包括角質(zhì)層的結(jié)構(gòu)組成、結(jié)構(gòu)分子活性、氧運(yùn)輸活性、用于轉(zhuǎn)移糖基的轉(zhuǎn)移酶活性、用于轉(zhuǎn)移氨?；霓D(zhuǎn)移酶活性，這些通路來源于中華大蟾蜍耳后腺組織的特異性和對(duì)應(yīng)的相關(guān)功能通路。

細(xì)胞組分多顯著富集于細(xì)胞骨架、中間纖維、角蛋白纖維和胞外區(qū)域條目，而其中角蛋白纖維、中間纖維、高分子細(xì)胞骨架纖維、超分子纖維是屬于同一條GO通路，這2條通路與細(xì)胞骨架有關(guān)，推測(cè)在細(xì)胞物質(zhì)運(yùn)輸中，各類小泡和細(xì)胞器可沿著細(xì)胞骨架定向轉(zhuǎn)運(yùn)，因此推測(cè)在中華蟾蜍的耳后腺組織中激活了相應(yīng)基因調(diào)控細(xì)胞骨架功能以更好完成其耳后腺的分泌功能。

生物過程顯著富集到的通路大部分是肽酶活性的負(fù)調(diào)控、蛋白水解的負(fù)調(diào)節(jié)、肽酶活性的調(diào)節(jié)、水解酶活性的負(fù)調(diào)控、蛋白水解作用的調(diào)控、蛋白代謝過程的負(fù)調(diào)控、大分子代謝過程的負(fù)調(diào)控、細(xì)胞代謝過程的負(fù)調(diào)控、水解酶活性的調(diào)節(jié)、代謝過程的負(fù)調(diào)控，而這些通路多是與蛋白質(zhì)水解、代謝相關(guān)過程的調(diào)控有關(guān)，且負(fù)調(diào)控過程居多。猜測(cè)中華蟾蜍耳后腺組織，除了需要維持著一種蛋白質(zhì)大量分泌的環(huán)境，還需要儲(chǔ)存其分泌的各類化學(xué)物質(zhì)，因此雖然存在大量差異表達(dá)基因富集于蛋白質(zhì)的合成、分泌與運(yùn)輸相關(guān)的通路，但也存在部分差異表達(dá)基因富集在維持耳后腺內(nèi)分泌蛋白質(zhì)穩(wěn)定的相關(guān)功能，這更進(jìn)一步表明中華蟾蜍耳后腺組織存在一定的機(jī)制維持這些化合物的分泌與穩(wěn)定。在生物體中，酶參與藥物生物轉(zhuǎn)化，其代謝產(chǎn)物的活性通常低于母體藥物或不活躍，但一些生物轉(zhuǎn)化的產(chǎn)物可能具有增強(qiáng)的毒性作用，因此，生物轉(zhuǎn)化一般包含“脫毒”和“中毒”的過程，決定生物體處理藥物和化學(xué)物質(zhì)能力的主要酶系統(tǒng)之一是細(xì)胞色素P450單加氧酶，其他酶系統(tǒng)包括脫氫酶、氧化酶、酯酶、還原酶，以及一些結(jié)合酶系統(tǒng)，包括葡萄糖醛酸轉(zhuǎn)移酶、硫代轉(zhuǎn)移酶、谷胱甘肽s-轉(zhuǎn)移酶等也參與其中(Meyer et al.，1996)。這也能夠解釋在此次分析中，中華蟾蜍耳后腺組織的差異表達(dá)基因顯著富集到了大量蛋白質(zhì)合成、酶的調(diào)節(jié)相關(guān)、各類分子跨膜轉(zhuǎn)運(yùn)體以及各類前體物質(zhì)的代謝合成通路，結(jié)果顯示富集顯著度最高的通路為藥物代謝中的其他酶通路，其次為硝基甲苯降解，4-硝基甲苯為農(nóng)藥的中間體，推測(cè)這與農(nóng)藥污染相關(guān)，而中華蟾蜍有相應(yīng)的降解硝基甲苯的適應(yīng)性進(jìn)化，這也與生物轉(zhuǎn)化有一定的關(guān)聯(lián)，環(huán)境和遺傳因素導(dǎo)致藥物代謝的個(gè)體間和個(gè)體內(nèi)差異，并可能改變中毒反應(yīng)和解毒反應(yīng)之間的平衡。酪氨酸、亞油酸、煙酸和煙酰胺、花生四烯酸、鞘脂類等化合物作為蟾蜍分泌物蟾酥中各種物質(zhì)的前體物質(zhì)、修飾物質(zhì)，其對(duì)應(yīng)的相關(guān)代謝通路鞘糖脂的生物合成通路及核糖體相關(guān)通路出現(xiàn)在結(jié)果中，可以看出中華蟾蜍耳后腺作為外分泌腺體，從生物信息學(xué)分析角度上也證實(shí)其組織細(xì)胞中與代謝合成相關(guān)的基因高度活躍。特別注意的是，研究表明蟾酥的主要化學(xué)成分為各類甾族化合物，作為生物體內(nèi)類固醇一類的化合物，是蟾酥的代表性成分(辛秀蘭等，2012)，這在差異基因富集分析中也得到一定的驗(yàn)證，例如，顯著富集的甾類激素生物合成通路，而富集在這條通路上的基因，主要包括了羥基類固醇脫氫酶Ⅱ、細(xì)胞色素P450家族的各種基因。其中細(xì)胞色素P450作為一類末端加氧酶，能夠參與生物體內(nèi)的甾醇類激素合成過程(王斌，李德遠(yuǎn)，2009)。這些結(jié)果對(duì)于分析顯著富集通路中差異表達(dá)基因，以及后續(xù)驗(yàn)證有重要意義。這些數(shù)據(jù)也為中華蟾蜍功能基因組分析和耳后腺各類化合物分泌相關(guān)研究提供了有價(jià)值的資源，有助于完善中華蟾蜍產(chǎn)物蟾酥、蟾衣的藥用價(jià)值的研究。本研究只是一次對(duì)同一個(gè)體轉(zhuǎn)錄組差異的探索，可從不同個(gè)體方向進(jìn)行下一階段研究。