γ射线辐射诱导水稻突变体的基因组变异分析

张 莉，张 勇，张渊海，张逸妍，解 涛，彭选明，杨 震*

（1. 湖南省核农学与航天育种研究所，长沙 410125；2. 湖南省农业生物辐照工程技术研究中心，长沙 410125）

水稻作为重要的粮食作物，它的丰产、稳产对国民经济发展具有重要的战略意义。诱变育种技术在农作物新材料创制和优良新品种的培育中已经发挥了显著作用［1-3］。常规水稻品种的诱变策略通常是改进一个或者一些容易鉴别的性状，并利用突变体植株扫描突变位点和定位相关基因来研究突变的分子机制［4-5］。植物辐射诱变通过改变染色体数目、染色体结构或DNA的核苷酸构成从而引起表型变异［6］。在重离子辐射诱变DRF缺失的烟草突变体中，NtDRF2整个基因缺失，NtDRF1发生1个碱基缺失使肽链编码发生移框突变，最终形成白花烟草［7］。Wang等［8］通过12C6+辐射月牙藻获得了5个叶绿素a缺陷的突变体，研究发现突变体中捕光色素复合体相关的Lhcb5、Lhcbm5和Lhcbm1基因的表达发生显著变化。目前基因组水平上诱导突变的分子机制仍不清楚，随着生物信息技术的高速发展，通过全基因组测序（whole-genome sequencing，WGS）技术可以大规模检测许多植物的个体突变［9］，WGS的高效率和低成本使其成为发掘基因组突变的一种好方法。

遗传变异包括序列变异和结构变异。序列变异包括单核苷酸多态性（single nucleotide polymorphisms，SNPs）、插入缺失（insertions and deletions，In-Dels）、微卫星或简单序列重复（microsatellites or simple sequence repeats，SSRs）和转座因子。这些序列多态性由于其低成本、稳定性和高通量的应用已被广泛应用于基因组选择、数量性状位点（quantitative trait locus，QTL）定位、单倍型和家系分析等领域［10］。如Zhang等［11］通过鉴定早熟三叶柑桔与其野生型的全基因组遗传变异，开发出新的遗传标记进行柑桔重要性状研究。结构变异（structural variations，SVs）包括大片断插入（large-scale insertions，INSs）、缺失（deletions，DELs）和倒置（inversions，INVs）以及基因组水平的染色体内部和染色体间易位，这些类型的结构变异统称为拷贝数变异（copy number variants，CNVs），这些变异在形成基因组多样性方面起着重要作用。Muñoz-Amatriaín等［12］采用8个大麦栽培品种和6个野生大麦品种进行基因组杂交比较。与栽培大麦相比，野生大麦中存在更高水平的拷贝数变异多样性。基因拷贝数占基因阵列编码序列的9.5%，被拷贝数影响的基因标记为抗病蛋白和蛋白激酶。栽培大麦Barke和Morex品种的CNVs序列比较表明，单链退火和合成依赖链退火的双链断裂DNA修复机制在大麦CNVs的发生过程中起着重要作用。Wallace等［13］对玉米5 000个近交系中的41种不同表型进行全基因组关联分析，检测到28 900 000个SNPs和800 000个CNVs，发现基因区和基因间区具有相反的富集模式、较小的等位基因频率和效应大小，同时全基因组关联分析（genomewide association study，GWAS）标记的基因具有丰富的调控功能，这表明基因调控和基因复制是表型变异的强大驱动因素。

目前对作物序列多态性的研究越来越多，不同作物中已经开发了几种模型的全基因组SNPs和InDels数据库［14-15］，并广泛应用于遗传多样性研究［16］、功能和进化研究［17-18］以及驯化和基因组进化研究［9，19-20］。本研究利用60Co-γ辐照水稻品种黄华占获得遗传稳定突变体湘辐1821，该突变体较黄华占株高增高，叶宽增宽，单产显著提高。同时，利用高通量深度测序来研究其基因组变异，以期为辐射育种技术提供理论支撑。

1 材料与方法

1.1 试验材料

水稻突变体湘辐1821（Xiangfu 1821，MT_HHZ，Oryza satiuaL.）系350 Gy60Co-γ射线辐照亲本黄华占（Huang Hua Zhan，WT_HHZ，Oryza satiuaL.）干种子，经多代大田选育而成。水稻品种黄华占购买于湖南农丰种业有限公司。

1.2 诱变处理

2016年在湖南省核农学与航天育种研究所辐照中心利用60Co-γ射线对黄华占进行诱变处理，辐照剂量为350 Gy，剂量率为5.25 Gy/min。

