仿刺参疣足数量SNP 遗传力评估*

倪 萍 任 强 王 静 丁 君 常亚青王扬帆① 胡景杰 包振民

(1. 中国海洋大学 海洋生物遗传学与育种教育部重点实验室 青岛 266003；2. 大连海洋大学 农业农村部北方海水增养殖重点实验室 大连 116023)

仿刺参(Apostichopus j aponicus)是中国传统的海珍品之一，营养价值高，药用价值广泛(常亚青等,2004)。海参养殖业在我国经过十几年的迅猛发展，已成为海水养殖的主要产业之一(王秀利等, 2006)。然而，由于刺参养殖的过速发展和不规范运作，基础研究滞后于刺参养殖产业的发展，加之基础设施、养殖工艺和方式的滞后，出现刺参病害频发、种质资源匮乏、成品参品质下降、单位面积产量和产值下降等问题(李成林等, 2010; 张春云等, 2004; 常亚青等,2006; 赵帅等, 2016)；野生资源过度捕捞和养殖刺参长期近亲繁殖、累代养殖，也使刺参品种质量逐步下降(封岩, 2019)，刺参的养殖业面临严重威胁。随着人们对仿刺参需求量的增加，迫切需要对其进行遗传改良，选育具有优良性状的新品种。遗传参数的估计是水产动物选择育种的一项基础工作，其中，遗传力是衡量育种进展和育种方法的关键性指标(Falconer et al, 2000)。

遗传力反映亲属间的相似程度和对人工与自然选择的反应速度，往往是育种工作的第一步。传统方法利用亲缘关系较近的个体表型间的相关系数及亲缘系数计算遗传力(Lstibůrek et al, 2018; 文超良等,2019)，在仿刺参遗传力估计中被广泛应用，例如，对仿刺参前期浮游阶段耳状幼体体长(栾生等, 2006)、耳状幼体和稚参体长(李云峰等, 2009)、幼参4 个生长性状(孟思远等, 2010)的遗传力估计等。传统方法利用系谱信息构建亲缘相关矩阵(Numerator relationship matrix, A 阵)，并利用线性混合模型估计遗传力，对系谱信息的准确性和完整性要求较高(李晶等,2020)；由于孟德尔抽样误差(Hill, 2014)，基于系谱推断的个体间亲缘关系准确性有限，并且存在系谱记录错误、部分群体系谱信息不完整或缺失的情况，使得传统法估计遗传力的准确性存在一定的局限。

随着测序技术的发展，越来越多的学者利用单核苷酸多态性(SNP)分子标记推断个体间的亲缘关系(Thompson, 1975; Lynch, 1988; Queller, 1989)，通过构建基因组亲缘关系矩阵(Genomic relationship matrix,GRM)(Li et al, 2014; 文超良等, 2019)，进行基于SNP分子标记的遗传力估计(SNP 遗传力)(Yang et al ,2017)，并在动植物以及人类群体中得到广泛应用(Ritland et al, 1996; Mousseau et al, 1998; Thomas et al, 2002; Visscher et al, 2006)。例如，Yang 等(2010)利用REML 分析方法估计了SNP 数据集解释的人类身高表型的变异程度；Guo 等(2018)利用GCTA 的REML 估计了栉孔扇贝(Chlamys fa rreri)的生长性状(扇贝长、宽、高和湿重)的SNP 遗传力；Fishback 等(2002)使用 REML 方法估计了虹鳟(Oncorhynchus mykiss)总长和重量的遗传力。Benjamin 等(2012)将这种基于GRM 的REML 分析方法称为GREML 法。GREML 法可以利用SNP 数据估计无亲缘关系的个体之间的遗传关系，然后，由全基因组SNP 解释的表型变异比例推断遗传力。与传统方法相比，REML 分析法适用于非均衡资料的方差组分估计，可大大提高方差组分的估计准确性。本研究利用已完成测序的仿刺参基因组(Li et al, 2018)，对不同养殖地理位置仿刺参群体进行全基因组重测序，并针对仿刺参重要经济性状疣足数量，进行全基因组水平的SNP 遗传力估计，为丰富仿刺参生长相关的分子育种理论、加快仿刺参优良品种的选育提供一定的理论依据。

