江東，王旭，李仁靜, 趙曉東, 戴祥生，柳正葳

基于GBS技術(shù)開展柚類資源群體遺傳評(píng)價(jià)并發(fā)掘酸含量相關(guān)基因

江東1，王旭1，李仁靜1, 趙曉東2, 戴祥生2，柳正葳3

1西南大學(xué)柑橘研究所，重慶 400712；2井岡山農(nóng)業(yè)科技園管理委員會(huì)，江西吉安 343016；3井岡山大學(xué)，江西吉安 343016

【目的】弄清我國(guó)柚類種質(zhì)資源的起源與演化、群體遺傳結(jié)構(gòu)和多樣性水平，高效利用柚類自然資源群體發(fā)掘與果實(shí)重要品質(zhì)性狀相關(guān)的基因，為柚類種質(zhì)資源群體的遺傳結(jié)構(gòu)研究、起源演化和柚類種質(zhì)創(chuàng)新提供重要的理論及實(shí)踐依據(jù)。【方法】利用國(guó)家果樹種質(zhì)重慶柑橘資源圃中保存的具有廣泛遺傳多樣性和不同地理來源的282份柚類資源作為材料，采用R I限制性內(nèi)切酶消化基因組DNA后構(gòu)建GBS文庫，然后進(jìn)行Illumina HiSeq PE150二代測(cè)序獲得短讀序列，通過BWA軟件將序列映射到柚參考基因組上，利用SAMTOOLS軟件鑒定SNP位點(diǎn)。依據(jù)SNP的基因分型結(jié)果，進(jìn)行主成分分析和群體結(jié)構(gòu)分析,并采用鄰接法構(gòu)建系統(tǒng)演化樹，分析亞群的遺傳多樣性水平。選擇23個(gè)低酸種質(zhì)和32個(gè)高酸種質(zhì)構(gòu)成兩個(gè)群體進(jìn)行Fst、XP-CLR分析，同時(shí)利用282個(gè)柚類種質(zhì)的基因型數(shù)據(jù)與果實(shí)可滴定酸含量的表型數(shù)據(jù)進(jìn)行GWAS全基因組關(guān)聯(lián)分析?！窘Y(jié)果】利用GBS簡(jiǎn)化基因組測(cè)序技術(shù)對(duì)282份柚類種質(zhì)進(jìn)行測(cè)序，獲得201.66 Gb的原始測(cè)序數(shù)據(jù)，每個(gè)樣本的平均測(cè)序數(shù)據(jù)量為0.72 Gb，經(jīng)過測(cè)序深度為5×次要等位基因頻率＞0.01、Miss0.2的篩選條件過濾，最終獲得121 726個(gè)高質(zhì)量的SNP位點(diǎn)。主成分分析、群體結(jié)構(gòu)分析和聚類分析均表明282份柚類資源可被分為6個(gè)亞群，其中柚與‘清見’雜交群體、葡萄柚等橘柚雜種群體可明顯區(qū)別于其他真正柚類種質(zhì)。不同地理來源和特定類型的柚類種質(zhì)在遺傳水平上存在明顯差異，來源于泰國(guó)、越南等東南亞國(guó)家的柚類資源形成了較為獨(dú)特的類群，與國(guó)內(nèi)的沙田柚類、文旦柚類、墊江柚類等群體能夠明顯區(qū)分開。在中國(guó)南方的大部分柚類種質(zhì)中有來自于越南柚的基因滲入，表明越南是柚的原始起源中心。另外，一些來源不清的柚類種質(zhì)也通過GBS技術(shù)得以準(zhǔn)確鑒定。Fst、XP-CLR選擇性清除分析發(fā)現(xiàn)在7號(hào)染色體上的強(qiáng)烈選擇信號(hào)區(qū)域中包含有丙酮酸脫氫酶復(fù)合物二氫硫辛酰賴氨酸殘基乙酰轉(zhuǎn)移酶（PDH-E2）和鋁激活蘋果酸轉(zhuǎn)運(yùn)蛋白（ALMT9），涉及檸檬酸的合成和運(yùn)輸。GWAS全基因組關(guān)聯(lián)分析在2號(hào)染色體上鑒定了一個(gè)與果實(shí)酸度高度關(guān)聯(lián)的區(qū)域?！窘Y(jié)論】GBS技術(shù)為研究柚的系統(tǒng)發(fā)育和演化提供了一種可靠和高效的方法。本研究表明人工雜交育種、長(zhǎng)期人工選擇和馴化、地理隔離是形成不同類型柚類種質(zhì)的驅(qū)動(dòng)力，同時(shí)也是導(dǎo)致柚類種質(zhì)遺傳分化和多樣性增加的重要原因。通過Fst、XP-CLR選擇性清除和全基因組關(guān)聯(lián)分析鑒定了多個(gè)與柚類果實(shí)檸檬酸含量相關(guān)的候選基因，為柚類的遺傳改良和育種提供了重要的基因資源。

GBS；柚系統(tǒng)演化；果實(shí)酸度調(diào)控基因；GWAS

0 引言

【研究意義】柚是柑橘屬的4個(gè)基本種之一[1-2]，其原生起源中心在南亞和東南亞國(guó)家的馬來西亞、印尼、泰國(guó)、越南等地，中國(guó)是世界柚類資源的次生演化中心。柚在傳播引入到我國(guó)后，與我國(guó)原有的寬皮橘、香橙、宜昌橙等資源進(jìn)行天然雜交而逐漸演化形成了甜橙[1-3]、橘柚、酸橙[4-5]、香圓[6]等不同類型的柑橘種質(zhì)，尤其對(duì)寬皮柑橘中雜柑的起源演化產(chǎn)生了重要影響[1,7-8]。由于柚種子為單胚，其實(shí)生后代的變異十分豐富，不同地區(qū)在引種過程中常通過實(shí)生選種獲得新的類型，在長(zhǎng)期的引種栽培馴化過程中，不同地區(qū)的柚類種質(zhì)由于地理隔離而保留了明顯的遺傳印跡，逐漸形成了獨(dú)具特色的地方品種。我國(guó)有豐富多樣的柚類種質(zhì)資源，這為研究柚的起源和演化提供了重要材料；但在以往的引種過程中，常用種子進(jìn)行實(shí)生繁殖，也造成了柚類品種同物異名、同名異物的現(xiàn)象十分突出，對(duì)資源的保存、鑒定帶來一定難度，對(duì)這些種質(zhì)進(jìn)行遺傳多樣性分析，有助于了解柚類種質(zhì)的遺傳背景和差異，從而為柚類資源的高效保存、基因資源挖掘提供科學(xué)的指導(dǎo)。長(zhǎng)期以來，園藝學(xué)上對(duì)柚類品種的認(rèn)識(shí)和劃分多依據(jù)形態(tài)學(xué)特征，由于表型數(shù)據(jù)易受外界環(huán)境因素的影響，不利于柚類種質(zhì)資源的精準(zhǔn)鑒定與評(píng)價(jià)。同時(shí)柚類樹體高大，田間資源保存占用土地面積較多，為加強(qiáng)柚類核心種質(zhì)的構(gòu)建，有必要弄清柚類種質(zhì)群體的遺傳背景和其中的特異資源，這對(duì)提高柚類種質(zhì)的高效保存與利用都具有重要意義。在柚類果實(shí)中，檸檬酸含量是影響果實(shí)品質(zhì)的重要因子，在柚類種質(zhì)資源中存在著低酸和高酸兩類種質(zhì)，這為利用自然資源群體開展基因型和表型的關(guān)聯(lián)分析，進(jìn)而發(fā)掘與果實(shí)酸代謝相關(guān)的基因奠定了基礎(chǔ)。本研究利用GBS簡(jiǎn)化基因組測(cè)序技術(shù)開展柚類種質(zhì)的基因型鑒定和遺傳多樣性分析，將對(duì)柚類種質(zhì)的收集保存、精準(zhǔn)鑒定、起源演化、基因挖掘以及資源的高效利用提供重要信息?！厩叭搜芯窟M(jìn)展】前期國(guó)內(nèi)對(duì)柚類種質(zhì)資源常采用分子標(biāo)記技術(shù)開展資源多樣性的研究，獲得了一些初步結(jié)果。比如Yu等[9]利用31個(gè)SSR分子標(biāo)記對(duì)274份國(guó)內(nèi)柚類種質(zhì)進(jìn)行了遺傳多樣性分析，發(fā)現(xiàn)國(guó)內(nèi)柚類種質(zhì)可分為3個(gè)亞群，這3個(gè)亞群與地理起源分布存在一定相關(guān)性。劉勇等[10]利用SSR分子標(biāo)記對(duì)122份我國(guó)柚類資源及近緣種進(jìn)行了遺傳多樣性分析，發(fā)現(xiàn)柚類資源主要由沙田柚品種群、文旦柚品種群和雜種柚品種群組成。GBS技術(shù)是近幾年來在第二代深度測(cè)序基礎(chǔ)上發(fā)展起來的簡(jiǎn)化基因組測(cè)序技術(shù)，通過DNA酶切后加上代表樣品的分組標(biāo)簽，實(shí)現(xiàn)了多樣本的高通量平行測(cè)序，不僅降低了測(cè)序成本，而且能夠?qū)Υ罅繕颖具M(jìn)行同時(shí)全基因組分型，為開展資源的遺傳研究提供了重要手段和方法[11-12]，利用GBS技術(shù)深入了解資源的遺傳背景和群體遺傳結(jié)構(gòu)、多樣化水平，對(duì)指導(dǎo)資源的收集、保存、利用都具有重要意義[13]。目前利用GBS測(cè)序技術(shù)進(jìn)行SNP發(fā)掘和基因分型已經(jīng)廣泛應(yīng)用在油橄欖[12]、柑橘[13]、棗[14]、油菜[15]、小麥[16-17]、菜豆[18]等作物的遺傳多樣性和群體結(jié)構(gòu)研究中。若聯(lián)合精準(zhǔn)的表型鑒定和全基因組關(guān)聯(lián)分析（GWAS），該技術(shù)還能夠加快對(duì)園藝作物重要性狀基因的挖掘[19-21]。另外，GBS技術(shù)還在全基因組選擇[22]、遺傳圖譜構(gòu)建[23-25]、遺傳圖譜加密[26]、基因組的輔助組裝[27]、品種資源鑒定[28-29]等研究領(lǐng)域發(fā)揮著重要作用。在柑橘類果實(shí)中，前期發(fā)現(xiàn)液泡P-ATP酶citPH1和citPH5對(duì)果實(shí)酸度有重要影響[30]，在柚類中是否還有其他的調(diào)控果實(shí)酸度的基因存在，值得深入研究?！颈狙芯壳腥朦c(diǎn)】作物的遺傳多樣性是品種改良的基礎(chǔ)，我國(guó)有豐富的柚類種質(zhì)資源，對(duì)這些種質(zhì)資源的精準(zhǔn)鑒定和深入評(píng)價(jià)將促進(jìn)資源的深度利用，對(duì)提高柚類品種的育種效率具有重要意義。目前利用GBS技術(shù)開展柚類資源的遺傳多樣性研究還未見報(bào)道。【擬解決的關(guān)鍵問題】利用GBS技術(shù)對(duì)國(guó)家柑橘種質(zhì)資源圃中保存的282份柚類資源進(jìn)行簡(jiǎn)化基因組測(cè)序，獲得不同柚類品種的基因型數(shù)據(jù)，利用SNP基因型數(shù)據(jù)開展柚類群體和不同亞群的多樣性分析，為闡明我國(guó)柚類種質(zhì)的來源、遺傳多樣性水平，進(jìn)而指導(dǎo)柚類核心種質(zhì)的構(gòu)建和柚類品種資源的保護(hù)提供重要的理論參考和依據(jù)。同時(shí)利用基因型和表型關(guān)聯(lián)分析發(fā)掘出與果實(shí)酸含量相關(guān)的基因，為柚類的低酸品種創(chuàng)新提供育種輔助選擇分子標(biāo)記。

