蠶豆農藝性狀的SSR標記關聯分析

劉玉玲， 張紅巖，韓雪梅，周仙莉，侯萬偉

(1.青海大學，西寧 810016；2.青海省農林科學院，西寧 810016；3.農業農村部作物基因資源與種質創制青海科學觀測實驗站，西寧 810016；4.國家作物種質復份庫，西寧 810016)

【研究意義】蠶豆(ViciafabaL.)屬豆科(leguminous)野豌豆屬(ViciaL.)，一年生或越年生草本植物[1]。世界蠶豆種植主要集中在亞洲和非洲，種植面積和產量較多的國家主要有中國、埃塞俄比亞、摩洛哥等，在我國蠶豆主產區主要集中在長江流域各省和春播區的青海、甘肅[1-2]。蠶豆為糧食、蔬菜、飼料、綠肥兼用作物，具有豐富的營養和商業價值，在土壤結構改良、人類健康及保障畜牧業優質飼草(料)等方面均具有重要的作用，發展蠶豆生產，對于改善人們食物結構、維護農田生態平衡和振興經濟都有重要作用[3]。【前人研究進展】數量性狀基因座(Quantitative trait locus, QTL)是指控制數量性狀的基因在基因組中的位置，QTL定位的研究方法主要有連鎖定位分析和關聯分析兩種，其中以構建遺傳連鎖圖譜為基礎的連鎖定位分析是解析植物數量性狀遺傳的主要方法[4-5]，并在蠶豆中已被廣泛應用[6-7]。關聯分析(Association analysis)是一種以連鎖不平衡為基礎，可直接鑒定出與表型變異密切相關且具有特定功能的標記位點或基因位點[8]。相較于連鎖分析(Linkage analysis)，關聯分析具有精確度高，不需專門構建群體，且可同時分析同一位點的多個等位基因等優勢[9]。現已廣泛應用于玉米[10-11]、水稻[12-13]、棉花[14-15]、大麥[16-17]等眾多重要經濟作物的分子標記輔助育種研究中，并獲得了良好的效果。在蠶豆上，Ali等[18]利用單核苷酸多態性(SNP)和擴增片段長度多態性(AFLP)對蠶豆耐旱性和抗凍性進行了關聯分析。簡單重復序列(Simple sequence repeat, SSR)標記具有高度多態性、基因組特定性、重復性好等特點，是目前在蠶豆研究中應用最為廣泛的一類分子標記，如遺傳多樣性分析[19]、指紋圖譜構建[20]、QTL定位[21]等。【本研究切入點】目前，蠶豆的育種研究主要以傳統的選擇育種和雜交育種為主，而真正意義上的分子標記輔助育種研究開展甚少。初莢高度、單株莢數、株高、百粒重等農藝性狀是蠶豆傳統育種的主要選擇標準，利用SSR分子標記對這些農藝性狀進行關聯分析，尋找與之顯著關聯的標記，并將其應用于蠶豆分子標記輔助選擇，可極大提高選擇效率，加快定向育種進程。【擬解決的關鍵問題】本文以321份國內外蠶豆種質資源為材料，選用76對具有顯著多態性的SSR標記對供試群體進行基因分型，在群體遺傳結構分析的基礎上，尋找與蠶豆主要農藝性狀相關聯的SSR位點，并對其優異等位變異和種質進行發掘，以期為蠶豆種質資源的評價、分子標記輔助選擇育種、雜交組合選配和品種合理布局提供理論基礎。

1 材料與方法

1.1 試驗材料

321份蠶豆種質資源由中國農業科學院提供。國內資源214份，來源于19個省份；國外資源107份，來源于26個國家。

1.2 方法

1.2.1 田間試驗及農藝性狀的鑒定 2019年，321份蠶豆材料種植在青海農林科學院作物所試驗基地，每份材料種1行，每行10株，行距為35 cm，株距10 cm，按常規田間管理。在成熟收獲時隨機抽樣，每份材料取5株，分別測定初莢高度、初莢節位、成莢節數、生殖節數、有效籽粒數、有效分枝、有效莢數、單莢粒數、莢長、莢寬、籽粒面積、籽粒厚、籽粒周長、籽粒長、籽粒寬、株高、百粒重，最后取平均值。

