魏文景，孫學(xué)良，占思遠(yuǎn)，陳 瑜，張國俊，苗 斌，張紅平*，郭家中*

（1.四川農(nóng)業(yè)大學(xué)動物科技學(xué)院，成都 611130；2.巴中市南江黃羊科學(xué)研究所，四川 南江 635600）

受益于高通量測序技術(shù)的迅速發(fā)展，目前主要畜禽物種全基因組序列均已發(fā)布（http://asia.ensembl.org/index.html）[1]；盡管如此，不斷完善各物種全基因組序列的注釋才能為相關(guān)物種重要經(jīng)濟(jì)性狀遺傳分析、基因組選擇等研究提供堅實的數(shù)據(jù)基礎(chǔ)。伴隨著 Hi-C[2]、CHIP-seq[3]和 RNA-seq[4]等功能基因組學(xué)技術(shù)的發(fā)展，人們得以在基因表達(dá)的不同水平上深入注釋畜禽全基因組序列，尤其是各種調(diào)控元件[5-7]。例如，與Illumina測序平臺為代表的短讀段測序技術(shù)相比，PacBio和Nanopore平臺實現(xiàn)了高通量長讀段測序，在轉(zhuǎn)錄水平上能夠捕獲全長轉(zhuǎn)錄本[8-10]，從而更準(zhǔn)確地探究可變剪接[11-12]、多聚腺苷酸化[13-14]等問題。

山羊是重要的家畜物種之一，山羊遺傳育種的深入研究依賴于完善的山羊基因組序列及注釋。目前被廣泛使用的山羊參考基因組（ARS1）具有較高組裝質(zhì)量（缺口總長度38 Mb，contig N50值為26.2 Mb等特點）[15]；但功能注釋并不完善，這些未注釋的大量基因可能具有重要的功能。例如，本課題組利用比較基因組學(xué)分析，初步確定山羊7號染色體45.05～59.76 Mb區(qū)域是影響山羊小耳表型的一個重要位點，但該區(qū)域內(nèi)的大部分序列缺少注釋，從而阻礙了因果突變鑒定的深入開展[16]。另一方面，已報道的山羊不同發(fā)育階段、不同組織的轉(zhuǎn)錄組研究主要使用短讀段測序技術(shù)，側(cè)重于基因差異表達(dá)分析。例如，Zhan S.Y.對不同階段山羊胎兒背最長肌組織的RNA-seq進(jìn)行轉(zhuǎn)錄組圖譜分析，以識別差異表達(dá)基因[17]；C.Muriuki根據(jù)表達(dá)模式將預(yù)測轉(zhuǎn)錄本分組，并進(jìn)行聚類分析以完善山羊基因表達(dá)圖譜數(shù)據(jù)集[18]。而有關(guān)山羊基因功能注釋方面的研究偏少[17-20]。本研究收集了南江黃羊的5種組織，利用PacBio測序技術(shù)開展全長轉(zhuǎn)錄組測序和分析，旨在完善功能注釋，從而為山羊遺傳育種研究提供理論參考。

1 材料和方法

1.1 組織樣品收集

在四川省南江黃羊原種場，選擇1只約7月齡健康的雄性斷奶山羊。在禁食12 h后進(jìn)行屠宰，分別收集羊的肝、脾、腎周脂肪、背最長肌和睪丸等組織樣品，隨后立即置于液氮中，帶回實驗室保存于-80℃冰箱中。

1.2 測序文庫構(gòu)建和高通量測序

選擇肝、脾、脂肪、肌肉和睪丸共5個組織的樣品送至北京諾禾致源科技股份有限公司，使用Trizol試劑進(jìn)行總RNA提取。待RNA樣品檢測合格后，將上述5個組織總RNA進(jìn)行等量混勻，形成1個總RNA混合樣本。按照PacBio測序說明書，利用帶有Oligo(dT)的磁珠富集總RNA混合樣本中的mRNA后，進(jìn)行大規(guī)模PCR擴增并進(jìn)行各種修復(fù)，最終獲得SMRT bell文庫。使用PacBio Sequel測序儀進(jìn)行高通量測序。

1.3 序列質(zhì)控和全長轉(zhuǎn)錄本序列的鑒定

