川芎基因組survey 測序及其特征分析

毛常清，沙秀芬, ，黃 靜，陶 珊，彭 芳，李 群，張 超，袁 燦*

1.四川省農(nóng)業(yè)科學(xué)院 經(jīng)濟(jì)作物育種栽培研究所，四川 成都 610300

2.四川師范大學(xué) 生命科學(xué)學(xué)院，四川 成都 610101

川芎Ligusticum chuanxiongHort.是傘形科藁本屬草本植物，以干燥的根莖入藥，始載于《神農(nóng)本草經(jīng)》，已有1500 多年的種植歷史，是我國傳統(tǒng)的大宗中藥材。川芎廣泛分布于我國四川彭州、什邡、眉山等地，在云南、貴州、廣西、湖北、江西等省也少量引進(jìn)種植。川芎性溫，味辛、微苦，具有活血行氣、祛風(fēng)止痛的功效，其主要化學(xué)成分包含揮發(fā)油、生物堿、有機(jī)酸和多糖等，在臨床上廣泛地用于治療心臟病、腦梗死及尿路結(jié)石等疾病[1-3]。

基因組序列是研究一個(gè)物種遺傳背景的基礎(chǔ)，隨著高通量測技術(shù)的逐漸成熟，許多植物基因組序列相繼被發(fā)表，包括大量藥用植物基因組。如通過高通量測序技術(shù)成功測算出靈芝、丹參、人參、三七、天麻、穿心蓮、黃花蒿、廣藿香、鐵皮石斛等幾十種重要藥用植物的基因組大小和特征[4-5]。但我國藥用植物種類豐富，約占中藥材資源總數(shù)的87%[6]，同大多數(shù)藥用植物一樣，現(xiàn)報(bào)道川芎的研究主要集中在化學(xué)成分[2]和藥理藥效機(jī)制上[7]，分子生物學(xué)方面僅開展了川芎轉(zhuǎn)錄組分析和利用通用引物分析其遺傳多樣性[8-10]，在分子遺傳學(xué)系統(tǒng)研究上存在較大空白。雖然藥用植物基因測序技術(shù)的應(yīng)用為川芎全基因組測序提供了技術(shù)基礎(chǔ)，但由于川芎基因組結(jié)構(gòu)龐大，遺傳背景復(fù)雜，直接進(jìn)行全基因組測序存在一定困難，因此在進(jìn)行全基因組測序之前，有必要對川芎基因組大小進(jìn)行調(diào)研。

本研究采用流式細(xì)胞術(shù)和Illumina Hiseq 2500高通量測序技術(shù)相結(jié)合的方式對川芎基因組大小進(jìn)行估算，對所得的基因組數(shù)據(jù)進(jìn)行K-mer 分析、基因組預(yù)測及注釋，關(guān)注阿魏酸等成分合成途徑的基因注釋，為進(jìn)一步全基因組精細(xì)測序和藥效成分合成分子機(jī)制研究提供參考依據(jù)和基因資源。

1 材料

分別采集2 個(gè)月左右的綠豆Vigna radiata(Linn.) Wilczek.、陸地棉Gossypium hirsutumL.幼嫩葉片和1 個(gè)月左右川芎幼嫩葉片（種植于四川省農(nóng)科院經(jīng)濟(jì)作物育種栽培研究所基地）用于流式細(xì)胞分析，并采集川芎葉片用于基因組測序。綠豆由四川省農(nóng)業(yè)科學(xué)院經(jīng)濟(jì)作物育種栽培研究所葉鵬盛研究員鑒定為豆科豇豆屬植物綠豆，陸地棉由中國農(nóng)業(yè)科學(xué)院棉花研究所杜雄明研究員鑒定為錦葵科棉屬植物陸地棉，川芎由四川農(nóng)業(yè)大學(xué)陳興福教授鑒定為傘形科藁本屬植物川芎。

2 方法

2.1 樣品的制備

