枯葉蛺蝶中華亞種基因組Survey分析

洪芳，徐堉峰，范博

(1.泉州師范學(xué)院 海洋與食品學(xué)院，福建 泉州 362000;2.臺(tái)灣師范大學(xué) 生命科學(xué)系，臺(tái)灣 臺(tái)北 116)

枯葉蝶，學(xué)名枯葉蛺蝶(Kallimainachus)，隸屬于蛺蝶科Nymphalidae，是大型蝴蝶，以前后翅相疊其翅形及斑紋似枯葉而著稱，具有較高的仿生學(xué)研究價(jià)值[1].枯葉蛺蝶[2]是枯葉蛺蝶屬中分布最為廣泛的物種，為該屬的代表種.目前，枯葉蛺蝶已知有6個(gè)亞種 (不含海南亞種)，包括指名亞種K.inachusinachus[2]、中華亞種K.inachuschinensisSwinhoe、臺(tái)灣亞種K.inachusformosanaFruhstorfer、日本亞種K.inachuseucercaFruhstorfer、暹羅亞種K.inachussiamensisFruhstorfer和北印度亞種K.inachushuegeli(Kollar).枯葉蝶一年多代，成蟲常棲息于闊葉林，出現(xiàn)在濕潤繁茂雨林里的灌木叢中和河床兩岸，在南方幾乎全年可見[3].枯葉蛺蝶棲息于天然植被覆蓋較好的山區(qū)，其幼蟲寄主植物有蕁麻科Urticaceae、蓼科Polygonaceae、薔薇科Rosaceae和爵床科Acanthaceae[4]；在我國僅記載有爵床科Acanthaceae蘆莉花族Ruellieae、鱗花草族Lepidagathideae和爵床族Justicieae的部分植物種類[5].幼蟲能以爵床科多種馬蘭屬植物為寄主.成蟲則以腐爛水果、闊葉樹干蟲蛀傷口流出的樹液和動(dòng)物糞便為食，是典型的食腐蝶類.垂直棲息海拔根據(jù)環(huán)境而定，一般可達(dá)到1 800 m的高度[6]，盡管Mark Alexander Wynter-Blyth記錄了它在海拔2 400 m森林茂密的山區(qū)和丘陵地區(qū)會(huì)遇到[4].在喜馬拉雅山脈印度庫曼地區(qū)，枯葉蛺蝶棲息在400～1 400 m的熱帶落葉森林和1 200 m以上的亞熱帶常綠森林[6].1998-2004年在中國重慶市進(jìn)行的一項(xiàng)調(diào)查發(fā)現(xiàn)，枯葉蝶棲息于潮濕闊葉林中[7].枯葉蛺蝶廣泛分布于東南亞的國家和地區(qū)，如印度、尼泊爾、不丹、孟加拉國、緬甸、中國、泰國、老撾、越南、巴基斯坦和日本[6].在我國分布于除海南外的秦嶺以南地區(qū).

當(dāng)前，第二代測序技術(shù) (next-generation sequencing technologies,NGS)和PacBio的單分子實(shí)時(shí)(single molecular real time,SMRT) 測序技術(shù)被廣泛開發(fā)和應(yīng)用，并迅速占領(lǐng)基因組測序領(lǐng)域的主流行業(yè).Survey評(píng)估中的K-mer分析能得到更為全面和準(zhǔn)確的信息，是因?yàn)镵-mer分析是從數(shù)學(xué)角度進(jìn)行分析的，克服了傳統(tǒng)Feulgen分光光度法[8]和流式細(xì)胞法用來測定基因組大小[9-10]的劣勢.由于內(nèi)參選擇不同，實(shí)驗(yàn)評(píng)估基因組大小與實(shí)際會(huì)有一些偏差的不足.第二代測序技術(shù)因其快速發(fā)展及推廣，進(jìn)而加速了全基因組的研究進(jìn)程[11].通過對(duì)枯葉蝶的全基因組測序，可為今后仿生學(xué)領(lǐng)域的研究提供科學(xué)理論依據(jù).

1 材料與方法

1.1 材料

