基于高通量測序的煙草甲基因組調(diào)查分析

許亞龍，金靜靜，趙艷珍，魏 攀，奚家勤，楊 軍，曹培健，張劍鋒

中國煙草總公司鄭州煙草研究院，鄭州高新技術(shù)產(chǎn)業(yè)開發(fā)區(qū)楓楊街2 號 450001

煙 草 甲（Lasioderma serricorne），屬 鞘 翅 目（Coleoptera）竊蠹科（Anobiidae），是一種雜食性倉儲害蟲，其寄主范圍廣泛，主要為害儲藏的煙草及其制品、儲藏的糧食以及藥材等，在全世界均有分布［1-2］。由于煙草甲主要是通過幼蟲潛居在寄主的體內(nèi)進(jìn)行蛀食，其發(fā)生為害具有較強的隱蔽性，對儲藏物的品質(zhì)造成了嚴(yán)重的影響。針對煙草甲為害損失的調(diào)查表明，我國卷煙工業(yè)企業(yè)的蟲害直接損失率約為 0.215%［3］。

DNA 測序技術(shù)的快速發(fā)展為昆蟲基因組學(xué)的興起奠定了堅實的技術(shù)基礎(chǔ)。5000 種昆蟲全基因組測序計劃（5000 Insect Genome Project, i5K）［4］和千種昆蟲轉(zhuǎn)錄組進(jìn)化項目（1K Insect Transcriptome Evolution,1KITE）［5］相繼啟動，目標(biāo)涵蓋昆蟲全基因組測序、轉(zhuǎn)錄組學(xué)、功能基因組學(xué)、進(jìn)化生物學(xué)、比較基因組學(xué)、生物信息學(xué)分析技術(shù)等研究方向，標(biāo)志著昆蟲學(xué)研究已全面進(jìn)入基因組學(xué)時代。2008 年，赤擬谷盜（Tribolium castaneum）成為首個完成基因組測序的鞘翅目昆蟲［6］，此后相繼對多種鞘翅目昆蟲進(jìn)行了全基因組測序［7-9］。煙草甲屬鞘翅目竊蠹科，但是在NCBI 已公布的29 種鞘翅目昆蟲基因組中未發(fā)現(xiàn)竊蠹科昆蟲相關(guān)數(shù)據(jù)。煙草甲目前僅見線粒體基因組的全序列報道［10］，其基因組信息的匱乏嚴(yán)重制約了其分子生物學(xué)的研究進(jìn)展。因此，針對煙草甲進(jìn)行全基因組測序就能夠獲得其基因組信息，支撐煙草甲生長發(fā)育、生理習(xí)性以及生物防治等分子機制的研究。

本研究中基于高通量測序的基因組survey 對煙草甲基因組大小及復(fù)雜程度等重要的基因組特征進(jìn)行了初步分析，結(jié)合煙草甲基因組特點對后續(xù)測序提出合理方案，旨在為進(jìn)一步解析高質(zhì)量煙草甲全基因圖譜奠定基礎(chǔ)。

1 材料與方法

1.1 實驗材料

煙草甲在中國煙草總公司鄭州煙草研究院煙草倉貯實驗室采集并在人工培養(yǎng)箱飼養(yǎng)。飼養(yǎng)食料：90%全麥粉+10%酵母粉；飼養(yǎng)條件：溫度28 ℃± 2 ℃，相對濕度70% ± 5%，暗處理。取雌雄成蟲各一只，連續(xù)繁殖純化4 代以上。取后代幼蟲10頭，利用超純水沖洗去除蟲體上粘附的飼料，經(jīng)液氮速凍后保存于超低溫冰箱中備用。

1.2 實驗方法

1.2.1 DNA 提取及檢測

利用DNA 提取試劑盒（Insect gDNA Isolation Kit, 美國Biomiga 公司）提取煙草甲基因組DNA，利用瓊脂糖凝膠電泳對所提DNA 的完整性、純度、片段大小進(jìn)行檢測，利用Qubit 熒光計（Invitrogen Qubit 2.0, 美國Thermo Fisher Scientific公司）進(jìn)行DNA 濃度測定。確保DNA 質(zhì)量達(dá)到建庫測序要求。

