农林育种研究

高通量技术在农林领域的广泛应用大大地加速了农林研究的进程，通过基因芯片及测序等技术可以帮助科研工作者解决农林领域研究中的各类问题，如：性状相关变异筛选、群体进化与选择、转录后基因调控、微生物群体结构多样性等等。

伯豪生物可以为您提供多组学的农林研究解决方案，包括：基因组层面全基因组重测序、Affymetrix SNP芯片、Illumina SNP芯片；转录组层面mRNA测序、全转录组测序/lncRNA测序、circRNA测序、miRNA测序；微生物层面16S rDNA/18S rDNA/ITS扩增子测序、宏基因组测序（元基因组测序）等。伯豪生物用科学的项目方案设计，严谨的实验数据质控，专业的生信分析团队，丰富的项目合作经验，为您的研究和数据带来扎实的服务，加速您的研究。

 

农林特色服务——BSA测序定位

一、BSA定位——快速定位性状基因

1.1 BSA方法简介

在动植物性状定位研究中，将家系后代中具有极端相反性状的个体分为两组，寻找出与性状呈共分离的分子标记，从而实现性状定位的方法，称为混合分组分析法（Bulked Segregant Analysis，BSA）。BSA的思想早在1991年就由Michelmore在莴苣的抗病性性状定位研究中提出1。在高通量测序普及的近十年内，BSA再次成为了性状定位研究的热点，NGM2、SHORE-Map3、QTL-seq4、MutMap5, 6等方法不断提出，并得到了广泛应用。

在实验设计上，BSA测序使用了极端性状的两个群体进行Pool-seq，对比变异位点基因频率的差异，找出那些具有明显差异的位点作为性状候选位点。根据应用的差异，BSA测序已产生了多种分析方法。QTL-seq可用于应对自然发生的数量性状和质量性状。MutMap是一种特定的应用，分析寻找的是EMS诱变导致的突变性状，而非自然产生的突变。下面将详细介绍这两种方法。

 

 1.2 QTL-seq

QTL-seq是Hiroki（2013）提出的QTL定位方法4。与传统QTL连锁定位方法一样，QTL-seq也是基于遗传群体来进行实验的，但它更为快速和节省成本，并可达到较高的准确性。下面是取样设计和分析原理的示意图（图1.2）。

图1.2. QTL-seq流程示意图（Hiroki 2013）4

 

1.3 MutMap

MutMap是用于寻找EMS诱变性状的方法。在取样设计和分析原理上与QTL-seq十分类似，同样是基于遗传群体的取样设计，在混池测序后比较基因频率，找出SNP-index接近1的候选位点。下面是原理示意图（图1.4）。

图1.4. MutMap流程示意图（Abe 2012）6

 提出QTL-seq分析方法的Ryohei Terauchi教授研究组还针对分析中可能遇到的各种情况提出了解决方法（表1.1）。MutMap+7和MutMap-Gap8是MutMap的特殊形式，适用于突变致死、无法异交和性状位点不在参考基因组的情况。

表1.1. MutMap的衍生应用

 

 二、技术路线

 

 

三、BSA测序实验设计

3.1 实验参数

QTL-seq和MutMap除了实验材料的准备上略有不同，其它几乎没有区别（表2.1）。

表2.1. BSA测序实验参数

 

3.2 注意事项

取样策略：对于QTL-seq取样，有时需要在极端性状和样本数量之间做均衡。理想情况下，取极端各5%的群体30个样本，则需群体大小为 30/0.05 = 600个。在群体数量不足时（尤其是RILs群体），则可放宽极端性状标准，并适当减少混池样本数目。

误差来源：由于使用了Pool-seq的测序策略，在估计基因频率时很难排除由于测序错误导致的偏差，需要在分析中使用sliding windows的策略来降噪；样本DNA质量不均、混样浓度不统一将会为频率估计带来偏差，但这可由增加群体样本数量来进行一定程度的弥补。

 

3.3 BSA测序应用案例

案例一 MutMap方法快速培育耐盐水稻品系Kaijin5