本研究所用材料為枯葉蛺蝶中華亞種，2017年6月購于四川省樂山市枯葉蝶養(yǎng)殖基地.所用活體樣本帶回實(shí)驗(yàn)室，液氮速凍后置于超低溫冰箱(-80 ℃)保存，備用.

1.2 方法

1.2.1 基因組DNA提取、檢測及測序 首先，采用改良的SDS法提取枯葉蛺蝶樣本的基因組DNA，用紫外分光光度計(jì) (日本島津,UV-1800)檢測濃度，電泳的Marker為1 kb DNA Ladder(TaKaRa)和λ-HindⅢdi-gest(TaKaRa,日本)；然后，隨機(jī)打斷枯葉蝶的DNA樣品，構(gòu)建長度為250 bp的小片段文庫，再通過末端修復(fù)、加測序接頭盒A尾、PCR擴(kuò)增、純化等處理后完成文庫制備.構(gòu)建好的文庫委托北京諾禾致源科技股份有限公司完成.利用Illumina Hiseq 2500測序平臺(tái)進(jìn)行雙末端測序.為了保證分析數(shù)據(jù)的準(zhǔn)確性，首先對(duì)測序產(chǎn)生的原始序列(raw reads)進(jìn)行過濾，去除接頭污染等低質(zhì)量的片段，其后將得到的高質(zhì)量序列(clean reads)用于每個(gè)物種的基因組大小預(yù)測以及雜合度和GC含量分析.

1.2.2 基因組特征預(yù)估 在基因組組裝前，通過測序得到的序列估計(jì)基因組特征.采用基于K-mer的分析方法和結(jié)果來估計(jì)基因組大小和雜合率等，這里取K=17[12]來進(jìn)行分析.

(1)根據(jù)Lander_waterman算法，基因組大小(G)滿足如下公式：

Cbase=Ck-mer×L/(l-k+1)，G=nk-mer/Ck-mer=nbase/Cbase.

其中：nbase和nk-mer為序列堿基總數(shù)和K-mer數(shù)，Cbase和Ck-mer為覆蓋堿基的期望深度和K-mer期望覆蓋深度.

(2)K-mer深度頻率分布服從泊松分布：

其中：a1/2為雜合K-mer種類數(shù)的百分比，為所有K-mer的種類數(shù).

1.2.3 基因組初步組裝 通過所有小片段庫測序得到的reads截?cái)喑筛〉男蛄衅沃g的重疊關(guān)系構(gòu)建de Brujin圖.去掉無法繼續(xù)連接的分支、低覆蓋度的分支，并且利用reads信息化簡單重復(fù)序列在de Brujin圖的分叉通路，才能簡化非常復(fù)雜的de Brujin圖.采用隨機(jī)選擇策略，合并少量的雜合位點(diǎn)，得到一個(gè)簡化后的de Brujin圖.只有在每個(gè)分叉位點(diǎn)將序列截?cái)?，才能得到最初的contigs，是因?yàn)楹喕蟮膁e Brujin圖仍然會(huì)有很多分叉位點(diǎn)無法確定真正的連接關(guān)系.根據(jù)插入片段大小信息和reads之間的連接關(guān)系，將所有文庫測序得到的reads比對(duì)回初步得到的contigs，將contigs組裝成scaffolds.構(gòu)建contig 和scaffold采用K=41，進(jìn)行SOAP de novo組裝[13]要根據(jù)高質(zhì)量數(shù)據(jù)，在得到scaffold序列后，要先用SOAP將過濾后的reads 比對(duì)到該組裝序列上，再直接拼接，從而獲得原始基因組序列及堿基深度.以10 kb為窗口，在序列上無重復(fù)前進(jìn)，從而做出GC_depth點(diǎn)圖.計(jì)算每個(gè)窗口的平均深度與GC含量，利用GC 聚成塊的分層來判斷基因組的雜合率和重復(fù)的分布情況.

2 結(jié)果與分析

2.1 測序數(shù)據(jù)量

