CN113151542A

CN113151542A - 一种华山松基因组snp的开发方法和应用

Info

Publication number: CN113151542A
Application number: CN202110290767.2A
Authority: CN
Inventors: 辛培尧; 王正德; 叶鹏; 刘成; 王飞; 辛静; 董章宏; 瞿绍宏
Original assignee: Southwest Forestry University
Current assignee: Southwest Forestry University
Priority date: 2021-03-18
Filing date: 2021-03-18
Publication date: 2021-07-23
Anticipated expiration: 2041-03-18
Also published as: CN113151542B

Abstract

本发明涉及基因领域，具体涉及一种华山松基因组SNP的开发方法及应用，以SLAF‑seq技术为核心，对209份华山松优良单株进行多态性SNP标记开发，共获得了3469074个SNP标记；经GWAS关联分析后，获得了与胸径相关的SNP位点6个，与树高相关的SNP位点7个，利用SNP的序列在NCBI数据库进行Blast比对之后，发现2个基因可能与胸径生长相关，分别为SNP Marker227806位点相近的与纤维素合成相关的CesA2基因，以及Marker227806位点相近的与木质素合成相关的CCoAOMT基因，为华山松育种性状的遗传机理及分子标记选择奠定了理论基础。

Description

一种华山松基因组SNP的开发方法和应用

技术领域

本发明涉及基因领域，具体涉及一种华山松基因组SNP的开发方法和应用。

背景技术

华山松(Pinus armandii Franch.)系松科(Pinaceae)，松属(Pinus)单维管束亚属的高大乔木植物，是我国特有的五针松高大乔木、同时也是高海拔地区重要的造林树种之一。由于其生长周期较长、基因组高度杂合以及多种因素的影响，使得华山松在基因组学分子标记辅助选择方面的研究相对滞后。

发明内容

为解决上述问题，本发明提供了一种华山松基因组SNP的开发方法和应用。

为实现上述目的，本发明采取的技术方案为：

一种华山松基因组SNP的开发方法，包括如下步骤：

S1、以 SLAF-seq技术为核心，对209份华山松优良单株进行多态性 SNP标记开发，用EcoRV-HF®+ScaI-HF®酶来构建简化基因组测序文库，最后在 IlluminaHiSeq2500system上测序，得到SLAF标签，用软件BWA把测序得到的reads进行比对，以次等位基因频率（MAF）和位点完整度为依据，过滤掉MAF<5%和完整度小于80%的低质量SNPs，而后结合GATK和SAMTOOLS软件作为开发 SNP标记的方法，共获得了 3469074 个 SNP标记；

S2、利用开发的华山松基因组 SNP，对表型树高（Height）和胸径（DBH）性状进行GWAS分析，采用GLM模型获得了与胸径相关的 SNP位点 6 个，与树高相关的SNP位点7个，利用SNP的序列在NCBI数据库进行 Blast比对之后，发现2 个基因可能与胸径生长相关，分别为SNPMarker227806 位点相近的与纤维素合成相关的CesA2基因（登录号：AC241331.1），以及Marker227806 位点相近的与木质素合成相关的CCoAOMT基因（登录号：AC241331.1）。

所述 SNPMarker227806可用于挑选华山松优良单株。

本发明首次开发了华山松基因组 SNP 3469074 个，获得与胸径、树高性状相关联的 SNP 标记 13 个，2 个候选基因CesA2和CCoAOMT，为华山松育种性状的遗传机理及分子标记选择奠定了理论基础。

附图说明

图1为本发明的技术路线图。

图2为信息分析流程图。

图 3为DNA 提取部分结果图。

图 4为华山松测序质量分布图。

图 5为华山松碱基分布图。

图 6为树高性状的全基因组关联分析 QQ 图。

图 7为RNA 提取结果图。

图 8为CesA2基因在各个组织的表达量结果图。

图 9为CCoAOMT基因在各个组织的表达量结果图。

具体实施方式

下面结合具体实施例对本发明进行详细说明。以下实施例将有助于本领域的技术人员进一步理解本发明，但不以任何形式限制本发明。应当指出的是，对本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进。这些都属于本发明的保护范围。

1.1材料来源

材料选自云南省楚雄市紫金山林场华山松无性系种子园，具体见表1。

表1 209份华山松试验材料来源。

