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

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

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

DNA 測序技術(shù)的快速發(fā)展為昆蟲基因組學(xué)的興起奠定了堅實(shí)的技術(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í)性以及生物防治等分子機(jī)制的研究。

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

1 材料與方法

1.1 實(shí)驗(yàn)材料

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

1.2 實(shí)驗(yàn)方法

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)超聲破碎后隨機(jī)打斷，構(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 點(diǎn)圖，并進(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)行補(bǔ)全，得到原始基因組序列。

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 的互補(bǔ)堿基數(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%的基因?yàn)閱慰截悾?8.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)的全基因組測序研究方案。基因組大小的預(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 個。