深度重测序揭示等位基因变异胡麻属indicum

背景

等位基因变异和连锁不平衡的全基因组模式的特征可以用来检测可靠的表型-基因型关联和分子选择的特征。然而，使用胡麻属indicum由于缺乏多态性数据，育种的种质资源受到限制。

结果

在这里，我们描述了来自12个国家的29个芝麻菌株的大规模并行重测序，每个测试样本的覆盖深度≥13倍。每个样本平均检测到127347个SNPs, 17961个小InDels和9266个结构变异。种群SNP率、种群多样性(π)和分离位点的Watterson估计量(θw)估计为8.6 × 10³2.5 × 10³3.0 × 10³英国石油公司^－1,分别。这些snp中，23.2%位于编码区内。多态性模式在基因家族中是非随机的，与生物或非生物环境相互作用的基因表现出高水平的多态性。连锁不平衡(LD)衰减距离估计为150kb，在种群中未观察到明显的结构。29个芝麻品系之间的系统发育关系符合芝麻起源于印度次大陆的假设。此外，我们还提出了腺苷酸异戊烯基转移酶(ITP)基因通过介导玉米素生物合成来决定芝麻叶腋花数的新作用。

结论

该研究首次报道了芝麻遗传变异的全基因组模式。这里描述的高LD距离和丰富的多态性增加了我们对影响芝麻全群体序列变异的力量的理解，这将是未来基因-表型和全基因组关联研究(GWAS)的宝贵资源。

胡麻属indicum(芝麻)是一种古老的作物，基因组大小约为357 Mb，含有高浓度的油脂和蛋白质，具有药用价值。然而，这种植物容易遭受内涝，特别容易受到许多真菌和细菌疾病的影响，如茎腐和根腐病，镰刀菌素枯萎病，白粉病等。这些生物和非生物胁迫会导致总体产量下降，产量与生长条件密切相关。为了克服环境压力和提高产量，植物育种项目需要丰富的种质资源和遗传信息[1]，等位基因变异和连锁不平衡的全基因组模式的表征确保了检测可靠的表型-基因型关联和分子选择的特征[2］．印度、中国和韩国是芝麻种质收集、保存和研究的领先国家[3.］．在中国，约6000株芝麻已存放在中国国家基因库(武汉，中期基因库;北京,长期资源库)。在韩国，位于韩国水原的农村振兴厅(RDA)的基因库中保存了> 7,698个变体[4]，在印度，国家植物遗传资源局(印度新德里)已存档了> 1万个变体。然而，由于基因组信息的缺乏和合适的生物标记物的缺乏，在全基因组尺度上对芝麻种质资源遗传多样性的研究较少。1]、[5] - [7］．

芝麻是该属中最常见的品种胡麻属，园内有超过20种开花植物。与芝麻不同的是，该属的大多数物种都没有被驯化，在多倍体水平上存在显著差异[1］．由于这些野生物种大多原产于撒哈拉以南的非洲，芝麻最初被认为起源于非洲;然而，驯化的芝麻被证明起源于印度次大陆。8]、[9］．由于缺乏多个芝麻品系的详细分子数据，对芝麻进化的进一步研究受到了阻碍。

参考基因组的完成为探索野生和驯化芝麻的遗传变异提供了必要的资源美国indicum（http://www.ocri-genomics.org/Sinbase/)．我们分析了来自12个国家的29个重测序芝麻品系，覆盖≥13倍。根据这些数据，我们构建了首个芝麻单倍型图谱，为了解芝麻跨多品系的遗传多样性提供了依据。这些数据可用于全基因组关联研究的发展，进而促进与简单和复杂性状相关的基因的定位。

重测序芝麻品系的表型多样性

我们手工选择了29个芝麻品系进行基因组重测序，其中6个来自印度和非洲，16个来自中国，2个来自美国，1个来自阿富汗、阿拉伯联合酋长国、韩国、缅甸、菲律宾和越南1:数据S1)。这些品系表现出广泛的表型，包括确定和不确定的生长习性、株高和株高的高低、开花的早和晚、种皮颜色的不同、每叶腋开一朵花和三朵花、单株和分枝风格等。遥远的地理关系和广泛的表型变异使这些菌株成为探索栽培芝麻遗传多样性的理想模型(图1)．

芝麻品种遗传多样性景观研究

为了识别大规模的多态性，更好地了解芝麻种质的遗传结构，对29个芝麻品系进行了重新测序，每个品系的覆盖深度≥13×，过滤后的数据超过120gb(附加文件2:图S1;额外的文件1:数据S1)。所有的序列都与“中子13号”的参考基因组进行比对，该基因组的有效长度为274 Mb (http://www.ocri-genomics.org/Sinbase/)，使用BWA软件[10］．不同品系间的映射率在88.8% ~ 95.2%之间，平均为91.4%。映射结果与GATK软件的映射结果一致1:数据S1)。

使用严格的流水线，我们使用SAMtools软件平均识别出每个菌株127,347个单核苷酸多态性(SNPs) [11]，范围为40,925至392,5441；额外的文件2:图S2)。总体而言，SNP率在1.5 × 10之间⁴到14.3 × 10⁴，其中G:A, A:G, C:T和T:C替换是最常见的(附加文件2:图S3)。通过组合所有菌株的SNP，我们鉴定出2,348,008个独特SNP，群体SNP率为8.6 × 10³英国石油公司^－1．接下来我们使用GATK软件调用并计算SNP种群，结果共生成2,003,821个种群SNP。GATK和SAMtools的符合率在80.0% ~ 89.9%之间，16个连锁组(假分子染色体)的符合率平均为85.4%1:数据S2)。Sanger测序结果显示，SNP呼叫准确率约为93.7%2:图S4;额外的文件1:数据S3)。这些结果表明，这里检测到的大多数snp是可靠的。在这些snp中，25.1%位于基因编码区，1.5%、7.5%和1.7%分别位于5 '非翻译区(UTR)、编码序列(CDS)和3'UTR2:图S5;额外的文件1:数据S4)。分别来自印度、缅甸和阿联酋的菌株20、24和26的SNP阳性率居前三;因此，这些地区可能蕴藏着更多样化的芝麻种质资源。