1.1.1 华山松树高与胸径的测定

在云南省楚雄市紫金山林场种子园内，进行同年种植的 6 个华山松无性系种源实地调查，利用测高仪测树高，测树围尺测量胸径，共测量了 1453 个单株。

1.1.2 材料

以前期树高、胸径基础数据及实测数据为依据，按照从大到小进行排列，最终筛选出 30 株树高大于 15m 且胸径大于 30.6cm 的高生长型单株及 31 株树高小于8.7m 且胸径小于 17.5cm 的低生长型单株，之后进行树高、胸径以及球果量的平均值计算分析之后，以树高大于 10.9m 且胸径大于 22.2cm 为参考，筛选出树高、胸径均处中上水平上的148 株，共计 209 株样本数。此外，选取 3 个优良单株用于qPCR 试验，分别取其幼嫩部位，包括针叶（ZY 1-3）、树皮（SP 1-3）、根（SG1-3）、韧皮部（RP 1-3）以及木质部（MZ 1-3）共15 个样本，每株以 3 个生物学重复为标准采样。

2.2 方法

2.2.1 样品总 DNA 提取方法

采集的 209 个叶片样品都是当年无病虫害的新鲜幼嫩针叶，并按每株植物进行编号，然后放入自封袋，标记后放于-80℃超低温冰箱保存以备提取 DNA。实验开始后，加入液氮经研磨后，用改良的 CTAB 法来进行华山松基因组 DNA 的提取，所用到试剂及具体步骤如下：

（1）1mol/L Tris-HCl

称取 12.11g 固体 Tris 缓慢倒于 100 mL 烧杯，然后量取 50 mL 蒸馏水倒入100 mL 烧杯中来溶解固体 Tris，而后用容量瓶定容至 100 mL，调节 pH 到 8.0即可，最后放入高压灭菌锅进行 25 min 灭菌处理。

（2）0.5mol/L EDTA

在 100 mL 烧杯中准确称量 18.61g 固体 EDTA，然后加入 50 mL 蒸馏水和少许固体氢氧化钠进行溶解，最后用容量瓶补足定容至 100 mL，调节 pH 到8.0 即可，最后放入高压灭菌锅进行 25 min 灭菌处理。

（3）5mol/L NaCl

在 200 mL 烧杯中准确称量 58.44 g 固体 NaCl，然后倒入 100 mL 蒸馏水以溶解固体 NaCl，之后用容量瓶补足定容至 200 mL，最后将其放入高压灭菌锅中进行 25min灭菌处理。

（4）2×CTAB

准确称量 4g 固体 CTAB 和 16.36g 固体 NaCl 缓慢倒入 200 mL 烧杯中，然后添加 100 mL 蒸馏水溶解，最后加入 20 mL 1 mol / L Tris-HCl 和 8 mL 0.5 mol/ LEDTA 混合，然后补足定容至 200 mL 容量瓶，之后放入高压灭菌锅进行 25min 灭菌处理。

（5）苯酚-氯仿-异戊醇（25∶24∶1）

严格依据各组分 25∶24∶1 的比例来进行配制混合液，而后缓慢倒入棕色瓶中进行密封保存。

（6）预处理液

量取上述配制的 50 mL 的 1mol/L Tris-HCl，20 mL 的 0.5mol/L EDTA，140mL的 5mol/L NaCl，混合后用容量瓶定容至 500 mL，高温灭菌后常温保存备用。

（7）70%乙醇

（8）PVP（分析纯）

（9）RNaseA 10mg/mL

（10）100%无水乙醇

（11）TE Buffer

（12）

2.2.1.2 步骤

（1）将幼嫩针叶经清水洗干净后晾干，取 50 mg 用剪刀剪碎放到经消毒处理的研钵中，加入一勺液氮快速研磨，待液氮快挥发完的时候加速研磨，有利于研磨充分，之后重复加入液氮，直至研磨彻底；

（2）拿经消毒处理过的药匙小心将粉末状样品装入至 2 mL 离心管中，而后迅速加入 PVP 一药匙进行防止氧化处理，之后再加入 1 mL 预处理液，快速摇匀，放置于冰盒上；

（3）经过 12000 r/min 转速进行为时 1 min 的离心之后，把上清液倒掉；

