毛常清,沙秀芬, ,黃 靜,陶 珊,彭 芳,李 群,張 超,袁 燦*
1.四川省農(nóng)業(yè)科學(xué)院 經(jīng)濟(jì)作物育種栽培研究所,四川 成都 610300
2.四川師范大學(xué) 生命科學(xué)學(xué)院,四川 成都 610101
川芎Ligusticum chuanxiongHort.是傘形科藁本屬草本植物,以干燥的根莖入藥,始載于《神農(nóng)本草經(jīng)》,已有1500 多年的種植歷史,是我國(guó)傳統(tǒng)的大宗中藥材。川芎廣泛分布于我國(guó)四川彭州、什邡、眉山等地,在云南、貴州、廣西、湖北、江西等省也少量引進(jìn)種植。川芎性溫,味辛、微苦,具有活血行氣、祛風(fēng)止痛的功效,其主要化學(xué)成分包含揮發(fā)油、生物堿、有機(jī)酸和多糖等,在臨床上廣泛地用于治療心臟病、腦梗死及尿路結(jié)石等疾病[1-3]。
基因組序列是研究一個(gè)物種遺傳背景的基礎(chǔ),隨著高通量測(cè)技術(shù)的逐漸成熟,許多植物基因組序列相繼被發(fā)表,包括大量藥用植物基因組。如通過高通量測(cè)序技術(shù)成功測(cè)算出靈芝、丹參、人參、三七、天麻、穿心蓮、黃花蒿、廣藿香、鐵皮石斛等幾十種重要藥用植物的基因組大小和特征[4-5]。但我國(guó)藥用植物種類豐富,約占中藥材資源總數(shù)的87%[6],同大多數(shù)藥用植物一樣,現(xiàn)報(bào)道川芎的研究主要集中在化學(xué)成分[2]和藥理藥效機(jī)制上[7],分子生物學(xué)方面僅開展了川芎轉(zhuǎn)錄組分析和利用通用引物分析其遺傳多樣性[8-10],在分子遺傳學(xué)系統(tǒng)研究上存在較大空白。雖然藥用植物基因測(cè)序技術(shù)的應(yīng)用為川芎全基因組測(cè)序提供了技術(shù)基礎(chǔ),但由于川芎基因組結(jié)構(gòu)龐大,遺傳背景復(fù)雜,直接進(jìn)行全基因組測(cè)序存在一定困難,因此在進(jìn)行全基因組測(cè)序之前,有必要對(duì)川芎基因組大小進(jìn)行調(diào)研。
本研究采用流式細(xì)胞術(shù)和Illumina Hiseq 2500高通量測(cè)序技術(shù)相結(jié)合的方式對(duì)川芎基因組大小進(jìn)行估算,對(duì)所得的基因組數(shù)據(jù)進(jìn)行K-mer 分析、基因組預(yù)測(cè)及注釋,關(guān)注阿魏酸等成分合成途徑的基因注釋,為進(jìn)一步全基因組精細(xì)測(cè)序和藥效成分合成分子機(jī)制研究提供參考依據(jù)和基因資源。
分別采集2 個(gè)月左右的綠豆Vigna radiata(Linn.) Wilczek.、陸地棉Gossypium hirsutumL.幼嫩葉片和1 個(gè)月左右川芎幼嫩葉片(種植于四川省農(nóng)科院經(jīng)濟(jì)作物育種栽培研究所基地)用于流式細(xì)胞分析,并采集川芎葉片用于基因組測(cè)序。綠豆由四川省農(nóng)業(yè)科學(xué)院經(jīng)濟(jì)作物育種栽培研究所葉鵬盛研究員鑒定為豆科豇豆屬植物綠豆,陸地棉由中國(guó)農(nóng)業(yè)科學(xué)院棉花研究所杜雄明研究員鑒定為錦葵科棉屬植物陸地棉,川芎由四川農(nóng)業(yè)大學(xué)陳興福教授鑒定為傘形科藁本屬植物川芎。
分別選取綠豆、棉花、川芎幼嫩葉片各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%瓊脂糖凝膠電泳檢測(cè)DNA 濃度及完整性。
將上述制備好的細(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ì)胞儀檢測(cè)PI 在488 nm 激發(fā)光下發(fā)出的熒光,CellQuest 軟件捕捉熒光信號(hào)數(shù)據(jù),ModFit 軟件分析結(jié)果。測(cè)定基因組大小。
待測(cè)樣品DNA 量=對(duì)照DNA 量×待測(cè)樣品的熒光強(qiáng)度/對(duì)照品的熒光強(qiáng)度
構(gòu)建270 bp 和500 bp 的小片段文庫,利用Hiseq2500 測(cè)序技術(shù)對(duì)文庫進(jìn)行雙端測(cè)序。從文庫中隨機(jī)取10 000 條單端read 與NCBI 數(shù)據(jù)庫中的核苷酸數(shù)據(jù)庫(NT)進(jìn)行BLAST[11]比對(duì),判斷樣本是否被外源物種污染。數(shù)據(jù)測(cè)序由北京百邁客生物科技有限公司完成。
選取K值為21 對(duì)基因組大小、重復(fù)序列、雜合率進(jìn)行預(yù)測(cè)。用程序Jellyfish 計(jì)算K-mer 分布,并通過K-mer 分布曲線初步評(píng)估基因組重復(fù)序列含量和雜合度。
基因組大小=K-mer 總數(shù)/K-mer 期望深度值
利用SOAPdenovo 進(jìn)行組裝得到contig,利用雙末端信息進(jìn)行g(shù)ap 填充,將無overlap 關(guān)系的contig 拼接組裝成scaffold,獲得含有N(重復(fù)序列)的初級(jí)基因組序列。過濾后的read 比對(duì)到已組裝好的基因序列上,獲得堿基深度,以10 kb 為窗口,在序列上無重復(fù)前進(jìn),計(jì)算每個(gè)窗口的平均深度與GC 含量,做出GC_depth 圖[12]。
對(duì)測(cè)序所得的數(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ù)序列。
基因從頭預(yù)測(cè),在濾過掉小于1000 bp 大小的scaffold 后,用Genscan 和Augustus 軟件通過擬南芥的訓(xùn)練集預(yù)測(cè)川芎基因。將預(yù)測(cè)到的基因比對(duì)到非冗余蛋白序列(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 分析來對(duì)預(yù)測(cè)基因進(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 分類。
使用ORTHOMCL v2.0.9 將預(yù)測(cè)到的川芎蛋白序列與丹參Salvia miltiorrhizaBunge、胡蘿卜Daucus carotavar.sativaHoffm.、葡萄Vitis viniferaL.和擬南芥Arabidopsis thalianaL.等4 種植物中氨基酸數(shù)目大于50 的序列匯集到一個(gè)蛋白數(shù)據(jù)庫中,通過blastp 比對(duì)獲得所有物種蛋白序列之間的相似性關(guān)系,E值為1×10-5,去掉序列一致度小于30%,覆蓋率小于30%的序列;并對(duì)比對(duì)結(jié)果進(jìn)行聚類,默認(rèn)膨脹系數(shù)為1.5。在ORTHOMCL 結(jié)果中檢索單拷貝基因家族,利用單拷貝同源基因基于最大似然法(maximum likelihood,ML)進(jìn)行進(jìn)化樹構(gòu)建。
采用綠豆和陸地棉作為內(nèi)標(biāo)植物,用流式細(xì)胞儀檢測(cè)其混合樣品的熒光強(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ù)測(cè)定Fig.1 FCM determination for mixed samples of V.radiata and L.chuanxiong
圖2 陸地棉花與川芎混合樣品流式細(xì)胞術(shù)測(cè)定Fig.2 FCM determination for mixed samples of G.hirsutum and L.chuanxiong
使用川芎基因組DNA 構(gòu)建270 bp 的文庫,通過Illumina Hiseq2500 測(cè)序平臺(tái)測(cè)序并過濾得到222.04 Gb 高質(zhì)量的數(shù)據(jù),測(cè)序深度為72.59 X,測(cè)序數(shù)據(jù)Q20 比例均在97.59%以上,Q30 比例均在94.74%以上。隨機(jī)篩選的1000 條單端read 能夠比對(duì)上NT 核酸數(shù)據(jù)庫的read 占總read 的7.62%,其中比對(duì)上野胡蘿Daucus carotaL.、細(xì)葉藁本Ligusticum tenuissimum(Nakai) Kitagawa 的read 數(shù)分別占比對(duì)上NT 庫reads數(shù)的61.81%、3.01%,且未發(fā)現(xiàn)動(dòng)物、微生物等異常比對(duì),說明樣本不存在污染。
通過270 bp 文庫數(shù)據(jù)構(gòu)建K=21 的K-mer 分布圖(圖3),K-mer 深度62 X 為主峰(由于雜合度較高,本研究對(duì)主峰判斷參考流式細(xì)胞結(jié)果),測(cè)序得到K-mer 的總數(shù)為189 618 848 178,估算出川芎基因組大小為3 058.37 Mb,與流式細(xì)胞實(shí)驗(yàn)中以陸地棉為內(nèi)標(biāo)植物測(cè)定的結(jié)果相近。在主峰對(duì)應(yīng)深度的1/2 處出現(xiàn)明顯的雜合峰,深度為31 X,說明川芎基因組具有較高的雜合度。進(jìn)一步根據(jù)變異數(shù)目占基因組大小的比例即雜合度估算該基因組雜合率,在所有有效數(shù)據(jù)中檢測(cè)到每139 個(gè)堿基對(duì)中就有3 個(gè)SNP,估算出川芎基因組具有較高雜合率,約為2.16%。K-mer 分布存在較長(zhǎng)拖尾,暗示川芎基因組存在較高的重復(fù)率,估算重復(fù)序列的基因組大小估計(jì)為 2 440.13 Mb,約為川芎基因組的79.79%。說明川芎屬于高重復(fù)、高雜合、大基因組等基因組特征的復(fù)雜物種。
圖3 K-mer 分布估算基因組大小Fig.3 K-mer distribution to estimate genome size
使用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,最大長(zhǎng)度達(dá)到26 842 bp,總長(zhǎng)度為3 198 598 874 bp;scaffold N50 為493 bp,N90 為191 bp,最大長(zhǎng)度為29 779 bp,總長(zhǎng)度為3 284 536 081 bp。其中,contig 和scaffold 的N50 值相對(duì)較短,可能是由于川芎高雜合率引起的。通過對(duì)組裝的contig 進(jìn)行GC 含量的統(tǒng)計(jì),結(jié)果顯示川芎基因組的GC 含量約為36.32%(圖4),說明測(cè)序不具有明顯的GC 偏向性,不影響測(cè)序分析的準(zhǔn)確性。
圖4 川芎基因組GC 含量和平均測(cè)序深度Fig.4 GC content and average sequencing depth of L.chuanxiong genome
重復(fù)序列檢測(cè)顯示其總長(zhǎng)度為2 250 901 770 bp,約為基因組大小的73.60%,低于K-mer 分析估算的重復(fù)序列含量,原因可能是組裝效果的限制,導(dǎo)致組裝過程中重復(fù)序列損失6.19%。注釋上,能夠找到明確重復(fù)序列元件的總長(zhǎng)度約為2 026.25 Mb。其中,長(zhǎng)末端重復(fù)序列(long terminal repeated,LTR)是最豐富的重復(fù)元件,占基因組的14.14%,其次是長(zhǎng)散在重復(fù)元件(long interspersed nuclear elements,LINE),占基因組的0.49%。SSR 重復(fù)序列總長(zhǎng)度約49.41 Mb,占基因組的1.62%,占重復(fù)序列的2.44%。
總共預(yù)測(cè)到34 864 個(gè)基因,總長(zhǎng)度為27 532 625 bp。在預(yù)測(cè)的基因中,查找到79 130 個(gè)外顯子,總長(zhǎng)度為18 557 406 bp;79 129 個(gè)內(nèi)含子,總長(zhǎng)度為8 975 219 bp。有30 737 個(gè)基因在功能數(shù)據(jù)庫中能比對(duì)到注釋信息,Nr 具有最高的注釋率(30 584,87.72%)KEGG 具有最低的注釋率(9623,27.60%)。在預(yù)測(cè)的基因中,15 184 個(gè)基因被分類為GO 功能類別;15 598 個(gè)基因被分類為KOG 功能類別;9623個(gè)基因注釋到125 個(gè)KEGG 代謝途徑。在GO 功能類別中,包含分子功能、細(xì)胞組成和生物過程(圖5);在KOG 功能類別中,通用功能預(yù)測(cè)的基因最多,其次是翻譯后修飾、蛋白質(zhì)周轉(zhuǎn)、分子伴侶和信號(hào)轉(zhuǎn)導(dǎo)機(jī)制(圖6);在參與基因數(shù)目最多的前10條KEGG 代謝途徑中,注釋基因分別參與核糖體代謝、植物激素信號(hào)轉(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
通過與胡蘿卜、丹參、葡萄和擬南芥等物種蛋白序列比對(duì)(圖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ì)。對(duì)2001 個(gè)單拷貝基因利用ML 構(gòu)建系統(tǒng)發(fā)育樹(圖7-C)。從川芎系統(tǒng)發(fā)育樹分析,川芎、胡蘿卜、丹參、葡萄和擬南芥來自于共同祖先。其中,川芎與胡蘿卜最晚與其他物種發(fā)生分歧,二者進(jìn)化分支長(zhǎng)度差異最小,分歧時(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
阿魏酸是川芎的主要藥效成分,屬于苯丙烷生物合成途徑,是木質(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ù),對(duì)阿魏酸合成過程的相關(guān)基因進(jìn)行了檢索,共有53 個(gè)功能基因覆蓋到阿魏酸合成途徑中(表1)。其中,編碼咖啡酸-O-甲基轉(zhuǎn)移酶、4-香豆酸輔酶A 連接酶、莽草酸O-羥基肉桂酰轉(zhuǎn)移酶、咖啡酰輔酶A-O-甲基轉(zhuǎn)移酶等酶的基因均有多個(gè)同源拷貝存在,預(yù)測(cè)川芎基因組在進(jìn)化的過程中某一時(shí)間點(diǎn)發(fā)生過基因擴(kuò)張。
基因組包含了一個(gè)生物體所有基因的總和,了解生物體基因組信息有助于深入了解生物體遺傳、進(jìn)化、生物合成、次生代謝等全部過程。通過基因組測(cè)序技術(shù)可以對(duì)特定物種基因組進(jìn)行測(cè)序,利用生物信息學(xué)方法對(duì)測(cè)序序列進(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ù)和高通量測(cè)序相結(jié)合的方式對(duì)川芎基因組進(jìn)行測(cè)算。通過流式細(xì)胞術(shù)估算出川芎基因組大小分別為3 897.12 Mb和3 034.28 Mb,通過高通量測(cè)序的K-mer 分析,綜合2 個(gè)分析結(jié)果估算川芎基因組大小為3 058.37 Mb,屬于基因組較大的物種。本研究測(cè)得川芎基因組GC 含量為36.32%,處于植物基因組GC 含量應(yīng)介于25%~65%的合理范圍[27],說明川芎基因組測(cè)序的結(jié)果和組裝是正確可靠的。通過K-mer 分析,估算出川芎重復(fù)序列含量為79.79%,雜合率約為2.16%,與地黃[28]、黃芪[29]等藥用植物類似,呈現(xiàn)高重復(fù)、高雜合的基因組特征,進(jìn)一步說明川芎基因?qū)儆诟咧貜?fù)、高雜合、基因組龐大的物種。
在基因預(yù)測(cè)、注釋和基因家族的鑒定中,本研究共預(yù)測(cè)到川芎編碼蛋白基因34 864 個(gè),遠(yuǎn)高于其傘形科親緣關(guān)系較近的胡蘿卜(32 113 個(gè))[30],可能是因?yàn)榻M裝的都是短片段測(cè)序文庫,川芎中的基因數(shù)量可能被高估了。此外,本研究中測(cè)序完成后從頭組裝中產(chǎn)生的contig N50 為286 bp,scaffold N50 為493 bp,明顯較預(yù)期短,這與廣藿香[4]、羅漢果[12]等藥用植物的全基因組調(diào)研結(jié)果一致。提示對(duì)于川芎這種具有復(fù)雜基因組特征的物種來說,利用二代高通量測(cè)序?qū)ζ淙蚪M進(jìn)行精確測(cè)序仍然存在技術(shù)難度。因此,提高川芎基因組的測(cè)序深度和組裝質(zhì)量,建議后續(xù)的研究可采用二代和三代測(cè)序技術(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ì)測(cè)序奠定基礎(chǔ)。本研究中組裝產(chǎn)生的大量川芎基因組序列和注釋基因?yàn)楹罄m(xù)分子標(biāo)記的開發(fā)和基因功能研究提供了大量的資源。同時(shí),本研究挖掘了阿魏酸合成途徑中的參與基因,對(duì)川芎阿魏酸生物合成途徑潛在分子機(jī)制的初步研究,為進(jìn)一步研究川芎生物學(xué)和選育具有優(yōu)良藥用性狀的品種奠定了基礎(chǔ)。
利益沖突所有作者均聲明不存在利益沖突