單葉省藤和黃藤基因組組裝研究

王思寧 趙韓生 高志民

(國際竹藤中心 北京 100102 )

棕櫚藤是世界上最重要的非木材森林植物資源之一，屬于棕櫚科(Arecaceae)，是自然進化中出現(xiàn)的攀爬植物的代表。近年來的研究表明，在棕櫚科省藤族省藤亞科中棕櫚藤植物可分為11屬，共有631種，其中省藤屬(NCBITaxonID:4711)和黃藤屬(NCBITaxonID:93268)種類最多，分別占藤種的約65%和20%。這2個屬也是最重要的棕櫚藤藤材來源，為工業(yè)生產提供了95%以上的藤條。全球有超過500萬人在經濟上依賴棕櫚藤，藤工業(yè)產值每年約70億美元，包括國內工業(yè)生產、國際貿易，以及藤編織的籃子、座椅和家具等。同時，對棕櫚藤遺傳育種的關注在持續(xù)增長，并且預計將來在幾年內棕櫚藤的種植面積將逐漸超過天然棕櫚藤面積。

單葉省藤(Calamussimplicifolius)(NCBITaxonID:746888)是中國原產、種植面積最廣的藤種之一，莖長約50m，直徑約15mm。單葉省藤是海南島特有的棕櫚藤，可生產高品質、中等直徑的藤條，用于捆扎和編織。黃藤(Daemonoropsjenkinsiana)(NCBITaxonID:1510057)，是常綠攀援棕櫚藤的代表，自然生長在海拔1000m以下的低地雨林，分布于孟加拉國、不丹、柬埔寨、印度、老撾、緬甸、尼泊爾、泰國、越南和中國的東南部。黃藤叢生生長，莖可達50m, 節(jié)間長40cm，直徑30mm。這2種產值最高的藤種現(xiàn)種植面積超過1000hm2，栽培在中國小于北緯23°30′的地區(qū)，例如海南島、廣東、廣西、云南、福建等華南地區(qū)。

單葉省藤和黃藤由于其藤條韌性強和耐久度好而具有廣泛的用途和巨大的開發(fā)潛力。需要應用分子育種技術開展棕櫚藤育種，來滿足對棕櫚藤質量和數(shù)量不斷增長的需求。然而，棕櫚藤遺傳背景不清，嚴重阻礙了對其分子生物學的研究和實際生產的應用，更影響了相關物種之間比較基因組分析的深入開展。因此，我們使用最新的高通量測序技術，完成了單葉省藤和黃藤的測序和染色體水平基因組組裝。隨著這2個棕櫚藤基因組的完成，使許多比較基因組分析和其他下游分析成為可能，包括分子標記的開發(fā)、功能基因的鑒定和分子輔助育種等。此外，高質量的基因組將有助于對棕櫚藤生物性狀進行基因組學、轉錄組學和代謝組學分析。本文選取并鑒定了棕櫚藤中改良藤條性能的一些重要候選基因——木質素生物合成基因家族，為今后利用木質素生物合成基因來改良藤條的品質以及研究物種多樣性奠定了基礎。

1 材料方法與數(shù)據(jù)處理

1.1 DNA分離、文庫構建和高通量測序

單葉省藤和黃藤材料取自中國林業(yè)科學研究院熱帶林業(yè)研究所(北緯23°11′29″，東經113°22′40″，海拔87m)，于2015年春季采集處于營養(yǎng)生長階段的幼嫩葉片組織用于總DNA提取。使用DNeasyPlantMiniKits(Qiagen)分離和提取總DNA。根據(jù)高分子量核DNA分離法純化基因組DNA獲得多個DNA文庫，并在IlluminaHiSeq4000和PacBioSequel平臺上進行測序(原文表1)。即：在基于Illumina的方法上構建了3個不同插入片段大小的文庫(270bp、500bp和800bp)，用于PE測序，以及4個不同插入片段大小的文庫(2kb、5kb、10kb和20kb)用于MP測序。另外，按照PacBio方法，構建了5個插入片段大小為20kb的PacBioSequel文庫，進行測序。通過對測序數(shù)據(jù)清理和預處理，最終得到約494.08Gb的單葉省藤高質量數(shù)據(jù)(322.3GbPEreads,93.4GbMPreads,78.38GbPacBiodata)，覆蓋率為單葉省藤基因組的252倍， 約426.17Gb的黃藤高質量數(shù)據(jù)(244.58GbPEreads,103.21GbMPreads,and78.38GbPacBiodata)，覆蓋率為黃藤基因組的266倍。

