González-Martínez SC, Ersoz E, Brown GR, et al (2006) DNA Sequence Variation and Selection of Tag Single-Nucleotide Polymorphisms at Candidate Genes for Drought-Stress Response in Pinus taeda L. Genetics 172:1915–1926. doi: 10.1534/genetics.105.047126

摘要

遗传关联研究正迅速成为解剖复杂性状的首选实验方法，包括对干旱胁迫的耐受性，干旱胁迫是森林树木死亡和产量损失的最常见原因。关联作图的优化需要知道核苷酸多样性和连锁不平衡的模式以及用于基因分型的合适多态性的选择。此外，应用于DNA序列变异数据的标准中性测试可用于选择推定选择用于关联作图的候选基因或氨基酸位点。在本文中，我们研究了18种候选基因的多态性模式，这些基因是一种重要的树木作物 - 火炬松（P. taeda L.）的干旱胁迫响应。基于一组21个推定中性核微卫星的数据分析未显示群体遗传结构或全基因组偏离中性【为什么不直接用基因序列数据呢？】。候选基因在沉默位点具有中等的平均核苷酸多样性（=0.00853），在单个基因中变化100倍。基因内LD水平较低，平均成对r2为0.30，在800bp时从0.50迅速衰减至0.20。基因间未发现明显的LD。虽然群体扩张也可以解释我们的结果并且选择的证据不是决定性的，但是在早期响应干旱-3（erd3）基因中可能已经发生了选择性清除。另一个基因，ccoaomt-1，一种参与木质化的甲基化酶，显示出二态性（即两个频率相等的高度不同的单倍型谱系），这通常与平衡选择的长期作用有关。最后，选择一组单倍型标记SNP（htSNP）。使用htSNPs，基因分型工作量减少30％-40％，同时采样最常见的等位基因变异，可以在我们正在进行的松树抗旱性相关研究中获得。

前言

