一种利用Δ(SNP-index)进行性状定位的QTL-seq方法及其应用

摘要

一种利用|Δ(SNP‑index)|进行性状定位的QTL‑seq方法及其应用，该方法以参考基因组代替亲本之一的基因组作为参照计算子代池SNP‑index，用Δ(SNP‑index)的绝对值|Δ(SNP‑index)|代替Δ(SNP‑index)用于滑动窗口法性状定位分析，通过对双亲具有相同表型的分离群体的极端表型池进行QTL‑seq分析，并利用|Δ(SNP‑index)|的分布进行目标性状基因或主效QTL的定位预测，改变了常规QTL‑seq方法只能应用于双亲纯合且具有相对性状的分离群体的状况，使QTL‑seq方法能够应用于双亲具有相同性状的分离群体，极大地拓展了QTL‑seq方法的应用范围，具有很好的应用前景。应用该方法对双亲均为红皮东方梨的“满天红”ד红香酥”杂交组合红皮表型池和绿皮表型池进行分析，将东方梨红皮/绿皮性状位点定位在1.86Mb的区间内。

说明书

技术领域

本发明属于生物信息学和基因组学领域，涉及一种利用|Δ(SNP-index)|进行性状定位的QTL-seq方法及其应用。

背景技术

QTL-seq是近年发展起来的采用高通量测序技术进行集群分离分析(BSA)定位数量性状主效基因或质量性状基因的一种方法(Takagi et al.,2013)，适用于双亲具有一对相对性状的子代群体。该方法选择常用作图群体中20-50个表现极端表型的单株组成极端表型池进行全基因测序，以双亲之一的基因组为参照对两个极端表型池进行分析，寻找两个池中的SNP，计算每个SNP的SNP-index，根据SNP-index的分布进行QTL或基因预测。

多年生果树大多基因组高度杂合，许多育种群体通常利用在目标性状具有相同表型的材料作为亲本配置杂交组合以获得更大比例具有目标性状表型的杂种后代。利用这样的育种群体进行目标性状的QTL-seq定位分析，需要对QTL-seq方法进行改进。

发明内容

解决的技术问题

为了解决QTL-seq方法无法在双亲具有相同目标性状表型的果树分离群体中进行应用的问题，本发明提供了一种利用|Δ(SNP-index)|进行性状定位的QTL-seq方法及其应用。

技术方案

一种利用|Δ(SNP-index)|进行性状定位的QTL-seq方法，其中|Δ(SNP-index)|是Δ(SNP-index)的绝对值，本发明方法以参考基因组作为参照，对双亲具有相同表型的分离群体的两个相对表型极端表型池进行QTL-seq分析，并利用|Δ(SNP-index)|的分布进行目标性状基因或主效QTL的定位预测。

进一步地，本发明方法利用|Δ(SNP-index)|的分布进行目标性状基因或主效QTL的定位预测时，|Δ(SNP-index)|规避了用参考基因组代替亲本基因组作为参照的情况下，SNP位点间连锁关系不清导致的滑动窗口内SNP位点的平均Δ(SNP-index)为0，无法进行目标性状基因或主效QTL定位预测的情况。

另一方面，本发明利用|Δ(SNP-index)|进行性状定位的QTL-seq方法在东方梨红皮/绿皮性状位点定位中的应用，利用本发明方法对双亲均为红皮梨的“满天红”ד红香酥”杂交组合红皮池和绿皮池进行分析，能够将所述“满天红”ד红香酥”杂交组合的东方梨红皮/绿皮性状位点定位在1.86Mb的区间内。

进一步地，上述“满天红”ד红香酥”杂交组合的东方梨红皮/绿皮性状位点所在的1.86Mb的区间包括如下14个子区间：

Scaffold NW_008988041.1上3800001-3880001bp区间；

Scaffold NW_008988076.1上640001-660001bp区间；

Scaffold NW_008988091.1上1-140001bp区间；

Scaffold NW_008988126.1上560001-720001bp区间；

Scaffold NW_008988130.1上1-300001bp区间；

Scaffold NW_008988130.1上500001-580001bp区间；

Scaffold NW_008988130.1上760001-960001bp区间；

