基于Survey分析的瀕危植物四合木基因組研究

王子寅劉秉儒*李子豪王繼飛楊鵬斌

（1.北方民族大學生物科學與工程學院，銀川 750021；2.國家民委黃河流域農牧交錯區(qū)生態(tài)保護重點實驗室，銀川 750021；3.寧夏特殊生境微生物資源開發(fā)與利用重點實驗室，銀川 750021；4.寧夏賀蘭山國家級自然保護區(qū)管理局，銀川 750021）

四合木（Tetraena mongolica）是蒺藜科（Zygophyllaceae）四合木屬（Tetraena）植物，是我國一級重點瀕危保護植物［1］，東阿拉善—西鄂爾多斯地區(qū)典型的單種屬群系，四合木有著植物界“活化石”和“大熊貓”的美稱，在內蒙古、寧夏等地區(qū)有分布且零散分布于俄羅斯、烏克蘭等國家［2］。四合木是一種強旱生植物，由于長期生長在干旱以及荒漠化的惡劣環(huán)境，使得四合木在形態(tài)結構及生理功能上產生了獨有的特征，如生活型為矮小灌木、根系十分發(fā)達、葉片肉質化且體積小以及植株被毛等［3］，在防風固沙、保持水土和維持荒漠生態(tài)系統(tǒng)功能等方面有重要價值［4］。

自20 世紀90 年代以來，由于全球氣候變化、過度砍伐、工業(yè)破壞以及害蟲侵蝕等多種自然和人為因素的綜合影響，四合木的生存環(huán)境被嚴重破壞，其生存面積減少了19.36%，棲息面積僅為2 700 km2，因此亟需開展對四合木種群的保護［5］。國內外眾多學者對于四合木的研究主要集中于群落結構及多樣性［6］、繁育技術［7］、葉綠體基因組［8-9］、病蟲害調查與保護措施［10］、化學成分分析［11-13］及其根際微生物［14-16］等方面。例如Wu 等［17］研究發(fā)現(xiàn)四合木莖葉中的提取物均顯示出一定的殺蟲潛力，促進了植物性來源的植物殺蟲劑的開發(fā)與利用。Usman等［18］通過研究四合木Pb的生物積累和抗氧化反應，發(fā)現(xiàn)四合木能夠合成關鍵的酶（如CAT、SOD 等），提高酶活性來應對高濃度的Pb，為四合木的重金屬脅迫的研究工作提供了一定的理論基礎。眾多國內外學者針對瀕危植物四合木展開了大量深入的研究，并且取得了豐富的研究成果。但由于關于四合木基因組的研究還較為缺乏，因此有必要測定四合木基因組大小，為四合木的生物遺傳和多樣性保護、種質資源利用等研究提供一定科學依據(jù)。

本研究采用第二代高通量測序技術的方法，對四合木全基因組的大小序列進行檢測分析并進行評價，全面了解其基因組特征，旨在深入挖掘四合木獨特的超旱生和繁殖的調控機理，深入挖掘四合木基因特色，進而有效地掌握保護和發(fā)揮四合木潛在巨大的基因價值，為后續(xù)四合木的分子遺傳育種奠定一定的基礎。

1 材料與方法

1.1 材料

四合木材料采自寧夏賀蘭山國家級自然保護區(qū)核心區(qū)四合木保護站，該保護站位于石嘴山市惠農區(qū)落石灘（39°22′N，106°45′E），海拔1 100～1 200 m。

1.2 實驗方法

1.2.1 樣品提取與檢測

采用改良的Cetyltrimethylammonium Bromide（CTAB）法提取的四合木DNA，在0.5%～2.0%的凝膠瓊脂糖條件下測定基因組DNA 模板完整性［19］。

1.2.2 測序

使用超聲波破碎儀（Covaris M220，美國）隨機破碎合格的四合木DNA 樣本，破碎片段長為350 bp，經過末端修復、加A 尾、加測序接頭、純化、PCR 擴增等步驟制備出整個文庫；構建完畢的文庫在Illumina NovaSeq 平臺進行高通量雙末端（Paired-End）測序［20］。為了保證高分析質量，需要對下機數(shù)據(jù)進行質控，將下機reads 的接頭序列、對后續(xù)干擾信息分析的reads、N（無法確定堿基信息）比大于10%的reads 與低質量reads 過濾，得到高質量的clean reads，為后續(xù)的基因組大小、雜合度、GC 含量等信息分析做準備［21］。該項工作由北京諾禾致源生物信息科技有限公司完成。

