牦牛皺胃全長轉錄組測序分析

謝書瓊，馬士龍，唐 嬌，劉益麗，江明鋒

(1.西南民族大學,青藏高原動物遺傳資源保護與利用教育部和四川省重點實驗室，成都 610041；2.西南民族大學青藏高原研究院，成都 610041；3.西南民族大學畜牧獸醫學院，成都 610041)

牦牛(Bosgrunniens)是高寒地區特有牛種，主要分布于青海、西藏、甘肅西北部、四川西部、新疆南部和云南香格里拉等地[1]。牦牛素有“高原之舟”的美稱，一直以來為高原地區的牧民提供了吃、住、行等有效生活保障和經濟來源[2]。近些年，有研究發現牦牛在長期自然選擇過程中為適應高海拔、高寒低氧獨特的生存環境，其感官知覺、能量代謝和缺氧相關的組織器官和基因調控都出現了明顯的特異性[3]。鑒于此，國內外學者應用RNA-Seq二代轉錄組測序對牦牛心臟、肝臟、脾臟、肺臟、腎臟、大腦[4]、肌肉[5]、睪丸[6]、卵巢[7]等多個組織器官進行研究發現，心臟和肺臟是適應高原缺氧轉錄變化的關鍵器官；對于營養代謝器官復胃的研究中絕大多數傾向于牦牛瘤胃微生物群落的分析；而關于牦牛皺胃研究較少。目前，對牦牛皺胃的研究集中在其形態和組織學結構的研究[8-9]，如魯延剛[10]對牦牛皺胃組織結構的特點及其黏膜免疫相關細胞的動態分布進行了研究，發現牦牛皺胃黏膜中存在比其他非高原反芻動物更豐富的淋巴組織，使得牦牛有更強的黏膜免疫功能和抵抗外來病原侵襲的能力。

牦牛皺胃不僅能分泌胃酸以消化食物中的蛋白質和脂肪，還分泌具有抗酸、抗胃蛋白酶能力的溶菌酶[11-12]，以消化從瘤胃和網胃進入皺胃的微生物。皺胃的功能實現是一個復雜的過程，涉及到大量基因的轉錄調控，因此探究牦牛皺胃分子機制需要進行轉錄組研究。二代測序技術雖然通量很高，成本低廉，但是測序時間長、讀長短、末端質量差讓其應用有所局限[13]。隨著高通量技術的發展，以PacBio測序為代表的第3代基因測序技術逐漸興起，美國太平洋生物科學公司研發的單分子實時(single molecule real time，SMRT)測序技術應用了邊合成邊測序的原理，并以SMRT芯片為測序載體[14]。PacBio Seque1系統相比前代RS Ⅱ系列測序儀，不僅大大簡化了建庫步驟，在測序通量上也大為提高，周期更短，測序成本更低[15]。本研究應用PacBio SMRT技術對牦牛皺胃全長轉錄組進行分析，以期為后續深入了解牦牛皺胃生物功能、分子機制及相關功能基因提供理論支持和數據支撐。

1 材料與方法

1.1 樣品采集

試驗選用四川省阿壩藏族羌族自治州紅原縣的1頭4歲健康麥洼公牦牛，飼養方式為傳統放牧；試驗動物集中屠宰后，用無菌手術剪剪下小塊牦牛皺胃組織，DEPC水沖洗后，對材料進行編號，并分裝于特定的2 mL凍存管中，立即投入液氮中保存。

1.2 主要試劑與儀器

mirVanaTMmiRNA ISOlation Kit(Ambion-1561)購自賽默飛世爾科技(中國)有限公司；SMARTer PCR cDNA Synthesis Kit購自寶日醫生物技術有限公司；SMRTbellTMTemplate Prep Kit 1.0-SPv3、PacBio Sequel Sequencing Kit 3.0測序系統均購自太平洋生物科學公司。

1.3 方法