分別選取綠豆、棉花、川芎幼嫩葉片各2 份，每份120 mg，洗凈置于預(yù)冷的培養(yǎng)皿中，向培養(yǎng)皿中加入1 mL 預(yù)冷的OttoI 細(xì)胞裂解液，快速切碎葉片后，用移液槍上下吹打混勻（避免氣泡），所得的提取液用42 μm 尼龍膜濾過到離心管中，低速冷凍離心后，棄上清液，向沉淀中加入1 mL冰浴的OttoII緩沖液重懸細(xì)胞，放置4 ℃?zhèn)溆谩?/p>

采用改良CTAB 法提取川芎基因組DNA，使用NanoDrop 2000C 超微量分光光度計(jì)和1%瓊脂糖凝膠電泳檢測DNA 濃度及完整性。

2.2 流式細(xì)胞術(shù)測定川芎基因組大小

將上述制備好的細(xì)胞懸浮液樣品中加50 μL 1 mg/mL RNase，50 μL、50 μg/mL PI（DNA 熒光染料碘化丙啶，預(yù)先經(jīng)0.22 μm 微孔濾膜濾過，-20 ℃保存），混勻，4 ℃避光染色 10 min。隨后用FACSCalibur 流式細(xì)胞儀檢測PI 在488 nm 激發(fā)光下發(fā)出的熒光，CellQuest 軟件捕捉熒光信號數(shù)據(jù)，ModFit 軟件分析結(jié)果。測定基因組大小。

待測樣品DNA 量=對照DNA 量×待測樣品的熒光強(qiáng)度/對照品的熒光強(qiáng)度

2.3 基因組測序和質(zhì)量評估

構(gòu)建270 bp 和500 bp 的小片段文庫，利用Hiseq2500 測序技術(shù)對文庫進(jìn)行雙端測序。從文庫中隨機(jī)取10 000 條單端read 與NCBI 數(shù)據(jù)庫中的核苷酸數(shù)據(jù)庫（NT）進(jìn)行BLAST[11]比對，判斷樣本是否被外源物種污染。數(shù)據(jù)測序由北京百邁客生物科技有限公司完成。

2.4 基因組大小、重復(fù)序列和雜合率預(yù)測

選取K值為21 對基因組大小、重復(fù)序列、雜合率進(jìn)行預(yù)測。用程序Jellyfish 計(jì)算K-mer 分布，并通過K-mer 分布曲線初步評估基因組重復(fù)序列含量和雜合度。

基因組大小=K-mer 總數(shù)/K-mer 期望深度值

2.5 基因組初步組裝和GC-Depth 分布分析

利用SOAPdenovo 進(jìn)行組裝得到contig，利用雙末端信息進(jìn)行g(shù)ap 填充，將無overlap 關(guān)系的contig 拼接組裝成scaffold，獲得含有N（重復(fù)序列）的初級基因組序列。過濾后的read 比對到已組裝好的基因序列上，獲得堿基深度，以10 kb 為窗口，在序列上無重復(fù)前進(jìn)，計(jì)算每個(gè)窗口的平均深度與GC 含量，做出GC_depth 圖[12]。

2.6 重復(fù)序列分析

對測序所得的數(shù)據(jù)進(jìn)行重復(fù)序列分析。使用4個(gè)互補(bǔ)程序 LTR_FINDER[13]、MITE-Hunter[14]，RepeatScout[15]和PILER-DF[16]構(gòu)建川芎重復(fù)序列文庫，隨后由PASTEClassifier[17]分類，并與Repbase[18]轉(zhuǎn)座因子庫結(jié)合起來作為最終的庫，然后運(yùn)行軟件Repeat-Masker[19]在最終文庫中找到同源重復(fù)序列。

2.7 基因預(yù)測和注釋

