白邊側足海天牛全基因組Survey分析

李歆毓，丁夢瑩，馮爾輝，張源源，章可蘭，萬迎朗

（1.海南東寨港國家級自然保護區管理局，?？?，570100; 2.海南大學 海洋學院，海口，570228;3.海南大學 熱帶作物學院，?？冢?70228; 4.中國熱帶農業科學院 橡膠研究所，海口，571101;5.農業農村部橡膠樹生物學與遺傳資源利用重點實驗室，海口，571101）

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

1 材料與方法

1.1 材料供試材料白邊側足海天牛（以下簡稱為海天牛）收集于海南省?？谑袞|寨港紅樹林保護區（110°38′26″ E, 19°56′31″N），并于海南大學熱帶作物學院飼養，84K楊樹由北京林業大學細胞生物學實驗室王鑫偉提供。

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

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

1.4 測序數據的質量控制比較原始數據與過濾后數據的質量控制數據后，使用 FastQC（v 0.11.9）軟件對過濾后的數據展開質量控制，包括對數據量的概覽，并統計了讀長每個位置測序質量，總體reads測序質量趨勢，A、T、G、C堿基含量估計測序是否存在偏差，是否存在污染，數據處理時是否需要去冗余；從而實現對前期數據處理時，盡量高標準，嚴格質量控制。

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

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

2 結果與分析

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

圖1 海天牛流式細胞分析

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

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

表 1 基因組序列數據量統計

表 2 過濾數據的基本信息

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

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

表 3 白邊側足海天牛的 K-mer 數據統計

圖3 白邊側足海天牛過濾數據情況圖

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

表 4 白邊側足海天牛預組裝結果統計

3 討 論

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

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

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