1.2.3 K-mer分析及基因組大小估計

通過K-mer 分析，利用過濾后的clean reads，取K=17，對四合木基因組大小、雜合率等基因組特征，進行初步分析。根據(jù)K-mer深度分布曲線獲得K-mer 深度估計值，可計算出該物種基因組大小。雜合率通過計算序列中雜合位點的比例得到。GC 圖可用于考察序列中是否存在污染或GC偏移，其分層可用于確定基因組的雜合度和重復序列情況。

2 結果與分析

2.1 測序數(shù)據(jù)統(tǒng)計及質量評估

本次測序使用高通量測序技術，得到四合木的raw reads 150.4 Gb，由于原始數(shù)據(jù)會含有一系列低質量的reads，會干擾后續(xù)的研究分析，為保證信息分析質量，過濾raw reads 中低質量數(shù)據(jù)，最終得到四合木clean reads 132.1 Gb，占原始數(shù)據(jù)的87.86%。有效測序堿基數(shù)占測序總數(shù)的99.89%。在第二代測序技術中，測序一個堿基會得到一個相對的質量值（Q），質量值是衡量測序準確度的重要指標。質量值越高，則該次測序中堿基排序錯誤的可能性越小。原始數(shù)據(jù)質量分布在Q30以上，表明測序質量較好（見圖1）。Illumina 平臺規(guī)定Q20、Q30應至少分別達到90%、85%，而此次測序中Q20、Q30分別為97.17%和92.22%（見表1），表明此次四合木基因組高通量測序也具有較高的準確性。圖2顯示了單堿基比例，用于檢測AT和CG是否存在分離情況，從中可以看出A、G、C、T 含量相近，N 含量幾乎為零。在實際操作過程中，由于DNA 模板擴增的偏差以及前幾個堿基測序質量值較低等因素，會導致每個read 前幾個堿基產生較大波動，這屬于正常情況，結果同樣表明測序質量較高。圖3顯示了測序錯誤率，錯誤率試驗可用于檢驗測序長度區(qū)域范圍內有無可能出現(xiàn)異常錯誤的堿基對序列，如在中間部位出現(xiàn)的錯誤堿基對序列錯誤率明顯地超過了其他中間部位，即有可能會出現(xiàn)異常錯誤的堿基對序列。通常對所有堿基位置的測序的失敗率都應小于1.00%。圖中所有點的堿基錯誤率都是小于0.04%，因此不存在異常堿基。

圖1 原始數(shù)據(jù)的質量分布Fig.1 Quality distribution of the original data

圖2 GC含量分布圖圖中堿基位置150 左半部分是read-1 GC 內容分布，150 右半部分是read-2 GC 內容分布；不同的顏色代表不同的堿基，用于檢測是否存在AT和GC分離Fig.2 Distribution map of GC content The left half of base position 150 in the figure is the read-1 GC content distribution，and the right half of 150 is the read-2 GC content distribution；Different colors represent different bases and are used to detect the presence of AT and GC separation

圖3 測序錯誤率分布圖中堿基位置150 左半部分為read-1 錯誤率分布，右半部分為read-2錯誤率分布Fig.3 Distribution map of sequencing error rate The left half of base position 150 in the figure is the read-1 GC error rate distribution，and the right half of 150 is the read-2 error rate distribution

表1 四合木測序數(shù)據(jù)統(tǒng)計Table 1 Data statistics of T.mongolica

2.2 污染評估及處理

