基于簡化基因組測序分析云南小桐子遺傳多樣性

鄭周宜　劉雨婷　陳春　王莎莎　唐唯

關(guān)鍵詞：小桐子；遺傳多樣性；分子標(biāo)記；全基因組關(guān)聯(lián)分析；單核苷酸多態(tài)性

中圖分類號：S565.9 文獻(xiàn)標(biāo)識碼：A

小桐子（Jatropha curcas L.）屬大戟科麻瘋樹屬的多年生木本油料植物，原產(chǎn)美洲，現(xiàn)廣泛分布于全世界熱帶及亞熱帶地區(qū)，在我國尤以云南分布面積最廣、資源數(shù)量最大[1-2]。小桐子具有生物產(chǎn)量大、適應(yīng)范圍廣、耐貧瘠等特性，是維持生態(tài)[3]、醫(yī)藥開發(fā)、生物防治及生產(chǎn)有機(jī)肥料的理想資源[4-5]。由于其種子含油率高達(dá)60%以上，也被認(rèn)為是具有巨大開發(fā)潛力的多用途能源植物樹種[6]。盡管近年來國內(nèi)外有關(guān)小桐子的研究報道數(shù)量大增，但大多集中在藥理毒理活性及新藥開發(fā)[4-5]、生理適應(yīng)機(jī)制[7-8]、生物柴油產(chǎn)業(yè)化[6]、經(jīng)濟(jì)用途開發(fā)[3]等方面，關(guān)于其種質(zhì)資源的遺傳多樣性研究起步相對較晚。另一方面，種質(zhì)資源不清、品種良莠不齊、基因資源庫保存不完善等問題已成為限制國內(nèi)小桐子產(chǎn)業(yè)開發(fā)的重要瓶頸[9]。

迄今為止，來自中國、巴西、印度、墨西哥等多國的研究人員已先后采用RAPD、AFLP 等多種分子標(biāo)記對國內(nèi)外小桐子種質(zhì)進(jìn)行了遺傳多樣性的鑒定和評價[10-11]，發(fā)現(xiàn)我國范圍內(nèi)收集的小桐子種質(zhì)資源遺傳背景雖普遍低于它國（南美、東南亞地區(qū)等）[12-13]，但云南地區(qū)的小桐子居群具有較為豐富的遺傳多樣性，且其居群內(nèi)的遺傳分化高于居群間[14-16]。上述研究在一定程度上印證了向振勇等[17]提出的云南小桐子種質(zhì)來源較為復(fù)雜的推測。

基于上述溯源，為了進(jìn)一步了解近年云南地區(qū)小桐子種群結(jié)構(gòu)特點，本研究首先從種子和幼苗性狀特點入手，結(jié)合簡化基因組測序和分子標(biāo)記技術(shù)解析云南各地區(qū)小桐子的遺傳多樣性特點；利用GWAS 分析方法，挖掘地下部分鮮重這一重要性狀差異基因及地區(qū)獨(dú)特性基因，預(yù)期研究結(jié)果可為小桐子種質(zhì)資源的保護(hù)利用、遺傳改良以及雜交育種親本選配提供重要遺傳基礎(chǔ)。

1 材料與方法

1.1 材料

供試材料來自云南省不同地區(qū)的9 個居群，種源地信息見表1。每個地區(qū)根據(jù)植株分布情況隨機(jī)選取50 個具有代表性的樣株，每個樣株至少間隔5 m 以上，每個樣株分別取3 顆種子，編號整理。

1.2 方法

1.2.1 幼苗培養(yǎng) 在各種源地的種子材料中挑選顆粒飽滿、大小均一的種子50 粒，清洗凈后以1% CuSO4 溶液浸泡消毒20 min，經(jīng)蒸餾水充分漂洗、浸泡、吸漲后于26℃萌發(fā)。將萌發(fā)14 d 后的種子移栽至裝有等量營養(yǎng)土的花盆中，置于人工氣候箱中繼續(xù)培養(yǎng)14 d，用于幼苗生物學(xué)性狀的測定。

1.2.2 性狀測定 （1）百粒重及種子形狀測定。

在各種源地隨機(jī)選取種子，取3 個重復(fù)，每個重復(fù)30 粒，以電子天平分別稱重（精確至0.01 g），計算種子百粒重。以游標(biāo)卡尺測定長徑、寬徑及側(cè)徑（精確至0.01 mm）。

