基于PacBio三代測序的香瓜茄（人參果）基因組的組裝

司誠 鐘啟文 楊世鵬

（1. 青海大學農林科學院 青海省蔬菜遺傳與生理重點實驗室，西寧 810016；2. 青海大學農牧學院，西寧 810016）

香瓜茄（Solanum muricatum）是原產于南美洲的二倍體（2n=24）茄科作物［1］，因其具有抗氧化、抗糖尿病、抗炎和抗腫瘤活性的作用而聞名［2-4］。果實通常呈圓形、橢圓形或細長形，成熟時黃色的果皮上覆蓋著紫色條紋，香氣濃郁，果肉黃色帶有甜味，芳香多汁，富含大量的維生素C［5］，具有較大的營養價值、商業價值以及藥用價值。

近幾年，香瓜茄作為被消費者熟知的水果作物，其大部分研究工作主要集中在抗病機理［6］、營養價值［7］、藥用代謝物［8］、茄科近緣物種比較［9-10］等方面。目前，根據測序技術的不斷更迭，番茄（Solanum lycopersicum）［11］、 辣 椒（Capsicum annuum）［12-13］、馬 鈴 薯（Solanum tuberosum）［14］、 茄 子（Solanum melongena）［15］、煙草（Nicotiana tabacum）［16］等物種的基因組已有多個版本的完整基因組，不同番茄的近緣種（醋栗番茄 Solanum pimpinellifolium［17］、潘那利番茄Solanum pennellii［18］）基因組也通過測序被解析。與香瓜茄親緣關系最近的馬鈴薯，基于組裝了雜合二倍體馬鈴薯RH89-039-16（RH）之后［19］，四倍體馬鈴薯的第一個高質量單倍型組裝也被報道［20］。

植物基因組包含重要的遺傳信息?；蚪M序列可用于植物比較基因組學的研究，同時也是研究植物進化的資源。高質量的參考基因組有利于選擇改進農藝性狀的基因，對研究其分子機理、加速植物育種至關重要。而Hi-C技術是染色體構象捕獲技術結合高通量測序衍生的一種技術，主要研究染色體的三維結構［21］。多種作物基因組序列的公布及解析，為植物重要性狀（如果肉性狀［22］、抗性水平［23-24］等）以及豐富的基因遺傳資源的挖掘等提供了有力工具，還可據此推測基因組的進化演變，促進對關鍵農藝性狀候選基因的篩選和分子標記的開發［25］，帶動CRISPR等技術的發展及其在作物育種中的應用［26-27］，已經成為作物育種改良的重要資源和工具。

目前尚未報道香瓜茄基因組，限制了該物種的各項研究。本研究利用PacBio和Hi-C測序技術，獲得香瓜茄的基因組序列，解釋香瓜茄與近緣茄科作物的進化關系，為豐富茄科作物基因組信息及進化發育歷程，同時為香瓜茄相關分子研究奠定堅實的基礎。

1 材料與方法

1.1 材料

材料為香瓜茄甜圓形果實類型（sweet-round friut，SRF）栽培種，采集自青海大學農林科學院園藝創新基地（36°38'N，101°55'E，海拔 2 200 m），經莖尖脫毒處理后，將組培苗新鮮葉片用蒸餾水清洗干凈后，擦干，-80℃保存，送北京百邁客生物科技有限公司進行測序。

1.2 方法

1.2.1 DNA的提取 采用CTAB法提取香瓜茄植物組織DNA。

1.2.2 基因組大小的預測

1.2.2.1 文庫構建與測序 將檢測合格的DNA樣品通過Covaris超聲波破碎儀隨機打斷成片段，經末端修復、加Ploy A尾、加測序接頭、純化、PCR擴增等步驟完成整個文庫制備。構建好的文庫通過Illumina Hiseq進行PE測序。

1.2.2.2 基因組組裝 通過對raw read質控得到clean read，采用SOAPdenovo軟件進行拼接。SOAPdenovo拼接的基本過程，利用K-mer頻數表數據糾錯。對于有低頻K-mer出現的reads進行糾錯，經過糾錯之后的數據用于后續的基因組組裝。將糾錯后的小片段庫的reads截斷成更小的序列片段，構建de Brujin圖［28-29］，獲得拼接的contigs。將所有文庫測序得到的reads比對回拼接的contigs，利用reads之間的連接關系和插入片段大小信息，將contigs組裝成scaffolds。