1.2.2 基因組DNA的提取和SSR分析 田間取幼嫩葉片，采用DS法提取基因組DNA，利用超微量核酸蛋白濃度測定儀及瓊脂糖凝膠電泳法對DNA濃度和純度進行檢測。選用214對蠶豆引物對表型差異較大的12份供試材料進行多態性引物篩選，篩選出76對條帶差異清晰的SSR引物。PCR反應體系(10 μL)：1 μL模板DNA、5 μL 2×Mix、1.5 μL引物、2.5 μL ddH2O。PCR擴增程序：95 ℃預變性5 min；95 ℃變性30 s，48～55 ℃(根據引物而定)退火30 s，72 ℃延伸45 s，35個循環；72 ℃延伸10 min。擴增結束后，每個PCR管加入3 μL變性緩沖液(98%甲酰胺，10 mmol/L EDTA(pH 8.0)，0.25%溴酚藍，0.25%二甲苯青)，95 ℃變性8 min。擴增產物在預變性的6%變性聚丙烯酰胺凝膠上80 W恒功率電泳1.5 h，利用銀染法顯色，隨后統計條帶。

1.3 數據統計分析

采用人工讀膠的方法，按照帶型自下而上“1”～“n”記錄。各農藝性狀的相關性分析采用SPSS 26.0軟件。利用Power Marker 3.25軟件計算引物的多態性信息(PIC)。以Structure 2.3.4軟件進行群體遺傳結構的分析，估計最佳群體組群數K，其取值范圍為1～10，每個K值重復運行10次，參照Evanno等[22]提出的ΔK值方法確定合適的K值，并計算出Q參數。親緣關系系數(Kinship)用SPAGeDi軟件計算。關聯分析采用Tassel 2.1軟件，選擇一般線性模型(General linear model, GLM)和混合線性模型(Mixed linear model, MLM)進行分析。表型效應計算參照文自翔等[23]SSR位點等位變異表型效應計算法。

2 結果與分析

2.1 表型性狀統計分析

對321份蠶豆材料的17項表型性狀進行統計分析(表1)，變異系數最大的性狀是有效籽粒數(51.64%)，其他性狀變異系數較大的有成莢節數(46.13%)、有效分枝(39.26%)、初莢高度(37.25%)和有效莢數(30.51%)，籽粒厚變異系數最小(10.06%)。17個性狀變異系數的變幅為10.06%～51.64%，平均值為25.07%，這表明供試群體具有較為豐富的表型多樣性。除初莢高度和籽粒周長外，其他各性狀的偏度系數值和峰度系數值分別為0.00～0.95和-0.29～1.57，偏度系數和峰度系數絕對值均較小，說明321份蠶豆材料的17個性狀符合正態分布，可用于關聯分析。

表1 蠶豆材料農藝性狀的統計分析

續表1 Continued table 1

供試群體的17個表型性狀相關性分析結果表明(表2)，各性狀之間存在顯著或極顯著的正或負相關性。如構成蠶豆產量的幾個主要因素與其他性狀間的相關性，有效籽粒數與有效分枝、有效莢數、單莢粒數、株高呈極顯著正相關，與莢長呈顯著正相關，與莢寬呈極顯著負相關；有效莢數與單莢粒數和株高呈極顯著正相關；百粒重與籽粒面積、籽粒厚、籽粒周長、籽粒長和籽粒寬呈極顯著正相關。

表2 豆材料農藝性狀的相關性分析

2.2 SSR標記多態性分析

76對多態性顯著的SSR標記在321份蠶豆材料中共檢測到389個多態性位點，平均每個位點等位基因數為5.12個，變幅為2～11。多態性位點最豐富的是T19標記，有11個等位變異位點；SSR-1914、SSR-11967和EST-475標記均只檢測到2個等位變異位點。多態性信息量(PIC)的變幅為0.1903(CAASES-V2888)～0.7766(CAASES-V3031)，平均值為0.5156(表3)，PIC值大于平均數占比為50%，表明本研究選用標記基因多樣性較高。

表3 76對SSR標記多態性統計

續表3 Continued table 3

2.3 群體遺傳結構分析