在獲得原始測序數(shù)據(jù)（polymerase reads）后，使用PacBio官方軟件SMRTlink（v.7.0）去除接頭序列并過濾掉長度小于等于50 bp的序列，獲得子序列（subreads）。對subreads進(jìn)行校正糾錯后，獲得環(huán)形一致性序列（circular consensus sequence，CCS）；再根據(jù)CCS是否包含5′端引物、3′端引物以及polyA序列分為全長非嵌合序列（full-length non-chimeric reads，F(xiàn)LNC）與非全長非嵌合序列（non-full-length reads，NFL）；然后利用同種類型聚類（ICE）對全長序列進(jìn)行聚類，獲得Cluster consensus序列；最后使用非全長序列對獲得的一致序列進(jìn)行校正；最終獲得高質(zhì)量一致性全長序列（polished consensus reads，PC），被用于后續(xù)分析。

1.4 生物信息學(xué)分析

使用GMAP[21]（v.2017-06-20）軟件（默認(rèn)參數(shù)），將高質(zhì)量一致性全長序列比對到山羊參考基因 組（ARS1，GCA_001704415.1）。 然 后 ，使 用TAPIS[12,22]（v.1.2.1）軟件對比對上基因組的PC序列進(jìn)行進(jìn)一步的校正、聚類去冗余得到轉(zhuǎn)錄本。利用GMAP軟件對FLNC序列進(jìn)行飽和度曲線分析，評估測序數(shù)據(jù)量。將轉(zhuǎn)錄本與山羊基因組比對，得到比對上的轉(zhuǎn)錄本與未比對的轉(zhuǎn)錄本。使用SUPPA[23]（v.2.3）軟件對比對上的轉(zhuǎn)錄本進(jìn)行可變剪接事件鑒定；該軟件將可變剪接分為7類：外顯子跳躍（skipped exon，SE）、外顯子互斥 (mutually exclusive exon，MX)、5′端 可 變 剪 接 (alternative 5′splice site，A5)、3′端可變剪接 (alternative 3′splice site，A3)、內(nèi)含子滯留 (retained intron，RI)、起始外顯子可變剪接(alternative first exon，AF)、終止外顯子可變剪接(alternative last exon，AL)。使用KEGG公共數(shù)據(jù)庫，對比對上與未比對到基因組的轉(zhuǎn)錄本序列進(jìn)行功能注釋。

2 結(jié)果與分析

2.1 全長轉(zhuǎn)錄本的鑒定和比對

對原始序列進(jìn)行質(zhì)控，共獲得21 078 721個subreads（36.97 Gb），subreads平均長度為 1 755 bp（N50=2 107 bp）。根據(jù)全長序列的判定標(biāo)準(zhǔn)，進(jìn)一步獲得341 023條FLNC序列（包含5′引物、3′引物和polyA尾巴），平均長度為1 940 bp（N50=2 242 bp）。對FLNC序列進(jìn)行校正及聚類去冗余，最終獲得28 432條PC序列，平均長度為2 045 bp（N50=2 398 bp）（表1）。上述PC序列被認(rèn)為是潛在的轉(zhuǎn)錄本，被用于后續(xù)的基因組比對。

表1 校正前后的轉(zhuǎn)錄本長度分布統(tǒng)計表Table 1 The summary of the length distribution of transcripts before and after correction

基因組比對結(jié)果表明，99.49%的PC序列（28 288條）被成功比對到山羊基因組，僅有144條序列未比對上（圖1A）。在比對上的序列中，93.87%的序列（26 689條）比對到基因組的唯一位置——13 532條序列（47.59%）比對到基因組“+”鏈，13 157條序列（46.28%）比對到基因組“-”鏈。另外，飽和度曲線表明，本研究獲得的全長序列能全面涵蓋相關(guān)組織表達(dá)的轉(zhuǎn)錄本（圖1B）。

圖1 PC序列的比對結(jié)果總結(jié)（A）和全長非嵌合序列的飽和度曲線（B）Figure 1 Mapping summary of PC sequences（A）and Saturation curve for full-length non chimera reads（B）

2.2 轉(zhuǎn)錄本可變剪接及新基因功能注釋分析

