棉籽大小与形状关联标记发掘及候选基因筛选

柯会锋，孙正文，王国宁，孟成生，吴立强

（1.河北农业大学农学院，华北作物改良与调控国家重点实验室，河北 保定 071001；2.河北农业大学生命科学学院，河北 保定 071001）

棉花（Gossypiumspp.）是重要的经济作物，棉纤维是重要的纺织工业原材料[1]；同时，棉籽也是棉花重要器官，因含蛋白质、脂肪等营养物质，在油脂、饲料加工和生物原料等领域具有重要价值[2-4]。然而，棉纤维是由胚珠表皮细胞分化而来，因而与棉籽发育之间存在一定竞争关系[5]。在当前棉花育种中，常注重纤维性状遗传改良，而忽略棉籽性状，导致品种纤维产量不断提高而棉籽质量逐渐下降，影响发芽和出苗，同时伴随棉苗长势减弱、抗性降低，并最终影响纤维产量和品质[6-7]。由此可见，在研究棉纤维性状遗传改良的同时，发掘控制棉籽性状遗传位点与基因对于协调棉纤维与棉籽竞争关系，实现二者同步遗传改良意义重大。

目前，关于籽粒形状大小遗传位点和候选基因发掘在水稻、小麦、玉米等作物中已有报道[8-13]。Song等[14]获得水稻粒宽和粒重数量性状位 点（quantitative trait locus，QTL）候选基因GW2，该基因编码E3泛素连接酶，其功能缺失可使水稻籽粒明显变宽；程瑞如[15]利用小麦重组自交系（recombinant inbred line，RIL）群体定位到粒长、粒宽和粒厚QTLs，贡献率9.3%~41.5%；耿庆河等[16]获得4个在2种环境下控制菜豆粒长、粒厚和百粒重的一因多效QTLs；秦伟伟等[17]对玉米RIL群体5种环境下的籽粒形状和大小进行鉴定，获得30个QTLs，其中粒长、粒长/粒宽主效QTL在3种环境条件下共定位；耿青青[18]对大豆籽粒形状大小进行关联分析，获得33、86和69个多环境条件下与粒长、粒宽和粒厚关联SNPs。

与上述研究相比，目前关于棉籽形状和大小的研究甚少。棉花籽指与种子活力指数以及棉苗鲜重正相关，且与籽粒油分和可溶性糖含量正相关，籽指越高，其发芽率越高，棉苗生长势越强[19]。棉纤维产量性状与种子品质性状呈现负相关，高产品种尤其高衣分品种通常籽指较低[20]。覃珊[6]认为，当前棉纤维产量的提高通常是以牺牲棉籽质量为代价，因而使得棉纤维产量进一步提升受到限制。

针对上述棉籽性状遗传位点及候选基因研究较少，难以实现纤维、棉籽同步遗传改良，本研究利用419份具有代表性棉花品种资源构成的自然群体，分别在2年鉴定其棉籽性状，结合群体重测序SNP基因型，发掘与棉籽形状和大小关联分子标记，并筛选其附近候选基因，为实现棉籽性状分子遗传改良提供标记和基因，并为今后实现棉纤维和棉籽性状同步遗传改良奠定基础。

1 材料与方法

1.1 供试材料

以419份具有代表性棉花种质资源构成的自然群体为材料，该群体具有丰富遗传变异，由河北农业大学棉花遗传育种课题组构建并提供，前期曾利用该群体对棉花纤维品质、产量性状和耐低磷特性等开展关联分析[21-22]。

1.2 研究方法

1.2.1 供试棉花自然群体田间种植与管理供试棉花自然群体分别于2015和2017年种植于海南省三亚市崖州区南滨农场，中性土壤，完全随机区组试验设计，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，雅苒商贸国际有限公司），开花期追施150 kg·hm-2氮肥（氮≥46%，河北正元氢能科技有限公司），棉花生长期间浇水4次，12月20日人工打顶，翌年3月20日左右收获。

1.2.2供试自然群体棉籽相关性状测定待供试棉花自然群体成熟后进行人工收获，经轧花、脱绒获得棉籽；选取健康、饱满、整齐一致的棉籽100粒，利用SC-G种子自动检测系统（杭州万深检测有限公司）对棉籽进行拍照；随后利用该系统对图像进行分析，具体性状包括棉籽的粒长（seed length,SL）、粒宽（seed width,SW）、籽粒面积（seed area,SA）、籽粒周长（seed perimeter,SP）和籽粒长宽比（ratio of seed length to width,RSLW）[23-24]。