（4）加入 0.8 mL 提前预热好的 2×CTAB 于 2 mL 离心管中，充分摇匀后，置于泡沫垫上，缓慢放在 65℃的水浴锅里 1 h，期间每隔 6min 摇一次；

（5）时间到之后，取出来冷却至室温，加入 35 μL

和 0.8 mL 酚-氯仿-异戊醇（25∶24∶1），然后进行 10 次左右的摇匀，而后在离心机上以 12000r/min 的转速离心 10 min，最后把上清液分批吸至离心管中；

（6）重复步骤 5 两次；

（7）小心吸取上层清液至新的离心管，加入上清体积的 4 倍 1×CTAB，在室温下放置 1 h；

（8）在 12000 r/min 转速下进行为时 5min 的离心处理，倒掉上层清液，此时白色沉淀出现；

（9）之后加入 2.4 倍体积的已预冷的-20℃无水乙醇，静置 1 h，而后在离心机上以 12000 r/min 的转速离心 5 min；

（10）加入 1 mL 70%的乙醇漂洗沉淀 2 次，而后在离心机上以 12000 r/min 离心的转速 1 min，倒掉乙醇，再以 12000 r/min 的转速离心 1 min，之后吹干或晾干；

（11）加入 3 μLRNaseA 去除 RNA，

（12）溶于 50 μLTE Buffer 中，放置于-20℃冰箱中以供后续使用。

2.2.2 DNA 浓度和纯度检测

华山松基因组DNA浓度和纯度的检测利用0.8%的琼脂糖电泳凝胶和超微量紫外分光光度计。其中0.8%的琼脂糖电泳凝胶是较直观的显示DNA的提取结果，以明亮、无拖尾的单条带为最佳，超微量紫外分光光度计检测的结果可以较准确的反映出DNA的浓度与纯度，其中波长230nm为提取样品中的一些污染物，如碳水化合物、盐（胍盐）等的吸收峰；260nm为核酸吸收峰；280nm为蛋白质的吸收峰，在260 nm和280 nm处的吸光度比值在1.8~2.0之间为较好。

对存在轻度杂质污染的 DNA 样品，可以通过添加 1mol/L 的 NaCl 溶液进行杂质去除。

测序

2.3.1 酶切方案

依托于北京百迈克生物科技有限公司，利用其开发的酶切预测软件，以华山松基因组大小和 GC 含量等信息为依据，采用无参考基因组的方法来最终确定最佳的酶切方案，主要有以下 3 个原则：

（1）使得酶切后的片段在重复序列中的比例降低；

（2）保证酶切片段在目标物种的基因组上呈均匀排列；

（3）得到的酶切片段大小、数量、长度应与实验预期需求一致。

2.3.2 实验流程

以日本晴水稻为对照来进行测序参考依据，根据上述方案，把检测合格的 209个华山松基因组 DNA 全部进行酶切处理（Kozich et al, 2013），然后进行以下操作：

（1）在 SLAF 标签 3′端加 A 尾处理

（2）连接测序接头-Dual-index

（3）PCR 扩增

（4）切胶获取目的条带

（5）最后进行质量检测合格后测序等。

2.3.3 信息分析流程

通过 Dual-index 来辨别测序之后的原始数据，然后对低质量的 reads 进行筛选，最后使用（Q）值评估该实验的测序质量。 Q 值水平是衡量单碱基错误率的重要标准，它对应于测序的单碱基错误率，可以从以下公式可以看出，即 Q 值越高，错误率越低，反之，则越高。

Q_-score=-10×log₁₀P

EcoRV-HF®+ScaI-HF®的酶切效率用日本晴水稻对照数据来评估，从而获得较高质量的 SNP标记，以期达到群体分析的目标，其分析过程如图1所示。

2.3.4SNP 标记的开发

本实验以 SLAF-seq技术为核心，对209份华山松优良单株进行多态性 SNP标记开发，用上述确定的EcoRV-HF®+ScaI-HF®酶来构建简化基因组测序文库，最后在 IlluminaHiSeq2500system上测序，得到SLAF标签，用软件BWA把测序得到的reads进行比对，以次等位基因频率（MAF）和位点完整度为依据，过滤掉MAF<5%和完整度小于80%的低质量SNPs，而后结合GATK和SAMTOOLS软件作为开发 SNP标记的方法（Li,2009；McKennaetal,2010；Lietal, 2009）。

2.3.5全基因组关联分析

