基于GBS技術開展柚類資源群體遺傳評價并發(fā)掘酸含量相關基因

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

基于GBS技術開展柚類資源群體遺傳評價并發(fā)掘酸含量相關基因

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

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

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

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

0 引言

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

1 材料與方法

1.1 試驗材料

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

1.2 試驗方法

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

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

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

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

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

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

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

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

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

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

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

2 結果

2.1 測序質量

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

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

2.2 主成分分析

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

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

2.3 群體遺傳結構分析

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

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

2.4 群體進化樹分析

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

2.5 群體的遺傳多樣性

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

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

圖3 利用282份柚類種質的121726 SNP采用鄰接法構建群體進化樹

表2 6個柚類亞群的遺傳多樣性

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

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

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

Fst最大值對應的SNP位于7號染色體上The corresponding maximum position of Fst and XP-CLR is on thesame location of chromosome 7

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

2.6 未知親本材料的鑒定

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

3 討論

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

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

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

3.2 越南柚在我國柚類起源演化中產生重要影響

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

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

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

4 結論

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

[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標記對中國柚類資源及近緣種遺傳多樣性研究. 農業(yè)生物技術學報, 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技術研究240份寬皮柑橘的系統(tǒng)演化. 中國農業(yè)科學, 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] 何天富. 中國柚類栽培. 北京: 中國農業(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

國家重點研發(fā)計劃（2018YFD1000101，2019YFD1001401）、種質資源精準鑒定（19211142）、江西科技計劃項目（20161BBF60048）

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

（責任編輯 趙伶俐）