基于流式細(xì)胞術(shù)和基因組survey分析的地黃基因組研究

趙 樂，朱畇昊，王 敏，韓永光，馬利剛，馮衛(wèi)生，鄭曉珂*

1.河南中醫(yī)藥大學(xué)藥學(xué)院，河南 鄭州 450046

2.呼吸疾病診療與新藥研發(fā)河南省協(xié)同創(chuàng)新中心，河南 鄭州 450046

3.北京工商大學(xué) 化學(xué)與材料工程學(xué)院，北京 100048

地黃Rehmannia g lutinosaLibosch 是玄參科地黃屬多年生草本植物，為《中國藥典》2015年版藥材地黃的唯一來源植物[1]，是我國傳統(tǒng)的大宗中藥材，具有很高的藥用價值和經(jīng)濟價值。地黃以其塊根入藥，塊根中含有環(huán)烯醚萜、黃酮、木脂素及多糖等多種化合物[2]。地黃藥理活性顯著，對心腦血管系統(tǒng)、神經(jīng)系統(tǒng)、免疫系統(tǒng)均有一定保護作用[3]。梓醇是地黃的主要活性成分，從化學(xué)結(jié)構(gòu)分析，梓醇屬于環(huán)烯醚萜苷類化合物，具有降血糖、神經(jīng)保護、抗炎、抗氧化等多種生物學(xué)活性[4]。Yan 等[5]發(fā)現(xiàn)梓醇能夠緩解2 型糖尿病中的胰島素抵抗，起到降血糖的作用，還能減緩糖尿病晚期的并發(fā)癥。此外，梓醇還可以改善小鼠阿爾茨海默癥的病理特征，同時改善相關(guān)的行為障礙[6-7]。

目前關(guān)于地黃的研究，主要集中于化學(xué)成分的分離、梓醇和地黃其他活性成分的藥理作用、以及地黃連作障礙等方面，關(guān)于地黃基因組方面的研究報道[8]較少。Zeng等[9]對地黃屬6種地黃植物的葉綠體基因組進行了測序，同時進行了系統(tǒng)發(fā)育關(guān)系分析，發(fā)現(xiàn)地黃葉綠體基因組大小為153 622 bp，而且地黃屬6 種地黃的葉綠體基因組在基因結(jié)構(gòu)方面高度一致。閆坤等[10]通過觀察地黃屬6 種地黃的小孢子母細(xì)胞減數(shù)分裂，發(fā)現(xiàn)地黃和茄葉地黃屬于四倍體植物（2n＝56），天目地黃、高地黃、裂葉地黃和湖北地黃屬于二倍體植物（2n＝28），說明地黃屬植物存在多倍體現(xiàn)象。這種多倍體現(xiàn)象可能是由于全基因組復(fù)制形成的，王姣等[11]利用系統(tǒng)發(fā)育生物學(xué)方法，通過分析地黃屬5 種地黃的轉(zhuǎn)錄組數(shù)據(jù)，發(fā)現(xiàn)地黃屬植物經(jīng)歷了2 次全基因組復(fù)制，加速了地黃屬植物的物種分化，通過二倍體祖先的基因組加倍形成了四倍體植物。閆婧[12]用二代高通量測序?qū)Φ攸S屬的天目地黃基因組進行了survey 分析，發(fā)現(xiàn)天目地黃基因組大小估算為1.14 Gb，GC 含量為40.64%，重復(fù)序列比例為 76.94%，雜合度為0.34%，從基因組基本結(jié)構(gòu)特征分析，天目地黃基因組屬于高重復(fù)、低雜合的大型基因組。

地黃的許多栽培種在長期的栽培過程中受到了選擇，這些栽培品種之間的遺傳關(guān)系不清楚，而且這些栽培品種在葉和花的顏色和形態(tài)、以及根的形狀方面變異較大[13]。由于地黃栽培品種之間的遺傳關(guān)系不清楚，因此獲得地黃的全基因組精細(xì)圖譜對于研究其主要藥用成分—梓醇的生物合成及代謝調(diào)控途徑具有重要的意義，同時也為研究地黃株高、葉形、塊根膨大等性狀變異奠定了基礎(chǔ)。根據(jù)報道地黃為四倍體（2n＝56）[10,13]，本研究以地黃栽培品種“沁懷1 號”為材料，通過二代高通量測序技術(shù)對地黃基因組進行survey 分析，同時結(jié)合流式細(xì)胞技術(shù)來評估地黃基因組大小，最終獲得地黃全基因組重復(fù)序列比例、基因組雜合度以及GC 含量等信息。這些信息將為今后繪制地黃全基因組的精細(xì)圖譜提供依據(jù)，也為闡明梓醇的生物合成及代謝調(diào)控途徑奠定基礎(chǔ)。