基因從頭預(yù)測，在濾過掉小于1000 bp 大小的scaffold 后，用Genscan 和Augustus 軟件通過擬南芥的訓(xùn)練集預(yù)測川芎基因。將預(yù)測到的基因比對到非冗余蛋白序列（non-redundant，Nr）、真核生物蛋白直系同源簇（clusters of euKaryotic orthologous groups，KOG）、基因本體（gene ontology，GO）、swiss-prot 和京都基因與基因組百科全書（Kyoto encyclopedia of genes and genomes，KEGG）等數(shù)據(jù)庫進(jìn)行BLAST 分析來對預(yù)測基因進(jìn)行注釋。然后，通過KEGG 在線網(wǎng)站（http://www.genome.jp/kegg/）檢索KEGG 途徑，使用Blast2GO 軟件處理得到的Nr 注釋結(jié)果進(jìn)行 GO 分類，使用在線網(wǎng)站（http://www.ncbi.nlm.nih.gov/COG/）處理KOG 注釋結(jié)果進(jìn)行KOG 分類。

2.8 基因家族鑒定及系統(tǒng)發(fā)育樹構(gòu)建

使用ORTHOMCL v2.0.9 將預(yù)測到的川芎蛋白序列與丹參Salvia miltiorrhizaBunge、胡蘿卜Daucus carotavar.sativaHoffm.、葡萄Vitis viniferaL.和擬南芥Arabidopsis thalianaL.等4 種植物中氨基酸數(shù)目大于50 的序列匯集到一個(gè)蛋白數(shù)據(jù)庫中，通過blastp 比對獲得所有物種蛋白序列之間的相似性關(guān)系，E值為1×10-5，去掉序列一致度小于30%，覆蓋率小于30%的序列；并對比對結(jié)果進(jìn)行聚類，默認(rèn)膨脹系數(shù)為1.5。在ORTHOMCL 結(jié)果中檢索單拷貝基因家族，利用單拷貝同源基因基于最大似然法（maximum likelihood，ML）進(jìn)行進(jìn)化樹構(gòu)建。

3 結(jié)果與分析

3.1 流式細(xì)胞術(shù)檢測基因組大小

采用綠豆和陸地棉作為內(nèi)標(biāo)植物，用流式細(xì)胞儀檢測其混合樣品的熒光強(qiáng)度。已知綠豆基因組大小為579 Mb[20]，陸地棉基因組為2173 Mb[21]，通過2 個(gè)內(nèi)標(biāo)植物計(jì)算出川芎的基因組大小分別為525/78×579=3 897.12 Mb，變異系數(shù)（coefficient of variation，CV）為3.40%（圖1）；613/439×2173=3 034.28 Mb，CV 為2.01%（圖2）。

圖1 綠豆與川芎混合樣品流式細(xì)胞術(shù)測定Fig.1 FCM determination for mixed samples of V.radiata and L.chuanxiong

圖2 陸地棉花與川芎混合樣品流式細(xì)胞術(shù)測定Fig.2 FCM determination for mixed samples of G.hirsutum and L.chuanxiong

3.2 川芎基因組測序數(shù)據(jù)統(tǒng)計(jì)及質(zhì)量評估

使用川芎基因組DNA 構(gòu)建270 bp 的文庫，通過Illumina Hiseq2500 測序平臺測序并過濾得到222.04 Gb 高質(zhì)量的數(shù)據(jù)，測序深度為72.59 X，測序數(shù)據(jù)Q20 比例均在97.59%以上，Q30 比例均在94.74%以上。隨機(jī)篩選的1000 條單端read 能夠比對上NT 核酸數(shù)據(jù)庫的read 占總read 的7.62%，其中比對上野胡蘿Daucus carotaL.、細(xì)葉藁本Ligusticum tenuissimum(Nakai) Kitagawa 的read 數(shù)分別占比對上NT 庫reads數(shù)的61.81%、3.01%，且未發(fā)現(xiàn)動(dòng)物、微生物等異常比對，說明樣本不存在污染。

3.3 川芎基因組特征