Scaffold NW_008988141.1上1-500001bp区间；

Scaffold NW_008988141.1上880001-940001bp区间；

Scaffold NW_008988461.1上240001-400001bp区间；

Scaffold NW_008988478.1上300001-400001bp区间；

Scaffold NW_008988581.1上300001-320001bp区间；

Scaffold NW_008989660.1上1-20001bp区间；

Scaffold NW_008989715.1上1-20001bp区间。

有益效果

“满天红”和“红香酥”是中国农业科学院郑州果树研究所选育的两个著名的红皮东方梨品种，在我国推广应用面积较大。这两个品种常被育种单位用作红皮梨育种的亲本，表现出良好的亲本遗传特点。本发明利用改进的QTL-seq方法对“满天红”ד红香酥”杂交组合群体红皮极端表型个体组成的红皮表型池和绿皮极端表型个体组成的绿皮表型池进行分析，成功地将东方梨红皮/绿皮性状位点定位到1.86Mb的区间内，为进一步开发可以大幅提高红皮梨育种选择效率的实用分子标记奠定了坚实的基础。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图，其中：

图1是现有技术中SNP-index计算原理图(Takagi et al.,2013)；

图2-1表示Scaffold NW89880130.1Δ(SNP-index)曼哈顿图(横坐标是Scaffold物理位置，纵坐标是Δ(SNP-index))；以及

图2-2表示Scaffold NW89880130.1|Δ(SNP-index)|曼哈顿图(横坐标是Scaffold物理位置，纵坐标是|Δ(SNP-index)|)。

具体实施方式

下面结合具体实施例对本发明做进一步的说明，但本发明不受以下实施例的限制。

根据本发明的示例性实施例，SNP-index根据reads测序深度信息计算，利用测序reads对每个碱基位点的碱基进行统计，计算与作为参照的基因组不相同的reads条数占总条数的比例，即为该碱基位点的SNP-index。完全与参照基因组相同的SNP位点其SNP-index为0；完全与参照基因组不同的SNP-index为1。Δ(SNP-index)是两个极端表型池中SNP位点的SNP-index差值。

通常，在双亲纯合且具有一对相对性状的群体中，利用亲本信息可以判断 子代SNP等位基因型的连锁状态，Δ(SNP-index)的分布是QTL-seq方法进行基因或主效QTL区间预测的最核心的指标，其符号的正负方向分别对应目标性状的一对相对性状。但若双亲杂合或不以亲本基因组作为计算SNP-index时的参照，则Δ(SNP-index)符号的正负方向无法与目标性状的一对相对性状分别对应。因而在双亲具有相同表型，子代群体中出现性状分离的情况下，一般也不建议以亲本之一的基因组数据为参照进行QTL-seq分析。

实施例1

1 材料来源

供试材料为中国农业科学院郑州果树研究所新乡试验基地定植的“满天红”ד红香酥”杂交组合群体极端表型单株中28个红皮单株组成的红皮池和27个绿皮单株组成的绿皮池。

2 方法

2.1 建库测序

对红皮池DNA样品和绿皮池DNA样品分别由池中单株DNA等量混合而成。两池各取1.5ug DNA通过Covaris破碎机随机打断成长度为350bp的片段，用于制备测序样品文库。测序文库的制备采用Truseq Nano DNA HT Sample preparation Kit(Illumina USA)试剂盒，按照说明书进行DNA片段末端修复、加ployA尾、加测序接头、纯化、PCR扩增等步骤后完成整个文库制备。文库先用Qubit2.0初步定量并稀释至1ng/μl，用Agilent 2100对文库插入片段大小(insert size)进行检测。插入片段大小符合预期后用Q-PCR方法对文库的有效浓度进行准确定量，确保待测序文库有效浓度>2nM，以保证文库质量。测序文库利用Illumina HiSeq2500测序平台进行125bp双端(paired-end)测序，测序读长(reads)约350bp。

为了保证信息分析质量，对fastq格式的原始测序数据进行过滤，包括去除带接头(adapter)的reads pair；去除单端测序read中未识别碱基(N)含量超过该条read长度比例10％的paired reads；以及去除单端测序read中低质量(Q≤5)碱基数超过该条read长度比例50％的paired reads等步骤。