1 材料與方法

1.1 試驗(yàn)材料

為使本研究的柚類品種資源具有更廣泛的遺傳多樣性和代表性，采用了282份柚類資源材料進(jìn)行研究，這些種質(zhì)是在不同時(shí)期，從越南、泰國(guó)以及國(guó)內(nèi)不同地區(qū)的資源調(diào)查活動(dòng)中收集，并保存在國(guó)家柑橘種質(zhì)資源圃（重慶）中；另外，材料中還包括了柚的天然雜種和人工雜交材料。選擇的282份柚類種質(zhì)中包括真正的柚類種質(zhì)194份、‘清見’與柚的人工雜交群體材料48份、與柚類親緣關(guān)系較近的40份葡萄柚及類似雜種。另外，還有來源于越南的柚類種質(zhì)14份，來源于泰國(guó)的柚類種質(zhì)10份，研究材料名單見附表1。

1.2 試驗(yàn)方法

1.2.1 文庫構(gòu)建及GBS測(cè)序 田間采集柚類資源無病蟲斑的健康春梢嫩葉，每份材料收集0.3—0.5 g，立刻保存于離心管置于液氮中帶回實(shí)驗(yàn)室，將凍干的葉樣組織用不銹鋼珠研磨成粉末。DNA提取采用Magen Plant DNA Mini Kit試劑盒，提取的DNA樣品送諾禾致源公司進(jìn)行后期的GBS建庫測(cè)序。DNA樣品采用R I酶切后加接頭標(biāo)簽，擴(kuò)增后回收350—400 bp長(zhǎng)度的DNA條帶，這些條帶經(jīng)過純化后上Illumina HiSeq測(cè)序平臺(tái)上進(jìn)行雙末端150 bp測(cè)序。為驗(yàn)證GBS測(cè)序的準(zhǔn)確性和可靠性，在測(cè)序樣本中加入生物學(xué)重復(fù)，同時(shí)一并建庫進(jìn)行分析。

1.2.2 單核苷酸多態(tài)性（SNP）鑒定 依據(jù)建庫樣品和對(duì)應(yīng)條形碼接頭序列聚合相同的樣本，利用fastqc軟件進(jìn)行質(zhì)控，去除低質(zhì)量序列，剩下的序列采用BWA[31]程序?qū)⑵淦ヅ鋵?duì)齊到柚參考基因組上（ftp://ftp.bioinfo.wsu.edu/www.citrusgenomedb.org/Citrus_maxima/），采用參數(shù)為：mem-t4-k32-M。SNP調(diào)用采用Samtools軟件進(jìn)行，最終獲得包含所有SNP變異及位點(diǎn)的VCFs文件。為獲得高質(zhì)量的SNP進(jìn)行后續(xù)分析，采用vcftools軟件設(shè)置篩選條件[32]，以測(cè)序深度大于5、缺失率小于0.2、最小等位基因頻率大于0.01為篩選條件，獲得的SNP位點(diǎn)用于進(jìn)一步的分析。

1.2.3 群體進(jìn)化樹分析 為估測(cè)282份柚類資源的系統(tǒng)演化，根據(jù)SNP基因型采用TreeBest軟件計(jì)算距離矩陣，利用獲得的遺傳矩陣采用鄰接法構(gòu)建系統(tǒng)進(jìn)化樹，設(shè)置自展值初始值為1 000。

1.2.4 主成分分析 通過GCTA軟件計(jì)算特征向量以及特征值，并利用R軟件繪制PCA分布圖。

1.2.5 群體遺傳結(jié)構(gòu)分析 利用PLINK（https://zzz. bwh.harvard.edu/plink/）[33]將vcf文件轉(zhuǎn)換成bed輸入文件，然后利用admixture（http://dalexander.github.io/ admixture/）軟件交叉驗(yàn)證選擇最優(yōu)的K值，設(shè)置K值初始值為3—10，根據(jù)ΔK劃分群體，最終選擇K=3和K=6分別進(jìn)行作圖。

1.2.6 群體遺傳多樣性分析 利用Plink對(duì)選擇的SNP等位基因數(shù)、基因型頻率、基因多樣性（GD）或預(yù)期雜合度（He）、次等位基因頻率（MAF）等進(jìn)行計(jì)算。為確定基因是否存在中性進(jìn)化，采用vcftools軟件計(jì)算6個(gè)亞群的核酸多態(tài)性值和Tajima’s D值，使用“--site-pi”選項(xiàng)估計(jì)每個(gè)SNP位點(diǎn)的核苷酸多樣性π值，取平均值作為群體的π值。設(shè)置窗口長(zhǎng)度為1 000計(jì)算Tajima’s D值，利用Plink計(jì)算群體間的分化系數(shù)Fst。

1.2.7 利用自然選擇和GWAS進(jìn)行果實(shí)酸含量候選基因的定位 為開展與果實(shí)酸含量相關(guān)的基因挖掘，根據(jù)282份柚類種質(zhì)在重慶北碚進(jìn)行多年的果實(shí)酸含量的分析數(shù)據(jù)，選擇低酸（23份，可滴定酸含量＜0.50%）和高酸（32份，可滴定酸含量＞1.00%）種質(zhì)構(gòu)成兩類亞群進(jìn)行Fst、XP-CLR計(jì)算。計(jì)算群體間的Fst值是目前最基礎(chǔ)、也是最常用的檢測(cè)自然選擇的方法，F(xiàn)st可以對(duì)分層群體中不同亞群間基因相關(guān)性進(jìn)行度量，F(xiàn)st的取值范圍為0—1，數(shù)值越大，表明群體間分化程度越大。利用plink軟件逐位點(diǎn)計(jì)算Fst，獲得Fst值后，利用R軟件包qqman完成曼哈頓圖的可視化。利用柑橘基因組數(shù)據(jù)庫（https://www.citrusgenomedb.org/）對(duì)篩選獲得的顯著性較高的SNP位點(diǎn)進(jìn)行基因注釋。利用XP-CLR方法對(duì)Fst鑒定出的7號(hào)染色體進(jìn)行選擇清除分析，利用plink軟件從高酸和低酸兩個(gè)群體中提取基因型和基因位點(diǎn)為輸入文件，選擇掃描窗口為50 kb，步長(zhǎng)為10 kb進(jìn)行XP-CLR分析，計(jì)算出XP-CLR正則化值后利用R軟件包進(jìn)行可視化。對(duì)篩選獲得的XP-CLR最大值對(duì)應(yīng)的位點(diǎn)利用柑橘基因組數(shù)據(jù)庫進(jìn)行注釋。

