赤點(diǎn)石斑魚(yú)抗神經(jīng)壞死病的基因組選擇評(píng)估

王卓標(biāo) ，方 銘，鄭樂(lè)云，葛 輝，陳欣欣，羅輝玉，王志勇*

(1.集美大學(xué)水產(chǎn)學(xué)院，農(nóng)業(yè)農(nóng)村部東海海水健康與養(yǎng)殖重點(diǎn)實(shí)驗(yàn)室，福建 廈門(mén) 361021；2.福建省水產(chǎn)研究所，福建省海洋生物增養(yǎng)殖與高值化利用重點(diǎn)實(shí)驗(yàn)室，福建 廈門(mén) 361013)

赤點(diǎn)石斑魚(yú)(Epinephelusakaara)是一種高經(jīng)濟(jì)價(jià)值的海洋魚(yú)類，主要分布于西太平洋地區(qū)包括中國(guó)東海、南海以及日本和韓國(guó)的南部海域[1]。20世紀(jì)60年代日本學(xué)者就研究了赤點(diǎn)石斑魚(yú)的產(chǎn)卵習(xí)性和早期生活史[2]，我國(guó)在20世紀(jì)80年代取得赤點(diǎn)石斑魚(yú)的人工育苗成功[3-4]，隨后在福建、廣東等沿海地區(qū)開(kāi)始了網(wǎng)箱養(yǎng)殖。然而多年來(lái)赤點(diǎn)石斑魚(yú)的苗種生產(chǎn)受到神經(jīng)壞死病(Viral nervous necrosis，VNN)的嚴(yán)重困擾，導(dǎo)致育苗的成活率與成功率都很低，VNN成為制約赤點(diǎn)石斑魚(yú)人工養(yǎng)殖業(yè)發(fā)展的一個(gè)重要因素。

魚(yú)類病毒性神經(jīng)壞死病是由神經(jīng)壞死病毒(Nervous necrosis virus，NNV)引起的一種全球范圍的魚(yú)類流行性傳染病。NNV為海水魚(yú)中較常見(jiàn)、危害嚴(yán)重的傳染病之一，至今已報(bào)告的受害魚(yú)類有鰻鱺目(Anguilliformes)、鱸形目(Perciformes)、鰈形目(Pleuronectiformes)、鲀形目(Tetraodontiformes)、鱈形目(Gadiformes)中的40余種[5-6]，其中被感染的種類集中在石斑魚(yú)、鱸魚(yú)，在中國(guó)受影響最大的是赤點(diǎn)石斑魚(yú)和斜帶石斑魚(yú)。對(duì)石斑魚(yú)而言，神經(jīng)壞死病的高發(fā)期為魚(yú)類的幼魚(yú)階段，孵化后的1～3周為發(fā)病高峰期，嚴(yán)重時(shí)發(fā)病率可達(dá)100%[7]，幼魚(yú)成活率低于10%。神經(jīng)壞死病的病癥通常表現(xiàn)為魚(yú)苗在水中以螺旋狀旋轉(zhuǎn)為主的異常游動(dòng)，伴隨著食欲減退，靜止時(shí)腹部向上。組織學(xué)檢測(cè)發(fā)現(xiàn)細(xì)胞空泡化主要集中發(fā)生于中樞神經(jīng)系統(tǒng)細(xì)胞以及視網(wǎng)膜[8]。受感染的幼魚(yú)絕大多數(shù)在短期內(nèi)死亡，因此常常導(dǎo)致人工育苗失敗，或成活率極低。近年來(lái)，隨著養(yǎng)殖密度的不斷提高和受感染魚(yú)類種類增加，其危害程度愈發(fā)嚴(yán)重。