通過270 bp 文庫數(shù)據(jù)構(gòu)建K=21 的K-mer 分布圖（圖3），K-mer 深度62 X 為主峰（由于雜合度較高，本研究對主峰判斷參考流式細(xì)胞結(jié)果），測序得到K-mer 的總數(shù)為189 618 848 178，估算出川芎基因組大小為3 058.37 Mb，與流式細(xì)胞實(shí)驗(yàn)中以陸地棉為內(nèi)標(biāo)植物測定的結(jié)果相近。在主峰對應(yīng)深度的1/2 處出現(xiàn)明顯的雜合峰，深度為31 X，說明川芎基因組具有較高的雜合度。進(jìn)一步根據(jù)變異數(shù)目占基因組大小的比例即雜合度估算該基因組雜合率，在所有有效數(shù)據(jù)中檢測到每139 個(gè)堿基對中就有3 個(gè)SNP，估算出川芎基因組具有較高雜合率，約為2.16%。K-mer 分布存在較長拖尾，暗示川芎基因組存在較高的重復(fù)率，估算重復(fù)序列的基因組大小估計(jì)為 2 440.13 Mb，約為川芎基因組的79.79%。說明川芎屬于高重復(fù)、高雜合、大基因組等基因組特征的復(fù)雜物種。

圖3 K-mer 分布估算基因組大小Fig.3 K-mer distribution to estimate genome size

3.4 川芎基因組組裝和GC 含量分析

使用222.04 Gb 高質(zhì)量數(shù)據(jù)，基于K=41 組裝產(chǎn)生12 973 787 個(gè)contig 和8 119 089 個(gè)scaffold。contig N50 為286 bp，N90 為131 bp，最大長度達(dá)到26 842 bp，總長度為3 198 598 874 bp；scaffold N50 為493 bp，N90 為191 bp，最大長度為29 779 bp，總長度為3 284 536 081 bp。其中，contig 和scaffold 的N50 值相對較短，可能是由于川芎高雜合率引起的。通過對組裝的contig 進(jìn)行GC 含量的統(tǒng)計(jì)，結(jié)果顯示川芎基因組的GC 含量約為36.32%（圖4），說明測序不具有明顯的GC 偏向性，不影響測序分析的準(zhǔn)確性。

圖4 川芎基因組GC 含量和平均測序深度Fig.4 GC content and average sequencing depth of L.chuanxiong genome

3.5 川芎基因組重復(fù)序列分析

重復(fù)序列檢測顯示其總長度為2 250 901 770 bp，約為基因組大小的73.60%，低于K-mer 分析估算的重復(fù)序列含量，原因可能是組裝效果的限制，導(dǎo)致組裝過程中重復(fù)序列損失6.19%。注釋上，能夠找到明確重復(fù)序列元件的總長度約為2 026.25 Mb。其中，長末端重復(fù)序列（long terminal repeated，LTR）是最豐富的重復(fù)元件，占基因組的14.14%，其次是長散在重復(fù)元件（long interspersed nuclear elements，LINE），占基因組的0.49%。SSR 重復(fù)序列總長度約49.41 Mb，占基因組的1.62%，占重復(fù)序列的2.44%。

3.6 基因預(yù)測和注釋

總共預(yù)測到34 864 個(gè)基因，總長度為27 532 625 bp。在預(yù)測的基因中，查找到79 130 個(gè)外顯子，總長度為18 557 406 bp；79 129 個(gè)內(nèi)含子，總長度為8 975 219 bp。有30 737 個(gè)基因在功能數(shù)據(jù)庫中能比對到注釋信息，Nr 具有最高的注釋率（30 584，87.72%）KEGG 具有最低的注釋率（9623，27.60%）。在預(yù)測的基因中，15 184 個(gè)基因被分類為GO 功能類別；15 598 個(gè)基因被分類為KOG 功能類別；9623個(gè)基因注釋到125 個(gè)KEGG 代謝途徑。在GO 功能類別中，包含分子功能、細(xì)胞組成和生物過程（圖5）；在KOG 功能類別中，通用功能預(yù)測的基因最多，其次是翻譯后修飾、蛋白質(zhì)周轉(zhuǎn)、分子伴侶和信號轉(zhuǎn)導(dǎo)機(jī)制（圖6）；在參與基因數(shù)目最多的前10條KEGG 代謝途徑中，注釋基因分別參與核糖體代謝、植物激素信號轉(zhuǎn)導(dǎo)、內(nèi)質(zhì)網(wǎng)蛋白質(zhì)、剪接體、碳代謝、氨基酸生物合成、RNA 轉(zhuǎn)運(yùn)、氧化磷酸化、植物-病原菌互作、淀粉和蔗糖代謝等途徑。