尽管芝麻传统上被认为是自花授粉的植物，但它也受到蝴蝶和蜜蜂等昆虫的交叉授粉。这与观察到的杂合度为12.37 ~ 49.30%，平均25.39%的杂合度率相一致2:图S6;额外的文件1:数据S5)。16个中国品种中有5个(品系1 ~ 5)杂合子率低于平均值(16.82 ~ 23.25%)和其他非中国品系。菌株16、22和26异常高的杂合子率表明这些菌株对昆虫更有吸引力;然而，更琐碎的解释，如测序和对齐错误，不能排除，特别是在重复区域。

种群snp用于计算两个常用的种群遗传统计数据:种群多样性(π)和分离位点的Watterson估计(θw)．平均π而且θw所有29个品系的值分别为2.5和3.0 / kb，低于水稻[12但比鹰嘴豆高(中投arietinum)[13)、西瓜(Citrullus lanatus)[14和大豆[15)(附加文件1:数据S6)。我们观察到许多连续滑动窗口伴随着假分子(LG1-LG16)，这些假分子包含的snp比正常snp少，反过来也更低π而且θw值(图2)，表明遗传多样性沿芝麻伪分子分布不均匀。我们检查了基因组中最高和最低的多态区域(位于顶部和底部5%的块)π值)，并发现最高多态区域的基因数量小于最低区域(524 vs. 1308)(附加文件1:资料S7和S8)，类似于其他物种，如鹰嘴豆[13］．在最高多态区域的许多基因与环境适应性有关，包括应激反应途径1:数据S9)。这些基因可能为研究芝麻的生物和非生物胁迫提供有价值的资源。值得注意的是，尽管在最低多态性区域的基因数量更多，但只有5个基因在两个基因本体(GO)术语中富集，所有这些基因都与基本生物功能有关;即核糖体结合(附加文件1:数据S10)。

接下来，我们使用满足所有对端要求的映射读数，但在contig的一端包含对齐间隙，以检测每个菌株的短InDels (1 - 5 bp)。检测到的InDel的总数与InDel(附加文件)的长度成反比2:图S7)。在所有29个品系中鉴定出的数量从4,495 - 56,594(平均为17,961)，总共有520,880个独特的InDels1:数据S11)。与SNP相似，InDels在基因组中的分布并不均匀，高密度区域与SNP高密度区域密切相关(图2)．在这些indel中，插入和删除的数量相似(分别为48.8%和51.2%)。纯合子内嵌子的比例是杂合子内嵌子的1.5倍以上。其中71.7%位于基因间区，1.5%(8221)位于CDS区，5.0%位于utr区。