1.3 种植方法

M1代编号B16314当季大田多本粗插种植，混收，只收获其主穗种子，以期获得更大的变异率。2017年随机选取B16314种子，单本种植，群体量20 000株左右，在其田间选择出综合性状良好的突变体，单株收获，编号B1716。2017年冬在海南三亚繁殖B1716种子，并观察其遗传稳定性，收获其种子，编号D71157。2018年继续种植D71157观察其遗传稳定性，并进行联合品比试验，收获其种子，编号B81247。2019年进行联合品比试验，并进行米质分析，种子编号B1821，命名为湘辐1821。

1.4 测序幼苗样品取样

将黄华占和湘辐1821种子常规浸种3 d，发芽后的种子播种至带土塑料桶于人工气候箱中培养。至3叶期对地上部分取样，将样品送至北京诺禾致源生物信息科技有限公司进行后续工作。

1.5 DNA文库构建、测序和数据处理及分析

十六烷基三甲基溴化铵（cetyltriethylammnonium，CTAB）法提取DNA，质检合格的DNA样品通过Covaris破碎机随机打断成长度为350 bp的片段。采用TruSeq Library Construction Kit 试剂盒进行建库。DNA片段经末端修复、加ployA尾、加测序接头、纯化、PCR扩增等步骤完成整个文库制备。构建好的文库通过illumina进行测序。首先，对测序得到的原始序列进行过滤，去除带接头的序列对。其次，检查测序错误率分布，每个碱基位置的测序错误率都应低于1%。最终得到的过滤序列为有效测序数据。通过BWA［21］软件比对到参考基因组，比对结果经 SAMTOOLS［22］软件去除重复 ；采用 SAMTOOLS软件进行个体 SNPs 的检测［SNPs的 reads 支持数不低于4；SNPs 的质量值（mapping quality，MQ）不低于20］和InDels的检测；通过CNVnator［23］检测CNVs；利用 BreakDancer［24］软件检测SVs。基因本体（gene ontology，GO）分析首先把所有候选基因向Gene Ontology数据库（http://www.geneontology.org/）的各个分类条目映射，计算映射到每个term的靶基因数目，然后应用超几何检验，找出与整个参考基因背景相比，在候选靶基因中显著富集的GO条目，通过GO功能显著性富集分析能确定候选靶基因行使的主要生物学功能。

2 结果与分析

2.1 测序数据质量情况汇总

通过测序数据分析，本次测序共产生原始数据82.89 Gb，过滤后的有效数据 82.73 Gb，各样品有效数据大小在39 884.828～43 004.680 Mb之间。GC含量在43.25%～43.63%之间，所有样品的数据量足够，测序质量合格，GC分布正常，符合要求，可以进行后续分析。本文选用籼稻蜀恢498（R498）作为参考基因组，其大小为390 983 850 bp，所有样本的比对率在96.25%～97.18%之间， 对参考基因组的平均覆盖深度在87.79倍（X）～88.28倍（X）之间， 1倍（1X）覆盖度（至少有1个碱基的覆盖）在95.46% 以上。比对结果正常，可用于后续的变异检测及相关分析。表1是测序数据和比对结果小结。

表1 黄华占（WT_HHZ）和湘辐1821（MT_HHZ）测序数据小结Tab. 1 Summary of illunima data in Huanghuazhan（WT_HHZ）and Xiangfu 1821（MT_HHZ）

2.2 SNPs变异分析

采用SAMTOOLS软件进行个体SNPs的检测。黄华占（WT_HHZ）和湘辐1821（MT_HHZ）分别得到758 215和799 434个SNPs，基因间隔区分布的SNPs最多，多达481 765和508 409个。其次为内含子区域有69 806和72 355个，基因上游1 kb区域有64 274和68 140个。外显子则只有47 603和50 561个，以非同义突变为主要类型。湘辐1821的SNPs总数较黄华占多了41 219个，且不同分布区域的SNPs数值都高于黄华占。湘辐1821的SNPs的杂合率为0.236%，比黄华占高了0.021 %。根据碱基替换不同，SNPs分为转换（transitions，Ts）和颠换（transversions，Tv）。湘辐1821转换与颠换的比率为2.546，黄华占的转换与颠换的比率为2.542，二者基本持平（表2）。这2种基因型的SNPs变异基本以转换为主要类型。转换即是G/A 和 C/T的2种变异，颠换则有A/C、C/G、G/T和T/A 4种变异类型。在黄华占中，转换变异达到544 204个，颠换只有214 011个，湘辐1821中类似，转换达到574 019个，颠换则只有225 415个，二者都以转换为主要类型（图1）。

