基于海甘藍RNA-Seq 序列開發(fā)EST-SSR 分子標記

戚維聰， 程計華， 黃邦全， 李 坦， 林 峰

(1.江蘇省農(nóng)業(yè)科學院農(nóng)業(yè)生物技術研究所，江蘇 南京210014;2.荷蘭瓦赫寧根大學植物育種系，荷蘭 瓦赫寧根6708PB;3.中國科學院武漢植物園，湖北 武漢430074;4.湖北大學生命科學院，湖北 武漢430062)

海甘藍(Crambe abyssinica)是異源六倍體(2n=6x=90)，原產(chǎn)于地中海地區(qū)，屬于十字花科(Brassicaceae)兩節(jié)薺屬(Crambegenus)［1-2］。生產(chǎn)芥酸是海甘藍經(jīng)濟價值的主要體現(xiàn)，海甘藍的芥酸含量達到60%。芥酸(順-13-二十二碳一烯酸)是一種長碳鏈不飽和脂肪酸，在化工中應用廣泛，其衍生物芥酸酰胺是塑料生產(chǎn)的必須添加劑，主要作用是潤滑，防止粘連。除此之外芥酸還可以被用在制藥、生產(chǎn)船舶表面涂料和尼龍等多個方面［1，3］。全球對芥酸年需求量大且增長迅速，從1990 年到2010 年，總消費量從1.8 ×107t 增長至3.5 ×107t，大約翻了一番。近年來由于大面積推廣種植雙低油菜，高芥酸資源出現(xiàn)短缺。為了滿足工業(yè)上對芥酸這一重要化工原料的需求，美國、荷蘭等國相繼培育、釋放、種植了一系列高芥酸含量的海甘藍新品種。海甘藍中度耐鹽，適合在鹽堿等邊際性土地上生長。中國擁有近2.2 ×106hm2海洋灘涂，利用這些邊際性土地種植高芥酸海甘藍將產(chǎn)生顯著經(jīng)濟效益［4-5］。

海甘藍作為一個“新”的經(jīng)濟作物，其遺傳研究較為薄弱。目前，NCBI 數(shù)據(jù)庫只搜集了37 條海甘藍EST 序列，基于海甘藍遺傳信息開發(fā)的分子標記尚未見報道。二代測序技術的推廣和應用為快速有效地發(fā)展“新”經(jīng)濟作物的分子標記提供了契機［6-7］。通過對轉錄組RNA-Seq 數(shù)據(jù)的從頭拼接(De novo)和組裝，海量測序短片段被連接成一類較長的重疊群(Contigs)，這些重疊群包含了特定發(fā)育時期的部分轉錄本信息［8］。轉錄本含有大量的簡單重復序列(Simple sequence repeat，SSR)，由于簡單重復序列變異度較高，這些位點已成為開發(fā)功能分子標記的理想靶位點［8-9］。利用RNA-Seq 組裝的重疊群序列開發(fā)EST-SSR 分子標記的技術已在蓖麻、草坪草、芝麻，紅薯和咖啡等植物中廣泛應用［10-12］。但是由于海甘藍研究中尚未獲得足夠數(shù)量的EST 和RNA-Seq 數(shù)據(jù)信息，開發(fā)EST-SSR 分子標記仍是空白。

海甘藍種子發(fā)育階段是合成芥酸的主要時期，關鍵基因在這段時間內相繼表達?；诖?，我們將對發(fā)育時期種子進行高通量深度測序，利用從頭拼接軟件獲得重疊群，對這些重疊群序列進行SSR 位點掃描，并對海甘藍基因組DNA 進行PCR 擴增，以期為揭示海甘藍遺傳背景的多樣性以及遺傳連鎖圖譜的構建提供第一代SSR 分子標記。

1 材料與方法

1.1 供試材料及全基因組測序

海甘藍品種為C.abyssinicacv.Galatic。去殼種子置于鋪有3 層潤濕濾紙的培養(yǎng)皿中，發(fā)芽后移入溫室土壤栽培，待植物開花后掛牌標記。取開花后20 d 的種子提取mRNA(10 粒為一組)，試劑盒為iScriptTMcDNA Synthesis Kit (Bio-rad，USA)。DNA提取材料種植條件如前所述，取0.5 g 新鮮葉片提取基因組DNA。

1.2 組裝拼接

采用Illumina HiSeq 2000 測序儀進行測序。進一步利用SolexaQA［13］軟件包對原始測序數(shù)據(jù)進行質量(Q20，Phred-Score≥20 即1%的錯誤率)和測序長度(L40，長度≥40 bp)過濾。隨后，利用Trinity轉錄組組裝軟件［14］按默認參數(shù)對清理后的序列進行轉錄組從頭de novo組裝。