對(duì)282份柑橘種質(zhì)的可滴定酸表型數(shù)據(jù)與基因型數(shù)據(jù)采用混合線性模型MLM進(jìn)行GWAS全基因組關(guān)聯(lián)分析，MLM公式如下：

y=Xα+Zβ+Wμ+e

通過關(guān)聯(lián)的顯著度篩選出潛在的候選SNP，以-log10(P)＞5為設(shè)定閾值，篩選出候選的SNP位點(diǎn)后，再對(duì)關(guān)聯(lián)位點(diǎn)兩端的20—30 kb連鎖不平衡區(qū)段篩選出候選基因。

1.2.8 芽變材料和同物異名材料的檢測(cè) 利用king軟件（http://people.virginia.edu/～wc9c/KING/Download. htm）檢查不同種質(zhì)直接的親緣關(guān)系，從而推導(dǎo)出芽變材料和同一家族姊妹系材料，進(jìn)而清楚地掌握種質(zhì)材料的遺傳背景和親緣關(guān)系，該結(jié)果除了用于GWAS中親緣關(guān)系的預(yù)測(cè)[34]，還可有效預(yù)測(cè)芽變材料和同物異名的材料。

2 結(jié)果

2.1 測(cè)序質(zhì)量

282份柚類種質(zhì)經(jīng)過DNA質(zhì)量檢測(cè)、酶切建庫和GBS測(cè)序共產(chǎn)生201.66 Gb的測(cè)序數(shù)據(jù)，平均每個(gè)樣本為0.72 Gb，去除低質(zhì)量的測(cè)序數(shù)據(jù)后，保留下的高質(zhì)量序列數(shù)據(jù)量為201.65 Gb，樣本產(chǎn)生的高質(zhì)量測(cè)序數(shù)據(jù)量為317.59—1 237.88 Mb，平均每個(gè)樣本產(chǎn)生的高質(zhì)量測(cè)序數(shù)據(jù)為715.05 Mb。本研究獲得的GBS測(cè)序質(zhì)量較高，Q20≥93.31%，Q30≥83.39%且GC分布正常。將高質(zhì)量的測(cè)序數(shù)據(jù)通過BWA軟件比對(duì)到柚參考基因組，該參考基因組大小為345.78 Mb，群體樣本與柚參考基因組的平均比對(duì)率為97.67%，變異范圍為94.52%—99.03%，說明本研究中的282份柚種質(zhì)的測(cè)序質(zhì)量較高，且不同柚種質(zhì)資源間的遺傳相似度較高。樣本的平均測(cè)序深度為8.31×，序列覆蓋度為17.74%—27.56%，至少有一個(gè)堿基的平均覆蓋度為22.93%，至少有4個(gè)堿基的覆蓋度為8.49%—17.94%，平均覆蓋度為13.44%。其中彭縣柚（131）的比對(duì)率最高，‘清見’雜交柚（282）的比對(duì)率最低。經(jīng)samtools軟件檢測(cè)，在282個(gè)柚類種質(zhì)中獲得了3 007 528個(gè)SNP變異位點(diǎn)，經(jīng)過測(cè)序深度為5、缺失率＜0.2、次要等位基因頻率＞0.01的條件篩選后，最終獲得121 726個(gè)高質(zhì)量的SNP位點(diǎn)。

檢測(cè)發(fā)現(xiàn)SNP在染色體上呈不均勻分布，在著絲粒附件的SNP數(shù)量相對(duì)較少。每條染色體上的SNP平均數(shù)量為12 841，最少為8號(hào)染色體，其上的SNP數(shù)量為8 581；2號(hào)染色體的SNP數(shù)量最多，為21 237；7號(hào)染色體的SNP數(shù)量為9 168。對(duì)獲得的121 726個(gè)SNP位點(diǎn)進(jìn)行注釋，位于基因上游的SNP有7 428個(gè)，基因下游的SNP有6 746個(gè)，在外顯子中導(dǎo)致密碼子提前終止的SNP有197個(gè)，終止子丟失的SNP有45個(gè)，同義突變的SNP有12 046個(gè)，非同義突變的SNP有13 446個(gè)；而存在于內(nèi)含子中的SNP有25 443個(gè)，基因間區(qū)的SNP有48 794個(gè)。單堿基替換產(chǎn)生的轉(zhuǎn)換更多，其中轉(zhuǎn)換的SNP有79 124個(gè)，顛換的SNP有42 602個(gè)，轉(zhuǎn)換/顛換的比值為1.857。觀察到A→G（23 665）和T→C（23 814）的轉(zhuǎn)換較G→A（15 733）和C→T（15 912）的轉(zhuǎn)換多。

2.2 主成分分析

利用PCA分析法對(duì)282份柚類種質(zhì)進(jìn)行主成分分析，其中前3個(gè)主成分能解釋大部分變異。根據(jù)第一主成分和第二主成分可將282份柚類資源分為3個(gè)亞群，分別為真柚群體、葡萄柚群體、‘清見’與柚雜交群體。而根據(jù)第一主成分和第三主成分，又可以將真柚亞群細(xì)分為4個(gè)亞群，分別是沙田柚亞群、越南柚和泰國(guó)柚亞群、文旦柚亞群、墊江柚亞群。綜合前3個(gè)主成分，可將282份資源分為6個(gè)亞群，分別是葡萄柚亞群、‘清見’與柚雜交群體、沙田柚亞群、文旦柚亞群、越南泰國(guó)柚亞群、墊江柚亞群（圖1）。

左圖依據(jù)第1和第2主成分構(gòu)圖，右圖依據(jù)第1和第3主成分構(gòu)圖。紅色圓圈代表‘清見’與柚雜交樣本，黑色圓圈代表真柚類樣本，藍(lán)色和其他顏色圓圈代表葡萄柚樣本

2.3 群體遺傳結(jié)構(gòu)分析

為評(píng)價(jià)282份柚類的群體結(jié)構(gòu)，使用ΔK來劃分亞群，當(dāng)K=3時(shí)觀察到最大的ΔK，表明282份柚類種質(zhì)可分為3組，分別對(duì)應(yīng)是葡萄柚亞群、‘清見’與柚雜種亞群和真柚亞群（圖2，K=3），結(jié)果與PCA的結(jié)果吻合。若要細(xì)分真柚亞群，可設(shè)置K=6時(shí)，將282份柚類種質(zhì)分為6個(gè)群體（圖2，K=6）。沙田柚類群（紅色）為單獨(dú)一個(gè)大類，其中包括真龍柚、合江柚、嶺南沙田柚、長(zhǎng)壽沙田柚等類型以及沙田柚的雜種如脆柚、沙暹柚、沙田柚×甜橙、正形沙田柚、杭晚蜜柚實(shí)生2號(hào)、衛(wèi)寺蜜柚、梁沙柚、舒化柚、早熟柚、金沙柚等，該類群果實(shí)形態(tài)多為梨形，具有明顯的沙田柚的果形特征；第二類為墊江柚類群，包括墊江白柚、墊江紅心柚、墊江周家白心柚、墊江曾家白柚、墊江黃沙白心柚、龍安柚以及灌香柚、北碚柚、蓬溪柚、貢水紅柚、江津紅心柚、金堂綠柚等，果形多為長(zhǎng)橢圓形或卵圓形；第三類為以琯溪蜜柚為代表的文旦柚類型，包括玉環(huán)文旦、三紅蜜柚、福建文旦、紅肉琯溪、紅綿蜜柚、黃肉蜜柚、通賢柚以及文旦柚雜種，包括左氏柚、尖頂梁平柚2號(hào)、四季拋、坪山柚等，果形多為錐狀卵圓形；第四類是‘清見’與柚的雜種群體，其中包括雞尾葡萄柚等材料，這類材料中除了有柚類的基因外，還含有其他柑橘種類的基因滲入，比如橙類的基因，果形比一般真正柚類果實(shí)偏小，果肉顏色多為橙色，果肉細(xì)軟多汁；第五類包括來源于東南亞國(guó)家越南和泰國(guó)的柚類品種，越南柚類包括光皮柚、裴紅柚、平陽柚、囊內(nèi)柚、越南小柚、越南小甜柚、綠柚；泰國(guó)柚包括泰國(guó)柚、暹羅低酸柚、高浦柚、永安柚、泰國(guó)柚kitik、強(qiáng)德勒柚、彭縣暹羅柚，以及與東南亞國(guó)家鄰近地域的品種，如勐侖早柚、晚白柚等，其中國(guó)內(nèi)的一些柚類品種也包括其中，表明這些品種中含有越南、泰國(guó)柚等外源基因的滲入，可能最早是通過品種引進(jìn)再與國(guó)內(nèi)品種雜交后馴化而成。這些品種如川渝地區(qū)的梁平柚、江北無核柚、納溪櫻桃柚、梅塆柚、蒲蓮柚、金堂綠柚、‘清見’ב強(qiáng)德勒’，湖南安江地區(qū)的安江香柚、安江石榴柚、安農(nóng)1號(hào)香柚、大庸菊花心、安農(nóng)無核蜜柚、白玉霜、鍋魁柚，浙江、江西等地的永嘉早香柚、金蘭柚、水晶文旦、盤谷文旦以及一些雜柚品種如8088柚、奇龍嘉、臍柚等。這個(gè)類群中的基因來源復(fù)雜，但主要以越南柚中的基因?yàn)橹饕z傳構(gòu)成。這個(gè)類群多數(shù)柚類品種果形為高扁圓形，成熟期普遍偏晚，表明起源于南亞熱帶地區(qū)的柚類品種在中亞熱帶氣候條件下，成熟期有推遲趨勢(shì)；第六類包括葡萄柚類型品種，主要有葡萄柚、胡柚、溫嶺高橙、建陽橘柚、‘清見’與柚的雜種、‘日本4’×夏橙、奧羅勃朗柯、世界蜜柚等，在這個(gè)類群中，比如像建陽橘柚具有橘類的基因滲入，‘日本4’×夏橙具有更多橙類基因的滲入。

