基于SLAF-seq 的华山松球果数量性状全基因组关联分析*

李启少，叶 鹏，赵文植，马路遥，王正德，曹正英，王 飞，辛培尧

(1.西南林业大学，国家林业和草原局西南风景园林工程技术研究中心，云南 昆明 650224；2.西南林业大学，西南地区生物多样性保育国家林业和草原局重点实验室，云南 昆明 650224；3.内蒙古浑善达克规模化林场，内蒙古 呼和浩特 010020)

华山松(Pinus armandiiFranch.)是松科(Pinaceae)松属(Pinus)大乔木，也是中国特有的五针松树种，其耐寒能力强，喜好温和凉爽、湿润的气候，适应性较广[1]。华山松是优良的建筑和家具材料，同时也可以用于制浆造纸，种子可供食用和榨油，利用价值较高。近年来，由于种质资源退化、良种选育滞后，导致华山松人工林陆续出现生长量下降、长势衰弱、病虫害严重、利用水平低和产品附加值低等问题[2-4]。虽然国内相关科研单位已开展了对华山松的遗传改良工作[5-9]，但由于其育种周期长和主要性状的优劣难以在早期进行选择，因此，找到准确、可靠且有效的早期选择标记是解决良种问题的关键技术之一。

随着分子生物学的不断进步，林木育种的方法与理论也在不断发展，尤其是分子标记辅助选择技术的不断更新，涌现出多种分子标记方法，如RAPD、AFLP、SRAP、SSR 和SNP 等标记，为植物遗传改良提供了更多的选择，也提高了选择的可靠性与准确度[10]。其中，单核苷酸多态性(single-nucleotide polymorphism，SNP)标记是由单碱基突变(缺失、插入、转换和颠倒)所引起的基因组水平上的DNA 序列多态性，具有稳定性高、位点丰富、分布较广和检测迅速等特点，目前已广泛应用于植物育种[11]。全基因组关联分析(genome-wide association study，GWAS)由RISCH等[12]首先提出。该技术能够在全基因水平上发掘与复杂性状相关联的遗传变异，以连锁不平衡为基础，通过识别和定位目标群体中的高密度分子标记获得与复杂性状表现型变异有关联的SNP 标记[13]。GWAS 最早应用于人类视网膜黄斑变性与人类年龄的相关性分析上[14]。近年来，由于引入混合线性模型和Bonferroni 校正统计方法，能够将自然群体特征性状和非家系群体作为研究对象开展GWAS 研究。因此，GWAS 也被大量应用于植物的关联分析与遗传育种工作。目前，不仅在水稻[15]、小麦[16]和玉米[17]的产量、株高、穗长、成花时间及抽穗期等性状关联出多个相关的基因位点，还在植物抗病基因的筛选[18]和代谢物合成[19]等方面得到了应用。GWAS 相关研究结果为阐明生物复杂性状的遗传结构提供了理论基础。生产中可将检测到的关联位点运用于分子标记辅助选择，对植物的遗传改良工作具有重要的理论和实践意义。

本研究基于前期通过SLAF-seq 测序获得的高通量SNP[6]构建华山松样品的系统发育树与主成分聚类图，继而进行全基因组关联分析，筛选与华山松球果数量性状显著关联的SNP 位点，并对显著关联位点区域的基因进行分析，挖掘与球果数量性状相关的候选基因，分析候选基因在华山松6 个部位的qPCR 基因组织特异性表达。研究结果可为华山松球果数量相关性状的早期选择提供理论依据。

1 材料与方法

1.1 采样地概况

采样地位于云南省楚雄市紫溪山华山松无性系种子园(N24°58′58″～25°4′，E101°22′29″～101°26′7″，海拔2 200～2 400 m)，于1990 年建成，园区内80%的林地坡度小于10°。种子园地处北亚热带温凉湿润气候，年降水量约1 000 mm，相对湿度80%～85%，年平均气温12 ℃，平均最低气温-6 ℃，平均最高气温20 ℃。

1.2 供试材料