1.3.1 總RNA提取和高通量測序 使用mirVanaTMmiRNA ISOlation Kit提取樣品總RNA，1.0%瓊脂糖凝膠電泳檢測RNA降解程度及污染；利用NanoDrop檢測RNA的純度(D260 nm/D280 nm值)；采用Qubit對RNA濃度進行精確定量；利用Agilent 2100 Bioanalyzer精確檢測RNA的完整性。使用SMARTer PCR cDNA Synthesis Kit將總RNA反轉錄為cDNA，再利用SMRTbellTMTemplate Prep Kit 1.0-SPv3構建文庫，使用PacBio Sequel Sequencing Kit 3.0測序系統進行測序，得到樣本雙端數據。牦牛皺胃全長轉錄組測序工作主要由歐易生物醫學科技有限公司協助完成。

1.3.2 全長轉錄組數據處理和校正評估 對庫檢合格文庫運用PacBio Sequel系統進行測序，PacBio SMRT為單條核苷酸鏈多次滾環測序，測序過程中產出的Polymerase reads去掉接頭(adapter)序列后的有效插入片段稱為子序列(subreads)，三代測序的原始子序列數據(raw subreads)為bam格式序列。利用SMRT Link v 6.0.0[16]分析包中的Isoseq3[17]分析流程對原始子序列進行質控分析。Isoseq3分析步驟包括：①通過分析原始序列數據進行相互比對和質量校正，得到循環一致性序列(circular consensus sequence,CCS)；②利用Isoseq3調用去接頭軟件lima以去除測序引物和標簽序列(barcodes)；③利用Isoseq3程序的聚類(cluster)以去除poly A尾和嵌合體結構，將相似的序列進行聚類，生成全長非嵌合(full-length non-chimera,FLNC)的一致性序列；④利用Isoseq3的拋光(polish)命令對FLNC序列進行打磨糾錯產生非冗余的一致性序列，根據準確度閾值分為高質量(high quality，HQ)和低質量 (low quality，LQ)的Isoform序列。

將獲得的高質量Isoform序列利用SQANTI[18]軟件對非冗余轉錄本進行質控和結果注釋優化，SQANTI軟件分析主要包括：參考基因組比對、去冗余分析、SQANTI質控分析、SQANTI過濾分析、新基因預測和基因結構優化。 通過軟件minimap2[19]將高質量Isoforms序列與參考基因組進行比對。將比對結果通過Cupcake ToFU在identify=0.85下進一步去冗余得到非冗余轉錄本序列。利用SQANTI軟件的sqanti_filter2.py程序包對校正后的非冗余轉錄本進一步過濾，以刪除一些建庫測序中產生的假陽性序列信息，獲取過濾后的Unigenes。

1.3.3 轉錄組序列功能注釋 采用Diamond和HMMER3等軟件進行比對，對牦牛皺胃全長轉錄組進行功能注釋，從而得到非冗余數據庫(non-redundant protein database,NR)、蛋白質真核同源(clusters of eukaryotic orthologous groups,KOG)、基因本體論(gene ontology,GO)、蛋白質序列(swiss-prot protein database,Swiss-Prot)、直系同源蛋白分組比對(evolutionary genealogy of genes：Non-supervised Orthologous Groups,eggNOG)、東京基因與基金組百科全書(kyoto encyclopedia of genes and genomes,KEGG)、Pfam數據庫等數據庫注釋。

1.3.4 轉錄組序列的生物信息學分析 將非冗余轉錄本序列比對(BLAST)到動物轉錄因子數據庫(animal transcription factor database,AnimalTFDB)，選取比對最好的一條作為該非冗余轉錄本的轉錄因子(transcription factor,TF)家族注釋信息。用軟件ESTScan[20]預測其編碼區，得到其編碼序列(CDS)(序列方向5′→3′)和氨基酸序列。使用MISA[21]和Primer 3.0[22]軟件進行簡單重復序列(SSR)分析，對單核苷酸最少重復次數為10；二核苷酸最少重復次數為6；三核苷酸、四核苷酸、五核苷酸和六核苷酸的最少重復次數均為5，獲得牦牛皺胃轉錄組SSR位點。通過軟件Astalavista[23]檢測樣本中存在的可變剪切 (alternative splicing,AS)。

2 結 果

2.1 全長轉錄組數據統計

PacBio SMRT測序統計結果顯示，牦牛皺胃全長轉錄組共獲得14 467 420條子序列，總堿基數為48 382 430 732 bp，最長序列長度為270 809 bp，最短序列長度為51 bp，子序列平均序列長度為3 344.23 bp(表1)。