圖5 川芎基因的GO 注釋Fig.5 GO annotations of L.chuanxiong genes

圖6 川芎基因功能注釋KOG 功能分類Fig.6 KOG functional classification of L.chuanxiong genes

3.7 基因家族鑒定及系統(tǒng)發(fā)育分析

通過與胡蘿卜、丹參、葡萄和擬南芥等物種蛋白序列比對（圖7-A、B），川芎中常見的基因家族中的基因數(shù)量小于同科的胡蘿卜，每個(gè)基因家族的平均基因數(shù)量與其他物種相當(dāng)，但川芎獨(dú)特基因家族的數(shù)量比其他物種的獨(dú)特基因家族中的基因數(shù)量要大得多，共計(jì)2058 個(gè)基因家族。所有5 個(gè)物種共有的基因家族為7112 個(gè)，其中2001 個(gè)基因是單拷貝基因，即每個(gè)基因家族中只存在一個(gè)直系同源基因，可用于系統(tǒng)發(fā)育推斷和發(fā)散時(shí)間估計(jì)。對2001 個(gè)單拷貝基因利用ML 構(gòu)建系統(tǒng)發(fā)育樹（圖7-C）。從川芎系統(tǒng)發(fā)育樹分析，川芎、胡蘿卜、丹參、葡萄和擬南芥來自于共同祖先。其中，川芎與胡蘿卜最晚與其他物種發(fā)生分歧，二者進(jìn)化分支長度差異最小，分歧時(shí)間更短，親緣關(guān)系更近。從遺傳變異度上來看，胡蘿卜所在的分支最短，遺傳變異度最小，進(jìn)化距離最近，川芎的遺傳變異度和進(jìn)化距離僅次于胡蘿卜，丹參的遺傳變異度最高，進(jìn)化距離最遠(yuǎn)。

圖7 川芎基因家族鑒定及系統(tǒng)發(fā)育樹Fig.7 Gene family identification and phylogenetic tree of L.chuanxiong

3.8 參與阿魏酸生物合成的基因

阿魏酸是川芎的主要藥效成分，屬于苯丙烷生物合成途徑，是木質(zhì)素合成的中間體。近年來研究者證實(shí)了在川芎、當(dāng)歸等傘形科物種中，阿魏酸的合成途徑主要包括COMT 途徑和CCoAMT 途徑[9,22-24]，主要參與的酶有苯丙氨酸解氨酶、肉桂酸-4-羥化酶、香豆酸-3-羥化酶、咖啡酸-O-甲基轉(zhuǎn)移酶、4-香豆酸輔酶A 連接酶、莽草酸O-羥基肉桂酰轉(zhuǎn)移酶、奎寧O-羥基肉桂酰轉(zhuǎn)移酶、咖啡酰輔酶A-O-甲基轉(zhuǎn)移酶、肉桂酰輔酶A 還原酶、醛脫氫酶家族成員C4。為挖掘川芎中參與阿魏酸生物合成的主要基因，本研究從注釋的KEGG 代謝通路中提取苯丙烷生物合成參考途徑（map00940）相關(guān)數(shù)據(jù)，對阿魏酸合成過程的相關(guān)基因進(jìn)行了檢索，共有53 個(gè)功能基因覆蓋到阿魏酸合成途徑中（表1）。其中，編碼咖啡酸-O-甲基轉(zhuǎn)移酶、4-香豆酸輔酶A 連接酶、莽草酸O-羥基肉桂酰轉(zhuǎn)移酶、咖啡酰輔酶A-O-甲基轉(zhuǎn)移酶等酶的基因均有多個(gè)同源拷貝存在，預(yù)測川芎基因組在進(jìn)化的過程中某一時(shí)間點(diǎn)發(fā)生過基因擴(kuò)張。

4 討論

