四川山鷓鴣的轉(zhuǎn)錄組組裝和注釋

劉乙， 何樂為， 蘭月， 周闖， 陳本平， 岳碧松, 3， 孟楊, 3*(.四川大學(xué)生命科學(xué)學(xué)院，生物資源與生態(tài)環(huán)境教育部重點(diǎn)實(shí)驗(yàn)室，成都60065; . 四川老君山國家級自然保護(hù)區(qū)管理局，四川屏山645350; 3. 四川大學(xué)生命科學(xué)學(xué)院，四川省瀕危動物保護(hù)生物學(xué)重點(diǎn)實(shí)驗(yàn)室，成都60065)

四川山鷓鴣隸屬雞形目Galliformes雉科Phasianidae，是中國特有的國家一級重點(diǎn)保護(hù)野生動物。由于四川山鷓鴣的分布區(qū)域狹窄，種群數(shù)量稀少，且長期受人類開發(fā)活動造成的棲息地喪失和片段化等因素的影響，其野生種群生存受到較為嚴(yán)重的威脅，被世界自然保護(hù)聯(lián)盟(IUCN)列為瀕危(EN)物種(BirdLife International，2016)。此外，IUCN與國際鳥盟、世界雉類協(xié)會將四川山鷓鴣納入鶉類保護(hù)組行動計劃(IUCN，2014)。四川山鷓鴣僅在中國的四川省和云南省被發(fā)現(xiàn)，分布于大相嶺山系的南緣、小相嶺山系的東緣、涼山山脈的東北部、烏蒙山系的西部。整個區(qū)域由10個棲息地斑塊和36個潛在棲息地斑塊組成，總面積約5 869 km(戴波等，2014)。

分子生態(tài)學(xué)方面，本課題組研究了四川山鷓鴣線粒體基因組、全基因組和該物種的系統(tǒng)發(fā)育(He.，2009；Zhou.，2019)，分析了基因組逆轉(zhuǎn)錄子和內(nèi)源性逆轉(zhuǎn)錄病毒(Cui.，2016；鄭帥等，2019)。系統(tǒng)發(fā)育關(guān)系研究方面，大部分基于線粒體基因組構(gòu)建的系統(tǒng)發(fā)育樹都認(rèn)同山鷓鴣屬是位于雉科的基部位置，是比較原始的類群，山鷓鴣屬的單系性得到了支持(Shen.，2009；Kan.，2010)。Chen等(2020)根據(jù)山鷓鴣屬14個物種的超保守元素、外顯子和線粒體基因組構(gòu)建的系統(tǒng)發(fā)育樹顯示，山鷓鴣屬是從非洲遷徙到東南亞定殖的。根據(jù)化石和地理模型估算，山鷓鴣屬的祖先在中新世早期到達(dá)印度支那，但是直到1 000萬年前全球冷卻加劇時才開始發(fā)散為 2個主要分支——“中國分支”和“東南亞分支”。四川山鷓鴣屬于“中國分支”的“橫斷山脈分支”，與環(huán)頸山鷓鴣的親緣關(guān)系最近，分化于330萬年前。

隨著測序技術(shù)的發(fā)展和成本的下降，單組學(xué)和多組學(xué)結(jié)合的分析將是今后一段時間內(nèi)的研究重點(diǎn)。本研究對1只成年雄性四川山鷓鴣個體的心臟、肝臟和腎臟進(jìn)行了轉(zhuǎn)錄組測序、組裝和注釋，為進(jìn)一步挖掘四川山鷓鴣的功能基因提供基礎(chǔ)數(shù)據(jù)，并為研究該瀕危物種提供遺傳數(shù)據(jù)。

1 材料與方法

1.1 組織采樣和測序