選育抗病品種是解決養(yǎng)殖魚(yú)類病害問(wèn)題的一個(gè)有效途徑[9-10]。但是傳統(tǒng)的選育方法進(jìn)展慢，效果較差。2001年Meuwissen T H E等[11]提出了基因組選擇(Genomic selection，GS)的方法，該方法具有不需構(gòu)建家系、育種值估計(jì)的準(zhǔn)確性高、育種效率高，并可以有效控制近交等多方面優(yōu)點(diǎn)。如今，隨著高通量DNA測(cè)序成本不斷降低，已有越來(lái)越多的水產(chǎn)動(dòng)物育種開(kāi)始使用基因組選擇的方法[12-18]。目前已經(jīng)有一些研究者將基因組選擇應(yīng)用于魚(yú)類抗病育種研究，如Tsai H Y等[19]報(bào)道了對(duì)大西洋鮭抗海虱、Liu Y等[20]報(bào)道了對(duì)牙鲆抗愛(ài)德華氏菌的基因組選擇研究。Palaiokostas C等[21]對(duì)歐洲鱸魚(yú)(Dicentrarchuslabrax)的神經(jīng)壞死病毒病抗性進(jìn)行評(píng)估，發(fā)現(xiàn)使用基因組選擇的方法與比系譜選育預(yù)測(cè)能力增加13%。福建省水產(chǎn)研究所石斑魚(yú)研究團(tuán)隊(duì)已經(jīng)完成了赤點(diǎn)石斑魚(yú)全基因組測(cè)序組裝[22]，華南農(nóng)業(yè)大學(xué)Yang M等[23]通過(guò)對(duì)100尾赤點(diǎn)石斑魚(yú)神經(jīng)壞死病毒(Red-spotted grouper nervous necrosis virus，RGNNV)易感與抗性石斑魚(yú)進(jìn)行全基因組關(guān)聯(lián)分析，找到了一些抗病相關(guān)的單核酸多態(tài)性(Single nucleotide polymorphisms,SNPs)位點(diǎn)；但迄今還沒(méi)有見(jiàn)到對(duì)赤點(diǎn)石斑魚(yú)抗神經(jīng)壞死病基因組選擇研究的報(bào)道。本研究對(duì)460尾赤點(diǎn)石斑魚(yú)神經(jīng)壞死病易感(染病死亡，230尾)和抗性(最終健康存活，230尾)魚(yú)苗，通過(guò)基因組重測(cè)序獲得高密度的SNPs集，進(jìn)行抗病遺傳力評(píng)估和基因組選擇預(yù)測(cè)，以期為后續(xù)的抗病育種實(shí)踐提供必要的理論參考。

1 材料與方法

1.1 樣本的采集與鑒定

實(shí)驗(yàn)材料采集于廈門(mén)劉五店，該育苗場(chǎng)在2019年赤點(diǎn)石斑魚(yú)育苗生產(chǎn)中遭遇了神經(jīng)壞死病，分別在發(fā)病時(shí)期采集已染病的瀕死個(gè)體作為易感組(共230尾)，渡過(guò)發(fā)病期后，采集同樣數(shù)量存活的健康個(gè)體作為抗病組。對(duì)460尾魚(yú)苗采集全魚(yú)，保存于95%乙醇中，用于DNA提取。發(fā)病魚(yú)苗病癥除了根據(jù)其病狀判斷外，還提取30尾病魚(yú)頭部組織DNA制成10個(gè)混合樣品，用RGNNV特異性檢測(cè)引物(正向引物：5’-CACCGCTTTGCAATCACAATG-3’，反向引物：5’-GTCATCAACGATACGCACTAGG-3’)進(jìn)行了PCR擴(kuò)增驗(yàn)證(結(jié)果均為陽(yáng)性)。

1.2 基因分型和質(zhì)量控制