GWAS是结合表型数据在全基因组内对SNP进行全体关联的一种方法，将表型和基因型结合起来进行群体的分析，筛选重要部位的阈值P=10^-5（1/N，N为标记数），而在环境的影响下，结果中会出现没有达到第一选择阈值（P <10^-4），但对于那些重复性好的标记位点，也有必要来统计分析，而后来挖掘与表型性状相关的功能基因。

GWAS的技术流程如下：（1）依据育种目的，建立研究群体，选择尽可能大的群体作为研究的样本，然后建立目标性状的数据库。（2）提取样本的DNA并进行质量控制以满足基因分型的要求，检测和质量控制基因型数据以满足后续关联分析的要求。（3）使用适当的统计模型分析SNP与目标性状之间的关联。4）对关联分析后的结果进行进一步分析和验证。

本研究在开发与生长量相关的 SNP分子标记的基础上，利用TASSEL软件进行关联分析，其公式如下所示：

y=Χα+Qβ+Kμ+e

首先使用 admixture软件来进行Q值（样本总体结构）的计算，然后使用SPAGeDi软件进行K值（样本之间的亲属关系）的计算，一般线性模型仅需要计算 Q值，然后使用混合线性模型需要计算Q和K，其中，X表示为基因型，Y表示为表型，由公式统计后可知，每一个SNP标记的位点都可以关联到一个结果（Initiative, 2000）。

候选基因在不同组织中特异性表达的 qPCR分析

2.4.1华山松总 RNA提取方法

用于qPCR试验的材料为 3 个优良单株的幼嫩部位，包括针叶（ZY1-3）、树皮（SP1-3）、根（SG1-3）、韧皮部（RP1-3）以及木质部（MZ 1-3）共15 个样本，每株以 3 个生物学重复为标准采取。

参照昆明硕阳科技有限公司提供的植物RNA小量提取试剂盒进行华山松总RNA的提取，具体操作步骤见下：

（1）将无病虫害的幼嫩针叶、根、经清水洗净，树皮用小刀轻轻刮取，称取50-70mg样品于消毒处理的研钵中，加入一勺液氮快速研磨，待液氮快挥发完的时候加速研磨，有利于研磨充分，之后重复加入液氮，直至针叶彻底研磨成粉末状；

（2）把样本粉末用药匙缓慢分批次加入到1.5mL离心管中，之后立刻加入

混合液，立刻摇匀；

（3）放入离心机，调整为14,000 r/min转速在室温下进行为时5 min离心处理，而后将gDNAFilter过滤柱嵌入2mL的收集管中，吸取上层清液至其中，在室温14,000 r/min转速下进行为时 2min 离心处理；

（4）往滤液中加入等体积的BufferRCB，上下摇匀混匀 10 次左右。

（5）将HiBind®RNAMini结合柱嵌入2mL收集管中，转移500μL混合液缓慢添加至HiBind®RNAMini结合柱上，在室温10,000r/min转速下进行 1min的离心处理，倒掉滤液；

（6）将HiBind®RNAMini结合柱嵌入同一2mL收集管中，转移余下混合液缓慢添加至HiBind®RNAMini结合柱上，在室温10,000r/min 转速下进行 1min 的离心处理，倒掉滤液；

（7）将HiBind®RNAMini结合柱嵌入同一2mL收集管中，往HiBind®RNAMini 结合柱缓慢加入400μLRWCWashBuffer，在室温10,000r/min转速下进行 1min的离心处理，倒掉滤液；

（8）将HiBind®RNAMini结合柱嵌入同一2mL收集管中，往 HiBind®RNAMini结合柱缓慢加入500μLRNAWashBufferII，在室温 10,000r/min转速下进行 1min的离心处理，倒掉滤液；

（9）重复步骤（8）；

（10）将HiBind®RNAMini结合柱套入同一2mL收集管中，在室温10,000r/min转速下进行 2min 的离心处理，甩干 HiBind®RNAMini结合柱基质；

（11）将HiBind®RNAMini结合柱嵌在一个新的1.5mL离心管中，吸取30-50μLDEPC水，准确添加在HiBind®RNAMini结合柱膜中央，常温放置 2min，在室温10,000r/min 转速下进行 1min 的离心处理，洗脱RNA，二次洗脱可增加洗脱效率，将所得到的RNA溶液存放于超低温冰箱中，以供后续实验备用。