（2）生長指標(biāo)及生物量測定。每個種源隨機(jī)抽樣10 株苗，于移栽后14 d 進(jìn)行株高測定（精確至0.01 cm）。生物量測定采用全稱重法。樣株在流動自來水下充分洗凈、吸干表面水分，分別對地上部分、地下部分進(jìn)行鮮重稱量。

（3）數(shù)據(jù)處理。以SPSS 軟件（v26.0）進(jìn)行數(shù)據(jù)統(tǒng)計、方差分析、多重比較（Duncans 法或非參數(shù)檢驗，以不同小寫英文字母表示顯著差異，P<0.05）及相關(guān)性分析。

利用變異系數(shù)（CV）分析種子性狀變異性，計算公式為：CV=（S/X）×100%（S 為標(biāo)準(zhǔn)差，X 為平均值）。計算根冠比=地下部分（鮮重）/地上部分（鮮重）×100%。

1.2.3 RAD-seq 基因組測序及分子標(biāo)記開發(fā) （1）基因組DNA 的提取。綜合考慮幼苗性狀（鮮重、根冠比）及平均變異系數(shù)，取差異顯著的2 個居群材料（各10 株）為BSA 分離群體，其中版納勐海BNMH 為BSA1，楚雄元謀CXYM 為BSA2，采用CTAB 法提取基因組DNA（以第2 片真葉為材料），經(jīng)瓊脂糖凝膠電泳（0.8%～1.0%）和超微量分光光度計檢測質(zhì)量、純度和完整性。

（2）基因組測序及分子標(biāo)記挖掘。構(gòu)建插入片段350 bp 的文庫，測序在ILLUMINA Hiseq2000 平臺進(jìn)行（Novogenes Inc.），測序讀長為雙端150 bp，測序序列用GenomeScope 2.0[18]進(jìn)行K-mer 分析確定染色體倍性。SNP 標(biāo)記挖掘由CLC Genomics Workbench 12.0.1（QIAGEN， Inc.）完成，工具為Basic Variation Detection（參數(shù)默認(rèn)）。

（3）SSR 引物設(shè)計。取上述2 個居群的小桐子材料進(jìn)行RAD-seq 基因組測序（5×覆蓋度），以測序基因組為基礎(chǔ)，利用基于Perl 的MISA 程序挖掘SSR 引物[19][（X）n；（XX）n-2；（XXX）n-3；（XXXX）n-4；（XXXXX）n-6；（XXXXXX）n-7，n=10] ， 挖掘到的SSR 標(biāo)記以軟件primer3（https：//primer3.org/）設(shè)計引物。引物位于SSR標(biāo)記2 個側(cè)翼，長度在20～22 bp 之間，正向引物的5端加上熒光標(biāo)記FAM，PCR 目標(biāo)產(chǎn)物預(yù)期大小在50～150 bp 范圍內(nèi)。引物用QC filter 進(jìn)行篩選，去掉純合位點，引物由昆明擎科生物科技有限公司合成。

（4）分子標(biāo)記的電泳檢測（瓊脂糖凝膠、毛細(xì)管電泳）及遺傳多樣性分析。對檢測和純化好的不同種群小桐子基因組DNA 進(jìn)行篩選，配制反應(yīng)體系后進(jìn)行PCR 擴(kuò)增，Tm 以引物合成報告單為準(zhǔn)。對瓊脂糖凝膠電泳初步確定擴(kuò)增效果理想的引物進(jìn)行毛細(xì)管電泳分析（ABI 3500），用GeneMapper 5.0（https：//tools. thermofisher.com/content/sfs/manuals/4476603A.pdf）讀取條帶，用MEGA 5.0[20] 構(gòu)建UPGMA 系統(tǒng)發(fā)育樹， 用PopGene 1.31（https：//sites. ualberta.ca/～fyeh/popgene.pdf）進(jìn)行遺傳多樣性計算。

（5）基于重要性狀差異的GWAS 分析。以地下部分鮮重統(tǒng)計數(shù)據(jù)為性狀表型，以全基因組SNPs 為基因型，使用R 包GAPIT 混合線性模型（MLM）進(jìn)行GWAS 分析（全基因組關(guān)聯(lián)分析，Genome Wide Association Study），以CMplot 包進(jìn)行GWAS 結(jié)果輸出，分別生成SNPs 密度分布圖和柱狀曼哈頓圖。

