Genetic diversity analysis and core collection of pinewood nematodiasis-resistant Pinus massoniana germplasm resources

2018年，在国家马尾松种质资源库中依据松脂和抗性的关系，从1 207株树中选择120份候选抗性马尾松无性系，同年在松材线虫病重灾区严重感染林分中收集健康树242份，并将所有无性系嫁接保存于浙江省临海市林业技术推广和场圃旅游服务总站(28°88′N，121°01′E)。2018—2020年连续3 a高强度对收集区内的候选抗性马尾松无性系进行人工接种松材线虫测定，筛选存活健康株系，共获得103份抗性马尾松无性系，另取安徽省林业科学研究院审定的抗松材线虫病良种(为目前中国所有的抗松材线虫病良种)作为对照。合计114份，其中浙江37份、江西26份、安徽42份、福建6份、贵州1份、广西2份。采集抗性样株幼嫩针叶，在−80 ℃冰箱中保存，各试验材料编号及来源见表1。

马尾松针叶基因组DNA采用北京艾德莱生物科技有限公司生产的改良CTAB植物基因组DNA快速提取试剂盒提取。利用超微量分光光度计，根据在260 nm的吸光度D(260)和D(260)/D(280)确定DNA样品的纯度和质量浓度，经检验合格后的DNA用TE溶液或去离子水(灭菌后)稀释至20~50 mg·L⁻¹，置于−20 ℃的冰箱保存备用。

参考董虹妤等^[18]和YANG等^[19]使用的马尾松 SSR 引物，用质量浓度为1.2%的琼脂糖凝胶电泳检测，筛选条带清晰且多态性较高的SSR引物(图1)，最终筛选出具有高多态性引物14对(表2)，引物序列由浙江尚亚生物技术有限公司合成。

Figure 1.  Primer polymorphism detection diagram

利用选取的14对 SSR 引物对114份马尾松DNA样本进行PCR扩增，设定本研究的反应体系为25.0 μL：2×Taq Plus Master Mix (Nazyme Code：P211-03)12.5 μL，引物F和R (10 μmol·L⁻¹)各1.0 μL，DNA模板2.0 μL，ddH₂O 8.5 μL。PCR扩增反应在Takara PCR Thermal cycler上进行，其程序为：95 ℃预变性3 min；95 ℃变性15 s，60 ℃退火15 s，72 ℃延伸30 s，35个循环；72 ℃延伸5 min，4 ℃保存。所得PCR产物采用Qsep100全自动毛细管电泳核酸分析仪进行检测和基因分型。

利用GenAlex 6.5^[20]计算各位点的遗传多样性参数，包括等位基因数(N_a)、有效等位基因数(N_e)、观测杂合度(H_o)、期望杂合度(H_e)、Shannon’s多样性指数(I)、近交系数(F_is)、固定指数(F)和基因流(N_m)，并对收集的抗松材线虫病马尾松种质进行主坐标分析(PCoA)。采用Cervus 2.0^[21]获得微卫星位点的多态信息含量(P_IC)。利用Structure 2.3.4^[22]对抗性马尾松种质资源进行分类，推测最佳聚类组群数目，确定群体结构，并对各亚群的遗传多样性参数进行分子方差分析(AMOVA)。

对抗松材线虫病马尾松种质进行群体结构分析后，利用Power Core^[23]和Core Finder^[24]对各亚群基于不同算法构建核心种质。其中，Power Core是利用随机取样策略和启发式算法，寻找从初始阶段到最终阶段的最优路径进行核心集合选择；Core Finder基于NP-完全覆盖问题，采用拉斯维加斯式(LasVegas-style)随机算法的M策略进行样本选择。在构建核心种质过程中，Power Core和Core Finder会自动生成合理的取样比例。最后利用SPSS 26.0和Excel中t检验对原有种质与核心种质各遗传多样性指标进行差异显著性分析，同时运用PCoA和Powermarker^[25]对所构建的核心种质构建UPGMA进化树，确认构建的核心种质的代表性。

由表3可知：14对 SSR 引物在114份抗松材线虫病马尾松样本中共检测出115个等位位点，N_a为5~11，平均为8.17。所有材料的N_e为3.45~8.05，平均为5.54。H_o和H_e的分别为0~0.57和0.51~0.73，H_e均值(0.64)高于H_o均值(0.06)，F_is平均为0.92，综合表明该抗性马尾松种质中存在杂合子缺失现象。P_IC为0.82~0.95，平均为0.9，I为1.15~1.90，平均为1.51，表明所选引物多态性高。

