西伯利亞白刺基因組信息初探

朱禮明，黎夢娟，張景波，楊秀艷,4，成鐵龍＊

(1. 南京林業大學，林木遺傳與生物技術省部共建教育部重點實驗室，江蘇 南京 210037；2. 南京林業大學南方現代林業協同創新中心，江蘇 南京 210037；3. 中國林業科學研究院沙漠林業實驗中心，內蒙古 磴口 015200；4. 中國林業科學研究院國家林業和草原局鹽堿地研究中心，北京 100091)

西伯利亞白刺（Nitraria sibirica Pall）系蒺藜科白刺屬植物，為第三紀孑遺植物，分布于蒙古、中亞以及我國西北、華北、東北的沙地、鹽堿地地區[1]。西伯利亞白刺具耐鹽堿、抗風沙等特性，能在沙漠鹽堿等惡劣環境下生存，是一種優良的沙地、鹽堿地改良物種，其果實富含多種氨基酸、糖類、黃酮等物質[2-4]，營養價值豐富，其地上部分也可作為牲畜飼料。因此，西伯利亞白刺兼有生態和經濟價值，有較好的開發利用前景。

目前，關于西伯利亞白刺的研究主要集中在繁殖技術優化[5-6]、果實成分測定[7-8]及生理生化測定[9-11]等方面，有關西伯利亞白刺的分子生物學方面的研究較少[12]，基因組學方面的研究也尚未見報道。宏觀的研究只能從表層揭示西伯利亞白刺抗逆適應現象 ，并不能從內部機制、進化等層面解釋西伯利亞白刺抗逆機理，而全基因組測序可以獲取典型基因組特征并獲得大量基因序列，對于剖析其生長、發育、抗逆等機理，發掘西伯利亞白刺的生態和經濟價值有積極意義[13-14]。

全基因組調查通過了解待測生物基因組的基本特征，可以對全基因組測序組裝難度、組裝時間和成本等作出大致的評估并作出相應的測序策略調整，是基因組測序前必不可少的步驟之一。

流式細胞術是一種快速預測基因組大小的技術，它通過比較待測植物和標定植物細胞懸液熒光吸收峰相對比值，再根據標定植物的基因組大小來計算待測植物基因組大小[15]。而隨著基因組測序技術的成熟及成本的下降，通過全基因組survey來探究待測植物的基因組基本特征不失為一種有效的方法，作為近年來發展較快的基因組預測技術，全基因組survey可以對生物的基因組基本特征測定評估[16-17]，相比于流式細胞術等基因組大小預測方法，不僅可以精準預測基因組大小，還可以對基因組復雜程度、雜合率、重復序列比例等有相應的評估，更能切合生物的基因組特征，因而有更好的參考價值。

SSR分子標記以其高重復性、高多態性、共顯性遺傳、豐度高等優良特性成為了研究群體遺傳學、遺傳變異和標記輔助選擇的有力工具，對于了解西伯利亞白刺的進化有積極的作用。

本研究基于流式細胞術和全基因組survey測序的方法對西伯利亞白刺基因組大小、復雜程度、雜合率等基因組特征有一個較為詳細的評估，同時也對其測序方案的制定提出建議，為后續西伯利亞白刺基因組組學研究奠定了良好的基礎。

1 材料和方法

1.1 材料

將取自內蒙古磴口的野生西伯利亞白刺種子置于4℃下沙藏30 d，置于萌發盒上進行萌發，再將發芽的種子定植于7 cm×7 cm的塑料花盆中(基質配方為河沙∶營養土=1∶1，并在其中摻入少量珍珠巖和蛭石)，幼苗生長2個月后取嫩葉備用。流式標定植物為 Jaroslav Dolezˇel博士惠贈的番茄‘Stupicke′poln?′ rane′’ 32 品種。

1.2 方法

1.2.1 流式細胞分析 使用BD公司influx型號流式細胞儀對西伯利亞白刺基因組大小進行分析，選用mG解離液對植物葉片進行解離，使用碘化丙啶（PI）溶液為熒光染料，采用本番茄作為內標，使用Influx自帶分析軟件FACSTM分析基因組大小。

操作步驟：于塑料皿上滴加1.5 mL mG解離液，分別取0.5 g西伯利亞白刺、番茄新鮮葉片用刀片迅速切碎后過400目濾網，將收集的濾液1 500 rpm，離心6 min，吸除上清液后重新加入500 μL預冷的mG解離液，加入PI染色液,最后加入10 μg·mL-1的 Rnase,避光 4℃ 孵育 5 min 后低速上機檢測。