通過對(duì)群體結(jié)構(gòu)的進(jìn)一步分析，發(fā)現(xiàn)相同地理來源的柚類品種通常聚在一起，比如來自墊江及其周邊的柚類品種多聚在一起；來自越南和泰國(guó)的柚類品種通常聚集在越南和泰國(guó)柚類群下，不同地理來源的品種存在較明顯的遺傳差異。

2.4 群體進(jìn)化樹分析

采用鄰接法構(gòu)建了群體進(jìn)化樹，經(jīng)過自展值達(dá)1 000次的計(jì)算獲得。群體進(jìn)化樹結(jié)果表明，282份柚類資源種質(zhì)聚為6個(gè)大類。其中第一個(gè)類群為川渝地區(qū)的地方種質(zhì)，包括墊江柚、墊江柚雜種、北碚柚、蓮蒲柚、梅灣柚等；第二類為沙田柚及其雜種類型，另外與沙田柚地理分布位置較近的越南柚、勐侖早柚、綠肉柚、暹羅蜜柚也與其聚在一起。第三類由琯溪蜜柚等文旦柚類型組成，另外還包括來自于湖南的部分柚類品種，如大庸菊花芯柚等。第四類由泰國(guó)柚類，如暹羅低酸柚、強(qiáng)德勒柚、泰國(guó)柚及一些與真正柚類遺傳距離較遠(yuǎn)的種質(zhì)構(gòu)成，包括水晶文旦柚、枳雀等柚類雜種?？梢娫谡嬲蔫诸惙N質(zhì)資源中，大致可以劃分為泰國(guó)柚類群、文旦柚類群、越南柚與沙田柚類群、墊江柚類群。第五類由‘清見’與柚的雜交群體構(gòu)成。第六類大部分由葡萄柚及其類似雜種組成。另外，從群體的聚類圖中看，水晶文旦、枳雀、285號(hào)柚品種在遺傳水平上存在明顯的異質(zhì)性，這些種質(zhì)應(yīng)該含有更多外源基因的滲入，可以作為核心種質(zhì)加以重點(diǎn)保存。從群體分化的結(jié)果來看，地理隔離、種間雜交、引種馴化是引起柚類種質(zhì)遺傳分化的重要驅(qū)動(dòng)力。

2.5 群體的遺傳多樣性

對(duì)282份柚類種質(zhì)進(jìn)行群體的遺傳多樣性評(píng)估，基因多樣性GD值（即期望雜合度He）在282份柚類種質(zhì)中變化范圍從0.02—0.5，平均值為0.2103；PIC值變化范圍從0.02—0.38，平均值為0.178；MAF值的變化范圍在0.01—0.50，平均值為0.139。利用Vcftools對(duì)6個(gè)不同亞群的p值和Tajima’s D值進(jìn)行計(jì)算（表2）。

從表2可見，6個(gè)亞群的p值變化范圍處于0.01—0.51，葡萄柚亞群的平均核酸多態(tài)性p?值最大，其次是‘清見’雜交柚類群，這主要?dú)w因于兩個(gè)群體均來源于柚與其他柑橘的種間雜交，造成基因組遺傳分化和多樣性水平的增加，因而后代個(gè)體在遺傳組成上與真正柚類種質(zhì)的差異較大。沙田柚亞群的平均p值最低，其次是文旦柚亞群，表明在沙田柚群體和文旦柚群體中可能經(jīng)歷了較大的人工選擇壓力或者人為繁殖導(dǎo)致基因組核酸多態(tài)性水平的降低?；蚪M的Tajiama’s D值往往與外源基因的滲入呈正相關(guān)，本研究發(fā)現(xiàn)‘清見’雜交柚群體的Tajima’s D值最大，平均值為1.12，表明雜交群體中個(gè)體間的異質(zhì)性較大，群體的多樣性水平較高。葡萄柚群體Tajima’s D值為0.92，文旦柚群體的Tajima’s D平均值為0.41，沙田柚群體的Tajima’s D平均值為-0.15，表明沙田柚群體由于長(zhǎng)期的大量引種、栽培面積的擴(kuò)大和人工選擇壓力的作用，亞群內(nèi)個(gè)體間的異質(zhì)性較小，基因組上存在選擇性清除效應(yīng)。利用Plink計(jì)算群體間的分化系數(shù)Fst，越南、泰國(guó)柚群體與葡萄柚群體的Fst值最大為0.24，其次是墊江柚群體和‘清見’雜交柚群體的Fst值為0.22，沙田柚群體和文旦柚群體、墊江柚群體的Fst值接近，分別為0.18和0.16。墊江柚群體與越南、泰國(guó)柚群體的Fst值最小，為0.06，其次是文旦柚群體與越南、泰國(guó)柚群體的Fst值為0.08，‘清見’雜交柚群體與葡萄柚群體的Fst值為0.09。一般認(rèn)為Fst大于0.2，群體間即存在明顯的遺傳分化，葡萄柚群體和‘清見’雜交柚群體與真正柚類亞群的Fst值均大于0.2。Fst值大小可反映親緣關(guān)系的遠(yuǎn)近，可見國(guó)內(nèi)真正柚類品種群與越南、泰國(guó)柚類群的分化遠(yuǎn)小于柚類與葡萄柚、‘清見’雜交柚的遺傳分化。

圖3 利用282份柚類種質(zhì)的121726 SNP采用鄰接法構(gòu)建群體進(jìn)化樹

表2 6個(gè)柚類亞群的遺傳多樣性

2.6 利用自然選擇法對(duì)可滴定酸候選基因的初步定位

果實(shí)中有機(jī)酸的含量是決定果實(shí)品質(zhì)和風(fēng)味的一個(gè)重要性狀，柚類果實(shí)中的可滴定酸絕大部分為檸檬酸[35]。通過對(duì)282份柚類種質(zhì)在重慶北碚地區(qū)多年的果實(shí)可滴定酸含量分析，利用Plink計(jì)算低酸（23）和高酸（32）兩個(gè)群體的Fst值，根據(jù)最大的Fst值，鑒定了4個(gè)可滴定酸的候選基因，均位于7號(hào)染色體上。在7號(hào)染色體上與可滴定酸含量顯著關(guān)聯(lián)的SNP位點(diǎn)，位于Cg7g014530.1和Cg7g014550.1的基因間區(qū)，其靠5′端的鄰近基因?yàn)镃g7g014530.1，該基因注釋為丙酮酸脫氫酶復(fù)合物（PDC）的一個(gè)亞基—二氫硫辛酰乙酰轉(zhuǎn)移酶（PDH-E2），該酶是檸檬酸合成中的一個(gè)重要酶，PDH-E2通過丙酮酸生成乙酰輔酶A，而乙酰輔酶A、草酰乙酸和H2O在ATP-檸檬酸合酶的催化下，在線粒體中以雙向反應(yīng)形成CoA和檸檬酸。由于該酶是參與乙酰輔酶A生成的重要蛋白，因此對(duì)檸檬酸的合成十分重要。對(duì)該SNP位點(diǎn)的基因型分析表明，基因型為CC、TC的為高酸品種，而基因型為TT的為低酸品種。Fst分化系數(shù)最大的SNP位于7號(hào)染色體16 299 216處的Cg7g014580，該基因注釋為絲氨酸/蘇氨酸蛋白激酶PRP4同源物。在該基因附近有Cg7g014630.1，注釋為鋁激活蘋果酸轉(zhuǎn)運(yùn)蛋白4，在這個(gè)基因的下游是Cg7g014640.1，基因注釋為絲氨酸/蘇氨酸蛋白磷酸酶2A 65 kDa調(diào)節(jié)亞基PP2A。利用XP-CLR方法在低酸（23）和高酸（32）兩個(gè)類群中進(jìn)行選擇性清除分析，同樣鑒定出PDH-E2所在基因位置有最強(qiáng)的選擇信號(hào)。