通過76個SSR標記獲得的基因型數據，利用Structure 2.3.4軟件對321份蠶豆材料進行群體遺傳結構分析。由圖1可知，對數似然函數值lnP(D)隨亞群數增多而增大，且當K=2時，ΔK值最大，因此可將321份材料分成2個亞群。如圖2所示，橫坐標代表材料編號。通過Q值分析，發現其中有274份種質材料(占85.98%)的Q值≥0.6，被劃分在2個群組中，其中群組1包含103份材料(國內材料96份，國外材料7份，共占32.09%)；群組2包含171份材料(國內材料87份，國外材料84份，共占53.27%)，其余47份種質材料(國內材料31份，國外材料16份，共占14.64%)的Q值<0.6，沒有明確的群類歸屬特性，形成了一個混合群體。上述結果表明，該供試群體結構單一，這可有效地降低群體結構對關聯分析產生的影響。研究表明，群體結構的劃分情況不僅能對關聯分析結果的準確性產生影響，而且可在一定程度上反應亞群分類與材料特性的關系，該結果顯示亞群分類與材料的來源存在一定的對應關系，其中78.51%國外材料歸于群組2中。

A. K值與lnP(D)值的折線圖；B. K值與ΔK的折線圖A. Lines chart of K with lnP(D); B. Lines chart of K with ΔK圖1 K值與lnP(D)、ΔK值的折線圖Fig.1 Lines chart of K with lnP(D) and ΔK

對每一個個體，若紅色部分超整體的60%，則歸于POP1亞群；若綠色部分高于整體的60%，則歸于POP2亞群；若紅色或綠色介于整體的40%～60%，則歸于混合類群For each individual, if the red part is more than 60% of the whole, it belongs to POP1 subgroup; If the green part is more than 60% of the whole, it belongs to POP2 subgroup; If the red or green part is between 40% and 60% of the whole, it belongs to the mixed group圖2 321份蠶豆材料的群體結構Fig.2 The population structure of 321 faba bean materials

2.4 連鎖不平衡分析

基因間的連鎖不平衡是關聯分析的前提和基礎。76對SSR引物共獲得了2850種組合，其中R2≥0.01的組合占18.32%，在P<0.01的概率支持下，成對存在的不平衡組合有268個(圖3)，占總組合數的9.4%。從分析結果來看，本研究所選用的76對SSR引物在321份蠶豆材料中存在連鎖不平衡，可與農藝性狀進行關聯分析。

2.5 SSR標記與農藝性狀的關聯分析

在GLM_Q和MLM_Q+K模型中，分別以蠶豆材料的對應Q值和Q+K值作為協變量，將76個SSR分子標記與17個農藝性狀進行關聯分析，尋找與之相關聯的標記，并確定其解釋率(表5)。GLM_Q分析結果顯示，共檢測到43個與17個農藝性狀顯著(P<0.05)關聯的標記，解釋率為0.0146～0.1764；其中有16個標記分別關聯2個性狀，9個標記分別關聯3個性狀，3個標記分別關聯4個性狀，1個標記分別關聯5個性狀，2個標記分別關聯6個性狀；有效莢數和籽粒周長關聯的引物最多，均與9個標記關聯，解釋率分別為0.0146～0.1101和0.0322～0.0916。11個標記與12種農藝性狀極顯著(P<0.01)關聯，其中CAASES-V2698和CAASES-V2937標記均同時與有效籽粒數和有效莢數極顯著關聯，解釋率分別為0.0826、0.1101、0.0787和0.0595；CAASES-V2986標記同時與籽粒面積、籽粒厚、籽粒長、籽粒寬和百粒重極顯著關聯，解釋率分別為0.0707、0.0773、0.0714、0.0668和0.0644；SSR-12751標記同時與成莢節數、有效籽粒數和單莢粒數極顯著關聯，解釋率分別為0.0547、0.0474和0.0404。