1 材料

地黃所用材料為地黃栽培品種“沁懷1 號”組培苗，經(jīng)河南中醫(yī)藥大學(xué)董誠明教授鑒定為地黃Rehmannia glutinosaLibosch，以已知基因組大小的大豆組培苗（中國農(nóng)業(yè)科學(xué)院）和辣椒組培苗（中國農(nóng)業(yè)科學(xué)院）作為流式細(xì)胞技術(shù)測定的對照樣品，采集這3 種植物組培苗新鮮的葉片，經(jīng)液氮速凍后，放置于-80 ℃超低溫冰箱保存。

BD FACSCalibur 流式細(xì)胞儀（美國BD 公司），離心機5810R（德國Eppendorf 公司）；碘化丙啶（propidium iodide，PI）批號340242，購自美國BD公司；吐溫20（批號P2287）購自美國Sigma-Aldrich公司。Otto I 緩沖液：0.1 mol/L 檸檬酸、0.5%聚山梨酯20，用0.22 μm 濾器濾過后4 ℃保存。Otto II緩沖液：0.4 mol/L Na2HPO4·12H2O，用0.22 μm 濾器濾過后室溫（18～25 ℃）保存。

2 方法

2.1 細(xì)胞核懸液制備及基因組大小測定

根據(jù)Dole?el 等[14]的方法，各取地黃、大豆和辣椒組培苗新鮮葉片100 mg，放入培養(yǎng)皿中，加入1 mL 預(yù)冷的Otto I 緩沖液，用鋒利的刀片將組織切碎，用移液器上下吹打幾次混合均勻（注意避免氣泡），然后將勻漿液通過42 μm 尼龍篩網(wǎng)濾過到離心管中，4 ℃、1000 r/min離心5 min，小心的棄去上清，收集細(xì)胞核沉淀。再用1 mL Otto II 緩沖液重懸細(xì)胞，然后加入DNA 熒光染料碘化丙啶（PI）和 RNA 酶，使二者的終濃度均在 50 μmol/mL，置于冰上，避光染色20 min，再將地黃、大豆和辣椒的染液混合形成混合樣品，然后上機測定這3 種植物的單獨樣品和混合樣品。

將已知基因組大小的大豆和辣椒作為對照樣品，首先用流式細(xì)胞儀單獨測定地黃、大豆和辣椒基因組DNA 的相對熒光強度，然后檢測混合樣品的相對熒光強度，并根據(jù)不同樣品相對熒光強度峰值的大小，同時參考對照樣品的基因組，計算地黃的基因組大小。用流式細(xì)胞儀測定時，每個樣品收集10 000 個顆粒，檢測PI 的熒光強度，收集FL2 通道的熒光，使用FACSComp 軟件對數(shù)據(jù)進行分析和處理，并通過FL2-A通道計算得到每個樣品G1 期的熒光強度[15]。

2.2 地黃基因組survey 分析

采用改良CTAB 法提取地黃葉片的基因組DNA，用紫外分光光度計和瓊脂糖凝膠電泳分別檢測地黃基因組DNA 的濃度、純度和完整性。

檢測合格后，用超聲波震蕩方法將地黃基因組DNA 破碎至目的片段（350 bp），然后經(jīng)過末端修復(fù)、加A、加接頭、目標(biāo)片段選擇和PCR等步驟構(gòu)建小片段測序文庫；用安捷倫2100 和定量PCR檢測文庫片段大小和文庫定量，確定文庫是否符合測序標(biāo)準(zhǔn)；通過橋式PCR 將文庫固定到測序芯片上；利用北京百邁客Illumina Hiseq 4000 測序儀對構(gòu)建好的文庫進行雙端150 bp（PE 150）測序，測序所產(chǎn)生的數(shù)據(jù)經(jīng)過質(zhì)控后用于下一步信息分析。雙端測序數(shù)據(jù)經(jīng)過評估（GC 分布統(tǒng)計、質(zhì)量值Q20、Q30評估）、過濾后得到高質(zhì)量的數(shù)據(jù)，用于地黃基因組大小的評估、基因組的組裝、GC 含量的統(tǒng)計、雜合度的統(tǒng)計（以及組裝后的評估）。根據(jù)地黃基因組大小，結(jié)合流式細(xì)胞技術(shù)測定結(jié)果來估算測序深度。

