柯會鋒，孫正文，王國寧，孟成生，吳立強

（1.河北農(nóng)業(yè)大學農(nóng)學院，華北作物改良與調(diào)控國家重點實驗室，河北 保定 071001；2.河北農(nóng)業(yè)大學生命科學學院，河北 保定 071001）

棉花（Gossypiumspp.）是重要的經(jīng)濟作物，棉纖維是重要的紡織工業(yè)原材料[1]；同時，棉籽也是棉花重要器官，因含蛋白質(zhì)、脂肪等營養(yǎng)物質(zhì)，在油脂、飼料加工和生物原料等領域具有重要價值[2-4]。然而，棉纖維是由胚珠表皮細胞分化而來，因而與棉籽發(fā)育之間存在一定競爭關系[5]。在當前棉花育種中，常注重纖維性狀遺傳改良，而忽略棉籽性狀，導致品種纖維產(chǎn)量不斷提高而棉籽質(zhì)量逐漸下降，影響發(fā)芽和出苗，同時伴隨棉苗長勢減弱、抗性降低，并最終影響纖維產(chǎn)量和品質(zhì)[6-7]。由此可見，在研究棉纖維性狀遺傳改良的同時，發(fā)掘控制棉籽性狀遺傳位點與基因?qū)τ趨f(xié)調(diào)棉纖維與棉籽競爭關系，實現(xiàn)二者同步遺傳改良意義重大。

目前，關于籽粒形狀大小遺傳位點和候選基因發(fā)掘在水稻、小麥、玉米等作物中已有報道[8-13]。Song等[14]獲得水稻粒寬和粒重數(shù)量性狀位 點（quantitative trait locus，QTL）候選基因GW2，該基因編碼E3泛素連接酶，其功能缺失可使水稻籽粒明顯變寬；程瑞如[15]利用小麥重組自交系（recombinant inbred line，RIL）群體定位到粒長、粒寬和粒厚QTLs，貢獻率9.3%~41.5%；耿慶河等[16]獲得4個在2種環(huán)境下控制菜豆粒長、粒厚和百粒重的一因多效QTLs；秦偉偉等[17]對玉米RIL群體5種環(huán)境下的籽粒形狀和大小進行鑒定，獲得30個QTLs，其中粒長、粒長/粒寬主效QTL在3種環(huán)境條件下共定位；耿青青[18]對大豆籽粒形狀大小進行關聯(lián)分析，獲得33、86和69個多環(huán)境條件下與粒長、粒寬和粒厚關聯(lián)SNPs。

與上述研究相比，目前關于棉籽形狀和大小的研究甚少。棉花籽指與種子活力指數(shù)以及棉苗鮮重正相關，且與籽粒油分和可溶性糖含量正相關，籽指越高，其發(fā)芽率越高，棉苗生長勢越強[19]。棉纖維產(chǎn)量性狀與種子品質(zhì)性狀呈現(xiàn)負相關，高產(chǎn)品種尤其高衣分品種通常籽指較低[20]。覃珊[6]認為，當前棉纖維產(chǎn)量的提高通常是以犧牲棉籽質(zhì)量為代價，因而使得棉纖維產(chǎn)量進一步提升受到限制。

針對上述棉籽性狀遺傳位點及候選基因研究較少，難以實現(xiàn)纖維、棉籽同步遺傳改良，本研究利用419份具有代表性棉花品種資源構(gòu)成的自然群體，分別在2年鑒定其棉籽性狀，結(jié)合群體重測序SNP基因型，發(fā)掘與棉籽形狀和大小關聯(lián)分子標記，并篩選其附近候選基因，為實現(xiàn)棉籽性狀分子遺傳改良提供標記和基因，并為今后實現(xiàn)棉纖維和棉籽性狀同步遺傳改良奠定基礎。

1 材料與方法

1.1 供試材料

以419份具有代表性棉花種質(zhì)資源構(gòu)成的自然群體為材料，該群體具有豐富遺傳變異，由河北農(nóng)業(yè)大學棉花遺傳育種課題組構(gòu)建并提供，前期曾利用該群體對棉花纖維品質(zhì)、產(chǎn)量性狀和耐低磷特性等開展關聯(lián)分析[21-22]。

1.2 研究方法