1.2.3 供试自然群体棉籽性状数据分析采用Microsoft Excel 2017软件对供试群体棉籽相关性状的平均值、标准差、变异系数、最大值、最小值、偏度系数和峰度系数等进行描述统计分析；同时采用该软件计算上述性状间的相关系数，并采用R语言corrplot软件绘制性状间相关系数的热图。

1.2.4 供试群体棉籽性状关联分析与候选基因筛选利用供试群体5个棉籽相关性状，结合群体已有重测序SNP基因型（含3 665 030个高质量SNPs，最小等位基因频率≥0.05，缺失率≤0.2），采用GEMMA（genome-wide efficient mixed-model association）软件进行全基因组关联分析，利用前3个主成分建立S矩阵对群体结构进行校正，利用简单匹配系数矩阵建立K矩阵，采用R软件的LDheatmap实现LD block可视化[22]。在此基础上，利用该群体进行关联分析以及候选基因挖掘[21]，在关联标记附近50 kb范围内筛选可能与棉籽大小和形状相关的候选基因，其中基因的注释参考TM-1基因组[25]。随后，利用已公布的棉花胚珠不同发育时期基因转录组数据（品种为Texas Marker-1，基因表达量用FPKM值表示，取样时间包括棉花开花后0、1、5、10、20、25和35 d），分析上述候选基因在胚珠中的表达量，并采用HemI 1.0软件绘制候选基因表达量的热图，制作热图时，候选基因表达量采用log2(FPKM+1)方法进行数据转换[25]。

2 结果与分析

2.1 自然群体棉籽性状遗传变异与相关性分析

2.1.1 棉籽性状遗传变异分析通过分析供试自然群体棉籽大小和形状相关性状（表1，图1）发现，棉籽的粒长平均8.46 mm，粒宽平均4.67 mm，长宽比平均1.82，面积与周长平均29.79 mm2和22.12 mm；同时发现，5个性状在供试群体中存在较大遗传变异，其变异系数为3.91%~9.55%。另外，棉籽性状在2年度间均呈正态分布，年度间分布趋势基本一致，表明棉籽性状由多基因控制。

图1 供试棉花品种的棉籽相关性状次数分布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）存在极显著相关性，并以粒长与面积、周长，粒宽与面积、周长，面积与周长间的相关系数更高。同时发现，同一性状在不同年度间的相关系数较高，粒长、粒宽、长宽比、面积和周长在年度间的相 关系数 分别为0.80、0.84、0.81、0.83和0.81，暗示可获得不同年度间共同关联的分子标记。

图2 供试棉花自然群体棉籽性状相关性Fig.2 Correlation of seed related traits in cotton natural population

2.2 自然群体棉籽性状关联位点及单倍型分析

2.2.1 棉籽性状关联SNPs分析通过分析供试自然群体棉籽相关性状关联SNPs（表2）发现，7个SNPs同时在2年与棉籽大小和形状关联，分别位于D11、D12、D05和A04染色体；23个SNPs在单年与3个以上棉籽性状关联，分别位于A07、A08、A01、A02、A05、A06和D08染色体；其中D11染色体SNP标记（D11∶50561153）同时在2年与3个相关性状（粒宽、面积和周长）关联，D12染色体SNP标记（D12∶56031535）在2年与2个相关性状（粒长和周长）关联，并在单年与籽粒面积关联。D05染 色 体4个SNPs（D05∶1490099、D05∶1490516、D05∶1496927和D05∶1512276，物理距离约22 kb）在2年同时与籽粒的长宽比显著关联；A04染色体SNP标记A04∶54398482在2年与籽粒周长关联，并在单年与籽粒面积显著关联。A07染色体12个SNPs与棉籽性状相关联，其中10个SNPs标记（物理距离约28 kb）在单年同时与粒长、籽粒面积和周长关联。

表2 供试自然群体棉籽大小和形状关联SNP标记Table 2 Associated SNP markers of seed size and shape in cotton natural population 续表Continued

表2 供试自然群体棉籽大小和形状关联SNP标记Table 2 Associated SNP markers of seed size and shape in cotton natural population

