基于BSA-seq的拟南芥辐射诱变突变体的基因定位

汪泽民 何 曦 张宇涵 周 明 洪丽兰，*

（1浙江大学原子核农业科学研究所/农业农村部和浙江省核农学重点实验室，浙江 杭州 310058；2浙江大学生命科学学院植物生物学研究所，浙江 杭州 310058）

从人工诱导或自然发生的突变体中挖掘出突变信息既是分子遗传研究的重要内容，也是培育高产优质新品种的重要途径［1-2］。传统的遗传学基因定位方法往往需要构建渐渗系（recombinant inbred line，RIL）群体和复杂的遗传图谱［3-6］，但这种方式会耗费大量的人力和时间成本。

分离群体分组分析方法（bulked segregant analysis，BSA）常用于鉴定与突变体表型紧密连锁的基因区域或开发连锁的分子标记［7］。在传统BSA方法中，得到突变体后，将突变体与亲本杂交，在后代群体中根据突变表型对个体分组得到两个分离群体，分别构建两个具有表型差异的基因池。在不同基因池中，与突变表型连锁的分子标记会表现出多态性，由此可通过比较分子标记在不同基因池中的差异实现对突变基因的定位［8］。近年来，随着第二代测序技术（next-generation sequencing，NGS）的逐渐成熟，衍生出众多将BSA与NGS相结合的基因组定位策略BSA-seq，如MutMap［9］、QTLseq［10］、QTG-seq［11］和GradedPool-Seq（GPS）［12］等。这些方法大大缩短了基因定位的周期，同时降低了成本。

目前，BSA-seq定位策略已被应用于挖掘黄瓜［13］、水稻［14］、小麦［15］、玉米［16］、花生［17］、大豆［18］等农作物的优良基因。程凤等［19］利用MutMap 定位策略，以黄瓜中短果突变体为试验材料，在黄瓜1号染色体中发现影响果长的基因CsaV3_1G044310。在小麦抗条锈病研究中，研究者们利用BSA-seq 定位到多个相关的数量性状基因座（quantitative trait locus，QTL）位点［20-22］。

辐射诱变具有简单快捷、突变效率高、突变类型多样等优势，常用于正向遗传学研究。在已有的报道中，采用BSA-seq定位策略的研究中突变体材料大多来源于自然突变或化学诱变，对辐射诱变得到的突变体采用BSA-seq 方法定位突变基因的报道较少。as2-7D突变体是ASYMMETRIC LEAVES（AS2）基因的突变体［23］。AS2 基因编码一个含有LATERAL ORGAN BOUNDARIES（LOB）结构域的转录因子，该蛋白质与ASYMMETRIC LEAVES 1（AS1）蛋白互作形成蛋白复合体，共同参与调节拟南芥叶片所有三个轴的正常形态［24］。as2-7D突变体的萼片背面（远离内部花器官的一面）具有凸起，为明显的表型特征，易于进行表型抑制突变体的筛选。鉴于此，本研究以拟南芥突变体as2-7D为材料，对其进行辐射诱变，从后代中筛选萼片背面较为平整光滑的抑制突变体。通过MutMap、QTL-seq 和GPS 三种BSA-seq 方法对其辐射诱变产生的突变体的基因进行定位鉴定，探索BSA-seq 方法用于辐射诱变产生的突变体的基因定位的可能，旨在为加快辐射诱变基因的定位过程提供新思路。

1 材料与方法

1.1 试验材料

拟南芥野生生态型Columbia-0 （Col-0）和Landsbergerecta（Ler）由浙江大学原子核农业科学研究所洪丽兰课题组（本课题组）保存。拟南芥突变体as2-5D和SALK_127850 订购于拟南芥生物资源保藏中心（Arabidopsis Biological Resource Center，ABRC），as2-7D为本课题组前期筛选得到的材料，as2-7D和as2-5D具有相同的点突变，均在AS2 基因启动子上游1 484 bp 处发生了一个G 到A 的点突变，造成AS2 基因的异位表达［23］。as2-7D是Ler 背景，as2-5D是Col-0背景。因为as2-7D比as2-5D具有更显著的萼片背面凸起的表型，故以as2-7D为材料进行抑制突变体筛选。

取as2-7D自交系种子，在浙江省农业科学院作物与核技术利用研究所用60Co-γ 射线诱变，辐射剂量率10 Gy·min-1，辐射剂量600 Gy，M1代植株自交一代，在M2代植株群体中筛选到若干萼片背面凸起表型显著受到抑制的突变体，选取其中一个抑制突变体进行下一步试验。