在日本的Ryohei Terauchi团队提出MutMap方法两年后，他们成功使用该方法在两年之内快速培育出了一个耐盐水稻品系Kaijin并成功推广（图2.1.）。

图2.1. a. SNP-index图，蓝色点为每个SNP的SNP-index数值，红色线为以窗口（Sliding windows= 4Mb，步移10kb）为单位的SNP-index分析，橙色线为99%置信区间；b. OsRR22/HST1基因结构，绿色框为Tos17插入突变位置；c. Sanger测序验证突变位点导致提前终止；d. Tos17插入突变验证OsRR22的功能，不耐和耐盐分离比为3:1。

 

研究使用了6,000株EMS诱变过的水稻品系Hitomebore，经过7天的1.5% NaCl条件培养，选择出耐盐突变品系hst1。经过与已有耐盐品系比较，发现hst1具有更强的耐盐性，适合作为优良品系育种的亲本。接下来，hst1与野生型Hitomebore进行杂交并构建F2群体，经过0.75% NaCl处理，发现不耐盐和耐盐分离比符合3:1的关系。将耐盐F2子代20株进行混池测序，并与野生型亲本进行比较，得到1,005个突变位点。对这些突变位点进行SNP-index的计算，统计出显著偏离0.5的点，发现2个SNP-index =1的位点，其中一个导致OsRR22基因的第三外显子发生提前终止，另一个位于非编码区。作者使用了水稻Tos17突变体库的OsRR22突变株验证了该突变结果确实将导致水稻耐盐，并通过回交亲本和Sanger测序验证培育出耐盐品系Kaijin。此项工作证明了MutMap方法在育种领域的高效性。

 

原文出处：Takagi H, Tamiru M, Abe A, et al. MutMap accelerates breeding of a salt-tolerant rice cultivar[J]. Nature Biotechnology, 2015, 33(5):445-9.

 

案例二 使用QTL-seq快速定位水稻的数量性状4

研究者使用了抗稻瘟病和易感稻瘟病品系杂交构建了RILs群体，并对直链淀粉低含量和高含量品系、稻苗低活力和高活力品系杂交构建了F2群体。随后，这些性状在RILs和F2群体中出现了性状分离，作者挑选了具有极端表型（极抗病/极易感病，极高/极低淀粉含量，以及极高/极低活力）的个体分别混为两个DNA池并进行建库，并进行全基因组重测序，通过群体间SNP频率的差异来进行基因定位。结果表明，QTL-seq与传统QTL连锁定位方法得到的QTLs位点相符（图2.2）。

图2.2.在水稻RILs群体中用QTL-seq寻找抗稻瘟病位点。a. 抗稻瘟病品系Nortai和易感病品系Hitomebore在2周真菌处理后的情况；b. 对RILs不同个体抗病表型进行分级；c. 4和5级作为抗病组，8 – 10 级作为易感组进行混池；d. 绿色和黄色分别代表抗病和易感组，蓝色代表两者SNP频率差值；底部图为QTL连锁分析结果。

 

原文出处：Takagi H, Abe A, Yoshida K, et al. QTL‐seq: rapid mapping of quantitative trait loci in rice by whole genome resequencing of DNA from two bulked populations [J]. Plant Journal for Cell & Molecular Biology, 2013, 74(1):174–183.

 