若四合木基因組DNA 樣品存在污染，不僅會導致數(shù)據(jù)的有效性產生誤差，還會影響Survey 分析結果的準確性，具體影響有：基因組大小、雜合度、重復序列以及GC 含量等多項基因組指標結果均會出現(xiàn)較大誤差［22］。為了判斷提取到的四合木DNA 是否受到污染，從350 bp 文庫中，隨機取10 000條reads，通過BLAST軟件比對NCBI核苷酸數(shù)據(jù)庫（NT 庫），將比對次數(shù)排名前6 的物種展示（見表2）。其中比對結果比率最高的物種是喬木維蠟木（Bulnesiaarborea），占比為0.19%，查看其分類地位發(fā)現(xiàn)該物種與四合木同屬蒺藜科。本次比對結果中沒有發(fā)現(xiàn)動物等異常物種的高比例情況，因此判斷此次測序樣品無污染，可用于進一步的基因組圖譜分析。與進化相對較遠的物種相比若出現(xiàn)高比例讀數(shù)，則判斷DNA樣品可能被污染，應進一步檢查原因或重新提取DNA 樣品，否則不能用于后續(xù)基因組分析。

表2 原始數(shù)據(jù)文庫與NT庫比對情況Table 2 Comparison between the original data and NT data

2.3 K-mer分析及基因組大小估計

使用四合木132.1 Gb的clean reads 數(shù)據(jù)，通過K-mer方法來預估四合木基因組的大小、雜合率和重復序列等基因組特征。K-mer 個數(shù)頻率分布和K-mer 種類數(shù)頻率分布如圖4 所示，顯示K-mer 分布曲線的成峰情況良好。圖4B 中可觀察到明顯的3 個峰，第1 個峰的Kmer 深度分布為11，第2 個峰的Kmer深度分布為23，第3個峰的Kmer深度分布為45，這3個峰的Kmer深度分布符合1∶2∶4的趨勢，因此推測四合木為四倍體。其中在depth=23附近是主峰值，即泊松分布對應的期望值。表3中根據(jù)SOAP軟件預測得到K-mer總數(shù)是24 822.75 Mb，由公式（基因組大小=K-mer總數(shù)/K-mer期望深度）計算得到了四合木基因組大小為1 079.25 Mb，經過修正后的四合木基因組大小為為1 065.84 Mb。四合木基因組雜合率為0.76%，重復序列比率為75.25%。

圖4 K-mer數(shù)量分布（A）和種類分布（B）Fig.4 K-mer number distribution（A）and species distribution（B）map

表3 K-mer分析所得基因組特征Table 3 Statistical analysis of genomic characteristics by K-mer analysis

2.4 基因組初步組裝統(tǒng)計

利用K-mer=41，對132.1 Gb 的四合木clean reads 進行初步組裝，構建Contig 和Scaffold，得到了最優(yōu)的組裝結果（見表4）。由表4 可知一共得到contigs 共3 502 126 條，總計682 Mb，其N50 為187 bp，最長contig 長22 300 bp，大于2 000 bp 的contig 有8 771 條。由于本研究只對大于100 bp 的Scaffold 進行統(tǒng)計，且Contig 統(tǒng)計是針對組裝長度大于等于100 bp 的Scaffold 內部的Contig 進行統(tǒng)計，故二者統(tǒng)計數(shù)據(jù)相等。圖5～6 顯示了顯著峰值，確定36×左右峰值為純合子峰，位于純合子峰橫坐標一半左右的峰為雜合子峰。故可以判斷四合木基因組是復雜的雜合子峰。

圖5 contig覆蓋深度和長度分布Fig.5 Distribution map of the depth and length of contig coverage

圖6 contig覆蓋深度和數(shù)量分布Fig.6 Distribution map of the depth and quantity of contig coverage

表4 四合木基因組組裝序列統(tǒng)計Table 4 Genomic assembly sequence statistics of T.mongolica

2.5 GC含量和分布情況

以10 kb 為窗口對四合木基因組的GC 堿基含量及測序深度統(tǒng)計，結果如圖7所示，橫坐標為GC含量，縱坐標為測序深度，右方表示測序深度分布，上方表示GC 含量分布，圖中紅色部分表示散點分布密度大。戚行江等［23］研究表明GC 含量過高（＞65%）、GC 含量過低（＜25%）均會導致基因組測序的偏向性，最終影響基因組的組裝效果。本研究結果表明幾乎所有窗口的GC 含量都在20%～60 %，通過對圖7 綜合分析，四合木基因組GC 含量約為33.57 %，結果較為適中，因此不會影響分析的準確度。其測序深度在15×～75×范圍內集中分布，說明測序過程中無污染。從GC_depth 分布圖中可以看出分為2層，結合雜合模擬實驗顯示四合木高雜合度可能是由于雜合引起?；蚪MGC含量對基因組進行測序的隨機性因素具有較大影響。