C值計算公式：C待測樣本=C標定×（G0/G1待測樣本/G0/G1標定）

式中：G0/G1為流式熒光吸收強度。

mG解離液配方：

45 mmol·L-1MgCl2，20 mmol·L-1MOPS， 30 mmol·L-1Na3C6H5O7·2H2O， 1%（ w/v） PVP-40，0.2%（v/v）TritonX-100，10 mmol·L-1Na2EDTA，20 μL·mL-1β-巰基乙醇，調節 pH 至 7.0，-20℃ 下保存。PI為碘化丙啶，使用時至終濃度為50 μg·μL-1，4℃保存。

1.2.2 DNA的提取以及質量檢測 采用CTAB法對西伯利亞白刺的新鮮葉片進行DNA提取，得到的DNA樣品用紫外分光光度計檢測其濃度、OD260/OD280，再經1%瓊脂糖凝膠電泳檢測其完整性(電泳條件為：電壓180 V,電泳時間：30 min)。

1.2.3 文庫制備及測序方法 檢測合格的DNA樣品通過Covaris超聲波破碎儀打斷成片段，并進行末端修復，加poly-A尾，加測序接頭，純化，PCR擴增等步驟后，構建出350 bp雙端PE150待測序文庫。文庫通過Illumina Hiseq平臺進行雙端PE測序。

1.2.4 K-mer分析 采用K-mer分析策略，若每條序列的長度為L，K-mer長度為K，可以得到LK+1個K-mer，再通過這些數據來對基因組大小進行預估，通過Lander-waterman算法對西伯利亞白刺基因組大小進行估計，滿足公式：

式中：Nbase和NK-mer為序列的堿基總數和K-mer數，Cbase和CK-mer為覆蓋堿基的期望深度和K-mer期望覆蓋深度。

對預估的基因組大小進行修正，將K-mer深度為1的情況認為是錯誤情況，計算錯誤率，并用于修正基因組大小，修正公式為

式中：Grevised為修正后的基因組大小，E為測序錯誤率。

通過K-mer數學分析模型，基因組雜合率公式為：

式中：a1/2為雜合K-mer種類數的百分比，nK為所有K-mer的種類數。

另外，計算標準泊松分布和實際數據曲線峰值后的面積差值，可得到重復序列百分比，在這里我們計算純合峰深度1.8倍后面的K-mer個數所占的比例來估計重復序列比例。

1.2.5 基因組組裝 由于西伯利亞白刺基因組重復序列較多，我們選擇K-mer=41將打斷的DNA序列拼接組裝到Scaffold，通過reads之間的overlap關系構建de Bruiji圖并對其簡化，在重復區域邊界位置進行剪切，得到contig序列，再根據大片段數據的Pair-end關系，構建Scaffold序列，最后用reads對Scaffold的gap區域進行填補，完成組裝過程，具體配置參數為

pregraph : -K 41 -R -d 1

-K kmer: K value in kmer

-R (optional): unsolve repeats by reads (default no)

-d KmerFreqCutoff(optional): delete kmers with frequency no larger than (default 0)

contig : -D 1 -M 1 -R

-D EdgeCovCutoff(optional): delete edges with coverage no largert than (default 1)

-M mergeLevel (default 1,min 0, max 3): the strength of merging similar sequences during contiging

-R solve_repeats (optional): solve repeats by read paths(default: no)

map : -K 41

-K kmer (default: the same as in pregraph): k value in kmer

scaff : -F 1 -L 43

-F (optional) fill gaps in scaffold. (default 0;1:normally; -1:only fill nonrepeat gap; 2:radically)

-L minLen : shortest contig (minus K value) for scaffolding

再根據組裝結果統計其contig分布情況，統計測序長度大于500 bp的測序深度和GC含量并做GC含量分布圖。

1.2.6 SSR分布特征分析 運行MISA腳本(pgrc.ipk-gatersleben.de/misa)對過濾后數據SSR位點鑒定并統計其類型、數量。篩選標準為單核苷酸SSR位點≥16次,雙核苷酸SSR位點≥6次,三四核苷酸SSR位點≥5次。

2 結果與分析

2.1 流式細胞基因組大小分析