2.4.2RNA浓度和质量的检测

（1）RNA浓度检测

用超微量紫外分光光度计进行浓度的测定，把刚提取出来的RNA溶液吸取 1μL，放到仪器上开始操作，在260nm和 280nm处的吸光度比值作为依据判断是否存在杂质污染，即OD₂₆₀/OD₂₈₀值在 1.8~2.0 之间为合格。

（2）RNA质量检测

配制 1％的琼脂糖凝胶：量取80mL1×MOPS电泳缓冲液、称取0.8g 琼脂于烧杯中，放入微波炉中加热溶解完后加入 3μL核酸染料，缓慢顺时针摇匀，对刚提取出来的 RNA进行电泳质量检测，首先先把一次性手套进行消毒处理，而后平铺于冰面上，然后吸取1 μL6×RNALoadingDye按顺序放于其上，再吸取3 μL待测样品RNA与之混合后，用10μL的移液枪将其混合液缓慢吸取添加至凝胶点样孔中，打开电泳仪器，调整电压为120 V，进行为时20min左右的电泳，最后用凝胶成像仪对胶板进行拍照保存，得到的 RNA质量均可满足后续反转录的要求。

2.4.3合成第一链 cDNA

合成第一链cDNA采用昆明硕阳生物科技有限公司反转录试剂盒来进行合成，具体操作步骤如下：

（1）将模板mRNA和其他试剂解冻置于冰上，在1.5mL离心管中依次添加2μL模板RNA、2μLAccuRT ReactionMix(4×)、4μL无酶水，而后在42℃下反应2min或者室温反应5min；

（2）然后依次添加 2μLAccuRTReactionStopper(5×)、4μL5×All-in-OneRTMasterMix、6μL无酶水，而后在 25℃下反应 10min 或者 42℃下反应 15min。

2.4.4第一链cDNA的检测

（1）第一链cDNA浓度检测

在核酸蛋白测定仪下操作进行浓度、纯度测定，把刚提取出来的DNA溶液吸取 1 μL，放到仪器上开始操作，在260 nm和 280 nm处的吸光度比值在 1.8-2.0之间为较好。

（2）第一链cDNA质量检测

采用 0.8%的琼脂糖凝胶对刚提取出来的DNA进行电泳DNA质量检测，具体操作如下：

①清洗容量为 250 mL的三角瓶，静置晾干；

②称取 0.64g琼脂糖，之后量取80mL的1×TAE，缓慢倒入容量瓶中，用微波炉加热4-5 次，期间拿出顺时针摇匀，直到琼脂糖溶解为透明状即可；

③取出后，冷却到 60℃左右，吸取2μL核酸染料添加至三角瓶底部，立即顺时针摇匀，倒入载胶板中，放好梳子，冷却 30 min；

④取掉梳子，将胶板取出，放入已倒好的1×TAE缓冲液电泳槽中，依据实际情况适量即可；

⑤将待检测的DNA溶液吸取2μL和2μL的DNA上样液（LoadingBuffer）混合，用量程为10 μL的移液枪缓慢吸取加入凝胶孔中，并加入2000bp 的 Marker作为标记物；

⑥打开电泳仪器，调整电压为120V，进行为时30min左右的电泳，之后用凝胶成像系统进行观察，检测其质量，最后保存，以备用。

2.4.5CesA2和CCoAOMT基因的引物设计

在 NCBI 中进行逐个序列比对之后，查找与树高、胸径性状相关的基因，利用NCBI网页在线版Primer3.0进行qPCR引物的设计，具体原则如下：引物长度在18~25bp左右即可，退火温度在 55℃~65℃之间较好，前后引物的Tm值不超过5℃，GC 含量在40%~60%之间较好，中间的PCR扩增片段一般在 100~400bp 最好，引物见下表（2）；之后选择与松树相关的AT5G42190作为内参基因来进行 qPCR，其序列为：

ATGCTGGACAGGCTTTGAAC

GAGTTGCTCCGAGATCTTTACA

经筛选后 2个基因各自得到一对在各个组织中稳定表达的引物，分别是引物CE1和 CC5。

表2qPCR分析所用引物

基因名称引物编号前引物后引物

2.4.6对CesA2和CCoAOMT基因的 qPCR验证