圖7 GC含量與測序深度（Depth）關聯(lián)分析Fig.7 Statistical chart of correlation analysis between GC content and sequencing Depth

3 討論與結論

基因組的雜合率和重復序列是決定基因組組裝質量的關鍵，評估測序數(shù)據(jù)的雜合率和重復序列，有助于找到合適的組裝策略［24］。研究發(fā)現(xiàn)細胞大小［25］、生物體壽命［26］、物種瀕危程度［27］等均與基因組的大小有一定關系。基于K-mer 分析的生物信息學方法已經成功應用于一串紅（Salvia splendens）［28］、西番蓮（Passifloracaerulea）［29］、樟樹（Cinnamomumcamphora）［30］等植物的基因組預測。植物基因組大小的檢測方法多種多樣，如孚爾根光密度測量法、脈沖場凝膠電泳、實時熒光定量PCR技術以及流式細胞儀等［29］。目前已有桫欏（Alsophilaspinulosa）［31］、泡桐（Paulowniafortunei）［32］、鳶尾（Iristectorum）［33］等多種植物通過流式細胞儀方法測得基因組大小。本研究對四合木基因組大小基于K-mer分析的預測結果為1 065.84 Mb，后續(xù)可采用常用的流式細胞儀測定法對其基因組大小進行測定，以檢驗K-mer分析結果的準確性。

不同植物的基因組大小有很大差異。在已知的被子植物基因組中，基因組最小的植物為螺旋貍藻（Genliseamaragretae）［34］，其基因組大小為63 Mb，而百合科（Liliaceae）的四倍體貝母（Fritillaria assyriaca）基因組［25］最大則有127 Gb，二者相差約2 050 倍。目前報道了多種旱生植物的基因組大小，如新疆沙冬青（Ammopiptanthusmongolicus）［35］的基因組大小約為787 Mb，鹽節(jié)木（Halocnemum strobilaceum）［36］的基因組大小約為663.25 Mb，好好芭（Simmondsiachinensis）［37］的基因組大小約為774 Mb等。本研究中預測得到的四合木基因組大小約為1 065.84 Mb，大于上述3 種旱生植物，推測可能是在四合木漫長的進化過程中其基因組朝著增大的方向進化，故四合木的基因組大小較高于一般常見的旱生植物。通過K-mer 分析得到四合木的雜合度為0.76%、重復序列比例為75.25%，從四合木上述基因組特征來看，四合木屬于高重復、高雜合的復雜基因組，這可能由于四合木種群位于東阿拉善—西鄂爾多斯狹長分布區(qū)中間，且為單種屬群系，與其他種群密切接觸，導致更多的基因介入，造成四合木雜合度和重復序列比例偏高，與Cheng等［38］研究結果一致。此外，由于二代測序中重復序列相似性較高、拷貝數(shù)量范圍較大，也會造成四合木重復序列比例偏高的現(xiàn)象［39］。

本研究對珍惜孑遺植物四合木進行基因組調研測序，預測基因組大小約為1 065.84 Mb，重復序列比例為75.25%；分析得出其基因組雜合度為0.76%，GC含量為33.57%，初步探索了四合木的基因組特征。隨著測序技術的發(fā)展，第二代測序技術讀長較短、組裝結果的連續(xù)性無法保證，以及由于基因組的高雜合導致過度組裝或組裝不徹底等缺點也日益突出。對于四合木高雜合、重復系列比例高的特點［39］，后續(xù)可采用第三代高通量測序技術（單分子測序）同時結合染色質區(qū)域捕獲技術進一步完善四合木組裝結果，提高重復序列的區(qū)分度，使其組裝結果達到染色質水平，二者結合不僅能夠確保準確性，而且周期短，性價比更高［40］，有望最終獲得高質量的四合木全基因圖譜，為研究四合木的遺傳進化關系提供相應的依據(jù)。