（6）基因組注釋及比較基因組分析。對上述重測序2 份材料的clean data 進(jìn)行全基因組比對（CLC Genomics Workbench 12.0.1）得到共線性結(jié)果。對2 個基因組進(jìn)行de novo 組裝（SOAPdenovo2[21]），得到的組裝序列用Augustus（http：//bioinf.uni-greifswald.de/augustus/）進(jìn)行CDS 預(yù)測。對預(yù)測到的基因以NCBI 中protein databank 和swiss-prot protein databank 進(jìn)行同源注釋和比較分析，找出同時具有以下2 個特征的特異蛋白：①僅在上述2 份材料中存在的；②與地下部分鮮重性狀緊密連鎖。

2 結(jié)果與分析

2.1 云南地區(qū)不同種源小桐子種子性狀及其變異分析

百粒重測量結(jié)果（表1）顯示其變幅在46.06～80.40 g，不同種源的差異可達(dá)顯著差異水平（P<0.05）。長徑、寬徑、側(cè)徑變幅依次為13.00～19.50、7.50～12.00、5.20～10.10 mm，最大平均值均為雙江大文鄉(xiāng)種源，最小平均值為麗江華坪種源。對種子性狀的變異性分析發(fā)現(xiàn)，供試種子變異系數(shù)差別較為明顯， 其中百粒重變異系數(shù)最大（13.93%），長徑、寬徑與側(cè)徑的變異系數(shù)依次為3.04%、5.34%和6.20%。

2.2 云南地區(qū)不同種源小桐子幼苗性狀及其變異分析

從表1 可以看出，株高性狀表現(xiàn)最好的是臨滄耿馬種源，平均苗高為11.46 cm，版納勐海種源具有相對最高的生物量積累（ 全株鮮重6.60 g），麗江華坪、楚雄祿豐和楚雄元謀種源在苗高和生物量積累上均表現(xiàn)最差，盡管楚雄元謀種源的根冠比表現(xiàn)最佳（0.25）。方差分析表明不同種源的苗高、地上部分鮮重、全株鮮重及根冠比均存在顯著差異。不同種源地小桐子表型性狀變異程度依次為：地上部分鮮重＞全株鮮重＞地下部分鮮重＞根冠比＞株高＞百粒重＞種子側(cè)徑＞種子寬徑＞種子長徑（表1），尤以生物量積累的性狀變異系數(shù)相對較高，暗示這些表型性狀存在遺傳不穩(wěn)定性，易受環(huán)境壓力的選擇。

2.3 RAD-seq 基因組測序及結(jié)果分析

測序結(jié)果顯示BNMH（BSA1）及CXYM（BSA2）的染色體為2 倍體（圖1A），得到具有EcoRⅠ酶切位點的Reads 數(shù)分別為64 734 412 條和60 187 766 條。平均Q20 為98.6%，平均Q30為93.9%，對reads 進(jìn)行有參組裝（參考基因組RJC1）后分別得到2783 和3033 條高質(zhì)量的contigs（表2）。

2.4 SSR 引物開發(fā)及遺傳多樣性分析

利用MISA 隨機(jī)開發(fā)出8 對SSR 引物，PCR擴(kuò)增后利用瓊脂糖凝膠電泳確認(rèn)有效擴(kuò)增引物為4 對（表3）；以4 對引物、9 個居群共95 個DNA樣本進(jìn)行毛細(xì)管電泳檢測，結(jié)果發(fā)現(xiàn)引物多態(tài)性擴(kuò)增條帶最多為6 條，最少為2 條，各居群間利用4 個毛細(xì)管電泳引物擴(kuò)增遺傳多樣性水平均較高（表4）。以每個地區(qū)所有個體多態(tài)性條帶為居群進(jìn)行UPGMA 聚類分析，結(jié)果表明保山地區(qū)和西雙版納地區(qū)可能作為云南省小桐子的親本或起源，各居群間親緣關(guān)系和地理位置有一定相關(guān)性（圖1B）。

2.5 基于地下部分鮮重差異的GWAS分析

基于在種子及幼苗性狀調(diào)查分析中發(fā)現(xiàn)，地下部分鮮重（underground fresh weight， UFW）為9 個居群間差異最大，因此本研究針對UFW，首先在全基因組水平挖掘到相關(guān)InDel 標(biāo)記8182個，SNPs 標(biāo)記82 155 個。因小桐子為雜合二倍體，因此經(jīng)純合和等位基因頻率（≥50%）過濾后，共挖掘到有效SNPs 22 756 個。GWAS 分析進(jìn)一步發(fā)現(xiàn)，SNPs 分布于11 對染色體（圖2A），UFW 高頻SNPs 集中于9 號染色體（圖2B）。

2.6 注釋及比較基因組分析