1.2.2.3 基因組大小預估 利用Illumina HiSeq測序得到測序結果，選取Kmer=41組裝到Scaffold，通過K-mer分析初步判斷樣品的基因組大小、雜合情況、重復序列信息等評估基因組大小。

1.2.3 三代基因組文庫構建及組裝 打斷DNA樣品后對打斷的DNA樣品進行損傷修復及末端修復，連接啞鈴型接頭，進行核酸外切酶消化，使用BluePippin進行目的片段篩選，獲得測序文庫。

PacBio測序數據通過初級分析評估、過濾低質量的reads、去除接頭后得到reads，進一步堿基糾錯后得到高準確性的CCS數據，用于基因組組裝、組裝后評估等信息分析。

1.2.4 Hi-C技術輔助組裝 將香瓜茄組培苗活體植株取樣后，利用甲醛將樣品固定，將細胞內蛋白與DNA、DNA與DNA之間進行交聯，利用限制性內切酶將DNA進行酶切，利用末端修復機制引入生物素標記的堿基，將末端修復后的DNA片段進行環化、DNA解交聯及純化后，破碎為300-700 bp的片段，利用鏈親和素磁珠捕獲含有互作關系的DNA片段進行文庫構建。

文庫構建完成后，分別使用Qubit（2.0）和Agilent 2100對文庫的濃度和插入片段大小進行檢測，使用qPCR方法對文庫的有效濃度進行準確定量，以保證文庫質量。庫檢合格后，用Illumina平臺進行高通量測序，測序讀長為PE150。

1.2.5 基因組注釋 對組裝完的基因組進行基因組注釋，包括重復序列、編碼基因及功能注釋、假基因、非編碼RNA注釋等。

1.2.5.1 重復序列注釋 采用RepeatModeler2（v2.0.1）［30］、LTR_retriever（v2.8）［31］進行從頭預測（從頭 預 測 軟 件 RECON（v1.0.8）［32］和 RepeatScout（v1.0.6）［33］），RepeatClassifier［30］借 助 于 repbase（v19.06）［34］、REXdb（v3.0）［35］和 Dfam（v3.2）［36］3個已知數據庫對預測結果進行分類。將上述從頭預測結果和已知數據庫合并去冗余后得到該物種特定的重復序列數據庫，最后使用RepeatMasker（v4.1.0）［37］基于構建好的重復序列數據庫對該基因組進行轉座子序列（TE）的預測。

1.2.5.2 編碼基因預測及評估 使用Augustus（v2.4）和 SNAP［38］進行從頭預測，使用 GeMoMa（v1.7）進行基于同源物種的預測。有參的轉錄本主要使用 Hisat（v2.0.4）和 Stringtie（v1.2.3）獲得，并利用GeneMarkS-T（v5.1）進行基因預測。無參轉錄本主要通過Trinity（v2.11）［39］組裝獲得，然后使用PASA（v2.0.2）進行基因預測。最后利用EVM（v1.1.1）［40］整合上述3種方法得到的預測結果。

1.2.5.3 非編碼RNA和假基因的預測 非編碼RNA包括microRNA、rRNA和tRNA等多種已知功能的RNA，針對不同非編碼RNA的結構特點，采用了不同的策略來預測。利用tRNAscan-SE（v1.3.1）識別 tRNA，rRNA 預測主要基于 Rfam（v12.0）［41］數據庫并采用barrnap（v0.9）進行預測，miRNA通過miRBase數據庫鑒定，snoRNA和snRNA基于Rfam（v12.0）數據庫并利用Infenal（1.1）進行預測。通過GenBlastA（v1.0.4）比對，在屏蔽完真基因座的基因組上尋找同源的基因序列，然后利用GeneWise（v2.4.1）［42］尋找基因序列中的不成熟的終止密碼子及移碼突變。

1.2.6 基因功能注釋 對預測得到的基因序列進行 NR（ftp：//ftp.ncbi.nlm.nih.gov/blast/db）、KEGG（http：//www.genome.jp/kegg）、SWISS-PROT（http：//ftp.ebi.ac.uk/pub/databases/swissprot）和 Pfam［41］等數據庫的注釋分析。

