基于ddGBS的桉樹SNP挖掘和系統(tǒng)進(jìn)化分析

賈翠蓉，朱顯亮，王 莉，周長品，翁啟杰，甘四明，李發(fā)根

（1.中國林業(yè)科學(xué)研究院 a.熱帶林業(yè)研究所；b.熱帶林業(yè)研究國家林業(yè)和草原局重點(diǎn)實(shí)驗(yàn)室，廣東 廣州 510520；2.南京林業(yè)大學(xué) 林學(xué)院，江蘇 南京 210037）

桉樹隸屬于桃金娘科Myrtaceae，主要包括三大屬：桉屬Eucalyptus、杯果木屬Angophora和傘房屬Corymbia，主要分布于澳大利亞，有808 個(gè)種和137 個(gè)變種[1]。大部分桉樹具有雄蕊先熟、花期不遇和地理隔離的特征，這使其自交或近交概率極低，但在重疊鄰接的桉樹自然分布區(qū)或引種栽培區(qū)有可能產(chǎn)生自然或人工雜種[2-4]。桉樹經(jīng)過世代的遺傳變異和自然選擇，逐步演化出種類龐雜的桉樹體系[5]，因此，需要借助分子標(biāo)記的手段來研究其系統(tǒng)進(jìn)化歷程。

單核苷酸多態(tài)性（Single nucleotide polymorphism, SNP）由于其在基因組上分布廣泛、易于高通量檢測和分型等優(yōu)點(diǎn)廣泛應(yīng)用于資源評(píng)價(jià)和遺傳圖譜構(gòu)建等研究[6-7]。隨著高通量測序技術(shù)的發(fā)展，測序成本大幅降低。對(duì)于群體遺傳研究而言，僅需要分布均勻且一定數(shù)量的SNP 標(biāo)記即可，基于測序的基因分型（Genotyping by sequencing,GBS）技術(shù)成本低、操作簡便和高通量等特點(diǎn)使其成為SNP 挖掘的有效方法[8]。ddGBS（Double digest GBS）主要利用雙酶切降低基因組的復(fù)雜度，使用標(biāo)簽序列區(qū)別樣本，實(shí)現(xiàn)了大量樣本的測序分型[9]，廣泛應(yīng)用于種質(zhì)資源的遺傳結(jié)構(gòu)和系統(tǒng)進(jìn)化分析[10-11]。ddGBS 技術(shù)在有或無參考基因組的情況下，均可通過對(duì)測序結(jié)果進(jìn)行拼接組裝，開發(fā)度SNP 標(biāo)記[12-13]。

多種分子標(biāo)記技術(shù)已應(yīng)用于桉樹群體遺傳和系統(tǒng)進(jìn)化研究，如隨機(jī)擴(kuò)增多態(tài)性（Random amplified polymorphic DNA, RAPD）、簡單重復(fù)序列（Simple sequence repeat, SSR）和擴(kuò)增片段長度多態(tài)性（Amplified fragment length polymorphism,AFLP）標(biāo)記以及多樣性陣列技術(shù)（Diversity array technology, DArT）[14-18]。盡管利用分子標(biāo)記對(duì)桉樹的研究已開展多年，但仍有許多系統(tǒng)發(fā)育問題和分類問題尚待解決。而準(zhǔn)確的系統(tǒng)發(fā)育對(duì)于理解進(jìn)化具有重要的作用，應(yīng)用也非常廣泛。例如將進(jìn)化多樣性納入生物多樣性測量指標(biāo)[19]，為生態(tài)學(xué)研究提供系統(tǒng)發(fā)育觀點(diǎn)[20]，預(yù)測物種對(duì)害蟲和病原體的敏感性[21]等。本研究利用ddGBS 技術(shù)對(duì)14 種不同桉樹（28 份樣品）進(jìn)行建庫測序、開發(fā)高質(zhì)量SNP 標(biāo)記并進(jìn)行基因分型，通過開展系統(tǒng)進(jìn)化樹的構(gòu)建研究，為桉樹的分類、分子鑒定及種質(zhì)資源評(píng)價(jià)提供大量的標(biāo)記資源。

1 材料與方法

1.1 試驗(yàn)材料