2.2.2 A07染色体棉籽性状关联SNPs单倍型分析进一步对A07染色体关联SNPs进行单倍型分析（图3）发现，在该染色体的70.68~70.70 Mb范围内，与棉籽性状关联的9个SNPs在不同品种间呈现出2种单倍型Hap1和Hap2，其中Hap1包含28份材料，Hap2包含186份材料；进一步对两种单倍型之间的籽粒长度、周长和面积进行分析发现，均存在极显著差异，并以单倍型Hap1的各相关性状表现更优，为揭示棉籽大小和形状遗传基础提供了依据，同时也为今后实现纤维和棉籽性状同步遗传改良奠定了基础。

图3 棉花A07染色体棉籽性状关联SNPs单倍型分析Fig.3 Associated SNPs and haplotypes on chromosome A07 for cotton seed related traits

2.3 棉籽大小和形状候选基因筛选

2.3.1 棉花胚珠不同发育时期候选基因综合分析通过分析棉籽大小和形状关联SNPs附近候选基因，结合其在棉花胚珠不同发育时期表达模式（图4）发现，有21个候选基因随棉花胚珠发育进程表达量增加，其中包括4个在胚珠发育前期（开花后1~10 d）表达量增加的候选基因、8个在胚珠发育后期（开花后20~35 d）表达量增加的基因和9个随胚珠整个生长发育期（开花后1~35 d）表达量增加的候选基因。

图4 棉花胚珠发育不同时期关联SNPs附近候选基因的表达Fig.4 Expressions of candidate genes associated to the adjacent SNPs at different developmental stages in cotton ovule

2.3.2 候选基因的筛选 进一步分析上述候选基因（表3）发现，D05染色体候选基因Gh_D05G0148编码蛋白为EBF蛋白（EIN3-binding F box protein），该蛋白已被证实可与乙烯信号调控因子EIN3（ethylene insensitive 3）发生互作，通过参与乙烯信号的调控反应进而影响植物生长发育[26]；分析Gh_D05G0148在棉花胚珠不同发育时期转录组数据发现，与0 d相比，该基因在开花后1~35 d始终保持较高表达水平，其不同时期表达量分别较0 d高出3.1~14.7倍，并以开花后10 d表达量最高，说明该基因对棉花种子生长发育有重要作用。

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

同时发现，D05染色体Gh_D05G0144编码蛋白属于YABBY转录因子家族，该类转录因子是种子植物特有的一类锌指蛋白超家族，已被证实参与植物多种组织和器官生长发育[27]；分析Gh_D05G0144在棉花胚珠不同发育时期表达发现，该基因在棉花开花后20~35 d表达量升高，推测其可能与棉籽后期发育密切相关。候选基因Gh_D05G0146编码蛋白属于三/四氨基酸重复基序（tetratricopeptide repeat，TPR）家族蛋白，该类蛋白也参与植物组织器官生长发育，且与生长素极性运输密切相关[28]；分析Gh_D05G0146表达发现，该基因在棉花开花后1~3和25~35 d表达量升高，推测其在棉籽生长发育中具有重要作用。

另外，A08染色体候选基因Gh_A08G0768属于

蛋白激酶家族基因，目前已有很多试验证明蛋白激酶在植物生长发育过程中具有重要作用，诸如植物细胞的伸长与扩张、器官发生与形态建成、果实生长发育以及多种植物激素反应等[29]；分析Gh_A08G0768表达量发现，其在棉花开花后10~35 d的胚珠中的表达量分别为0 d的2.7~3.7倍，可能参与棉花种子生长发育。由此可见，上述候选基因为进一步开展棉籽发育相关基因功能鉴定奠定了基础。

3 讨论

有学者利用分离群体或自然群体，对棉籽的籽指、种仁重量、种壳重量、蛋白质、脂肪和氨基酸含量等开展研究。Wang等[30]通过鉴定棉花RIL群体3种环境条件下的籽指、百粒种仁重、种仁长度和宽度等性状，结合群体SLAF-seq SNP遗传连锁图谱，获得6个种仁长度QTLs，位于A05、A11、A12、D02、D08和D13染 色 体，贡 献 率11.2%~20.8%，其中3个QTLs为多环境检测；同时获得6个种仁宽度QTLs，位于A01、A08、A11、D01、D03和D09染色体，贡献率10.7%~56.9%，其中3个QTLs为多环境检测。本研究利用构建的棉花自然群体，分别在2年鉴定其粒长、粒宽、长宽比、面积和周长等性状，结合群体重测序SNP基因型，获得30个可在2年间共同关联或1年与3个以上性状关联SNPs，分别位于棉花11条染色体，其中7个SNPs标记位于D11、D12、D05和A04染色体，可在2年间共同检测到。