原始子序列質量低，冗余程度高,且存在嵌合體。首先，通過質控分析得到CCS reads 296 840條，帶有5′-和3′-引物接頭的全長序列(FL)為278 053條，帶有5′-和3′-引物接頭且不含嵌合體的全長序列277 712條，包含5′-引物、3′-引物和3′-端poly A尾且不存在嵌合體序列的FLNC序列277 402條，占CCS序列的93.45%。對得到的FLNC序列進行聚類分析，然后借助arrow算法進行打磨糾錯分析，最終獲得高質量的Isoform序列。質控得到高質量序列15 669條。 經過SQANTI分析去冗余后得到非冗余轉錄本序列11 104條，進一步對校正后的非冗余轉錄本過濾以刪除一些建庫測序中產生的假陽性序列信息，獲取過濾后的8 556條Unigenes。

2.2 非冗余轉錄本數據庫注釋

將獲得的8 556條Unigenes通過7個數據庫進行注釋，與NR庫比對，共有8 544條(99.86%)Unigenes被注釋，占比最多；其次是eggNOG數據庫，共有8 491條(99.24%)Unigenes被注釋；得到注釋結果數量最少的是KEGG數據庫，只有1 721條(20.11%)；而在Swiss-Prot、KOG、GO、Pfam數據庫中，分別有8 475(99.05%)、6 572(76.81%)、7 725(90.29%)、8 162(95.40%)條Unigenes被注釋(表2)。

表2 各數據庫注釋統計表

2.2.1 NR注釋庫 在NR數據庫共有8 544條Unigenes被注釋，其中有34.7%被注釋到野牦牛(Bosmutus)上，22.93%被注釋到黃牛(Bostaurus)上，6.38%被注釋到美洲草原野牛(Bisonbisonbison)上，有3.91%、2.97%、2.45%、1.52%、1.25%、1.17%、1.17%的Isoforms分別在瘤牛×黃牛(Bosindicus×Bostaurus)、亞洲水牛(Bubalusbubalis)、瘤牛(Bosindicus)、中華白海豚(Sousachinensis)、綿羊(Ovisaries)、智人(Homosapiens)、東歐馬鹿(Cervuselaphushippelaphus)上被注釋(圖1)。

圖1 NR庫前十物種分布圖Fig.1 Distribution map of top 10 species in NR library

2.2.2 KOG數據庫注釋 在KOG數據庫中有6 572條Unigenes被注釋，分別被注釋到25個分類中，其中在一般功能預測中被注釋的Unigenes最多(1 044條)，其次是翻譯后修飾、蛋白折疊和分子伴侶(880條)，信號轉導機制(816條)，細胞內運輸、分泌和囊泡轉運(555條)，而最少的是細胞運動性功能(17條)(圖2)。

2.2.3 GO功能注釋 在GO數據庫中共有7 725條Unigenes被注釋，可分為三大類：生物過程、細胞組分、分子功能。3個大類又注釋到52個小類。在生物過程中的23個小類里以細胞過程(6 276條)、單生物過程(5 254條)、代謝過程(4 583條)、生物調節(4 045條)、生物過程調節(3 754條)定位的居多，細胞殺傷功能定位最少(62條)；細胞組分中以細胞(7 056條)、細胞部分(7 049條)、細胞器(5 906條)、細胞器部分(4 325條)、膜結構(3 598條)定位居多，以細胞外基質(35條)、細胞外基質部分(35條)定位最少；在分子功能中以結合活性(5 091條)、催化活性(2 918條)居多，以通道調節器活性(4條)、受體調節活性(4條)最少(圖3)。

圖2 KOG功能分類圖Fig.2 KOG functional classification diagram

圖3 GO功能分類圖Fig.3 GO function classification diagram

2.2.4 KEGG功能注釋 在KEGG數據庫中共有1 721條Unigenes被注釋，包含了細胞過程、環境信息加工、遺傳信息加工、人類疾病、新陳代謝、生物體系統6條一級通路和42條二級通路。其中基因數量多的有遺傳信息加工里的翻譯(284條)、折疊分類和降解(176條)、轉錄(116條)；細胞過程中的運輸和分解代謝(153條)居多；生物體系統中的消化系統(149條)居多，這可能與皺胃是牦牛的消化器官有關。基因數量最少的是傳染病：寄生蟲病(1條)(圖4)。