1 材料与方法

1.1 仿刺参样本与表型

本研究215 个仿刺参样本分别来自于8 个不同的养殖地理位置：平岛、西小磨、山东栖霞口、蚆蛸岛、中俄杂交仿刺参、旅顺、黄龙尾和俄罗斯海参崴。首先，对所有仿刺参样本的疣足数量进行统计(表1)。

表1 仿刺参样本表型统计Tab.1 Statistics of papillae number for A. japonicus

1.2 全基因组重测序与参考SNP 选择

按照动物组织基因组DNA 提取试剂盒(天根，北京)说明书提取仿刺参基因组DNA。使用NanoDrop-2000 超微量紫外分光光度计进行定量分析。

检验合格的DNA 片段经末端修复、加polyA 尾、加测序接头、纯化和PCR 扩增完成整个文库制备。库检合格的文库通过Illumina Hiseq Xten 的PE150bp模式进行测序，每个样品测序深度为10×。对测序获得的测序数据进行质量过滤，使用 Cutadapt 软件(Martin, 2011)去除接头序列，使用SolexaQA 软件(Cox et al, 2010)去除质量值低于20 的碱基，将得到的高质量测序数据利用BWA 软件(Li et al, 2009)比对到仿刺参参考基因组。使用Samtools(Li et al, 2009)进行去重复，GATK(Ye et al, 2009)进行局部重比对，碱基质量值校正等处理以及 SNP 小片段插入缺失(Small INDEL)的检测，按照条件“QD＜2.0||MQ＜40.0||FS＞60.0||MQRankSum＜-12.5||ReadPosRankSum＜-8.0”对SNP 进行过滤，得到最终的SNP 位点集。

筛选SNP 是为了在重测序获得的高密度SNP 中减少或删除处于高度连锁不平衡的SNP，特别是在有限群体中有相当数量的SNP 之间存在高度的连锁不平衡关系，会降低对动物个体遗传参数的估计准确度。本研究采用 SNP 基因频率的平均欧式距离(Average Euclidean distance, AED) (Wu et al, 2016)来筛选SNPs，挑选出次要等位基因频率(Minor allele frequency, MAF)大于0.05 和0.1 的50K SNP，在筛选出来的 50K SNP 中随机均匀挑选不同密度的SNP(SNP 数量分别为5K、10K、15K、20K、25K、30K、35K、40K、45K 和50K)。

1.3 基于SNP 的遗传力估计

本研究使用 GCTA(Yang et al , 2011)(gcta_1.93.2beta 版 本) 的 GREML(Genomic relatedness matrix restricted maximum likelihood)中的期望最大约束似然法(EM-REML)、平均信息约束似然法(AI-REML)对仿刺参疣足数量的SNP 遗传力进行评估。EM 算法是由Dempster 等(1977)提出的一种迭代算法，可针对有缺失数据的数据资料，获得参数的最大似然估计值。缺失数据在实际的统计分析中非常普遍，EM 算法的基本思想是用缺失数据在给定参数的某个估计值的条件下的期望代替缺失数据，从而将不完全数据变成完全数据，使得似然函数的最大化变得相对简单。REML 的 AI 算法(平均信息算法)由Gilmour 等(1995)、Johnson 等(1995)和Jensen 等(1997)提出，它是将求似然函数最大值的 2 种常见的Newton-Raphson 算法和Fisher's scoring 方法结合起来的一种算法(张勤, 2007)。

GCTA 软件的核心即如下所示线性混合模型(Yang et al, 2011)：

式中，y 表示表型，b 表示固定效应协变量的系数，u 表示随机效应自变量的系数，这里的随机效应指的是所有SNP 位点对表型的效应，e 表示随机误差。

其中，u 和e 服从如下正态分布：I为一个n × n 的单位矩阵，n 表示样本量。

表型方差用如下公式表示：

W 为第ij 元素的标准化基因型矩阵，xij为第j 个个体的第i 个SNP 的参考等位基因的拷贝数，pi为参考等位基因频率。