分子进化的中性理论指出核苷酸多样性由种群突变参数控制，其中是每代、每点的突变率。在过去的二十年中，在自然群体中选择下的候选基因的鉴定依赖于物种内和物种之间的核苷酸多样性模式的分析以及来自中性期望的等位基因（单倍型）分布的偏离（即，中性测试；见综述文章Kreitman 2000; Ford 2002; Rosenberg和Nordborg 2002）。这些分析在广泛的基因和生物体中出现了两种主要模式。一种类型的基因座显示出过多的中频单倍型，通常排列在两个高度不同的谱系周围（例如，Filatov和Charlesworth 1999; Tian等2002），另一种特征是过量的稀有单倍型（例如，Olsen等 2002;参见松树研究Pot et al 2005）。只要不是由于群体统计学或群体结构，这些的频谱与中性期望的偏差分别与平衡选择和由正选择引起的净化选择或选择性清除相关联。
等位基因变异与性状之间的遗传关联在种群尺度上的差异是鉴定导致适应性状变异的基因或等位基因的强大且相对较新的方法（Long和Langley 1999;参见针对针叶树的Neale和Savolainen 2004）。群体分层是关联研究中最常见的系统偏差来源（Buckler和Thornsberry，2002; Hirschhorn和Daly，2005）。假设中性分子标记，例如核微卫星，通常用于检测群体结构和其他可能在关联研究中产生假阳性的群体和群体统计过程（Rosenberg等2002）。关联作图的优化需要知道每个特定物种和候选基因组的核苷酸多样性模式和连锁不平衡。另外，应用于单个或几个基因的DNA序列的标准中性测试可用于选择候选基因或氨基酸位点，其被推定选择用于关联作图。
森林树木在陆地生态系统中发挥着至关重要的作用，在气候控制，碳固定和野生动物维护方面提供了重要的生态效益。干旱胁迫是树木死亡的最常见原因，并且是造成商业物种严重年度产量损失的原因（Pinus taeda L。中的高达65％; Burns和Honkala 1990）。了解干旱胁迫耐受性的生理机制和遗传基础一直是植物生物学家长期关注的问题（例如，Ingram和Bartels 1996; Seki等2003;参见Newton等1991的森林树木）。然而，在鉴定干旱相关基因和开发森林树木表达研究方面的进展是相对较新的（Chang等1996; Dubos和Plomion 2003; Watkinson等2003）。树木中脱水耐受性的分子基础极其复杂，并且已经提出了多种表达候选基因。在种子发育（Jarvis等1996）和干旱胁迫（Richard等2000; Watkinson等2003）期间，在不同的针叶树中发现了脱水蛋白的表达增加。 Chang等（1996）使用消减杂交方法鉴定了四个在P.taeda中具有干旱诱导表达的cDNA克隆：lp2，与乙烯合成中的中间体S-腺苷甲硫氨酸合成酶（sams）具有高度同源性; lp3，主要在根中表达，后来发现属于ABA诱导基因的小家族（Padmanabhan等1997）; lp4，类似于I型含铜糖蛋白;lp5，几乎完全在根中表达并编码与细胞壁蛋白类似的富含甘氨酸的蛋白质。在树中鉴定的干旱 - 胁迫反应的其他主要表达候选基因编码蛋白激酶（Dubos和Plomion 2003; Dubos等2003），半胱氨酸蛋白酶（Tranbarger和Misra 1996），铁贮存蛋白（Li等1998），抗氧化剂（Li等1998; Karpinska等2001）和发病相关蛋白（Dubos和Plomion 2001; Dubos等2003）。
针叶树是长寿的，广泛分布的生物，通常表现出高水平的杂合性和大的有效种群大小。因此，有人提出针叶树可能表现出高水平的核苷酸变异（Dvornyk等2002）。然而，针对针叶树的DNA序列变异的第一个结果最多显示出对核苷酸多样性的中等估计（例如，Kado等人2003; Brown等人2004; Pot等人2005）。针叶树的平均种群分化也是中等的（Kado等2003;但参见Pot等2005的korrigan和pp1基因），即使采样极端表型（Garc磇a-Gil等2003）。例如，Garc磇a-Gil等（2003）未发现沿纬度克隆采样的两个光敏色素位点的光敏结构域的任何功能分化，这与生长物候学的显着差异相关（如常见的园林实验所示）。核苷酸多样性和/或群体分化的模式偏离中性期望，可能表明自然选择的作用，仅描述了少数基因和树种[日本柳杉cl5（ Kado等2003）; Pseudotsuga menziesii中的f3h1、4cl1和mt-like（Krutovsky和Neale 2005）;和松树中的pp1、korrigan和CesA3（Pot等2005）]。由于群体水平的高重组率，针叶树中的大的有效群体大小将导致低连锁不平衡（LD）。该预测与针叶树中的经验数据一致，其中观察到基因中缺乏LD和基因内的LD相对快速衰减（200-500bp）（Brown等2004; Rafalski和Morgante 2004）。然而，在欧洲和美国，在更新世冰川避难所经历过历史瓶颈的特定树种或种群中，是否存在更广泛的LD是可能的，但目前尚不清楚。
由于遗传漂移、群体统计学、种群结构和自然选择之间的相互作用，自然种群的常规变异被模式化。在本文中，我们使用21个核微卫星的数据集来检测群体结构和群体统计过程，这些过程可能在关联研究和偏倚中性测试中引起虚假关联，并对18个候选基因的全部或部分进行测序，用于一种重要的树木作物P. taeda的干旱 - 应激反应。我们的样本覆盖了东南部本地的P. taeda，包括佛罗里达，这个物种的推定的更新世冰川避难所（Schmidtling等1999; Al-Rabab'ah和Williams 2002），在我们之前的研究中没有广泛采样（见Brown等，2004）。我们使用DNA序列来估计核苷酸多样性和连锁不平衡的水平，以确定选择的候选基因（通过中性测试），并选择单倍型标记单核苷酸多态性（htSNP）用于我们当前的遗传关联研究。

材料和方法

植物材料：