MLM_Q+K分析結果顯示，共檢測到33個與17個農藝性狀顯著(P<0.05)關聯的標記，解釋率為0.0175～0.1706，且相校于GLM_Q分析結果出現了2個分別與有效分枝和籽粒面積新顯著關聯的標記(CAASES-V3058、SSR-11885)。總體上，MLM_Q+K分析檢測到與表型性狀關聯的SSR標記比GLM_Q分析結果少了10個，除解釋率不同外，剩余標記和性狀的關聯情況與GLM_Q分析結果基本相同。8個標記與8個農藝性狀極顯著(P<0.01)關聯，其中除CAASES-V2986不再與籽粒寬和百粒重以及SSR-12751不再與單莢粒數極顯著關聯外，CAASES-V2698、CAASES-V2937、CAASES-V2986和SSR-12751標記與性狀的關聯結果GLM_Q分析完全一致，且解釋率均有所下降。另外，CAASES-V2731、CAASES-V2745、CAASES-V3035和SSR-11448標記分別與株高、籽粒周長、單莢粒數和籽粒周長極顯著關聯，解釋率為0.0471、0.0583、0.0464和0.0897。

圖3 SSR標記之間的連鎖不平衡(LD) Fig.3 Linkage imbalance(LD) between SSR markers

表5 SSR標記與農藝性狀的關聯分析結果

續表5 Continued table 5

續表5 Continued table 5

2.6 優異等位變異發掘

為了獲得農藝性狀關聯標記位點的優異等位基因，分析了MLM_Q+K模型中達到極顯著(P<0.01)水平的等位變異表型效應及表型性狀。共發掘到增效的等位變異有26個，減效的等位變異53個(表6)。

由表6可知，成莢節數增效和減效效應最大的估計值分別為0.4966和-0.7459，對應的等位基因分別是SSR-12751-1和SSR-12751-4，其典型材料之一分別是湖北—青皮蠶豆和云南—金鐘豆；單莢粒數增效和減效效應最大的估計值分別為0.1376和-0.2164，對應的等位基因分別是CAASES-V3035-4和CAASES-V3035-1，其典型材料之一分別是江蘇—青皮和浙江—闊板青；有效莢數和有效籽粒數增效和減效效應最大的估計值分別為14.7732、35.7356和-5.5416、-12.768，對應的等位基因均分別是CAASES-V2937-4和CAASES-V2937-2，典型材料之一則均分別為甘肅—武山蠶豆和青海—青海五號；株高增效和減效效應最大的估計值分別為11.1656和-6.9988，對應的等位基因分別是CAASES-V2731-2和CAASES-V2731-3，其典型材料之一分別是內蒙古—當地樹豆和阿爾及利亞-372；籽粒厚增效和減效效應最大的估計值分別為0.3276和-0.2636，對應的等位基因分別是CAASES-V2986-4和CAASES-V2986-2，其典型材料之一分別是四川—金川本地胡豆和俄羅斯-zpyzue；籽粒長增效和減效效應最大的估計值分別為0.9784和-1.0496，對應的等位基因分別是CAASES-V2986-5和CAASES-V2986-3，其典型材料之一分別是日本—引 98-22和蘇丹-991；籽粒周長增效和減效效應最大的估計值分別為35.0313和-10.563，對應的等位基因分別是SSR-11448-5和SSR-11448-4，其典型材料之一分別是敘利亞-23705和青海—青海1號。

表6 優異等位變異發掘結果

續表6 Continued table 6

3 討 論

關聯分析在數量性狀研究中具有廣泛的應用前景，然而，關聯分析結果的準確性卻受多種因素的影響。其中，群體結構是影響關聯分析結果的重要因素，亞群的混合使整個群體的LD強度增強，可能導致不連鎖的基因多態性位點與性狀的關聯，從而得出假陽性結果，因此群體結構越簡單，關聯分析結果則更可靠[24]。由本研究群體遺傳結構分析結果可見，該供試材料可分為2個亞群，且僅14.64%的材料形成混合群體，群體結構相對單一，在一定程度上保證了關聯分析結果的準確性。另外，群體遺傳結構分析是基于亞群是否達到哈德溫伯格平衡的數學模型并計算Q值，不受人為因素對亞群劃分的影響，將所獲得的Q值作為協變量帶入回歸分析中，可有效提高關聯分析的可靠性[25-26]。但僅通過Q值作為矯正條件，卻并不能完全避免假陽性關聯的出現。張吉順等[27]的研究表明，以群體結構Q值和親緣系數值共同作為協變量的混合線性模型，即利用MLM_Q+K值模型進行關聯分析，更有利于提高關聯分析的準確性。鑒于此，本研究選用了GLM和MLM兩種模型進行關聯分析，相較于GLM模型，MLM模型關聯到的位點變少，但分析結果卻更為可靠。