Structure分析得出的结果经Structure Haevester处理后发现(图2A)：lnP(D)(马尔科夫链不作数迭代后所得的对数似然值)持续增大，没有拐点，可用∆K [L(K)在连续K之间变化率的增量]来确定最优分组数(K)。由图2B可知：当K=4时，∆K取得最大值，表明参试的114份马尾松样本被划分为4个亚群。亚群Ⅰ(绿色)包含了12份种质资源，其中2份来自广西，1份来自贵州，3份来自福建，2份来自安徽，4份来自浙江；亚群Ⅱ(黄色)包含24份种质资源，其中8份来自安徽，7份来自江西，9份来自浙江，该亚群有12份抗性马尾松种质掺杂有来自红色部分的基因；亚群Ⅲ(红色)包含了29份种质资源，其中10份来自安徽(均为皖马抗良种)，8份来自江西，11份来自浙江，共有12份抗性马尾松种质掺杂有来自黄色部分的基因；亚群Ⅳ(蓝色)包含49份种质资源，其中2份来自福建，21份来自安徽(其中1份为皖马抗良种)，11份来自江西，15份来自浙江(图3)。由Structure输出的分组结果可知(表4)：筛选获得的抗性马尾松种质中有75.4%的个体样本隶属不同亚群的比例(Q)＞0.8，仅有12.3%的个体Q＜0.6，表明抗性马尾松种质群体结构较强，不同来源的个体被聚在相同的遗传结构中表明不同省间的马尾松存在基因渗入现象。

Figure 2.  Line plots of lnP(D) and ∆K as a function of K

Figure 3.  Structure cluster diagram of 114 samples of resistant P. massoniana

基于遗传距离矩阵对抗性马尾松材料进行主坐标分析(PCoA)。由图4可知：全部马尾松种质资源可大致分为3组，主坐标1和2分别解释了位点信息数据中10.45%和6.74%的变异。基于抗性马尾松114份样本主坐标分析图可知：大部分地理来源一致的抗性马尾松种质分布比较集中，但来自安徽、江西、浙江的部分抗性马尾松种质在3个组中均有分布。结合群体结构分析结果可知：114份抗性马尾松种质资源的分布情况基本符合群体结构分析。

Figure 4.  Principal coordinate analysis of 114 samples

由表5可见：抗性马尾松种质各亚群的遗传多样性整体水平较高。亚群Ⅳ的等位基因数在4个亚群中最大(13.43)，亚群Ⅰ具有最低的等位基因数(6.07)。抗性马尾松4个亚群所包含的马尾松个体数不同，但亚群间遗传多样性水平差异不明显。亚群Ⅳ的N_e最高，为7.14，而最低的是亚群Ⅰ，为4.87。亚群Ⅳ的H_o和H_e均最高，但与具有最小值的亚群Ⅲ差异不显著，分别高0.03和0.15。

AMOVA方差分析表明(表6)：抗性马尾松遗传变异14%存在于亚群间，而亚群内变异为80%，表明亚群内的遗传变异是抗性马尾松变异的主要来源。4个亚群间F_ST为0.143，表明亚群间存在中等程度遗传分化。

利用M策略和随机取样策略构建核心种质的遗传多样性指标见表7。结果显示：基于随机取样策略的Power Core和基于M策略的Core Finder分别筛选得到的核心种质均保留了原有种质100%的等位基因数。但Core Finder核心种质包含的种质(72份)少于Power Core (79份)，且两者间各遗传多样性参数差异不显著(P＞0.05)。综合以上分析表明：选择M策略抽取72份抗性马尾松种质能够以最小的种质份数最大程度地代表收集区内抗性马尾松种质资源的遗传多样性。利用M策略构建的核心种质中包含广西2份(100.0%)、贵州1份(100.0%)、福建4份(66.7%)、安徽33份(其中4份为皖马抗良种)(78.6%)、江西8份(30.8%)和浙江24份(64.8%)种质材料。

