基于InDel标记的谷子株高QTL定位

杜晓芬 钱枰励 唐楚楚 杜德杰 韩康妮 李禹欣 王智兰 王 军，*

（1山西农业大学谷子研究所/山西省后稷实验室，山西 长治 046011；2山西农业大学农学院，山西 太谷 030801；3中国农业大学农学院，北京 100193）

谷子（Setariaitalica）具有耗水少、抗旱耐瘠、营养丰富均衡、粮饲兼用等特点，是旱作生态农业绿色发展的主栽作物［1］。随着谷子参考基因组的公布［2-3］，谷子基础研究在近十年飞速发展，并逐渐发展为禾本科作物功能基因组研究的模式作物［4］。

分子标记是作物遗传研究的基础，其中，插入/缺失（insertion/deletion，InDel）多态性标记是指等位基因的核苷酸序列在个体间因插入或缺失引起的长度变异，具有分布广、密度高、变异稳定、重复性好、性价比高等优点［5］，广泛应用于种质资源分析、遗传多样性分析、亲缘关系鉴定、种质鉴定、遗传图谱构建、重要性状QTL 定位与图位克隆等方面。赵庆英等［6］对名优谷子品种晋谷21 号进行重测序，与参考基因组相比发现了169 037 个InDel 位点，进一步验证了68 个InDel 位点，并开发了一个晋谷21 号的特异InDel 标记，为分子标记辅助育种提供了分子标记资源。Zhang 等［7］基于谷子品种SSR41 的重测序结果，开发了156 个适合琼脂糖凝胶电泳的InDel 标记，并利用多态性InDel 标记对谷子白叶鞘基因进行了图位克隆。王智兰等［8］根据SiGRAS23开发了一个与株高性状紧密连锁的InDel 标记D8-1，可用于谷子种质资源株高变异的筛选。张硕等［9］利用7 个InDel 标记将谷子白条纹叶基因定位于20.13 Mb 内，为该基因的精细定位和克隆奠定了基础。由此可见，InDel标记在谷子基础研究方面发挥了重要的作用。

株高是谷子的重要农艺性状之一，主要通过调节株型和抗倒伏性影响谷子产量。前人遗传分析结果表明，谷子株高是典型的数量性状，由多基因控制［10］。近十年来，研究人员定位了许多株高相关的数量性状位点（quantitative trait loci，QTL）［11-16］，这些QTL 在谷子的9条染色体上均有分布。同时，前人精细定位和克隆了一些候选基因，例如编码DELLA蛋白的SiD1和编码细胞色素P450 蛋白的SiD2［17-18］。尽管如此，相比水稻、小麦等主要农作物，谷子株高QTL 或基因的鉴定仍相对较少，这严重制约了谷子株高遗传机理研究的进展。

鉴于此，本研究对2 个株高存在极显著差异的谷子品种（衡谷12号和长农35号）进行深度重测序，比较两者之间的InDel差异，开发InDel标记，并利用多态性InDel标记扫描“衡谷12号×长农35号” F2分离群体，构建遗传连锁图谱、定位株高QTL，进一步在重组自交系群体中验证株高QTL，并对候选基因进行预测。本研究旨在为谷子基础研究提供一套好用的InDel标记，同时为株高基因的克隆奠定基础。

1 材料与方法

1.1 试验材料及性状调查

本研究中，用于株高QTL 定位的杂交组合的母本为衡谷12 号，由河北农林科学院旱作农业研究所选育，具有矮秆、极早熟、分蘖性强等特点［19］，父本长农35 号由山西农业大学谷子研究所选育，属于高秆、晚熟、无分蘖品种［20］。2016年，将衡谷12号和长农35号进行杂交；2017年，鉴定F1植株；2018年，将F1自交获得F2分离群体，在谷子成熟期调查株高，其中，双亲调查5株，取平均值作为性状值，F2子代逐株挂牌调查株高。

采用单粒传法，将F2分离群体繁殖到F6代。2022年，将双亲和重组自交系群体（recombinant inbred line，RIL，F6代共283个株系）种植2行，行长2 m，株距0.33 m。在谷子成熟期调查株高，每个株系调查5株，以平均值作为性状值。所有材料均种植于山西农业大学谷子研究所试验田，田间管理、水肥等条件一致。

1.2 DNA提取

