文獻(xiàn)解讀|利用Pacbio Iso-Seq測序技術(shù)發(fā)現(xiàn)水稻非生物脅迫下的新轉(zhuǎn)錄本和新基因

文章題目：Utilizing PacBio Iso-Seq for Novel Transcript and Gene Discovery of Abiotic Stress Responses in Oryza sativa L.

中文題目：利用Pacbio Iso-Seq測序技術(shù)發(fā)現(xiàn)水稻非生物脅迫下的新轉(zhuǎn)錄本和新基因

發(fā)表期刊：International Journal of Molecular Sciences

發(fā)表時間：2020年10月31日

影響因子：4.556

研究背景

全球氣候變化導(dǎo)致高溫、干旱和夜間高溫等非生物脅迫條件的嚴(yán)重程度和頻率增加，這些都造成了作物產(chǎn)量的降低。隨著世界人口的增長，植物育種專家面臨著開發(fā)高產(chǎn)、優(yōu)質(zhì)、減少環(huán)境污染的新品種的艱巨任務(wù)。水稻是世界上一半以上人口的主要卡路里來源，特別是對亞洲最貧窮的人來說。世界各地的基因庫中保存著23萬多份水稻及其野生近緣種的廣泛自然遺傳多樣性的種質(zhì)資源，是一種無價的可用于作物改良的重要基因庫。

雖然世界上近80%的水稻種植是基于indica（秈稻亞種）品種，但目前的標(biāo)準(zhǔn)基因組及其注釋來自粳稻亞種Nipponbare。由于缺乏適當(dāng)?shù)幕蚪M，不同水稻亞種的研究大多都基于Nipponbare基因組。例如，在3000水稻基因組計(jì)劃中將測序序列比對到Nipponbare基因組上，丟棄了不能比對到該參考基因組的序列。這可能會導(dǎo)致非粳稻亞種特有的遺傳信息的丟失。另外，最近已經(jīng)對其它水稻亞種的栽培品種的基因組進(jìn)行了測序，例如indica（Shuhui498，Zhenshan 97, Minghui 63）、aus（Kasalath，N22），但其完整性和注釋程度仍存在差異。值得注意的是aus亞種是抗病、耐磷酸鹽缺失、耐澇、耐厭氧發(fā)育和抗旱等潛在性狀的寶貴基因來源。例如，在aus品種基因組中發(fā)現(xiàn)了耐磷酸鹽缺失相關(guān)基因OsPSTOL1、耐澇相關(guān)基因OsSNORKEL1/2和OsSUB1A。值得注意的是，這些基因在粳稻的Nipponbare亞種基因組序列中是不存在的。

在過去的幾年里，RNA測序（特別是基于illumina的短序列RNA-seq）已經(jīng)成為分析轉(zhuǎn)錄組的有力工具，用來識別在非脅迫控制和各種環(huán)境脅迫條件下差異表達(dá)的基因。然而，需要基于參考基因組或轉(zhuǎn)錄組序列對RNA-seq數(shù)據(jù)進(jìn)行比對和注釋來確定轉(zhuǎn)錄水平。在水稻中，參考基因組決定了可以鑒別的差異表達(dá)基因和轉(zhuǎn)錄本亞型。顯然，參考基因組/轉(zhuǎn)錄組中沒有的基因的表達(dá)信息在分析過程中會丟失。這在研究耐脅迫的外來品種、陸地品種或野生稻種時尤其相關(guān)，因?yàn)樗鼈兛赡芎袇⒖计贩NNipponbare不存在的耐受基因。這將嚴(yán)重限制識別支持作物改良計(jì)劃的新候選基因的可能性。
解決這一問題的一個顯而易見的辦法是對所需的基因組進(jìn)行測序、組裝和注釋。但是，這種方法比較昂貴和耗時。在這篇文章里，我們探索了一種更有針對性的RNA-Seq序列方法來測序和重建了三個不同亞種的水稻品種的部分轉(zhuǎn)錄本作為參考，Pacific Bioscience（PacBio）屬于提供高通量全長轉(zhuǎn)錄本序列的新一代測序方法。該方法已成功應(yīng)用于對現(xiàn)有植物轉(zhuǎn)錄本及注釋的探索和擴(kuò)展，如高粱、小麥、甘蔗、野生棉花、不同的穗型草、苜蓿等。

取樣材料

樣品取自三個水稻品種的10個不同亞種的不同組織部位（表1）：

分析結(jié)果

1.重構(gòu)轉(zhuǎn)錄本