本試驗(yàn)用的14 種桉樹來自課題前期收集的材料，分別是：檸檬桉C.citriodora、方格皮桉C.tessellaris、粗皮桉E.pellita、尾葉桉E.urophylla、亮果桉E.nitens、鄧恩桉E.dunnii、藍(lán) 桉E.globulus、赤 桉E.camaldulensis、細(xì) 葉桉E.tereticornis、加拉桉E.marginata、銀頂山岑桉E.sieberi、克蒂斯桉E.curtisii、四齒桉E.tetrodonta和王桉E.regnans。每種桉樹選擇2 棵，均來自不同的家系（詳細(xì)的種源信息及分類學(xué)狀態(tài)列于表1）。物種名稱和分類學(xué)來自Hill 和Johnson 于1995年公布的分類系統(tǒng)[1]。

1.1.1 DNA 樣品制備

采用改良的CTAB 法提取基因組DNA[22]，保存于-80℃冰箱備用，使用瓊脂糖凝膠電泳檢測DNA 質(zhì)量，并利用超微量分光光度計(jì)（Nanodrop 2000, Thermo fisher scientific）測定DNA 濃度和純度（即OD260 與OD280 的比值分布在1.8～2.0之間），再利用Qubit 3.0 熒光定量儀（Thermo Fisher scientific）對(duì)樣品DNA 進(jìn)行定量。

1.1.2 GBS 方案選擇及測序文庫構(gòu)建

建庫主要參照Poland 等[23]方案，選取在基因組上酶切位點(diǎn)分布較多的MspI 酶（C|CGG）和分布較少的PstI 內(nèi)切酶（CTGCA|G）以減少基因組的復(fù)雜性。首先利用Restriction Digest 酶切軟件（https://github.com/JINPENG-WANG/RestrictionDigest）分析巨桉參考基因組序列（https://phytozome.jgi.doe.gov），預(yù)估MspI 和PstI 的 電子酶切效果，推測其在目標(biāo)文庫范圍內(nèi)的酶切片段數(shù)。參照Poland 等[23]的建庫流程（酶切—連接—混樣—PCR 擴(kuò)增—純化），本研究將混樣環(huán)節(jié)改進(jìn)為：等納摩爾質(zhì)量混樣，雙輪磁珠分選350 bp 片段，分選磁珠采用諾唯贊VAHTS DNA Clean Beads。使 用Thermo Fisher Scientific Qubit 3.0 檢測文庫濃度，之后利用PekinElmer Labchip GXII Touch 微流控毛細(xì)管電泳系統(tǒng)對(duì)文庫進(jìn)行質(zhì)檢。質(zhì)檢合格后送至北京貝瑞基因進(jìn)行雙端測序（PE 150），測序平臺(tái)為Illumina HiSeq X Ten。

1.2 數(shù)據(jù)分析

1.2.1 SNP 開發(fā)與注釋

利用Stacks V2.41（http://catchenlab.life.illinois.edu/stacks/）軟件拆分下機(jī)數(shù)據(jù)，按照如下標(biāo)準(zhǔn)去除低質(zhì)量序列：序列中N 的比例≥10%及序列中＞50%的堿基質(zhì)量低于5 的序列。去除低質(zhì)量序列后通過Bowtie2 v2.2.2.9（https://github.com/BenLangmead/bowtie2）軟件與巨桉基因組進(jìn)行比對(duì)，然后利用GATK 3.8（https://software.broadinstitute.org/gatk/）軟件開發(fā)SNP 位點(diǎn)。

表1 14 種桉樹的分類學(xué)狀態(tài)及種源家系信息?Table 1 The information of fourteen Eucalyptus used for analysis

使用Vcftools v4.2（http://vcftools.sourceforge.net/）軟件對(duì)檢測到的SNP 位點(diǎn)數(shù)據(jù)進(jìn)行了過濾，過濾參數(shù)為：最小平均測序深度為4，無缺失數(shù)據(jù)，最小雜合比為0.05，從而獲得高質(zhì)量的SNP。使用SnpEff v4.3（http://snpeff.sourceforge.net/）軟件對(duì)這些SNP 位點(diǎn)進(jìn)行注釋。

1.2.2 系統(tǒng)進(jìn)化樹的構(gòu)建及核苷酸多樣性分析

基于過濾后的高質(zhì)量SNP 位點(diǎn)，利用IQTREE v1.6.11（http://www.iqtree.org/） 軟件采用最大似然法（Maximum Likelihood, ML） 對(duì)28個(gè)樣品進(jìn)行聚類分析，構(gòu)建系統(tǒng)進(jìn)化樹，再利用Figtree v1.4.4（http://tree.bio.ed.ac.uk/software/figtree/）軟件生成可視化樹文件。

