李歆毓，丁夢(mèng)瑩，馮爾輝，張?jiān)丛?，章可蘭，萬迎朗

（1.海南東寨港國(guó)家級(jí)自然保護(hù)區(qū)管理局，海口，570100; 2.海南大學(xué) 海洋學(xué)院，?？?，570228;3.海南大學(xué) 熱帶作物學(xué)院，?？?，570228; 4.中國(guó)熱帶農(nóng)業(yè)科學(xué)院 橡膠研究所，?？?，571101;5.農(nóng)業(yè)農(nóng)村部橡膠樹生物學(xué)與遺傳資源利用重點(diǎn)實(shí)驗(yàn)室，?？?，571101）

白邊側(cè)足海天牛（Elysia leucolegnote）屬于軟體動(dòng)物門囊舌目（Sacoglossa）海天牛超科（Elysioidea）海天牛屬（Elysia），該種分布于中國(guó)的香港[1]、?？赱2]以及泰國(guó)[3]、菲律賓[4]等地紅樹林。目前，在我國(guó)已報(bào)道有8種海天牛的分布記錄[1，5]。海天牛屬部分種具有保留藻類葉綠體并使其發(fā)揮光合作用的能力[6-7]。這種被吸收后且繼續(xù)發(fā)揮功能的葉綠體被稱為盜質(zhì)體[8]。不同物種的盜質(zhì)體壽命的決定因素和盜質(zhì)體發(fā)揮功能的分子機(jī)制等相關(guān)的研究[9-13]對(duì)葉綠體的內(nèi)共生理論的揭示[14]以及指導(dǎo)植物抗逆方案的設(shè)計(jì)都具有重要的意義[15-17]，但受限于光合軟體動(dòng)物材料的稀缺性，此相關(guān)研究的規(guī)模十分有限。白邊側(cè)足海天牛在我國(guó)的分布有利于科研人員就地取材，為我國(guó)在此方向的深入研究提供了獨(dú)特的材料。分子機(jī)制的揭示離不開物種基因組的信息。海天牛屬目前尚未有染色體水平上的高質(zhì)量基因組的報(bào)道?；蚪M已經(jīng)成為深入研究分子機(jī)制的關(guān)鍵內(nèi)容?；蚪M含量在同一個(gè)物種里是保持穩(wěn)定的[18-20]。高通量測(cè)序?yàn)橐环N廣泛用于多種動(dòng)植物基因組測(cè)序的技術(shù)[21]，將高通量獲得的原始數(shù)據(jù)通過質(zhì)控以后，利用K-mer分析法評(píng)估基因組大小、雜合度和重復(fù)率等特征[22-23]。基因組測(cè)序大大促進(jìn)了動(dòng)植物的遺傳進(jìn)化及功能基因研究，但深度測(cè)序之前的低覆蓋度的全基因組調(diào)查尤為重要，因?yàn)槠淇梢詻Q定基因組測(cè)序中的最合適的測(cè)序、拼裝方式[24]。本實(shí)驗(yàn)旨在采用高通量測(cè)序技術(shù)結(jié)合K-mer分析[25]，對(duì)白邊側(cè)足海天牛基因組進(jìn)行測(cè)定及評(píng)估，為后續(xù)深度測(cè)序提供基礎(chǔ)信息。

1 材料與方法

1.1 材料供試材料白邊側(cè)足海天牛（以下簡(jiǎn)稱為海天牛）收集于海南省?？谑袞|寨港紅樹林保護(hù)區(qū)（110°38′26″ E, 19°56′31″N），并于海南大學(xué)熱帶作物學(xué)院飼養(yǎng)，84K楊樹由北京林業(yè)大學(xué)細(xì)胞生物學(xué)實(shí)驗(yàn)室王鑫偉提供。