在谷子拔节期，分别采集新出叶片，液氮速冻后置于-80 ℃备用。采用改良的十二烷基硫酸钠（sodium dodecyl sulfate，SDS）法提取总DNA。SDS提取液配制：1 mol·L-1Tris-HCl（pH 值8.5）100 mL，0.5 mol·L-1乙二胺四乙酸（ethylene diamine tetraacetic acid，EDTA）（pH 值8.0）100 mL，5 mol·L-1NaCl 20 mL，20% SDS 100 mL，加去离子水定容至1 L。提取步骤为：将叶片放入研钵中加入液氮研磨成粉状，转移至2.0 mL 离心管中；加入SDS 提取液750 µL，混匀后置于65 ℃水浴锅中1 h；加入750 µL苯酚/氯仿（1∶1），混匀15 min后，4 ℃，12 000 r·min-1离心15 min；吸取上清液600 µL，加等体积氯仿，混匀15 min 后，4 ℃，12 000 r·min-1离心15 min；吸取上清液450 µL，加0.6 倍体积的异丙醇混匀，静置30 min；4 ℃，12 000 r·min-1离心10 min 后弃上清，用70%乙醇清洗2 次；室温干燥后加入100 µL 的Tris-EDTA（TE）缓冲液备用。

1.3 RNA提取和cDNA合成

分别在谷子孕穗期和开花期，对衡谷12 号和长农35 号的倒二茎进行取样，液氮速冻研磨，参照RNAiso Plus 试剂（9108，北京宝日医生物技术有限公司）说明书提取总RNA，参照PrimeScript Ⅱ 1st strand cDNA synthesis kit 试剂盒（6210A，北京宝日医生物技术有限公司）说明书合成cDNA。

1.4 重测序

在苗期，分别随机采集衡谷12号和长农35号的叶片，采用高效植物基因组DNA 提取试剂盒（DP350，北京天根生化科技有限公司）提取DNA，送至北京百迈客生物科技有限公司进行重测序，重测序按照Illumina公司（美国）提供的标准步骤执行，简要流程为：对提取的基因组DNA 进行质检，DNA 检测合格后用超声波将DNA片段化；对片段化的DNA进行末端修复、3'端加A和连接测序接头；用琼脂糖凝胶电泳对片段大小进行选择，并进行PCR 扩增形成测序文库；对测序文库进行文库质检，质检合格的文库用Hiseq 2500 测序仪（Illumina公司，美国）进行测序。

1.5 基因组变异检测

测序得到的原始序列（raw reads）进行过滤，获得高质量的序列（clean reads），过滤的标准为：去除带接头的reads；过滤N 含量超过10%的reads；去除质量值低于10、碱基超过50%的reads。利用Burrows-Wheeler Transform 软件［21］将获得的clean reads 比对到谷子参考基因组豫谷1号上（Setariaitalicav2.2，https：//phytozome.jgi.doe.gov/pz/portal.html#%21info？alias=Org_Sitalica_er）。使用Picard（https：//sourceforge.net/projects/picard/）去掉重复，Genome Analysis ToolKit 软件［22］用于检测单核苷酸多态性（single nucleotide polymorphisms，SNP）和InDel，并用SnpEff 软件［23］对检测到的SNP 和InDel 进行注释，Breakdancer（http：//breakdancer.sourceforge.net/）用于检测结构变异（structral variation，SV），对于检测到的SV，通过比对变异位点在参考基因组上的位置对其进行注释。

1.6 InDel标记开发和多态性鉴定

比较衡谷12号和长农35号差异InDel位点，选择衡谷12号和长农35号插入或缺失碱基数量差异大于等于3 bp且间隔大于5 kb的位点作为引物设计的目标位点。分别提取目标位点上下游侧翼序列150 bp，整理成Fasta格式，利用BatchPrimer3 v2.0（https：//probes.pw.usda.gov/batchprimer3/）批量设计引物，产物大小设置为100～300 bp，引物由南京金斯瑞生物科技股份有限公司合成。

PCR 反应体系为10 µL：10×PCR buffer 1.0 µL，2 µmol·L-1上下游引物各1 µL，rTap DNA聚合酶0.1 µL，2.5 mmol·L-1dNTP 0.8 µL，模板DNA 1.0 µL，灭菌水5.1 µL。PCR 扩增程序为：94 ℃预变性5 min，94 ℃变性30 s，58 ℃退火30 s，72 ℃延伸30 s，35个循环；72 ℃终延伸10 min。PCR 扩增产物利用8%聚丙烯酰胺凝胶电泳进行分离，经硝酸银染色后，拍照记录。

1.7 遗传图谱构建和QTL定位