1.3 SSR 全基因組挖掘及引物設計

對Trinity 程序組裝獲得的轉錄組重疊群序列，用MISA 程序(http://pgrc.ipk-gatersleben.de/misa/misa.html)掃描轉錄組SSR 位點(http://pgrc.ipkgatersleben.de/misa/)，滿足單堿基重復超過10 次或雙堿基重復超過8 次或三堿基重復超過5 次或四堿基重復超過4 次或五堿基重復超過2 次或六堿基重復超過3 次被定義為SSR 位點。通過Perl 語言抽取SSR 位點側翼序列后利用Primer 3(http://biotools.umassmed.edu/bioapps/primer3_www.cgi)進行PCR 特異引物設計，退火溫度為(58 ±3)℃。

1.4 海甘藍SSR 位點在白菜基因組的定位及注釋

為了進一步揭示海甘藍EST-SSRs 相關的基因組信息及變異效應，利用所開發(fā)的SSRs 分子標記保守側翼序列，結合白菜基因組序列信息(http://www.ncbi.nlm.nih.gov/nuccore/AENI00000000)，通過Blast［14］同源比對，去除冗余以及進行全基因組物理定位。SSR 位點多態(tài)性比較則主要通過電子ePCR 軟件(http://www.ncbi.nlm.nih.gov/projects/e-pcr/)完成，利用本地數(shù)據(jù)庫預測白菜基因組中PCR 擴增產(chǎn)物的大小。

1.5 PCR 擴增和電泳

PCR 反應總體系為25.0 μl，包括10 ×Buffer(含Mg2+)2.5 μl，dNTP 2.5 μl，Taq酶0.1 μl，引物2.0 μl，模板DNA 2.0 μl，ddH2O 15.9 μl，上覆20.0 μl礦物油。PCR 反應程序為:94 ℃預變性3 min;94℃變性30 s，58 ℃退火1 min，72 ℃延伸80 s，共進行38 個循環(huán)。用2%瓊脂糖膠加壓100 V 進行電泳，紫外透射儀上觀察、拍照。

2 結果與分析

2.1 組裝拼接

利用Illumina Paired-End 100 bp×2 雙端測序，經(jīng)Q20 質量過濾和L40 長度清理，獲得了86 224 256 條平均長度約93 bp，共計7 999 663 632個堿基的重疊群。數(shù)據(jù)進一步利用Trinity 組裝程序進行從頭de novo組裝，共獲得長度等于200 個堿基的重疊群序列234 622條，重疊群平均長度為956 bp，其中N50 為1 428 bp，最長的重疊群為16 475 bp。由于組裝后的重疊群存在選擇性剪切產(chǎn)生的冗余，利用Blastn (1e-50)去除高相似度的序列后得到186 778條重疊群，這些重疊群序列用于后續(xù)的SSR 挖掘。

2.2 EST-SSR 位點掃描

MISA 軟件掃描重疊群后共發(fā)現(xiàn)在38 601條序列中含有47 073個SSR 位點，其中有6 816條重疊群含有1 個以上的SSR 位點，平均每4 500 bp 含有1 個SSR位點。除了單堿基重復10 次的SSR 位點外(45%)，三堿基重復SSR 位點占了總數(shù)的29%，雙堿基重復占10%，五堿基重復占7%，六堿基重復占5%，四堿基重復占2%。三堿基重復中AAG/CTT 基序存在于4 993條序列中，豐度最高(表1、表2、圖1)。

表1 海甘藍EST-SSR 搜索結果Table 1 Summary of EST-SSR searching resulting results

表2 海甘藍EST-SSR 位點長度信息Table 2 Length distribution of EST-SSR based on the number of repeat units

2.3 EST-SSR 引物設計及過濾

利用SSR 位點側翼保守序列，共設計了16 355對SSR 引物。PCR 產(chǎn)物擴增范圍為150 ～400 bp。利用Blastn 對電子PCR 產(chǎn)物進行了冗余性評估，發(fā)現(xiàn)只有6 639條引物可以專化擴增出唯一特異性條帶。利用已測序白菜基因組信息對6 639對引物進行篩選，在5'端允許3 個堿基錯配，3'端不允許錯配，PCR 產(chǎn)物長度大于200 bp，小于500 bp 的條件下，獲得了1 206條物種間可通用的EST-SSR 引物。根據(jù)白菜基因組每10 kb DNA 區(qū)間內只允許存在一個分子標記的條件下，共獲得了688 個分子標記。圖2 列出了這些分子標記在白菜基因組的物理位置。除個別位置外，ESTSSR 分子標記的位置在染色體上基本呈均勻分布。

圖1 三核苷酸不同基序排列分布Fig.1 Frequency distribution of trinucleotide EST-derived SSRs based on motif sequences

2.4 EST-SSR 分子標記驗證

基于白菜基因組信息，在每條染色體臂兩端各挑選了一個SSR 位點，共計20 對EST-SSR 引物在海甘藍中進行引物驗證(表2)。通過瓊脂糖凝膠電泳發(fā)現(xiàn)，在檢測的20 對SSR 引物中，擴增成功率達100%。20 對PCR 產(chǎn)物片段大小均與預測結果完全吻合(圖3)。

3 討論

通過對已公開EST 序列進行SSR 位點掃描開發(fā)分子標記的方法已被廣泛地運用。海甘藍作為一個“新”作物，其EST 信息在NCBI 數(shù)據(jù)庫中幾乎為空白。我們通過二代測序技術，首次獲得了大量海甘藍EST 序列。其中約有3%左右的序列含有SRR位點，這些位點側翼序列為設計引物提供了模板。雖然通過以上的方法可以快速高效設計EST-SSR引物，但是在研究過程中發(fā)現(xiàn)，二代測序拼接結果中含有大量冗余(Redundancy)序列，基于這些冗余序列開發(fā)的一系列EST-SSR 引物只能被認定為一個位點的分子標記。導致大量冗余序列的主要原因是RNA 選擇性剪切和測序過程中的污染。因此，在設計海甘藍EST-SSR 引物前，我們首先對Trinity 拼接結果進行了過濾，通過序列之間相互對比，去除了30%高相似度的重疊群。其次對余下的序列進行了功能注釋，刪除明顯來自人類和細菌的污染序列。最后，利用海甘藍近緣物種白菜基因組序列，淘汰了物理位置在10 kb 之內的分子標記。通過以上3 種過濾，本研究最大程度地獲得了688 個全基因組較為分散的SSR 分子標記。

由于EST 主要編碼功能基因，受選擇壓力大，序列保守性高，導致EST-SSR 所揭示的多態(tài)性相對基因組SSR 低［13，15-17］。研究表明，在EST 的3'區(qū)存在較高的三核苷酸變異率。如控制人類亨丁頓舞蹈癥的HTT基因是在第四對染色體，HTT基因重復排序數(shù)目CAG 太多時會造成亨丁頓舞蹈癥，重復數(shù)28 次以下沒有亨丁頓舞蹈癥表現(xiàn)，而40 以上有亨丁頓舞蹈癥，基于CAG 重復序列開發(fā)的SSR 標記在人類不同個體間具有較高的多態(tài)性［9］。玉米Prolamin 蛋白結合因子(PBF)在3'區(qū)也存在一個高變異率的ACC 三核苷酸序列，基于該重復序列開發(fā)的umc1065 EST-SRR 分子標記在玉米不同自交系間揭示極高的多態(tài)性［18］。盡管根據(jù)EST 設計的SSR 較根據(jù)基因組設計的多態(tài)性偏低，但由于與功能基因緊密相連，可被開發(fā)成功能性分子標記，顯著提高分子標記對目的功能基因的選擇效率。

EST-SSR 引物較根據(jù)基因組信息設計的SSR 引物擴增成功效率偏低，原因主要一方面引物結合位點落在外顯子與外顯子的連接區(qū)，該區(qū)段在基因組上被內含子隔開，導致引物無特定結合區(qū)域;另一方面，PCR 擴增區(qū)域若含有較長的內含子，引物在基因組DNA 為模板時由于區(qū)段太長，無法獲得有效的擴增。由于缺乏海甘藍基因組信息，無法僅通過生物信息學手段預測內含子的位置和大小。為了最大程度解決EST-SSR 引物擴增效率偏低的問題，我們利用海甘藍的近緣種白菜的基因組測序信息來過濾外顯子連接區(qū)和大內含子區(qū)，最終獲得了近688 個高質量的引物信息。這些白菜和海甘藍通用分子標記(Transferable EST-SSR markers)可以被用于這兩個物種間的比較基因組研究，利用白菜的遺傳信息加速海甘藍關鍵基因的鑒定和克隆。

注:本研究獲得的EST-SSR 引物及對應的具體信息可向作者索取。

［1］ KANEKO Y，BANG S，MATSUZAWA Y.Wild crop relatives:genomic and breeding resources［M］.Berlin:Springer，2011:247-258.

［2］ WARWICK S，GUGEL R.Genetic variation in theCrambe abyssinica-C.hispanica-C.glabratacomplex［J］.Genetic Resources and Crop Evolution，2003，50(3):291-305.

圖2 海甘藍EST-SSR 分子標記在白菜基因組物理分布Fig.2 The physical localization of EST-SSR on Brassica napa genome

圖3 20 對SSR 引物基因組擴增帶型Fig.3 PCR banding patterns of Crambe abyssinice cv.Galatic by 20 SSR primer pairs

［3］ ABADA Y S K，NGUYEN H P，SCHREIBER R，et al.Assessment of motor function，sensory motor gating and recognition memory in a novel BACHD transgenic rat model for Huntington disease［J］.PLoS ONE，2013，8(7):e68584.

［4］ METTERNICHT G，ZINCK J.Remote sensing of soil salinity:potentials and constraints ［J］.Remote Sensing of Environment，2003，85(1):1-20.

［5］ RHOADES J，MANTEGHI N，SHOUSE P，et al.Soil electrical conductivity and soil salinity:new formulations and calibrations［J］.Soil Science Society of America Journal，1989，53(2):433-439.

［6］ VARSHNEY R K，NAYAK S N，MAY G D，et al.Next-generation sequencing technologies and their implications for crop genetics and breeding［J］.Trends in Biotechnology，2009，27(9):522-530.

［7］ DAVEY J W，HOHENLOHE P A，ETTER P D，et al.Genomewide genetic marker discovery and genotyping using next-generation sequencing［J］.Nature Reviews Genetics，2011，12(7):499-510.

［8］ RAHMAN M，SUN Z，MCVETTY P B.High throughput genomespecific and gene-specific molecular markers for erucic acid genes inBrassica napus(L.)for marker-assisted selection in plant breeding［J］.Theoretical and Applied Genetics，2008，117(6):895-904.

［9］ RIESS O，NOERREMOELLE A，SOERENSEN S A，et al.Improved PCR conditions for the stretch of (CAG)n repeats causing Huntington's disease［J］.Human Molecular Genetics，1993，2(6):637.

［10］ AGGARWAL R K，HENDRE P S，VARSHNEY R K，et al.Identification，characterization and utilization of EST-derived genic microsatellite markers for genome analyses of coffee and related species［J］.Theoretical and Applied Genetics，2007，114(2):359-372.

［11］ FEINGOLD S，LLOYD J，NORERO N，et al.Mapping and characterization of new EST-derived microsatellites for potato (Solanum tuberosumL.)［J］.Theoretical and Applied Genetics，2005，111(3):456-466.

［12］ ZHANG T，GE M，YE X，et al.Construction of a linkage map for quantitative trait loci associated with economically important traits in creeping bentgrass (Agrostis stoloniferaL.)［J］.Euphytica，2012，188(3):347-360.

［13］ GORBACH D M，HU Z L，DU Z Q，et al.Mining ESTs to determine the usefulness of SNPs across shrimp species［J］.Animal Biotechnology，2010，21(2):100-103.

［14］ ALTSCHUL S F，MADDEN T L，SCH?FFER A A，et al.Gapped BLAST and PSI-BLAST:a new generation of protein database search programs［J］.Nucleic Acids Research，1997，25(17):3389-3402.

［15］ MUCHERO W，DIOP N N，BHAT P R，et al.A consensus genetic map of cowpea［Vigna unguiculata(L.)Walp.］and synteny based on EST-derived SNPs［J］.Proceedings of the National Academy of Sciences，2009，106(43):18159-18164.

［16］ SEHGAL D，RAJARAM V，ARMSTEAD I P，et al.Integration of gene-based markers in a pearl millet genetic map for identification of candidate genes underlying drought tolerance quantitative trait loci［J］.BMC Plant Biology，2012，12(1):9.

［17］ XUE S，ZHANG Z，LIN F，et al.A high-density intervarietal map of the wheat genome enriched with markers derived from expressed sequence tags ［J］.Theoretical and Applied Genetics，2008，117(2):181-189.

［18］ VICENTE-CARBAJOSA J，MOOSE S P，PARSONS R L，et al.A maize zinc-finger protein binds the prolamin box in zein gene promoters and interacts with the basic leucine zipper transcriptional activator Opaque2［J］.Proceedings of the National Academy of Sciences，1997，94(14):7685-7690.