1.2.2 建庫測序

質(zhì)檢合格的DNA 樣本委托北京諾禾致源科技股份有限公司進(jìn)行建庫測序。DNA 樣品經(jīng)超聲破碎后隨機打斷，構(gòu)建小片段（180 bp、300 bp、500 bp）文庫。通過Illumina Hiseq 2000 進(jìn)行PE150 雙末端測序。針對測序數(shù)據(jù)進(jìn)行質(zhì)控，過濾去除掉低質(zhì)量的數(shù)據(jù)，利用有效數(shù)據(jù)進(jìn)行煙草甲基因組特征評估及初步組裝。

1.2.3 K-mer 分析

采用K-mer 分析法［11］估算基因組大小。選取K 值為17 進(jìn)行預(yù)測分析，統(tǒng)計K-mer 頻數(shù)分布，計算獲得K-mer 深度估計值，作K-mer 分布曲線。估算基因組大小，計算公式為基因組大小= K-mer數(shù)量/峰深度。將 Kdepth=1 的情況認(rèn)為是錯誤情況，計算錯誤率，并用于修正基因組大?。塾嬎愎剑盒拚蚪M大小=預(yù)估基因組大小×（1-錯誤率）］。以計算出的純合峰深度1.8 倍后面的K-mer個數(shù)所占比例來估算重復(fù)序列比例。通過雜合峰值和純合峰值比例來確定基因組的雜合率。

1.2.4 GC 含量分布及分析

針對組裝的contigs 進(jìn)行GC 含量的統(tǒng)計，利用contigs 覆蓋深度分布與GC 含量分布構(gòu)建GC-depth 點圖，并進(jìn)行關(guān)聯(lián)分析。

1.2.5 基因組初步組裝

利用SOAPdenovo2［12］對測序數(shù)據(jù)進(jìn)行拼接組裝，將測序所得reads 進(jìn)行比對得到的contigs。根據(jù)雙末端數(shù)據(jù)之間的配對關(guān)系連接contigs，將contigs 組裝成 scaffolds，并對 contigs 之間的空隙進(jìn)行補全，得到原始基因組序列。

1.2.6 基因預(yù)測及評估

利用Augustus［13］基于赤擬谷盜的基因訓(xùn)練集對初步組裝的基因組進(jìn)行基因位置的注釋。利用BUSCO［14］選擇真核模式生物中的255 個保守基因作為參考數(shù)據(jù)庫對基因組組裝質(zhì)量進(jìn)行評估。

2 結(jié)果與分析

2.1 測序數(shù)據(jù)量及質(zhì)量評估

基于Illumina Hiseq 平臺進(jìn)行雙端PE150 測序，過濾掉無效或低質(zhì)量的reads 數(shù)據(jù)，共獲得煙草甲中reads 數(shù)量為163 929 635 條，測序總數(shù)據(jù)量為49.18 GB。Q20 與Q30 均為衡量測序質(zhì)量優(yōu)劣的指標(biāo)，本研究中煙草甲高通量測序Q20 比率達(dá)97.00%、Q30 比率達(dá)93.10%，表明煙草甲基因組高通量測序數(shù)據(jù)質(zhì)量較高（圖1）。煙草甲基因組測序數(shù)據(jù)中位置堿基 N 基本為零，A 與 T、C 與 G 的互補堿基數(shù)基本一致，表明本研究中煙草甲基因組的測序質(zhì)量較好。

2.2 K-mer 分析

利用K-mer 的分析方法來預(yù)測煙草甲基因組的大小、雜合率和重復(fù)序列等基因組特征。當(dāng)取K=17 時，根據(jù)SOAP de novo 軟件預(yù)測得到K-mer總數(shù)為43 906 084 422（表1）。根據(jù)圖2 中K-mer的深度分布，根據(jù)公式估算出煙草甲基因組大小為245.29 Mb，經(jīng)修正后的基因組大小為242.25 Mb。根據(jù)計算公式，煙草甲基因組雜合率為0.77%，重復(fù)序列比率為42.95%。