选择能够良好分型的多态性InDel标记，分别扫描F2群体（120 个单株）和RIL 群体（283 个株系），对扩增结果进行统计，其中与衡谷12号带型一致的记为“A”，与长农35 号带型一致的记为“B”，杂合带型记为“H”，条带缺失时记为“-”。采用JoinMap 4［24］构建遗传连锁图谱，选择回归算法（regression mapping）和Kosamb，s函数，其他参数为默认参数。采用WinQTLCart 2.5［25］对株高进行QTL 定位，使用复合区间作图（composite interval mapping，CIM）法，在P＜0.05 条件下，进行1 000 次排列测验（permutation test）确定QTL 的似然比的自然对数（likelihood of odd，LOD）阈值。QTL命名规则按照q+性状简称+染色体进行命名，当同一条染色体出现多个QTL，按照-1、-2、-3的顺序命名。

1.8 候选基因预测和分析

利用Phytozome数据库（https：//phytozome-next.jgi.doe.gov/），以豫谷1 号为参考基因组，对定位区间的基因及其功能注释进行分析。同时根据衡谷12 号和长农35 号重测序数据，分析定位区间内编码区发生突变的基因，预测候选基因。对预测的候选基因进行克隆测序，采用DNAMAN 7 软件对测序结果进行比对分析。参考王智兰等［8］的方法进行候选基因的荧光实时定量表达分析。

1.9 数据统计分析

利用SPSS 17 分析群体性状的最大值、最小值、平均数、标准差、峰度、偏度以及进行独立样本t检验。利用GraphPad Prism 5进行频率分布分析和作图。

2 结果与分析

2.1 衡谷12号和长农35号重测序及变异分析

衡谷12 号和长农35 号的clean data 为31.82 Gbp，Q30 达到88.38%，两个样品与参考基因组平均比对率为98.63%，平均覆盖深度为34×，其中衡谷12 号的覆盖深度为30×，长农35 号的覆盖深度为38×。变异检测结果表明，衡谷12 号中单核苷酸多态性（SNP）、插入/缺失（InDel）和结构变异（SV）数量分别为524 486、103 442 和22 295 个，长农35 号中SNP、InDel 和SV 数量分别为1 131 571、196 477 和24 509 个（图1-A），比较而言，检测到的SNP 数量最多，SV 数量最少，且长农35 号中这三种变异类型都高于衡谷12 号，同时发现，长农35 号和衡谷12 号在相同位点共同检测到SNP 变异位点和InDel 变异位点分别有259 259 和55 224 个（图1-E、F）。

SNP 变异类型分为碱基转换（T/A 和C/G）和碱基颠换（A/C、G/T、C/T 和G/A），衡谷12 号中碱基颠换的数量占总变异类型的87.7%，长农35 号中碱基颠换的数量占总变异类型的88%，同时发现，长农35号中的碱基转换和碱基颠换的数量均高于衡谷12号（图1-B）。

SV 变异类型分为大片段插入（insertion，INS）、大片段缺失（deletion，DEL）、反转（inversion，INV）、染色体内易位（intra-chromosomal translocation，ITX）和染色体间易位（inter-chromosomal translocation，CTX）等，衡谷12 号和长农35 号的变异类型主要是INS 和DEL，分别占总变异类型的84.76%和79.17%（图1-C）。

InDel变异类型分为小片段插入（INS）和小片段缺失（DEL），分析表明，在衡谷12号和长农35号中，主要发生均为1 bp的INS或DEL，同时发现大于等于39 bp且小于等于50 bp的INS数量高于DEL数量（图1-D）。

2.2 InDel标记的开发和多态性分析

比对衡谷12 号和长农35 号中InDel 变异位点，选择两者差异大于等于3 bp 且间隔大于5 kb 的位点，最后筛选了3 961 个位点，并设计引物。首先，筛选在双亲间具有多态性的标记，初步获得1 392 对多态性标记，多态性率为35.14%，其中第3 染色体的多态性标记最多，第4 染色体的多态性标记最少（图2-A）。然后，用F2群体的部分单株进一步筛选，获得591 个稳定性好，且在后代群体中能够良好分型的多态性标记（图2-B），为下一步遗传图谱构建奠定了基础。

图2 InDel标记的多态性分析Fig.2 Development and polymorphism analysis of InDel markers

2.3 遗传图谱构建及F2群体株高QTL定位

利用591个InDel标记对包含120个单株的F2群体进行基因分型，采用JoinMap 4 进行连锁分析，最后获得一张包括467 个InDel 标记的谷子遗传连锁图谱（电子附图1），该图谱总图距448.45 cM，平均图距0.96 cM，其中第4 和第8 染色体标记数量最少（16 个），第3 染色体标记数量最多（160 个），第1 染色体图距最短（9.67 cM），第9染色体图距最长（105.22 cM）（表1）。