作為與Illumina和PacBio文庫構建測序，進行了另一平行分析，即使用相同組織材料為單葉省藤和黃藤建立了2個高通量測序文庫。使用MboI限制性內切酶在甲醛固定后消化DNA，然后用生物素標記5個殘基。在原位連接平末端片段后，對分離出的DNA進行反向交聯(lián)、純化和過濾，得到包含生物素標記的片段。隨后，按此順序進行DNA片段端修復、適配體連接和聚合酶鏈式反應。然后，根據(jù)BGISEQ-500的使用說明，對其進行100PE的測序。在分別獲得了大約148Gb和154Gb的原始數(shù)據(jù)后，經過低質量數(shù)據(jù)的過濾和其他預處理，最終獲得約6.7Gb和13.1Gb的高質量數(shù)據(jù)，并使用HiC-Pro(version2.8.0devel)軟件分別對單葉省藤和黃藤進行了評價和分析(原文表1)。

1.2 基因組調查

對給定物種基因組大小、雜合性等基因組特征的了解解，將有助于開發(fā)一套適合的測序方案和組裝策略。因此，本研究采用4種獨立的方法來估算基因組大小，即KmerSpectrumPlot.plinALLPATHS-LG(版本r52488)、GCE(GenomeCharacteristicsEstimation,released7Jan.2015)、JELLYFISH(版本2.0)和流式細胞術(原文表S1、S2，圖S1、S2)。在調研中，我們測序獲得單葉省藤約為98Gb，黃藤約為60Gb的二代原始測序數(shù)據(jù)作為調研數(shù)據(jù)。經過數(shù)據(jù)預處理和組裝，結果顯示(原文表S1)，最終預測單葉省藤和黃藤基因組大小分別約為1.98Gb和1.61Gb，相關雜合度分別為1.32%～1.52%和1.19%～1.31%。因此，基因組分析結果表明，這2種棕櫚藤的基因組比較適用于Illumina和PacBio測序數(shù)據(jù)的聯(lián)合組裝。

1.3 利用Illumina和PacBio高通量測序數(shù)據(jù)聯(lián)合組裝基因組