在四川老君山國家級自然保護(hù)區(qū)發(fā)現(xiàn)受傷的雄性四川山鷓鴣1只，因救助無效死亡。收集心臟、肝臟和腎臟組織，液氮研磨。將組織粉末加入TRIzol試劑(Invitrogen，Carlsbad，CA，USA)提取總RNA。對提取的RNA進(jìn)行2%瓊脂糖凝膠電泳以評估RNA的降解和污染。使用NanoPhotometer分光光度計(Implen，Los Angeles，USA)和帶有Qubit 2.0熒光計的Qubit RNA分析試劑盒(Life Technologies，Carlsbad，USA)檢查總RNA的純度和濃度。另外，使用Bioanalyzer 2100上的RNA Nano 6000分析試劑盒(Agilent Technologies，Santa Clara，USA)評估RNA完整性。使用Epicentre Ribo-zero去除rRNA試劑盒(Epicentre，Madison，USA)去除rRNA后，使用Illumina的NEBNext RNA文庫制備試劑盒(NEB，Ipswich，USA)構(gòu)建3個組織的cDNA文庫，并在北京諾禾致源公司的Illumina Novaseq 6000平臺上以150 bp的配對末端測序長度進(jìn)行測序。轉(zhuǎn)錄組數(shù)據(jù)以登錄號PRJNA638287保存在Gene Expression Omnibus(GEO)數(shù)據(jù)庫中(https://www.ncbi.nlm.nih.gov/bioobject/PRJNA638287)。

1.2 質(zhì)量控制和轉(zhuǎn)錄組組裝

使用NGS QC Toolkit v2.3.3(Patel & Jain，2012)和HISAT2 v2.1.0(Kim.，2015)對測序獲得的原始序列進(jìn)行數(shù)據(jù)質(zhì)量控制，過濾低質(zhì)量和帶接頭的數(shù)據(jù)。使用MultiQC(Ewels.，2016)評估數(shù)據(jù)質(zhì)量。使用Trinity v2.1.1(Grabherr.，2011)組裝得到的非冗余干凈序列。統(tǒng)計每個基因和轉(zhuǎn)錄本的原始讀取計數(shù)。并使用BUSCO v 5.1.2 (Sim?o.，2015)評估組裝質(zhì)量。

1.3 轉(zhuǎn)錄組注釋

轉(zhuǎn)錄組注釋分為結(jié)構(gòu)注釋和功能注釋2個部分。使用TransDecoder-v5.5.0(Grabherr.，2011)預(yù)測unigenes的蛋白編碼區(qū)，每個開放閱讀框(open reading frame，ORF)的長度最少為100個氨基酸。注釋功能基因時，通過BLAST搜索NCBI Non-Redundant Protein Sequences(NR)數(shù)據(jù)庫、Gene Ontology(GO)數(shù)據(jù)庫、KyotoEncyclopedia of Genes and Genome(KEGG)數(shù)據(jù)庫、Swiss數(shù)據(jù)庫和Cluster of Orthologous Group(COG)數(shù)據(jù)庫進(jìn)行比對，將閾值(E-value)設(shè)定為1×10。

2 結(jié)果

2.1 測序結(jié)果

心臟、肝臟和腎臟的原始序列過濾后分別產(chǎn)生了5.70 G、4.60 G和5.16 G數(shù)據(jù)(表1)。3個文庫的GC含量在48%左右。MultiQC序列質(zhì)量報告顯示，每條序列各位置堿基的測序質(zhì)量均位于綠色區(qū)間，具有平均質(zhì)量分?jǐn)?shù)的序列數(shù)量也是位于綠色區(qū)間，說明數(shù)據(jù)質(zhì)量很好。每條序列的重復(fù)水平均在20%以下的可接受水平。3個文庫中沒有檢測到過表達(dá)的序列(1條序列占比超過總數(shù)據(jù)的1%視為過表達(dá))，說明樣品未被污染。每條序列各位置的N堿基(無法識別的堿基記為N)比例可忽略不計，除了心臟文庫在97位堿基處有0.42%未識別的堿基。接頭序列的含量非常低，說明在對原始數(shù)據(jù)質(zhì)控時，接頭序列被徹底清除。

286 661條轉(zhuǎn)錄本經(jīng)過Trinity組裝并去掉冗余后共得到234 488個基因(表2)。使用BUSCO評估組裝的完整性為97.6%，在255個核心基因中檢測到249個完整的基因(122個單拷貝基因和 127個多拷貝基因)，另外6個基因包括4個片段化基因和2個未檢出基因。

2.2 基因結(jié)構(gòu)注釋

選取每個基因中最長的轉(zhuǎn)錄本作為代表序列(unigene)，對unigenes的結(jié)構(gòu)分析統(tǒng)計發(fā)現(xiàn)，81.29%不包含ORF，16.43%的包含1個ORF，而2.28%的包含2個及以上ORFs。

2.3 基因功能注釋

將組裝好的基因在主要的5個數(shù)據(jù)庫中進(jìn)行注釋，總共有70 737個基因獲得注釋結(jié)果，占30.17%，其中NR數(shù)據(jù)庫的注釋結(jié)果最多，COG數(shù)據(jù)庫的注釋結(jié)果最少(表3)。5個數(shù)據(jù)庫共同注釋到的基因6 998個，NR數(shù)據(jù)庫中單獨(dú)被注釋到的基因最多，為23 773個(圖1)。

表1 四川山鷓鴣RNA測序數(shù)據(jù)結(jié)果統(tǒng)計Table 1 Statistics of RNA sequencing data results of Arborophila rufipectus

表2 轉(zhuǎn)錄本和基因序列信息統(tǒng)計表Table 2 Statistics of transcripts and genes sequences

注： Contig N50： 將所有的Contigs按照從長到短進(jìn)行排序, 并將Contig按照這個順序依次相加, 當(dāng)相加的長度達(dá)到Contig總長度的一半時, 最后一個加上的Contig長度即為Contig N50

Notes: Contig N50: sort all Contigs from the longest to the shortest, and add Contigs in this order, when the added length reaches half of the total length of Contig, the last added Contig length is Contig N50