表1 “衡谷12号×长农35号” F2群体遗传连锁图谱信息Table 1 Characteristics of the genetic map in F2 population of ‘Henggu12 × Changnong35’

对双亲及120 个F2子代株高进行分析，发现双亲之间差异极显著（P＜0.01，图3-A），株高在F2群体中呈双向超亲分离（表2），频率分布结果显示，株高呈双峰分布（图3-A）。

表2 双亲及群体株高统计分析Table 2 Phenotypic analysis of plant height in the parents，F2 population and RIL population /cm

图3 “衡谷12号×长农35号” F2群体和RIL群体株高频率分布直方图Fig.3 Histogram of frequency distribution of plant height in F2 population and RIL population of ‘Henggu 12×Changnong 35’

利用WinQTLCart 2.5 对F2群体的株高进行QTL定位，共检测到4个QTL 位点，除qPH5-1外，加性效应均为负值，说明来自父本长农35 号的等位位点对株高贡献率大。在4 个QTL 中，第5 染色体检测到2 个QTL，分别为qPH5-1和qPH5-2，其中qPH5-1可解释表型变异率最大，为8.41%；第9 染色体检测到2 个QTL，分别为qPH9-1和qPH9-2，其中qPH9-2可解释表型变异率为16.47%，为主效QTL（表3）。

表3 “衡谷12号×长农35号” F2群体株高QTL定位Table 3 Putative QTL detected for plant height in F2 population of ‘Henggu12×Changnong35’

2.4 RIL群体株高QTL的验证

首先，对双亲和283 个株系的株高进行分析，发现株高在RIL 群体中呈单向超亲分离，变异范围比F2群体小，但是变异系数比F2群体大（表2）。频率分布结果显示，与F2群体的结果相似，株高呈双峰分布（图3-B）。然后在RIL 群体中对F2群体中鉴定的效应较大的QTL（qPH5-1和qPH9-2）进行了验证。结果表明，在RIL 群体中，仅检测到qPH9-2，LOD 值为93.6，加性效应为-46.15，解释表型贡献率达91.06%（表4）。

表4 “衡谷12号×长农35号”重组自交系群体株高QTL定位Table 4 Putative QTL detected for plant height in RIL population of ‘Henggu12 × Changnong35’

2.5 qPH9-2候选基因分析

利用Phytozome 数据库分析qPH9-2定位区间（2 776 961～7 257 031 bp），该区间包含706 个基因，分析衡谷12号和长农35号在该区间的重测序结果，发现70个基因在编码区、5'UTR和3'UTR发生突变。其中，27 个基因由于SNP 变异发生非同义突变，4 个基因由于InDel 变异发生移码突变（电子附表1）。分析发生非同义突变的基因功能注释，发现编码bHLH 转录因子的基因Seita.9G080400，在其5'UTR 处和编码区分别存在2 处SNP。由于转录因子在植物生长发育中发挥重要作用，因此对其进行了重点分析。

根据参考基因组Seita.9G080400的基因组序列，设计特异引物（电子附表2），分别对衡谷12 号和长农35 号进行克隆测序，测序结果验证了其5'UTR 的2 处SNP（SNP1和SNP2）和第一外显子的2处SNP（SNP3和SNP4），同时还发现在其第一内含子、第二内含子和第三内含子分别存在1 处SNP（SNP5、SNP6 和SNP7）（电子附图2）。外显子中的SNP3在衡谷12中编码缬氨酸（V），在长农35号中编码丙氨酸（A）；外显子中的SNP4在衡谷12 中编码组氨酸（H），在长农35 号中编码天冬氨酸（D）（图4-A）。表达分析结果表明，孕穗期茎中，Seita.9G080400相对表达量在衡谷12号中是长农35号的10倍；开花期茎中，Seita.9G080400相对表达量在衡谷12号中是长农35号的42倍（图4-B）。

图4 Seita.9G080400序列分析和表达分析Fig.4 Sequences and expression level analysis of Seita.9G080400 between Henggu12 and Changnong35

3 讨论