在Illumina數(shù)據(jù)的預處理過程中，過濾掉了低質量的讀長和接頭序列，獲得單葉省藤和黃藤的有效數(shù)據(jù)分別為416Gb和348Gb。對于PacBio數(shù)據(jù)，使用了MECAT(發(fā)布于2017年6月27日)軟件，得到單葉省藤和黃藤分別為52Gb和32Gb的高質量數(shù)據(jù)。隨后，使用FALCON(版本0.3)軟件進行第1次contig組裝。如原文表S3所示，使用不同參數(shù)為單葉省藤基因組生成了2個組裝結果，一個ContigN50為67.2kb，基因組約為1.59Gb(約為預測基因組大小的80%)，另一個ContigN50為66.7kb，基因組約為1.53Gb(約為預測基因組大小的77%)。此外，獲得一個黃藤組裝結果ContigN50為81.5kb，基因組約為1.27Gb(約為預測基因組大小的79%)。通過以上結果顯示，僅用三代數(shù)據(jù)和MECAT軟件不能獲得對2種高質量棕櫚藤基因組，可能是由于2種棕櫚藤的高雜度(單葉省藤1.32%～1.52%，黃藤1.19%～1.31%)、高比例重復序列(單葉省藤為54.15%，黃藤為70%；有關詳細信息請參閱原文中的后續(xù)分析)和使用了測序深度不足的三代數(shù)據(jù)(PacBio數(shù)據(jù)經錯誤校正后分別為26×和20×)造成的。因此結合以上結果，我們嘗試利用Illumina和PacBio測序數(shù)據(jù)，對單葉省藤和黃藤進行了聯(lián)合基因組組裝。首先，利用Platanus(版本1.2.4，是一種用于高雜合數(shù)據(jù)的全新基因組匯編器，可以通過構建具有自動優(yōu)化的k-mer大小的DeBruijn圖)將PE片段整合到contig序列中。其次，用DBG2OLC(2015年7月11日發(fā)布，參照以下參數(shù)：DBG2OLCContigscontig.faLD0K17KmerCovTh4MinOverlap25AdaptiveTh0.007RemoveChimera1fscaffold.fa.) 矯正PacBio測序數(shù)據(jù)并組裝scaffolds，以分別獲得單葉省藤和黃藤約1.92Gb和1.56Gb的初始組裝結果。再次，在進行SSPACE分析前，參照DBG2OLC(原文表S4)進行了平行分析，這一步驟有助于提高基因組組裝的質量，減少SSPACE過程中的錯誤。然后，使用SSPACE(版本3.0)和MP大片度測序數(shù)據(jù)延長scaffolds，并通過使用GapCloser(版本1.12)和PBJelly(2015年8月24日發(fā)布) 進行補洞。最終，獲得了1.96Gb的單葉省藤組裝數(shù)據(jù)，包含5116個Contigs，ContigN50長度為107kb，ScaffoldN50長度為803kb，同時獲得了一個1.60Gb的黃藤組裝數(shù)據(jù)，包含ContigN50和ScaffoldN50長度分別為108kb和784kb(原文表2)。

隨后，通過3D-DNA分析流程(版本170123)對有效的高通量測序數(shù)據(jù)與上述裝配數(shù)據(jù)一起進行處理，得到了一個合理的連接模式，準確的染色體組裝結果。如原文圖2所示，使用Juicebox(版本1.5.2)顯示連接模式圖譜。原文表S5展示了單葉省藤12個最長的染色體水平Scaffold和黃藤13個最長的染色體水平Scaffold。預測染色體總長度分別占單葉省藤和黃藤基因組總長度的92.08%和92.01%，ScaffoldN50分別為169Mb和119Mb。

1.4 基因組評估

采用3種獨立的方法評估了單葉省藤和黃藤基因組組裝的準確性和完整性。首先，對基因組中不確定堿基的百分比(Ns)和鳥嘌呤、胞嘧啶含量(GC)進行了分析。結果顯示在每個基因組中低豐度的Ns(單葉省藤約為0.6%和黃藤約為0.7%)和總GC含量(單葉省藤為41.78%和黃藤為41.07%)與相關轉錄組數(shù)據(jù)相似(單葉省藤為41.68%和黃藤為41.89%)。此外，利用BasicLocalAlignmentSearchTool(BLAST)-like比對工具(版本1.0)(默認參數(shù))進行RNA測序(RNAseq)數(shù)據(jù)所組裝出來的功能基因序列與組裝得到的基因組序列比對，結果顯示，在一個大片段上90%以上的序列可以與基因組的組裝序列匹配(單葉省藤為92.89%和黃藤為81.81%)(原文表S6)。最后，利用BUSCO(版本3.0)對2種棕櫚藤基因組的完整性進行了評估，這是一種從進化的角度定量對單拷貝直系同源基因含量進行計算，以評估基因組完整性的方法。BUSCO結果顯示，在單葉省藤基因組組裝結果中檢測到保守性為96.4%的BUSCO蛋白(基于embryophytaodb9庫文件)(包括3.8%的BUSCO蛋白片段)。同時，在黃藤基因組組裝結果中檢測到保守性為87.3%的BUSCO蛋白(包括4.0%的BUSCO蛋白片段)(原文表S7)。

1.5 重復注釋