使用第二代高通量測序技術(shù)測序，共獲取28.25 Gb原始數(shù)據(jù).若預(yù)計(jì)樣品基因組大小為545.05 M，則枯葉蝶樣品的測序深度為52×.先建立一個(gè)DNA文庫， 插入片段大小350 bp，嚴(yán)格過濾后得到clean data.產(chǎn)出初始數(shù)據(jù)28 249 103 100 bp，過濾后28 183 572 900 bp，GC含量33.94%，測序錯(cuò)誤率0.02%.

2.2 K-mer分析

枯葉蛺蝶的原始數(shù)據(jù)用于17-mer分析使用，得出其頻率分布(圖1).K-mer的總數(shù)是21 588 M.以17-mer出現(xiàn)的次數(shù)作為橫坐標(biāo)，出現(xiàn)的頻率作為縱坐標(biāo)，可以看出，17-mer分布曲線成峰情況較好.K-mer 深度分布受基因組中雜合子和重復(fù)序列存在的影響.因?yàn)镵-mer曲線呈現(xiàn)明顯拖尾(圖1)，說明枯葉蝶基因組重復(fù)序列含量較高.從survey分析的結(jié)果可以看出，在depth=39附近是主峰值，由K-mer-number/depth)得到的基因組大小為553.55 Mb左右，修正后的基因組大小為545.05 Mb(圖2).經(jīng)過計(jì)算和分析得出，枯葉蝶的基因組雜合率為1.22%，重復(fù)序列比例為48.03.%.

圖1 K-mer=17 Depth和K-mer個(gè)數(shù)頻率分布 圖2 K-mer=17 Depth和K-mer種類數(shù)頻率分布Fig.1 Frequency distribution of K-mer =17 Fig.2 The frequency distribution of K-mer=17 Depth and K-mer number Depth and K-mer species number

2.3 基因組初步組裝

用K-mer 41進(jìn)行初步組裝，統(tǒng)計(jì)得到Contig N50為555 bp，總長為653 208 907 bp，Scaffold N50為972 bp，總長為708 125 193 bp(表1).由各分析指標(biāo)來看，測序枯葉蝶基因組為復(fù)雜基因組，后續(xù)可以采用相應(yīng)策略進(jìn)行基因組組裝.

圖3 contig GC含量和覆蓋深度Fig.3 GC content and average sequencing depth of the Kallima inachus genome data 注: 紅色的部分代表散點(diǎn)圖中點(diǎn)的密度比較大的部分

表1 組裝結(jié)果統(tǒng)計(jì)Tab.1 Statistics of assembly results in the Kallima inachus genome bp

2.4 GC含量及分布

在全基因組尺度上，基因組中的堿基組成，腺嘌呤(A)、鳥嘌呤(G)、胸腺嘧啶(T)、胞嘧啶(C)的相對(duì)含量，一般是用GC含量(GC-content)來表示[14-15].因GC含量的變化跟重復(fù)元件分布、基因密度、同義密碼子的使用頻率等密切相關(guān)，故其含量能體現(xiàn)生物體核酸序列組成的重要特征[16].

經(jīng)過計(jì)算預(yù)估，枯葉蛺蝶基因組GC含量為35.24%.對(duì)組裝的contigs進(jìn)行GC含量統(tǒng)計(jì)的結(jié)果見圖3.由圖中可知，枯葉蛺蝶的GC深度分布被分成了2層，出現(xiàn)了一個(gè)低深度區(qū)域，基因組的雜合率(1.97%)高是可能的原因.因?yàn)殡s合會(huì)導(dǎo)致兩條同源染色體雜合的部位只組裝出了1條，或2條都沒有裝出；同時(shí)，GC含量出現(xiàn)較低的一層，正是由于該部位以上的read乘數(shù)是整個(gè)基因組乘數(shù)的一半[17].根據(jù)圖3中GC含量的深度信息，發(fā)現(xiàn)GC圖發(fā)生輕微的分離現(xiàn)象；根據(jù)驗(yàn)證比對(duì)后推測，枯葉蛺蝶樣品中含有外源污染.

2.5 蛺蝶科昆蟲基因組比較