本研究通過對基因組測序材料BNMH（版納勐海）和CXYM（楚雄元謀）進(jìn)行注釋和功能預(yù)測，分別得到27 133 和26 282 個候選功能蛋白。Blast2GO 分析對預(yù)測蛋白進(jìn)行GO 功能注釋，發(fā)現(xiàn)在生物學(xué)途徑（biological process）中，獲得注釋條目最多的是蛋白修飾、跨膜運(yùn)輸和轉(zhuǎn)錄調(diào)節(jié)；在分子功能（molecular function）和細(xì)胞組件（cellular component）中分別為催化活性、結(jié)合功能、膜整合組分和胞質(zhì)（圖3）。

為了初步明確云南本地小桐子的獨(dú)特性，進(jìn)一步比較分析發(fā)現(xiàn)，僅在BNMH 和CXYM 兩份測序材料中存在的候選蛋白共40個，其中未注釋的功能蛋白有16 個，和UFW 性狀連鎖的有3 個（表5）。

3 討論

小桐子（J. curcas L.）是一種在全世界（亞）熱帶地區(qū)廣泛種植的多用途油料作物，其作為一種新興能源植物的開發(fā)潛力而引起了全世界的廣泛關(guān)注，然而，現(xiàn)有的育種種質(zhì)在遺傳多樣性上的匱乏成為小桐子完全商業(yè)化的主要阻礙，而種質(zhì)資源的收集和評估對于培育高產(chǎn)、優(yōu)質(zhì)和高適應(yīng)性的小桐子新品種至關(guān)重要。迄今為止在研究多樣性的各種分子標(biāo)記中，基于PCR 的標(biāo)記如RAPD（隨機(jī)擴(kuò)增多態(tài)性DNA）、ISSR（簡單序列間重復(fù)）、SSR（簡單序列重復(fù)）、EST-SSR 和AFLP（擴(kuò)增片段長度多態(tài)性）均已被用于探究小桐子的遺傳多樣性和系統(tǒng)發(fā)生。PECINAQUINTERO等[22]借助AFLP 技術(shù)對來自墨西哥中部和東南部9 個州的175 份種質(zhì)（包括有毒性、無毒性、蛋白質(zhì)含量和種油含量對比大的基因型）展開全面的遺傳多樣性分析，發(fā)現(xiàn)墨西哥小桐子具有很高的遺傳多樣性（明顯高于GRATIVOL 等[23]對巴西來源小桐子的相關(guān)報道），并伴隨著不同州的種群結(jié)構(gòu)。SUDHEER 等[24]使用RAPD 和AFLP技術(shù)對42 個印度小桐子種質(zhì)的分子多樣性分析，指出其多態(tài)性的平均百分比為26.47，顯著低于墨西哥（33.18）。這些結(jié)果都支持了目前的主流觀點，墨西哥和中美洲可能是麻瘋樹起源和多樣化中心[23-24]，而非洲和亞洲的小桐子種質(zhì)存在低的遺傳變異性[11， 25-26]。WEN 等[13]早期利用36 個EST-SSR 和20 個G-SSR 從來自南美洲、中國和印度尼西亞的45 個種質(zhì)中檢測到183 個多態(tài)性等位基因，平均遺傳多樣性指數(shù)為0.5572，發(fā)現(xiàn)中國云南的種質(zhì)依然顯示出比其他地理區(qū)域（如格林納達(dá)和中國海南）更高水平的遺傳變異，可作為豐富小桐子遺傳背景的重要種質(zhì)資源。在本研究中，小桐子種子采集自云南省的不同地域，并首先用于多項形態(tài)學(xué)或農(nóng)藝性狀分析，包括種子徑圍、百粒重、幼苗株高和生物量積累等。這些性狀對作物育種改良研究很重要，但應(yīng)用于小桐子遺傳多樣性的系統(tǒng)研究還很少見到報道。本研究還發(fā)現(xiàn)小桐子種群內(nèi)存在顯著性狀差異，尤以生物量相關(guān)性狀的變異系數(shù)最大（ 21.90%～25.87% ）， 這與PESTANA- CALDAS 等[27] 和LAVIOLA 等[28]的研究相似，均認(rèn)為單株植物種子產(chǎn)量（包括種子百粒重）、植株高度對小桐子遺傳變異性貢獻(xiàn)最小，而生長有關(guān)的特征則表現(xiàn)出最大的表型變異，這可能和地區(qū)適應(yīng)性和馴化有密切關(guān)系。