进行qPCR时，首先先进行避光处理，而后采用昆明硕阳的配套试剂盒，然后在ABI7500FastReal-TimePCR System仪器上，进行避光qPCR反应。具体操作步骤如下：

确定为 20 μL的反应体系，配套的反应液配制方法见表 3；

表3qPCR 反应体系

具体的qPCR扩增标准程序如下：

第一步：预变性：95℃ 30 s；

第二步：PCR反应：95℃ 5 s、60℃ 30 s40cycles；

第三步：溶解曲线：95℃15 s、60℃ 1min、95℃ 15 s。

经实验结束后，统计Ct值于Excel表中进行计算，利用公式2-△△Ct=2-[(处理组Ct-内参Ct)-(对照组Ct-内参Ct)]计算出各组织的Ct值（Schmittgen＆Livak,2008），而后用 origin 来绘制相应图片。

结果分析

3.1华山松DNA提取

采用改良CTAB法提取的华山松基因组DNA的效果如图 5所示，经1%的琼脂糖凝胶电泳后，在209 个样品中全部提取成功，其条带不仅亮度清晰、杂质少，而且还没有拖带现象，而后用核酸蛋白测定仪对其浓度进行检测（附表1），可以满足后续构建文库的要求。由附表1 可以看出，所提取的华山松总DNA浓度范围在 6.5-944ng/μL之间，相差 145 倍之多，反应了样品间的浓度差异较大，DNA总量平均值为5.54μg，其中样品XSQG-131DNA总量最少为0.1μg，样品XS-30的最多，为34.2μg，此外，209个华山松样品中OD₂₆₀/OD₂₈₀值在 1.60-1.99之间的样品有 139 个，占总量的66.5%，可满足后续实验的要求；而OD₂₆₀/OD₂₈₀值大于2.00 的样品有 70 个，这些样品中存在轻度杂质污染，经处理后可满足后续要求。

3.2SLAF-seq测序数据统计评估

经文库构建后，进行 SLAF-seq测序，然后对数据进行处理分析。

3.2.1测序质量检查

如果某个碱基在测序时出现错误的概率为 0.001，则该碱基的质量值Q应该为30。继而用Q30 来检测，209 个华山松样品中 Q30 值在95.27 -97.39%之间，平均值为96.53%，可以看出测序单碱基错误率很低，华山松样品的测序质量值分布图如5所示：从碱基分布图 5 中能看出，在华山松 209 个样品中， A 和 T，G 和 C 的含量比例呈现一致趋势。

本研究中对209个华山松样品测序共得到2332Mbreads数据，平均每个样11.16Mb，GC含量在 37.15-38.67%之间，平均含量为37.74%，结合以上分析结果部分样品的reads数量、Q30 和GC含量数据结果，见表 4。

表4部分样品测序数据统计

3.3华山松SLAF 标签开发

3.3.1华山松SLAF标签统计

本项目共开发 12952676 个SLAF标签，标签的平均测序深度为 20.24x，部分样品SLAF标签的统计结果见下表 5：

表5部分SLAF标签统计

华山松SNP标记开发

本实验利用北京百迈客生物科技有限公司自主研发的SLAF-seq技术，对209份华山松优良单株进行全基因组SNP标记开发，209份华山松材料共获得了12952676 个SLAF标签，其中，多态性SLAF标签有1456486 个，然后用软件BWA把经过测序所得的reads数据与火炬松（Pinus taeda)的基因组进行完整比对，以次等位基因频率（MAF）和位点完整度为依据，过滤掉MAF<5%和完整度小于80%的低质量 SNPs，而后结合GATK和SAMTOOLS软件作为开发SNP标记的方法，共获得了 3469074 个 SNP标记，结合以上分析结果，部分样品的SNP信息，如表6 所示：

表 6 部分样品 SNP 信息统计表

3.4 全基因组关联分析

3.4.1 多重假设检验校正

本研究把群体结构加入当做协变量，从而建立混合线性模型，这样做的目的是降低了全基因组关联分析中的假阳性，再把 Bonferroni 方法用在校对全基因组关联分析中，这样既保证了关联分析之后，精准定位的结果正确性，又以此来降低了多重假设检验校对后的 P 值假阳性的概率。从下列 QQ 图（图 6）中可以看出：假设在标记位点和表型性状之间没有关联的背景下，观测 P 值和假设没有关联的期望 P 值在图中表现为在右端发生了偏离，这可以说明树高这一性状的差异不是群体分层导致的结果，表明了全基因组关联分析所用的 GLM 模型是适合本研究的。