利用Vcftools 軟件計(jì)算核苷酸多態(tài)性，以200 kb 為滑動(dòng)窗口，100 kb 為步長，進(jìn)行核苷酸多樣性（兩條序列中不同核苷酸的比例）分析。

2 結(jié)果與分析

2.1 測序質(zhì)量及SNP 分布

基因組DNA 經(jīng)雙酶切建庫測序，去除低質(zhì)量序列后，共獲得30.49 Gb 數(shù)據(jù)，平均每個(gè)樣品的數(shù)據(jù)量為1.09 Gb。如表2所示：測序共產(chǎn)生210 276 036 條序列，樣品平均序列條數(shù)為7 509 858，14 個(gè)樹種產(chǎn)生的測序片段數(shù)存在差異，其中最多的為加拉桉，平均測序片段數(shù)為10 913 288 條，最少為藍(lán)桉，平均產(chǎn)生4 741 137 條。測序Q30 值在97.7%～98.3%范圍內(nèi)，平均為98.0%；GC 含量為43%～47%，平均的GC 含量為44%。樣品間比對(duì)匹配率在46.0%～88.9%之間，平均為77.8%，其中，比對(duì)匹配率最高為尾葉桉（88.9%），最低為方格皮桉（45.0%）。各樹種與巨桉平均對(duì)比匹配率在47%～49%之間的是方格皮桉和檸檬桉；在66%～84.5%范圍的包括加拉桉、銀頂山岑桉、四齒桉、王桉和克蒂斯桉；在86%～89%范圍的是粗皮桉、赤桉、藍(lán)桉、細(xì)葉桉、尾葉桉、鄧恩桉和亮果桉。

利用GATK 軟件共檢測到2 606 828 個(gè)SNP位點(diǎn)，經(jīng)Vcftools 軟件對(duì)SNP 位點(diǎn)數(shù)據(jù)進(jìn)過濾后，保留下42 222 個(gè)高可信度SNP。根據(jù)繪制的SNP 位點(diǎn)在巨桉參考基因組上的分布（圖1）可見，SNP 位點(diǎn)在染色體上分布較為均勻，約99.6%（42 059）的SNP 成功定位在巨桉參考基因組11條染色體上，僅有約0.4%（163 個(gè)）的SNP 無法定位到11 條染色體上。對(duì)高可信度SNP 進(jìn)行功能注釋，統(tǒng)計(jì)結(jié)果（表3）表明：SNP 標(biāo)記約有29.0%位于外顯子區(qū)域，內(nèi)含子區(qū)、上游和下游區(qū)域分別有11.3%、13.2%和12.4%，5′和3′端的UTR 區(qū)含有5.1%的SNP 位點(diǎn)，剪切位點(diǎn)為2.2%，基因間區(qū)僅有1.9%。

表2 桉樹GBS 測序數(shù)據(jù)統(tǒng)計(jì)Table 2 GBS sequencing data statistics in fourteen eucalyptus

2.2 進(jìn)化樹構(gòu)建及多樣性分析

基于42 222 個(gè)高質(zhì)量SNP，統(tǒng)計(jì)各樹種的特有SNP 位點(diǎn)（即僅在該樹種存在的SNP 位點(diǎn)）數(shù)目（表4），分別包括樣品特有SNP 位點(diǎn)數(shù)目和物種特有SNP 位點(diǎn)數(shù)目。14 種桉樹的SNP 數(shù)目存在差異，其中克蒂斯桉（3 348 個(gè)）、方格皮桉（1 509 個(gè)）和王桉（1 610 個(gè)）的特有SNP 較多，檸檬桉（841 個(gè)）、加拉桉（593 個(gè)）次之，粗皮桉（57個(gè)）和細(xì)葉桉（56 個(gè)）最少。14 種桉樹的核苷酸多樣性值主要分布在1.03×10-6～1.07×10-4之間，最大值為1.95×10-4，平均值為2.82×10-5（圖2）。

圖1 SNP 位點(diǎn)在巨桉基因組染色體上的分布Fig.1 Distribution of SNP loci on the chromosome of E.grandis genome