表3 五大數(shù)據(jù)庫注釋結(jié)果Table 3 Annotation results of five major databases

2.4 基因的GO注釋結(jié)果

GO數(shù)據(jù)庫將注釋得到的基因分為細(xì)胞成分(cellular component)、分子功能(molecular function)和生物過程(biological process)三大類，在細(xì)胞成分大類中，細(xì)胞(cell)和細(xì)胞組分(cell part)注釋到的基因最多；在分子功能大類中，結(jié)合(bingding)注釋到的基因最多；在生物過程大類中，細(xì)胞過程(cellular process)注釋到的基因最多(圖2)。

圖1 四川山鷓鴣轉(zhuǎn)錄組組裝的基因的數(shù)據(jù)庫注釋韋恩圖Fig. 1 Venn diagram of database annotation of Arborophila rufipectus transcriptome assembled genes

圖2 四川山鷓鴣轉(zhuǎn)錄組GO數(shù)據(jù)庫注釋圖Fig. 2 GO database annotation diagram of Arborophila rufipectus transcriptome

3 討論

本研究采用讀長短但準(zhǔn)確性高的二代Illumina測序。數(shù)據(jù)質(zhì)量控制檢測報告顯示，每條序列各位置堿基的測序質(zhì)量均位于代表數(shù)據(jù)質(zhì)量很好的綠色區(qū)間，每條序列各位置中無法識別的N堿基比例可以忽略不計，這體現(xiàn)了二代測序的準(zhǔn)確性。但二代測序的讀長短(50～500 bp)，如基于150 bp的配對末端測序長度進(jìn)行測序組裝出來的基因的Contig N50只有964 bp，說明組裝出來的基因碎片化，從大量的短序列中準(zhǔn)確組裝出完整的轉(zhuǎn)錄組仍然充滿挑戰(zhàn)。三代測序解決了讀長問題，其測序長度高達(dá)上百kbp，但錯誤率高(1%～15%)，需要二代測序輔助降低錯誤率。

從頭組裝的方法是基于序列之間的重疊部分完成的轉(zhuǎn)錄組組裝，使用從頭組裝的方法進(jìn)行組裝的應(yīng)用價值是對四川山鷓鴣原有基因組的補(bǔ)充，因其已有參考基因組也是通過二代Illumina測序組裝的(Zhou.，2019)，其組裝過程中存在錯誤和組裝不完全的情況。相對于基于DNA的基因組組裝來說，基于RNA的可以將那些由多個轉(zhuǎn)錄本加工而成的成熟mRNA組裝出來。在研究目標(biāo)是轉(zhuǎn)錄本序列，且物種無參考基因組或者參考基因組質(zhì)量不高的轉(zhuǎn)錄組組裝中，Trinity是所有從頭組裝軟件中準(zhǔn)確度最高且組裝更完整的軟件，能在不損失準(zhǔn)確性和運(yùn)行速度的前提下得到最完整的組裝版本(盧戌，2013)。

基因的注釋是整個流程中最重要的一個環(huán)節(jié)，尤其是在高通量測序成本日益下降的情況下。組裝拼接出來的草圖需要進(jìn)行注釋后才能體現(xiàn)價值，后續(xù)的功能基因挖掘也建立在注釋結(jié)果上。基因功能的注釋是建立在結(jié)構(gòu)注釋的基礎(chǔ)上的，根據(jù)結(jié)構(gòu)注釋的結(jié)果提取具有翻譯功能的區(qū)域，并與主要的數(shù)據(jù)庫比對。結(jié)構(gòu)注釋中可以看到只有18.71%的基因具有ORF框，即這部分序列有編碼翻譯成蛋白的潛能。這種現(xiàn)象普遍存在于真核生物的基因組中，大部分的基因被認(rèn)為是“垃圾”序列，但具有重要的轉(zhuǎn)錄后調(diào)控作用，涉及轉(zhuǎn)錄起始、轉(zhuǎn)錄中和轉(zhuǎn)錄后的每一個過程，如lncRNA、MicroRNA、piRNA、siRNA等(Kapranov.，2007)。功能注釋中NR數(shù)據(jù)庫的結(jié)果最豐富。

四川山鷓鴣的近緣物種海南山鷓鴣的棲息地主要是海南島的熱帶雨林和山地常綠林。海南山鷓鴣在海南霸王嶺國家級自然保護(hù)區(qū)海拔350～1 560 m(韋鋒等，2008)、海拔80～260 m的南味嶺林區(qū)的次生林均有分布(邱勝榮，丁長春，2014)。四川山鷓鴣在四川麻咪澤自然保護(hù)區(qū)內(nèi)的生境(2 100～2 300 m)(趙成等，2015)與海南山鷓鴣有明顯的海拔梯度差異，后續(xù)可利用本研究數(shù)據(jù)與海南山鷓鴣進(jìn)行比較轉(zhuǎn)錄組學(xué)研究，探討近緣物種對不同海拔生境適應(yīng)的遺傳學(xué)機(jī)制。