GWAS鑒定了與果實(shí)可滴定酸含量顯著關(guān)聯(lián)的一系列關(guān)聯(lián)位點(diǎn)。其中2號(hào)染色體上1 168 947—1 170 318有最大的關(guān)聯(lián)信號(hào)。對(duì)前后20 Kb區(qū)段進(jìn)行基因注釋，其中最強(qiáng)關(guān)聯(lián)信號(hào)的基因注釋為線粒體DGP2_ARATH DAR GTPase 2，該區(qū)段內(nèi)還包括編碼絲氨酸/蘇氨酸蛋白磷酸酶PP2A的同源基因Cg2g000890以及F-box/kelch-重復(fù)蛋白的基因。而絲氨酸/蘇氨酸蛋白磷酸酶PP2A在Fst和XP-CLR選擇清除的顯著區(qū)段也被發(fā)現(xiàn)與可滴定酸的含量顯著相關(guān)。Fst、XP-CLR和GWAS結(jié)果均表明PP2A與果實(shí)中可滴定酸含量有顯著的相關(guān)性。

Fst最大值對(duì)應(yīng)的SNP位于7號(hào)染色體上The corresponding maximum position of Fst and XP-CLR is on thesame location of chromosome 7

圖5 利用282份柚種質(zhì)采用GWAS分析鑒定柚果實(shí)可滴定酸的曼哈頓圖（A）和QQ圖（B）

2.6 未知親本材料的鑒定

本研究采用King軟件對(duì)柚種質(zhì)個(gè)體間親緣關(guān)系進(jìn)行研究，不僅是利用自然群體進(jìn)行GWAS定位的前提，也為揭示材料間的親緣關(guān)系提供了重要線索。研究發(fā)現(xiàn)King能夠準(zhǔn)確對(duì)生物學(xué)重復(fù)、芽變材料、全同胞家系個(gè)體進(jìn)行準(zhǔn)確鑒定。本研究中，加入了一些生物學(xué)重復(fù)材料，比如2份晚白柚、2份胡柚、2份東風(fēng)早柚、3份斐紅柚都被準(zhǔn)確鑒定出來，表明GBS的基因分型結(jié)果準(zhǔn)確可靠。同時(shí)King能夠?qū)σ恍┐嬉善贩N進(jìn)行區(qū)分，比如64、150、208三個(gè)品種被鑒定為生物學(xué)重復(fù)，其中64、208是龍安柚，而150是緬甸柚的概率應(yīng)該很小，更大可能是龍安柚。另外，泰國(guó)柚3-Kitik與永安柚、潼南柚與琯溪蜜柚、安柚與越南小甜柚、貢水紅柚和江津紅心柚的相似度都較高。墊江白柚和墊江紅心柚的遺傳相似度也很高，其實(shí)生后代也多數(shù)聚集在一起。梁平柚與鳳凰柚、虎蜜柚的親緣系數(shù)很高，表明這些品種間有較近的親緣關(guān)系。一些同名的品種如段氏柚，被劃分到3個(gè)亞群中，明顯屬于同名異物，需要對(duì)這些材料進(jìn)一步核實(shí)。福建文旦柚雖然被劃分到1個(gè)亞類下，但兩份福建文旦柚仍不完全相同，表明這兩份材料盡管親緣關(guān)系較近，可能仍然具有一定遺傳差異。本研究表明利用GBS技術(shù)結(jié)合King軟件來鑒定柚類品種的親緣關(guān)系是可行的。從系統(tǒng)發(fā)育樹和King分析的結(jié)果來看，即使柚類品種是從種子實(shí)生繁殖得來，其全同胞家系個(gè)體間的遺傳組成并沒有太大的變異，實(shí)生后代常與親本材料聚在一起，暗示了柚類種質(zhì)的基因純合水平可能較高，而雜合水平可能較低。

3 討論

3.1 不同地理區(qū)域柚類種質(zhì)資源存在顯著遺傳差異和形態(tài)特征

不同地理區(qū)域由于氣候生態(tài)環(huán)境的不同，對(duì)群體的遺傳特異性塑造產(chǎn)生了重要影響，地理的隔離會(huì)造就獨(dú)特的種群，這在許多植物和動(dòng)物中都已經(jīng)得到證實(shí)[36-37]。本研究表明不同地域的柚類種質(zhì)具有明顯不同的遺傳特征，來源于泰國(guó)、越南和國(guó)內(nèi)的柚類資源存在明顯的遺傳差異，而相同或鄰近區(qū)域內(nèi)的柚類種質(zhì)在遺傳組成上的差異往往較小，比如泰國(guó)的暹羅低酸柚、強(qiáng)德勒柚、泰國(guó)柚等幾個(gè)品種聚在一起，而越南的光皮柚、囊內(nèi)柚、綠柚等明顯聚在一起。墊江白柚類群起源于重慶，由于群體后代較多，也形成了具有特色的地方品種群。沙田柚是華南地區(qū)的重要柚類品種，受人為選擇壓力較大或者人為的廣泛繁殖，因此，該品種類群個(gè)體中保持了明顯的沙田柚特征。綜合來看，柚類品種存在明顯的地理分布特征，柚起源于東南亞地區(qū)，在引入到國(guó)內(nèi)后，為適應(yīng)當(dāng)?shù)氐臍夂蚝腿藶檫x擇壓力，演化出了新的類型，如沙田柚、文旦柚和墊江白柚等類型。

不同類型的柚類亞群具有獨(dú)特的形態(tài)特征和地理起源，如沙田柚亞群起源于華南地區(qū)，樹勢(shì)較強(qiáng)，樹形較直立，果形為梨形，果頸明顯，果肉白色，果實(shí)酸度較低，成熟往往較晚，被引入到不同地區(qū)，形成不同的品系，如長(zhǎng)壽沙田柚、墊江沙田柚、合江柚、真龍柚等。而文旦柚類的典型代表是琯溪蜜柚，起源于東南沿海的福建等地，樹勢(shì)偏弱，樹形扁圓形，果皮光滑，果形尖卵圓形或扁圓形，中心柱多空虛，果實(shí)酸度較高，成熟較早。而墊江柚亞群的種質(zhì)，樹勢(shì)強(qiáng)，樹形較圓頭或開張形，果實(shí)橢圓形，中心柱多空虛，成熟較早。這3個(gè)類群具有明顯的地域特征和遺傳構(gòu)成。這與國(guó)內(nèi)Yu等[9]、劉勇等[10]利用SSR分子標(biāo)記和何天富[38]利用形態(tài)學(xué)特征對(duì)國(guó)內(nèi)柚類種質(zhì)進(jìn)行多樣性分析的結(jié)果相吻合。利用GBS的群體遺傳分析結(jié)果，提出了不同地理區(qū)域的代表性核心種質(zhì)，其中華南地區(qū)的代表柚類種質(zhì)是沙田柚，東南沿海的代表柚類種質(zhì)是琯溪蜜柚，而西南地區(qū)代表柚類種質(zhì)是墊江柚。

3.2 越南柚在我國(guó)柚類起源演化中產(chǎn)生重要影響