2 結果

2.1 基因組預估及初步組裝

2.1.1 測序數據量統計 利用Illumina HiSeq測序得到54.26 Gb的raw reads，經質控后獲得54.11 Gb clean reads，測序深度31 X。Clean reads Q20=97.31%，Q30=92.75%，均大于90%，測序錯誤率（0.04%）＜0.05%，也在容錯范圍內，表明測序質量較好。

2.1.2 17-mer分析及基因組大小估計 通過對香瓜茄過濾得到的54.11 Gb的有效數據進行17-mer分析（圖1），根據survey分析結果，在主峰值前約1/2處（depth=15）出現一個較為明顯的小峰，說明香瓜茄基因組的雜合程度較高。主峰后depth=62處同樣也出現一個小峰，并且與主峰成倍數關系，但由于其峰值較低，峰形不明顯，應是重復序列所導致，而非同源多倍體。Depth＞62之后的拖尾則是由于香瓜茄基因組重復導致。由公式Kmer-number/depth計算得到的基因組大小約為1 252.41 Mb，修正后的基因組大小為1 238.06 Mb，基因組雜合率為0.84%，重復序列比例為65.87%（表1）。

表1 香瓜茄基因組特征Table 1 Pepino genomic characteristics

圖1 Depth和K-mer個數及種類頻率分布圖Fig. 1 Depth and number of K-mer as well species frequency distribution

2.1.3 基因組組裝結果統計 采用Soapdenovo軟件對香瓜茄序列進行拼接（表2），以Kmer=41組裝到Scaffold，contig N50為2 049 bp，總長為1 141 353 553 bp，scaffold N50為3 185 bp，總長為1 169 596 440 bp。圖2-a及2-b展示contig分布情況。

圖2 Contig覆蓋深度、長度和數量分布圖Fig. 2 Contig coverage depth, length and number distribution map

表2 組裝結果統計Table 2 Statistics of assembly results

2.1.4 GC含量及其分布 通過對組裝的contig進行GC含量統計，根據contigs的GC分布以及覆蓋深度信息繪制散點圖（圖3）。發現大多分布在20%-50%，主要集中在36%左右，經計算得到基因組GC含量為36.30%。

圖3 GC含量與測序深度（depth）關聯分析統計圖Fig. 3 Statistical analysis of GC content and sequencing depth

2.2 PacBio三代測序組裝結果及評估

2.2.1 PacBio測序結果 在PacBio測序平臺獲得香瓜茄基因組的raw reads及組裝結果（圖4）。使用該樣品的基因組DNA構建PacBio文庫，測序獲得約55 080 918 774 bp（55.08 Gb）的CCS數據，總測序深度約為47.64 X，reads N50為14 640 bp，平均讀長為14.179 bp。過濾低質量數據得到的質控后讀數共包含3 884 556條reads。

圖4 Reads長度分布統計Fig. 4 Reads length distribution statistics

2.2.2 組裝結果 PacBio數據進行質控后得到高準確性的CCS數據，然后基于CCS數據使用hifiasm（v0.12）軟件進行初步組裝，得到基因組序列（表3）。基因組序列總長度為1.15 Gb，contig N50為22.63 Mb，其中，1 kb以上contig數目1 813個，contig N90為596 645 bp，最長的contig為83 851 337 bp，GC含量為35.83%。

表3 組裝結果的統計信息Table 3 Statistical information of assembly results

2.2.3 組裝結果評估 利用bwa軟件將二代高通量測序（如Illumina HiSeq測序平臺）得到的短序列與參考基因組比對，統計比對率（99.85%），可評估組裝基因組的完整性。使用CEGMA（v2.5）來評估最終基因組組裝的完整性，有97.16%的CEGMA基因存在香瓜茄基因組中。使用BUSCO軟件評估基因組組裝完整性，有98.20%的BUSCO基因存在香瓜茄基因組中，表明基因組組裝完整性較高。

2.3 Hi-C輔助基因組組裝及評估