结构变异(SV)最初定义为插入、删除、DNA倒置和其他大于1kb大小的序列重排[16］．测序现在已成为常规[17]，结构变异的操作范围已扩大到包括小得多的事件[18]、[19］．在本研究中，我们使用软件包Breakdancer v1.2检测了10 bp到1 Mb之间的SVs [20.设置为默认参数。与参考基因组相比，我们在所有29个菌株中发现了7220 - 12,458个SVs(平均= 9,266)1:数据S12)。对于这些sv，删除事件的数量几乎是插入事件的二比一(附加文件2:图S8)。除InDels外，其他sv(包括DNA倒置、染色体内易位和染色体间易位)的发生率相对较低，在739 ~ 2360之间(平均= 1140)。大多数sv的大小在100 - 1000 bp之间，较长的变异(>1 kb)较少，尤其是大于10 kb的变异(附加文件)2:图S9)，与高粱的情况一致[21］．

变异对基因的影响

基因内DNA序列的变化对植物的形态和进化起着至关重要的作用。芝麻基因组中有27148个注释基因(http://www.ocri-genomics.org/Sinbase/)，与参考基因组相比，74.8%(20,311)包含一个或多个snp。此外，62.6%(16,997)、15.5%(4,218)和18.0%(4,892)的基因在CDS、5 ' utr和3 ' utr中分别含有SNPs。这些基因被分为43个分子功能基团，其中30%与蛋白质结合、水解酶活性和ATP结合有关;然而，所有预测水解酶活性的基因仅在其CDS区域内包含snp2:图S10)。进一步的分析发现，在其CDS区域有258个snp基因显著富集(P< 0.01)为生物过程细胞死亡和凋亡过程(附加文件1:数据向)。在编码区识别出的136,130个非同义snp和142,103个同义snp表示非同义与同义替换比为0.99(附加文件1: S4数据;额外的文件2:图S11)，与高粱(1.0)相似[22，但高于拟南芥(0.83) (23低于大豆(1.38)[15]和大米(1.2)[24］．具有非同义SNPs基因的氧化石墨烯术语富集与细胞死亡、凋亡和防御反应密切相关(附加文件1:数据S14)，特别是编码抗病蛋白、udp -糖基转移酶或含有富含亮氨酸重复序列和NB-ARC结构域的蛋白(附加文件)2:图S12;额外的文件1:数据S15)。这些结果表明，参与生物胁迫反应的基因突变率较高，这与植物-病原体相互作用导致病原体相关分子模式识别受体多样化的理论一致[25]、[26］．

位于关键结构位置的编码区snp可以导致蛋白质形态的显著变化，进而导致蛋白质整体功能的变化。在研究的29个芝麻菌株中，我们发现了1281个与过早停止密码子形成相关的snp和246个与停止密码子到不停止密码子突变相关的snp。在186个基因中观察到起始密码子到非起始密码子突变，以及额外的404个剪接位点突变2:图向)。这些大效应snp大多位于假分子(LG)的近端(图)2)．这四种大效应snp的注释显示了不同的功能富集模式。例如，起始密码子到非起始密码子突变主要出现在涉及转运、凋亡和防御反应的基因中，而剪接位点突变则更常见于与细胞代谢、氧化还原、有机物质代谢和氮化合物代谢相关的基因中2:图S14系列)。在四种大效应snp中，过早停止密码子尤其有趣，因为这些突变通常与功能丧失有关。大多数突变是在与氧化石墨烯生物过程相关的基因中发现的，这些生物过程与逆境有关，包括细胞死亡、凋亡和防御反应1:数据S16)。