使用P.taeda的32个种子巨型配子体【胚乳】（针叶树种子的单倍体，母系衍生营养组织）的样品（来自30棵树中的1颗和来自1棵树中的2颗）用于SNP开发。种子供体包括22个不相关的第一代选择（优树）来自未受干扰的自然林分，覆盖东南部的P. taeda[大西洋沿海平原（ACP），佛罗里达州中部，佛罗里达州北部，马里恩县和墨西哥湾沿岸种源;参见http://www.genetics.org/supplemental/上的补充表S1]和大西洋沿岸平原种源中第一代选择中受控杂交产生的9个第二代选择。这些树目前是森林生物学研究合作社（FBRC）树木改良计划（佛罗里达大学，佛罗里达州盖恩斯维尔）的一部分。由于包含四对半同胞和三棵树作为亲本的第一代选择，第二代树可能引入轻微偏差（参见补充表S1）。然而，由于该物种中的杂合性水平较高且减数分裂偏析，因此偏倚被认为可以忽略不计。

候选基因选择：

根据（1）公共数据库中的P.taeda表达序列标签（ESTs）的重叠群组成（DDBJ / EMBL / GenBank）与模型物种中的干旱胁迫响应基因的同源性选择干旱胁迫响应的候选基因; （2）基于6个木质部EST文库（通过http：// pinetree.ccgb.umn.edu/访问）与干旱胁迫响应在北卡罗来纳州立大学组装的单基因组（20,500个非冗余基因）序列的同源性模型物种中的基因; （3）干旱胁迫下来自P. taeda树的根库中EST的过多，与使用MAGIC Gene Discovery工具的“电子”Northerns所示的对照树相比（佐治亚大学，http：// fungen.org） /Projects/Pine/Pine.htm）。还选择了另外两个基因ppap12和lp3-3，因为它们在干旱处理中显示差异表达，如分别在P. pinaster（Dubos等2003）和P.taeda（Padmanabhan等1997）中的反向Northerns所示。

DNA分离，扩增和测序：

在种子萌发后使用Plant DNeasy试剂盒（QIAGEN，Valencia，CA）从巨型配子体中提取单倍体基因组DNA。设计PCR引物以扩增9个核基因座中的400至1000bp片段，并且先前公开的引物用于另外9个基因（参见补充表S2，在http://www.genetics.org/补充/）。设计引物以扩增lp3-3，dhn-1和dhn-2的全长基因。使用BigDye Terminator v.3.1循环测序试剂盒（Applied Biosystems，Foster City，CA），在ABI 377自动测序仪上直接从PCR产物获得序列数据。所有样品均从两端测序。使用phred和phrap程序（Ewing等1998; Gordon等1998; http：// bozeman.mbt.washington.edu/phredphrapconsed.html）在Unix环境下进行碱基调用和正向和反向读取的组装。来自基因座的多个等位基因在多重比对扩展（MACE）程序（B.Gilliland和C.Langley，University of California，Davis，CA）中比对。肉眼检查所有色谱图，并且只有当所有序列的phred得分在该位点超过25时才接受推定的序列变体。根据需要进行重新测序以维持该质量标准。由于DNA样品是单倍体，单倍型（即等位基因）的鉴定是明确的。

候选基因的定位：

18个候选基因中的6个先前已被定位（Brown等，2003）。使用P.taeda的两个参考作图群体，qtl和基础谱系（Brown等人2001中的细节）尝试绘制剩余的12个基因座。根据Temesgen等人（2001）使用变性梯度凝胶（DGGE）绘制五个候选基因（lp3-1，dhn-1，rd21Alike，cpk3和ppap12），并使用模板指导绘制1（lp3-3）。具有荧光偏振（FP）检测的染料终止掺入测定（TDI）（TDI 59→9引物：TTGCCAGTAGCATACACA TCTG）。使用AcycloPrime-FP SNP检测试剂盒和Wallac VICTOR2荧光板读数器（Perkin-Elmer Life and Analytical Sciences，Torrance，CA）进行FP + DI。
其他6个候选基因未连接（sod-chl）或缺乏合适的多态性（即，由于SNP附近存在重复区域，因此无法为任何SNP分离谱系的亲本或FP的DI引物无法设计：铁蛋白，erd-3，dhn-2，lp5-like和ug-2_498）。在Brown等人（2001）之后获得了共识图以及其他标记。

群体结构和群体统计过程：

群体分层是最常见的系统偏差，在关联研究中产生假阳性关联（Marchini等2004; Hirschhorn和Daly，2005）。此外，群体遗传结构或群体统计过程的存在，例如范围扩展或撤退，可能会产生类似于自然选择行为产生的等位基因频谱的签名，并误导对标准中性测试的解释，例如Tajima's D。我们使用21个高度多态性（平均每个基因座15个等位基因）核微卫星（nuSSRs），覆盖大多数P. taeda连锁群，以测试种群结构或群体统计过程。 nuSSR数据由C. Dana Nelson（美国农业部南方森林遗传研究所）友情提供，其中包括94棵树，其采样范围与此处提供的序列数据大致相同（参见补充表S3）。
为了测试种群结构，我们首先使用基于模型的聚类算法（Structure软件; Pritchard等2000; Rosenberg等2002），其构建没有任何先前地理信息的群体组。具有推定数量的簇（K参数）的模型，从1到4，不相关的等位基因频率，以及两个burn-in，以最小化起始配置的影响，并且运行长度为10^6。其次，我们计算了样本中包括的三个地理区域（墨西哥湾沿岸，东北和东南部）的遗传分化估计（F-统计，基于Weir和Cockerham 1984之后的嵌套ANOVA）。使用置换测试（10,000个排列）和基因座上的jiackkifed估计来测试区域之间群体遗传结构的显着性。
为了测试全基因组偏离中性，例如由群体统计过程产生的那些，Ewens-Watterson中立性检验（Watterson 1978,1986），基于纯合性和Fisher精确检验计算的概率（Ewens璚atterson璖latkin's exact; Slatkin 1994,1996），使用Arlequin v.2000程序（Schneider等2000）进行。考虑到在基因座处发现的等位基因的数量，Ewens璚atterson测试使得能够检测到与中性模型的偏差，作为相对于中性平衡期望的缺陷或过量的纯合性。应该注意的是，纯合性过量是种群扩增的典型全基因组特征（Payseur等2002; Luikart等2003）。一旦计算了21个nuSSR基因座中的每一个的测试，使用Mann-Whitney U-检验来检测预期的和观察到的纯合性值是否来自相同的分布。必要时应用Bonferroni校正多次测试。

核苷酸变异和中性测试：

使用DnaSP v.4.0（Rozas等人，2003）进行序列数据的分析。通过Watterson的（Watterson 1975）和估计核苷酸多样性，样品中序列之间的成对核苷酸差异的平均数量（Nei和Li 1979）。使用没有重组的聚结模拟评估基因座上序列变异的异质性。进行了许多统计分析以鉴定偏离标准中性进化模型的基因或氨基酸位点。针对完整序列和滑动窗口（窗口长度和步长分别为100和25个站点）计算每个基因座的田岛（1989）D-统计量。田岛的D统计量反映了和之间的差异。在突变 - 漂移平衡时，D的预期值接近于零。还计算了基于以值（由估计）为条件的单倍型（基因）频率分布的中性Fs检验统计量（Fu 1997）。 Tajima的D和Fu的Fs检验统计数据也可以反映群体变化（Fu 1997; Sano和Tachida 2005）。为了计算需要来自外群的数据的测试，14个基因的推定直向同源物【外群】来自P. pinaster，这是一种欧洲物种，可能与1.2亿年前的P. taeda不同（Krupkinet al.1996）。对于8个基因，我们使用来自GenBank的序列（登录号：AL751338，lp3-1; BX255067，dhn-1; BX677401，lp5样; BX252032，sod-chl; BX681838，sams-2; AY641535，pal-1; CR393126，ccoaomt-1;和AJ309112，ppap12），对于其他6个，我们使用直接从P. pinaster megagametophyte DNA获得的序列，使用相同的引物对进行测序，如P. taeda（基因dhn-2，rd21A-like） ，pp2c，Aqua-MIP，erd-3和ug-2_498; A. Soto和MT Cervera，未发表的数据）。然后，我们计算：（1）Fay和Wu的H检验（Fay和Wu 2000），基于选择性清除后立即预期的高频衍生等位基因的相对过量; （2）Hudson-Kreitman-Aguadé（HKA）检验（Hudsonet al.1987），检验特定区域的多态性和分歧之间的解耦; （3）McDonald-Kreitman（MK）测试（McDonald和Kreitman 1991），基于物种内和物种之间的同义和非同义替换的比较。进行了HKA测试，比较每个基因与其他每个基因的比较。最后，为了检测单个氨基酸位点的阳性选择，我们使用在线DataMonkey软件包中实现的基于似然的方法估计序列比对中每个位点的非同义和同义变化率（Kosakovsky-Pond和Frost 2005a， b）。对于这些分析，我们使用保守的单可能性祖先计数（SLAC）方法，与Suzuki-Gojobori（Suzuki和Gojobori 1999）相关，以及固定效应可能性（FEL）方法，它直接估计非同义和每个站点的同义替换率更适合具有中等数量序列的数据集（n = 20-40; Kosakovsky-Pond和Frost 2005a）。

LD，单倍型多样性和用于关联作图的htSNP的选择：

使用Tassel软件（http://www.maizegenetics.net/index.php?page=）计算LD描述性统计量r2（Hill和Robertson 1968），仅基于信息位点（频率=2/32= 0.063）。r2统计量总结了重组和突变历史，并且它对样本量的敏感性低于其他常见的LD统计量，如D'（Flint-Garcíaetal.2003）。 r2的统计显着性用单尾Fisher精确检验计算，并应用Bonferroni校正进行多次检验。利用r2估计的多态性位点之间的LD的非线性回归和位点之间的碱基对的距离来估计连锁不平衡与物理距离的衰减（Remingtonet al.2001; Ingvarsson 2005）。为了调整非线性函数，我们使用Hill和Weir（1988）提供的r2期望，用于具有低水平突变和样本大小n的调整的漂移 - 重组平衡，

（1）
其中C是总体重组参数。使用在SAS v.8.0统计软件包（SAS Institute，Cary，NC）的proc nlin中实现的Gauss-Newton算法拟合方程1。在Nei（1987）之后计算了单倍型多样性（He）。我们使用HaploblockFinder软件（Zhang和Jin 2003; http://cgi.uc.edu/cgi-bin/kzhang/haploBlockFinder.cgi/）在线确定htSNP，即代表常见等位基因变体的htSNP，阈值为r2 = 0.2来定义LD块。随着等位基因频率的降低，关联研究（对于固定样本量）的能力显着降低（Wanget al.2005）。然后，考虑到仅对应于常见（MAF> 5％）和频繁（MAF> 15％）SNP的次要等位基因频率（MAF）来选择htSNP。鉴于松树中发现的LD水平较低，导致LD区块较短，其他识别htSNP的方法，如LD块内LD子群的识别（详见Takeuchiet al.2005），效果不佳，就没展示。