3.4.2树高和胸径性状的全基因组关联分析

利用开发的华山松基因组 SNP，对表型树高（Height）和胸径（DBH）性状进行 GWAS分析，采用GLM模型获得了与胸径相关的 SNP位点 6 个，与树高相关的SNP位点7个，获得的13个SNP所在SLAF标签上的位置以及对应的碱基序列，见表7，利用SNP的序列在NCBI数据库进行 Blast比对之后，发现2 个基因可能与胸径生长相关，分别为SNPMarker227806 位点相近的与纤维素合成相关的CesA2基因（登录号：AC241331.1），以及Marker227806位点相近的与木质素合成相关的CCoAOMT基因（登录号：AC241331.1）。

表7SNP所在的序列

3.5与胸径生长相关的SNP 标记群体分析

对华山松 209 个样品胸径表型数据结合关联到的 SNP（Marker227806）基因型进行统计后发现，在此群体中有两种基因型，分别是（GG）175 份，（GT）34 份（表8），其中，GT杂合基因型，结合胸径平均值来看，基因型为GG的胸径平均值为 23.80cm，基因型为GT的胸径平均值为 26.10cm，经t检验分析后发现 P =0.037 <0.05，说明SNP（Marker227806）2 种基因型的样本胸径存在显著差异，因此，进一步说明此 SNP与胸径性状相关，有利基因型为GG，可以用此SNP标记来进行挑选优良单株。

表8209株华山松基因型统计

3.6与胸径生长相关基因的 qPCR 组织特异性分析

3.6.1华山松RNA提取

使用昆明硕阳科技有限公司植物RNA小量提取试剂盒进行15个华山松总RNA的提取，效果如图7 所示，在琼脂糖凝胶上展现出清晰明亮且无拖带的2 条带，满足后续qPCR的最低要求。

3.6.2候选基因的 qPCR 组织特异性表达分析

经关联分析后，选用与松树相关的 AT5G42190作为内参基因来进行qPCR结果验证分析华山松在不同组织中的表达情况，所用材料，包括针叶（ZY）、树皮（SP）、根（SG）、韧皮部（RP）以及木质部（MZ）共15 个样本，包含三个生物学重复，得到Ct值，从而利用上述所用公式计算出各组织的 2^-△△Ct值，进而得到各组织的表达量图，从图9中可以看出，与纤维素合成相关的基因CesA2在各个组织中都有表达，在RP组织中的表达量明显比ZY、SP、SG、MZ的高，其中RP组织中的表达量最高，SP组织中的表达量最低；而CesA2基因被证实与次生壁生长相关，进一步证明很可能在胸径的生长发育过程中发挥着重要作用；在图 9中的与木质素合成相关的基因CCoAOMT在各个组织中均有表达，在RP组织中的表达量明显比ZY、SP、SG、MZ 的高，其中在RP组织中的表达量也是最高，SP组织中的表达量相对来说较低。此基因在植物体内表达的量严重影响木质素单体的含量，从而影响其在木质素总量中的比例，这进一步说明在胸径生长发育中起到关键作用。

以上对本发明的具体实施例进行了描述。需要理解的是，本发明并不局限于上述特定实施方式，本领域技术人员可以在权利要求的范围内做出各种变化或修改，这并不影响本发明的实质内容。在不冲突的情况下，本申请的实施例和实施例中的特征可以任意相互组合。

Claims

1.一种华山松基因组SNP的开发方法，其特征在于：包括如下步骤：

S2、利用开发的华山松基因组 SNP，对表型树高（Height）和胸径（DBH）性状进行 GWAS分析，采用GLM模型获得了与胸径相关的 SNP位点 6 个，与树高相关的SNP位点7个，利用SNP的序列在NCBI数据库进行 Blast比对之后，发现2 个基因可能与胸径生长相关，分别为SNPMarker227806 位点相近的与纤维素合成相关的CesA2基因，以及Marker227806位点相近的与木质素合成相关的CCoAOMT基因。

2.一种华山松基因组SNP的应用，其特征在于： SNPMarker227806可用于挑选华山松优良单株。