随着高通量测序技术的快速发展，基于重测序技术，在全基因组水平发掘的分子标记在作物中得到广泛应用。大量研究表明，InDel标记是植物中的第二大类多态性高的标记，具有准确性高、共显性、易于基因分型等诸多优点［5］，已在遗传图谱构建、重要性状QTL定位、种质资源分析与鉴定等方面广泛应用。谷子中相关文献证实，InDel标记数量仅次于SNP标记［6-7，26-28］。例如，Li等［27］对312份谷子种质进行重测序，共鉴定了3 027 706 个SNP 和256 826 个InDel，Bai 等［28］在农家种十里香中检测到915 434 个SNP 和28 546 个InDel。本研究中，衡谷12 号中的SNP 和InDel 分别为524 486和103 442个，长农35号中的SNP和InDel分别为1 131 571和196 477 个。由于测序深度和材料的不同，导致检测到的SNP 和InDel 数量有所不同。本研究进一步验证了InDel 标记是一类具有丰富多态性的标记，此外，本研究鉴定的大量InDel位点为开发标记提供了分子基础。

目前已有研究中，主要利用简单重复序列（simple sequence repeat，SSR）标记［29-31］进行种质资源群体结构分析［32］、图位克隆［33］、优异等位变异检测［34］，而InDel标记应用的相关报道较少。尽管谷子中InDel变异丰富，但对其进行开发并验证的标记数量同样相对较少，谷子中可用的InDel 标记资源有限，因此在全基因组水平开发InDel 标记非常必要。赵庆英等［6］在晋谷21 中开发验证了68 个InDel 标记，Zhang 等［7］在SSR41中开发验证了156个InDel标记。本研究选择了3 961 个在衡谷12 号和长农35 号之间有差异的InDel位点，经过初步筛选，获得了1 392 个多态性InDel 标记，多态性率为35.14%。与前人研究结果相比，本研究验证的InDel标记数量较多，并能够有效应用于株高QTL定位。

数量性状QTL 定位结果易受定位群体和环境条件影响，尽管初级定位群体（F2、F2：3和DH）表型数据不能重复，但是一些大效应的QTL 还是会被准确定位。例如，水稻“绿色”基因DEP1的克隆，最初利用DH 群体检测到控制水稻直立穗型的主效QTL（qPE9-1，贡献率为41.72%）［35］，随后利用2 个F2群体及1 521 个F2子代将其精细定位并确定候选基因DEP1［36］。本研究利用开发的多态性InDel 标记，构建了“衡谷12 号×长农35 号” F2群体的遗传连锁图谱，初步定位了4 个株高QTL，并对其中2 个效应值较大的QTL 在RIL 群体中进行了验证，相比F2群体的定位结果，qPH9-2在RIL 群体中能重复检测到，贡献率达91.06%，且物理区间（2 776 961～7 257 031 bp）更小。前人在谷子第9染色体上相同区间也检测到株高QTL（H9a：3 452 395～9 964 754 bp），贡献率为7.6%～15.5%［12］，暗示该区间存在控制株高的主效QTL。

bHLH 是植物中重要的一类转录因子，参与植物生长发育与形态建成。研究表明，小麦TabHLH123调控小麦根系发育，TabHLH123-6B与根重、千粒重和株高相关［37］，TabHLH123-6A与根系构型相关［38］。本研究中，Seita.9G080400编码bHLH 转录因子，测序结果表明，在衡谷12 号和长农35 号中，其第一外显子有2处SNP，导致氨基酸发生非同义突变。其中，SNP3在衡谷12 号中编码缬氨酸，长农35 号中编码丙氨酸，这两种氨基酸都为疏水性氨基酸；SNP4 在衡谷12 号中编码组氨酸，长农35 号中编码天冬氨酸，由碱性氨基酸变为酸性氨基酸，这些氨基酸的改变是否影响蛋白的功能，需要进一步验证。另外，在不同发育时期的茎中，Seita.9G080400相对表达量在衡谷12 号和长农35号之间存在极显著差异。综合以上结果，推测Seita.9G080400可能是控制株高的关键候选基因，但是仍然需要通过进一步精细定位和功能鉴定进行验证。

4 结论

本研究利用2 个谷子品种的重测序数据，在全基因组水平开发了1 392 个分布于谷子9 条染色体的InDel 标记。利用多态性InDel 标记，构建了一张包含467 个InDel 标记的谷子遗传连锁图谱，定位了1 个控制株高的主效QTL，预测了1 个关键候选基因，验证了InDel标记的有效性和应用性，为今后谷子种质资源鉴定、群体结构分析和亲缘关系分析等提供了良好的分子标记资源。

电子附表1 突变基因的信息Electronic Table S1 Information of mutant genes

电子附图1 “衡谷12号×长农35号”F2群体InDel标记连锁遗传图Electronic Fig.S1 The genetic linkage map constructed using 467 InDel markers in F2 population derived from ‘Henggu12 ×Changnong35’