對初步組裝的基因組序列利用有效Hi-C數據進行進一步的組裝，包括初步組裝基因組序列的分群、排序和排序后的定向，最終獲得染色體水平上基因組序列。共產生143 362 025 128 bp數據。通過Hi-C文庫質量評估分析共獲得404 890 166對唯一比對到基因組上的reads（unique paired alignments），其中，252 453 038對是有效的Hi-C數據（valid interaction pairs），占唯一比對到基因組數據的62.35%；116 720 788對無效數據中末端懸掛類型的Hi-C數據（dangling end pairs），占唯一比對到基因組數據的28.83%；6 442 178對無效數據中屬于相鄰連接類型的（re-ligation pairs），占比1.59%；1 535 133對無效數據中屬于為自連類型的（self-circle ligation pairs），占比0.83%；27 739 029對無效數據中屬于其他未定義的（dumped pairs），占比6.85%。

經過Hi-C組裝和人工調整后，共有1 123 245 570 bp的序列長度的基因組序列被定位到12條染色體上，占比97.16%；在定位到染色體上的序列中，能夠確定順序和方向的序列長度為1 079 553 436 bp，占定位染色體序列總長度的96.11%。對Hi-C糾錯和組裝后得到的基因組序列進行統計（表4），獲得最終版本的基因組組裝統計結果，Contig N50為22 628 432 bp，Scaffold N50為87 253 278 bp。

表4 香瓜茄Hi-C組裝的基因組信息Table 4 Hi-C assembly information of the pepino genome

對于Hi-C組裝到染色體的基因組等長切割成500 000 bp一個bin，然后任意2個bin之間覆蓋Hi-C Read Pairs的數目作為2個bin之間交互的強度信號（圖5），可以明顯區分出12個染色體分組；在每一分組內部可以看出位于對角線位置的交互的強度要高于非對角線的位置，說明Hi-C組裝的染色體結果中鄰近的序列間（對角線位置）交互強度高，而非鄰近的序列之間（非對角線位置）的交互信號強度弱，與Hi-C輔助基因組組裝的原理一致，證明香瓜茄基因組序列掛載率高。

圖5 香瓜茄基因組Hi-C組裝染色體交互熱圖Fig. 5 Hi-C assembly chromosome interaction heat map of pepino genome

2.4 基因預測結果

對組裝完的基因組進行基因組注釋，包括重復序列、編碼基因及功能注釋、假基因、非編碼RNA注釋等。重復序列注釋主要包括串聯重復序列（tandem repeats）和散在重復序列（interspersed repeats），其中，第二類主要是轉座子序列（transposable elements，TE）是研究的主要對象。將從頭預測結果和已知數據合并去冗余后得到該物種特定的重復序列數據庫，最后基于構建好的重復序列數據庫對香瓜茄基因組進行TE的預測。最終得到約742 491 882 bp的TE，占比64.22%，最終得到約201 341 835 bp的串聯重復序列，占比17.42%。

對香瓜茄基因采用同源預測、從頭預測和轉錄組預測，基因預測結果（表5）顯示，編碼基因預測最終得到41 571個基因；非編碼RNA即不編碼蛋白質的RNA，包括miRNA、rRNA和tRNA等多種已知功能的RNA，針對不同非編碼RNA的結構特點，采用了不同的策略來預測不同的非編碼RNA，總共預測得到4 360個tRNA、5 677個rRNA、154個miRNA、202個snRNA、287個snoRNA；假基因預測得到449個。利用擬南芥、辣椒、番茄、潘那利番茄以及馬鈴薯等開展同源預測香瓜茄基因信息，其中，香瓜茄與近緣作物馬鈴薯預測得到的基因數量最多，有51 586個。

表5 香瓜茄基因預測結果Table 5 Prediction results of pepino gene

BUSCO中embryophyta數據庫包含1 614個保守的核心基因。使用BUSCO（v4.0）軟件來評估基因預測的完整性，其中，有98.64%的BUSCO基因存在預測的基因中，說明基因預測的完整性高。

2.5 基因功能注釋

香瓜茄中99.06%的基因可以注釋到所有數據庫中（表6）。通過GO注釋分析（圖6），共有30 713個基因具有GO注釋預測的功能，占預測到總基因數的73.88%。GO注釋結果顯示整個分類中基因分布在細胞組分（cellular component）的相較于分子功能（molecular function）和生物學過程（biological process）較少，生物學過程最多。其中，二級功能分布在細胞內（intracellular）、細胞結構體（cellular anatomical entity）、催化活動（catalytic activity）、結合（binding）、代謝過程（metabolic process）、細胞過程（cellular process）的基因數目相對較多。