表3 SNP 注釋結(jié)果Table 3 The results of SNP annotation

表4 桉樹樹種的特有SNP 位點(diǎn)統(tǒng)計(jì)?Table 4 The statistics of unique SNP in fourteen eucalypts

圖2 14 種桉樹的核酸多樣性Fig.2 Nucleotide diversity value of 14 species in eucalypts

ML 法構(gòu)建的進(jìn)化樹結(jié)果（圖3）顯示：14 個(gè)樹種可以聚為2 大類，聚類一是傘房屬布萊克亞屬的檸檬桉和方格皮桉，聚類二是其他12 個(gè)樹種。聚類二分為2 個(gè)分支：第一分支包括單蒴蓋亞屬的加拉桉、王桉和銀頂山岑桉，高伯亞屬的克蒂斯桉、紋蒴亞屬的四齒桉。第二分支含有7 個(gè)來自雙蒴蓋亞屬的樹種：粗皮桉、藍(lán)桉、細(xì)葉桉、尾葉桉、赤桉、亮果桉以及鄧恩桉，其中粗皮桉和尾葉桉同屬橫脈組脂桉系，細(xì)葉桉和赤桉隸屬于窿緣組，鄧恩桉、亮果桉和藍(lán)桉同屬藍(lán)桉組多枝桉系。

圖3 14 種桉樹的系統(tǒng)進(jìn)化樹Fig.3 Systematic evolutionary tree of fourteen species in eucalypts

3 結(jié)論與討論

本研究利用ddGBS 技術(shù)對(duì)14 個(gè)桉樹樹種進(jìn)行測序分型，發(fā)掘42 222 個(gè)高質(zhì)量SNP 位點(diǎn)，可廣泛應(yīng)用于桉樹SNP 開發(fā)。并在此基礎(chǔ)上進(jìn)行系統(tǒng)進(jìn)化樹分析，結(jié)果與Hill 和Johnson 的桉屬雙蒴蓋亞屬及傘房屬的分類結(jié)果一致，同時(shí)也和新的分類系統(tǒng)保持一致，這為桉樹系統(tǒng)分類提供新思路，對(duì)于進(jìn)一步開展桉樹關(guān)聯(lián)遺傳學(xué)研究和重要性狀QTL 定位等方面具有重要的應(yīng)用價(jià)值。

GBS 技術(shù)受內(nèi)切酶種類的限制，有研究表明物種酶切效率因選擇的內(nèi)切酶不同而有所差異[24]，良好的酶切組合有利于提高文庫質(zhì)量，增加標(biāo)簽數(shù)量和SNP 位點(diǎn)數(shù)[25]。目前在植物中，ApeKI、PstI 和EcoRI 通常與MspI、MseI 和HpaII 等酶結(jié)合使用，這取決于每個(gè)物種的基因組特異性。前期陳升侃等[26]利用EcoRI 和HindIII 酶切組合開發(fā)尾葉桉×細(xì)葉桉雜種SNP 位點(diǎn)，并應(yīng)用于材性性狀的關(guān)聯(lián)分析；Timothy 等人[27]利用GBS 技術(shù)和桉樹60K 芯片對(duì)藍(lán)桉無性系子代群體進(jìn)行SNP開發(fā)及分型，比較發(fā)現(xiàn)GBS 跟桉樹60K 芯片有13 645 個(gè)不同的SNP 位點(diǎn)；Collins 等[31]采用單酶PstI 從基因組學(xué)的角度證明E.magnificata是一個(gè)獨(dú)特的分類群。本研究采用PstI 和MspI 內(nèi)切酶組合來開發(fā)桉樹高質(zhì)量SNP 標(biāo)記，并應(yīng)用于桉樹系統(tǒng)發(fā)育研究。

形態(tài)學(xué)、植物化學(xué)和基因組學(xué)分析均在桉樹分類研究中發(fā)揮重要作用[28-30]，其中基因組數(shù)據(jù)是闡明物種進(jìn)化和分化過程中的有力工具。二代測序數(shù)據(jù)在很大程度上與形態(tài)學(xué)和植物化學(xué)分析一致，同時(shí)也能進(jìn)一步推斷出不同種群之間的變異[31]。基于ddGBS 技術(shù)的進(jìn)化樹構(gòu)建為桉樹系統(tǒng)發(fā)育研究提供了良好且低成本的技術(shù)支持。如Collins 等[31]整合形態(tài)學(xué)、植物化學(xué)和GBS 的方法發(fā)現(xiàn)并鑒定了一個(gè)桉樹新的稀有種，并且GBS的結(jié)果也與形態(tài)學(xué)和植物化學(xué)的結(jié)果相互印證。前人的研究表明，檸檬桉和方格皮桉同屬于傘房屬的布萊克亞屬[1]，本研究構(gòu)建的系統(tǒng)發(fā)育樹中，傘房亞屬的檸檬桉和方格皮桉聚為一類，這也完全與桉樹的最新版分類系統(tǒng)一致，證實(shí)了本研究的可靠性[32]。