试验材料选自种子园内6 个种源地的1 453株华山松单株，以连续多年的年份内球果数量为基础数据，共筛选出210 株华山松作为试验材料，其中高、低球果产量的华山松单株各30株，球果数量及生长量均处中上水平的华山松单株150 株。另选3 株不同华山松优良单株的幼嫩且无病害的球果、树皮、韧皮部、木质部、根和针叶，用于提取华山松RNA 以进行q-PCR 验证。

1.3 试验方法

1.3.1 系统发育树与主成分分析

利用前期已开发的SNP 标记，采用MGEA X软件[20]，基于邻接法的Kimura2-parameter 模型构建华山松样本的系统发育树，bootstrap 重复1 000次。利用EIGENSOFT 软件[21]将SNP 标记用于主成分分析，获得华山松的主成分聚类图，每1 种颜色代表1 个种源地，每1 个点代表1 个样品，样品间存在广泛的差异，样品距离越近，则代表亲缘关系越近，主成分分析的结果可以辅助进化分析，充分体现群体中个体亲缘关系的离散程度。

表1 qPCR 引物序列Tab.1 qPCR primer sequence

1.3.2 球果数量的GWAS 分析

(1) GWAS 分析：将球果数量作为鉴定指标，利 用TASSEL[22]、FaST-LMM[23]和EMMAX[24]软件的光谱变换线性混合模型(factored spectrally transformed linear mixed models，FaST-LMM)、一般线性模型(general linear models，GLM)、混合线性模型(mixed linear models，MLM)和压缩混合线性模型(compressed linear mixed model，CMLM)，在充分考虑210 份华山松单株的总体结构与亲缘关系的情况下进行关联分析。利用admixture 软件[25]获得Q值(样本总体结构)；样品间亲缘关系K值则是利用SPAGeDi 软件[26]获得；X表示基因型，y表示表现型。所有SNP 位点均可得到1 个关联结果。其中，TASSEL 软件的混合线性模型公式为：y=Xα＋Qβ＋Kμ＋e。式中：α 为固定效应向量；β 为SNP 替代效应；μ 为服从分布N (0，Gσ2)的随机加性遗传效应的向量，G 为从SNP 标记导出的亲属关系矩阵，σ2为加性方差；e 为残差。GWAS 结果基于R v4.2.0 中的CMplot 程序包绘制，Q-Q 图(quantile-quantile plot)用于评估关联分析结果的准确性。

(2) 候选基因筛选：关联分析得到与球果数量性状相关的SNP 位点，并在SLAF 标签上找出对应的位置以及碱基序列，利用SNP 序列在NCBI数据库进行BLAST 比对，筛选与球果数量性状相关的基因。

1.3.3 荧光定量PCR 反应

(1) 华山松总RNA 提取

使用植物RNA 提取试剂盒(OMEGA R6827)提取华山松针叶、树皮、韧皮部、木质部、球果和根的总RNA。

(2) mRNA 反转录为cDNA

参照周延清等[27]的方法得到第1 链cDNA 并储存于-20 ℃冰箱中，备用。

(3) 引物设计

从NCBI 中检索相关的基因序列，用Primer 3设计荧光定量PCR 引物(表1)。引物条件为：引物长度在18～24 bp 之间；GC 含量在40%～60%之间；5′端或中间区域为G 或C；扩增片段的长度在100～300 bp 之间，退火温度在55～60 ℃之间。从设计的7 对引物中筛选1 对能让候选基因在各个组织中稳定表达的引物，对检测到的目标基因进行q-PCR 分析。

(4) 荧光定量PCR 反应

使 用EvaGreen 2×qPCR MasterMix 试剂盒(上海爱必梦生物科技有限公司)进行荧光定量PCR 反应，反应体系为20 μL，包括EvaGreen2×qPCR MasterMix 10 μL，上、下游引物各1 μL，模板DNA 1 μL，无酶水7 μL。q-PCR 扩增程序为：95 ℃预变性 10 min；95 ℃ PCR 反应15 s，60 ℃ 1 min，共40 个循环；95 ℃溶解曲线15 s，60 ℃ 1 min，95 ℃ 15 s。选择和松树相关的AT5G42190作为内参基因，其引物序列为：F：ATGCTGGACAGGCTTTGAAC，R：GAGTTGCTCCGAGATCTTTACA。