圖1 數(shù)據(jù)質(zhì)量分布Fig.1 Data quality distribution

表1 K-mer 分析所得基因組特征統(tǒng)計分析Tab.1 K-mer analysis of genome features

圖 2 Depth 和 K-mer 頻率分布圖Fig.2 Depth and K-mer frequency distribution

2.3 基因組初步組裝

利用SOAP de novo 軟件對煙草甲測序數(shù)據(jù)進(jìn)行初步組裝?？紤]到基因組存在的雜合和重復(fù)的情況，以K=41 作為初步組裝的K-mer 值，首先組裝成為contigs，得到比較理想的組裝結(jié)果（表2）。針對組裝好的長度大于等于100 bp 的contigs 進(jìn)行統(tǒng)計，N50 長度為1 309 bp，組裝得到最長的序列長度為678 872 bp，組裝的contigs 總數(shù)量為461 378條，總長度為206.74 Mb。進(jìn)一步將所有文庫測序得到的reads 比對到初步得到的contigs，利用reads之間的連接關(guān)系和插入片段大小信息，最終將contigs 組裝成scaffolds。統(tǒng)計分析發(fā)現(xiàn)，scaffold的N50 長度為1 864 bp，最長序列長度為1 500 785 bp，組裝的 scaffolds 總量為 418 693 條，總長度為211.10 Mb。

表2 基因組組裝結(jié)果Tab.2 Genomic assembly results

2.4 GC 含量分布分析

GC 含量是基因組核酸序列組成的重要特征，GC 含量-測序深度關(guān)聯(lián)分析可以用于檢測樣本基因組是否存在GC 分布偏好以及是否存在外源的污染等。針對組裝的contigs 進(jìn)行GC 含量的統(tǒng)計，進(jìn)行了GC 含量與測序深度的關(guān)聯(lián)分析。如圖3所示，橫坐標(biāo)表示GC 含量，縱坐標(biāo)表示測序深度，右方是contigs 覆蓋深度分布，上方是GC 含量分布。GC 含量主要集中在窗口的30%～50%之間，表明煙草甲基因組沒有顯著的GC 偏好性。GC 含量也沒有顯著的分層現(xiàn)象，表明煙草甲基因組的雜合率不高。計算分析發(fā)現(xiàn)，煙草甲基因組初步組裝版本GC 含量為44.61%。圖中低深度區(qū)出現(xiàn)了少量的GC 聚集，經(jīng)NCBI 核苷酸數(shù)據(jù)庫blast 比對分析發(fā)現(xiàn)，低深度區(qū)域部分存在少量小麥和細(xì)菌等污染。這可能與樣本采集前飼喂全麥粉有關(guān)。此外，煙草甲基因組中雜合度為0.77%。由于在組裝過程中同源染色體上雜合部位只能被識別出一半，導(dǎo)致在低測序深度區(qū)域也出現(xiàn)了GC 富集的現(xiàn)象。

圖3 GC 含量與測序深度關(guān)聯(lián)分析Fig.3 GC content and sequencing depth correlation analysis

2.5 基因預(yù)測及組裝質(zhì)量評估

為了在煙草甲基因組中較為準(zhǔn)確地預(yù)測基因序列，利用Augustus 選擇鞘翅目模式昆蟲赤擬谷盜為基因模型物種，預(yù)測煙草甲基因組中基因，并對初步組裝的基因組進(jìn)行基因注釋。分析發(fā)現(xiàn)，煙草甲中預(yù)測基因數(shù)量為38 401 個。以真核生物中保守基因作為參考數(shù)據(jù)庫，以單拷貝基因拼接的完整性和準(zhǔn)確性來評價煙草甲基因組組裝質(zhì)量。BUSCO 分析發(fā)現(xiàn)，初步組裝的煙草甲基因組可完整覆蓋89.8%的BUSCO 核心基因（n=255）（圖4），其中71.37%的基因為單拷貝，18.43%的基因包含多拷貝，7.45%的基因部分覆蓋，僅有2.75%的基因未能比對上。以上結(jié)果表明煙草甲基因組具備較高的完整性和準(zhǔn)確性，基因組組裝質(zhì)量較好。

