Preliminary Localization of a Rice Lesion Mimic Gene Based on BSA-seq

本实验中所用到的水稻材料有籼稻（Oryza sativa ssp.）‘黄华占’（‘HHZ’），以及经过EMS诱变筛选出‘HHZ’背景的LMM43突变体，已经过多代自交获得了稳定遗传的表型。

将LMM43突变体及‘HHZ’野生型（Wild Type, WT）水稻种植于海南省三亚市崖州区独村实验中心基地（109º18′N, 18º36′E），进行常规水肥管理。观察在全生育期内LMM43和WT植株的表型差异，记录LMM43突变体类病斑的出现情况，关注其斑点出现的时间点、形状、颜色、位置、密度等。在拔节期（60 d）取叶片用于BSA-seq测序的极端性状样本池构建。

为确定该突变体的遗传模式，以‘HHZ’作为母本，LMM43作为父本进行杂交，该工作于陵水试验基地完成。收获的F₁代种子，种植后观察F₁代植株表型，将杂株排除后自交获得F₂代种子。继续将F₂代群体种植于独村实验中心基地获得遗传分离群体，对其进行观察和取样工作，统计性状分离情况。

选取了白叶枯病菌的7个生理小种：PXO61、PXO86、PXO79、PXO71、PXO112、PXO99和PXO145，使用剪叶法分别对野生型和突变体进行接菌处理^[21]，每个生理小种接种3株水稻，20 d后对接种叶片的病斑长度进行拍照及测量。

在F₂代群体处于拔节期时，分别选择类病斑表型和野生型表型各30株以及LMM43突变体亲本和HHZ亲本各2株的叶片，提取DNA后等量混合构建为Mu-F₂、Mu-QB、WT-F₂、WT-QB 四个混池。然后送至杭州联川生物技术股份有限公司进行样品检测、文库构建、使用Illumina HiSeq 2500平台进行测序，测序深度为30×。

通过测序，得到最初的原始序列即Raw Reads。该数据里面含有带接头的和低质量的reads，通过去除接头、过滤掉reads上难确定的碱基类型占比大于10%的以及Q≤10碱基数占整reads50%以上的低质量reads，将低质量的数据过滤掉后得到Clean Reads。以明恢63（MH63RS3，大小为396 Mb）作为参考基因组，用BWA软件将测序获得的Reads定位至基因组上，通过比对定位Clean Reads在参考基因组上的位置，统计各样品的测序深度、基因组覆盖度等信息，并进行变异的检测。

将上一步比对到的结果使用samools继续过滤掉冗余reads，然后使用GATK对SNP和InDel进行变异检测，4个样本各自生成gVCF后，再进行群体间的基因型比对过滤，最终可以得到变异的位点集。将以上获得的位点使用SnpEff进行变异注释和预测变异影响。

在开展关联分析之前，先对SNP和InDel进行筛选，筛选条件包括：剔除具有多个基因型的位点；去除reads支持度低于4的位点；排除混池间基因型相同的位点以及隐性混池基因不是来源于隐性亲本的SNP位点。

将经上述条件过滤后的变异位点进行关联分析，有2种方法分别为欧式距离（Euclidean Distance, ED）算法和Index算法^[13]。ED算法是利用测序数据寻找混池间存在显著差异标记，并以此评估与性状关联区域的方法。理论上，BSA项目构建的两个混池间除了目标性状相关位点存在差异，其他位点均趋向于一致，因此非目标位点的ED值应趋向于0。ED值越大表明该标记在两混池间的差异越大^[22]。SNP-index近年来发表的一种通过混池间的基因型频率差异进行标记关联分析方法，主要是寻找混池之间基因型频率的显著差异。利用两亲本的SNP数据，分别计算两混池的SNP-index，并通过ΔSNP-index观测可能与性状分离相关的位点。通过以上两种算法分别对获得的SNP和InDel分别进行关联分析，得到候选区域。