圖4 Unigenes的KEGG功能分類Fig.4 KEGG functional classification of Unigenes

2.3 轉錄因子注釋分析

根據注釋結果，本次共注釋得到的943個轉錄因子，分布在51個轉錄因子家族。其中轉錄因子較多的家族有zf-C2H2(157個)、bHLH(76個)、Homeobox(67個)、MYB(62個)，而CSL、E2F、HTH、NCU-G1、NF-YA、Nrf1、P53、RFX、zf-BED家族轉錄因子最少，均只有1個(圖5)。

2.4 編碼區預測

利用軟件ESTScan預測其蛋白編碼區，結果顯示共有8 544個基因片段可視為蛋白編碼區，片段大小主要集中在400～1 600 bp(圖6)。

2.5 SSR分析

SSR分析結果表明，共檢索出3 596個SSR位點；其中，單核苷酸重復占大多數(67.38%)，雙核苷酸重復、三核苷酸重復、四核苷酸重復、五核苷酸重復和六核苷酸重復分別占11.62%、20.22%、0.42%、0.17%和0.19%(圖7)。

2.6 可變剪接

通過可變剪接類型分析,牦牛皺胃轉錄本共有1 825個可變剪接事件，其中外顯子跳躍最多(876個)，其次是內含子保留(382個)，外顯子互斥最少(88個)(圖8)。

圖5 轉錄因子家族Unigenes分布圖Fig.5 Distribution map of transcription factor family Unigenes

圖6 編碼區序列長度分布圖Fig.6 CDS sequence length distribution diagram

Mono-，單核苷酸重復；Di-，雙核苷酸重復；Tri-，三核苷酸重復；Tetra-，四核苷酸重復；Penta-，五核苷酸重復；Hexa-，六核苷酸重復；5～8/9～12/13～16/17～20/21～24，重復單元的重復次數Mono-，Mono nucleotide motifs;Di-，Di nucleotide motifs;Tri-，Tri nucleotides;Tetra-，Tetra nucleotides;Penta-，Penta nucleotides;Hexa-，Hexa-nucleotide motifs;5-8/9-12/13-16/17-20/21-24，Repetition times of repetition unit圖7 SSR類型統計結果圖Fig.7 SSR motif statistical results

A3SS，可變3′-端剪切；A5SS，可變 5′-端剪切；ES，外顯子跳躍；IR，內含子保留；MXE，外顯子互斥A3SS，Alternative 3′-splice sites；A5SS，Alternative 5′-splice sites；ES：Exons skipping；IR，Intron retention；MXE，Mutually exclusive exons圖8 不同剪接事件頻數統計圖Fig.8 Statistical diagram of frequency of different splicing class

3 討 論

2015年10月，PacBio公司推出全新升級的PacBio Sequel測序系統，該系統具有讀長長、通量高、準確率高等特點，無需打斷擴增，直接讀取反轉錄的全長cDNA，能夠有效獲取高質量的轉錄本全長序列，從而可準確辨別產生可變剪接的外顯子和轉錄起始位點的準確信息[24]。本研究利用PacBio Sequel系統對牦牛皺胃轉錄組全長進行測序分析，共獲得14 467 420條子序列，平均子序列長度為3 344.23 bp，結果顯示轉錄組所獲得的組裝序列完整性較好。對比牦牛參考基因組后，利用NR、Swiss-Prot、KEGG、KOG、eggNOG、GO、和Pfam七大公共數據庫進行功能注釋分析，在8 556條Unigenes中，有8 544條獲得NR數據庫注釋；有6 572條獲得KOG數據庫注釋，分為25個功能組分；有7 725條被注釋到GO數據庫，有1 721條注釋到42個代謝通路。

2018年，王明成[25]采用PacBio公司的RS Ⅱ系統對家牦牛的心臟、肝臟、脾臟、肺臟、腎臟、皮膚和睪丸組織進行全長轉錄組測序，得到約381 884條CCS序列，FLNC的轉錄本為251 712條，并將得到的全長轉錄組數據用于基因注釋分析。而本研究對牦牛皺胃組織進行全長轉錄組測序與分析，得到296 840條CCS序列，相較于王明成的結果，CCS序列少。由于單個組織樣本較多組織樣本混合測序得到序列更少，本研究得到277 402條FLNC并用于后續分析，對比顯示本研究得到更多的FLNC的轉錄本，說明新一代的Sequel系統更優于前代RS Ⅱ測序系統，序列完整度較高，產生較少的嵌合體序列。