對成功比對到基因組的序列進(jìn)行聚類、去冗余，共鑒定出來源于9 879個基因的19 115個轉(zhuǎn)錄本，平均長度為2 069 bp（N50=2 452 bp）。發(fā)生轉(zhuǎn)錄的基因數(shù)量占山羊參考基因組（ARS1）注釋基因總數(shù)（27 271個）的36.23%。根據(jù)基因組注釋，9 879個被表達(dá)基因包括8 672個已知基因和1 207個新基因。在上述新基因中，1 120個基因位于常染色體上；其中18（65）、19（64）、3（60）和1號（60）染色體上數(shù)量較多；其余87個新基因位于Scaffold上；另外，這些新基因總共表達(dá)1 790個轉(zhuǎn)錄本。

統(tǒng)計分析結(jié)果表明，每個基因平均轉(zhuǎn)錄形成1.93個轉(zhuǎn)錄本（中位數(shù)為2個）；表達(dá)多個（≥2）轉(zhuǎn)錄本的基因（4 019）占被轉(zhuǎn)錄總基因數(shù)的40.68%，如ODF2、BAG6、ALDOA等基因在內(nèi)的114個基因形成的轉(zhuǎn)錄本數(shù)目大于或等于10個（表2）。由表3可知，共有4 089個基因發(fā)生了可變剪接，占表達(dá)基因總數(shù)的41.39%。在9 721個可變剪接事件中，AF類型占比最高（36.47%），其次是SE類型的可變剪接（25.44%），最少的是MX類型的可變剪接（1.78%）。

表2 轉(zhuǎn)錄本種類和可變剪接事件較多的5種基因Table 2 Five genes with more transcript types and alternative splicing events

表3 可變剪接事件統(tǒng)計結(jié)果Table 3 Alternative splicing events statistics results

功能注釋結(jié)果表明，在1 207個新基因中共有577個被成功注釋到KEGG數(shù)據(jù)庫（圖2A）。按照生物體系統(tǒng)（organismal systems）分類，其中23個與免疫系統(tǒng)相關(guān)，22個與消化系統(tǒng)相關(guān)，19個與內(nèi)分泌系統(tǒng)相關(guān)。依據(jù)人類疾病（human disease）分類，與各類癌癥相關(guān)的轉(zhuǎn)錄本共19個；與內(nèi)分泌和代謝相關(guān)的轉(zhuǎn)錄本共26個，其中注釋到轉(zhuǎn)錄本數(shù)目最多的基因為GP1BA。上述結(jié)果表明，相較于其他功能或表型相關(guān)基因，當(dāng)前的山羊參考基因組中關(guān)于免疫和內(nèi)分泌系統(tǒng)的注釋相對較少。

2.3 未比對轉(zhuǎn)錄本功能注釋分析

由圖3A可知，144個未比對到基因組的轉(zhuǎn)錄本長度的最小值為388 bp、最大值為6 083 bp、平均長度值為2 186 bp（中位數(shù)為1 884 bp）。GC含量分析結(jié)果表明（圖3B），這些序列GC含量最低值、最高值分別為36.17%、67.86%，平均值為56.50%。對未比對轉(zhuǎn)錄本序列與比對上轉(zhuǎn)錄本序列兩者進(jìn)行CG含量分析，得到未比對轉(zhuǎn)錄本GC含量顯著高于比對上的其他轉(zhuǎn)錄本這一結(jié)果（Wilcoxon秩和檢驗，P-value=1.12-13）。功能注釋結(jié)果表明，共有141個轉(zhuǎn)錄本序列被成功注釋到KEGG數(shù)據(jù)庫（圖2B）。按照生物體系統(tǒng)（organismal systems）分類，其中28個與免疫系統(tǒng)相關(guān)。依據(jù)人類疾?。╤uman disease）分類，與病毒性、細(xì)菌性和寄生蟲性傳染病例相關(guān)的轉(zhuǎn)錄本共79個，與免疫疾病相關(guān)的轉(zhuǎn)錄本共26個，其中GC含量最高（67.86%）的轉(zhuǎn)錄本被注釋到了NUDT14基因。