Kmer 是從測序數(shù)據(jù)中提取出的長度為K的寡聚核苷酸序列，本研究中取K＝19 進行分析，在測序reads 均勻分布的前提下，根據(jù)公式，計算基因組大小（基因組大小=總堿基數(shù)/平均測序深度=總Kmer 數(shù)/平均Kmer 深度）。由于測序片段是隨機打斷的，標(biāo)準(zhǔn)的Kmer 深度分布曲線呈正態(tài)分布，根據(jù)實際曲線偏離正態(tài)分布的程度，通過Jellyfish和GenomeScope 軟件[16]對測序數(shù)據(jù)進行分析，統(tǒng)計產(chǎn)生Kmer 頻數(shù)分布數(shù)據(jù)，再使用該數(shù)據(jù)進行擬合作圖，得到Kmer 分布圖，進一步評估地黃基因組雜合度和重復(fù)序列比例。

3 結(jié)果與分析

3.1 流式細(xì)胞技術(shù)測定地黃基因組大小

以已知基因組大小的物種為對照品：大豆（基因組大小為1.1 Gb）[17-18]和辣椒（基因組大小為3.2 Gb）[19]，進行細(xì)胞流式實驗，來估算地黃基因組的大小。首先用流式細(xì)胞技術(shù)將大豆、地黃和辣椒單獨測定，檢測每個樣品的相對熒光強度（圖1-A～C），然后采用內(nèi)標(biāo)法，以大豆和辣椒為對照品，用流式細(xì)胞技術(shù)對混合樣品進行相對熒光強度的檢測，大豆、地黃和辣椒混合樣品的基因組DNA 相對熒光強度峰值分別為41.99、73.25、110.61（圖1-D）。根據(jù)流式細(xì)胞實驗結(jié)果，地黃基因組大小介于大豆和辣椒之間，估算地黃基因組大小應(yīng)在2.00～2.12 Gb。

圖1 不同樣品流式細(xì)胞術(shù)檢測結(jié)果Fig.1 Peak value image of different species

3.2 測序結(jié)果統(tǒng)計

使用地黃“沁懷1 號”樣品的基因組DNA 構(gòu)建350 bp 文庫，在Illumina Hiseq 4000 測序平臺測序并過濾得到132.79 Gb 高質(zhì)量的數(shù)據(jù)，總測序深度約為66.40×，測序數(shù)據(jù)Q20比例均在97.15%以上，Q30比例均在92.86%以上，滿足60×以上的測序數(shù)據(jù)量。

3.3 樣品污染評估

如果地黃基因組DNA 樣品存在污染不僅會降低有效數(shù)據(jù)量，同時還會影響基因組survey分析結(jié)果的準(zhǔn)確性，導(dǎo)致基因組大小、重復(fù)序列比例、雜合度和GC 含量等基因組特征評估結(jié)果出現(xiàn)較大偏差，使得基因組組裝建庫策略出現(xiàn)偏差，最終影響后續(xù)的基因組組裝效果。為了判斷提取的地黃基因組DNA 是否受到污染，從測序得到的350 bp 文庫中，隨機取10 000 條單端reads，與NT 庫（Nucleotide Sequence Database）進行BLAST 比對。BLAST 使用ncbi-blast＋2.2.29 版本，參數(shù)設(shè)置為“-num_descriptions 100-num_alignments 100-evalue 1×10-5”。能夠比對上NT 庫的reads 占提取reads數(shù)的8.74%，其中比對到芝麻Sesamum i ndicumL.、寬葉溝酸漿Erythranthe guttatusDC.、可可Theobroma cacaoL.、番茄Solanum lycopersicumL.和地黃的reads 數(shù)分別占比對上NT 庫reads 數(shù)的12.47%、4.57%、4.11%、3.08%和2.17%。芝麻是地黃的近緣物種，于2014年完成芝麻基因組測序[20]，在NT 庫中基因注釋信息豐富，所以在比對物種所占比例中最高為12.47%，在前期地黃基因的進化分析中也發(fā)現(xiàn)地黃和芝麻的親緣關(guān)系較近[21-24]，而地黃基因組信息未知，在NT 庫中基因注釋極少，所以在比對物種所占比例中僅為2.17%。比對結(jié)果中未發(fā)現(xiàn)動物、微生物等異常比對，因此該地黃基因組DNA 樣品測序數(shù)據(jù)不存在污染，可用于基因組survey分析。如果有一定比例的reads 比對到進化距離較遠(yuǎn)的物種如動物、微生物等，則判斷樣品可能存在污染，需要進一步檢查原因，不能用于基因組survey 分析。