使用PacBio Sequel I測序平臺對每個品種進(jìn)行SMRT測序，得到15.49~24.51GB的轉(zhuǎn)錄本數(shù)據(jù)。用IsoSeq3軟件對原始測序數(shù)據(jù)進(jìn)行ccs和lima處理，每個品種SMRT cell分別得到460340~736747條全長非嵌合序列（full-length non-chimeric reads簡稱FLNC，包含3 ‘ 引物、5 ‘ 引物以及polyA尾）。全長非嵌合經(jīng)過IsoSeq3聚類和polish分別得到37951~54684高質(zhì)量轉(zhuǎn)錄本（HQ）以及1233 ~2170低質(zhì)量轉(zhuǎn)錄本（LQ）。先將HQ與NCBI核苷酸數(shù)據(jù)庫進(jìn)行blastn比對（E<=1e-10），再將上一步為比對上的轉(zhuǎn)錄本序列與NCBI蛋白數(shù)據(jù)庫進(jìn)行blastx比對（E<=1e-10），去除未比對上兩個數(shù)據(jù)庫的轉(zhuǎn)錄本序列，最終得到37535~54594條HQ用于后續(xù)分析（表2）。

Pacbio RSII平臺聲明使用RNA-seq二代測序數(shù)據(jù)對轉(zhuǎn)錄本數(shù)據(jù)進(jìn)行矯正，可以得到更多的HQ序列，因?yàn)長Q序列中含有大量的插入和缺失。然而與RSII相比，PacBio Sequel I測序平臺的測序結(jié)果更好。為了驗(yàn)證這一結(jié)果，我們用minimap2將未矯正的HQ比對到相應(yīng)亞種的基因組，結(jié)果表明缺失的比例很小，所以進(jìn)一步的分析中沒有包含LQ序列。

2.轉(zhuǎn)錄本去冗余

在文庫準(zhǔn)備過程中，會產(chǎn)生5 ‘ RNA降解產(chǎn)物，并進(jìn)行測序。這些降解產(chǎn)物具有相同的外顯子結(jié)構(gòu)，但缺乏5 ‘序列信息，因此產(chǎn)生與技術(shù)偏差或生物學(xué)背景無關(guān)的冗余異構(gòu)體。我們測試了三種不同的去冗余方法，包括cogent、cDNA cupcake和TAMA，其中cDNA cupcake和TAMA需要基于參考基因組，而cogent不需要基于參考基因組。cogent基于pacbio全長轉(zhuǎn)錄本序列重構(gòu)一個參考基因組，然后將相同的序列比對到重建的基因組，基于比對結(jié)果利用cDNA cupcake算法對轉(zhuǎn)錄本去冗余。cDNA cupcake和TAMA直接用minimap2軟件和各自的亞種參考基因組進(jìn)行比對?；谶@三種方法，只有很少的轉(zhuǎn)錄本不能回比到基因組上（表3）?？偟膩碚f，這三種去冗余方法均能顯著減少異構(gòu)體的數(shù)量，分別為47.6% （cDNA cupcake，Nipponbare）和68.3%（cogent，Dular）。

基于植物中430個高度保守的同源蛋白利用BUSCO軟件對TAMA算法去冗余前后的HQ進(jìn)行完整性評估（圖一），由于取樣不完全，缺失了54%~27%的重要蛋白，其中Nipponbare參考基因組（IRGSP）只缺失了6種。cDNA cupcake和TAMA的結(jié)果相似，而對于cogent，超過50%的蛋白缺失，最有可能的原因是轉(zhuǎn)錄本沒有回比到重建的基因組。

去冗余后轉(zhuǎn)錄本長度中值都有所增長，長度分布和轉(zhuǎn)錄本長度中值與Nipponbare參考基因組相似。統(tǒng)計(jì)了去冗余后10個品種每個基因相應(yīng)的轉(zhuǎn)錄本數(shù)量，其中基因只有一個轉(zhuǎn)錄本的比例，TAMA最高達(dá)到了75%，cDNA cupcake在60%左右，cogent只有50%。同時計(jì)算了Nipponbare參考基因組每個基因?qū)?yīng)的轉(zhuǎn)錄本數(shù)量進(jìn)行比較，該參考基因中基因只有一個轉(zhuǎn)錄本的比例達(dá)到了85%（圖2）。

來自同一亞種的不同品種親緣關(guān)系更近，我們使用系統(tǒng)發(fā)育樹評估亞種之間的遺傳距離。利用去冗余后的轉(zhuǎn)錄本序列基于Nipponbare參考基因組識別SNPs，使用SNPhylo繪制進(jìn)化樹（圖3）。SNPhylo提取高質(zhì)量并且具有代表性的SNPs進(jìn)行后續(xù)分析，cDNA cupcake算法大約30000個SNPs，cogent算法大約23200個SNPs，TAMA算法大約16000個SNPs。三種方法中，同一亞種的不同品種聚類在了一起，cDNA cupcake算法和cogent的聚類結(jié)果更相似。三種方法均能將aus與另外兩個亞種區(qū)分開，但cogent和TAMA對indica和japonica種間的區(qū)分不如cDNA cupcake明顯。

3.評估重構(gòu)的轉(zhuǎn)錄本