1.2.1 供試棉花自然群體田間種植與管理供試棉花自然群體分別于2015和2017年種植于海南省三亞市崖州區(qū)南濱農(nóng)場，中性土壤，完全隨機區(qū)組試驗設計，3次重復，行長6.00 m，行距0.70 m，株距0.25 m，密度5.7萬株·hm-2，10月25日足墑起壟覆膜播種（膜寬1.1 m，膜厚0.012 mm），播種前施300 kg·hm-2基肥（復合肥，氮磷鉀比例15∶15∶15，雅苒商貿(mào)國際有限公司），開花期追施150 kg·hm-2氮肥（氮≥46%，河北正元氫能科技有限公司），棉花生長期間澆水4次，12月20日人工打頂，翌年3月20日左右收獲。

1.2.2供試自然群體棉籽相關性狀測定待供試棉花自然群體成熟后進行人工收獲，經(jīng)軋花、脫絨獲得棉籽；選取健康、飽滿、整齊一致的棉籽100粒，利用SC-G種子自動檢測系統(tǒng)（杭州萬深檢測有限公司）對棉籽進行拍照；隨后利用該系統(tǒng)對圖像進行分析，具體性狀包括棉籽的粒長（seed length,SL）、粒寬（seed width,SW）、籽粒面積（seed area,SA）、籽粒周長（seed perimeter,SP）和籽粒長寬比（ratio of seed length to width,RSLW）[23-24]。

1.2.3 供試自然群體棉籽性狀數(shù)據(jù)分析采用Microsoft Excel 2017軟件對供試群體棉籽相關性狀的平均值、標準差、變異系數(shù)、最大值、最小值、偏度系數(shù)和峰度系數(shù)等進行描述統(tǒng)計分析；同時采用該軟件計算上述性狀間的相關系數(shù)，并采用R語言corrplot軟件繪制性狀間相關系數(shù)的熱圖。

1.2.4 供試群體棉籽性狀關聯(lián)分析與候選基因篩選利用供試群體5個棉籽相關性狀，結(jié)合群體已有重測序SNP基因型（含3 665 030個高質(zhì)量SNPs，最小等位基因頻率≥0.05，缺失率≤0.2），采用GEMMA（genome-wide efficient mixed-model association）軟件進行全基因組關聯(lián)分析，利用前3個主成分建立S矩陣對群體結(jié)構(gòu)進行校正，利用簡單匹配系數(shù)矩陣建立K矩陣，采用R軟件的LDheatmap實現(xiàn)LD block可視化[22]。在此基礎上，利用該群體進行關聯(lián)分析以及候選基因挖掘[21]，在關聯(lián)標記附近50 kb范圍內(nèi)篩選可能與棉籽大小和形狀相關的候選基因，其中基因的注釋參考TM-1基因組[25]。隨后，利用已公布的棉花胚珠不同發(fā)育時期基因轉(zhuǎn)錄組數(shù)據(jù)（品種為Texas Marker-1，基因表達量用FPKM值表示，取樣時間包括棉花開花后0、1、5、10、20、25和35 d），分析上述候選基因在胚珠中的表達量，并采用HemI 1.0軟件繪制候選基因表達量的熱圖，制作熱圖時，候選基因表達量采用log2(FPKM+1)方法進行數(shù)據(jù)轉(zhuǎn)換[25]。

2 結(jié)果與分析

2.1 自然群體棉籽性狀遺傳變異與相關性分析

2.1.1 棉籽性狀遺傳變異分析通過分析供試自然群體棉籽大小和形狀相關性狀（表1，圖1）發(fā)現(xiàn)，棉籽的粒長平均8.46 mm，粒寬平均4.67 mm，長寬比平均1.82，面積與周長平均29.79 mm2和22.12 mm；同時發(fā)現(xiàn)，5個性狀在供試群體中存在較大遺傳變異，其變異系數(shù)為3.91%~9.55%。另外，棉籽性狀在2年度間均呈正態(tài)分布，年度間分布趨勢基本一致，表明棉籽性狀由多基因控制。

圖1 供試棉花品種的棉籽相關性狀次數(shù)分布Fig.1 Distributions of seed related traits of cotton varieties

表1 自然群體棉籽相關性狀遺傳變異Table 1 Genetic variation of seed related traits in cotton natural population