地理隔離雖然塑造了品種亞群，但也存在亞群間的基因交流。本研究發(fā)現(xiàn)沙田柚、文旦柚、墊江白柚亞群間通過基因漸滲形成新的種質(zhì)。從柚類起源演化的角度來看，來源于起源中心的柚類品種，更廣泛地將其基因滲入到周邊的柚類品種中。本研究中納入了起源中心的越南和泰國(guó)的柚類資源，從群體結(jié)構(gòu)分析來看，越南柚所提供的基因滲入在真正柚類群體中的作用非常明顯，這可以從墊江柚群體與越南柚群體的Fst值最小為0.06，文旦柚群體與越南柚群體的Fst值為0.08得以證實(shí)，因此，越南的柚類資源對(duì)我國(guó)柚類的起源演化具有重要影響。由于越南柚的基因更多滲入我國(guó)的柚類種質(zhì)中，并與沙田柚、文旦柚和墊江柚存在基因交流，由此推測(cè)，越南是柚的原初起源中心，而次生中心所演化形成的獨(dú)特柚類種質(zhì)，如沙田柚、文旦柚、墊江柚的基因也滲入到國(guó)內(nèi)多數(shù)柚品種中，但不如越南柚的貢獻(xiàn)更大。在沙田柚系中沙暹柚、正形沙田柚、杭晚蜜柚實(shí)生、梁砂柚等資源中都有越南柚基因的滲入，貢水紅柚、江津紅心柚、金堂綠柚等資源中也具有較大比例的越南柚基因的貢獻(xiàn)；而文旦柚系中，除了坪山柚、四季拋、左氏柚、楚門文旦等幾份柚外，其余品種資源中越南柚的基因貢獻(xiàn)很少，表明地理位置的鄰近對(duì)柚類種質(zhì)的形成具有重要的影響和促進(jìn)作用。

3.3 柚類果實(shí)中決定檸檬酸含量基因的初步定位

對(duì)柚類果實(shí)中檸檬酸含量的基因定位，以往采用混池分離法或構(gòu)建雜交群體進(jìn)行QTL定位方式來發(fā)掘。由于GBS能夠獲得高密度的SNP位點(diǎn)，利用性狀迥異的兩類自然群體通過計(jì)算Fst和XP-CLR來發(fā)掘與性狀相關(guān)點(diǎn)的基因在動(dòng)、植物中都已經(jīng)證明是高效可靠的[39-40]。本研究通過對(duì)低酸和高酸柚群體的Fst和XP-CLR選擇性清除分析，定位到7號(hào)染色體上與可滴定酸含量顯著關(guān)聯(lián)的SNP位點(diǎn)，其鄰近的基因編碼丙酮酸脫氫酶復(fù)合物（DHC）的一個(gè)亞基——二氫硫辛酰賴氨酸殘基乙酰轉(zhuǎn)移酶（PDH-E2），該酶主要負(fù)責(zé)乙酰-CoA的生成。由于GBS是通過對(duì)酶切位點(diǎn)附近的序列進(jìn)行測(cè)序分型，SNP的密度相對(duì)重測(cè)序來說要低很多，在基因上有可能缺乏酶切位點(diǎn)而未能鑒定到SNP差異。該研究發(fā)現(xiàn)的SNP雖未直接定位于基因上，但這種結(jié)構(gòu)變異所帶來的“搭車效應(yīng)”也反映出高酸品種和低酸品種群體間，在該SNP位點(diǎn)附近的基因有可能存在其他結(jié)構(gòu)上的變異，進(jìn)而對(duì)基因的功能產(chǎn)生影響。而該SNP位點(diǎn)的基因分型結(jié)果表明，基因型為CC、TC的為高酸品種，而基因型為TT的為低酸品種，證實(shí)該SNP位點(diǎn)可作為低酸柚品種選育的重要分子標(biāo)記加以利用。檸檬酸在細(xì)胞內(nèi)一般通過三羧酸循環(huán)（TCA）途徑不斷產(chǎn)生和循環(huán)。在線粒體內(nèi)，丙酮酸和硫辛酰胺-E首先形成S-乙酰基二氫硫辛酰胺-E，其在PDH-E2的作用下，與(R)-硫辛酸、CoA催化形成二氫硫辛酰胺-E和乙酰-CoA，草酰乙酸、H2O和乙酰-CoA進(jìn)一步在ATP-檸檬酸合酶的催化下，在線粒體中以雙向反應(yīng)形成CoA和檸檬酸。乙酰-CoA是檸檬酸循環(huán)的最重要的起始底物，因此，PDH-E2編碼基因結(jié)構(gòu)或者活性的改變，會(huì)影響其蛋白酶活性或者含量，通過調(diào)節(jié)乙酰-CoA的生成，進(jìn)而影響檸檬酸的合成。另外，DHC的E1亞基還受磷酸化調(diào)控，具有無活性（磷酸化）和有活性（去磷酸化）兩種形式，丙酮酸脫氫酶激酶（PDK）通過磷酸化PDH可抑制其酶活，使乙酰-COA的含量減少，進(jìn)而影響檸檬酸的合成[41]。在其他果樹中，如桃的低酸品種‘歐亨利’，PDK的表達(dá)量增高，減少了乙酰-COA的生產(chǎn)，進(jìn)而影響了果實(shí)中檸檬酸的含量[42]。本研究發(fā)現(xiàn)7號(hào)染色體16 299 216處的絲氨酸/蘇氨酸蛋白激酶PRP4，和2號(hào)染色體上的絲氨酸/蘇氨酸蛋白磷酸酶PP2A都與果實(shí)中檸檬酸的高低存在一定關(guān)聯(lián)，這些激酶和磷酸酶是否能對(duì)DHC的E1亞基進(jìn)行磷酸化或去磷酸化，進(jìn)而調(diào)控果實(shí)中酸的含量還需在后續(xù)試驗(yàn)中進(jìn)一步研究。另外，鋁激活蘋果酸轉(zhuǎn)運(yùn)蛋白4，也可能是影響檸檬酸含量的重要候選基因。在番茄中調(diào)控蘋果酸積累的主效基因編碼一個(gè)鋁激活蘋果酸轉(zhuǎn)運(yùn)蛋白（）[43]，是導(dǎo)致番茄果實(shí)蘋果酸含量變異的重要基因，其啟動(dòng)子區(qū)的GTC序列插入/缺失引起基因表達(dá)差異是影響蘋果酸積累的主要因素。此外，Li等[44]還報(bào)道調(diào)控蘋果酸積累的與具有高同源性，葡萄中不僅具備蘋果酸轉(zhuǎn)運(yùn)能力，還可介導(dǎo)酒石酸轉(zhuǎn)運(yùn)[45]；而在桃中，的下調(diào)，會(huì)引起蘋果酸轉(zhuǎn)運(yùn)到液泡中的含量減少從而導(dǎo)致果實(shí)酸度降低[42]。本研究中，該基因在柚類中是否通過介導(dǎo)檸檬酸的轉(zhuǎn)運(yùn)，進(jìn)而影響果實(shí)的酸度，也還需要繼續(xù)深入研究。

4 結(jié)論

本研究利用GBS技術(shù)對(duì)282份柚類資源進(jìn)行了簡(jiǎn)化基因組測(cè)序，根據(jù)系統(tǒng)演化樹、群體遺傳結(jié)構(gòu)和遺傳多樣性分析，結(jié)合King軟件能夠準(zhǔn)確地對(duì)282份柚類資源進(jìn)行準(zhǔn)確劃分，證明GBS技術(shù)可高效用于資源材料的精準(zhǔn)鑒定和親緣關(guān)系的研究，對(duì)提高柚類資源的保存、育種效率和科學(xué)管理都具有重要價(jià)值。同時(shí)，研究證明了地理隔離、雜交育種、人工選擇是柚種質(zhì)資源遺傳分化和多樣性增加的內(nèi)在驅(qū)動(dòng)力，越南是柚類起源中心，在我國(guó)柚類的起源演化中扮演重要角色。中國(guó)作為柚類資源的次生演化中心，也具有獨(dú)特的核心骨干種質(zhì)，可根據(jù)中國(guó)柚類種質(zhì)的品種特點(diǎn)和地域分布劃分為沙田柚系、文旦柚系、墊江柚系。本研究利用Fst、XP-CLR和GWAS技術(shù)鑒定了柚類果實(shí)中與可滴定酸性狀相關(guān)的重要候選PDH-E2基因和鋁激活蘋果酸轉(zhuǎn)運(yùn)蛋白4，為今后柚類品質(zhì)遺傳改良提供了重要的可候選基因。

[1] WU G A, TEROL J, IBANEZ V, LóPEZ-GARCíA A, PéREZ- ROMáN E, BORREDá C, DOMINGO C, TADEO F R, CARBONELL- CABALLERO J, ALONSO R, CURK F, DU D L, OLLITRAULT P, ROOSE M L, DOPAZO J, GMITTER F G, ROKHSAR D S, TALON M. Genomics of the origin and evolution of. Nature, 2018, 554(7692): 311-316.

[2] WANG X, XU Y T, ZHANG S Q, CAO L, HUANG Y, CHENG J F, WU G Z, TIAN S L, CHEN C L, LIU Y, YU H W, YANG X M, LAN H, WANG N, WANG L, XU J D, JIANG X L, XIE Z Z, TAN M L, LARKIN R M, CHEN L L, MA B G, RUAN Y J, DENG X X, XU Q. Genomic analyses of primitive, wild and cultivated citrus provide insights into asexual reproduction. Nature Genetics, 2017, 49(5): 765-772.