1.2 表型观察

试验所用拟南芥的生长温度为22 ℃、光周期为16 h 光/8 h 暗，光照强度约为100 µmol·m-2·s-1。植物开花后一周，在DF295 体视镜（Leica，德国）下观察比较拟南芥花萼的表型差异并拍照。

1.3 遗传学分析

将抑制突变体与Ler 和as2-7D分别进行杂交，观察F1代表型并分析突变基因的显隐性；F1代自交得到F2，统计F2代植株群体的表型分离比，通过卡方测验分析抑制突变体萼片凸起受到抑制的表型是否由单基因所调控，从而确定突变基因的遗传模式。

1.4 BSA-seq定位抑制基因

1.4.1 混池测序 取100 株亲本as2-7D的幼嫩叶片构建亲本混池，为亲本池。从抑制突变体和as2-7D杂交得到F2代群体中挑选as2-7D表型和抑制表型的植株各100株，取幼嫩叶片构建子代混池，分别为as2-7D池和抑制突变体池。使用十六烷基三甲基溴化铵（cetyltrimethylammonium bromide，CTAB）提取液分别对三个混池样本提取DNA，送往北京诺禾致源科技股份有限公司测序，DNA 样本质检合格后构建文库，在NovaSeq 6000 测序平台（Illumina，美国）进行双端重测序，测序深度为30×。使用fastp v0.23.1 软件［25］对原始测序结果进行过滤筛选，得到净数据（clean data）。

1.4.2 测序结果分析 对测序结果的分析在服务器上进行，服务器的操作系统为Ubuntu 20.04.4 LTS。从1001Genomes（https：//1001genomes. org/data/MPIPZ/MPIPZJiao2020/releases/current/strains/Ler/）上下载野生型Ler 的参考基因组序列以及基因注释信息文件。使用fastqc v0.11.8 软件查看测序质量，质量合格后分别使用MutMap、QTL-seq和GPS分析测序结果。

使用MutMap 方法分析测序结果：从github 上下载由Sugihara等［26］开发的MutMap pipeline（https：//github.com/YuSugihara/MutMap），使用该pipeline 默认参数和实际混池样本数，分析亲本池和抑制突变体池的clean data，得到抑制基因的候选区间。

使用QTL-seq 方法分析测序结果：从github 上下载由Sugihara 等［26］开发的QTL-seq pipeline（https：//github.com/YuSugihara/QTL-seq），使用该pipeline 默认参数和实际混池样本数，分析亲本池、as2-7D池和抑制突变体池的clean data，得到抑制基因的候选区间。

使用GPS方法分析测序结果：利用GPS方法对as2-7D池和抑制突变体池的clean data进行分析。首先利用bwa v0.7.17软件［27］将clean data与Ler的参考基因组比对定位。使用Picard软件（https：//broadinstitute.github.io/picard）标记重复序列。利用GATK v4.2.4.0软件［28］实现对单核苷酸多态性（single nucleotide polymorphisms，SNP）和插入缺失（insertion-deletion mutations，InDel）变异位点的检测。过滤掉低质量和低测序深度的SNP和InDel 位点，使用Ridit（relative to an identified distribution unit）分析法计算-lg（P）值，得到与参考基因组序列具有显著差异的突变位点及抑制基因的候选区间。

1.4.3 筛选可信的突变位点 在BSA-seq 的分析流程中得到基因序列变异信息（variant call format，vcf）文件，使用SnpEff v4.3t 软件［29］对所有突变位点进行注释，结合三种BAS-seq 分析策略和Integrative Genomics Viewer （IGV）软件［30］的基因组可视化结果，实现对抑制基因的精细定位。

2 结果与分析

2.1 as2-7D抑制突变体的获得及表型特征

在as2-7D的辐射诱变群体中筛选到一株萼片背面凸起表型显著受到抑制的突变体。野生型Ler 萼片背面光滑平整（图1-A），as2-7D萼片背面具有大量的凸起和褶皱（图1-B），而抑制突变体的萼片背面较为平整，没有凸起和褶皱，更接近于Ler野生型（图1-C）。对该突变体的AS2基因测序，发现其仍然含有as2-7D的突变位点，且该基因区域不具有第二个突变位点，说明该突变体是由AS2基因之外的位点突变导致的对as2-7D萼片表型的抑制，选择该抑制突变体进一步研究。

图1 as2-7D ssp1抑制突变体的表型Fig.1 The phenotypes of as2-7D ssp1 mutant

2.2 as2-7D抑制突变体的遗传分析