圖6 GO二級節點注釋分類統計圖Fig. 6 Statistical chart of GO secondary node annotation classification

表6 香瓜茄基因功能注釋統計信息Table 6 Statistical information of pepino gene function annotation

eggNOG注釋結果（圖7）顯示，香瓜茄的蛋白序列功能主要集中在復制、重組和生物生成（L：replication，recombination and repair），占比 10.94%，轉錄（K：transcription）占比7.58%，信號轉導機制（T：signal transduction mechanisms）占比7.1%，翻譯后修飾、蛋白質周轉、分子伴侶（O：posttranslational modification，protein turnover，chaperones） 占 比6.82%，能量的產生和轉換（C：energy production and conversion）占比5.36%。eggNOG采用了COG，KOG和arCOG中引入的20個功能類別，在功能層面上對基因進行分類。eggNOG結果反映在不同的功能類別中，通過基因數目的多少能夠展示出該物種在進化過程中對環境的適應性。

圖7 eggNog功能注釋分類統計圖Fig. 7 eggNog functional annotation classification

2.6 物種進化分析

選擇8個已知基因組信息的物種，構建系統進化樹。結果表明，香瓜茄與馬鈴薯的進化時間大約在12.82 MYA（圖8）。從進化時間上來看，香瓜茄的進化時間稍晚于煙草、辣椒和茄子。對比于基因預測結果（表5），近緣物種馬鈴薯作為香瓜茄同源預測物種，預測到的香瓜茄上的基因個數也最多，揭示了馬鈴薯與香瓜茄較近的進化關系。

圖8 物種間分化時間Fig. 8 Differentiation time between species

3 討論

連續性和完整性是基因組組裝的重要指標，PacBio基因組組裝和Illumina數據的糾錯可以大大提高測序數據連續性和完整性［43］。本研究通過這種策略對香瓜茄基因組組裝顯示出高度分辨的結果，N50=22.62 Mb（megabases），與近期測序的芒苞草 N50 = 6.96 Mb［44］，黑麥 N50＞29 Mb［45］和板藍根N50=36.16 Mb［46］結果相近。在定位到染色體上的序列中，能夠確定順序和方向的序列長度為1 079 Mb，占定位染色體序列總長度的96.11%。BUSCO和基因預測結果分析進一步證實了香瓜茄基因組的高質量和完整性。

根據基因組數據構建了香瓜茄的系統發育關系，同時對比基因預測結果，香瓜茄比對到馬鈴薯上的基因個數最多，發現在茄科作物中，香瓜茄與馬鈴薯進化關系最近，這與先前發表的系統發育分析的結果一致［9］。茄科作為雙子葉植物中最重要的果蔬類群，其包含的較多物種的基因組或基因組草圖已經被測序完成［47-48］?；蚪M測序技術及生物信息技術的不斷發展［49］，顯著推動了香瓜茄這種非模式作物的測序研究，香瓜茄基因組的測序完成是對茄科基因組研究的又一補充。我們的研究結果將為茄科的起源，進化和多樣化分析增加了功能見解。

本研究通過對香瓜茄基因組的預估，對比PacBio第三代測序技術測序以及Hi-C輔助基因組組裝結果，首次揭示了香瓜茄基因組的大小。根據各分析指標，推測香瓜茄基因組為高雜合基因組，針對一些物種基因組重復序列偏多的特征，可以采用三代測序或者HiFi測序等兼顧長讀長和高精度的測序手段開展基因組研究。隨著長度測序的出現和完善，基因組組裝的數量和質量正在不斷提升，但一些具有顯著生態價值和較低經濟價值的植物中參考基因組的數量和質量仍然較低，在137個植物目中，有76個植物目缺乏代表性的參考基因組，62個目至少有1個參考基因組。例如，十字花科目有83個種的參考基因組，禾本目和唇形目分別有80個種和67個種的參考基因組。伴隨著技術的進步和越來越多其他物種的關注度提升，未來完整的植物基因組測序數據庫構建將成為可能。

4 結論

獲得香瓜茄高質量染色體水平參考基因組，推測該測序香瓜茄基因組為高復雜基因組。香瓜茄與馬鈴薯具有較近的進化關系。