1.2 流式細(xì)胞分析取海天牛和 84K 楊樹幼嫩新鮮葉片 1 g，置于盛有 500 μL細(xì)胞裂解液Galbraith中，并迅速將其切碎。細(xì)胞核通過30 μm濾膜后，與 1 mL 濃度為 10 μg·mL-1的 PI（碘化丙啶）染色液混合，室溫孵育30 s后上機(jī)檢測(cè)。使用CyFlow?Cube8（希森美康,日本）流式細(xì)胞儀對(duì)海天?；蚪M大小進(jìn)行評(píng)估，變異系數(shù)控制在5%以內(nèi)。每個(gè)樣品重復(fù)3次。

1.3 文 庫 構(gòu) 建 及 基 因 組 測(cè) 序粉 碎 合 格 的DNA樣品為350 bp左右的目的片段，構(gòu)建文庫，經(jīng)過末端修復(fù)、加A、加接頭、目標(biāo)片段選擇和PCR等步驟，用安捷倫2 100 和定量PCR的方法檢測(cè)文庫片段大小和文庫定量，以確定文庫是否符合測(cè)序標(biāo)準(zhǔn)，通過橋式PCR的方法將文庫固定到測(cè)序芯片上；將這些兩端的片段在 Illumina Hiseq Xten（Illumina, U.S）測(cè)序儀上進(jìn)行雙末端（Paired-End）測(cè)序，獲得全基因組測(cè)序數(shù)據(jù)，去除未成對(duì)匹配的讀長(zhǎng)（reads）、低質(zhì)量讀長(zhǎng)、接頭遭受污染以及過濾掉接頭重復(fù)（duplication）的讀長(zhǎng)等對(duì)測(cè)序所產(chǎn)生的數(shù)據(jù)進(jìn)行預(yù)處理。

1.4 測(cè)序數(shù)據(jù)的質(zhì)量控制比較原始數(shù)據(jù)與過濾后數(shù)據(jù)的質(zhì)量控制數(shù)據(jù)后，使用 FastQC（v 0.11.9）軟件對(duì)過濾后的數(shù)據(jù)展開質(zhì)量控制，包括對(duì)數(shù)據(jù)量的概覽，并統(tǒng)計(jì)了讀長(zhǎng)每個(gè)位置測(cè)序質(zhì)量，總體reads測(cè)序質(zhì)量趨勢(shì)，A、T、G、C堿基含量估計(jì)測(cè)序是否存在偏差，是否存在污染，數(shù)據(jù)處理時(shí)是否需要去冗余；從而實(shí)現(xiàn)對(duì)前期數(shù)據(jù)處理時(shí)，盡量高標(biāo)準(zhǔn)，嚴(yán)格質(zhì)量控制。

1.5 K-mer 分 析 以 及 基 因 組 特 征 估 計(jì)通 過jellyfish-2（v10.7.7）軟件對(duì)序列文件進(jìn)行 K-mer 的計(jì)數(shù)和統(tǒng)計(jì)；隨后，利用負(fù)二項(xiàng)式模型（Negative binomial model）對(duì)應(yīng)的軟件 GenomeScope（v1.0）對(duì)基因組大小及其雜合度的評(píng)估，并生成最終基因組評(píng)估結(jié)果。選用K-mer值17、19兩種條件對(duì)評(píng)估結(jié)果進(jìn)行比較。

1.6 基 因 組 初 步 組 裝使 用 SOAPdenovo2（v2.03）軟件對(duì)過濾后的數(shù)據(jù)進(jìn)行拼接，拼接出Contigs序列，然后組裝基因組。SOAPdenovo2的K-mer參數(shù)設(shè)置29，其他參數(shù)選擇默認(rèn)值。再將Contigs序列構(gòu)圖形成Scaffolds序列，并利用不同插入片段估計(jì)出 Contigs間的距離，用N 堿基填起來。最后，再利用測(cè)序的雙末端數(shù)據(jù)之間的配對(duì)關(guān)系（Paired-End）以及短Reads數(shù)據(jù)對(duì)已組裝的Contigs的覆蓋信息，對(duì) Contigs間空隙（“N”）進(jìn)行局部組裝，補(bǔ)充Contigs信息，適當(dāng)延長(zhǎng) Contigs序列。有效數(shù)據(jù)與原始序列進(jìn)行對(duì)比后獲得堿基深度，在序列上以5 kb為窗口，無重復(fù)前進(jìn)，從而得到GC depth點(diǎn)圖，對(duì)組裝后的基因組進(jìn)行評(píng)估。