将抑制突变体分别与as2-7D和Ler 杂交，观察F1代的表型，抑制突变体 ×as2-7DF1代的萼片表型与as2-7D萼片表型一致（图1-B、F），说明抑制突变体对as2-7D萼片表型的抑制是由隐性突变控制的。统计F1代自交得到的F2代植株群体的表型分离比，群体共43 个植株，其中类as2-7D32 株，类抑制突变体11 株，经卡方测验分析得P值约为0.93，大于0.05，符合3∶1的分离比（表1），说明抑制突变体抑制as2-7D萼片表型的性状是由一个单基因控制的，将该抑制基因命名为SUPPRESSORSofas2-7D SEPAL PHENOTYPE 1（SSP1），由此抑制突变体命名为as2-7D ssp1。

表1 as2-7D ssp1 × as2-7D F2代表型分离统计表Table 1 Phenotypic segregation of the as2-7D ssp1 × as2-7D F2 population

as2-7D ssp1× Ler F1代植株萼片背面仍然出现凸起，但凸起的数量较as2-7D有所减少（图1-E），与as2-7D× Ler F1代具有相似的表型（图1-D）。由于该表型凸起的数量介于Ler 和as2-7D之间，将该表型称为中间型。as2-7D ssp1× Ler F2代群体中，类as2-7D∶中间型∶类Ler∶类as2-7D ssp1∶新表型的个体数量比符合3∶6∶3∶3∶1 的分离比（表2）。综合上述试验数据推测，ssp1是一个单基因隐性突变，且不与AS2基因连锁；as2-7D ssp1× Ler F2代群体中的新表型可能由SSP1单基因突变导致。

表2 as2-7D ssp1 × Ler F2代表型分离统计表Table 2 Phenotypic segregation of the as2-7D ssp1 × Ler F2 population

2.3 ssp1突变基因定位

对as2-7D亲本池及as2-7D ssp1×as2-7DF2代群体的两个表型池进行全基因组测序，分别用MutMap、QTL-seq、GPS 方法分析测序结果，限定SSP1基因所在的位置，同时用基因组可视化软件IGV 筛选候选基因。

图2、3 分别为MutMap 和QTL-seq 的分析结果，蓝色散点表示某些突变位点的SNP-index 值，绿色折线表示95%置信区间的阈值，黄色折线表示99%置信区间的阈值，红色折线为2 Mb 内的平均SNP-index。MutMap 法引入了SNP-index，表示混池中非参考基因组类型的SNP频率。抑制突变体池中的某SNP标记与抑制表型的关联度越高，该标记的SNP-index 越接近于1。查看MutMap 分析结果图，在2 号染色体大概9 Mb 以及13～16 Mb 的位置范围内，平均SNP-index 已超过99%的置信阈值，但是在9～11 Mb 区间内无可信的SNP 位点，未能绘制折线，使得后续2 Mb 范围内的平均SNP-index 值偏低。因此，根据MutMap 的结果，将SSP1基因限定在2号染色体9～16 Mb区间内（图2）。理论上，as2-7D ssp1×as2-7DF2代群体中两个极端表型混池在SSP1基因所在区域的SNP-index 具有显著差异，而其他区域差异较小。因此，可使用QTL-seq定位SSP1基因。QTL-seq 分别计算子代两个极端表型混池的SNP-index，将两个混池的SNP-index 做减法得到∆SNP-index。∆SNP-index越高的基因区域，与SSP1基因的关联度也越高。从QTL-seq 分析图同样得出SSP1基因位于2号染色体9～16 Mb区间内的结论（图3）。GPS 分析结果显示，2号染色体9～15 Mb区间内富集了高-lg（P）值的散点（图4），与MutMap 和QTL-seq 的分析结果基本一致。

图2 MutMap分析定位SSP1基因的结果Fig.2 Results of MutMap analysis to locate the SSP1 gene

图3 QTL-seq分析定位SSP1基因的结果Fig.3 Results of QTL-seq analysis to locate the SSP1 gene

图4 GPS分析定位SSP1基因的结果Fig.4 Results of GPS analysis to locate the SSP1 gene

在2 号染色体的9～16 Mb 区间内，共有391 个SNP和Indel 突变位点，根据质量进行初步筛选，得到139个可信的突变位点。以SNP-index 和ΔSNP-index 为指标进一步过滤候选位点，共得到2 个SNP-index 或ΔSNP-index高于99%置信度的候选位点。snpEff软件的注释信息表明这两个突变位点位于基因的上游或下游，影响其附近基因功能的概率较低，排除这两个突变位点为候选突变位点。