1.3.4 数据统计与分析

试验结束后，统计Ct 值于Excel 表中，利用公式ΔΔCt=(Ct处理组-Ct内参)-(Ct对照组-Ct内参)]计算各组织的Ct 值[28]，之后利用Excel 软件作图；采用Excel 2010 和SPASS 20.0 软件进行相关数据的描述统计与相关性分析。

2 结果与分析

2.1 系统发育树与主成分分析

系统发育树(图1)显示：华山松各材料之间存在着较为丰富的遗传变异，来自相近地理位置的材料都基本在相近的位置，但是来源于楚雄和巍山的材料分布相当广泛，在各群体中均有分布。将SNP 标记用于主成分分析(图2)发现：同一种源地的样品基本聚在一起，说明它们的亲缘关系较近，其中来自巍山的样品分布范围较广泛，表明巍山各材料间遗传差异较大，与系统发育树的分析结果基本一致。

图1 华山松系统发育树Fig.1 Phylogenetic tree of Pinus armandii

图2 华山松PCA 聚类图Fig.2 PCA cluster graph of P.armandii

2.2 华山松球果数量性状的GWAS 结果

2.2.1 华山松全基因组关联分析

GWAS 结果表明：球果数量性状能够在FaSTLMM、GLM 和MLM 模型中被定位。由图3 可知：FaST-LMM 和MLM 模型的观测值与期望值仅在最右侧发生了偏离现象，说明球果数量性状的差异并非因群体分层所造成，所选的模型为本试验适合的模型；而GLM 模型中SNP 位点出现在预测线上方，且不与对角线重合，说明观测值大大超过了预测值，表明该模型不合理。

图3 华山松球果数量Q-Q 图Fig.3 Quantile-quantile (Q-Q) plot of cones of P.armandii

共检测到12 个与球果数量性状显著相关的SNP 位点，获得12 个SNP 所在SLAF 标签上的位置和P值(表2)。12 个SNP 位点分别为Marker170924、Marker193307、Marker357784、Marker6337143、Marker17681543、Marker201128、Marker215196、Marker250033、Marker251795、Marker-6508087、Marker292650、Marker18895864。

表2 与华山松球果数量显著相关的SNP 位点Tab.2 SNP loci significantly related to the number of cones in Pinus armandii

2.2.2 候选基因的筛选

将12 个SNP 所在的序列片段在NCBI 中进行基因检索(表3)，在Marker193307 附近检测到1 个紧密关联的基因Korrigan(KOR1，编码跨膜内切β-1,4-D 葡聚糖酶)，其序列相似度较高(为74.74%)，NCBI 基因登录号为AC241332.1，该基因对植物纤维素合成具有重要作用。

表3 与球果数量相关的单核苷酸多态性(SNP)序列Tab.3 Single-nucleotide polymorphism (SNP) sequence related to the number of cones

2.3 基因组织特异性分析

由图4 可知：KOR1在韧皮部的表达量最高，在幼嫩球果中的表达量次之，在根中的表达量最低。已知KOR1基因是与纤维素合成相关的基因，该基因在韧皮部和幼嫩球果中的表达量较高，说明其参与了韧皮部与幼嫩球果的生长代谢。因此，推测球果数量性状与该基因相关。

图4 KOR1 在华山松各个部位中的表达量变化Fig.4 Expression changes of KOR1 gene in different tissus of P.armandii

3 讨论

3.1 华山松不同分子标记方法比较

近年来，分子标记辅助选择技术不断更新，越来越多的分子标记方法应用于华山松育种。赵杨等[29]对2 个不同年份的华山松子代材料进行ISSR 遗传多样性分析，结果显示：该种子园的遗传多样性处于较高水平，遗传变异多处于种源内，与种子园内亲代群体相比较，初期子代的遗传多样性指数略有降低，而后期子代的遗传多样性指数则为最低，可见，对于提高华山松球果产量和品质刻不容缓。辛静等[30]对紫溪山华山松无性系种子园进行SSR 遗传多样性分析，结果表明：种源间遗传变异不高，其主要遗传变异存在于种源内。徐剑等[31]利用SCoT 分子标记对紫溪山华山松种子园进行遗传多样性分析，表明种子园内不同无性系拥有较高水平的遗传多样性，且遗传分化主要存在于种源内。刘成等[32]同样利用紫溪山华山松无性系种子园的材料，通过SRAP分子标记技术对其进行遗传多样性分析，结果显示：该园的华山松群体遗传多样性处于相对较高的水平，遗传变异在种源内。因此，利用该种子园进行华山松杂交育种时，应多进行种源内杂交，兼顾遗传距离较远的种源间杂交。同时，需尽量开发较多的早期选择标记，以实现早期选择。