2.1.2 相關性分析供試自然群體5個棉籽性狀間（圖2）存在極顯著相關性，并以粒長與面積、周長，粒寬與面積、周長，面積與周長間的相關系數(shù)更高。同時發(fā)現(xiàn)，同一性狀在不同年度間的相關系數(shù)較高，粒長、粒寬、長寬比、面積和周長在年度間的相 關系數(shù) 分別為0.80、0.84、0.81、0.83和0.81，暗示可獲得不同年度間共同關聯(lián)的分子標記。

圖2 供試棉花自然群體棉籽性狀相關性Fig.2 Correlation of seed related traits in cotton natural population

2.2 自然群體棉籽性狀關聯(lián)位點及單倍型分析

2.2.1 棉籽性狀關聯(lián)SNPs分析通過分析供試自然群體棉籽相關性狀關聯(lián)SNPs（表2）發(fā)現(xiàn)，7個SNPs同時在2年與棉籽大小和形狀關聯(lián)，分別位于D11、D12、D05和A04染色體；23個SNPs在單年與3個以上棉籽性狀關聯(lián)，分別位于A07、A08、A01、A02、A05、A06和D08染色體；其中D11染色體SNP標記（D11∶50561153）同時在2年與3個相關性狀（粒寬、面積和周長）關聯(lián)，D12染色體SNP標記（D12∶56031535）在2年與2個相關性狀（粒長和周長）關聯(lián)，并在單年與籽粒面積關聯(lián)。D05染 色 體4個SNPs（D05∶1490099、D05∶1490516、D05∶1496927和D05∶1512276，物理距離約22 kb）在2年同時與籽粒的長寬比顯著關聯(lián)；A04染色體SNP標記A04∶54398482在2年與籽粒周長關聯(lián)，并在單年與籽粒面積顯著關聯(lián)。A07染色體12個SNPs與棉籽性狀相關聯(lián)，其中10個SNPs標記（物理距離約28 kb）在單年同時與粒長、籽粒面積和周長關聯(lián)。

表2 供試自然群體棉籽大小和形狀關聯(lián)SNP標記Table 2 Associated SNP markers of seed size and shape in cotton natural population 續(xù)表Continued

表2 供試自然群體棉籽大小和形狀關聯(lián)SNP標記Table 2 Associated SNP markers of seed size and shape in cotton natural population

2.2.2 A07染色體棉籽性狀關聯(lián)SNPs單倍型分析進一步對A07染色體關聯(lián)SNPs進行單倍型分析（圖3）發(fā)現(xiàn)，在該染色體的70.68~70.70 Mb范圍內(nèi)，與棉籽性狀關聯(lián)的9個SNPs在不同品種間呈現(xiàn)出2種單倍型Hap1和Hap2，其中Hap1包含28份材料，Hap2包含186份材料；進一步對兩種單倍型之間的籽粒長度、周長和面積進行分析發(fā)現(xiàn)，均存在極顯著差異，并以單倍型Hap1的各相關性狀表現(xiàn)更優(yōu)，為揭示棉籽大小和形狀遺傳基礎提供了依據(jù)，同時也為今后實現(xiàn)纖維和棉籽性狀同步遺傳改良奠定了基礎。

圖3 棉花A07染色體棉籽性狀關聯(lián)SNPs單倍型分析Fig.3 Associated SNPs and haplotypes on chromosome A07 for cotton seed related traits

2.3 棉籽大小和形狀候選基因篩選

2.3.1 棉花胚珠不同發(fā)育時期候選基因綜合分析通過分析棉籽大小和形狀關聯(lián)SNPs附近候選基因，結(jié)合其在棉花胚珠不同發(fā)育時期表達模式（圖4）發(fā)現(xiàn)，有21個候選基因隨棉花胚珠發(fā)育進程表達量增加，其中包括4個在胚珠發(fā)育前期（開花后1~10 d）表達量增加的候選基因、8個在胚珠發(fā)育后期（開花后20~35 d）表達量增加的基因和9個隨胚珠整個生長發(fā)育期（開花后1~35 d）表達量增加的候選基因。

圖4 棉花胚珠發(fā)育不同時期關聯(lián)SNPs附近候選基因的表達Fig.4 Expressions of candidate genes associated to the adjacent SNPs at different developmental stages in cotton ovule