基因組包含了一個(gè)生物體所有基因的總和，了解生物體基因組信息有助于深入了解生物體遺傳、進(jìn)化、生物合成、次生代謝等全部過程。通過基因組測序技術(shù)可以對特定物種基因組進(jìn)行測序，利用生物信息學(xué)方法對測序序列進(jìn)行拼接和組裝，最終獲得該物種基因組序列，進(jìn)而了解其基因組信息[25]。2009 年，陳士林團(tuán)隊(duì)[6]首次提出本草基因組計(jì)劃，此后，越來越多的學(xué)者開始在藥用植物基因組學(xué)研究上投入大量精力。藥用植物全基因組水平研究，有助于闡明藥用植物活性成分生物合成和代謝調(diào)控途徑之間的關(guān)系，為具有相似有效成分和藥理活性的近緣物種間的系統(tǒng)發(fā)育關(guān)系研究奠定了基礎(chǔ)，也為藥用植物的遺傳育種和基因資源保護(hù)提供了重要依據(jù)[26]。

本研究利用流式細(xì)胞術(shù)和高通量測序相結(jié)合的方式對川芎基因組進(jìn)行測算。通過流式細(xì)胞術(shù)估算出川芎基因組大小分別為3 897.12 Mb和3 034.28 Mb，通過高通量測序的K-mer 分析，綜合2 個(gè)分析結(jié)果估算川芎基因組大小為3 058.37 Mb，屬于基因組較大的物種。本研究測得川芎基因組GC 含量為36.32%，處于植物基因組GC 含量應(yīng)介于25%～65%的合理范圍[27]，說明川芎基因組測序的結(jié)果和組裝是正確可靠的。通過K-mer 分析，估算出川芎重復(fù)序列含量為79.79%，雜合率約為2.16%，與地黃[28]、黃芪[29]等藥用植物類似，呈現(xiàn)高重復(fù)、高雜合的基因組特征，進(jìn)一步說明川芎基因?qū)儆诟咧貜?fù)、高雜合、基因組龐大的物種。

在基因預(yù)測、注釋和基因家族的鑒定中，本研究共預(yù)測到川芎編碼蛋白基因34 864 個(gè)，遠(yuǎn)高于其傘形科親緣關(guān)系較近的胡蘿卜（32 113 個(gè)）[30]，可能是因?yàn)榻M裝的都是短片段測序文庫，川芎中的基因數(shù)量可能被高估了。此外，本研究中測序完成后從頭組裝中產(chǎn)生的contig N50 為286 bp，scaffold N50 為493 bp，明顯較預(yù)期短，這與廣藿香[4]、羅漢果[12]等藥用植物的全基因組調(diào)研結(jié)果一致。提示對于川芎這種具有復(fù)雜基因組特征的物種來說，利用二代高通量測序?qū)ζ淙蚪M進(jìn)行精確測序仍然存在技術(shù)難度。因此，提高川芎基因組的測序深度和組裝質(zhì)量，建議后續(xù)的研究可采用二代和三代測序技術(shù)相結(jié)合，并利用全基因組染色體構(gòu)象捕獲技術(shù)（high-through chromosome conformation capture，Hi-C），解析全基因組范圍內(nèi)整個(gè)染色質(zhì)DNA 在空間位置上的關(guān)系，獲得完整準(zhǔn)確的全基因組圖譜[31]，得到高質(zhì)量的川芎基因組序列。

本研究首次利用流式細(xì)胞技術(shù)和基因組survey分析，初步獲得川芎基因組大小和結(jié)構(gòu)特征，即基因組龐大、序列重復(fù)率高、序列雜合度高，為下一步進(jìn)行全基因組精細(xì)測序奠定基礎(chǔ)。本研究中組裝產(chǎn)生的大量川芎基因組序列和注釋基因?yàn)楹罄m(xù)分子標(biāo)記的開發(fā)和基因功能研究提供了大量的資源。同時(shí)，本研究挖掘了阿魏酸合成途徑中的參與基因，對川芎阿魏酸生物合成途徑潛在分子機(jī)制的初步研究，為進(jìn)一步研究川芎生物學(xué)和選育具有優(yōu)良藥用性狀的品種奠定了基礎(chǔ)。

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