將西伯利亞白刺和番茄的葉片混合解離液放入流式細胞運行并在480 nm波長下檢測其熒光吸收強度(圖1)，其中,P0為西伯利亞白刺的吸收峰，P1為番茄的吸收峰，番茄參考2C值為1.96 pg，實驗重復3次。將平均值代入C值計算公式得出：2C西伯利亞白刺=2C番茄×（G0/G1西伯利亞白刺）/（G0/G1番茄）=1.96 pg×0.534，得西伯利亞白刺C值大小為523.4 Mbp。

圖 1 流式細胞測定結果Fig. 1 Flow cytometry results

2.2 DNA的提取以及質量檢測

取1 μL DNA樣品于分光光度計的檢測，結果顯示 OD260/OD280為 1.89，濃度為 206.9 ng·μL-1。再利用1%瓊脂糖凝膠電泳檢測其條帶完整性，圖2 表明：電泳條帶單一，無明顯雜帶。綜合二者推測，此DNA完整度較高，可用于下游實驗。

圖 2 DNA瓊脂糖凝膠電泳圖Fig. 2 DNA agarose gel electrophoresis

2.3 測序數據產出及質控

2.3.1 測序數據統計 過濾掉無效或低質量的reads數據，再經圖像識別、去污染等步驟，得出最終的測序結果（表1）：其中，測序的總reads數為212 852 294個，測序的總數據大小為63 855.69 Mbp，按照536.16 Mbp的預估基因組大小得出本次測序深度為119.09×，測序的錯誤率為0.04%，Q20的含量為95.59%，Q30的含量為89.33%，GC含量為36.78%。

表 1 測序結果統計Table 1 Sequencing results statistics

2.3.2 測序質量檢測 測序數據的質量主要分布在Q30（≥80%）以上，這樣才能保證后續分析的正常進行，如圖3所示，實驗Q30含量為89.33%滿足后續分析要求。

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

此外，測序錯誤率也影響測序結果的準確性，對于下游分析至關重要，本實驗2個reads的測序錯誤率均低于1%（圖4），表明本次測序錯誤率控制良好。為進一步保證測序結果的可信性，還需對本次測序的堿基含量分布進行分析。GC含量分布檢查用于檢測有無AT、GC分離現象，理論上G和C含量以及A和T含量在每個測序循環上應分別相等，且整個測序過程中穩定不變，呈水平線。由于DNA模板擴增偏差等原因使測序前幾個堿基測序質量值較低，發生小幅度波動，屬于正常情況。本實驗中（圖5）測序的G和C的含量和A和T的含量接近也保證了測序的可信度。

圖 4 測序錯誤率分布Fig. 4 Sequencing error rate distribution

圖 5 GC含量分布圖Fig. 5 GC content distribution map

2.4 K-mer分析

利用K-mer分析法對西伯利亞白刺基因組大小進行估計，根據測序結果（表2、圖6）發現：當K-mer深度為89×時存在明顯的主峰，由K-mer相關公式計算得到的基因組大小為536.16 Mbp，并通過后續基因修正得修正后基因組大小為526.30 Mbp; 而在主峰前橫坐標二分之一處出現次峰。一般當目標序列存在雜合現象時，存在雜合位點的K-mer被分成2份，頻率變成原頻率的1/2，因此，此峰為雜合峰，并統計得出西伯利亞白刺基因組雜合率為0.90%，雜合率較高，屬于復雜基因組。此外，在約為主峰2倍depth的地方存在次峰，并有明顯的拖帶現象，該片段出現的期望值是大部分的2倍，這些片段為重復片段，由相關統計結果得重復序列數占總序列數的55.39%。

表 2 K-mer=17分析所得各項數據Table 2 K-mer=17 analysis of the data

圖 6 K-mer=17 Depth和K-mer種類數頻率分布圖Fig. 6 K-mer=17 Depth and K-mer species frequency distribution

2.5 基因組組裝

2.5.1 數據組裝結果 運用Soapdenovo軟件拼接上述測序數據，并對數據進行糾錯，構建contig、scaffold等優化過程，得到初步的基因組組裝信息(表3)：針對組裝好的長度大于等于100 bp的scaffold內部contig進行統計，得N50長度為1 076 bp，N90為 147 bp，組裝得到最長的序列長度為45 660 bp，組裝的contig總數量為917 423個，總長度為424 458 883 bp。進一步將所有文庫測序得到的reads比對回初步得到的contigs，利用reads之間的連接關系和插入片段大小信息，過濾掉長度＜100 bp的 contig序列，最終將 contigs組裝成scaffolds，結果顯示：N50的長度的1 889 bp，N90為189 bp，最長序列長度為89 063 bp，組裝總量為717 232個，總長度為443 258 576 bp。