由表8可见：构建的72份抗性马尾松核心种质占114份原有种质的63.16%，其N_a均值、N_e均值、I均值分别占原有种质相应遗传参数的100.00%、95.67%和94.96%。t检验结果表明：列入研究的抗性马尾松种质与入选的核心种质的遗传多样性无显著差异(P＞0.05)，说明本研究基于M策略构建的抗性马尾松核心种质能充分代表列入研究的原抗性马尾松种质资源的遗传多样性，剔除了与114份原有种质重复或亲缘关系较近的材料。

PCoA分析(图5)显示：抗性马尾松核心种质均匀分布在整个马尾松种质资源中，表明本研究所构建的抗性马尾松核心种质具有很好的代表性。UPGMA进化树(图6)显示：72份抗性马尾松可分为3个类群，其中类群Ⅰ的11份抗性马尾松种质全部位于亚群Ⅰ(91.7%)，类群Ⅱ的30份抗性马尾松种质全部位于亚群Ⅳ(61.2%)，类群Ⅲ的抗性马尾松种质18份位于亚群Ⅱ(75.0%)，其余13份位于亚群Ⅲ(44.8%)。结合结构分析可知：位于亚群Ⅱ中混杂有黄色和红色的核心种质与亚群Ⅲ中的核心种质聚为一个亚类，亚群Ⅱ中以黄色为主的核心种质聚为另一个亚类，但抗性马尾松核心种质的UPGMA聚类结果与列入研究的原有种质结构分析和PCoA分析结果基本一致，进一步对核心种质进行了确认。

Figure 5.  Principal coordinate analysis of 72 core germplasm and 114 original germplasm

Figure 6.  UPGMA clustering map of 72 samples of resistant P. massoniana

遗传多样性在生物育种工作中至关重要，是评价生物多样性的重要组成部分，物种的遗传多样性越高或遗传变异越丰富，其对环境变化的适应能力也就越强^[26]。本研究利用14对 SSR 引物对经高强度人工接种松材线虫后存活下来的103份马尾松种质和安徽省林业科学研究院审定的11份高抗良种进行遗传多样性和群体结构分析。N_a是衡量SSR位点多态性高低的重要指标，其中所有种质的平均N_a为8.17，高于沈敬理等^[27]利用SSR分子标记技术对131个马尾松无性系的研究(N_a=3.67)，可能是本研究选择的 SSR 引物具有较高多态性和基因分型检测手段不一样造成的。N_e、H_e、I等是评价遗传多样性的重要指标。高景斌等^[28]利用 SCAR 分子标记对抗性马尾松材料进行遗传多样性研究，其N_e和I分别为1.51和0.46；YANG等^[19]运用20对 SSR 引物分析马尾松二代亲本种质的遗传多样性，其N_e为3.01， I为1.26，P_IC为0.89；董虹妤等^[18]对选自马尾松第2代育种群体的12份马尾松材料运用13对 SSR 引物进行研究，其中N_a和I分别为3.17和0.47。与上述研究结果相比，本研究所选择的14对 SSR 引物能有效用于检测马尾松的遗传多样性(N_e=5.54、H_e=0.64、I=1.51、P_IC=0.90)，同时表明筛选获得的103份抗性马尾松种质和安徽省林业科学研究院审定的11份良种具有丰富的遗传多样性，在未来的抗性育种进程中具有很大的遗传增益潜力。

Structure结构分析、UPGMA聚类分析和PCoA分析被广泛应用于推断植物的群体结构，但三者的统计原理不同。PCoA分析是一种非约束性的数据降维方法，通过对特征值和特征向量进行排序来寻找主坐标，从而有效地绘制所分析的样本；Structure结构分析基于贝叶斯模型的计算方法，根据基因型信息对多组样本进行分类；UPGMA聚类分析利用非加权组平均法对参试样本进行层次聚类。PCoA和UPGMA聚类分析可直观看到材料间的分类关系，但Structure结构分析能清晰地反映材料间的遗传背景和基因交流情况^[29−30]。总体而言，Structure结构分析、UPGMA聚类分析和PCoA分析相辅相成，两两使用有助于全面了解不同样本间的遗传差异。本研究Structure结构分析结果表明：114份抗性马尾松种质资源被分成4个亚群，而PCoA分析将所有抗性马尾松种质划分为3组，将结构分析所得4个亚群中的亚群Ⅱ和亚群Ⅲ合并为一组，可能是因为这2个亚群中的部分抗性马尾松种质Q＜0.6，存在一定频率的基因交流和渗透。在结构分析中，每个亚群内都包含来自不同地区的个体，地理来源一致的抗性马尾松种质并没有完全聚在一起，这可能是由于不同地区育种者之间亲本交换或频繁的人工引种，导致新环境中种质的遗传信息发生了改变。安徽省林业科学研究院审定的11份抗松材线虫病良种是中国所有的抗松材线虫病马尾松优良种质，已在安徽省松材线虫病疫区进行了大面积的林相改造和造林，在松材线虫病防控中发挥了重要作用。这11份抗松材线虫病良种主要分布在亚群Ⅲ中，仅与19份抗性马尾松种质聚为一个亚群，这表明本研究筛选获得的大部分抗性马尾松优良品系与审定的良种遗传关系较远，具有较低的遗传相似性。因此，本研究筛选获得的103份抗性马尾松种质与审定的11份良种遗传信息同质性程度较低，且具有较高的遗传多样性，可保障今后营造的抗性林在生态系统上具有较强的稳定性。