通过将得到的SNP候选区域和InDel候选区域结果取交集，得到最终的候选区域。应用BLAST软件对候选区间内的编码基因进行多个数据库（NR、Swiss-Prot、GO、KEGG、COG）的深度注释。通过详细的注释，快速筛选候选基因。

在拔节期取LMM43和WT的叶片，使用天根生化科技有限公司的总RNA提取试剂盒（DP450）提取水稻叶片的RNA，送至公司进行转录组测序。

对WT和LMM43突变体进行全生育期观察，发现WT和LMM43在四叶期以前没有出现差异，从四叶期开始LMM43突变体开始从叶尖部位出现少量红褐色的坏死斑点（图1−A），随着叶片的生长，类病斑斑点数目逐渐增加从叶尖部位逐渐扩散至整个叶片及茎秆（图1−C—E），从幼苗到成熟都带有类病斑表型，类病斑在叶片上随机分布，为扩散的全生育期型突变体。通过对比剑叶、倒二叶、倒三叶发现类病斑在老叶先产生，随后在新叶产生。与WT相比，LMM43突变体叶片枯死，严重时导致整株植株死亡，表明类病斑严重影响了水稻的生长发育（图1−B）。

Figure 1.  Comparison of phenotypes between WT and LMM43

为了确定该突变体的遗传特性，本研究将野生型‘HHZ’作为母本， LMM43为父本进行杂交，通过观察F₁表型发现其均表现出类病斑，这说明了LMM43的类病斑性状受到显性基因控制。同时，对F₁进行自交后得到的F₂代水稻进行了观察，发现类病斑表型在F₂代群体中出现了分化，经统计，有斑单株与无斑单株分别为289株和86株，经过卡方检验（χ²=0.854<χ²_0.05=3.84，df=1），其分离比符合3∶1（表1）。综上确定了LMM43突变体表型是由单个核基因控制的显性性状。

对LMM43突变体进行白叶枯病菌的抗性鉴定，通过接种实验，发现突变体对白叶枯病菌的7个生理小种的抗病性均显著增强（图2）。

Figure 2.  Resistance characterization of WT and LMM43 leaf blight fungi

为了定位引起LMM43表型的基因，分别构建F₂类病斑表型（Mu-F₂）、LMM43亲本表型（Mu-QB）、F₂野生型（WT-F₂）、HHZ亲本（WT-QB）4个混池，通过BSA-seq重测序的方法，挖掘造成水稻类病斑突变体表型的基因位点。

使用Illumina HiSeq 2500测序平台进行双端测序， 经Raw Data过滤质控后共获得Clean Data共有为64.55 Gbp，Mu-F₂、Mu-QB、WT-F₂、WT-QB分别获得了16.32、18.40、16.90、12.93 Gbp的Clean Base，Q30平均值达到了90.96%。将得到的Clean Reads比对至MH63基因组上后，Mu-F₂、Mu-QB、WT-F₂、WT-QB样本池的比对效率分别为97.83%、97.50%、97.44%、97.32%，平均比对效率达97.52%（表2）。4个样本的平均覆盖深度为39×、44×、40×、31×，平均覆盖深度为38.50×（表2）。以上测序质量统计证明了此次测序质量高，一致性好，可以用于后续的生物信息学分析来定位候选基因。

用BAW进行比对，使用GATK和samtools（v1.9）进行检测和过滤，共发现1160877个SNP。在Mu-QB和WT-QB间检测，获得76489个SNP，造成非同义突变的有8120个；Mu-F₂和WT-F₂间检测SNP，获得128357个，非同义突变有12300个（表3）。在InDel检测中，亲本混池间共检测到20742个InDel，而F₂代混池间有32086个InDel（表4）。