[3] XU Q, CHEN L L, RUAN X A, CHEN D J, ZHU A D, CHEN C L, BERTRAND D, JIAO W B, HAO B H, LYON M P, CHEN J J, GAO S, XING F, LAN H, CHANG J W, GE X H, LEI Y, HU Q, MIAO Y, WANG L, XIAO S X, BISWAS M K, ZENG W F, GUO F, CAO H B, YANG X M, XU X W, CHENG Y J, XU J, LIU J H, LUO O J, TANG Z H, GUO W W, KUANG H H, ZHANG H Y, ROOSE M L, NAGARAJAN N, DENG X X, RUAN Y J. The draft genome of sweet orange (). Nature Genetics, 2013, 45(1): 59-66.

[4] WU G A, SUGIMOTO C, KINJO H, AZAMA C, MITSUBE F, TALON M, GMITTER F G, ROKHSAR D S. Diversification of mandarin citrus by hybrid speciation and apomixis. Nature Communications, 2021, 12(1): 4377.

[5] BARKLEY N A, ROOSE M L, KRUEGER R R, FEDERICI C T. Assessing genetic diversity and population structure in a citrus germplasm collection utilizing simple sequence repeat markers (SSRs). Theoretical and Applied Genetics, 2006, 112(8): 1519-1531.

[6] DEMARCQ B, CAVAILLES M, LAMBERT L, SCHIPPA C, OLLITRAULT P, LURO F. Characterization of odor-active compounds of ichang lemon (Tan.) and identification of its genetic interspecific origin by DNA genotyping. Journal of Agricultural and Food Chemistry, 2021, 69(10): 3175-3188.

[7] GONZALEZ-IBEAS D, IBANEZ V, PEREZ-ROMAN E, BORREDá C, TEROL J, TALON M. Shaping the biology of citrus: II. Genomic determinants of domestication. The Plant Genome, 2021, 14(3): e20133.

[8] WU G A, PROCHNIK S, JENKINS J, SALSE J, HELLSTEN U, MURAT F, PERRIER X, RUIZ M, SCALABRIN S, TEROL J,. Sequencing of diverse mandarin, pummelo and orange genomes reveals complex history of admixture during citrus domestication. Nature Biotechnology, 2014, 32(7): 656-662.

[9] YU H W, YANG X M, GUO F, JIANG X L, DENG X X, XU Q. Genetic diversity and population structure of pummelo () germplasm in China. Tree Genetics & Genomes, 2017, 13(3): 58.

[10] 劉勇, 劉德春, 吳波, 孫中海. 利用SSR標(biāo)記對(duì)中國(guó)柚類資源及近緣種遺傳多樣性研究. 農(nóng)業(yè)生物技術(shù)學(xué)報(bào), 2006, 14(1): 90-95.

LIU Y, LIU D C, WU B, SUN Z H. Genetic diversity of pummelo and their relatives based on SSR markers. Journal of Agricultural Biotechnology, 2006, 14(1): 90-95. (in Chinese)

[11] ELSHIRE R J, GLAUBITZ J C, SUN Q, POLAND J A, KAWAMOTO K, BUCKLER E S, MITCHELL S E. A robust, simple genotyping-by-sequencing (GBS) approach for high diversity species. PLoS ONE, 2011, 6(5): e19379.

[12] Islam A S M F, SANDERS D, MISHRA A K, JOSHI V. Genetic diversity and population structure analysis of the USDA olive germplasm using genotyping-by-sequencing (GBS). Genes, 2021, 12(12): 2007.

[13] 王小柯, 江東, 孫珍珠. 利用GBS技術(shù)研究240份寬皮柑橘的系統(tǒng)演化. 中國(guó)農(nóng)業(yè)科學(xué), 2017, 50(9): 1666-1673.doi: 10.3864/j.issn. 0578-1752.2017.09.012.

WANG X K, JIANG D, SUN Z Z. Study on phylogeny of 240 mandarin accessions with genotyping-by-sequencing technology. Scientia Agricultura Sinica, 2017, 50(9): 1666-1673. doi: 10.3864/ j.issn.0578-1752.2017.09.012. (in Chinese)

[14] CHEN W, HOU L, ZHANG Z Y, PANG X M, LI Y Y. Genetic diversity, population structure, and linkage disequilibrium of a core collection ofassessed with genome-wide SNPs developed by genotyping-by-sequencing and SSR markers. Frontiers in Plant Science, 2017, 8: 575.

[15] BIRD K A, AN H, GAZAVE E, GORE M A, PIRES J C, ROBERTSON L D, LABATE J A. Population structure and phylogenetic relationships in a diverse panel ofL. Frontiers in Plant Science, 2017, 8: 321.

[16] ELTAHER S, SALLAM A, BELAMKAR V, EMARA H A, NOWER A A, SALEM K F M, POLAND J, BAENZIGER P S. Genetic diversity and population structure of F3:6Nebraska winter wheat genotypes using genotyping-by-sequencing. Frontiers in Genetics, 2018, 9: 76.

[17] ALIPOUR H, BIHAMTA M R, MOHAMMADI V, PEYGHAMBARI S A, BAI G H, ZHANG G R. Genotyping-by-sequencing (GBS) revealed molecular genetic diversity of Iranian wheat landraces and cultivars. Frontiers in Plant Science, 2017, 8: 1293.

[18] DELFINI J, MODA-CIRINO V, DOS SANTOS NETO J, RUAS P M, SANT'ANA G C, GEPTS P, GON?ALVES L S A. Population structure, genetic diversity and genomic selection signatures among a Brazilian common bean germplasm. Scientific Reports, 2021, 11(1): 2964.

[19] LUO Z N, BROCK J, DYER J M, KUTCHAN T, SCHACHTMAN D, AUGUSTIN M, GE Y F, FAHLGREN N, ABDEL-HALEEM H. Genetic diversity and population structure of aspring panel. Frontiers in Plant Science, 2019, 10: 184.

[20] HYUN D Y, SEBASTIN R, LEE G A, LEE K J, KIM S H, YOO E, LEE S, KANG M J, LEE S B, JANG I, RO N Y, CHO G T. Genome-wide SNP markers for genotypic and phenotypic differentiation of melon (L.) varieties using genotyping- by-sequencing. International Journal of Molecular Sciences, 2021, 22(13): 6722.

[21] AKRAM S, ARIF M A R, HAMEED A. A GBS-based GWAS analysis of adaptability and yield traits in bread wheat (L.). Journal of Applied Genetics, 2021, 62(1): 27-41.

[22] POLAND J, ENDELMAN J, DAWSON J, RUTKOSKI J, WU S Y, MANES Y, DREISIGACKER S, CROSSA J, SáNCHEZ-VILLEDA H, SORRELLS M, JANNINK J L. Genomic selection in wheat breeding using genotyping-by-sequencing. The Plant Genome, 2012, 5(3): 103-113.

[23] KAUR G, PATHAK M, SINGLA D, SHARMA A, CHHUNEJA P, SARAO N K. High-density GBS-based genetic linkage map construction and QTL identification associated with yellow mosaic disease resistance in bitter gourd (L.). Frontiers in Plant Science, 2021, 12: 671620.

[24] ABED A, BADEA A, BEATTIE A, KHANAL R, TUCKER J, BELZILE F. A high-resolution consensus linkage map for barley based on GBS-derived genotypes. Genome, 2022, 65(2): 83-94.

[25] VERMA S, GUPTA S, BANDHIWAL N, KUMAR T, BHARADWAJ C, BHATIA S. High-density linkage map construction and mapping of seed trait QTLs in chickpea (L.) using Genotyping-by-Sequencing (GBS). Scientific Reports, 2015, 5(1): 17512.

[26] POLAND J A, BROWN P J, SORRELLS M E, JANNINK J L. Development of high-density genetic maps for barley and wheat using a novel two-enzyme genotyping-by-sequencing approach. PLoS ONE, 2012, 7(2): e32253.

[27] GAUR R, VERMA S, PRADHAN S, AMBREEN H, BHATIA S. A high-density SNP-based linkage map using genotyping-by-sequencing and its utilization for improved genome assembly of chickpea (L.). Functional & Integrative Genomics, 2020, 20(6): 763-773.

[28] WANG B Y, TAN H W, FANG W P, MEINHARDT L W, MISCHKE S, MATSUMOTO T, ZHANG D P. Developing single nucleotide polymorphism (SNP) markers from transcriptome sequences for identification of Longan () germplasm. Horticulture Research, 2015, 2: 14065.

[29] PEMBLETON L W, DRAYTON M C, BAIN M, BAILLIE R C, INCH C, SPANGENBERG G C, WANG J P, FORSTER J W, COGAN N O I. Targeted genotyping-by-sequencing permits cost- effective identification and discrimination of pasture grass species and cultivars. Theoretical and Applied Genetics, 2016, 129(5): 991-1005.

[30] STRAZZER P, SPELT C E, LI S J, BLIEK M, FEDERICI C T, ROOSE M L, KOES R, QUATTROCCHIO F M. Hyperacidification offruits by a vacuolar proton-pumping P-ATPase complex. Nature Communications, 2019, 10(1): 744.