GCTA 对于样本遗传相似度的定义公式：

定义所有SNP 位点的方差：

式中，N 为SNP 位点数，G 为个体间的亲缘关系矩阵(GRM)。

定义SNP 遗传力的计算公式：

对来自不同养殖地理位置的215 个仿刺参样本进行高通量重测序，获得一组高质量SNP。构建GRM矩阵后，使用 GCTA 的 2 种遗传参数估计方法(AI-REML 和EM-REML)，添加协变量的校正，针对仿刺参疣足数量，分别计算不同密度SNP 下的SNP遗传力估计值及染色体水平的SNP 遗传力估计值。

2 结果与分析

2.1 仿刺参疣足数量的表型数据分布

从215 个仿刺参样本疣足数量的频率直方图和概率密度曲线(图1a)可以看出，本研究不同地域的仿刺参疣足数量大致符合正态分布，表明仿刺参疣足数量性状不受地域环境的影响。考虑到图1a 中出现大量疣足数量大于60 的个体，且主要集中在中俄杂交仿刺参和俄罗斯海参崴群体中，因此，利用R MASS包中nlminb 函数，对215 个样本数据进行混合正态分布参数的最大似然估计(图1b)，结果显示，样本疣足数量符合2 个正态分布(均值39.8、方差6.79 和均值54.35、方差10.9)的最优拟合，说明样本中存在多疣足仿刺参群体情况。进一步把8 个不同地域的仿刺参群体分为正常群体(平岛、西小磨、山东栖霞口、蚆蛸岛、旅顺和黄龙尾)与多疣足群体(中俄杂交仿刺参和俄罗斯海参崴)，并利用R t.test 函数对上述2 个群体疣足数量进行t 检验，结果差异显著(P＜0.01)。根据上述表型数据分布结果，本研究将不同地域的仿刺参群体分为正常群体与多疣足群体，并在SNP 遗传力估计模型设计群体固定效应，以校群体分层对SNP 遗传力估计的影响。

图1 215 份仿刺参样本疣足数量的频率直方图和概率密度曲线(a)和双重高斯分布最优拟合(b)Fig.1 True frequency histogram and probability density curve (a) and optimal fitting of a double Gaussian distribution (b) of papillae number for the 215 A. japonicus samples

2.2 仿刺参疣足数量SNP 遗传力估计

本研究全基因组重测序质量较高，GC 分布正常，对基因组的覆盖度均超过90%，高质量SNPs 数量约为500 万个。使用GCTA 的2 种REML 分析方法(EM-REML 和AI-REML)对仿刺参疣足数量的遗传力进行评估，比较相同密度、相同次要等位基因频率筛选条件下，不同分析方法得到的仿刺参疣足数量SNP遗传力估计值。发现通过2 种方法得到的遗传力估计值差异极小(表2，表3)，证明本研究采用的分析方法有较好稳定性、可靠性。

除此之外，本研究对不同密度SNP 条件下的仿刺参疣足数量SNP 遗传力估计值进行比较，结果显示，当MAF＞0.05 时(表2)，不同密度SNP 条件下得到的仿刺参疣足数量SNP 遗传力估计均值为(0.566±0.022)～(0.612±0.003)，SNP 数量在5K、10K 时，得到的SNP 遗传力估计均值分别为0.567±0.022、0.588±0.011(AI-REML)和0.566±0.022、0.587±0.011(EMREML)；SNP 数量在15K～50K 时，得到的遗传力估计均值相较于低密度SNP 遗传力估计均值更大，且趋于稳定；SNP 数量为50K 时，得到的SNP 遗传力估计均值最高，分别为0.612±0.003(AI-REML)和0.611±0.003(EM-REML)。当MAF＞0.1 时(表3)，不同密度SNP 条件下得到的仿刺参疣足数量SNP 遗传力估计均值为(0.586±0.015)～(0.615±0.016)，SNP 数量在5K、10K、15K 时，得到的遗传力估计均值分别为0.593±0.026、0.596±0.017、0.587±0.014(AI-REML)和 0.592±0.026 、 0.595±0.017 、 0.586±0.015(EMREML)；SNP 数量在20K～50K 时，得到的遗传力估计均值相较于低密度SNP 遗传力估计均值更大，且趋于稳定；SNP 数量为30K 时，得到的SNP 遗传力估计均值最高，分别为0.615±0.016(AI-REML)和0.614±0.016 (EM-REML)，并且当SNP 数量在35K、40K、50K 时得到的遗传力估计均值均达到0.614。