前期Steane 等[16]和Jones 等[17]利用可覆蓋全基因組的DArT 標(biāo)記對(duì)雙蒴蓋亞屬桉樹進(jìn)行了大樣本多類群的系統(tǒng)發(fā)育研究，為進(jìn)化和生態(tài)研究提供了一個(gè)強(qiáng)有力的系統(tǒng)發(fā)育框架。Nicolle 等人[33]也在Jones 和Steane 研究的基礎(chǔ)上加入野外調(diào)查數(shù)據(jù)、植物標(biāo)本及其他系統(tǒng)發(fā)育的研究結(jié)果，對(duì)雙蒴蓋亞屬桉樹分類進(jìn)行了補(bǔ)充和修訂。本研究采用的ddGBS 技術(shù)以經(jīng)濟(jì)簡便的方法開發(fā)SNP 來構(gòu)建系統(tǒng)發(fā)育樹，其中雙蒴蓋亞屬中的樹種聚為平行的三個(gè)分支：橫脈組脂桉系（粗皮桉和尾葉桉）、窿緣組（細(xì)葉桉和赤桉）及藍(lán)桉組多枝桉系（鄧恩桉、亮果桉和藍(lán)桉），結(jié)果與上述的研究基本一致[17,33]，同時(shí)與Hill 等的分類高度一致。其他亞屬的分類結(jié)果，如高伯亞屬的克蒂斯桉、紋蒴亞屬的四齒桉、單蒴蓋亞屬的加拉桉、王桉和銀頂山岑桉，與Hill 等的分類系統(tǒng)及Nicolle 的最新分類系統(tǒng)一致性較高。

14 個(gè)桉樹樹種與巨桉參考基因組的平均比對(duì)匹配率存在差異，其分布規(guī)律與進(jìn)化樹存在相似性，故推測可能是緣于樹種與巨桉親緣關(guān)系遠(yuǎn)近不同有關(guān)，如來自桉屬的粗皮桉、藍(lán)桉、細(xì)葉桉、尾葉桉、赤桉、亮果桉以及鄧恩桉與巨桉同屬雙蒴蓋亞屬，其比對(duì)匹配率較高（86%～89%），而親緣關(guān)系較遠(yuǎn)的來自傘房屬的檸檬桉和方格皮桉，其匹配率則低于50%。

GBS 技術(shù)以其操作過程簡單且使用甲基化敏感酶回避了基因組重復(fù)區(qū)域成為了簡化基因組測序的主流技術(shù)，當(dāng)前也廣泛應(yīng)用于系統(tǒng)發(fā)育、全基因組關(guān)聯(lián)、QTLs 定位、連鎖圖譜構(gòu)建等領(lǐng)域[10-12]。本研究驗(yàn)證了ddGBS 方法構(gòu)建進(jìn)化樹的可靠性，在后續(xù)的試驗(yàn)中，將通過使用更多桉樹基因組特異性的酶切組合，以高效的開發(fā)桉樹SNP 標(biāo)記。Thornhill 等人[34]通過matK 和psbA-trnH 基因，及核內(nèi)轉(zhuǎn)錄間隔區(qū)、外轉(zhuǎn)錄間隔區(qū)的DNA 序列，來評(píng)估桉樹的分類與711 個(gè)物種的系統(tǒng)發(fā)育的匹配度，發(fā)現(xiàn)杯果木屬不能單獨(dú)成為一個(gè)屬，而應(yīng)歸于傘房屬，或者兩個(gè)屬親緣關(guān)系非常近。下一步研究可增加其他更多桉樹樹種的分析，例如杯果木屬，這對(duì)于探索杯果木屬與傘房屬的遺傳進(jìn)化關(guān)系具有重大的意義。