基于全長轉(zhuǎn)錄組測序的牛脂肪組織可變剪接事件分析

夏晗，李凡，彭玲偉，杜玉琴，郭愛珍，楊利國，周揚(yáng)

華中農(nóng)業(yè)大學(xué)動物科學(xué)技術(shù)學(xué)院、動物醫(yī)學(xué)院/農(nóng)業(yè)動物遺傳育種與繁殖教育部重點(diǎn)實(shí)驗(yàn)室，武漢 430070

肌內(nèi)脂肪的含量以及沉積分布影響牛肉的嫩度、系水力、剪切力值、風(fēng)味和多汁性，適當(dāng)提高牛肉肌內(nèi)脂肪含量可改變牛肉品質(zhì)[1-2］。但是，當(dāng)皮下脂肪的積累較多時(shí)會降低飼料報(bào)酬，造成資源損耗[3］。脂肪沉積受復(fù)雜的調(diào)控網(wǎng)絡(luò)影響，Utrera 等[4］針對包括脂肪沉積能力的牛胴體性狀遺傳力進(jìn)行了72 項(xiàng)研究，也未能完全了解其遺傳機(jī)制。

可變剪接是剪接體對前體mRNA 選擇性的修飾，使同一基因可以表達(dá)不同成熟mRNA 剪接異構(gòu)體的現(xiàn)象[5］?？勺兗艚颖徽J(rèn)為在基因中廣泛存在，有研究顯示，在高通量測序結(jié)果中95%以上的人類基因存在可變剪接，由此產(chǎn)生表達(dá)多種不同的剪接異構(gòu)體[6］。可變剪接使基因表達(dá)具有更多可能性，展示了基因功能的多樣性[7］。研究顯示可變剪接的發(fā)生會導(dǎo)致蛋白結(jié)構(gòu)的變化，導(dǎo)致同一基因的不同蛋白亞型在功能上出現(xiàn)差異[8］；且有研究表明，可變剪接會對牛的脂肪相關(guān)基因表達(dá)產(chǎn)生重要影響[9］。例如，NOVA 依賴性的可變剪接會調(diào)節(jié)白色脂肪組織褐變[10］、質(zhì)量減輕導(dǎo)致的皮下脂肪中TCF7L2 可變剪接調(diào)節(jié)可導(dǎo)致脂肪組織中的高血糖和胰島素作用受損[11］等。所以了解可變剪接事件的發(fā)生對于明晰脂肪沉積調(diào)控網(wǎng)絡(luò)有著重要的作用。

基于PacBio 平臺的全長轉(zhuǎn)錄組測序是第三代測序技術(shù)的代表，能直接測序得到RNA 轉(zhuǎn)錄本的全長片段，可以規(guī)避短片段帶來的缺失信息的影響，提供更完整的轉(zhuǎn)錄組信息，更好地進(jìn)行差異表達(dá)基因分析及功能注釋[12］。同時(shí)，長讀技術(shù)有利于更精確、更準(zhǔn)確、更多地發(fā)現(xiàn)可變剪接事件，有助于提高對脂肪沉積調(diào)控網(wǎng)絡(luò)的認(rèn)識[13］。目前，第三代測序技術(shù)已經(jīng)被廣泛應(yīng)用于植物的遺傳信息挖掘，并取得了較好的結(jié)果，例如：赪桐的種質(zhì)資源挖掘[14］、紫娟茶樹葉片呈色機(jī)理及物質(zhì)合成途徑[15］、掌葉魚黃草的遺傳信息挖掘[16］等。但是，目前基于全長轉(zhuǎn)錄組的研究較少，尤其是牛類動物繁殖周期長、生長緩慢更是限制了其研究的進(jìn)展。本研究利用基于PacBio 測序平臺的第三代測序技術(shù)對夷陵黃牛脂肪組織遺傳信息進(jìn)行挖掘，以期為正確認(rèn)識可變剪接對脂肪沉積調(diào)控網(wǎng)絡(luò)的機(jī)制提供理論參考。

1 材料與方法

1.1 試驗(yàn)材料

供試牛為湖北省宜昌地區(qū)、年齡在3 歲、發(fā)育正常、體態(tài)良好的3 頭健康夷陵黃牛。屠宰后，使用消毒手術(shù)刀，分別在每頭牛相同部位采集5 cm3大小的皮下脂肪、腹腔脂肪、肌間脂肪樣品。使用75%乙醇徹底清洗組織，PBS 清洗殘留在組織上的乙醇，吸干水分后于?80 ℃液氮保存，用于總RNA 提取。

1.2 總RNA提取

使用TRIZOL 法提取夷陵黃牛脂肪組織總RNA。提取后的RNA 需使用1%凝膠電泳檢測樣品純凈度，并經(jīng)過Nano Drop 儀器檢測，確定總RNA的OD260/OD280值在1.9~2.1、OD260/OD230值在1.4~1.8。并使用Qubit 對RNA 進(jìn)行定量，Agilent 2100、Bioanalyzer 等進(jìn)一步進(jìn)行質(zhì)量檢測，提取的總RNA質(zhì)量滿足要求后，低溫保存，待后續(xù)進(jìn)行文庫構(gòu)建。

1.3 文庫構(gòu)建及測序

構(gòu)建文庫所使用的RNA 需要使用Qiagen 試劑盒進(jìn)行回收純化預(yù)處理，以獲得更高質(zhì)量的RNA。隨后使用SMARTer?PCR cDNA Synthesis Kit 試劑盒反轉(zhuǎn)得到cDNA。依據(jù)cDNA 質(zhì)量優(yōu)化擴(kuò)增程式，使用KAPA HiFi PCR Kits 試劑盒進(jìn)行擴(kuò)增，對獲得的全長cDNA 片段使用SMRT bell template prep kit 1.0試劑盒構(gòu)建文庫。同時(shí)，對cDNA 進(jìn)行末端修復(fù)、損傷修復(fù)處理，并在全長cDNA 的雙端添加stem loop sequencing adapter 進(jìn)行篩選，文庫構(gòu)建成功后，利用PacBio Sequel 測序平臺開展全長轉(zhuǎn)錄組測序。

1.4 全長轉(zhuǎn)錄組數(shù)據(jù)組裝與篩選

對建庫測序得到的原始數(shù)據(jù)，首先使用SMRTLINK 5.1 軟件對其進(jìn)行質(zhì)控，然后使用SNR（Signal Noise Ratio）篩選剔除低質(zhì)量數(shù)據(jù)及接頭序列，得到環(huán)形一致序列（circular consensus sequence，CCS）。去除含5' primer、3' primer、poly A 片段，獲得的全長轉(zhuǎn)錄組序列經(jīng)過校正和去冗余后，可用于后續(xù)分析。

1.5 可變剪接檢測

試驗(yàn)使用ICE（iterative isoform-clustering）中的algorithm 模塊進(jìn)行聚類，使用DAGCon 得到一致性序列，并通過arrow 進(jìn)行校正。將二代數(shù)據(jù)校正后的全長序列與參考基因組比對，去除融合基因及非注釋數(shù)據(jù)后，將比對到參考基因組同一基因上的不同isoform 進(jìn)行提取，即為可變剪接。根據(jù)可變剪接不同的形成方式進(jìn)行分類。

1.6 全長轉(zhuǎn)錄組功能注釋與分類

獲得高質(zhì)量的全長轉(zhuǎn)錄組序列后，使用KOG 蛋白質(zhì)真核同源數(shù)據(jù)庫、GO 基因本體數(shù)據(jù)庫、KEGG代謝通路數(shù)據(jù)庫、NR 非冗余蛋白數(shù)據(jù)庫、Swiss Prot精準(zhǔn)蛋白數(shù)據(jù)庫對預(yù)測的Isoforms進(jìn)行注釋，用以評估脂肪組織中基因的功能注釋信息及其參與的代謝通路及調(diào)控網(wǎng)絡(luò)，并判斷基因的物種相似性。

2 結(jié)果與分析

2.1 測序質(zhì)量分析

對采集的夷陵黃牛的皮下脂肪、腹腔脂肪、肌間脂肪樣品使用PacBio 平臺第三代測序技術(shù)進(jìn)行全長轉(zhuǎn)錄組測序，共獲得19.81 G的原始數(shù)據(jù)，篩除接頭和長度小于50 bp的原始離線數(shù)據(jù)后，得到13 014 954條Subreads，平均長度為1 474 bp（Full Passes≥1、最小預(yù)測準(zhǔn)確度0.8），核驗(yàn)下機(jī)數(shù)據(jù)中的CCS（circular consensus sequence）的質(zhì)量，共獲得779 345 個(gè)CCS reads，CCS 總堿基數(shù)為1 812 898 181 bp，CCS Read平均長度2 320.43 bp，平均Full Pass 數(shù)量為11。通過驗(yàn)證5′引物、3′引物和poly-A尾的保留，獲得高質(zhì)量的FLNC（full-length nonchimeric reads）640 205條，占比82.15%，平均長度為2 091 bp，以上結(jié)果說明全長轉(zhuǎn)錄組測序數(shù)據(jù)質(zhì)量較好，可以用于后續(xù)試驗(yàn)。

2.2 全長轉(zhuǎn)錄組Reads 聚類與校正

利用來源于相同樣品的二代測序數(shù)據(jù)對PacBio平臺的第三代測序reads 進(jìn)行校正[17］。校正數(shù)據(jù)后，我們獲得了779 345 個(gè)consensus reads。使用ICE 去除缺失序列和冗余序列后，進(jìn)一步使用Lordec 軟件，以二代數(shù)據(jù)作為參照對356 772 條一致性序列結(jié)果進(jìn)行校正，獲得可用于進(jìn)一步分析的一致性序列總長為790 043 023 bp。

2.3 可變剪接分析

校正后的數(shù)據(jù)與參考基因組進(jìn)行比較，得到與2 095 個(gè)融合基因相關(guān)的2 768 條reads。將剩余數(shù)據(jù)進(jìn)一步過濾去冗余。使用Match Annot 軟件將比對后結(jié)果與注釋信息進(jìn)行比較，并將三代注釋結(jié)果和原基因組結(jié)果進(jìn)行合并，得到基因注釋的PB iso?form 63 400 條，非基因區(qū)序列有36 480 條。將比對到參考基因組cDNA 區(qū)域的序列進(jìn)行檢測，共識別到可變剪接事件50 520 個(gè)。其中2 406 個(gè)基因發(fā)生外顯子跳躍、1 184 個(gè)基因發(fā)生外顯子替換、1 262 個(gè)基因發(fā)生外顯子接受、4 085個(gè)基因發(fā)生內(nèi)含子保留、2 365個(gè)基因發(fā)生外顯子位置替換、4 143個(gè)基因發(fā)生其他可變剪接事件（表1）。

表1 可變剪接類型統(tǒng)計(jì)表Table 1 Statistical table of alternative splicing types

將檢測到可變剪接事件的基因進(jìn)行GO 富集分析，顯示與83 個(gè)條目顯著相關(guān)（P<0.05），其中富集到分子功能類的GO 條目有22 個(gè)，細(xì)胞組分類的GO條目14 個(gè)，生物過程類的GO 條目47 個(gè)。在細(xì)胞組分類別中最顯著是“薄膜涂層”，其次是“涂層膜”“黏著斑”條目。在分子功能類別中最顯著的是“分子功能”，其次是“連接”“蛋白質(zhì)結(jié)合”條目。在生物學(xué)過程類別中最顯著的是“生物調(diào)節(jié)”，其次是“細(xì)胞過程的調(diào)節(jié)”“生物過程的調(diào)節(jié)”，另外，我們還檢測到與脂質(zhì)合成及代謝相關(guān)過程的顯著富集，包括“糖異生”“葡萄糖代謝過程”“己糖生物合成過程”“單糖生物合成過程”“糖酵解過程”“丙酮酸代謝過程”“己糖代謝過程”等。 這說明可變剪接對于脂肪沉積調(diào)控網(wǎng)絡(luò)可能有著重要的影響（圖1）。

圖1 GO及KEGG功能富集圖Fig.1 Functional enrichment of GO and KEGG

對檢測到可變剪接事件的基因進(jìn)行KEGG 富集分析，以期研究可變剪接對基因之間互作關(guān)系的影響，結(jié)果顯示有27條通路顯著富集，其中直接與脂質(zhì)合成及代謝相關(guān)的有15條通路，例如“脂肪細(xì)胞中脂解的調(diào)節(jié)”“脂肪消化吸收”“脂肪細(xì)胞因子信號通路”“PPAR 信 號 通 路”“AMPK 信 號 通 路”“PI3KAkt 信號通路”等。這說明可變剪接的參與導(dǎo)致脂肪沉積調(diào)控網(wǎng)絡(luò)的復(fù)雜化。

2.4 新基因功能注釋分析

為了獲得更直觀的注釋結(jié)果，使用NR、KOG、GO、KEGG 和Swissport 5 個(gè)數(shù)據(jù)庫進(jìn)行注釋，在全部80 756 條序列中有69 259 條被注釋，如表2 所示。如圖2 顯示至少有16 572 個(gè)基因被5 個(gè)數(shù)據(jù)庫所注釋，有26 653個(gè)基因至少在1個(gè)數(shù)據(jù)庫中被注釋。其中占全部基因數(shù)量81.61%的65 902個(gè)基因被NR 數(shù)據(jù)庫所注釋，根據(jù)序列同源性比對結(jié)果，其中28 530條與普通牛同源、10 028條與牦牛同源、3 265條與平原野牛同源、2 970 條與綿羊同源、2 431 條與亞洲水牛同源以及與其他物種同源的18 678 條，如圖3 所示。這說明了夷陵黃牛在品種形成過程中主要受到普通牛和瘤牛的影響。

圖2 功能注釋韋恩圖Fig.2 Venn diagram of functional annotation

圖3 NR分類注釋統(tǒng)計(jì)圖Fig.3 NR classification annotation statistics

表2 功能注釋統(tǒng)計(jì)表Table 2 Function note statistics table

KOG 數(shù)據(jù)庫注釋到41 004 條序列在基因區(qū)域，占全部序列的50.78%。根據(jù)序列信息可以分為26組，其中一般功能預(yù)測條目注釋到8 406 條序列，信號轉(zhuǎn)導(dǎo)機(jī)制條目7 221 條，注釋最少的為未知蛋白組（僅為52 條）。從分組可以發(fā)現(xiàn)，注釋到細(xì)胞代謝發(fā)育等通用途徑的序列占絕大多數(shù)。例如：胞內(nèi)運(yùn)輸、分泌和囊泡運(yùn)輸、脂質(zhì)運(yùn)輸和代謝、能量生產(chǎn)和轉(zhuǎn)化分別注釋到2 822、1 620、920 條序列，這些途徑證明了在夷陵黃牛脂肪組織中進(jìn)行著高頻率的能量代謝以及物質(zhì)跨膜運(yùn)輸，這可能與脂肪細(xì)胞中脂質(zhì)積累相關(guān)（圖4）。

圖4 KOG 分類注釋統(tǒng)計(jì)圖Fig.4 KOG classification annotation statistics

GO 數(shù)據(jù)庫注釋到的基因有26 653 個(gè)，占全部基因數(shù)量的33%，富集到GO條目共1 280個(gè)，在分子功能條目上富集到559 個(gè)條目，30 672 條序列，其中以蛋白質(zhì)結(jié)合、核酸結(jié)合條目聚集到最多的序列，分別為5 123 條和2 706 條；生物過程富集到507 個(gè)條目15 039條序列，其中轉(zhuǎn)錄調(diào)控及蛋白磷酸化聚集序列最多，分別為1 495 條和1 423 條；細(xì)胞組分富集到214 個(gè)條目7 907 條序列，其中膜和膜的組成部分條目聚集的序列最多，分別為1 142 條和1 721 條。其中本研究檢測到與細(xì)胞能量儲存與代謝等相關(guān)條目的顯著富集，可能與脂肪組織復(fù)雜的生物學(xué)過程及活躍的細(xì)胞代謝有關(guān)（圖5）。

圖5 GO分類注釋統(tǒng)計(jì)圖Fig.5 GO classification annotation statistics

KEGG 的功能注釋結(jié)果發(fā)現(xiàn)，夷陵黃牛脂肪組織全長轉(zhuǎn)錄組數(shù)據(jù)中共注釋到基因45 984 個(gè)，占全部基因的56.94%，KEGG 將注釋得到的基因信息分為5 個(gè)大類：細(xì)胞過程、環(huán)境信息處理、遺傳信息處理、代謝、有機(jī)系統(tǒng)，這些類別分別富集到10 745、9 516、5 789、7 005、1 655 個(gè) 基 因，如 圖6 所 示。KEGG 除富集到“細(xì)胞生長和死亡”“細(xì)胞運(yùn)動”“免疫系統(tǒng)”“神經(jīng)系統(tǒng)”等通路外，還富集到如“能量代謝”“聚糖生物合成和代謝”“脂質(zhì)代謝”等與脂肪合成、能量代謝高度相關(guān)的路徑，這些顯著富集的通路可能參與脂肪沉積網(wǎng)絡(luò)的調(diào)控。

圖6 KEGG 分類注釋統(tǒng)計(jì)圖Fig.6 KEGG classification annotation statistics

將80 756 條序列與Swiss-Prot 數(shù)據(jù)庫蛋白庫進(jìn)行比對，共注釋得到69 259 個(gè)基因，占全部基因數(shù)量的85.76%。采用angel 軟件，對序列做CDS 和Pro?tein 預(yù)測，共發(fā)現(xiàn)94 458 條CDS 序列及其對應(yīng)的pep序列，并且比對到UTR 非編碼區(qū)域的155 214 條序列。將數(shù)據(jù)與基因注釋結(jié)果對比優(yōu)化后，得到3'or 5'延長區(qū)域共29 448 處，分別對應(yīng)21 145 條reads，7 176 個(gè)基因。

3 討 論

可變剪接已經(jīng)被證實(shí)在95%以上的人類基因中存在，可變剪接的存在使得基因的功能存在不確定性，脂肪調(diào)控網(wǎng)絡(luò)變得更為復(fù)雜[7］。三代測序技術(shù)可以彌補(bǔ)二代測序序列短、完整性差的缺陷，直接獲得全長序列，從而準(zhǔn)確地表征基因結(jié)構(gòu)特點(diǎn)，同時(shí)使用RNA-seq 對全長轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行校正，在進(jìn)一步提高PacBio 循環(huán)共有序列（CCS）質(zhì)量的同時(shí)，降低長序列的高堿基錯誤率，提升可變剪接檢出率和準(zhǔn)確率[18］。本研究利用三代測序技術(shù)對夷陵黃牛脂肪組織（皮下脂肪、腹腔脂肪、肌間脂肪）的總RNA 進(jìn)行全長轉(zhuǎn)錄組測序。通過數(shù)據(jù)質(zhì)控和分析，一共確定了779 345 個(gè)CCS reads，平均長度為2 320.43 bp，共識別可變剪接事件50 520 個(gè)，超過在牛上發(fā)現(xiàn)全部可變剪接事件的33.4%，說明可變剪接可能與脂肪沉積復(fù)雜的調(diào)控網(wǎng)絡(luò)緊密聯(lián)系[19］。前5%的發(fā)生可變剪接基因GO 和KEGG 富集分析結(jié)果顯示，83 個(gè)條目顯著富集，其中與脂質(zhì)代謝相關(guān)過程的的有“糖異生”“葡萄糖代謝過程”“糖酵解過程”“丙酮酸代謝過程”“己糖代謝過程”等。KEGG 結(jié)果顯示，這些可變剪接相關(guān)基因涉及到的通路共有27 條，與脂質(zhì)合成及代謝相關(guān)的15 條，其中“PPAR 信號通路[20］”“AMPK 信 號 通 路[21］”“PI3K-Akt 信 號 通 路[22］”“MAPK 信號通路[23］”“Ras 信號通路[24］”“Rap1 信號通路[25］”“ECM-受體相互作用通路[26］”“cGMP-PKG信號通路[27］”“Apelin信號通路[28］”等均有文獻(xiàn)報(bào)道。這些可變剪接事件的發(fā)現(xiàn)，是對現(xiàn)有脂肪可變剪接數(shù)據(jù)庫的補(bǔ)充，也為充分正確認(rèn)識脂肪沉積調(diào)控網(wǎng)絡(luò)提供了新的理論依據(jù)補(bǔ)充。

本研究共發(fā)現(xiàn)80 756 個(gè)新基因序列，KOG、KEGG、NR、Swiss Prot、GO 數(shù)據(jù)庫分別注釋到其中41 004、45 984、65 902、57 880、26 653條序列，在物種同源性比對中發(fā)現(xiàn)夷陵黃牛主要與普通牛和瘤牛同源，屬于南方黃牛的血統(tǒng)組成，與夏小婷等[29］的研究結(jié)果一致。但是小眾動植物普遍存在基因組信息研究匱乏等現(xiàn)象，在進(jìn)行同源比對時(shí)可參考的基因組序列信息有限，容易導(dǎo)致比對結(jié)果的偏離[30］。盡管已經(jīng)證明了全長轉(zhuǎn)錄組在缺少基因組信息時(shí)有獲取數(shù)據(jù)的優(yōu)勢，但長讀技術(shù)依舊受到原始樣品質(zhì)量的影響[31］。本研究中，KOG 數(shù)據(jù)庫注釋到的序列占全部序列的50.78%，其中聚集到一般性功能預(yù)測分組的序列數(shù)量最多，說明維持細(xì)胞一般性功能的序列占據(jù)主要部分。GO 數(shù)據(jù)庫注釋的序列占比33%，主要與維持細(xì)胞正常形態(tài)及功能相關(guān)，其中分子功能組以蛋白質(zhì)結(jié)合、核酸結(jié)合條目聚集到最多的序列，生物過程組中轉(zhuǎn)錄調(diào)控及蛋白磷酸化聚集序列最多，細(xì)胞組分中膜的組成部分聚集的序列最多。KEGG 將注釋得到的基因信息分為5個(gè)大類：細(xì)胞過程10 745 條、環(huán)境信息處理9 516 條、遺傳信息處理5 789 條、代謝7 005 條、有機(jī)系統(tǒng)1 655 條，除富集到“細(xì)胞生長和死亡”“細(xì)胞運(yùn)動”“免疫系統(tǒng)”、“神經(jīng)系統(tǒng)”等主要通路外，還富集到如“能量代謝”“聚糖生物合成和代謝”“脂質(zhì)代謝”等與脂質(zhì)合成及代謝高度相關(guān)的路徑。除此之外，我們還使用Swiss-Prot數(shù)據(jù)庫蛋白庫注釋到69 259 條序列，發(fā)現(xiàn)94 458 條可能的CDS 序列及其對應(yīng)的pep 序列。這些數(shù)據(jù)可以為現(xiàn)有脂肪沉積相關(guān)研究進(jìn)行數(shù)據(jù)補(bǔ)充，并為后續(xù)研究提供參考信息，為后續(xù)的研究提供寶貴資源。