表2 仿刺参疣足数量SNP 遗传力(MAF>0.05)Tab.2 SNP heritability estimates for papillae number in A. japonicus (MAF＞0.05)

表3 仿刺参疣足数量SNP 遗传力(MAF>0.1)Tab.3 SNP heritability estimates for papillae number in A. japonicus (MAF＞0.1)

2.3 仿刺参疣足数量染色体水平SNP 遗传力估计

根据已发表的仿刺参基因组和遗传连锁图谱，共发现22 个连锁群(Tian et al, 2015; Li et al, 2018)，本研究将Scaffold 拼接到染色体水平，并对仿刺参疣足数量进行染色体水平的SNP 遗传力估计，对染色体SNP 遗传力估计值和染色体长度进行回归分析(图2)，结果显示，7 号染色体的SNP 遗传力估计值最小，为0.085(MAF＞0.05)、0.094(MAF＞0.1)；2 号染色体的SNP 遗传力估计值最大，为0.598(MAF＞0.05)、0.599(MAF＞0.1)；14 号染色体长度最小，SNP遗传力估计值为 0.269(MAF＞0.05)、0.282(MAF＞0.1)；1 号染色体长度最大，SNP 遗传力估计值为0.522(MAF＞0.05)、0.526(MAF＞0.1)。在MAF＞0.05和MAF＞0.1 条件下，单个染色体的贡献和其长度显著相关，较长的染色体具有明显的线性趋势(P＜0.05)，且遗传力估计值最小的7 号染色体和长度最小的14 号染色体均偏离线性回归较远。

3 讨论

遗传力是反映性状遗传能力大小的重要遗传参数，准确合理的遗传参数估计有助于更好地理解遗传因素对特定群体某一性状的表型影响程度。估计动物个体的SNP 遗传力，首先需要筛选出一组适用的SNP。筛选这些SNP，一般需要满足两个条件：第一，筛选的 SNP 是该物种基因芯片中共同的SNP，这样可以很方便地将这些基因SNP 芯片用于分子育种，而不需要重新设计新的芯片或增加芯片中的SNP 位点；第二，选出的SNP 有较高的信息含量，其统计准确度要高。Hulsegge 等(2013)采用LD 的r2＞0.30 作为删除SNP 的尺度，结果表明，在保持相同准确性的前提下，使用这个尺度来筛选SNP，可以明显降低所需SNP 标记的数量。同时，需要筛选信息量高的SNP，筛选高信息量的SNP 可以依据不同的统计指标。例如，Hulsegge 等(2013)分别使用 Delta(群体间等位基因频率的绝对差异)、Wright(1978)的FST 以及Weir 等(1984)的FST 衡量标记信息量的效果。信息熵(Entropy)(Mitt et al, 2017)也是衡量SNP 信息量的指标。本研究采用SNP 基因频率的平均欧式距离(Average Euclidean distance,AED) (Wu et al, 2016)来筛选SNPs。后续工作将根据上述指标，挑选信息含量高并可进行准确的遗传参数估计的低密度SNP 标记。

图2 仿刺参疣足数量染色体水平SNP 遗传力估计Fig.2 Chromosome-wise heritability estimates for papillae number in A. japonicus