使用南京諾唯贊生物科技股份有限公司的快速組織基因組DNA試劑盒提取鰭條組織基因組DNA，進(jìn)行質(zhì)檢和建庫(kù)后，在Illumina NovaSeq 6000平臺(tái)(Illumina，USA)進(jìn)行WGS測(cè)序。其中450尾的目標(biāo)測(cè)序深度為4×，另隨機(jī)挑選10尾測(cè)序深度為20×，用于提供高質(zhì)量的SNPs變異參考。首先使用fastQC(https：//www.bioinformatics.babraham.ac.uk/projects/fastqc)對(duì)測(cè)序數(shù)據(jù)進(jìn)行質(zhì)量檢測(cè)，后使用fastp[24]對(duì)測(cè)序數(shù)據(jù)進(jìn)行過(guò)濾。最后使用MultiQC[25]對(duì)最終的質(zhì)控結(jié)果進(jìn)行匯總檢查。

通過(guò)BWA-MEM[26]將clean reads 比對(duì)到赤點(diǎn)石斑魚(yú)基因組[22]上，將產(chǎn)生的文件利用Samtools 進(jìn)行排序并轉(zhuǎn)化為bam文件格式。之后使用sambamba[27]對(duì)bam中構(gòu)建文庫(kù)時(shí)的PCR重復(fù)進(jìn)行標(biāo)記，對(duì)標(biāo)記重復(fù)后的bam文件使用GATK[28]中的“HaplotypeCaller”的模塊進(jìn)行單個(gè)樣本變異的檢測(cè)，再通過(guò)“CombineGVCFs”將單個(gè)的變異集整合為群體的VCF(Variant call format)文件。使用“SelectVariants”“VariantFiltration”模塊進(jìn)行硬過(guò)濾以及雙等位基因的提取，并使用BCFtools[29]對(duì)缺失率大于20%的變異進(jìn)行過(guò)濾，然后對(duì)缺失的基因型使用Beagle[30]進(jìn)行填充。最后通過(guò)PLINK[31]對(duì)VCF進(jìn)行過(guò)濾以及格式轉(zhuǎn)換，過(guò)濾標(biāo)準(zhǔn)為：1)次等位基因頻率MAF>0.05；2)HWE<1e-6。最終獲得460個(gè)樣本的高質(zhì)量SNPs數(shù)據(jù)用于后續(xù)分析。

1.3 群體結(jié)構(gòu)分析