有研究表明在小桐子中，分子標(biāo)記和性狀數(shù)據(jù)連鎖不緊密[28-29]。經(jīng)典的遺傳定位依靠基因型和表型的緊密連鎖[30]，結(jié)合基因組測序的BSA 分析法擴(kuò)大了群體內(nèi)一致性數(shù)據(jù)和群體間差異性數(shù)據(jù)[31]，本研究選擇地下鮮重這一極端差異性狀的2 個分離群體，BNMH（版納勐海）和CXYM（楚雄元謀）以BSA 分離法進(jìn)行全基因組測序，分別組裝獲得2783 和3033 條contigs 并通過GWAS對重要SNP 標(biāo)記的基因組區(qū)域進(jìn)行高精度的關(guān)聯(lián)和標(biāo)記[32]，在9 號染色體上確定了5 個UFW 高頻SNPs 的峰值區(qū)域，有3 個候選蛋白與UFW 性狀連鎖。從基因組測序和比較分析來看，云南省小桐子具有較高的單核苷酸多態(tài)性，具有特有蛋白40 個（含未注釋功能蛋白16 個），顯示了較好的多樣性、獨(dú)特性和地區(qū)適應(yīng)性。此外，地下鮮重主要涉及根系生長和發(fā)育，這一性狀的差異也明顯和地區(qū)馴化有關(guān)。

本研究利用MISA 開發(fā)并隨機(jī)選出8 對SSR引物（有效擴(kuò)增引物4 對），對9 個居群共95 個DNA 樣本進(jìn)行毛細(xì)管電泳檢測，結(jié)果發(fā)現(xiàn)引物多態(tài)性條帶最多為6 個，最少為2 個，引物位點的多態(tài)信息含量指數(shù)（PIC）變幅為0.041～0.711（平均為0.370），高于V?SQUEZ-MAYORGA 等[33]（PIC=0.274±0.165）和KUMAR 等[34（] PIC=0.25±0.16）的報道，可作為小桐子種質(zhì)親緣關(guān)系鑒定、遺傳變異檢測等分子生物學(xué)研究的重要可用資源。此外，群體結(jié)構(gòu)研究雖表現(xiàn)為中度多態(tài)性[35]，但本研究發(fā)現(xiàn)云南小桐子的雜合度（ 平均He=0.437）與目前世界范圍內(nèi)的小桐子報道相比十分可觀，明顯高于WEN 等[13]對南美洲和云南種質(zhì)資源的報道（He 分別為0.33、0.3473）。多數(shù)研究認(rèn)為小桐子起源于中美洲，但傳入我國時間和路線不詳[36]。本研究基于UPGMA 聚類分析發(fā)現(xiàn)各居群間親緣關(guān)系和地理位置有一定相關(guān)性，各居群間遺傳多樣性較高，這與WEN 等[13]、RAFII等[37]的結(jié)論相一致，暗示小桐子可能由多點多分化種傳入云南，傳入后在各生長地之間基因交流較少，變異主要為獨(dú)立群體演化事件，其主要結(jié)果為充分了解小桐子在中國的起源中心的遺傳多樣性和種群結(jié)構(gòu)提供參考，也為其作為可持續(xù)能源作物的復(fù)興奠定新的理論基礎(chǔ)。在后續(xù)研究中應(yīng)進(jìn)一步結(jié)合胞質(zhì)基因組多態(tài)性特點及相關(guān)歷史考據(jù)文獻(xiàn)對小桐子在云南的傳入、傳播及分布特點進(jìn)行更加深入的探討。此外，小桐子作為云南省的一個主要的外來物種，在將來的育種工作中有必要加大對云南?。ㄓ绕涫沁吘车貐^(qū)）的小桐子種質(zhì)資源的勘探與發(fā)掘，結(jié)合InDel 等大片段標(biāo)記[34， 38]，開發(fā)適用于種質(zhì)資源檢測的各類有效標(biāo)記，進(jìn)一步通過分子設(shè)計育種加速小桐子的種植資源改良。

4 結(jié)論

本研究采集了云南省10 個小桐子聚居地種子，調(diào)查了9 個重要農(nóng)藝性狀，利用地下鮮重差異顯著的2 個BSA 分離群體，通過RAD-seq 測序結(jié)合GWAS 分析挖掘到SNPs 22 756 個，分析發(fā)現(xiàn)高頻SNPs 集中在9 號染色體，進(jìn)一步定位到來自于版納和楚雄群體的特有候選蛋白40 個，篩選了3 個UFW 候選蛋白。此外，研究基于4對SSR 標(biāo)記對9 個地區(qū)小桐子種質(zhì)資源進(jìn)行了群體結(jié)構(gòu)分析，得到了較好的多態(tài)性和推測了可能的起源中心。