本研究在蛋白質編碼基因模型預測之前，完成了單葉省藤和黃藤基因組中轉座因子(TEs)和串聯(lián)重復序列的鑒定。采用2種獨立的方法來預測序列重復：同源注釋和從頭預測(denovo)法。同源注釋結果顯示，利用RepeatMasker(版本4.0.5)和RepeatProteinMasker(版本4.0.5)對Repbase文庫(發(fā)布于2017年12月01日)進行搜索，查找轉座因子序列。從頭預測結果顯示，利用Repeat-Modeler(版本1.0.8)和LTRFINDER方法在去除污染和多拷貝基因后，可以構建一個重復序列文庫。然后，使用RepeatMasker對基因組序列進行了分類，并與重復序列文庫進行了對比。此外，串聯(lián)重復序列由TandemRepeatFinder(版本4.09)識別，其參數(shù)為：Match=2,Mismatch=7,Delta=7,PM=80,PI=10,Minscore=50andMaxPeriod=2000。最終結果表明，長末端重復序列(longterminalrepeat,LTR)是最豐富的重復序列類型，2種非LTR逆轉座子(短穿插核元件和長穿插核元件)在2種棕櫚藤基因組中占比例最少(原文表S8)。TEs在單葉省藤和黃藤基因組中分別占54.15%和70%。TEs序列分歧度結果表明，從頭預測的重復序列比Repbase預測的重復序列更多(原文圖3)。

1.6 RNA樣本收集、文庫構建和轉錄組組裝

分別選取單葉省藤和黃藤3個不同發(fā)育階段的纖鞭樣本，每個樣本3個生物學重復(原文表S9)。為保證與2種棕櫚藤基因組研究的一致性，RNA取樣的位置和DNA取樣的位置一致。利用RNA提取試劑盒(TaKaRa, 日本)提取各樣本的總RNA，并使用NanoDrop2000分光光度計測定RNA的質量和濃度。使用RecombinantDNaseI于37℃孵育30min去除殘余DNA，并使用反轉錄試劑盒(Promega, 美國)合成cDNA第一條鏈，操作參照相關試劑盒說明書進行。然后，使用BGISEQ-500平臺對文庫池中的短100PE片段測序。在對轉錄組數(shù)據(jù)進行預處理時，使用SOAPnuke(版本1.5.6)對匹配序列和低質量的讀長序列進行過濾，過濾參數(shù)為：-n0.001-l20-q0.4-q2。所有樣本均使用Trinity(版本2.0.6)按以下參數(shù)組裝：(1)grouppairsdistance500, (2)mincontiglength200, (3)minkmercov2, (4)minglue2, (5)bflyopts-V5, (6)edge-thr=0.1, (7)stderr,and(8)SSlibtypeRF。然后，對Trinity的輸出結果進行聚類分析，使用TGI聚類工具(版本v2.0.6)生成一組非冗余集合，參數(shù)如下：(1)aminimumof95%identitybetweenthecontigs, (2)aminimumof35overlappingbases, (3)aminimumscoreof35,and(4)amaximumof20unmatchedoverhangingbasesatthesequenceends。最后，根據(jù)序列相似度將整合的轉錄本分為2類：成簇轉錄本和單一轉錄本。在成簇轉錄本中，轉錄本之間的序列相似區(qū)都超過了70%，由這些轉錄本被剪接成同源基因或旁系同源基因。此外，所有鑒定的基因將被用于后續(xù)的功能分析。

1.7 基因建模與預測