利用PLINK的PCA模塊進(jìn)行主成分分析(Principal component analysis，PCA)。根據(jù)PCA結(jié)果所占比重前2個(gè)主成分PC1～PC2的數(shù)據(jù)，通過(guò)在線繪圖工具bioinformatics(http：//www.bioinformatics.com.cn)繪制PCA圖。

1.4 遺傳參數(shù)估計(jì)

使用GCTA[32]和R語(yǔ)言的EMMREML包(https：//CRAN.R-project.org/package=EMMREML)對(duì)遺傳參數(shù)進(jìn)行估計(jì)以及后續(xù)基因組育種值(Genomic estimated breeding value,GEBV)的計(jì)算?；蚪M最佳線性無(wú)偏預(yù)測(cè)(Genomic best linear unbiased prediction，GBLUP)模型公式為：

y=Xβ+Zμ+e

(1)

其中y代表的是性狀，Z是根據(jù)SNPs效應(yīng)而設(shè)計(jì)的矩陣(基因型“AA”、“Aa”和“aa”分別對(duì)應(yīng)著SNPs基因型中的“0”、“1”和“2”)，β是固定效應(yīng)(地點(diǎn)效應(yīng))，混合線性模型也可以展開(kāi)表示為：

(3)

其中B代表的是等位基因的m×n共享矩陣，m代表總的標(biāo)記數(shù)量，n代表樣本的總數(shù)；P代表的是每個(gè)SNPs位點(diǎn)次等位基因的頻率(Minor allele frequency，MAF)組成的矩陣；pj為每個(gè)SNPs在第j個(gè)位點(diǎn)等位基因“a”的頻率。本研究的G矩陣通過(guò)460尾的SNPs位點(diǎn)的基因型數(shù)據(jù)，使用GCTA[32]的計(jì)算得出，相較于普通的BLUP模型，將系譜推測(cè)的個(gè)體遺傳關(guān)系的A矩陣改為由全基因組測(cè)序得到的SNPs標(biāo)記構(gòu)建的G矩陣，GEBV估算的準(zhǔn)確度更高[12]。

1.5 交叉驗(yàn)證

交叉驗(yàn)證用于測(cè)試預(yù)測(cè)的準(zhǔn)確度，本研究采用5折的交叉驗(yàn)證，將460尾個(gè)體隨機(jī)分為5組(參考群體與測(cè)試群體比例為4∶1)，使用GBLUP對(duì)460個(gè)體的GEBV進(jìn)行計(jì)算，測(cè)試群體的表型設(shè)置為空，通過(guò)參考群體的表型對(duì)測(cè)試群體的表型進(jìn)行預(yù)測(cè)，比較兩個(gè)值之間的相關(guān)系數(shù)來(lái)衡量模型的預(yù)測(cè)能力。重復(fù)以上步驟5次，使得每一組均有一次機(jī)會(huì)作為候選群體。交叉驗(yàn)證的重復(fù)次數(shù)為100次。

1.6 不同標(biāo)記密度對(duì)預(yù)測(cè)能力的影響

為考察不同SNPs數(shù)量GEBV的預(yù)測(cè)力，除了使用全部的SNPs(6 132 865個(gè)SNPs，6 132 k)外，還設(shè)計(jì)了13個(gè)不同個(gè)數(shù)的標(biāo)記子集，分別為0.5、1、3、5、10、30、50、100、250、500、1 000、2 500 k。為了降低抽樣的誤差，使用GATK SelectVariants包對(duì)每個(gè)數(shù)量的標(biāo)記集都進(jìn)行50次的隨機(jī)抽樣后，對(duì)每次抽樣的結(jié)果都進(jìn)行100次的5折交叉驗(yàn)證。

1.7 不同覆蓋度標(biāo)準(zhǔn)對(duì)預(yù)測(cè)的影響

根據(jù)本研究使用原始的測(cè)序深度，對(duì)每個(gè)個(gè)體在每個(gè)SNPs上的覆蓋度設(shè)置了6個(gè)過(guò)濾標(biāo)準(zhǔn)，即DP3、DP4、DP5、DP6、DP8和DP10，分別對(duì)應(yīng)于3×、4×、5×、6×、8×與10×的reads覆蓋度。使用過(guò)濾后的VCF文件，通過(guò)BCFtools對(duì)于每個(gè)SNPs位點(diǎn)大于80%的個(gè)體的基因型覆蓋度大于過(guò)濾標(biāo)準(zhǔn)，就保留該位點(diǎn)。過(guò)濾后使用BEAGLE對(duì)VCF文件進(jìn)行填充，后對(duì)每個(gè)過(guò)濾標(biāo)準(zhǔn)產(chǎn)生的SNPs標(biāo)記集進(jìn)行100次5折交叉驗(yàn)證。

2 結(jié)果

2.1 標(biāo)記分型及群體結(jié)構(gòu)分析

本研究對(duì)460尾赤點(diǎn)石斑魚(yú)苗進(jìn)行基因組重測(cè)序，共獲得2.67 Tb clean data數(shù)據(jù)挖掘SNPs，通過(guò)哈代-溫伯格平衡(HWE)>10-6以及次要等位基因頻率(MAF)>5%經(jīng)過(guò)質(zhì)控后，最終用于分析的SNPs共有5 412 683個(gè)(未進(jìn)行覆蓋度過(guò)濾)，平均標(biāo)記密度約為200.1 bp/SNPs。對(duì)460個(gè)樣品進(jìn)行主成分分析，提取前2個(gè)主成分(Principal components，PCs)的結(jié)果(圖1)，可以看出群體存在明顯的分層現(xiàn)象，分成了2個(gè)小的聚類群，但易感(Case)和抗性(Control)個(gè)體在兩個(gè)亞群中均勻分布。

注：Case組為易感組；Control組為抗性組。

2.2 遺傳參數(shù)估計(jì)及基因組選擇預(yù)測(cè)

利用全部460尾魚(yú)苗的表型(230尾易感、230尾抗性)和全部SNPs位點(diǎn)的基因型數(shù)據(jù)計(jì)算抗神經(jīng)壞死病性狀的遺傳力，使用R包EMMREML的結(jié)果為0.566 2，GCTA的AI-REML經(jīng)過(guò)1 000迭代的結(jié)果為0.566 6，兩種方法估算結(jié)果相近。對(duì)460尾赤點(diǎn)石斑魚(yú)進(jìn)行抗神經(jīng)壞死病基因組選擇的預(yù)測(cè)力分析，80%個(gè)體作參考群，20%個(gè)體作驗(yàn)證群，分別進(jìn)行了100次隨機(jī)抽樣的交叉驗(yàn)證，平均的基因組預(yù)測(cè)準(zhǔn)確度為(0.359±0.019)。

2.3 不同SNPs數(shù)量的基因組選擇預(yù)測(cè)力

從全部分型位點(diǎn)(5 412 683個(gè))中隨機(jī)取不同數(shù)量SNPs進(jìn)行基因組選擇分析，每組標(biāo)記數(shù)量隨機(jī)抽樣50次。隨機(jī)抽取的SNPs標(biāo)記數(shù)量與育種值估計(jì)準(zhǔn)確度的關(guān)系見(jiàn)圖2。由圖2可知，當(dāng)標(biāo)記密度由0.5 k升到5 k時(shí)基因組的預(yù)測(cè)能力迅速增加，標(biāo)記密度達(dá)到5 k以上時(shí)，育種值估計(jì)準(zhǔn)確度隨著標(biāo)記的增加趨于平穩(wěn)。

2.4 不同覆蓋度標(biāo)準(zhǔn)對(duì)預(yù)測(cè)的影響

對(duì)測(cè)序數(shù)據(jù)進(jìn)行不同覆蓋度的過(guò)濾，獲得的SNPs數(shù)目如表1所示。隨著覆蓋度過(guò)濾標(biāo)準(zhǔn)提高，保留下來(lái)的SNPs標(biāo)記數(shù)量急劇減少。用不同覆蓋度過(guò)濾得到的SNPs分別進(jìn)行基因組選擇預(yù)測(cè)能力評(píng)估，結(jié)果如圖3所示，當(dāng)覆蓋度標(biāo)準(zhǔn)為DP3(即≥3×)時(shí)，基因組選擇預(yù)測(cè)力最高；覆蓋度標(biāo)準(zhǔn)為DP4和DP5(保留的SNPs數(shù)量分別為5 982和11 291)時(shí)，基因組選擇預(yù)測(cè)力幾乎完全一樣，均為0.346。覆蓋度標(biāo)準(zhǔn)提高為DP6時(shí)，保留的SNPs只有3 695個(gè)，預(yù)測(cè)力下降到0.30左右，DP8和DP10過(guò)濾剩余的SNPs數(shù)分別只有1 817和830個(gè)，預(yù)測(cè)力下降為0.26左右。

表1 不同過(guò)濾標(biāo)準(zhǔn)剩余的SNPs數(shù)

3 討論

基因組選擇的首要工作建立參考群，利用參考群估算各個(gè)分子標(biāo)記位點(diǎn)對(duì)性狀表型的效應(yīng)。本文建立了460尾石斑魚(yú)的參考群，主成分分析顯示群體聚集為2個(gè)不同的亞群，表明群體間存在較分化的親緣關(guān)系。赤點(diǎn)石斑魚(yú)人工育苗中親本雌雄配比一般達(dá)到(20～30)∶1，使用的雄性親本數(shù)量非常有限，因此同一批、尤其是同一池的受精卵可能來(lái)自少數(shù)雄性親本。本研究發(fā)現(xiàn)魚(yú)苗親緣關(guān)系分化成2個(gè)亞群，推測(cè)可能采樣的魚(yú)苗來(lái)自于2個(gè)或幾個(gè)雄性親本及為數(shù)不多的母本。

候選群與參考群之間的親緣關(guān)系是影響基因組選擇效果的關(guān)鍵因素[34-35]。在應(yīng)用基因組選擇技術(shù)時(shí)，候選群親魚(yú)要盡可能地與參考群中個(gè)體存在親緣關(guān)系，如果是在同一家育苗場(chǎng)持續(xù)選育，并在魚(yú)苗繁育中有可追蹤的生產(chǎn)記錄，則可以較可靠地推斷候選親魚(yú)與參考群之間的親緣關(guān)系，這將更有利于基因組選擇技術(shù)的應(yīng)用。另外，如果經(jīng)費(fèi)允許，有必要進(jìn)一步擴(kuò)大參考群規(guī)模，使參考群包含更多家系，增加參考群的多樣性，對(duì)于擴(kuò)展所建立的基因組選擇技術(shù)的適用面、提高基因組選擇的預(yù)測(cè)準(zhǔn)確性和選育效率都具有重要意義[34]。

此外，基因組選擇準(zhǔn)確性還受基因組的標(biāo)記密度的影響[36]。本文通過(guò)隨機(jī)抽樣研究了不同標(biāo)記密度對(duì)育種值預(yù)測(cè)準(zhǔn)確性的影響，當(dāng)使用5 k個(gè)標(biāo)記時(shí)，即可使預(yù)測(cè)準(zhǔn)確性接近利用全基因組SNPs的水平，其后隨著標(biāo)記個(gè)數(shù)增加，預(yù)測(cè)準(zhǔn)確性變化趨于平緩；當(dāng)標(biāo)記個(gè)數(shù)為50 k時(shí)，預(yù)測(cè)準(zhǔn)確性與利用全基因組所有SNPs的預(yù)測(cè)準(zhǔn)確性幾乎一致，這與Tsai H Y等[19]對(duì)大西洋鮭抗海虱性狀基因組預(yù)測(cè)的研究結(jié)果非常相似，提示如果設(shè)計(jì)赤點(diǎn)石斑魚(yú)育種芯片，50 k的標(biāo)記密度可滿足育種要求。盡管如此，利用全基因組SNPs標(biāo)記有利于挖掘具有較大效應(yīng)的分子標(biāo)記及因果突變，將這些大效應(yīng)分子標(biāo)記和因果突變嵌入預(yù)測(cè)模型，建立基于主效-微效多基因效應(yīng)相結(jié)合的預(yù)測(cè)模型，能進(jìn)一步提高育種值預(yù)測(cè)準(zhǔn)確性。利用同一批數(shù)據(jù)，筆者在基因組上已經(jīng)發(fā)現(xiàn)了兩個(gè)效應(yīng)值極強(qiáng)的GWAS信號(hào)(結(jié)果未列出)，計(jì)劃在信號(hào)內(nèi)部挖掘可靠的分子標(biāo)記或因果突變，并將其作為協(xié)變量加入GBLUP模型中，通過(guò)主效標(biāo)記/因果突變進(jìn)一步吸收殘余誤差的方式進(jìn)一步提高預(yù)測(cè)準(zhǔn)確性。對(duì)測(cè)序數(shù)據(jù)進(jìn)行覆蓋度過(guò)濾能夠提高對(duì)挖掘的標(biāo)記位點(diǎn)基因分型的準(zhǔn)確性，從而在一定程度上提高育種值估算的準(zhǔn)確性(圖3)。但是如表1和圖3所示，在每個(gè)個(gè)體測(cè)序量不變的情況下，隨著覆蓋度過(guò)濾標(biāo)準(zhǔn)提高，可保留用于分析的標(biāo)記數(shù)量急劇減少，當(dāng)DP≥4時(shí)，盡管標(biāo)記分型的準(zhǔn)確度會(huì)明顯提升，但是由于保留的標(biāo)記數(shù)量減少，育種值預(yù)測(cè)準(zhǔn)確性已低于不進(jìn)行覆蓋度過(guò)濾。據(jù)Zhang W等[37]對(duì)大黃魚(yú)研究的結(jié)果，采用全基因組重測(cè)序，由于可挖掘到的標(biāo)記數(shù)多，即使測(cè)序覆蓋度低至0.5×，而且不進(jìn)行覆蓋度過(guò)濾，也能獲得與8×覆蓋度測(cè)序基本一致的基因組選擇效果，這將大大降低候選親本標(biāo)記分型費(fèi)用。因此，可以預(yù)期，隨著高通量DNA測(cè)序價(jià)格進(jìn)一步降低，基因組重測(cè)序?qū)⒃絹?lái)越多被用于基因組選擇，也將成為石斑魚(yú)全基因組選擇育種的主要工具。