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

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

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

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

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

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

1 材料與方法

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

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

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

1.2.1 DNA 提取及檢測(cè)

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

1.2.2 建庫(kù)測(cè)序

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

1.2.3 K-mer 分析

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

1.2.4 GC 含量分布及分析

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

1.2.5 基因組初步組裝

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

1.2.6 基因預(yù)測(cè)及評(píng)估

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

2 結(jié)果與分析

2.1 測(cè)序數(shù)據(jù)量及質(zhì)量評(píng)估

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

2.2 K-mer 分析

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

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

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

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

2.4 GC 含量分布分析

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

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

2.5 基因預(yù)測(cè)及組裝質(zhì)量評(píng)估

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

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

2.6 鞘翅目昆蟲(chóng)基因組比較

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

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

3 討論

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

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

4 結(jié)論

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