使用3種獨立的方法對完整的蛋白質編碼基因模型進行了綜合預測，包括從頭預測、同源蛋白預測和轉錄組數(shù)據(jù)預測。使用AUGUSTUS(版本3.3)的默認參數(shù)對重復區(qū)段進行第一次注釋，該程序是一個基于自我調試模型的從頭預測軟件。通過進行多次調試優(yōu)化數(shù)據(jù)，在單葉省藤和黃藤中分別預測獲得了85246和87613個基因模型。在基于同源性的預測中，以油棕(Elaeisguineensis)、長葉刺葵(Phoenixdactylifera)、二穗短柄草(Brachypodiumdistachyon)、水稻(Oryzasativa)、粟(Setariaitalica)、高粱(Sorghumbicolor)和玉米(Zeamays)等7個物種作為參考基因組。從ENSEMBL數(shù)據(jù)庫下載它們的蛋白序列并應用TBLASTN(版本2.2.26)對單葉省藤和黃藤中E-value不小于1e-5的組裝序列進行比對。然后，由GeneWise(版本2.0)生成拼接結果。轉錄組數(shù)據(jù)預測結果顯示，使用HISAT2(版本2.0.2)來識別外顯子-內含子剪接位置，并進一步進行轉錄組序列與基因組序列的比對。使用Cufflinks(版本2.2.1)分別在單葉省藤和黃藤中分別鑒定了56024和58134個蛋白質編碼基因模型(原文表S10)。最后，使用MAKER(版本2)整合了上述3個獨立預測的數(shù)據(jù)結果，隨著相同基因序列在單葉省藤和黃藤基因組中的出現(xiàn)，最終共分別鑒定出51235和53342個完整的蛋白編碼基因模型。

1.8 基因注釋和基因功能預測

使用2種獨立的方法對預測的基因集合進行評估：基因功能評估和BUSCO完整性評估。在基因功能評估中，比照基因在近緣物種中同源蛋白序列和預測的蛋白集合，評估了預測的注釋和蛋白質序列的一致性。參照NCBI非冗余蛋白質數(shù)據(jù)庫(2018年3月13日發(fā)布)、SWISS-PROT(2018年1月1日發(fā)布)、GeneOntology(GO)(2013年10月30日發(fā)布)、KyotoEncyclopediaofGenesandGenomes(KEGG)(數(shù)據(jù)集v81)和InterPro(數(shù)據(jù)集v.53)等5個權威蛋白數(shù)據(jù)庫進行注釋比對，結果表明(原文表S11)在單葉省藤和黃藤預測的基因模型中分別有5.34%和2.89%的基因為未注釋的基因。此外，BUSCO的評估顯示，在單葉省藤和黃藤中分別有88.7%和91.3%的保守BUSCO蛋白(embryophytaodb9)存在。在這些保守的BUSCO蛋白中，分別有76.2%和81.2%是完整的。此外，tRNA，rRNA，miRNA和snRNA4種非編碼RNA也被預測出來(原文表S12)。

1.9 基因家族結構和棕櫚藤特異性基因家族分析

本研究中，通過逐一序列比對的方法來預測基因組水平的直系同源基因。這種方法因速度快而常用于處理大數(shù)據(jù)量分析。采用BLAST-based的OrthoMCL(版本2.0.9)來識別單葉省藤和黃藤中E-value≤1e-5且匹配系數(shù)≥80的同源基因。Markovchain聚類(默認膨脹參數(shù))也共同用于與其他8個物種的全面BLASTP分析，包括無油樟(Amborellatrichopoda)、油棕、擬南芥、二穗短柄草、水稻、紫萍(Spirodelapolyrhiza)、長葉刺葵和高粱。在所有10個物種中鑒定得到30936個基因家族，其中單葉省藤和黃藤基因組中分別檢測到44700和44537個同源基因。10個物種中共有的基因家族近6132個(占19.8%)，在單葉省藤和黃藤基因組中別檢測到的2366和2707個特有的基因家族(原文圖4b)。此外，家族分析結果顯示其中637個基因家族是棕櫚藤所特有的。這些特異的基因家族大多富集在膜組成成分、轉錄因子活性相關的基因簇上(原文表S13)，且聚集在植物與病原體相互作用、植物激素信號轉導相關的KEGG通路中(原文表S14)。

1.10 系統(tǒng)發(fā)育進化和分化時間分析