尽管大多数SNPs在CDS区域被检测到，但CDS区域仅占12651个InDel突变的14.3%，低于5 '和3 ' utr(分别为18.4%和19.6%)。5’和3’UTRs中包含InDels的基因数量随着InDel大小从1到5 bp的增加而减少，但在3 bp的InDels中，CDS InDels显著增加，这与在高粱中观察到的相似[21)(附加文件2:图S15)。这种3-bp InDels的富集是意料之中的，因为InDels如果不是3-bp的倍数，就会导致帧移位，通常是致命的。最后，我们在基因本体的基础上分析了InDels的分布，发现它们与导致过早停止密码子的SNPs相似，且富集具有统计学意义(P< 0.001)的基因参与细胞死亡，凋亡和防御反应(附加文件1:数据肌力表现)。

29个芝麻品系间的亲缘关系

考虑到芝麻最初是在非洲还是在印度次大陆种植的时候[8]、[9因此，研究地理环境对芝麻遗传多样性的影响具有重要意义。采用邻居联结法构建了包含29个芝麻品系的系统发育树。该分析揭示了中国菌株之间的最高亲缘度，而来自其他国家的菌株则遍布全球(图3.a).基于主成分分析(PCA)，来自不同地理位置的芝麻品系的这种相互交织的特性也很明显(图)3.b和c).使用贝叶斯聚类软件STRUCTURE观察模糊组[27]， K由2 - 5逐渐变化(图4一个)。