利用F₂代的两个混池，将SNP和InDel分别通过ED和Index算法进行关联分析。在ED算法中先计算出每个位点的ED值再进行5次方处理以消除背景噪音，采用Distance法进行拟合后分别可以得到SNP-ED值分布图和InDel-ED值分布图（图3-A、图4−A），然后取所有位点拟合值的median+3SD作为分析的关联阈值（1.86、5.62）即得到8个SNP区域（图3−A中突起峰区域）和17个InDel区域（图4−A中突起峰区域）；在Index算法中计算出分别计算出WT-F₂和Mu-F₂每个位点的SNP-index值，然后两者相减得到ΔSNP-index值，Distance法进行拟合后得到SNP-index值分布图和InDel-index值分布图（图3B、图4B），随后取阈值线为置信度0.99后得到了9个SNP区域（图3−B ΔSNP-index值分布图中位于置信度0.99的阈值线上的突起峰）和26个InDel区域（图4−B ΔInDel-index值分布图中位于置信度0.99的阈值线上的突起峰）。将SNP经2种算法得到的8个和9个区域取交集后得到3个SNP关联区域572个位点（图3−C），InDel的17个和26个区域取交集后得到17个关联区域4 278个位点（图4−C）。

Figure 3.  SNP association analysis graphs

Figure 4.  InDel association analysis graphs

根据以上的关联结果再将获得的3个SNP关联区域和17个InDel关联区域取交集后得到了一个候选区域。在10号染色体13 830 000～15 200 000的1.37 Mb内（表5），共199个基因。其中根据这些变异位点的注释可知，亲本间存在的非同义突变20个，移码突变的InDel 2个；混池间非同义突变26个，移码突变的InDel 2个；所以在该候选区间内可能与性状密切相关的为这46个非同义突变的SNP位点和4个移码突变的InDel位点。通过对候选位点的编码基因进行多个数据库（NR、Swiss-Prot、GO、KEGG、COG）的深度注释，共注释到16个非同义突变基因，3个移码突变基因（表6）。

通过以上综合分析，得到了OsMH_10G0226700、OsMH_10G0225900等19个基因，其中，16个为非同义突变，3个为移码突变，鉴于EMS造成的主要为G到A和C到T的点突变，所以非同义突变是我们重点关注的。在这16个非同义突变中，其中5个基因可能与水稻的类病斑性状有关，OsMH_10G0223300, red chlorophyll catabolite reductase（PCCR）; OsMH_10G0206700, calmodulin-binding transcription activator 3 isoform X2（CAMTA3）; OsMH_10G0202100, ethylene-responsive transcription factor RAP2-13（RAP2-13）; OsMH_10G0201300, probable glucan 1,3-beta-glucosidase A（β-1,3-葡聚糖苷酶A）; OsMH_10G0201200, disease resistance protein RGA4-like（RGA4-like）（表6）。

为了解候选基因LMM43突变体中是否有表达量的变化，我们通过转录组测序结果分析5个候选基因的表达量变化，结果表明，只有OsMH_10G0223300和OsMH_10G0206700在水稻中的基因表达量较高，其中OsMH_10G0223300的表达量下降，但不显著；OsMH_10G0206700的基因表达量在LMM43突变体中显著上调。

Figure 5.  Candidate gene expression analysis

水稻类病斑自1965年被发现后，到目前为止已有近100多个基因被报道^[3]。研究表明，大部分水稻类病斑突变体对白叶枯病、稻瘟病或纹枯病具有抗性。水稻的抗病机制，特别是类病斑突变体的抗病机制，对于提高水稻的抗病性和产量具有重要意义。但是在水稻类病斑形成和激活免疫防御基因高表达的同时，类病斑突变体的生长发育同样受到了严重的影响，如突变体 lms1^[23]、spl28^[24]等不仅表现出类病斑的表型, 而且植株矮小，存在早衰现象，spl29^[25]，spl33 ^[26]等突变体的株高、穗长、结实率和千粒重与野生型相比均显著下降，本研究中的LMM43突变体类病斑属于扩散型的全生育期类病斑突变体，其生长发育和已报道过的类病斑突变体一样，类病斑的大面积的蔓延也会使得其农艺性状的下降，同时植株的矮缩也让LMM43的结实较为困难。