群体的分层效应会对遗传力估计值的准确性产生影响。本研究根据215 个仿刺参样本疣足数量的频率直方图和概率密度曲线(图1)初步将研究样本分为两大类，一类为来自中国的品种(正常群体)，另一类为俄罗斯以及中国与俄罗斯杂交的品种(多疣足群体)。计算时使用群体分层这一因素作为协变量进行校正，并得到了稳定结果，证明本研究使用的分组依据可靠。本研究使用了基于GCTA 的期望最大约束似然法(EM-REML)和平均信息约束似然法(AI-REML)对仿刺参疣足数量的SNP 遗传力进行估计，由估计值结果(表2、表3)可知，此方法得到的SNP 遗传力估计值均具有较高的稳定性和可信度。此前，孟思远等(2010)通过建立家系，使用传统方法估算仿刺参幼参阶段肉刺数目的遗传力，得到仿刺参肉刺(疣足)数量遗传力为0.191～0.404。本研究使用GCTA 的REML分析法得到了与传统方法相比更高、更稳定的遗传力估计值，说明疣足数量主要受加性效应的控制，对仿刺参的选择育种具有较大潜力。通过比较不同SNP密度下仿刺参疣足数量SNP 遗传力估计值，发现当SNP 数量在50K 时，得到的遗传力估计值均已达到稳定，说明50K SNP 的密度足够捕获数量性状基因座(QTL)大效应和小效应，证实了低密度SNP 用于准确估计遗传参数的可行性。同时，在全基因组选择方面，Wu 等(2016)通过设计牛的低密度SNP 芯片，对育种值估计的准确性进行评估，获得了较准确的结果，也证实了低密度SNP 芯片应用于分子育种的可靠性。由于全基因组分析对样本量要求较高，通常在500 以上(Yang et al, 2010)，本研究的样本量为215，统计量较少，容易导致统计结果出现一定程度的偏差，但所得结果与已报道数据(孟思远等, 2010; 和飞等, 2017)一致。后续工作将扩大样本统计量进行仿刺参数量性状的遗传力评估。

本研究在染色体水平对仿刺参疣足数量SNP 遗传力进行估计(图2)，提供性状遗传结构的信息。结果显示，仿刺参疣足数量是一个复杂的数量性状，不同染色体上的变异导致了仿刺参疣足数量的变异，与该性状相关的效应位点散布在各染色体上，说明了通过重测序的技术手段，在全基因组上寻找SNP 的重要性。单个染色体的贡献和其长度显著相关，表示在每条染色体上可能有多个小到中等效应的变异，而不是少数具有主要效应的变异。Goddard 等(2009)关于家畜复杂性状的分子标记辅助选择结果也支持这一假设，该报道称，性状变异有很大的遗传因素，但由于统计力量不足，大多数具有小到中度影响的变异尚未在遗传关联分析中被确认。而GCTA 可以克服分子标记辅助选择的缺陷，通过全基因组标记面板，所有QTL 都能和至少1 个marker 处于连锁不平衡状态。分析发现，7 号染色体虽不是最小染色体，但其遗传力估计值最低，与此相反，2 号染色体虽不是最大染色体，却具有最大遗传力估计值，说明较大的基因效应位点可能不在7 号染色体上，而位于2 号染色体，意味着QTL 出现在2 号染色体的概率较大，出现在7 号染色体上的概率较小。此结果的发现为将来开展仿刺参疣足数量的全基因组关联分析(GWAS)提供了支持。

4 总结

本研究对不同养殖地理位置仿刺参群体进行全基因组重测序，并针对仿刺参重要经济性状疣足数目，使用GCTA 软件的GREML 进行全基因组水平的SNP 遗传力估计，染色体水平SNP 遗传力估计，同时比较不同密度SNP 标记条件下的仿刺参疣足数量SNP 遗传力估计值；结果表明，仿刺参疣足数量SNP遗传力为中等遗传力(MAF＞0.05 时，0.566～0.612；MAF＞0.1 时，0.586～0.615)，当SNP 数量在50K 时，得到的遗传力估计值均已达到稳定，且染色体SNP遗传力与其长度显著相关；说明仿刺参疣足数量是复杂数量性状，由全基因组水平多基因共同作用，且50K SNP 密度足够捕获QTL 大效应和小效应。本研究通过分子学研究方法得到可靠的遗传力估计值，有效解决了传统方法系谱记录繁琐，遗传力估计准确性差等问题，为分子育种技术提供了一定理论依据；同时，证明低密度SNP 标记用于仿刺参分子育种的可能性，可为低密度SNP 芯片的设计开发提供支持。