由于本研究没有包括芝麻的任何亲缘种或野生种，因此不可能对芝麻的起源做出明确的结论。然而，观察到的29个芝麻品系之间的系统发育关系为芝麻的进化提供了一些线索。这三株来自印度、缅甸和阿拉伯联合酋长国的菌株相对于其他菌株表现出更高的遗传距离(图3.a).根据多样性理论的瓦维洛夫中心(Vavilov center of diversity theory)，该理论认为，在植物最初被驯化的地方，可以观察到更丰富的遗传多样性[28，这些结果表明芝麻起源于印度次大陆。

芝麻的高度连锁不平衡

在设计关联研究时，LD模式对于确定映射分辨率是必要的[29]、[30.和解释联想峰值[31］．为了估计芝麻的LD，我们计算了r²使用Haploview在snp对之间进行分析[32]，并发现在大约150 kb的过程中，它从初始值0.30衰减到~0.15(图4b和c)。芝麻的LD衰减估计与自花大豆相当(约150kb) [15，但比在答:芥(~ 4 kb) (29]，籼稻(约65 kb) [12)和谷子(约100 kb)。也显著高于高粱等异花授粉植物(1 kb) [33]和玉米(< 1kb) [34］．芝麻的高LD不仅使它成为研究极端LD对基因组和群体结构影响的好植物[15]，但也适用于多态标记相对较少的GWAS。

芝麻叶腋花数候选位点的聚类分离分析

大块分离分析(BSA)是一种快速检测特定基因组区域标记的方法[35]，并已成功应用于检测水稻各种性状的数量性状位点(QTL)或基因[36)、玉米(37和小麦[38］．结合高通量测序技术，BSA已被用于鉴定木糖利用的新基因酿酒酵母．在此，我们利用BSA方法探索可能负责每叶腋花数的候选基因。这种表型是芝麻的一个重要农艺性状，因为它在预测产量中起着作用。29株芝麻菌株根据单花与三花(13株与16株)分为两组1:数据S1)。我们在两个基因库中鉴定出695个snp一致的基因。在这些基因中，有181个、21个和31个分别在cd中含有SNPs, 5'UTR和3'UTR1S18:数据美国)。GO术语注释主要与ATP结合、锌离子结合、核酸结合和热休克蛋白结合有关。特别感兴趣的是6个腺苷酸异戊烯基转移酶(ITP)同源物(SIN_1002735;SIN_1000260;SIN_1000476;SIN_1000477;SIN_1016115和SIN_1001679)，它们在玉米素生物合成途径中显著富集(图5)．玉米素是细胞分裂素植物激素家族的一员，已知其参与多种与植物生长发育相关的过程，包括促进侧芽生长和刺激细胞分裂以产生更浓密的植物[39]、[40］．本研究结果表明，ITP基因可能通过介导玉米素的生物合成而影响芝麻叶腋花数。然而，还需要利用转基因模型或双亲杂交群体进行进一步的研究。

下一代测序正迅速增进我们对作物遗传变异的了解[41］．本研究首次对高油料作物芝麻进行了全面的重测序分析。这些数据来自12个国家的29个品系，提供了对芝麻种质基因组遗传变异的洞察，并促进了广泛的功能和进化研究，包括基因组进化、群体遗传学、标记辅助育种和基因鉴定。在芝麻基因组中发现高LD，表明标记辅助育种是芝麻改良的较好选择。这里提供的数据为芝麻起源于印度次大陆的假说提供了新的证据。在编码和非编码区域，我们发现了数十万个多态性，这为进化遗传学和功能研究提供了重要的资源。特别令人感兴趣的是含有非同义突变的基因，包括可能介导与环境相互作用的大效应snp。本研究还表明，ITP基因可能在决定芝麻叶腋花数方面起一定作用。然而，还需要进一步的研究来充分了解本研究中确定的遗传变异的功能相关性。

选取29个栽培芝麻品系进行基因组重测序，其中16个来自中国，2个来自美国，1个来自阿富汗、阿拉伯联合酋长国、韩国、缅甸、菲律宾和越南。

图书馆建设与排序

用CTAB法从各品系的鲜叶和黄化叶中提取基因组DNA。根据制造商说明书(Illumina)，使用5 μg基因组DNA构建每个菌株的成对端测序文库，嵌入大小约500 bp。测序使用Illumina Hiseq 2000平台。然后，原始测序读取经过一系列严格的过滤步骤，根据以下标准删除读取:

类型(1):对于短插入库和长插入库，分别具有≥10%和≥3%的未识别核苷酸。

类型(2):40%的基础阅读有>，质量分数< 7。

类型(3):>的读取10bp对齐适配器序列，允许≤2-bp的不匹配。

类型(4):与对应配对端重叠≥10bp的配对端读数。

Type(5):两个完全相同的对端reads的Read1和read2(认为是PCR复制的产物)。

在所有过滤步骤之后，在深度≥13倍处生成了> 120 Gb(附加文件1:数据S1)。

SNP打电话

利用BWA软件将Reads映射到“中植13号”芝麻基因组组装图上[10］．使用的详细参数如下:

" bwa aln -m 200000 -o 1 -e 30 -i 15 -l 35-L -i -t 4 -n 0.04 -R 20 -f "

" bwa sampe -a 800 "

考虑到所有菌株都是一组，我们使用SAMtools函数“mpileup”[11使用映射质量≥20的reads检测原始种群snp。所用的参数如下:

" samtools mmpeup -uf -b -D | bcftools view -bvcgI -p 0.99 "

使用SAMtools程序“vcfutils”，使用上述过程提取的snp首先被过滤，得到30到581之间的测序深度。所用的参数如下:

“perl vcfutils.pl varFilter -d 30 -d 581”

根据以下标准进一步筛选原始SNP位点:拷贝数≤2，间隔至少5 bp，除了少数等位基因频率(MAF≥0.05)，当SNP之间的距离< 5 bp时，SNP会被保留。多样性参数π而且θ_w使用10 KB的窗口和1 KB的滑动窗口进行测量[12]、[14］．

为了检验SAMtools的SNP调用准确性，我们随机选取4个大小为4.5 ~ 8.1 kb的片段，用重叠引物进行扩增，得到的PCR产物进行Sanger测序。两种方法检测snp的符合率为92.3 ~ 95.2(平均为93.7%)2:图S4;额外的文件1:数据S3)。

此外，GATK工具包[42也被用来调用snp，如下所示:

我们首先使用BWA软件将干净的reads映射到芝麻基因组，参数如下:

" bwa aln -m 200000 -o 1 -e 30 -i 15 -l 35-L -i -t 4 -n 0.04 -R 20 -f "

" bwa sampe -a 800 "。

SAMtools用于分离、排序、rmdup和合并SAM对齐结果，picard-tools用于对bam结果进行排序并标记为副本。接下来，我们使用GATK程序从统一的基因型原始VCF中使用以下参数重新排列和过滤snp:

-T selectvariables -R -variant -concordance -o .jar

java -jar GenomeAnalysisTK.jar -T variantfiltering -R——filterExpression " QD < 20.0 || ReadPosRankSum < -8.0 || FS > 10.0 || QUAL < $MEANQUAL "——filterName LowQualFilter——missingValuesInExpressionsShouldEvaluatAsFailing——logging_level ERROR -o .jar

-T CombineVariants -R -V sample1. jar GenomeAnalysisTK.jarvcf - v sample2。vcf - gentypemergeoptions unique ify -o *。

利用GATK从16个连锁群中共获得2,003,821个群体snp。

短InDel检测

使用软件SOAPInDel的默认参数[43]，根据满足对端要求并包含对齐间隙的映射读提取至多5 bp的主短插入或删除，所有间隙都由至少3个非冗余的对端读支持。然后过滤主要InDel集，包括读取质量值> 20和InDels < 5 bp。

结构变异检测

根据配对端测序的原理，配对端读取必须有一个对正向序列，另一个对反向序列。两个对准位置之间的距离应与插入件尺寸一致。因此，与基因组对齐的两个对端reads应该具有正常的方向和适当的跨度。通过聚类分析异常配对端对齐，并与之前使用Breakdancer软件定义的结构变化类型进行比较[20.使用默认参数运行。所得的SV数据集包括INS(插入)、DEL(删除)、ITX(染色体内易位)、INV(倒置)和CTX(染色体间易位)，范围从10 bp到1 Mb。

连杆不平衡的计算

为了测量种群中的LD，我们计算了相关系数(r²)使用Haploview软件[32),如下:

1. (1）

   Ped和info文件作为输入文件生成。

2. （2）

   对于每一条染色体，例如LG1，参数设置为“java -jar haploview.jar -n -log LG1.log -pedfile LG1. genetype”。ped -info lg1 . gene - type.info -dprime -minGeno 0.6 -minMAF 0.01 -hwcutoff 0.001 -memory 2000 -maxdistance 500 "。

3. （3）

   然后用R脚本绘制曲线，绘制平均(R²)对一对明智的标记距离。

群体遗传学分析

多样性参数π而且θ_w使用10 KB的窗口和1 KB的滑动窗口进行测量[12]、[14］．顶部和底部的5%块基于π值，将这些区块中的基因分别定义为高发散度基因和低发散度基因(附加文件1:数据S7和S8)。

单个snp被用来计算样本之间的距离。在带引导(1000)的p-distance模型下，使用TreeBest (http://sourceforge.net/projects/treesoft/files/treebest/)对29个芝麻品系进行筛选。系统发育树显示软件MEGA5 [44］．使用软件EIGENSOFT [45］．软件FRAPPE [46]被用来确定种群结构。

BSA:

的隔离分析

cd:

编码序列

CTX:

染色体间易位

InDel:

插入和删除

走:

基因本体论

GWAS:

全基因组关联研究

INS:

插入

财年:

染色体内易位

发票:

反演

国际旅游业伙伴关系:

Isopentenyltransferase

LD:

连锁不平衡

主成分分析:

主成分分析

RDA:

农村发展管理

SNP:

单核苷酸多态性

SV:

结构变化

UTR:

翻译区

中国农业科学院油料作物研究所的研究得到国家作物种质资源基础设施项目(NB2012-2130135)、国家自然科学基金项目(no:31271766,没有。31101182也没有。国家作物种质资源保护(2012-014);中国农业科学研究系统(no. 31201242);汽车−15)。

从属关系

1. 中国农业科学院油料作物研究所，农业部油料作物生物学与遗传改良重点实验室，武汉430062

   王林海，张艳欣，李东华，魏鑫，丁霞，张秀荣

2. 北京华大基因(BGI)， -深圳，中国深圳

   雪莲汉

3. 1基因,杭州,中国

   雪莲汉

相应的作者

对应到Xiurong张．

相互竞争的利益

作者声明他们没有竞争利益。

作者的贡献

XRZ和LHW对研究的设计做出了贡献。LHW和XRZ撰写了手稿。XLH和LHW分别进行测序实验和数据分析。YXZ和DHL制备了材料并进行了实验。JYY执行数据库构建。所有作者阅读并批准了最终稿件。

开放获取本文由BioMed Central Ltd.授权发布。这是一篇开放获取文章，根据创作共用授权协议(https://creativecommons.org/licenses/by/2.0)，它允许在任何媒体上不受限制地使用、分发和复制，只要原始作品的名称正确。创作共用公共领域奉献放弃书(https://creativecommons.org/publicdomain/zero/1.0/)除另有说明外，适用于本条提供的资料。

再版和权限