仿刺參疣足數(shù)量SNP 遺傳力評估*

倪 萍 任 強 王 靜 丁 君 常亞青王揚帆① 胡景杰 包振民

(1. 中國海洋大學(xué) 海洋生物遺傳學(xué)與育種教育部重點實驗室 青島 266003；2. 大連海洋大學(xué) 農(nóng)業(yè)農(nóng)村部北方海水增養(yǎng)殖重點實驗室 大連 116023)

仿刺參(Apostichopus j aponicus)是中國傳統(tǒng)的海珍品之一，營養(yǎng)價值高，藥用價值廣泛(常亞青等,2004)。海參養(yǎng)殖業(yè)在我國經(jīng)過十幾年的迅猛發(fā)展，已成為海水養(yǎng)殖的主要產(chǎn)業(yè)之一(王秀利等, 2006)。然而，由于刺參養(yǎng)殖的過速發(fā)展和不規(guī)范運作，基礎(chǔ)研究滯后于刺參養(yǎng)殖產(chǎn)業(yè)的發(fā)展，加之基礎(chǔ)設(shè)施、養(yǎng)殖工藝和方式的滯后，出現(xiàn)刺參病害頻發(fā)、種質(zhì)資源匱乏、成品參品質(zhì)下降、單位面積產(chǎn)量和產(chǎn)值下降等問題(李成林等, 2010; 張春云等, 2004; 常亞青等,2006; 趙帥等, 2016)；野生資源過度捕撈和養(yǎng)殖刺參長期近親繁殖、累代養(yǎng)殖，也使刺參品種質(zhì)量逐步下降(封巖, 2019)，刺參的養(yǎng)殖業(yè)面臨嚴重威脅。隨著人們對仿刺參需求量的增加，迫切需要對其進行遺傳改良，選育具有優(yōu)良性狀的新品種。遺傳參數(shù)的估計是水產(chǎn)動物選擇育種的一項基礎(chǔ)工作，其中，遺傳力是衡量育種進展和育種方法的關(guān)鍵性指標(biāo)(Falconer et al, 2000)。

遺傳力反映親屬間的相似程度和對人工與自然選擇的反應(yīng)速度，往往是育種工作的第一步。傳統(tǒng)方法利用親緣關(guān)系較近的個體表型間的相關(guān)系數(shù)及親緣系數(shù)計算遺傳力(Lstib?rek et al, 2018; 文超良等,2019)，在仿刺參遺傳力估計中被廣泛應(yīng)用，例如，對仿刺參前期浮游階段耳狀幼體體長(欒生等, 2006)、耳狀幼體和稚參體長(李云峰等, 2009)、幼參4 個生長性狀(孟思遠等, 2010)的遺傳力估計等。傳統(tǒng)方法利用系譜信息構(gòu)建親緣相關(guān)矩陣(Numerator relationship matrix, A 陣)，并利用線性混合模型估計遺傳力，對系譜信息的準(zhǔn)確性和完整性要求較高(李晶等,2020)；由于孟德爾抽樣誤差(Hill, 2014)，基于系譜推斷的個體間親緣關(guān)系準(zhǔn)確性有限，并且存在系譜記錄錯誤、部分群體系譜信息不完整或缺失的情況，使得傳統(tǒng)法估計遺傳力的準(zhǔn)確性存在一定的局限。

