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

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

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

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

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

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

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

SSR分子標(biāo)記以其高重復(fù)性、高多態(tài)性、共顯性遺傳、豐度高等優(yōu)良特性成為了研究群體遺傳學(xué)、遺傳變異和標(biāo)記輔助選擇的有力工具，對于了解西伯利亞白刺的進(jìn)化有積極的作用。

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

1 材料和方法

1.1 材料

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

1.2 方法

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

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

C值計(jì)算公式：C待測樣本=C標(biāo)定×（G0/G1待測樣本/G0/G1標(biāo)定）

式中：G0/G1為流式熒光吸收強(qiáng)度。

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β-巰基乙醇，調(diào)節(jié) pH 至 7.0，-20℃ 下保存。PI為碘化丙啶，使用時(shí)至終濃度為50 μg·μL-1，4℃保存。

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

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

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

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

對預(yù)估的基因組大小進(jìn)行修正，將K-mer深度為1的情況認(rèn)為是錯(cuò)誤情況，計(jì)算錯(cuò)誤率，并用于修正基因組大小，修正公式為

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

通過K-mer數(shù)學(xué)分析模型，基因組雜合率公式為：

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

另外，計(jì)算標(biāo)準(zhǔn)泊松分布和實(shí)際數(shù)據(jù)曲線峰值后的面積差值，可得到重復(fù)序列百分比，在這里我們計(jì)算純合峰深度1.8倍后面的K-mer個(gè)數(shù)所占的比例來估計(jì)重復(fù)序列比例。

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

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

再根據(jù)組裝結(jié)果統(tǒng)計(jì)其contig分布情況，統(tǒng)計(jì)測序長度大于500 bp的測序深度和GC含量并做GC含量分布圖。

1.2.6 SSR分布特征分析 運(yùn)行MISA腳本(pgrc.ipk-gatersleben.de/misa)對過濾后數(shù)據(jù)SSR位點(diǎn)鑒定并統(tǒng)計(jì)其類型、數(shù)量。篩選標(biāo)準(zhǔn)為單核苷酸SSR位點(diǎn)≥16次,雙核苷酸SSR位點(diǎn)≥6次,三四核苷酸SSR位點(diǎn)≥5次。

2 結(jié)果與分析

2.1 流式細(xì)胞基因組大小分析

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

圖 1 流式細(xì)胞測定結(jié)果Fig. 1 Flow cytometry results

2.2 DNA的提取以及質(zhì)量檢測

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

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

2.3 測序數(shù)據(jù)產(chǎn)出及質(zhì)控

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

表 1 測序結(jié)果統(tǒng)計(jì)Table 1 Sequencing results statistics

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

圖 3 數(shù)據(jù)質(zhì)量分布Fig. 3 Data quality distribution

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

圖 4 測序錯(cuò)誤率分布Fig. 4 Sequencing error rate distribution

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

2.4 K-mer分析

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

表 2 K-mer=17分析所得各項(xiàng)數(shù)據(jù)Table 2 K-mer=17 analysis of the data

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

2.5 基因組組裝

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

表 3 基因組組裝結(jié)果統(tǒng)計(jì)Table 3 Genomic assembly results statistics

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

圖 7 GC含量與測序深度關(guān)聯(lián)分析統(tǒng)計(jì)圖Fig. 7 GC content and sequencing depth correlation analysis

2.6 SSR位點(diǎn)分析

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

表 4 西伯利亞白刺SSR位點(diǎn)統(tǒng)計(jì)Table 4 SSR locus statistics of N. sibirica

3 討論

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

基因組大小預(yù)測常使用流式細(xì)胞術(shù)[21]、Feulgen圖像分析法[22]、全基因組survey調(diào)查[23]等方法。流式細(xì)胞術(shù)通過比較待測植物和內(nèi)標(biāo)植物細(xì)胞懸液熒光吸收峰比值，根據(jù)公式由內(nèi)標(biāo)植物的基因組來計(jì)算待測植物基因組的大小，是一種快速、便捷的基因組預(yù)估的方法，在測定動(dòng)植物體的基因組大小方面均有較廣的應(yīng)用。

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

西伯利亞白刺基因組GC含量為36.78%，沒有明顯的過高或過低的情況[27]，對NGS測序準(zhǔn)確性影響較?。欢潆s合率為0.9%，基因組重復(fù)序列比例達(dá)55.39%，屬于高雜合基因組。推測可能是由于西伯利亞白刺在地理分布上較廣，生態(tài)條件懸殊、植物形態(tài)變化也較大有關(guān)[28]。

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

4 結(jié)論

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