同时，本研究在A07染色体检测到10个相邻SNPs与籽粒长度、面积和周长显著关联，其中9个SNPs在课题组之前利用该群体进行棉花籽指关联分析中被检测到[21]。进一步分析供试群体的籽指与籽粒长度、宽度、长宽比、面积和周长性状间的相关关系，结果发现，除长宽比以外，其他4个籽粒性状与籽指之间的相关系数均达到极显著水平（r0.01=0.13），其中籽指与籽粒长度、宽度、面积和周长在2015年条件下的相关系数分别为0.79、0.89、0.92和0.85，在2017年条件下的相关系数则分别为0.79、0.87、0.89和0.84。并且，Wang等[30]研究也发现，棉花的籽指与棉仁重量、长度和宽度之间的相关系数亦达到显著或极显著水平，其中籽指与棉仁长度的相关系数为0.21~0.41，籽指与棉仁宽度之间的相关系数为0.64~0.79。由此可见，因棉花籽指与粒长、粒宽、面积和周长之间存在显著相关性，故本研究一方面证实了结果的可靠性，另一方面为揭示性状相关关系的分子机制奠定了理论基础。另外，通过比较本研究与Wang等[30]报道的种仁长度和宽度QTLs发现，在A01、A08、A05和D08染色体存在控制棉籽形状和大小遗传位点，但因Wang等[29]所用遗传图谱属于SLAF-seq SNP标记，未给出标记具体物理位置，因而无法进一步比较本研究关联标记与该文献QTL连锁标记之间的物理距离。

基于上述关联SNPs标记，本研究在其附近筛选到多个与棉籽形状和大小相关的候选基因，在这些基因中，D05染色体候选基因Gh_D05G0144编码种子植物特有转录因子YABBY家族蛋白，该家族属于锌指蛋白超家族，已报道与植物叶片、花和果实等生长发育密切相关[27]。拟南芥基因组中含有6个YABBY基因，其中4个在叶片及其叶源器官（子叶、萼片等）表达，另外2个则只在特异的花器官中表达，如INNER NO OUTER（INO）基因只在胚珠外珠被的最外层细胞层特异表达，并具有促进外珠被生长发育的功能[31]。最近，Sun等[32]研究发现，拟南芥YABBY转录因子INO可通过抑制铁转运蛋白基因（natural resistance associated macrophage protein 1,NRAMP1）的表达，从而减少拟南芥幼嫩种子中的铁元素积累，进而保证植株的生殖生长和种子发育。

同时，候选基因Gh_D05G0148可编码EBF蛋白，该类蛋白通常作为乙烯信号调控因子EIN3的互作蛋白来调控植物生长发育。Guo等[26]研究发现，EBF蛋白可通过介导EIN3/EIL蛋白的降解实现调控乙烯信号通路的功能，番茄中的EBF蛋白基因SlEBF1和SlEBF2通过乙烯信号途径影响植株衰老过程以及果实的生长发育与成熟，并且，番茄EBF基因SlEBF2-like的上游启动子区含有3个乙烯响应元件，在番茄中过表达SlEBF2-like不仅增加了番茄果实的长度，而且改变了果实形状。本研究通过分析EBF蛋白基因Gh_D05G0148在棉花胚珠中的表达量发现，以开花后10 d表达量最高，并在开花后1~35 d始终保持较高水平，暗示其可能参与棉花种子的生长发育。

另外，候选基因Gh_D05G0146编码含TPR结构的家族蛋白，该类蛋白已被发现可参与拟南芥根系生长发育，且与生长素极性运输密切相关，拟南芥TPR家族基因SSR1（short and swollen root 1）突变后，因细胞伸长和增殖严重受损，导致初生根生长被明显抑制，且根系生长素水平显著下降[28]。在本研究中，Gh_D05G0146随棉花胚珠发育表达量升高，因而推测其在棉花种子生长发育过程中发挥一定作用。

综上，本研究发掘的关联SNP标记以及候选基因，为今后开展棉籽大小和形状分子遗传改良提供了依据，也为进一步解析棉花籽粒大小和形状分子机理奠定了基础。