传统的核心种质构建方法是对比不同种质资源在植物学、形态学等性状上的差异来剔除亲缘关系相近的样本并且抽取核心种质，但表型性状容易受到外界因素影响而表现不稳定^[31]。基于分子标记数据构建核心种质可减少温度、气候、环境等外界因素的影响，能够很好地作为植物形态和生理特性的有效补充。本研究根据SSR分子标记基因分型数据，利用Core Finder和Power Core，基于M策略和随机取样策略构建抗松材线虫病马尾松核心种质。综合对比分析不同构建策略，得到72份抗性马尾松核心种质，其保留了列入研究的原有种质全部的等位基因数，且各遗传多样性参数如N_e、I、H_e、P_IC分别为5.30、1.43、0.63、0.92，均与114份原有种质的遗传多样性参数无显著差异(P＞0.05)。本研究构建的72份抗性马尾松核心种质在广西、贵州、福建、安徽、江西和浙江6省中均有分布，但江西抗性马尾松核心种质抽取比例为30.8%，低于其他省份抗性马尾松种质的入选比例，表明该省抗性马尾松样本存在较高的遗传相似性。结合结构分析可知：72份抗性马尾松核心种质全面涵盖了列入研究的抗性马尾松种质的4个亚群，亚群Ⅲ抗性马尾松核心种质入选比例最低(44.8%)，表明该亚群抗性马尾松种质遗传信息比较一致；亚群Ⅰ的抗性马尾松核心种质入选比例最高(91.7%)，表明该亚群可能存在最丰富的遗传多样性。

取样比例的大小是核心种质库能否构建成功的重要因素。WANG等^[32]和ESCEIBANO等^[33]认为核心种质取样比例为5%~40%，大多数为5%~15%。李魁鹏等^[34]利用SSR分子标记从864份杉木Cunninghamia lanceolata中筛选获得了80份核心种质，占原有种质的9.3%；黄小凤等^[35]利用SNP分子标记从221份荔枝Litchi chinensis品种(品系)中获得了44份核心种质，占原有种质的20%；张馨芳等^[36]基于SSR分子标记数据，从161份板栗Yanshan chestnut中构建了含46份种质的核心种质，取样比例为28.6%。但本研究核心种质取样比例为63.13%，高于其他木本植物的取样比例。李自超等^[37]认为核心种质取样比例通常由原有种质的数量、群体结构和遗传多样性水平等决定，对于遗传多样性丰富、原有种质数量较小的群体，可适当增大取样比例。

利用14对 SSR 引物对114份抗松材线虫病马尾松种质的遗传多样性分析表明：所选抗性马尾松种质具有丰富的遗传多样性(N_a=8.17，N_e=5.54，H_e=0.64，I=1.51，P_IC=0.90)，且本研究筛选的103份抗性马尾松材料遗传结构简单、遗传变异丰富，可用于开展马尾松重要性状的关联分析研究。根据SSR分子标记数据基于M策略构建了含有72份材料的抗性马尾松核心种质，保留了列入研究的原有种质100%的等位基因。经过t检验、PCoA分析和UPGMA聚类分析，构建的抗性马尾松核心种质与列入研究的原有种质之间无显著差异，剔除了遗传相近的材料，具有良好的代表性。