2.3.2 候選基因的篩選 進一步分析上述候選基因（表3）發(fā)現(xiàn)，D05染色體候選基因Gh_D05G0148編碼蛋白為EBF蛋白（EIN3-binding F box protein），該蛋白已被證實可與乙烯信號調(diào)控因子EIN3（ethylene insensitive 3）發(fā)生互作，通過參與乙烯信號的調(diào)控反應進而影響植物生長發(fā)育[26]；分析Gh_D05G0148在棉花胚珠不同發(fā)育時期轉(zhuǎn)錄組數(shù)據(jù)發(fā)現(xiàn)，與0 d相比，該基因在開花后1~35 d始終保持較高表達水平，其不同時期表達量分別較0 d高出3.1~14.7倍，并以開花后10 d表達量最高，說明該基因?qū)γ藁ǚN子生長發(fā)育有重要作用。

表3 棉籽大小和形狀相關候選基因Table 3 Candidate genes related to the seed size and shape in cotton

同時發(fā)現(xiàn)，D05染色體Gh_D05G0144編碼蛋白屬于YABBY轉(zhuǎn)錄因子家族，該類轉(zhuǎn)錄因子是種子植物特有的一類鋅指蛋白超家族，已被證實參與植物多種組織和器官生長發(fā)育[27]；分析Gh_D05G0144在棉花胚珠不同發(fā)育時期表達發(fā)現(xiàn)，該基因在棉花開花后20~35 d表達量升高，推測其可能與棉籽后期發(fā)育密切相關。候選基因Gh_D05G0146編碼蛋白屬于三/四氨基酸重復基序（tetratricopeptide repeat，TPR）家族蛋白，該類蛋白也參與植物組織器官生長發(fā)育，且與生長素極性運輸密切相關[28]；分析Gh_D05G0146表達發(fā)現(xiàn)，該基因在棉花開花后1~3和25~35 d表達量升高，推測其在棉籽生長發(fā)育中具有重要作用。

另外，A08染色體候選基因Gh_A08G0768屬于

蛋白激酶家族基因，目前已有很多試驗證明蛋白激酶在植物生長發(fā)育過程中具有重要作用，諸如植物細胞的伸長與擴張、器官發(fā)生與形態(tài)建成、果實生長發(fā)育以及多種植物激素反應等[29]；分析Gh_A08G0768表達量發(fā)現(xiàn)，其在棉花開花后10~35 d的胚珠中的表達量分別為0 d的2.7~3.7倍，可能參與棉花種子生長發(fā)育。由此可見，上述候選基因為進一步開展棉籽發(fā)育相關基因功能鑒定奠定了基礎。

3 討論

有學者利用分離群體或自然群體，對棉籽的籽指、種仁重量、種殼重量、蛋白質(zhì)、脂肪和氨基酸含量等開展研究。Wang等[30]通過鑒定棉花RIL群體3種環(huán)境條件下的籽指、百粒種仁重、種仁長度和寬度等性狀，結(jié)合群體SLAF-seq SNP遺傳連鎖圖譜，獲得6個種仁長度QTLs，位于A05、A11、A12、D02、D08和D13染 色 體，貢 獻 率11.2%~20.8%，其中3個QTLs為多環(huán)境檢測；同時獲得6個種仁寬度QTLs，位于A01、A08、A11、D01、D03和D09染色體，貢獻率10.7%~56.9%，其中3個QTLs為多環(huán)境檢測。本研究利用構(gòu)建的棉花自然群體，分別在2年鑒定其粒長、粒寬、長寬比、面積和周長等性狀，結(jié)合群體重測序SNP基因型，獲得30個可在2年間共同關聯(lián)或1年與3個以上性狀關聯(lián)SNPs，分別位于棉花11條染色體，其中7個SNPs標記位于D11、D12、D05和A04染色體，可在2年間共同檢測到。