第1个突变体sl^[27]，是自然变异引起的，后来随着科学技术手段的发展，人们通过物理、化学及生物诱变型的多种手段获得突变体。目前已克隆近50个水稻类病斑突变体基因，其中大多数是通过图位克隆获得的^[28]，其中大部分为单个核基因控制隐性遗传，极少数的如显性基因遗传突变体有Spl18^[7]、spl-D^[9]、Spl26^[10]，半显性遗传突变体LIL1^[29]、nls1-1D^[8]，其中nls1-1D类病斑突变纯合体后代在生长方面存在着严重的缺陷，发芽后3周内就会死亡。所以显性水稻类病斑突变体需要构建F2:3群体使用纯合体进行定位，在F₂代杂合基因型干扰及纯合后代材料的致死现象，使得我们的基因定位上存在着较大困难。

BSA-seq是一种结合了分离体混合分析和高通量测序技术的方法，可以用于快速定位与特定性状相关的基因，该技术的核心在于选择具有目标性状极端表型的个体^[13]。所以选择了SBA-seq来对LMM43类病斑突变体的基因进行初步定位，由于使用F₂代群体构建极端性状混池，其中有显性杂合体杂的存在，使得干扰较多，最终通过分析，将候选区域确定在了10号染色体的1.37 Mb的一个区间内，其中包括199个基因。根据目前报道过已经克隆的10号染色体上共有5个类病斑突变体，其中4个隐性突变，1个显性突变，都不在我们的候选区域内^[3]，证明了LMM43是由一个尚未报道全新的显性基因控制的类病斑突变体。在此基础上我们通过经候选区域内的基因经多个数据库深度注释后确定了5个候选基因，OsMH_10G0223300编码的PCCR，是植物中叶绿素降解的关键酶，Tang等发现敲除OsPCCR1会使水稻产生类病斑^[30]；OsMH_10G0206700编码的CAMTA3，是一类转录因子可感知来自Ca²⁺的应激信号，具有多种胁迫响应性^[31]；OsMH_10G0202100编码的RAP2-13，作为乙烯响应因子（ERF）家族，可能参与了水稻的抗病等过程；OsMH_10G0201300编码的β-1,3-葡聚糖苷酶A主要参与病原体防御、碳水化合物代谢和逆境响应等过程^[32]；OsMH_10G0201200编码的RGA4-like，是RGA4的同源蛋白，在Zhu等研究中发现了RGA4在水稻中通过与RGA5协同作用，识别稻瘟菌的效应蛋白并激活抗病反应^[33]。所以综上这些基因可能与类病斑性状的出现相关。通过对候选基因的表达量进行分析发现只有2个基因在水稻中的表达量较高，OsMH_10G0223300的表达量下降，OsMH_10G0206700显著上调，由于碱基突变造成的非同义突变不一定会影响该基因的表达量变化，所以后续需要继续通过互补实验进行验证。

本研究从籼稻‘黄华占’（‘HHZ’）经EMS诱变后的群体中分离了1个水稻类病斑突变体LMM43，通过对LMM43类病斑突变体的表型进行全生育期的观察，是1个全生育期扩散型的类病斑突变体，且LMM43突变体对白叶枯病菌的抗性显著增强。在F₂代群体中，通过统计性状分离情况，得知LMM43突变体表型是由单个核基因控制的显性性状。利用BSA-seq对LMM43类病斑突变体的基因进行初步定位。在‘HHZ’×LMM43杂交构建的F₂代分离群体中进行混池构建，将基因定位在水稻10号染色体的1个约1.37 Mb的区间内，该区间包含199个预测基因，综合分析预测出5个候选基因，为后期基因克隆和功能研究奠定了基础。