基因可视化软件IGV 发现，抑制突变体混池在2 号染色体11 950 417～11 975 606 bp 范围内发生了一个25 189 bp的缺失（图5）。为验证在as2-7D ssp1群体中存在该片段缺失，分别在缺失片段两端设计鉴定引物，对6 株as2-7D ssp1植株的基因组DNA 进行扩增，结果显示，as2-7D ssp1确实存在该大片段缺失。在该缺失片段中共包含AT2G28580、AT2G28590、AT2G28600、AT2G28605、AT2G28610、AT2G28620 这6 个基因。其中，AT2G28610编码一个转录因子PRESSED FLOWER（PRS），在拟南芥侧生器官（如萼片、叶等）边缘表达，具有调控细胞增殖的功能。根据前人的研究报道，PRS过表达会导致拟南芥萼片背面出现凸起的表型［31］，与as2-7D类似，因此将PRS作为SSP1基因的候选基因。

图5 as2-7D ssp1突变体基因组中的大片段缺失Fig.5 Large deletion in the genome of the as2-7D ssp1

2.4 PRS作为抑制基因的验证

为验证PRS基因是否是SSP1基因，从Col-0 背景的PRS基因T-DNA 插入突变体SALK_127850，分离出纯合的PRS基因突变体prs-3，构建了as2-5D与prs-3的双突变体as2-5D prs-3。观察Col-0、as2-5D、prs-3和as2-5D prs-3 的萼片表型发现，as2-5D的萼片背面具有少量的凸起（图6-B），prs-3 萼片背面平整光滑（图6-C），与Col-0（图6-A）相比无明显差异，而as2-5D prs-3 的萼片背面没有凸起，整体较为平整（图6-D），说明PRS基因的缺失可以抑制as2-5D萼片背面凸起的表型，验证了PRS基因是SSP1基因。

图6 PRS突变抑制as2-5D突变体的萼片表型Fig.6 PRS mutation suppresses the sepal phenotype of as2-5D mutant

3 讨论

本研究使用60Co-γ 射线辐射诱变as2-7D，在辐射诱变后代中筛选得到抑制突变体as2-7Dssp1，成功地运用BSA-seq 定位策略对抑制基因SSP1进行了定位。BSA-seq 虽然无法精确定位大片段的插入或缺失突变，但可以在全基因组范围内筛选出于目标基因高度连锁的SNP 位点，从而将目标基因的位置限制在较小的范围内。在得到目标基因的候选区间后，可通过基因组注释文件对该范围内的可信SNP 位点进行筛选，并结合IGV软件查看该区间内是否存在大片段的插入或缺失，从而实现对目标基因的精细定位。为使用BSA-seq 定位策略挖掘由辐射导致的突变基因提供了参考依据。

本研究通过MutMap、QTL-seq 和GPS 三种BSAseq 分析方法对辐射诱变突变体的基因进行定位。MutMap 基于子代突变体池的SNP-index 对目标性状的基因区域进行定位。在子代突变体群体的遗传背景与亲本一致或高度相似时，MutMap具有良好的定位效果。本研究中的抑制突变体由as2-7D辐射诱变得到，二者的杂交F2代群体具有相似的遗传背景，因此可用MutMap定位抑制基因。QTL-seq的技术思路与MutMap相类似。为减小背景噪音的影响，对两个子代混池分别求SNP-index，相减后得∆SNP-index，根据∆SNPindex 对抑制基因进行定位。本研究中MutMap 和QTL-seq 均将抑制基因定位在2 号染色体的9～16 Mb区间内。GPS 方法没有使用SNP-index 算法，而是使用Ridit 分析法来筛选与目标基因关联度较高的SNP位点，通过窗口smooth 来减小背景噪音的影响。本研究中，三种BSA-seq 分析方法均成功定位到了抑制突变体的变异区间，且取得了相似的定位效果。因此，MutMap、QTL-seq 和GPS 均可应用于辐射诱变突变体的基因定位。

4 结论

本研究在as2-7D的辐射诱变突变体后代中筛选到抑制突变体as2-7D ssp1，使用MutMap、QTL-seq 以及GPS 三种BSA-seq 定位策略对引起突变性状的SSP1基因进行定位，三种方法均成功定位到了抑制突变体的变异区间；并通过基因组注释信息和遗传学试验分析，确认变异区间内的PRS基因的缺失导致了as2-7D ssp1的突变表型。结果表明BSA-seq定位策略可用于定位辐射诱变引起的大片段缺失突变。