同時，本研究在A07染色體檢測到10個相鄰SNPs與籽粒長度、面積和周長顯著關聯(lián)，其中9個SNPs在課題組之前利用該群體進行棉花籽指關聯(lián)分析中被檢測到[21]。進一步分析供試群體的籽指與籽粒長度、寬度、長寬比、面積和周長性狀間的相關關系，結(jié)果發(fā)現(xiàn)，除長寬比以外，其他4個籽粒性狀與籽指之間的相關系數(shù)均達到極顯著水平（r0.01=0.13），其中籽指與籽粒長度、寬度、面積和周長在2015年條件下的相關系數(shù)分別為0.79、0.89、0.92和0.85，在2017年條件下的相關系數(shù)則分別為0.79、0.87、0.89和0.84。并且，Wang等[30]研究也發(fā)現(xiàn)，棉花的籽指與棉仁重量、長度和寬度之間的相關系數(shù)亦達到顯著或極顯著水平，其中籽指與棉仁長度的相關系數(shù)為0.21~0.41，籽指與棉仁寬度之間的相關系數(shù)為0.64~0.79。由此可見，因棉花籽指與粒長、粒寬、面積和周長之間存在顯著相關性，故本研究一方面證實了結(jié)果的可靠性，另一方面為揭示性狀相關關系的分子機制奠定了理論基礎。另外，通過比較本研究與Wang等[30]報道的種仁長度和寬度QTLs發(fā)現(xiàn)，在A01、A08、A05和D08染色體存在控制棉籽形狀和大小遺傳位點，但因Wang等[29]所用遺傳圖譜屬于SLAF-seq SNP標記，未給出標記具體物理位置，因而無法進一步比較本研究關聯(lián)標記與該文獻QTL連鎖標記之間的物理距離。

基于上述關聯(lián)SNPs標記，本研究在其附近篩選到多個與棉籽形狀和大小相關的候選基因，在這些基因中，D05染色體候選基因Gh_D05G0144編碼種子植物特有轉(zhuǎn)錄因子YABBY家族蛋白，該家族屬于鋅指蛋白超家族，已報道與植物葉片、花和果實等生長發(fā)育密切相關[27]。擬南芥基因組中含有6個YABBY基因，其中4個在葉片及其葉源器官（子葉、萼片等）表達，另外2個則只在特異的花器官中表達，如INNER NO OUTER（INO）基因只在胚珠外珠被的最外層細胞層特異表達，并具有促進外珠被生長發(fā)育的功能[31]。最近，Sun等[32]研究發(fā)現(xiàn)，擬南芥YABBY轉(zhuǎn)錄因子INO可通過抑制鐵轉(zhuǎn)運蛋白基因（natural resistance associated macrophage protein 1,NRAMP1）的表達，從而減少擬南芥幼嫩種子中的鐵元素積累，進而保證植株的生殖生長和種子發(fā)育。

同時，候選基因Gh_D05G0148可編碼EBF蛋白，該類蛋白通常作為乙烯信號調(diào)控因子EIN3的互作蛋白來調(diào)控植物生長發(fā)育。Guo等[26]研究發(fā)現(xiàn)，EBF蛋白可通過介導EIN3/EIL蛋白的降解實現(xiàn)調(diào)控乙烯信號通路的功能，番茄中的EBF蛋白基因SlEBF1和SlEBF2通過乙烯信號途徑影響植株衰老過程以及果實的生長發(fā)育與成熟，并且，番茄EBF基因SlEBF2-like的上游啟動子區(qū)含有3個乙烯響應元件，在番茄中過表達SlEBF2-like不僅增加了番茄果實的長度，而且改變了果實形狀。本研究通過分析EBF蛋白基因Gh_D05G0148在棉花胚珠中的表達量發(fā)現(xiàn)，以開花后10 d表達量最高，并在開花后1~35 d始終保持較高水平，暗示其可能參與棉花種子的生長發(fā)育。

另外，候選基因Gh_D05G0146編碼含TPR結(jié)構(gòu)的家族蛋白，該類蛋白已被發(fā)現(xiàn)可參與擬南芥根系生長發(fā)育，且與生長素極性運輸密切相關，擬南芥TPR家族基因SSR1（short and swollen root 1）突變后，因細胞伸長和增殖嚴重受損，導致初生根生長被明顯抑制，且根系生長素水平顯著下降[28]。在本研究中，Gh_D05G0146隨棉花胚珠發(fā)育表達量升高，因而推測其在棉花種子生長發(fā)育過程中發(fā)揮一定作用。

綜上，本研究發(fā)掘的關聯(lián)SNP標記以及候選基因，為今后開展棉籽大小和形狀分子遺傳改良提供了依據(jù)，也為進一步解析棉花籽粒大小和形狀分子機理奠定了基礎。