與本文枯葉蝶進(jìn)行比對(duì)的有33種蛺蝶的基因組組裝信息，這些蝴蝶的基因組組裝信息來自于NCBI-genome (https://www.ncbi.nlm.nih.gov/genome).檢索得到的33 條蛺蝶科昆蟲相關(guān)的全基因組信息見表2.由表中可知，尚有喙蝶亞科Libytheinae、螯蛺蝶亞科Charaxinae、閃蝶亞科Morphinae、絹蛺蝶亞科Calinaginae、釉蛺蝶亞科Heliconiinae、線蛺蝶亞科Limenitidinae、絲蛺蝶亞科Cyrestinae、苾蛺蝶亞科Biblidinae和小紫蛺蝶亞科Apaturinae等8亞科未見全基因組信息.在已報(bào)道的完成全基因組測序的8 種蛺蝶中，帕眼蝶(P.aegeria)的基因組為501.7 Mb，偏瞳蔽眼蝶(B.anynana)為475.4 Mb[24]，金斑蛺蝶(H.misippus) 為408.79 Mb[19]，黑妹袖蝶(Heliconiushimera)為399.74 Mb，Dionevanilla為390.75Mb[21]，慶網(wǎng)蛺蝶(M.cinxia)為389.908 Mb[18]，黑脈金斑蝶(D.plexippus)為241.87 Mb[20]，紅帶袖蝶(H.melpomene)為239.614 Mb[22].我們測定的枯葉蝶的基因組大小為545.05 Mb，除略小于Maniolajurtina(614.74 Mb)和翠袖鋸眼蝶(Elymniashypermnestra)(566.93 Mb)外，明顯大于上述8種蛺蝶及其余蛺蝶的基因組.

表2 枯葉蝶基因組組裝數(shù)據(jù)及與33個(gè)蛺蝶科昆蟲基因組比較Tab.2 Genome assembly statistics of K.inachus and comparisons to 33 available genomes of Nymphalidae

枯葉蛺蝶基因組的GC 含量為35.24%(圖3).枯葉蝶的GC含量低于蛺蝶科眼蝶亞科的Maniolahyperantus(37.89%)、帕眼蝶(37%)、偏瞳蔽眼蝶(36.6%)及Taenariscatops(36.6%)，高于其他29種蝴蝶的GC 含量(表2).已有研究表明GC 含量存在種的差異[25]，在蛺蝶科GC含量在30%～38%(表2). 如果GC 含量過高(>65%)或過低(<25%)，都可能導(dǎo)致Illumina 測序平臺(tái)上的序列偏倚，并嚴(yán)重影響基因組裝配[26].

3 討論

在全基因組測序中，一般要求Contig N50>20 kb，Scaffold N50>300 kb，才能拼接出比較滿意的基因組序列[17].二代測序技術(shù)中常用的組裝軟件包括 GS Assembler、ABySS、Velvet和SOAP 等.為了克服單個(gè)拼接軟件分析結(jié)果精確度的局限性，我們使用了多種序列的拼接軟件用于枯葉蛺蝶測序結(jié)果的拼接中，這樣序列拼接的精確度才能得到顯著的提高.這對(duì)于眾多的非模式昆蟲全基因組測序研究具有一定的借鑒意義.枯葉蛺蝶基因組的雜合率為1.22%，說明其基因組具有較高雜合率，其雜合率遠(yuǎn)遠(yuǎn)高于同科著名的遷徙性蝴蝶黑脈金斑蝶(0.55%).這有待于今后的進(jìn)一步研究探討.

本研究經(jīng)過對(duì)枯葉蛺蝶基因組大小的首次測定及相應(yīng)參數(shù)的初步評(píng)估，現(xiàn)得出如下結(jié)論：(1)目前枯葉蛺蝶基因組大小可以粗略估計(jì)為545 Mb；(2)通過對(duì)枯葉蛺蝶基因組28.25 G測序數(shù)據(jù)采用K-mer17進(jìn)行survey分析，預(yù)估得到基因組大小為553.55 Mb，修正后為545.05 Mb，重復(fù)序列比例為48.03%.因?yàn)榭萑~蛺蝶有較高的雜合度和重復(fù)，由各種分析指標(biāo)來看，可以推測該測序枯葉蛺蝶基因組為復(fù)雜基因組范疇，但用全基因組鳥槍法策略進(jìn)行組裝還是存在一定的風(fēng)險(xiǎn)和難度.