圖4 BUSCO 評估結(jié)果Fig.4 BUSCO assessment results

2.6 鞘翅目昆蟲基因組比較

在 NCBI 的 Genome 數(shù)據(jù)庫（https://www.ncbi.nlm.nih.gov/genome）中查詢已公布鞘翅目昆蟲基因組信息，得到29 種鞘翅目昆蟲的基因組信息，將其基因組信息與本研究中獲得的煙草甲基因組數(shù)據(jù)進(jìn)行比較，結(jié)果見表3。

從已公布的數(shù)據(jù)來看，鞘翅目昆蟲基因組大小介于 12.08～2 418.07 Mb 之間，相差近 200 倍。煙草甲基因組242.25 Mb 處于一個相對較小的水平，與蜂房小甲蟲（Aethina tumida，234.34 Mb）和沙漠鐵包甲蟲（Asbolus verrucosus，249.61 Mb）較為接近，約為模式昆蟲赤擬谷盜（Tribolium castaneum，165.94 Mb）的1.5 倍。煙草甲基因組的GC 含量為44.6%，高于大多數(shù)已知鞘翅目昆蟲，略低于北美眼斑叩甲（Alaus oculatus，45.6%）和Aenictocupidus jacobsonorum（44.8%），遠(yuǎn)高于赤擬谷盜（Tribolium castaneum，35.2%）。已公布的鞘翅目昆蟲基因數(shù)量介于11 990～27 558 個之間，本研究中預(yù)測煙草甲基因數(shù)量為38 401 個?？紤]到基因組survey 的測序深度不夠帶來的基因組片段化以及存在重復(fù)區(qū)域等問題，可能會導(dǎo)致煙草甲初步組裝后預(yù)測基因數(shù)量偏高。

3 討論

全基因組測序是破譯物種遺傳密碼的重要基礎(chǔ)。在啟動物種全基因組測序工作之前，有必要對其基因組大小及復(fù)雜程度進(jìn)行初步評估，從而確定對應(yīng)的全基因組測序研究方案?；蚪M大小的預(yù)測常使用流式細(xì)胞術(shù)［15］、Feulgen 圖像分析法［16］、基因組 survey 分析［17］等方法。相比于其他方法，基因組survey 分析是一種更為精確的分析未知基因組特征的方法。除此之外，通常認(rèn)為基因組雜合度越大，重復(fù)片段越多，該物種的組裝難度就越大。本研究中，煙草甲基因組雜合率為0.77%，重復(fù)序列比例高達(dá)42.95%，高質(zhì)量基因組組裝難度較大。隨著近年來測序成本的下降和三代測序技術(shù)的普及，采用二代Illumina 測序結(jié)合三代PacBio RSII 測序策略，輔以Hi-C 技術(shù)進(jìn)行煙草甲全基因組測序研究，有望獲得高質(zhì)量的煙草甲全基因組圖譜。

表3 煙草甲基因組組裝數(shù)據(jù)與鞘翅目29 種昆蟲基因組比較Tab.3 Genome assembly statistics of L. serricorne and comparisons to 29 genomes of Coleoptera

4 結(jié)論

通過對煙草甲進(jìn)行全基因組survey 分析，預(yù)估煙草甲基因組大小為242.25 Mb，GC 含量為44.61% ，雜 合 率 為 0.77% ，重 復(fù) 序 列 比 例 為42.95%；組裝后得到的 contig N50 為 1 309 bp，總長為 206.74 Mb，scaffold N50 為 1 864 bp，總長為211.10 Mb；預(yù)測基因數(shù)量為38 401 個。