赤點石斑魚抗神經壞死病的基因組選擇評估

王卓標 ，方 銘，鄭樂云，葛 輝，陳欣欣，羅輝玉，王志勇*

(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.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 結果

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數

3 討論

基因組選擇的首要工作建立參考群，利用參考群估算各個分子標記位點對性狀表型的效應。本文建立了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測序價格進一步降低，基因組重測序將越來越多被用于基因組選擇，也將成為石斑魚全基因組選擇育種的主要工具。