2 結(jié)果與分析

2.1 流式細(xì)胞結(jié)果預(yù)測(cè)基因組大小基于流式細(xì)胞術(shù)分析海天?；蚪M大小，當(dāng)變異系數(shù)控制在 5% 以內(nèi)時(shí)，以84K楊作為對(duì)照樣品（圖1-A）信號(hào)峰清晰集中，84K楊與海天牛的混合樣品的信號(hào)峰獨(dú)立分離且距離較近（圖1-B）。利用84K楊作為對(duì)照樣品，根據(jù)混合樣品PI 熒光強(qiáng)度以及峰值的倍數(shù)關(guān)系，計(jì)算海天?；蚪M是84K楊的1.69倍，84K楊的核DNA相對(duì)含量為1.129 20，基因組平均值為 470.155 Mb；估算出海天牛的核DNA相對(duì)含量為 2.218 71，基因組平均值為794.562 Mb。

圖1 海天牛流式細(xì)胞分析

2.2 建庫信息及數(shù)據(jù)量統(tǒng)計(jì)基因組調(diào)查利用第二代高通量測(cè)序技術(shù)進(jìn)行雙末端測(cè)序，獲得全基因組的序列結(jié)果。測(cè)序共得到海天牛原始數(shù)據(jù)約 25.8 Gb，共 171 847 064 條原始序列；過濾后約25.1 Gb，原始數(shù)據(jù) Q30比例為 91.33%，過濾后Q30 比例為91.78%，滿足基因組調(diào)查需要的測(cè)序數(shù)據(jù)量（表1）。比較原始數(shù)據(jù)與過濾數(shù)據(jù)（表2）的堿基的分布情況（圖2-A、 2-B），過濾前后除了測(cè)序時(shí)前幾個(gè)bp堿基含量略有波動(dòng)屬正常現(xiàn)象，其余每個(gè)測(cè)序位置A堿基和T堿比例相等，G堿基和C堿基比例相等，N堿基的數(shù)量為0。原始數(shù)據(jù)與過濾后數(shù)據(jù)的測(cè)序質(zhì)量分布在Q30到Q40之間，Q30序列占比高，表明測(cè)序結(jié)果質(zhì)量高可用于后續(xù)分析。

圖2 白邊側(cè)足海天牛過濾前后測(cè)序情況

表 1 基因組序列數(shù)據(jù)量統(tǒng)計(jì)

表 2 過濾數(shù)據(jù)的基本信息

過濾數(shù)據(jù)所有讀長(zhǎng)上的堿基質(zhì)量值大于30且波動(dòng)小，說明過濾后數(shù)據(jù)質(zhì)量穩(wěn)定（圖3-A）；實(shí)際G、C堿基含量與理論G、C堿基含量都在38%左右出現(xiàn)峰值，且沒有明顯的偏差，表明測(cè)序結(jié)果沒有偏向性（圖3-B）；過濾后所有的數(shù)據(jù)讀長(zhǎng)都為150 bp（圖3-C）；全部序列達(dá)到 Q20，超過 95% 序列達(dá)到Q30，且集中在Q36（表2）。以上結(jié)果表明，過濾后的數(shù)據(jù)讀長(zhǎng)長(zhǎng)，質(zhì)量高，沒有堿基偏好性適用于K-mer分析。