参考文献

  1. Michelmore, R.W., Paran, I. & Kesseli, R.V. Identification of markers linked to disease-resistance genes by bulked segregant analysis: a rapid method to detect markers in specific genomic regions by using segregating populations. Proceedings of the National Academy of Sciences88, 9828-9832 (1991).
  2. Austin, R.S. et al. Next-generation mapping of Arabidopsis genes. Plant Journal for Cell & Molecular Biology67, 715 (2011).
  3. Sun, H. & Schneeberger, K. SHOREmap v3.0: Fast and Accurate Identification of Causal Mutations from Forward Genetic Screens. Methods in Molecular Biology1284, pp 381-395 (2015).
  4. Takagi, H. et al. QTL-seq: rapid mapping of quantitative trait loci in rice by whole genome resequencing of DNA from two bulked populations. Plant Journal74, 174-183 (2013).
  5. Takagi, H. et al. MutMap accelerates breeding of a salt-tolerant rice cultivar. Nature Biotechnology33, 445 (2015).
  6. Abe, A. et al. Genome sequencing reveals agronomically important loci in rice using MutMap. Nature Biotechnology30, 174-178 (2012).
  7. Fekih, R. et al. MutMap+: genetic mapping and mutant identification without crossing in rice. 8, e68529 (2013).
  8. Takagi, H. et al. MutMap-Gap: whole-genome resequencing of mutant F2 progeny bulk combined with de novo assembly of gap regions identifies the rice blast resistance gene Pii. New Phytologist200, 276 (2013).

 

 

农林特色服务——全基因组关联分析

一、全基因组关联分析概述

全基因组关联分析（genome-wide association studies，GWAS）是一种在全基因组水平上将大量样本的表型和基因型进行对照分析或相关性分析，从而对复杂性状进行基因定位的方法。目前，GWAS已经成为动植物数量性状位点（quantitative trait locus，QTL）精细定位的重要手段。随着越来越多的物种全基因组测序完成，GWAS将更广泛地应用于育种和进化领域。

在动植物QTL定位的研究设计中，根据样本材料来源不同可分为杂交家系材料和自然群体材料。在家系材料中，主要使用连锁分析方法考察重组事件带来的连锁。由于家系传代次数和样本数受限制，可供观察的重组事件有限，因此定位的结果通常位于一段较长的染色体区间内。而在自然群体中，使用的是关联分析方法来考察群体中存在的连锁不平衡。相对于杂交家系有限的传代次数来说，自然群体中蕴藏着丰富的重组事件。因此，关联分析定位的精度要远高于连锁分析。但是，由于不同地理群体间存在着的遗传分化也会带来连锁不平衡，关联分析的假阳性将高于传统的连锁分析。

 

二. 技术路线

 

三. 实验参数

实验方法：全基因组重测序 / SNP基因芯片

样本数量：最少 > 200个，推荐 > 500个

测序深度： > 5X（样本数 < 500）；大于500个样本可适当降低测序深度

测序平台：Illumina Hiseq

 

四. 结果展示

 (a) 数据质控

对SNP位点的质控（例）：去除MAF < 0.01的数据，去除位点call-rate < 0.9的数据，并去除显著偏离哈温平衡的位点（HW-p < 10-6）。

对样本的质控（例）：个体call-rate < 0.9、个体常染色体异质性（FDR < 0.01）、剔除亲缘关系较近的个体（IBS > 0.95）。

 

(b) 数量性状正态分布检验

对于质量性状，采用case-control分析（性状只有0和1）；对于数量性状，数据需要验证是否符合正态分布，并将不符合正态分布的数据转换为正态分布的结果。样本量如小于2000，结果以Shapiro-Wilk检验为准，Sig<0.05则需要进行数据转换。

表3. 数量性状正态性检验（示例）

 

对Sig. < 0.05的Trait4和Trait7进行正态转换，得到以下结果：

表4. 数量性状转换后正态性检验（示例）

 

(c). 关联分析

使用线性回归模型进行全基因组关联分析，公式为y = Xα+ Zβ + Wμ+ e，其中y为表型性状；Xα为固定效应，即影响y的因素如群体分层、性别年龄等可预测因素；Zβ为标记效应；Wμ为随机效应，一般指个体亲缘关系。

对所得p值进行多重检验校正，以校正后的p值计算Q-Q plot和生成曼哈顿图。

 

(d). Q – Qplot

Q-Q plot 即Quantile-Quantile Plot，主要是直观的表示观测值与期望值之间的差异，可直观的看出显著位点情况。位点值为经过校正后的p值。理论上每个点在Q-Q plot中应该是观测值与预测值相等，这时可得到一条斜率为1的直线。如偏离直线则说明此SNP点出现了偏离。