表2 SNPs分布区间及数目Tab. 2 Summary of SNPs discovery via illunima sequencing

图1 SNPs突变频谱图Fig. 1 The mutation spectrum of SNPs

2.3 InDels变异分析

利用SAMTOOLs软件检测长度小于50 bp的小片段的插入缺失（InDels），黄华占（WT_HHZ）/湘辐1821（MT_HHZ）一共得到142 313/147 686个InDels。基因间隔区分布的InDels最多，多达77 888/80 604个，其次为内含子区域有19 174/19 977个，基因上游1 kb区域有15 632/16 319个。湘辐1821外显子中的变异以移码变异为主，共有2 468个，非移码的有2 357个（表3）。基因组InDels长度从1 bp到21 bp不等，2种基因型都是1 bp的InDels为主要变异类型，占50%（图2）。编码区的InDels也是1～21 bp，1 bp的InDels在2种基因型中约占27%，2 bp的InDels约占11%，3 bp的InDels约占24%（图3）。在黄华占和湘辐1821中基本以短序列InDels为主，说明点突变为主要突变方式。

表3 InDels分布区间及数目统计Tab. 3 Summary of InDels discovery via illunima sequencing

2.4 SVs分析

SVs指基因组水平上INSs、DELs、INVs、染色体内部迁移（intra-chromosomal translocations，ITXs）、染色体间的迁移（inter-chromosomal translocations，CTXs）。黄华占中一共得到16 775个SVs，包括29个 INSs、5 751个 DELs、1 275个 INVs、2 354个ITXs 和7 366个CTXs。大部分 SVs位于外显子区域（2 542）和基因间隔区（2 917），基因上游1 kb区域有 527个， 下游1 kb区域有447个，内含子区域有334个， 8个 SVs影响剪接位点。在湘辐1821中一共得到18 382个SVs，比黄华占多了1 607个，包括29个 INSs、6 251个 DELs、1 318个 INVs、2 596 个 ITXs 和 8 188 个 CTXs。多数 SVs 位于外显子区域（2 648）和基因间隔区（3 211），基因上游1 kb区域有598个， 下游1 kb区域有463个，内含子区域有345个， 14个 SVs影响剪接位点。SVs的长度一般以100 bp为单位，长度大小从0～100 bp至 ＞1 200 bp。黄华占中58%的SVs和湘辐1821中57% 的SVs都是大于1 200 bp的序列（图4）。

图2 基因组InDels长度分布Fig. 2 The length distribution of the genome InDels

图3 CDS区InDels长度分布Fig. 3 The length distribution of the CDS InDels

2.5 CNVs拷贝数变异分析

CNVs 可以在转录和翻译水平上影响基因表达。利用CNVnator 软件分别检测到黄华占和湘辐1821中一共有16 614 和17 658个 CNVs。这些CNVs几乎都分布在基因间隔区，黄华占11 430个，湘辐1821 12 203个。外显子区域中黄华占2 090个，湘辐1821 2 238个。基因上游1 kb区域中黄华占1 092个，湘辐1821 1 118个。下游1 kb区域中黄华占826个，湘辐1821 866个。内含子区域中黄华占有552个，湘辐1821有 614个（图5）。黄华占中拷贝数增加个数是1 758个，拷贝数减少个数是14 856个，拷贝数总增加长度达到8 275 400 bp，拷贝数减少长度达到44 603 000 bp。湘辐1821中拷贝数增加个数是1 827个，拷贝数减少个数是15 831个，拷贝数总增加长度达到8 469 400 bp，拷贝数减少长度达到445 569 500 bp。2种基因型均以拷贝数减少为主要变异方式。

图4 SVs长度分布图Fig. 4 The length distribution of the SVs

图5 CNVs位置分布图Fig. 5 The CNVs distribution

2.6 变异的差异比较分析

SNPs和InDels变异可以直接地对应到基因的变异，利于突变体性状分析。通过对黄华占和湘辐1821中的SNPs和InDels进行差异分类统计，共找出差异SNPs 86 163个。其在12条染色体均有分布，其中第1染色体最多，达到11 389个，其次是第11染色体和12染色体。从分布区间来看，基因间隔区最多，有65 553个，外显子区域只有3 375个（表4、表5）。对这些差异SNPs类型进行分类统计，湘辐1821中转换类型的SNPs为38 121个，占总差异SNPs的44%；颠换类型的SNPs为13 100个，占总差异SNPs的15%。黄华占中转换类型的SNPs为35 666个，占总差异SNPs的41%；颠换类型的SNPs为12 470个，占总差异SNPs的14%。差异InDels共88 777个，其在12条染色体上也均有分布，其中第1染色体最多，达到13 261个，其次是第5染色体和第2染色体。从分布区间来看，基因间隔区最多，有48 045个，外显子只有2 843个（表4、表5）。