有效测序数据通过BWA(Burrows-Wheeler Aligner)软件(Li et al.,2009)，比对到参考基因组(ftp://ftp.ncbi.nlm.nih.gov/genomes/Pyrus_x_bretschneideri/CHR_Un/)，参数设置：mem-t 4-k 32–M。比对结果文件用SAMTOOLS软件(Li et al.,2009)转换为BAM文件，参数设置：–bS–t，并利用SAMTOOLS软件的“rmdup”命令去除重复。如果多个read pairs具有相同的外部坐标，只保留比对质量最高的read pair。

2.2 SNP/InDel检测和注释

采用GATK3.3软件(McKenna et al.,2010)的UnifiedGenotyper模块进行样本SNP/InDel的检测，使用VariantFiltration进行SNP过滤的参数为—filterExpression"QD<4.0||FS>60.0||MQ<40.0"，-G-filter"GQ<20"，--clusterWindowSize 4。使用VariantFiltration进行InDel过滤的参数为 --filterExpression"QD<4.0||FS>200.0||ReadPosRankSum<-20.0||InbreedingCoeff<-0.8"。SNP注释用ANNOVAR软件(Wang et al.,2010)基于参考基因组GFF3注释文件进行(ftp://ftp.ncbi.nlm.nih.gov/genomes/Pyrus_x_bretschneideri/GFF/)。

2.3 SNP-Index计算

SNP-index根据其reads测序深度信息计算(Takagi,Abe et al.2013)，利用测序reads对每个碱基位点的碱基进行统计，因为双亲在目标性状上表现相同，均为红皮，所以本研究以参考基因组(Wu et al.,2013)为参照，统计子代池(Red-Pool和Green-Pool)中和参考基因组在某一个碱基位点相同或者不相同的reads条数，计算不相同reads条数占总条数的比例，即为该碱基位点的SNP-index。完全与参考基因组相同的SNP-index为0；完全与参考基因组不同的SNP-index为1。其原理如附图1所示。Δ(SNP-index)是每个SNP在绿皮池中的SNP-index减去红皮池中的SNP-index的差值。

为减少测序错误和比对错误造成的影响，我们将两个子代池中SNP-index都小于0.3，并且SNP深度都小于7的位点，以及一个子代池SNP index缺失的位点，均予以过滤去除。

2.4 曼哈顿图绘制

QTL-seq策略适用来源于纯合亲本的子代群体，利用亲本信息可以判断子代SNP等位基因型的连锁状态。若子代群体来源于杂合亲本，则相邻SNP位点等位基因的连锁状态无法根据亲本之一或参考基因组进行判断，利用滑动窗口方 法进行区间定位时可能会因为正链负链无法判断的情况而导致错误(注：正链和负链分别指2条同源染色体)，出现窗口Δ(SNP-index)趋向于0的“假交换效应”。

用滑动窗口法绘图，每个窗口中所有SNP位点的Δ(SNP-index)和|Δ(SNP-index)|的平均值即为该窗口的Δ(SNP-index)和|Δ(SNP-index)|。选择200Kb为滑动窗口，20Kb为步长，以Scaffold(梨尚未公开完整的染色体序列，以Scaffold代替)的物理位置为横坐标，窗口Δ(SNP-index)和|Δ(SNP-index)|为纵坐标绘图。

2.5 候选区间确定

按照大约1500cM的梨图谱平均长度和性状粗定位的要求，我们根据滑动窗口法绘制的|Δ(SNP-index)|曼哈顿图，以Top 0.5％作为阈值(此阈值之上SNP区间内的Δ(SNP-index)绝对值较大，数量占所有SNP数量的0.5％)，划定候选基因区间，跨度大约为1500cM×0.5％＝7.5cM。该区间范围内的SNP作为候选SNP位点。

2.6 候选区间验证

在候选区间内设计InDel和SSR引物，先在9个红皮梨和7个绿皮梨中进行筛选，能较好地区分红皮梨和绿皮梨的引物用于对“满天红”ד红香酥”群体所有单株进行连锁分析，以验证区间定位的准确性。

3 结果

3.1 测序结果

测序数据的Q20(％)和Q30(％)分别在93％和86％以上，红皮池和绿皮池的平均测序深度分别为19.65X和22.98X，基因组覆盖度分别为94.74％(1X)和95.15％(1X)，84.38％(4X)和86.39％(4X)。

3.2 SNP检测及注释

共检测注释4,134,489个SNP，其中2,812,500个SNP突变属于转换，1,321,989个SNP突变属于颠换，转换和颠换的比例为2.127。4,134,489个SNP中，位于外显子区的有556,583个SNP，包括使基因获得终止子的2,749个SNP，使基因失去终止子的401个SNP，以及288,384个同义突变SNP和265,049个非同义突变SNP。

3.3 SNP-index曼哈顿图绘制

以Scaffold NW89880130.1作为绘图Scaffold来说明按Δ(SNP-index)绘图和按|Δ(SNP-index)|绘图的效果。未对Δ(SNP-index)取绝对值绘制的附图2-1和对Δ(SNP-index)取绝对值绘制的附图2-2表现出明显的不同，图2-1中窗口Δ(SNP-index)均＜0，从前往后总的趋势是趋向于0，图2-2因为是取了绝对值后再做的滑窗分析，所以窗口|Δ(SNP-index)|均＞0，在窗口Δ(SNP-index)趋向0的位置，窗口|Δ(SNP-index)|反而距离0较远，说明对窗口内的Δ(SNP-index)进行平均时出现了正负值相互抵消的情况，这种抵消效应很好地证明了正负链不清、标记间连锁状态不明造成的影响。

3.4 候选区间确定

以200Kb为窗口，20Kb为步长绘制|Δ(SNP-index)|曼哈顿图，设置Top0.5％为阈值线，该阈值线为|Δ(SNP-index)|＝0.33。在|Δ(SNP-index)|＝0.33的阈值线之上划定东方梨红皮/绿皮性状候选区间，共发现14个子区间，合计1.86Mb。14个子区间如下：

Scaffold NW_008988041.1上3800001-3880001bp区间；

Scaffold NW_008988076.1上640001-660001bp区间；

Scaffold NW_008988091.1上1-140001bp区间；

Scaffold NW_008988126.1上560001-720001bp区间；

Scaffold NW_008988130.1上1-300001bp区间；

Scaffold NW_008988130.1上500001-580001bp区间；

Scaffold NW_008988130.1上760001-960001bp区间；

Scaffold NW_008988141.1上1-500001bp区间；

Scaffold NW_008988141.1上880001-940001bp区间；

Scaffold NW_008988461.1上240001-400001bp区间；

Scaffold NW_008988478.1上300001-400001bp区间；

Scaffold NW_008988581.1上300001-320001bp区间；

Scaffold NW_008989660.1上1-20001bp区间；

Scaffold NW_008989715.1上1-20001bp区间。

显然，图2-2中的Scaffold NW_008988130.1上1-300001bp区间、 500001-580001bp区间和760001-960001bp区间均不能从利用Δ(SNP-index)进行滑窗扫描的图2-1中被定位出来。

3.5 候选区间验证

候选区间内设计的InDel和SSR引物经筛选后，对“满天红”ד红香酥”群体所有单株进行连锁分析，发现区间内相关InDel和SSR引物与东方梨红皮/绿皮性状位点具有明显的连锁关系，其中位于ScaffoldNW_008988130.1上的In2130-12、In2130-16、In2130-19和ZFRI130-16等标记与目标性状位点遗传距离均在2.6cM以内。

连锁分析验证结果表明，以参考基因组代替亲本基因组作为参照，计算SNP-index，以|Δ(SNP-index)|代替Δ(SNP-index)进行滑动窗口扫描的改良QTL-seq方法在双亲均为红皮梨(亲本表型相同)的东方梨红皮/绿皮性状分离群体中定位的东方梨红皮/绿皮性状候选区间是准确的。而用Δ(SNP-index)进行滑动窗口扫描的传统QTL-seq方法则不能定位出相应区间。

以上所述，仅为本发明的优选实施方式，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，对本发明技术方案进行的简单变化或等效替换均落入本发明的保护范围。