图11. Q-Qplot图（示例）

 

(e). 曼哈顿图

使用R软件包进行曼哈顿图的绘制。横坐标为各SNP位点在染色体上的位置，纵坐标为p值的负对数值。设置阈值，筛选显著相关的位点作为候选位点。

图12. GWAS分析结果（示例）

 

(f). 显著位点功能注释

对以上筛选得到的显著位点进行功能注释，得到该位点对应的基因以及该基因的功能。

 

五. 研究案例

案例四 借助Affymetrix Axiom平台使用GWAS研究鸡攻击行为相关基因7

攻击行为是很多物种共有的进化属性，它可以帮助雄性物种快速建立生存空间或社会主导地位。然而可想而知，这种属性在禽类的饲养中是十分不利的。之前已有一些研究表明鸡的攻击行为可由遗传因素导致，并估计了遗传力大小约为0.57，说明群体中一半以上的攻击行为是由遗传导致的。作者试图使用GWAS方法来探究攻击行为背后的遗传因素，以期改善养殖效率和动物福利。

研究共使用了265只雄性家鸡，在出生第60天时开始测量攻击性状，共4个，分别为T1：在记录期间（16天）内打斗的次数；T2：每天打斗的次数（超过4次则被记录）；T3：打斗的天数；T4：超过4次打斗的天数。基因分型工作由上海伯豪完成，使用了Affymetrix Axiom平台的高密度芯片，型号为600K Affymetrix Axiom HD chicken genotyping array，共生成了599,898个SNP位点。经过质控后，保留468,020个SNP用于GWAS分析。结果共得到33个显著位点，涉及26个基因。

图13. T1-T4四个攻击行为性状的Manhattan图，蓝和红色线分别表示P=4.6E-6和P=2.3E-7

 

作者使用了Ingenuity Pathway Analysis（IPA）对26个相关基因进行代谢通路分析，根据基因互作网络进一步得到9个基因，其中SORCS2基因与多个基因互作，可能具有重要功能。在80个显著代谢通路中，“行为”相关通路排第26位。接下来，作者选择了SORCS2基因进行后续验证。

图14. Ingenuity Pathway Analysis（IPA）得到的9个基因。其中，SORCS2基因与多个基因互作，可能具有重要功能

 

在IPA分析中，SORCS2基因与NGF, NGFR, L-dopa and dopamine四个基因相关，而NGF已在小鼠中被证明其表达量与攻击行为相关。于是，研究在鸡成纤维细胞系DF-1中敲降了SORCS2基因的表达量，观察到了NGF, NGFR, L-dopa基因的表达量也随之降低，其中NGF, L-dopa和多巴胺受体基因家族DRD1–DRD4显著下降，但NGFR下降量并不显著。此结果显示SORCS2基因可能会影响多巴胺通路的表达，从而影响到鸡的攻击行为程度。

 

图15. 左图：NGF, NGFR, L-dopa基因由于SORCS2基因敲降而降低；右图：多巴胺受体基因家族DRD1–DRD4由于SORCS2基因敲降而降低

 

文章对比了攻击行为程度高的和低的鸡个体的SORCS2基因表达量，发现在攻击行为程度高的鸡个体中，SORCS2基因表达量显著高于攻击行为程度低的鸡个体。这个结果进一步说明了SORCS2基因可能在鸡攻击行为中扮演了重要的角色。

图16. 攻击行为程度高的鸡个体中SORCS2基因表达量显著升高

 

研究使用了265只雄性家鸡，通过Affymetrix Axiom基因分型平台进行SNP分型，得到大量高密度SNP位点进行GWAS分析，得到33个候选位点。通过Ingenuity Pathway Analysis（IPA）代谢通路分析，将研究对象聚焦到少数几个基因上，并通过使用DF-1细胞系进行表达量验证，结果说明了SORCS2基因与鸡攻击行为之间的关系。

 

原文出处：Li ZH, Zheng M, Abdalla A, et al. Genome-wide association study of aggressive behaviour in chicken. Scientific Reports, 2016.