SNPs突变位点为单一碱基，对所涉及到的基因进行注释，筛选候选基因可靠方便。通过对这些差异SNPs所涉及基因进行筛选，筛选得到的候选基因有3 092个，并对其进行GO分析。3 092个候选基因分为生物学过程（biological process）、分子功能（molecular function）和细胞成分（cellular component）3大类。如图6所示：在生物学过程的17个小类中，大分子生物合成（macromolecule biosynthetic process）、细胞大分子生物合成（cellular macromolecule biosynthetic process）、大分子代谢调控（regulation of macromolecule metabolic process）、基因表达（gene expression）和DNA复制（DNA replication）类型基因较多，约占整个GO分析的70%；在分子功能的8个小类中，转移酶活性（transferase activity）、铁离子结合（iron ion binding）和DNA解旋酶活性（DNA helicase activity）的基因占主要类型；细胞成分包含5个小类，转录因子和RNA聚合酶相关基因较多。对这些差异候选基因的GO分析对湘辐1821表型分析可能具有重要参考意义。

表4 差异SNPs和InDels在染色体上的分布Tab. 4 The distribution of differential SNPs and InDels on chromosome

表5 差异SNPs和InDels在基因组上不同区间的分布Tab. 5 The distribution of differential SNPs and InDels on different genomic regions

3 讨论

图6 差异SNPs所在基因的GO分类Fig. 6 Gene ontology categories of differential SNPs associated gene

二代测序技术的高速发展推动了水稻基因组学的研究，高通量、高精度、低成本的超高优势使其迅速应用于结构基因组学、功能基因组学以及转录组学研究［25］。水稻基因组重测序是指在已知水稻基因组序列的前提下，对某个水稻品种基因组进行测序，对照参考基因组序列分析重测序品种变异情况，可挖掘大量的SNPs、InDels和SVs等，从而阐明该水稻品种的遗传特征。结合二代测序技术，对辐射诱变机理的探索也广泛开展起来，空间诱变获得的突变体主要发生DNA水平上的突变［26］。重离子辐射改变植物体内基因序列及表达方式，使突变体DNA甲基化、转座子/反转座子活性发生改变［27］。Cheng等［28］对从9311辐照来的突变体Red-1进行重测序，发现Red-1中9.19%的基因组序列发生改变，有381 403个SNPs，50 116个1～5 bp的InDels，在Red-1基因组中点突变变异方式是主要变异，这与本研究中的湘辐1821的变异方式一致，表明SNPs和InDels变异是γ射线辐射诱变的主要特征。

水稻育种已经从高产育种目标转变为兼顾优质、多抗、高氮利用率等绿色性状。核诱发突变技术能够诱发产生自然界稀有的新基因，大量的实践证明利用人工诱发遗传变异是丰富水稻种质资源和选育新品种的重要手段之一［29］。辐射诱变可以改变一个或两个主要性状，并保持其他性状不变。如经γ辐照后，水稻和马铃薯中的淀粉含量发生改变［30-31］。本研究利用γ射线辐射处理黄华占干种子，获得了综合性状良好、遗传稳定的突变新品系，较野生型其产量显著提高，剑叶增宽，株高增加，且米质经农业农村部稻米检测中心检测达到经农业农村部稻米检测中心颁布的一等食用稻标准。通过二代高通量测序，在该突变体中发现了大量分子变异，湘辐1821中变异类型的个数均高于野生型黄华占。湘辐1821的SNPs转换与颠换的比率为2.546，黄华占的为2.542，这2种基因型的SNPs变异基本以转换为主要类型。在黄华占和湘辐1821的InDels中基本以短序列InDels为主，说明点突变为主要突变方式。拷贝数减少是黄华占和湘辐1821的主要拷贝数变异方式。通过比较分析，2种基因型的差异SNPs 共有86 163个，差异InDels共88 777个。差异SNPs所涉及基因筛选到的候选基因有3 092个。这些研究结果将为湘辐1821表型分析提供重要参考，同时为辐射诱变机理提供理论支撑。