表 3 基因組組裝結果統計Table 3 Genomic assembly results statistics

2.5.2 GC含量分布分析 GC含量是反映植物基因組成的重要指標之一，GC含量深度分析圖用于檢測測序是否存在GC分布偏向，樣品是否存在細菌的污染等。由圖7可得：西伯利亞白刺基因組測序沒有明顯的GC偏向。圖中有2處GC聚集處，為了確認低測序深度區域是否為細菌污染造成，將低測序深度序列比對到NCBI核苷酸數據庫，并沒有細菌序列被比對上，說明樣品沒有被細菌污染，推測這是由于西伯利亞白刺基因組高雜合度所造成的。由于在組裝過程中同源染色體上雜合部位只能被識別出一半，導致此部位的GC含量分布在低測序深度區域。

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

2.6 SSR位點分析

由MISA腳本分析西伯利亞白刺基因組數據并統計（表4），共搜尋到521 125個SSR位點，其中，單核苷酸位點出現比例最高，達342 883個，占總SSR位點的65.80%；二核苷酸位點146 312個，占比28.06%；三核苷酸位點26 133個，占5.02%；四個及以上核苷酸位點8 678個，占1.67%。所以，單核苷酸重復是西伯利亞白刺主要的SSR重復位點，同時單核苷酸重復中A/T占比最多，達到了63.94%。

表 4 西伯利亞白刺SSR位點統計Table 4 SSR locus statistics of N. sibirica

3 討論

基因組大小是指生物單倍體染色體中DNA的含量，也稱為C值[18]。目前為止已有數千種動植物的C值被檢測并收錄入相應的動植物C值庫[19-20]。DNA的C值是生物體重要的基因特征，是種群分類的證據之一，也是開展各項基因工作的基礎。了解基因組大小對于推測物種的演化趨勢、進化地位、種屬間進化關系、生物進化分類等具有深遠的意義。

基因組大小預測常使用流式細胞術[21]、Feulgen圖像分析法[22]、全基因組survey調查[23]等方法。流式細胞術通過比較待測植物和內標植物細胞懸液熒光吸收峰比值，根據公式由內標植物的基因組來計算待測植物基因組的大小，是一種快速、便捷的基因組預估的方法，在測定動植物體的基因組大小方面均有較廣的應用。

全基因組survey測序是基于小片段文庫的低深度從頭測序，通過對原始數據進行圖像識別，去污染、去接頭等步驟，再進行K-mer分析，Soapdenovo軟件組裝繼而完成整個分析過程，可對基因組的大小、GC含量、雜合率以及重復序列的含量等重要的基因組特征信息進行分析，相比于流式細胞儀、Feulgen圖像分析法等基因組大小預測方法更能切合所測生物體基因組特征，是一種更精確的分析未知基因組特征的途徑[24-26]。

西伯利亞白刺基因組GC含量為36.78%，沒有明顯的過高或過低的情況[27]，對NGS測序準確性影響較小；而其雜合率為0.9%，基因組重復序列比例達55.39%，屬于高雜合基因組。推測可能是由于西伯利亞白刺在地理分布上較廣，生態條件懸殊、植物形態變化也較大有關[28]。

一般來說，基因組雜合度越大，重復片段越多，該物種的組裝難度就越大。西伯利亞白刺屬于高雜合基因組植物，而同為高雜合基因組的胡楊利用全基因組鳥槍法結合Fosmid拼裝策略獲得了精度較高的基因組圖譜[29]具有一定參考意義，如果使用二代測序Platanus組裝軟件[30]可能更適合于西伯利亞白刺基因組的拼裝。隨著近年來測序成本的下降和3代測序技術的普及，二代llumina搭配三代Pacbio輔以Hi-C技術的方案將會是西伯利亞白刺全基因組測序更好的選擇，更有利于獲得高質量的全基因組圖譜。

4 結論

本實驗測得西伯利亞白刺基因組大小為536.16 Mbp，修正后為526.30 Mbp，雜合率為0.90%，重復序列比例為55.39%；西伯利亞白刺Contig N50為1 076 bp，總長為424 458 553 bp，Scaffold N50為1 889 bp，總長為443 258 576 bp。西伯利亞白刺有521 125個SSR位點，其中單核苷酸位點有342 883個，二核苷酸位點有146 312個，三核苷酸位點有26 133個，四個及以上為8 678個，單核苷酸為其主要的SSR特征。