各數量性狀之間常表現出某種特定的相關性，在本研究中，初莢高度與初莢節位、生殖節數、莢寬、籽粒周長、株高呈極顯著正相關，與成莢節數呈極顯著負相關。已有研究表明，這種表型性狀的相關性可能是由于某個基因的多效性或控制不同性狀的基因緊密連鎖造成，關聯分析中同一標記同時與多個相關性狀的關聯可解釋這種性狀之間的遺傳關聯性[28]。胡倩文等[29]在對大麥4個穗部性狀進行關聯分析時發現，小穗數與穗粒數呈極顯著正相關，并檢測到28個同時控制穗粒數和小穗數的關聯位點，從分子水平證實了小穗數與穗粒數的遺傳關聯性。本研究中，有效籽粒數和有效莢數之間呈極顯著相關，在MLM模型的顯著(P<0.05)條件下，共檢測到4個同時控制2種性狀的標記(CAASES-V2698、CAASES-V2937、SSR-11448和SSR-12751)；籽粒周長和百粒重之間呈極顯著相關，共檢測到2個標記同時控制2種性狀(CAASES-V2745和SSR-11885)；株高和有效分枝之間呈極顯著相關，共檢測到2個標記同時控制2種性狀(CAASES-V2698和CAASES-V3058)，可見相關性狀之間的遺傳關聯性在本研究中得到了進一步證實。在育種中利用這些關聯SSR標記進行輔助選擇，有助于實現多個性狀的協同改良[30-31]。

在本研究中，利用基于MLM模型極顯著(P<0.01)條件下的關聯分析結果，共發掘出7個分別對成莢節數、單莢粒數、有效莢數、有效籽粒數、株高、籽粒厚、籽粒長和籽粒周長增效潛力最大的優異等位變異(SSR-12751-1、CAASES-V3035-4、CAASES-V2937-4、CAASES-V2731-2、CAASES-V2986-4、CAASES-V2986-5和SSR-11448-5)和7份典型材料，CAASES-V2937-4等位基因對有效粒籽數的表型效應值高達35.7356，其中甘肅—武當蠶豆同時載有與有效莢數和有效籽粒數關聯的優異等位變異位點。通過對蠶豆產量構成因素的分析發現，構成蠶豆產量的主要因素有單株粒數、單株莢數、單株粒重、有效粒數、百粒重和有效莢數[32-33]。本研究的農藝性狀相關性分析結果顯示，籽粒長、籽粒周長和籽粒厚均與百粒重呈極顯著正相關，成莢節數和單莢粒數與有效粒數和有效莢數均呈極顯著正相關，株高與有效籽粒數呈極顯著相關，這6個農藝性狀可視為蠶豆產量的間接構成因素。鑒于此，本研究所挖掘到的增效等位變異可用于蠶豆產量的分子標記輔助選擇育種，攜帶增效基因的材料則可作為優異親本用于雜交育種。

由于目前有關蠶豆關聯分析的研究相對較少，本研究中極顯著關聯到的SSR標記，還沒有發現在其他文獻中報道，因此，還需進一步的分析來驗證本結果的準確性。因為影響關聯作圖準確性和精確性的因素還有群體大小、標記質量等，所以可進一步通過更高密度的標記，更大的群體規模以及多年多點的數據進一步提高關聯分析結果的可靠性。另外，本研究采用線性回歸方法進行性狀關聯分析，比較簡單直觀，但不能估計QTL的具體位置和加性、顯性或超顯性效應。相比之下，連鎖分析能檢測QTL的加性效應，且能克服關聯分析對稀有等位基因檢測效能低的缺點[34]。鑒于此，可采用連鎖分析和關聯分析相結合的策略，提高蠶豆農藝性狀QTL分析的效率和準確性。

4 結 論

321份蠶豆材料在17個農藝性狀間部分性狀存在顯著或極顯著的線性相關性，以MLM模型關聯分析共發現8個標記與8種農藝性狀極顯著(P<0.01)關聯，解釋率為0.0291～0.0897，基于這一關聯分析的結果，最終發掘出7個優異等位變異及7份優異種質，這為蠶豆育種親本的選擇和分子標記輔助選擇育種奠定了理論基礎。