2.3 K-mer分析以及基因組大小、雜合率的估計(jì)使用K-mer的分析方法可以預(yù)測(cè)白邊側(cè)足海天牛的基因組特征。選擇K-mer的條件為17和19展開分析，樣本17-mer和19-mer分布曲線為非正常泊松分布，呈現(xiàn)雙峰分布，在17×和27×附近各有1個(gè)峰值（圖3-A、圖3-B）。總測(cè)序深度約為30×，根據(jù)17-mer分析，預(yù)測(cè)海天?；蚪M大小約為724.8 Mb，基因組重復(fù)率為52.8%，雜合度為1.55%，模型擬合值為99.38%；19-mer分析預(yù)測(cè)海天牛基因組大小約為730.8 Mb，基因組重復(fù)率為35.1%，雜合度為1.68%，模型擬合值為99.72%（表3）。

表 3 白邊側(cè)足海天牛的 K-mer 數(shù)據(jù)統(tǒng)計(jì)

圖3 白邊側(cè)足海天牛過濾數(shù)據(jù)情況圖

2.4 白 邊 側(cè) 足 海 天 牛 基 因 組 的 預(yù) 組 裝選 用SOAPdenovo2軟件對(duì)海天牛樣本進(jìn)行預(yù)組裝，設(shè)置K-mer參數(shù)為29時(shí)，在scaffold尺度上，得到含 N 堿基的基因組大小 628 574 653 bp，不含 N 的基因組大小 627 289 254 bp；Scaffold N50 長(zhǎng)度為373 bp，共 405 072 條；Scaffold 數(shù)量 2 258 693 條，最長(zhǎng)的 scaffold長(zhǎng)度為 22 424 bp。在 contig的尺度上，以 contig N50 為 358 bp 數(shù)量有 419 361 條。得到含 N 堿基的基因組大小 624 854 764 bp，不含N 的基因組大小 624 854 764 bp，最長(zhǎng)的 contig 為22 424 bp（表4）。組裝成 Scaffold 的 contig 的數(shù)量為168 878條，每個(gè)scaffold的平均contig數(shù)目為1.5。除此，還得到scaffold尺度上的各堿基的含量，堿基 A 數(shù)量為 207 638 986 bp，占總的堿基數(shù)目的 33.03%；堿基 C 數(shù)量為 110 706 106 bp，占總的堿基數(shù)目的17.61%；堿基G數(shù)量為109 134 885 bp，占總的堿基數(shù)目的17.36%；堿基T與堿基A的數(shù)量及占比基本相同，堿基T數(shù)量為199 809 277 bp占總的堿基數(shù)目的31.79%；剩下的所有為堿基 N，數(shù)量為 1 285 399 bp，占總的堿基數(shù)目的0.20%。最后計(jì)算得到G、C堿基含量為35.05%。GC-depth分析顯示，測(cè)序無偏向性；平均深度集中在30×，GC 深度分布被分為2層。

表 4 白邊側(cè)足海天牛預(yù)組裝結(jié)果統(tǒng)計(jì)

3 討 論

目前，在我國(guó)已報(bào)道有8種海天牛的分布記錄[1，5]，其中部分海天牛具有利用藻類葉綠體進(jìn)行光合作用的能力。盜質(zhì)體壽命是不等的，有的盜質(zhì)體能維持長(zhǎng)達(dá)9個(gè)月[26]，而有的只能維持短短的幾個(gè)小時(shí)。根據(jù)葉綠體在海天牛中停留的時(shí)間將海天牛分為三類，第一類為長(zhǎng)時(shí)間保存葉綠體物種（long-term retention (LtR) slugs），已報(bào)道的包括E.chlorotica, E.timida, E.crispata, E.clarki, E.viridis, Plakobranchus ocellatusandCostasiella ocellifera[27-30]；第二類是短時(shí)間保存葉綠體物種（short-term retention species, StR），其對(duì)葉綠體的保留時(shí)間不超過兩周；第三類是不保存葉綠體物種（non-retention species, NR），在食用藻類后迅速分解葉綠體[31]。在本研究中發(fā)現(xiàn)，白邊側(cè)足海天牛至少能保持盜質(zhì)體活性2個(gè)月以上，屬于能長(zhǎng)時(shí)間保存葉綠體的海天牛。