利用獲得的962個種間高度保守的單拷貝直系同源基因，進行棕櫚藤和其他物種的進化關系分析。首先，通過MUSCLE(版本3.8.31)進行蛋白質序列的多序列比對；在此基礎上，基于蛋白質序列進行了DNA編碼序列(CDS)比對。隨后，將所有比對的CDS序列連接起來，使用Perl腳本為每個物種生成一個超級基因。然后，提取了每個密碼子第2個位置上(相位1)的核苷酸使用RAxML(版本8.2.3)構建系統(tǒng)發(fā)育樹，模型為“GTRGAMMA”。結果表明，4種棕櫚科植物分布在一支上，由2個獨立的子分支組成，其中一個分支上包含單葉省藤和黃藤，另一個分支上包含油棕和長葉刺葵(原文圖4a)。

此外，使用了PAML(版本4.5)的MCMCTree程序并按照參數(shù)“-nsample200000-burnin40000”估計單葉省藤和黃藤與其他8個物種的分化時間。時間標準是根據(jù)參考物種發(fā)生分化的時間推算得出的。結果表明，單葉省藤和黃藤之間的分化時間約在1930萬年前，而棕櫚科的油棕和長葉刺葵的分化時間約在4080萬年前(原文圖4c)。

1.11 木質素生物合成通路基因家族全基因組鑒定分析

木質素是一種由木質素單體組成的復雜芳香族聚合物，與纖維素和半纖維素基相互作用共同構成次生細胞壁。木質素通常由3種木質素單體對-羥基苯基(p-hydroxyphenyl，H)、香草醛(vanillin，G)和丁香醛(syringaldehyde，S)組成。本研究利用擬南芥、二穗短柄草、水稻、高梁、毛竹(Phyllostachysedulis)、毛果楊(Populustrichocarpa)、單葉省藤和黃藤等8個物種的基因組進行木質素生物合成通路13個基因家族的全基因鑒定。多數(shù)基因組序列(擬南芥、二穗短柄草、水稻、高梁和毛果楊)從ENSEMBL數(shù)據(jù)庫下載獲得，從毛竹基因組數(shù)據(jù)庫(BambooGenomeDatabase)中下載了毛竹的基因組序列。在前人研究的基礎上，收集得到了經過實驗驗證的140個木質素生物合成通路基因(原文表S15)，這些基因作為查詢序列進行進一步的基因功能鑒定。如上所述，在全基因組基因鑒定過程中使用了BLAST比對和保守結構域分析方法。簡單地說，使用已知基因的編碼序列，對所有基因組序列(包括兩種棕櫚藤的基因組序列)進行了目標蛋白的BLAST比對(版本2.2.26)，選取E-value<1e-10，相似度>40%且序列的覆蓋率>95%的序列。隨后，利用Pfam-A.hmm數(shù)據(jù)庫(2017年3月31日發(fā)布)對已過濾的序列進行hmmsearch(版本3.1b2)分析，通過人工校正的方法去除結果中結構域不完整的序列。結果表明，單葉省藤和黃藤中木質素相關基因數(shù)量均有明顯增多(原文表3)。單葉省藤和黃藤中每個木質素生物合成通路基因家族平均成員數(shù)量分別約為15和13個，總數(shù)分別為193和172個。在2種棕櫚藤中過氧化物酶基因家族是最常見的基因，單葉省藤的苯丙氨酸解氨酶基因和黃藤中的香豆酸3-羥化酶基因與肉桂酸4-羥化酶基因是最不常見的。棕櫚藤中木質素生物合成基因的數(shù)目擴增可能是由于全基因組復制(WGD)事件造成的，因為全基因組復制可以提供更多的基因拷貝，從而促進具有新功能基因的進化。

2 結論

本研究采用多種類型的測序數(shù)據(jù)，首次完成了單葉省藤和黃藤染色體水平基因組的組裝，2種棕櫚藤的基因組作為重要的參考基因組，將有助于促進其他藤種的從頭基因組測序組裝和重測序研究，同時通過與不同物種進行比較研究，為物種進化研究提供證據(jù)。借助兩種棕櫚藤高質量基因組數(shù)據(jù)，使木質素生物合成通路關鍵基因的鑒定更加便捷，這些候選基因對棕櫚藤的生長發(fā)育非常重要。本研究為進一步對棕櫚藤及相關物種基因組的研究奠定了基礎。