王卓标 ，方 铭，郑乐云，葛 辉，陈欣欣，罗辉玉，王志勇*
(1.集美大学水产学院，农业农村部东海海水健康与养殖重点实验室，福建 厦门 361021；
2.福建省水产研究所，福建省海洋生物增养殖与高值化利用重点实验室，福建 厦门 361013)
赤点石斑鱼(Epinephelusakaara)是一种高经济价值的海洋鱼类，主要分布于西太平洋地区包括中国东海、南海以及日本和韩国的南部海域[1]。20世纪60年代日本学者就研究了赤点石斑鱼的产卵习性和早期生活史[2]，我国在20世纪80年代取得赤点石斑鱼的人工育苗成功[3-4]，随后在福建、广东等沿海地区开始了网箱养殖。然而多年来赤点石斑鱼的苗种生产受到神经坏死病(Viral nervous necrosis，VNN)的严重困扰，导致育苗的成活率与成功率都很低，VNN成为制约赤点石斑鱼人工养殖业发展的一个重要因素。
鱼类病毒性神经坏死病是由神经坏死病毒(Nervous necrosis virus，NNV)引起的一种全球范围的鱼类流行性传染病。NNV为海水鱼中较常见、危害严重的传染病之一，至今已报告的受害鱼类有鳗鲡目(Anguilliformes)、鲈形目(Perciformes)、鲽形目(Pleuronectiformes)、鲀形目(Tetraodontiformes)、鳕形目(Gadiformes)中的40余种[5-6]，其中被感染的种类集中在石斑鱼、鲈鱼，在中国受影响最大的是赤点石斑鱼和斜带石斑鱼。对石斑鱼而言，神经坏死病的高发期为鱼类的幼鱼阶段，孵化后的1～3周为发病高峰期，严重时发病率可达100%[7]，幼鱼成活率低于10%。神经坏死病的病症通常表现为鱼苗在水中以螺旋状旋转为主的异常游动，伴随着食欲减退，静止时腹部向上。组织学检测发现细胞空泡化主要集中发生于中枢神经系统细胞以及视网膜[8]。受感染的幼鱼绝大多数在短期内死亡，因此常常导致人工育苗失败，或成活率极低。近年来，随着养殖密度的不断提高和受感染鱼类种类增加，其危害程度愈发严重。
选育抗病品种是解决养殖鱼类病害问题的一个有效途径[9-10]。但是传统的选育方法进展慢，效果较差。2001年Meuwissen T H E等[11]提出了基因组选择(Genomic selection，GS)的方法，该方法具有不需构建家系、育种值估计的准确性高、育种效率高，并可以有效控制近交等多方面优点。如今，随着高通量DNA测序成本不断降低，已有越来越多的水产动物育种开始使用基因组选择的方法[12-18]。目前已经有一些研究者将基因组选择应用于鱼类抗病育种研究，如Tsai H Y等[19]报道了对大西洋鲑抗海虱、Liu Y等[20]报道了对牙鲆抗爱德华氏菌的基因组选择研究。Palaiokostas C等[21]对欧洲鲈鱼(Dicentrarchuslabrax)的神经坏死病毒病抗性进行评估，发现使用基因组选择的方法与比系谱选育预测能力增加13%。福建省水产研究所石斑鱼研究团队已经完成了赤点石斑鱼全基因组测序组装[22]，华南农业大学Yang M等[23]通过对100尾赤点石斑鱼神经坏死病毒(Red-spotted grouper nervous necrosis virus，RGNNV)易感与抗性石斑鱼进行全基因组关联分析，找到了一些抗病相关的单核酸多态性(Single nucleotide polymorphisms,SNPs)位点；
但迄今还没有见到对赤点石斑鱼抗神经坏死病基因组选择研究的报道。本研究对460尾赤点石斑鱼神经坏死病易感(染病死亡，230尾)和抗性(最终健康存活，230尾)鱼苗，通过基因组重测序获得高密度的SNPs集，进行抗病遗传力评估和基因组选择预测，以期为后续的抗病育种实践提供必要的理论参考。
1.1 样本的采集与鉴定
实验材料采集于厦门刘五店，该育苗场在2019年赤点石斑鱼育苗生产中遭遇了神经坏死病，分别在发病时期采集已染病的濒死个体作为易感组(共230尾)，渡过发病期后，采集同样数量存活的健康个体作为抗病组。对460尾鱼苗采集全鱼，保存于95%乙醇中，用于DNA提取。发病鱼苗病症除了根据其病状判断外，还提取30尾病鱼头部组织DNA制成10个混合样品，用RGNNV特异性检测引物(正向引物：5’-CACCGCTTTGCAATCACAATG-3’，反向引物：5’-GTCATCAACGATACGCACTAGG-3’)进行了PCR扩增验证(结果均为阳性)。
1.2 基因分型和质量控制
使用南京诺唯赞生物科技股份有限公司的快速组织基因组DNA试剂盒提取鳍条组织基因组DNA，进行质检和建库后，在Illumina NovaSeq 6000平台(Illumina，USA)进行WGS测序。其中450尾的目标测序深度为4×，另随机挑选10尾测序深度为20×，用于提供高质量的SNPs变异参考。首先使用fastQC(https：//www.bioinformatics.babraham.ac.uk/projects/fastqc)对测序数据进行质量检测，后使用fastp[24]对测序数据进行过滤。最后使用MultiQC[25]对最终的质控结果进行汇总检查。
通过BWA-MEM[26]将clean reads 比对到赤点石斑鱼基因组[22]上，将产生的文件利用Samtools 进行排序并转化为bam文件格式。之后使用sambamba[27]对bam中构建文库时的PCR重复进行标记，对标记重复后的bam文件使用GATK[28]中的“HaplotypeCaller”的模块进行单个样本变异的检测，再通过“CombineGVCFs”将单个的变异集整合为群体的VCF(Variant call format)文件。使用“SelectVariants”“VariantFiltration”模块进行硬过滤以及双等位基因的提取，并使用BCFtools[29]对缺失率大于20%的变异进行过滤，然后对缺失的基因型使用Beagle[30]进行填充。最后通过PLINK[31]对VCF进行过滤以及格式转换，过滤标准为：1)次等位基因频率MAF>0.05；
2)HWE<1e-6。最终获得460个样本的高质量SNPs数据用于后续分析。
1.3 群体结构分析
利用PLINK的PCA模块进行主成分分析(Principal component analysis，PCA)。根据PCA结果所占比重前2个主成分PC1～PC2的数据，通过在线绘图工具bioinformatics(http：//www.bioinformatics.com.cn)绘制PCA图。
1.4 遗传参数估计
使用GCTA[32]和R语言的EMMREML包(https：//CRAN.R-project.org/package=EMMREML)对遗传参数进行估计以及后续基因组育种值(Genomic estimated breeding value,GEBV)的计算。基因组最佳线性无偏预测(Genomic best linear unbiased prediction，GBLUP)模型公式为：
y=Xβ+Zμ+e
(1)
其中y代表的是性状，Z是根据SNPs效应而设计的矩阵(基因型“AA”、“Aa”和“aa”分别对应着SNPs基因型中的“0”、“1”和“2”)，β是固定效应(地点效应)，混合线性模型也可以展开表示为：
(3)
其中B代表的是等位基因的m×n共享矩阵，m代表总的标记数量，n代表样本的总数；
P代表的是每个SNPs位点次等位基因的频率(Minor allele frequency，MAF)组成的矩阵；
pj为每个SNPs在第j个位点等位基因“a”的频率。本研究的G矩阵通过460尾的SNPs位点的基因型数据，使用GCTA[32]的计算得出，相较于普通的BLUP模型，将系谱推测的个体遗传关系的A矩阵改为由全基因组测序得到的SNPs标记构建的G矩阵，GEBV估算的准确度更高[12]。
1.5 交叉验证
交叉验证用于测试预测的准确度，本研究采用5折的交叉验证，将460尾个体随机分为5组(参考群体与测试群体比例为4∶1)，使用GBLUP对460个体的GEBV进行计算，测试群体的表型设置为空，通过参考群体的表型对测试群体的表型进行预测，比较两个值之间的相关系数来衡量模型的预测能力。重复以上步骤5次，使得每一组均有一次机会作为候选群体。交叉验证的重复次数为100次。
1.6 不同标记密度对预测能力的影响
为考察不同SNPs数量GEBV的预测力，除了使用全部的SNPs(6 132 865个SNPs，6 132 k)外，还设计了13个不同个数的标记子集，分别为0.5、1、3、5、10、30、50、100、250、500、1 000、2 500 k。为了降低抽样的误差，使用GATK SelectVariants包对每个数量的标记集都进行50次的随机抽样后，对每次抽样的结果都进行100次的5折交叉验证。
1.7 不同覆盖度标准对预测的影响
根据本研究使用原始的测序深度，对每个个体在每个SNPs上的覆盖度设置了6个过滤标准，即DP3、DP4、DP5、DP6、DP8和DP10，分别对应于3×、4×、5×、6×、8×与10×的reads覆盖度。使用过滤后的VCF文件，通过BCFtools对于每个SNPs位点大于80%的个体的基因型覆盖度大于过滤标准，就保留该位点。过滤后使用BEAGLE对VCF文件进行填充，后对每个过滤标准产生的SNPs标记集进行100次5折交叉验证。
2.1 标记分型及群体结构分析
本研究对460尾赤点石斑鱼苗进行基因组重测序，共获得2.67 Tb clean data数据挖掘SNPs，通过哈代-温伯格平衡(HWE)>10-6以及次要等位基因频率(MAF)>5%经过质控后，最终用于分析的SNPs共有5 412 683个(未进行覆盖度过滤)，平均标记密度约为200.1 bp/SNPs。对460个样品进行主成分分析，提取前2个主成分(Principal components，PCs)的结果(图1)，可以看出群体存在明显的分层现象，分成了2个小的聚类群，但易感(Case)和抗性(Control)个体在两个亚群中均匀分布。
注：Case组为易感组；
Control组为抗性组。
2.2 遗传参数估计及基因组选择预测
利用全部460尾鱼苗的表型(230尾易感、230尾抗性)和全部SNPs位点的基因型数据计算抗神经坏死病性状的遗传力，使用R包EMMREML的结果为0.566 2，GCTA的AI-REML经过1 000迭代的结果为0.566 6，两种方法估算结果相近。对460尾赤点石斑鱼进行抗神经坏死病基因组选择的预测力分析，80%个体作参考群，20%个体作验证群，分别进行了100次随机抽样的交叉验证，平均的基因组预测准确度为(0.359±0.019)。
2.3 不同SNPs数量的基因组选择预测力
从全部分型位点(5 412 683个)中随机取不同数量SNPs进行基因组选择分析，每组标记数量随机抽样50次。随机抽取的SNPs标记数量与育种值估计准确度的关系见图2。由图2可知，当标记密度由0.5 k升到5 k时基因组的预测能力迅速增加，标记密度达到5 k以上时，育种值估计准确度随着标记的增加趋于平稳。
2.4 不同覆盖度标准对预测的影响
对测序数据进行不同覆盖度的过滤，获得的SNPs数目如表1所示。随着覆盖度过滤标准提高，保留下来的SNPs标记数量急剧减少。用不同覆盖度过滤得到的SNPs分别进行基因组选择预测能力评估，结果如图3所示，当覆盖度标准为DP3(即≥3×)时，基因组选择预测力最高；
覆盖度标准为DP4和DP5(保留的SNPs数量分别为5 982和11 291)时，基因组选择预测力几乎完全一样，均为0.346。覆盖度标准提高为DP6时，保留的SNPs只有3 695个，预测力下降到0.30左右，DP8和DP10过滤剩余的SNPs数分别只有1 817和830个，预测力下降为0.26左右。
表1 不同过滤标准剩余的SNPs数
基因组选择的首要工作建立参考群，利用参考群估算各个分子标记位点对性状表型的效应。本文建立了460尾石斑鱼的参考群，主成分分析显示群体聚集为2个不同的亚群，表明群体间存在较分化的亲缘关系。赤点石斑鱼人工育苗中亲本雌雄配比一般达到(20～30)∶1，使用的雄性亲本数量非常有限，因此同一批、尤其是同一池的受精卵可能来自少数雄性亲本。本研究发现鱼苗亲缘关系分化成2个亚群，推测可能采样的鱼苗来自于2个或几个雄性亲本及为数不多的母本。
候选群与参考群之间的亲缘关系是影响基因组选择效果的关键因素[34-35]。在应用基因组选择技术时，候选群亲鱼要尽可能地与参考群中个体存在亲缘关系，如果是在同一家育苗场持续选育，并在鱼苗繁育中有可追踪的生产记录，则可以较可靠地推断候选亲鱼与参考群之间的亲缘关系，这将更有利于基因组选择技术的应用。另外，如果经费允许，有必要进一步扩大参考群规模，使参考群包含更多家系，增加参考群的多样性，对于扩展所建立的基因组选择技术的适用面、提高基因组选择的预测准确性和选育效率都具有重要意义[34]。
此外，基因组选择准确性还受基因组的标记密度的影响[36]。本文通过随机抽样研究了不同标记密度对育种值预测准确性的影响，当使用5 k个标记时，即可使预测准确性接近利用全基因组SNPs的水平，其后随着标记个数增加，预测准确性变化趋于平缓；
当标记个数为50 k时，预测准确性与利用全基因组所有SNPs的预测准确性几乎一致，这与Tsai H Y等[19]对大西洋鲑抗海虱性状基因组预测的研究结果非常相似，提示如果设计赤点石斑鱼育种芯片，50 k的标记密度可满足育种要求。尽管如此，利用全基因组SNPs标记有利于挖掘具有较大效应的分子标记及因果突变，将这些大效应分子标记和因果突变嵌入预测模型，建立基于主效-微效多基因效应相结合的预测模型，能进一步提高育种值预测准确性。利用同一批数据，笔者在基因组上已经发现了两个效应值极强的GWAS信号(结果未列出)，计划在信号内部挖掘可靠的分子标记或因果突变，并将其作为协变量加入GBLUP模型中，通过主效标记/因果突变进一步吸收残余误差的方式进一步提高预测准确性。对测序数据进行覆盖度过滤能够提高对挖掘的标记位点基因分型的准确性，从而在一定程度上提高育种值估算的准确性(图3)。但是如表1和图3所示，在每个个体测序量不变的情况下，随着覆盖度过滤标准提高，可保留用于分析的标记数量急剧减少，当DP≥4时，尽管标记分型的准确度会明显提升，但是由于保留的标记数量减少，育种值预测准确性已低于不进行覆盖度过滤。据Zhang W等[37]对大黄鱼研究的结果，采用全基因组重测序，由于可挖掘到的标记数多，即使测序覆盖度低至0.5×，而且不进行覆盖度过滤，也能获得与8×覆盖度测序基本一致的基因组选择效果，这将大大降低候选亲本标记分型费用。因此，可以预期，随着高通量DNA测序价格进一步降低，基因组重测序将越来越多被用于基因组选择，也将成为石斑鱼全基因组选择育种的主要工具。