對(duì)盜質(zhì)體活性的長(zhǎng)期保持依賴于宿主核基因編碼的功能基因與盜質(zhì)體基因的協(xié)調(diào)表達(dá)。例如E.chlorotica與E.timida食用藻類的葉綠體基因組中存在一種特定基因（ftsH，一種對(duì)光系統(tǒng)II修復(fù)至關(guān)重要的D1質(zhì)控蛋白酶）其中M41金屬蛋白酶結(jié)構(gòu)域是維持盜質(zhì)體長(zhǎng)期活動(dòng)的關(guān)鍵[32-33]。同時(shí)動(dòng)物內(nèi)源的脂肪酸合酶-（FAS）樣聚酮合酶（PKS）蛋白也可以提供光保護(hù)能力，盜質(zhì)體固定二氧化碳，固定碳被轉(zhuǎn)化為甲基丙二酰輔酶a，并被軟體動(dòng)物EcPKS1酶修飾，合成紫外線-氧化阻斷吡喃，保護(hù)軟體動(dòng)物及其葉綠體免受光合損傷[34]。另外，有觀點(diǎn)認(rèn)為吞食葉綠體后從植物中攝取的酶的豐度可能限制盜質(zhì)體發(fā)揮功能。也有一種觀念認(rèn)為，盜質(zhì)體壽命的維持是通過從藻類細(xì)胞核到動(dòng)物細(xì)胞核的廣泛水平基因轉(zhuǎn)移（HGT）來實(shí)現(xiàn)。但是對(duì)于該假說還存在很大的爭(zhēng)議，早期研究中，TORRES等證實(shí)了核編碼的基因在質(zhì)體核糖體抑制劑存在的條件下可以合成LHCⅠ，并提出可以通過病毒或逆轉(zhuǎn)錄病毒實(shí)現(xiàn)HTG的假說[35]。這些假說的討論，都必須基于對(duì)海天牛核基因組及其攝取的盜質(zhì)體基因組的分析研究。

E.chlorotica的全基因組測(cè)序與組裝是以二代為主，三代PacBio輔助的scafford的基因組組裝水平, 其全基因組大小為 557 Mb，scaffold N50 為442 kb，BUSCO 注釋率為 93.3%[36]。海天?？颇壳斑€沒有染色體水平的基因組組裝結(jié)果，所以對(duì)海天?；蚪M的檢測(cè)仍舊是一個(gè)具有新穎性的課題?；蚪M調(diào)查，也稱作Survey，基于深度達(dá)到20～30×以上的高質(zhì)量的二代測(cè)序數(shù)據(jù)對(duì)物種的基因組大小與特征進(jìn)行解讀，可以為物種基因組測(cè)序方案提供重要的指導(dǎo)[37-38]。在基因組調(diào)查的基礎(chǔ)上，結(jié)合流式細(xì)胞儀可提升基因組大小預(yù)測(cè)結(jié)果的精準(zhǔn)性[39-40]。本研究中，利用 84K 楊為對(duì)照，通過流式細(xì)胞術(shù)預(yù)測(cè)海天牛基因組大小均值為794.562 Mb，K-mer分析的結(jié)果顯示其基因組大小為724～730 Mb，兩者的結(jié)果偏差較小，不影響對(duì)基因組測(cè)序決策的判斷[41]。所有的結(jié)果顯示，白邊側(cè)足海天牛是一個(gè)高度雜合的物種，且基因組大小超過700 Mb。為了達(dá)到染色體級(jí)別的組裝水平，全基因組測(cè)序建議使用以三代測(cè)序技術(shù)為主，Hi-C 或 Hi-Fi技術(shù)相結(jié)合的測(cè)序手段[42-44]，測(cè)序量達(dá)到80×～100×的深度足夠完成海天?；蚪M的精細(xì)組裝。