基于高通量測序的煙草甲基因組調查分析

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

中國煙草總公司鄭州煙草研究院，鄭州高新技術產業開發區楓楊街2 號 450001

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

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

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

1 材料與方法

1.1 實驗材料

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

1.2 實驗方法

1.2.1 DNA 提取及檢測

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

1.2.2 建庫測序

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

1.2.3 K-mer 分析

采用K-mer 分析法［11］估算基因組大小。選取K 值為17 進行預測分析，統計K-mer 頻數分布，計算獲得K-mer 深度估計值，作K-mer 分布曲線。估算基因組大小，計算公式為基因組大小= K-mer數量/峰深度。將 Kdepth=1 的情況認為是錯誤情況，計算錯誤率，并用于修正基因組大小［計算公式：修正基因組大小=預估基因組大小×（1-錯誤率）］。以計算出的純合峰深度1.8 倍后面的K-mer個數所占比例來估算重復序列比例。通過雜合峰值和純合峰值比例來確定基因組的雜合率。

1.2.4 GC 含量分布及分析

針對組裝的contigs 進行GC 含量的統計，利用contigs 覆蓋深度分布與GC 含量分布構建GC-depth 點圖，并進行關聯分析。

1.2.5 基因組初步組裝

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

1.2.6 基因預測及評估

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

2 結果與分析

2.1 測序數據量及質量評估

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

2.2 K-mer 分析

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

圖1 數據質量分布Fig.1 Data quality distribution

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

表2 基因組組裝結果Tab.2 Genomic assembly results

2.4 GC 含量分布分析

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

圖3 GC 含量與測序深度關聯分析Fig.3 GC content and sequencing depth correlation analysis

2.5 基因預測及組裝質量評估

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

圖4 BUSCO 評估結果Fig.4 BUSCO assessment results

2.6 鞘翅目昆蟲基因組比較

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

從已公布的數據來看，鞘翅目昆蟲基因組大小介于 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%，高于大多數已知鞘翅目昆蟲，略低于北美眼斑叩甲（Alaus oculatus，45.6%）和Aenictocupidus jacobsonorum（44.8%），遠高于赤擬谷盜（Tribolium castaneum，35.2%）。已公布的鞘翅目昆蟲基因數量介于11 990～27 558 個之間，本研究中預測煙草甲基因數量為38 401 個。考慮到基因組survey 的測序深度不夠帶來的基因組片段化以及存在重復區域等問題，可能會導致煙草甲初步組裝后預測基因數量偏高。

3 討論

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

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

4 結論

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