基于TAMA算法得到的HQ進(jìn)行轉(zhuǎn)錄本的評估。由于TAMA只對比對到參考基因組上的轉(zhuǎn)錄本進(jìn)行去冗余，我們用cogent對沒比對上參考基因組的轉(zhuǎn)錄本進(jìn)行去冗余。合并結(jié)果后，每個品種最終得到10511（Dular）~15011（IR64）個基因，14255（Dular）~20803（Moroberekan）個轉(zhuǎn)錄本（表4）。與Nipponbare參考基因組相比，大約三分之一的基因位點(diǎn)和大約一半的轉(zhuǎn)錄模型被重建。每個品種每個基因的平均轉(zhuǎn)錄數(shù)約為1.4～1.5，略高于參考基因組的1.2。中位轉(zhuǎn)錄本長度為986 bp（Dular）~1394 bp（Nipponbare），與Nipponbare參考值1385 bp相似。平均GC含量在50.87%（Dular）~52.76%（IR64），與Nipponbare參考值51.24%相似。利用gffcompare軟件與Nipponbare參考基因組進(jìn)行比較識別新基因與轉(zhuǎn)錄本。

4.功能注釋

為了深入了解重建轉(zhuǎn)錄本的生物學(xué)信息，我們進(jìn)行了功能注釋。使用TransDecoder軟件預(yù)測開放閱讀框（ORFs），包括blast和PFAM，結(jié)果表明大約有60%~70%的完整ORFs（包括啟動和終止密碼子）。此外還發(fā)現(xiàn)了26%~38%的5 ‘ ORF、很少比例的3 ‘ ORF和中間ORF（既沒有起始密碼子也沒有終止密碼子）（圖4）。

使用Trinotate和Mercator4進(jìn)行功能注釋。Mercator4是專門為植物開發(fā)的，它使用了一種簡單的層次樹結(jié)構(gòu)，被稱為“容器”，用來描述生物學(xué)概念。主要的生物過程如光合作用，都是由頂層的容器來表示的，每個子容器描述的是一個更詳細(xì)的子過程。目前本體包括27個功能類別，代表了植物中不同的生物過程。N22、IR64和Nipponbare三個品種作為各自亞種的代表與植物中所有水稻基因的分類進(jìn)行比較分析，結(jié)果顯示三個品種的注釋結(jié)果分布相似（圖5）。超過28000個水稻已知基因在Mercator庫中沒有注釋分類信息，因此三個品種有8000~10000個轉(zhuǎn)錄本沒有分類注釋到Mercator庫。

5.品種間共有和特有的轉(zhuǎn)錄本

為了鑒定品種特異性轉(zhuǎn)錄本，以N22、IR64和Nipponbare三個品種的轉(zhuǎn)錄本作為blast比對庫，其它9個品種與其進(jìn)行比對（圖6）。識別到N22特有轉(zhuǎn)錄本652個，IR64特有轉(zhuǎn)錄本2426個，Nipponbare特有轉(zhuǎn)錄本349個。

6.aus特有轉(zhuǎn)錄本的差異表達(dá)分析

aus品種N22特別抗旱和耐熱脅迫，因此我們想知道在這些條件下是否有aus特異轉(zhuǎn)錄本受到調(diào)控。以N22為研究對象，分析干旱和熱脅迫下差異表達(dá)的基因。利用從發(fā)育種子中分離的RNA進(jìn)行Illumina測序，將測序數(shù)據(jù)回比到重構(gòu)的N22轉(zhuǎn)錄本。使用DESeq2基于參數(shù)FDR<0.1和|log2FC|>=1軟件識別出56個aus特異的差異表達(dá)基因。這56個差異基因進(jìn)行blast比對，其中46%比對上擬南芥，27%沒有任何注釋信息，11%僅描述了一個PFAM域或與其它植物物種的序列同源，而在水稻中僅有5%已知同源基因。

舉個例子，在高溫和干旱雙重脅迫下顯著上調(diào)的基因B12288。它在japonica和indica中均有同源基因RAB21，這個基因受干旱的誘導(dǎo)，其編碼的蛋白屬于LEA脫氫蛋白家族。與水稻其它脫氫蛋白進(jìn)行多序列比對研究（圖7），N22實(shí)際上的基因與野生稻、O. sativa ssp. japonica中其它4種脫氫酶的親緣關(guān)系密切。序列覆蓋率為89.5%，序列同源性86.0%，其中包含脫氫酶高度保守的重復(fù)區(qū)。相比japonica蛋白，N22蛋白序列與野生稻更接近。

總結(jié)
本文主要探討了Pacbio Iso-Seq獲得的轉(zhuǎn)錄本相比于Nipponbare參考基因組是否可以用于aus等水稻亞種的下游分析。此外通過這些轉(zhuǎn)錄本，我們希望發(fā)現(xiàn)水稻非生物脅迫下新的轉(zhuǎn)錄本和基因。我們的分析表明所有品種都可以鑒定出特異的轉(zhuǎn)錄本，還確定了aus亞種特異的差異表達(dá)基因。Pacbio Iso-Seq這種方法適用于其它沒有基因組或者基因組質(zhì)量不高的物種，相比對基因組組裝，這種方法更省時便宜。