3.4 核外DNA 含量評估

由于細(xì)胞核外也存在核酸序列，如果測序文庫中核外DNA 含量過高，會加大后期基因組組裝的難度，甚至產(chǎn)生錯誤。因此評估文庫中核外DNA 含量對判斷數(shù)據(jù)能否用于后續(xù)基因組組裝非常必要。為了評估測序數(shù)據(jù)中核外DNA 的含量，利用Illumina Hiseq 測序得到的350 bp 文庫與地黃的葉綠體序列（NC_034308.1，153 622 bp）進行SOAP 比對[25]。SOAP 使用2.21 版本，參數(shù)設(shè)置為“-m 260-x 440”。雙端比上的reads 數(shù)為17 415 302，占總reads 的1.96%，單端比上的reads 數(shù)為11 313 201，占總reads 的1.28%，這2個的比例都低于經(jīng)驗值5%。由此判斷測序數(shù)據(jù)的核外DNA 含量很低，不影響基因組survey 分析的準(zhǔn)確性和后期基因組的組裝。

3.5 地黃基因組大小、重復(fù)序列比例和雜合度評估

使用350 bp 文庫數(shù)據(jù)構(gòu)建K＝19 的Kmer 分布圖（圖2），進行地黃基因組大小、重復(fù)序列比率和雜合度的評估。由圖2 可觀察到明顯的3 個峰，第1個峰的Kmer 深度分布為27，第2 個峰的Kmer 深度分布為56，第3 個峰的Kmer 深度分布為102，這3個峰的Kmer深度分布符合1∶2∶4的趨勢，推測地黃為四倍體。假設(shè)地黃基因組為“AABB”，則代表“AA”一致，“BB”一致，但“AA” “BB”不同，根據(jù)地黃基因組呈現(xiàn)的Kmer 分布，從測序數(shù)據(jù)中得到的總Kmer 數(shù)為116 809 179 227 個，去除深度異常的Kmer 后，共113 026 782 538 個Kmer用于基因組大小估計，用Jellyfish 和GenomeScope軟件計算得到的基因組大小約2.03 Gb，與流式細(xì)胞技術(shù)測定結(jié)果一致。根據(jù)Kmer 分布情況，估計地黃基因組重復(fù)序列比例約 78.48%，雜合度約1.93%。

圖2 Kmer 分布圖預(yù)估地黃基因組大小Fig.2 Kmer analysis for estimating genome size of R.glutinosa

基因組GC 含量對二代基因組測序的隨機性有較大影響，過高（＞65%）或過低（＜25%）的GC含量會導(dǎo)致測序偏向性，嚴(yán)重影響基因組分析結(jié)果。物種GC 含量是評估基因組survey 分析準(zhǔn)確性和后續(xù)基因組組裝難度的重要指標(biāo)之一，通過對350 bp 文庫測序數(shù)據(jù)分析，地黃基因組的GC 含量約37.27%，較為適中，不會影響分析的準(zhǔn)確性。

綜上所述，地黃基因組大小約2.03 Gb，重復(fù)序列比例約78.48%，雜合度約1.93%，基因組的GC 含量約37.27%，從基因組基本結(jié)構(gòu)特征分析，地黃基因組屬于高重復(fù)、高雜合、大基因組的復(fù)雜基因組。

4 討論