[31] LI H, DURBIN R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics, 2009, 25(14): 1754-1760.

[32] DANECEK P, AUTON A, ABECASIS G, ALBERS C A, BANKS E, DEPRISTO M A, HANDSAKER R E, LUNTER G, MARTH G T, SHERRY S T, MCVEAN G, DURBIN R, GROUP 1 G P A. The variant call format and VCFtools. Bioinformatics, 2011, 27(15): 2156-2158.

[33] CHANG C C. Data management and summary statistics with PLINK. Methods in Molecular Biology, 2020, 2090: 49-65.

[34] MANICHAIKUL A, MYCHALECKYJ J C, RICH S S, DALY K, SALE M, CHEN W M. Robust relationship inference in genome-wide association studies. Bioinformatics, 2010, 26(22): 2867-2873.

[35] PAN T F, ALI M M, GONG J M, SHE W Q, PAN D M, GUO Z X, YU Y, CHEN F X. Fruit physiology and sugar-acid profile of 24 pomelo ((L.) osbeck) cultivars grown in subtropical region of China. Agronomy, 2021, 11(12): 2393.

[36] D'AGOSTINO N, TARANTO F, CAMPOSEO S, MANGINI G, FANELLI V, GADALETA S, MIAZZI M M, PAVAN S, DI RIENZO V, SABETTA W, LOMBARDO L, ZELASCO S, PERRI E, LOTTI C, CIANI E, MONTEMURRO C. GBS-derived SNP catalogue unveiled wide genetic variability and geographical relationships of Italian olive cultivars. Scientific Reports, 2018, 8(1): 15877.

[37] MEDINA R, WOGAN G O U, BI K, TERMIGNONI-GARCíA F, BERNAL M H, JARAMILLO-CORREA J P, WANG I J, VáZQUEZ- DOMíNGUEZ E. Phenotypic and genomic diversification with isolation by environment along elevational gradients in a neotropical treefrog. Molecular Ecology, 2021, 30(16): 4062-4076.

[38] 何天富. 中國(guó)柚類栽培. 北京: 中國(guó)農(nóng)業(yè)出版社, 1999.

HE T F. Cultivation of Pomelo in China. Beijing: China Agriculture Press, 1999. (in Chinese)

[39] BAAZAOUI I, MCEWAN J, ANDERSON R, BRAUNING R, MCCULLOCH A, VAN STIJN T, BEDHIAF-ROMDHANI S. GBS data identify pigmentation-specific genes of potential role in skin-photosensitization in two Tunisian sheep breeds. Animals, 2019, 10(1): 5.

[40] ZHANG M M, YANG L, SU Z C, ZHU M Z, LI W T, WU K L, DENG X M. Genome-wide scan and analysis of positive selective signatures in Dwarf Brown-egg Layers and Silky Fowl chickens. Poultry Science, 2017, 96(12): 4158-4171.

[41] NILO-POYANCO R, MORAGA C, BENEDETTO G, ORELLANA A, ALMEIDA A M. Shotgun proteomics of peach fruit reveals major metabolic pathways associated to ripening. BMC Genomics, 2021, 22(1): 17.

[42] ZHENG B B, ZHAO L, JIANG X H, CHERONO S, LIU J J, OGUTU C, NTINI C, ZHANG X J, HAN Y P. Assessment of organic acid accumulation and its related genes in peach. Food Chemistry, 2021, 334: 127567.

[43] YE J, WANG X, HU T X, ZHANG F X, WANG B, LI C X, YANG T X, LI H X, LU Y E, GIOVANNONI J J, ZHANG Y Y, YE Z B. An InDel in the promoter of Al-ACTIVATED MALATE TRANSPORTER9 selected during tomato domestication determines fruit malate contents and aluminum tolerance. The Plant Cell, 2017, 29(9): 2249-2268.

[44] LI C L, DOUGHERTY L, COLUCCIO A E, MENG D, EL-SHARKAWY I, BOREJSZA-WYSOCKA E, LIANG D, PI?EROS M A, XU K N, CHENG L L. Apple ALMT9 requires a conserved C-terminal domain for malate transport underlying fruit acidity. Plant Physiology, 2020, 182(2): 992-1006.

[45] DE ANGELI A, BAETZ U, FRANCISCO R, ZHANG J B, CHAVES M M, REGALADO A. The vacuolar channel VvALMT9 mediates malate and tartrate accumulation in berries of. Planta, 2013, 238(2): 283-291.

Population Genomic Structure of Pomelo Germplasm and Fruit Acidity Associated Genes Identification by Genotyping-by-Sequencing Technology

JIANG Dong1, WANG Xu1, LI RenJing1, ZHAO XiaoDong2, DAI XiangSheng2, Liu ZhengWei3

1Citrus Research Institute, Southwest University, Chongqing 400712;2Jinggangshan Agricultural Science and Technology Park Management Committee, Ji’an 343016, Jiangxi;3Jianggangshan University, Ji’an 343016, Jiangxi

【Objective】To reveal the phylogeny, population genetic structure and diversity level of pomelo ((L.) Osbeck) germplasm, and to efficiently utilize them to explore genes related to important fruit quality traits, this research provided an insight into the population genetic structure and phylogeny of pomelo germplasm and facilitated the pomelo varieties innovation.【Method】 A total of 282 pomelo accessions including landraces from different geographical regions and hybrid offspring of kiyomi tangor and pomelo were contained in this study. GBS library was constructed with genomic DNAs digested byR I restriction endonuclease and sequenced on Illumina HiSeq PE150 platform, the clean short reads were then mapped to pomelo reference genome by BWA, and SNPs were called out with SAMTOOLS pipeline. Based on 121 726 SNPs genotyping data, principal component analysis (PCA) and population genetic structure analysis were carried out and phylogenetic trees were constructed with Neighbor-joining method. Furthermore, two sub-populations containing high-acid accessions (32) and low-acid accessions (23) were used to identify candidate genes related to fruit acidity by Fst and XP-CLR selective sweeping analysis. Meanwhile, the genotype data of 282 pomelo accessions and the phenotypic data of titratable acid content in fruit were used for GWAS.【Result】A total of 201.66 Gb original reads were generated from 282 pomelo germplasm by GBS approach, in average each sample produced 0.72 Gb reads. After the screening conditions of sequencing depth of dp5, the miss less than 0.2 and minor alleles frequency (MAF)＞0.01, and a total of 121 726 SNPs were selected out for subsequent analysis. The PCA, structure and phylogenetic analysis all supported that the 282 pomelo germplasm could be divided into 6 subgroups, among which pomelo and kiyomi hybrid population, grapefruits and other pomelo hybrid populations could be obviously different from true-to-type pomelo populations, pomelos originated from different geographical region displayed unique genetic feature. The pomelo germplasm from Thailand and Vietnam formed a relatively unique group different from other domestic groups, such as ShaTian pomelo, Wen Dan pomelo, and Dian Jiang pomelo in China. The genetic introgression from Vietnam pomelos were exhibited in most pomelo germplasm in southern China, suggested that Vietnam was the origin center for pomelo. In addition, some pomelo germplasm with unknown origin have been identified accurately by GBS technology. This study showed that different geographical distribution and artificial selection pressure had great effect on the genomic composition of pomelo. Besides, Fst, XP-CLR selective sweeping analysis revealed a strong selection signal region on chromosome 7 contained genes annotated as dihydrolipoyl transacetylase (DLT-E2) of pyruvate dehydrogenase complex (PDC) and aluminum-activated malate transporter (ALMT9), which involved in the synthesis and transportation of citric acid. In additional GWAS genome-wide association analysis identified another region on chromosome 2, which was also highly associated with fruit acidity. 【Conclusion】GBS technology provided reliable and efficient method for studying the phylogeny and evolution of pomelo. The study showed that artificial cross breeding, long-term artificial selection, geography isolation and domestication were the major driving forces for the formation of different types of pomelo germplasm. In addition, it clearly showed that Southeast Asian was primary center for pomelo origin and China mainland was secondary evolutionary center. Several candidate genes related to citric acid content in pomelo fruits were identified by Fst, XP-CLR selective sweeping and GWAS. This study provided important gene resources for the further genetic improvement and breeding of pomelo fruits.

GBS; pomelo phylogeny; genes related to fruit acidity; GWAS

10.3864/j.issn.0578-1752.2023.08.010

2022-05-25；

2022-10-08

國(guó)家重點(diǎn)研發(fā)計(jì)劃（2018YFD1000101，2019YFD1001401）、種質(zhì)資源精準(zhǔn)鑒定（19211142）、江西科技計(jì)劃項(xiàng)目（20161BBF60048）

通信作者江東，Tel：13983194771；E-mail：jiangdong@cric.cn

（責(zé)任編輯 趙伶俐）