圖2 新基因轉(zhuǎn)錄本（A）和未比對轉(zhuǎn)錄本（B）的功能注釋Figure 2 Functional annotations for the transcripts of novel genes（A）and the unmapped transcripts（B）

圖3 未成功比對到基因組的轉(zhuǎn)錄本的長度（A）和GC含量（B）Figure 3 Length(A)and GC content(B)for unmapped transcripts

3 討論

山羊遺傳育種的深入研究依賴于完善的山羊基因組序列及注釋信息。目前，“黃金級”高質(zhì)量山羊基因組（ARS1）[15,24]被廣泛使用，但山羊基因組序列的功能注釋相對滯后[25]。本研究使用PacBio測序技術(shù)對山羊5個組織的混合樣本進(jìn)行轉(zhuǎn)錄組測序與分析。飽和度曲線分析結(jié)果表明，本研究獲得的測序數(shù)據(jù)全面地涵蓋了混合樣本中所表達(dá)的轉(zhuǎn)錄本。在約3萬個預(yù)測轉(zhuǎn)錄本中，超過60%是參考基因的新亞型，該結(jié)果與其他家畜物種以及人類基因組的注釋結(jié)果相似[26]，表明山羊基因組轉(zhuǎn)錄的復(fù)雜性。

本研究發(fā)現(xiàn)接近50%的基因轉(zhuǎn)錄出多個（≥2）轉(zhuǎn)錄本，表明了可變剪接在山羊基因表達(dá)過程中廣泛存在。已有研究表明，在動物體內(nèi)SE和AF兩種可變剪接類型事件占比最高[27-29]，這與本研究結(jié)果一致?？偟膩碚f，本研究結(jié)果顯著提高了每個基因平均表達(dá)的轉(zhuǎn)錄本數(shù)量（從1.53個增加到1.93個），這與在水牛中觀察到的相近(每個注釋基因表達(dá)1.91個轉(zhuǎn)錄本)[30]。但基于海福特牛32個組織全長轉(zhuǎn)錄本研究表明，被轉(zhuǎn)錄基因平均表達(dá)3.57個轉(zhuǎn)錄本[26]。造成上述結(jié)果差異主要原因是海福特牛研究包含的組織數(shù)量更多，種類更全。另一個原因是上述研究使用了不同的全長測序平臺；有研究發(fā)現(xiàn)，與PacBio平臺相比ONT平臺在轉(zhuǎn)錄本預(yù)測方面假陽性較高、可重復(fù)性相對較差[31-33]。

山羊的免疫基因是其抗病能力的遺傳學(xué)基礎(chǔ)；但由于免疫基因復(fù)雜的結(jié)構(gòu)以及測序技術(shù)的限制，目前相關(guān)研究僅在抗寄生蟲感染、血液免疫等部分領(lǐng)域有所進(jìn)展[34-35]。相較于蛋白質(zhì)編碼基因、rRNA、miRNA等其他類型的基因，目前山羊參考基因組關(guān)于免疫系統(tǒng)注釋相對匱乏。已有研究表明，較高的GC含量是導(dǎo)致免疫基因難以被鑒定的主要原因之一[36]。將脾臟組織加入混樣，旨在完善抗性基因相關(guān)注釋，最終在新基因轉(zhuǎn)錄本中鑒定到GP1BA等免疫相關(guān)基因。在未比對轉(zhuǎn)錄本功能注釋中，鑒定到免疫相關(guān)轉(zhuǎn)錄本數(shù)量最多，且GC含量顯著高于其他轉(zhuǎn)錄本。其中GC含量最高（67.86%）的轉(zhuǎn)錄本被注釋到了NUDT14基因（已被證明編碼蛋白影響病毒DNA復(fù)制[37]）。免疫基因Ig重鏈與輕鏈特殊的組合方式?jīng)Q定了漿細(xì)胞產(chǎn)生抗體的特異性；在該機制下，B細(xì)胞基因理論重排結(jié)果超過2×106個。而目前山羊參考基因組中IgV區(qū)段與IgC區(qū)段注釋共10個[38]，遠(yuǎn)小于理論值。綜上，本研究提高了山羊每個基因平均表達(dá)的轉(zhuǎn)錄本數(shù)量；改善了山羊免疫基因的功能注釋。