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

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

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

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

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

1 材料與方法

1.1 樣本的采集與鑒定

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

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

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

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

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

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

1.4 遺傳參數(shù)估計

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

y=Xβ+Zμ+e

(1)

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

(3)

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

1.5 交叉驗證

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

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

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

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

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

2 結(jié)果

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

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

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

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

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

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

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

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

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

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

3 討論

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

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

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