隨著測序技術(shù)的發(fā)展，越來越多的學(xué)者利用單核苷酸多態(tài)性(SNP)分子標(biāo)記推斷個體間的親緣關(guān)系(Thompson, 1975; Lynch, 1988; Queller, 1989)，通過構(gòu)建基因組親緣關(guān)系矩陣(Genomic relationship matrix,GRM)(Li et al, 2014; 文超良等, 2019)，進行基于SNP分子標(biāo)記的遺傳力估計(SNP 遺傳力)(Yang et al ,2017)，并在動植物以及人類群體中得到廣泛應(yīng)用(Ritland et al, 1996; Mousseau et al, 1998; Thomas et al, 2002; Visscher et al, 2006)。例如，Yang 等(2010)利用REML 分析方法估計了SNP 數(shù)據(jù)集解釋的人類身高表型的變異程度；Guo 等(2018)利用GCTA 的REML 估計了櫛孔扇貝(Chlamys fa rreri)的生長性狀(扇貝長、寬、高和濕重)的SNP 遺傳力；Fishback 等(2002)使用 REML 方法估計了虹鱒(Oncorhynchus mykiss)總長和重量的遺傳力。Benjamin 等(2012)將這種基于GRM 的REML 分析方法稱為GREML 法。GREML 法可以利用SNP 數(shù)據(jù)估計無親緣關(guān)系的個體之間的遺傳關(guān)系，然后，由全基因組SNP 解釋的表型變異比例推斷遺傳力。與傳統(tǒng)方法相比，REML 分析法適用于非均衡資料的方差組分估計，可大大提高方差組分的估計準(zhǔn)確性。本研究利用已完成測序的仿刺參基因組(Li et al, 2018)，對不同養(yǎng)殖地理位置仿刺參群體進行全基因組重測序，并針對仿刺參重要經(jīng)濟性狀疣足數(shù)量，進行全基因組水平的SNP 遺傳力估計，為豐富仿刺參生長相關(guān)的分子育種理論、加快仿刺參優(yōu)良品種的選育提供一定的理論依據(jù)。

1 材料與方法

1.1 仿刺參樣本與表型

本研究215 個仿刺參樣本分別來自于8 個不同的養(yǎng)殖地理位置：平島、西小磨、山東棲霞口、蚆蛸島、中俄雜交仿刺參、旅順、黃龍尾和俄羅斯海參崴。首先，對所有仿刺參樣本的疣足數(shù)量進行統(tǒng)計(表1)。

表1 仿刺參樣本表型統(tǒng)計Tab.1 Statistics of papillae number for A. japonicus

1.2 全基因組重測序與參考SNP 選擇

按照動物組織基因組DNA 提取試劑盒(天根，北京)說明書提取仿刺參基因組DNA。使用NanoDrop-2000 超微量紫外分光光度計進行定量分析。

檢驗合格的DNA 片段經(jīng)末端修復(fù)、加polyA 尾、加測序接頭、純化和PCR 擴增完成整個文庫制備。庫檢合格的文庫通過Illumina Hiseq Xten 的PE150bp模式進行測序，每個樣品測序深度為10×。對測序獲得的測序數(shù)據(jù)進行質(zhì)量過濾，使用 Cutadapt 軟件(Martin, 2011)去除接頭序列，使用SolexaQA 軟件(Cox et al, 2010)去除質(zhì)量值低于20 的堿基，將得到的高質(zhì)量測序數(shù)據(jù)利用BWA 軟件(Li et al, 2009)比對到仿刺參參考基因組。使用Samtools(Li et al, 2009)進行去重復(fù)，GATK(Ye et al, 2009)進行局部重比對，堿基質(zhì)量值校正等處理以及 SNP 小片段插入缺失(Small INDEL)的檢測，按照條件“QD＜2.0||MQ＜40.0||FS＞60.0||MQRankSum＜-12.5||ReadPosRankSum＜-8.0”對SNP 進行過濾，得到最終的SNP 位點集。