基因組是一個生物體中全部的DNA 信息，包括其所有基因，通過對藥用植物的全基因組進行測序，獲得藥用植物基因組的DNA 序列信息，就能夠從分子水平研究活性成分的生物合成和代謝調(diào)控途徑，以及藥用植物之間的進化關(guān)系。隨著高通量測序技術(shù)的迅速發(fā)展，目前已有幾十種藥用植物的全基因組被解析，如大麻、靈芝、鐵皮石斛、丹參、三七、罌粟、黃花蒿等重要的藥用植物[26]。這些藥用植物全基因組序列信息的發(fā)布，為藥用植物活性成分生物合成和代謝調(diào)控途徑的研究奠定了堅實的基礎(chǔ)，極大的推動了藥用植物分子生物學(xué)的研究水平，也為藥用植物的遺傳育種提供了重要依據(jù)[27]。

由于流式細(xì)胞技術(shù)方便、快捷而且結(jié)果可靠，是測定植物基因大小的首選方法[14]，目前應(yīng)用流式細(xì)胞技術(shù)，測定了黃芪[28]、柴胡[29]和人參[30]等藥用植物基因組的大小。除了流式細(xì)胞技術(shù)以外，利用高通量測序技術(shù)對植物基因組進行survey 分析，也是檢測植物基因組大小和特征的有效方法。通過對基因組進行survey分析，提前獲得基因組的大小、重復(fù)序列、GC 含量等信息，為后續(xù)進行全基因組深度測序時選擇合適的測序策略提供依據(jù)。目前利用高通量測序方法，已對羅漢果[31]、馬鞍藤[32]、三島柴胡[33]等藥用植物的基因組進行survey分析，獲得了這些植物基因組大小、復(fù)雜程度以及基因組的其他特征等信息。

玄參科有4500 種植物，包括多種藥用植物，如地黃、玄參、洋地黃、胡黃連、陰行草等，具有重要的藥用價值和經(jīng)濟價值[34]，但是對于這些藥用植物進行遺傳育種和基礎(chǔ)研究的主要障礙在于缺乏參考基因組序列，截至目前還沒有玄參科植物的參考基因組發(fā)布。準(zhǔn)確測定地黃基因組大小是后續(xù)進行地黃全基因組深度測序和遺傳分析的基礎(chǔ)，目前在許多基因組測序項目中都使用高通量測序的Kmer 分析估算基因組大小[35]。本研究應(yīng)用流式細(xì)胞技術(shù)估算地黃基因組大小應(yīng)在2.00～2.12 Gb，同時結(jié)合高通量測序的Kmer 分析，地黃基因組大小修正為2.03 Gb，而且通過Kmer分析還能獲得地黃基因組的重復(fù)序列78.48%、雜合度1.93%、GC 含量37.27%等基因組特征信息。通過將流式細(xì)胞技術(shù)和基因組survey分析2種方法結(jié)合在一起，對地黃基因組大小和特征進行分析，提高了實驗結(jié)果的可靠性。從基因組基本結(jié)構(gòu)特征上看，地黃基因組屬于高重復(fù)、高雜合、大基因組的復(fù)雜基因組。由于二代高通量測序技術(shù)限制，測序讀長太短，僅有200～500 bp，對于地黃這種大型復(fù)雜基因組的組裝難度大、花費多，而從頭組裝de novo組裝對算法的要求更高，所以針對地黃基因組雜合度高、重復(fù)序列比例高、基因組大的特點，可采用3 代高通量測序技術(shù)（PacBio Sequel）同時結(jié)合染色質(zhì)區(qū)域捕獲（chromosome conformation capture，Hi-C）技術(shù)來完善三代高通量測序基因組組裝結(jié)果，將其組裝到染色體水平[36]。PacBio Sequel 測序平均長度達(dá)12～15 kb，能夠跨越基因組高重復(fù)區(qū)域和高復(fù)雜區(qū)域，減少拼接成本，同時結(jié)合Hi-C 技術(shù)，能夠獲得高質(zhì)量染色體水平的基因組序列，可以闡明物種間基因組的演化歷程，揭示物種間的進化關(guān)系。本研究通過流式細(xì)胞技術(shù)和基因組survey 分析，獲得的地黃基因組大小和特征等信息，為今后繪制地黃全基因組的精細(xì)圖譜奠定了基礎(chǔ)，也為研究地黃屬植物的多倍體現(xiàn)象以及遺傳進化關(guān)系提供依據(jù)。

利益沖突所有作者均聲明不存在利益沖突