本研究將所得到的Unigenes與KOG數據庫進行比對，根據其功能可分為25類，結果顯示較多基因富集到一般功能預測，翻譯后修飾、蛋白折疊和分子伴侶，信號轉導機制，細胞內運輸、分泌和囊泡轉運，說明牦牛皺胃的功能涉及多種蛋白或酶的合成以及多種大分子的轉運。在GO數據庫共有7 725條Unigenes被注釋，可劃分到3個大類共52個分支，GO功能富集顯示，三大類中分別是細胞過程、細胞、結合活性所富集基因最多，這與張春蘭[26]在羊上的研究結果相似。KEGG注釋到6個大類42個小類，其中涉及轉錄、翻譯、折疊、運輸、分解以及消化系統的基因居多，這可能跟牦牛皺胃所具有的消化功能有關，在不斷地合成各種消化酶，可以消化食物、分解蛋白質脂肪以及瘤胃和網胃進入皺胃中的微生物[27]。

本次共注釋得到的943個轉錄因子，分布在51個轉錄因子家族。其中轉錄因子較多的家族有zf-C2H2(157個)、bHLH(76個)和Homeobox(67個)等。zf-C2H2鋅指基因是人類最大的一類轉錄因子，也是哺乳動物最大的基因家族之一[28]；牦牛bHLH轉錄因子參與轉錄調節、DNA結合、DNA模板化、肌肉器官發育以及神經元分化等[29]，這可能與牦牛皺胃肌肉發育和皺胃細胞分化有關；Homeobox家族既在消化道發育過程中的內胚層又在中胚層成分中表達，而且許多持續表達至成年，本研究正是從成年牦牛皺胃注釋到很多Homeobox家族轉錄因子，Hox基因家族的異常表達還可能導致先天性消化道畸形[30]。

牦牛皺胃編碼區預測發現共有8 544個基因片段可視為蛋白編碼區，而共有82條注釋到溶菌酶基因，其他的注釋到DNA聚合酶、RNA聚合酶、ATP酶家族等，這些基因序列可以為后期研究提供參考。從轉錄組中共檢測到3 596個SSR位點，其中單核苷酸重復最多67.38%，占大多數，二核苷酸、三核苷酸為重復單元的序列也較多，分別是11.62%、20.22%，而以四、五、六核苷酸為重復單元的SSR序列則較少，分別只有0.42%、0.17%、0.19%，說明要開發利用SSR序列作為分子標記應重點考慮單核苷酸重復SSR序列，其次為復雜SSR序列：二核苷酸和三核苷酸重復SSR序列。獲得的SSR分析結果可作為后期牦牛SSR分子標記的基礎。此外，本研究還獲得牦牛皺胃轉錄組1 825個可變剪接事件，其中外顯子跳躍最多。鮑晶晶等[31]通過研究綿羊不同發育階段背最長肌組織中可變剪接發現，外顯子跳躍最多，與本研究相似，說明可變剪接的發生對皺胃肌肉組織發育過程中形成的功能各異的蛋白質起著重要作用。

4 結 論

本研究利用PacBio Sequel測序系統獲得了可靠的牦牛皺胃全長轉錄組數據，共獲得14 467 420條子序列，平均子序列長度為3 344.23 bp，質控得到CCS序列有296 840條，FLNC序列有277 402條。對轉錄組數據進行了NR、Swiss-Prot、KEGG、KOG、eggNOG、GO和Pfam七大數據庫對比注釋，分別注釋到8 544、8 475、1 721、6 572、8 491、7 725、8 162條Unigenes，通過轉錄因子注釋分析得到943個轉錄因子，分布在51個轉錄因子家族，編碼區預測到8 544個基因片段可視為蛋白編碼區，SSR分析檢索出3 596個SSR位點，得到1 825個可變剪接事件。本研究可為進一步研究牦牛皺胃生物學特性、相關代謝途徑、信號通路及其分子機制提供數據支持。