篩選SNP 是為了在重測序獲得的高密度SNP 中減少或刪除處于高度連鎖不平衡的SNP，特別是在有限群體中有相當(dāng)數(shù)量的SNP 之間存在高度的連鎖不平衡關(guān)系，會降低對動物個體遺傳參數(shù)的估計準(zhǔn)確度。本研究采用 SNP 基因頻率的平均歐式距離(Average Euclidean distance, AED) (Wu et al, 2016)來篩選SNPs，挑選出次要等位基因頻率(Minor allele frequency, MAF)大于0.05 和0.1 的50K SNP，在篩選出來的 50K SNP 中隨機均勻挑選不同密度的SNP(SNP 數(shù)量分別為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)對仿刺參疣足數(shù)量的SNP 遺傳力進行評估。EM 算法是由Dempster 等(1977)提出的一種迭代算法，可針對有缺失數(shù)據(jù)的數(shù)據(jù)資料，獲得參數(shù)的最大似然估計值。缺失數(shù)據(jù)在實際的統(tǒng)計分析中非常普遍，EM 算法的基本思想是用缺失數(shù)據(jù)在給定參數(shù)的某個估計值的條件下的期望代替缺失數(shù)據(jù)，從而將不完全數(shù)據(jù)變成完全數(shù)據(jù)，使得似然函數(shù)的最大化變得相對簡單。REML 的 AI 算法(平均信息算法)由Gilmour 等(1995)、Johnson 等(1995)和Jensen 等(1997)提出，它是將求似然函數(shù)最大值的 2 種常見的Newton-Raphson 算法和Fisher's scoring 方法結(jié)合起來的一種算法(張勤, 2007)。

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

式中，y 表示表型，b 表示固定效應(yīng)協(xié)變量的系數(shù)，u 表示隨機效應(yīng)自變量的系數(shù)，這里的隨機效應(yīng)指的是所有SNP 位點對表型的效應(yīng)，e 表示隨機誤差。

其中，u 和e 服從如下正態(tài)分布：I為一個n × n 的單位矩陣，n 表示樣本量。

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

W 為第ij 元素的標(biāo)準(zhǔn)化基因型矩陣，xij為第j 個個體的第i 個SNP 的參考等位基因的拷貝數(shù)，pi為參考等位基因頻率。

GCTA 對于樣本遺傳相似度的定義公式：

定義所有SNP 位點的方差：

式中，N 為SNP 位點數(shù)，G 為個體間的親緣關(guān)系矩陣(GRM)。

定義SNP 遺傳力的計算公式：

對來自不同養(yǎng)殖地理位置的215 個仿刺參樣本進行高通量重測序，獲得一組高質(zhì)量SNP。構(gòu)建GRM矩陣后，使用 GCTA 的 2 種遺傳參數(shù)估計方法(AI-REML 和EM-REML)，添加協(xié)變量的校正，針對仿刺參疣足數(shù)量，分別計算不同密度SNP 下的SNP遺傳力估計值及染色體水平的SNP 遺傳力估計值。

2 結(jié)果與分析

2.1 仿刺參疣足數(shù)量的表型數(shù)據(jù)分布

從215 個仿刺參樣本疣足數(shù)量的頻率直方圖和概率密度曲線(圖1a)可以看出，本研究不同地域的仿刺參疣足數(shù)量大致符合正態(tài)分布，表明仿刺參疣足數(shù)量性狀不受地域環(huán)境的影響?？紤]到圖1a 中出現(xiàn)大量疣足數(shù)量大于60 的個體，且主要集中在中俄雜交仿刺參和俄羅斯海參崴群體中，因此，利用R MASS包中nlminb 函數(shù)，對215 個樣本數(shù)據(jù)進行混合正態(tài)分布參數(shù)的最大似然估計(圖1b)，結(jié)果顯示，樣本疣足數(shù)量符合2 個正態(tài)分布(均值39.8、方差6.79 和均值54.35、方差10.9)的最優(yōu)擬合，說明樣本中存在多疣足仿刺參群體情況。進一步把8 個不同地域的仿刺參群體分為正常群體(平島、西小磨、山東棲霞口、蚆蛸島、旅順和黃龍尾)與多疣足群體(中俄雜交仿刺參和俄羅斯海參崴)，并利用R t.test 函數(shù)對上述2 個群體疣足數(shù)量進行t 檢驗，結(jié)果差異顯著(P＜0.01)。根據(jù)上述表型數(shù)據(jù)分布結(jié)果，本研究將不同地域的仿刺參群體分為正常群體與多疣足群體，并在SNP 遺傳力估計模型設(shè)計群體固定效應(yīng)，以校群體分層對SNP 遺傳力估計的影響。