由于分子标记方法的不同，开发基因标记的试验结果会有一定差异。本研究所采用的分子标记方法是SNP 标记技术，与其他分子标记方法相比较，SNP 标记数量多、遗传稳定，且不需要凝胶电泳检测，为试验节省了大量的时间和成本。而利用如简单重复间序列分子标记(ISSR)、序列相关扩增多态性分子标记(SRAP)和简单重复序列分子标记(SSR)等方法开发华山松基因组的DNA 分子标记难度相对较大、需要的时间较长，且标记数量有限。赵杨等[29]利用ISSR 分子标记共获得13 条具有多态性的引物；刘成等[32]利用SRAP 分子标记共筛选出15 对多态性引物组合；祝娟[33]利用SSR 标记对华山松基因组DNA 进行PCR 扩增，筛选获得12 对具有多态性的引物。而本研究利用SLAF-seq 测序技术共得到3 469 074个SNP 标记。可见，SLAF-seq 测序技术开发分子标记的效率较其他分子标记技术效率较高。

3.2 华山松球果数量性状的相关基因

本研究共检测到与华山松球果数量显著相关的SNP 位点12 个，并在Marker13307 位点中检索到KOR1基因。KOR1基因是一种纤维素酶，它参与调控纤维素合成，对植物生长发育、组织脱落以及维持植物正常生理代谢有着极其重要的作用[34]。李萍[35]对新疆杏中早熟、中熟和晚熟品种的果实发育期进行了生理生化机制研究，结果表明：伴随着果实的发育，果实中纤维素以及半纤维素的含量先呈上升趋势，14 d 后急剧下降，而纤维素酶在果实生长发育过程中则一直呈上升趋势。另有学者以壶瓶枣作为试验材料，测定其果实的果肩果皮、果肩果肉、果底果皮和果底果肉在盛花后40、55、70 和85 d 的纤维素含量变化，结果表明：壶瓶枣果实整体的纤维素含量呈先上升后下降的趋势[36]，这与李萍[35]对新疆杏果实的研究基本一致。可见，纤维素在果实的生长发育中起到重要的作用，而参与调控纤维素合成的KOR1基因也在果实的生长发育中有着不可或缺的地位。KOR1是纤维素合成所必需的基因，在植物和其他生物(如细菌)质膜—细胞壁的合成中也发挥着重要作用[37-38]。目前，已经确定KOR1与纤维素合成酶复合物具有物理相互作用，并参与了纤维素伸长过程中甾醇连接引物的断裂或纤维素纤维结晶度的降低，KOR1等位基因突变会导致纤维素合成受阻以及生长迟缓，在KOR1基因点突变的irx2 突变体中，束间纤维和木质部的纤维素含量降低，导致导管塌陷，表明KOR1参与次生壁中纤维素的生物合成[39-41]。此外，KOR1的其他突变会导致细胞板形成缺陷[42]和初级细胞壁纤维素含量减少[43-44]，表明KOR1在初级细胞壁的纤维素生物合成中也起着重要作用。本研究的qPCR 结果显示KOR1基因在华山松幼嫩球果中的表达量较高，因此认为KOR1基因参与了幼嫩球果的生长发育。

4 结论

巍山种源的华山松较其他种源拥有较为广泛的遗传变异。结合华山松表型数据和SNP 标记进行全基因组关联分析，得到12 个与球果数量性状显著关联的SNP 标记，并通过检索得到1 个与植物纤维素合成相关的基因KOR1，该基因在幼嫩球果中的表达量较高，可能与球果数量性状相关。