圖1 215 份仿刺參樣本疣足數(shù)量的頻率直方圖和概率密度曲線(a)和雙重高斯分布最優(yōu)擬合(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 仿刺參疣足數(shù)量SNP 遺傳力估計

本研究全基因組重測序質(zhì)量較高，GC 分布正常，對基因組的覆蓋度均超過90%，高質(zhì)量SNPs 數(shù)量約為500 萬個。使用GCTA 的2 種REML 分析方法(EM-REML 和AI-REML)對仿刺參疣足數(shù)量的遺傳力進行評估，比較相同密度、相同次要等位基因頻率篩選條件下，不同分析方法得到的仿刺參疣足數(shù)量SNP遺傳力估計值。發(fā)現(xiàn)通過2 種方法得到的遺傳力估計值差異極小(表2，表3)，證明本研究采用的分析方法有較好穩(wěn)定性、可靠性。

除此之外，本研究對不同密度SNP 條件下的仿刺參疣足數(shù)量SNP 遺傳力估計值進行比較，結(jié)果顯示，當(dāng)MAF＞0.05 時(表2)，不同密度SNP 條件下得到的仿刺參疣足數(shù)量SNP 遺傳力估計均值為(0.566±0.022)～(0.612±0.003)，SNP 數(shù)量在5K、10K 時，得到的SNP 遺傳力估計均值分別為0.567±0.022、0.588±0.011(AI-REML)和0.566±0.022、0.587±0.011(EMREML)；SNP 數(shù)量在15K～50K 時，得到的遺傳力估計均值相較于低密度SNP 遺傳力估計均值更大，且趨于穩(wěn)定；SNP 數(shù)量為50K 時，得到的SNP 遺傳力估計均值最高，分別為0.612±0.003(AI-REML)和0.611±0.003(EM-REML)。當(dāng)MAF＞0.1 時(表3)，不同密度SNP 條件下得到的仿刺參疣足數(shù)量SNP 遺傳力估計均值為(0.586±0.015)～(0.615±0.016)，SNP 數(shù)量在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 數(shù)量在20K～50K 時，得到的遺傳力估計均值相較于低密度SNP 遺傳力估計均值更大，且趨于穩(wěn)定；SNP 數(shù)量為30K 時，得到的SNP 遺傳力估計均值最高，分別為0.615±0.016(AI-REML)和0.614±0.016 (EM-REML)，并且當(dāng)SNP 數(shù)量在35K、40K、50K 時得到的遺傳力估計均值均達到0.614。

表2 仿刺參疣足數(shù)量SNP 遺傳力(MAF>0.05)Tab.2 SNP heritability estimates for papillae number in A. japonicus (MAF＞0.05)

表3 仿刺參疣足數(shù)量SNP 遺傳力(MAF>0.1)Tab.3 SNP heritability estimates for papillae number in A. japonicus (MAF＞0.1)

2.3 仿刺參疣足數(shù)量染色體水平SNP 遺傳力估計

根據(jù)已發(fā)表的仿刺參基因組和遺傳連鎖圖譜，共發(fā)現(xiàn)22 個連鎖群(Tian et al, 2015; Li et al, 2018)，本研究將Scaffold 拼接到染色體水平，并對仿刺參疣足數(shù)量進行染色體水平的SNP 遺傳力估計，對染色體SNP 遺傳力估計值和染色體長度進行回歸分析(圖2)，結(jié)果顯示，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 條件下，單個染色體的貢獻和其長度顯著相關(guān)，較長的染色體具有明顯的線性趨勢(P＜0.05)，且遺傳力估計值最小的7 號染色體和長度最小的14 號染色體均偏離線性回歸較遠。

3 討論

遺傳力是反映性狀遺傳能力大小的重要遺傳參數(shù)，準(zhǔn)確合理的遺傳參數(shù)估計有助于更好地理解遺傳因素對特定群體某一性狀的表型影響程度。估計動物個體的SNP 遺傳力，首先需要篩選出一組適用的SNP。篩選這些SNP，一般需要滿足兩個條件：第一，篩選的 SNP 是該物種基因芯片中共同的SNP，這樣可以很方便地將這些基因SNP 芯片用于分子育種，而不需要重新設(shè)計新的芯片或增加芯片中的SNP 位點；第二，選出的SNP 有較高的信息含量，其統(tǒng)計準(zhǔn)確度要高。Hulsegge 等(2013)采用LD 的r2＞0.30 作為刪除SNP 的尺度，結(jié)果表明，在保持相同準(zhǔn)確性的前提下，使用這個尺度來篩選SNP，可以明顯降低所需SNP 標(biāo)記的數(shù)量。同時，需要篩選信息量高的SNP，篩選高信息量的SNP 可以依據(jù)不同的統(tǒng)計指標(biāo)。例如，Hulsegge 等(2013)分別使用 Delta(群體間等位基因頻率的絕對差異)、Wright(1978)的FST 以及Weir 等(1984)的FST 衡量標(biāo)記信息量的效果。信息熵(Entropy)(Mitt et al, 2017)也是衡量SNP 信息量的指標(biāo)。本研究采用SNP 基因頻率的平均歐式距離(Average Euclidean distance,AED) (Wu et al, 2016)來篩選SNPs。后續(xù)工作將根據(jù)上述指標(biāo)，挑選信息含量高并可進行準(zhǔn)確的遺傳參數(shù)估計的低密度SNP 標(biāo)記。

圖2 仿刺參疣足數(shù)量染色體水平SNP 遺傳力估計Fig.2 Chromosome-wise heritability estimates for papillae number in A. japonicus

群體的分層效應(yīng)會對遺傳力估計值的準(zhǔn)確性產(chǎn)生影響。本研究根據(jù)215 個仿刺參樣本疣足數(shù)量的頻率直方圖和概率密度曲線(圖1)初步將研究樣本分為兩大類，一類為來自中國的品種(正常群體)，另一類為俄羅斯以及中國與俄羅斯雜交的品種(多疣足群體)。計算時使用群體分層這一因素作為協(xié)變量進行校正，并得到了穩(wěn)定結(jié)果，證明本研究使用的分組依據(jù)可靠。本研究使用了基于GCTA 的期望最大約束似然法(EM-REML)和平均信息約束似然法(AI-REML)對仿刺參疣足數(shù)量的SNP 遺傳力進行估計，由估計值結(jié)果(表2、表3)可知，此方法得到的SNP 遺傳力估計值均具有較高的穩(wěn)定性和可信度。此前，孟思遠等(2010)通過建立家系，使用傳統(tǒng)方法估算仿刺參幼參階段肉刺數(shù)目的遺傳力，得到仿刺參肉刺(疣足)數(shù)量遺傳力為0.191～0.404。本研究使用GCTA 的REML分析法得到了與傳統(tǒng)方法相比更高、更穩(wěn)定的遺傳力估計值，說明疣足數(shù)量主要受加性效應(yīng)的控制，對仿刺參的選擇育種具有較大潛力。通過比較不同SNP密度下仿刺參疣足數(shù)量SNP 遺傳力估計值，發(fā)現(xiàn)當(dāng)SNP 數(shù)量在50K 時，得到的遺傳力估計值均已達到穩(wěn)定，說明50K SNP 的密度足夠捕獲數(shù)量性狀基因座(QTL)大效應(yīng)和小效應(yīng)，證實了低密度SNP 用于準(zhǔn)確估計遺傳參數(shù)的可行性。同時，在全基因組選擇方面，Wu 等(2016)通過設(shè)計牛的低密度SNP 芯片，對育種值估計的準(zhǔn)確性進行評估，獲得了較準(zhǔn)確的結(jié)果，也證實了低密度SNP 芯片應(yīng)用于分子育種的可靠性。由于全基因組分析對樣本量要求較高，通常在500 以上(Yang et al, 2010)，本研究的樣本量為215，統(tǒng)計量較少，容易導(dǎo)致統(tǒng)計結(jié)果出現(xiàn)一定程度的偏差，但所得結(jié)果與已報道數(shù)據(jù)(孟思遠等, 2010; 和飛等, 2017)一致。后續(xù)工作將擴大樣本統(tǒng)計量進行仿刺參數(shù)量性狀的遺傳力評估。

本研究在染色體水平對仿刺參疣足數(shù)量SNP 遺傳力進行估計(圖2)，提供性狀遺傳結(jié)構(gòu)的信息。結(jié)果顯示，仿刺參疣足數(shù)量是一個復(fù)雜的數(shù)量性狀，不同染色體上的變異導(dǎo)致了仿刺參疣足數(shù)量的變異，與該性狀相關(guān)的效應(yīng)位點散布在各染色體上，說明了通過重測序的技術(shù)手段，在全基因組上尋找SNP 的重要性。單個染色體的貢獻和其長度顯著相關(guān)，表示在每條染色體上可能有多個小到中等效應(yīng)的變異，而不是少數(shù)具有主要效應(yīng)的變異。Goddard 等(2009)關(guān)于家畜復(fù)雜性狀的分子標(biāo)記輔助選擇結(jié)果也支持這一假設(shè)，該報道稱，性狀變異有很大的遺傳因素，但由于統(tǒng)計力量不足，大多數(shù)具有小到中度影響的變異尚未在遺傳關(guān)聯(lián)分析中被確認。而GCTA 可以克服分子標(biāo)記輔助選擇的缺陷，通過全基因組標(biāo)記面板，所有QTL 都能和至少1 個marker 處于連鎖不平衡狀態(tài)。分析發(fā)現(xiàn)，7 號染色體雖不是最小染色體，但其遺傳力估計值最低，與此相反，2 號染色體雖不是最大染色體，卻具有最大遺傳力估計值，說明較大的基因效應(yīng)位點可能不在7 號染色體上，而位于2 號染色體，意味著QTL 出現(xiàn)在2 號染色體的概率較大，出現(xiàn)在7 號染色體上的概率較小。此結(jié)果的發(fā)現(xiàn)為將來開展仿刺參疣足數(shù)量的全基因組關(guān)聯(lián)分析(GWAS)提供了支持。

4 總結(jié)

本研究對不同養(yǎng)殖地理位置仿刺參群體進行全基因組重測序，并針對仿刺參重要經(jīng)濟性狀疣足數(shù)目，使用GCTA 軟件的GREML 進行全基因組水平的SNP 遺傳力估計，染色體水平SNP 遺傳力估計，同時比較不同密度SNP 標(biāo)記條件下的仿刺參疣足數(shù)量SNP 遺傳力估計值；結(jié)果表明，仿刺參疣足數(shù)量SNP遺傳力為中等遺傳力(MAF＞0.05 時，0.566～0.612；MAF＞0.1 時，0.586～0.615)，當(dāng)SNP 數(shù)量在50K 時，得到的遺傳力估計值均已達到穩(wěn)定，且染色體SNP遺傳力與其長度顯著相關(guān)；說明仿刺參疣足數(shù)量是復(fù)雜數(shù)量性狀，由全基因組水平多基因共同作用，且50K SNP 密度足夠捕獲QTL 大效應(yīng)和小效應(yīng)。本研究通過分子學(xué)研究方法得到可靠的遺傳力估計值，有效解決了傳統(tǒng)方法系譜記錄繁瑣，遺傳力估計準(zhǔn)確性差等問題，為分子育種技術(shù)提供了一定理論依據(jù)；同時，證明低密度SNP 標(biāo)記用于仿刺參分子育種的可能性，可為低密度SNP 芯片的設(shè)計開發(fā)提供支持。