摘 要： 為了解析半滑舌鰨（Cynoglossus semilaevis）性逆轉(zhuǎn)性狀的分子遺傳作用機制，定位篩選可用于性控育種的分子標(biāo)記或侯選基因，本研究提出了一種期望最大化算法（ Expectation-Maximization algorithm， EM），并基于該算法開展了半滑舌鰨性逆轉(zhuǎn)性狀的全基因組關(guān)聯(lián)分析。EM算法直接使用閾模型中隱含連續(xù)正態(tài)分布表型的期望作為因變量，用迭代最小二乘代替logit 回歸法的迭代重加權(quán)最小二乘，它具有比logit 回歸法更直觀、更易于編程的優(yōu)點。本研究采用顯著主成分控制群體分層后，使用EM 算法與logit 回歸對對半滑舌鰨數(shù)據(jù)進行GWAS（Genome-wideAssociation Study， GWAS）分析。結(jié)果顯示，EM算法結(jié)果無明顯的假陽性或假陰性，比logit 回歸法的檢測效力更高?；贓M算法的全基因組關(guān)聯(lián)分析共定位到13 個與性逆轉(zhuǎn)性狀顯著關(guān)聯(lián)的QTN （quantitative trait nucleotide，QTN），其中3 個QTN位于W染色體上，10 個QTN位于Z 染色體上。經(jīng)過基因注釋發(fā)現(xiàn)，上述定位獲得的QTN位于LOC103396896、MALT1、ADGRD2、FBXl17、DMXl1、SMARCA2、DMRT1、LOC103397760、NEUR13 和PDLIM5a 基因區(qū)段內(nèi)。當(dāng)進行檢索時發(fā)現(xiàn)，這些基因參與了其他物種中涉及性別決定或性腺發(fā)育等相關(guān)過程。本研究提供了一種基于EM算法的具有高檢測效力的全基因組關(guān)聯(lián)分析方法，同時也為半滑舌鰨的性逆轉(zhuǎn)遺傳機制解析和性控育種提供有效的理論指導(dǎo)。

關(guān)鍵詞： 全基因組關(guān)聯(lián)分析；半滑舌鰨；性逆轉(zhuǎn)；主成分；期望最大化算法；廣義線性模型

中圖法分類號： S917.4 文獻標(biāo)識碼： A 文章編號： 1000-2324（2024）04-0531-09

半滑舌鰨作為重要的海水養(yǎng)殖魚類，其生長發(fā)育表現(xiàn)出巨大的性別差異，雄性比雌性成熟早、生長慢，雌性成熟個體的體重可達雄性的2-4 倍［1， 2］。由于這一特性，性別差異是半滑舌鰨的一種極其重要的經(jīng)濟性狀［3， 4］。半滑舌鰨的遺傳性別決定是由ZW染色體控制的，即“遺傳雄魚”為ZZ 型染色體，“遺傳雌魚”為ZW型染色體［5］。雖然半滑舌鰨的偽雄魚現(xiàn)象的機制尚未有實質(zhì)進展，但當(dāng)前的研究主要集中在于確定遺傳因素和后天環(huán)境因素對此現(xiàn)象的影響程度［5-6］。由于這種性別轉(zhuǎn)變是不可逆的，造成了生理雄魚比例過高，嚴重影響了其養(yǎng)殖經(jīng)濟效益。所以了解半滑舌鰨的性別決定機制進從而挖掘性別關(guān)聯(lián)標(biāo)記將為單雌性品種選育或高雌性比例養(yǎng)殖提供理論基礎(chǔ)，有助于降低養(yǎng)殖成本，獲得更大的經(jīng)濟效益［6-8］。

性逆轉(zhuǎn)現(xiàn)象在昆蟲、爬行動物、兩棲動物和魚類中均有發(fā)現(xiàn)［9-13］，性逆轉(zhuǎn)這類僅由0/1 數(shù)據(jù)所表示的間斷性狀被稱為二元性狀，這類性狀一般不符合簡單的孟德爾遺傳定律，需要借助關(guān)聯(lián)分析方法進行基因定位［14］。為定位半滑舌鰨的性別決定位點，早期相關(guān)研究多采用擴增片段長度多態(tài)性（ AFLPs） 分子標(biāo)記、微衛(wèi)星標(biāo)記等大分子標(biāo)記進行基因分型［15-17］。隨著簡化基因組測序技術(shù)限制性位點關(guān)聯(lián)DNA技術(shù)（ Restrictionsite associated DNA， RAD）的應(yīng)用和高密度連鎖圖譜的繪制［18， 19］，高通量高密度的單核苷酸多態(tài)分子標(biāo)記（Single Nucleotide Polymorphisms，SNP）被獲得，在此基礎(chǔ)上進行全基因組關(guān)聯(lián)分析（Genome-wide association studies， GWAS）要比大分子標(biāo)記更能準(zhǔn)確的定位半滑舌鰨的性別相關(guān)的性狀［20， 21］。此前的GWAS 研究發(fā)現(xiàn)半滑舌鰨在Z 染色體上可能存在參與性逆轉(zhuǎn)發(fā)育過程的基因，例如：FBXl17 基因、si：deky-193c22.1基因、DAPK1（ death-associated protein kinase 1）基因、ADGRD2 （adhesion G protein-coupledreceptor D2）基因、DMXL1 （Dmx like 1）基因和LOC103396896 基因［22-26］。這些針對舌鰨性逆轉(zhuǎn)性狀的GWAS研究多采用線性模型，然而由于群體分層會導(dǎo)致違背線性模型（linear model， LM）的正態(tài)分布假設(shè)中的殘差方差，這可能導(dǎo)致關(guān)聯(lián)分析結(jié)果呈現(xiàn)假陽性，所以LM 通常不適合分析二元性狀［27］。對于二元性狀，有研究建議使用基于logit 回歸的廣義線性模型（Generalize LinearModel，GLM）［27］。廣義線性模型作為一種非線性回歸模型，所分析的數(shù)據(jù)可以具有非線性或者非恒定的方差結(jié)構(gòu)［14］。

雖然當(dāng)分子遺傳標(biāo)記的效應(yīng)很小且沒有群體分層時，LM計算的p 值和基于GLM的logit 回歸法結(jié)果可以近似等價［28］，但實施線性模型的前提是表型數(shù)據(jù)要滿足正態(tài)分布假設(shè)，否則會有假陽性的風(fēng)險［27］，所以應(yīng)采用基于logit 回歸的廣義線性模型去分析二元性狀。然而，每次求解廣義線性模型需要進行迭代重加權(quán)最小二乘運算，這將比線性模型消耗更多時間和計算機內(nèi)存，特別是當(dāng)協(xié)變量比較多或所分析的群體較小時，logit回歸還會產(chǎn)生嚴重的估計偏差［ 29， 30］。

本研究提出了一種基于期望最大化的廣義線性模型算法用于分析半滑舌鰨性逆轉(zhuǎn)性狀。EM法不需要去求解二元性狀閾模型［14］中隱含的表型，而是將閾模型函數(shù)產(chǎn)生的隱含連續(xù)正態(tài)分布的期望作為表型，此時用于求解計算的迭代方程與多元線性回歸中的正規(guī)方程組是相同的，直接采用簡單回歸求解遺傳標(biāo)記的效應(yīng)代替logit 回歸的加權(quán)最小二乘進行運算。該期望最大化算法具有比單純Logit 算法更直觀、更易于編程的優(yōu)點［31］，同時也降低了廣義線性模型的求解難度，避免了由權(quán)重引起的異常解。本研究在矯正群體分層后，再使用EM算法對半滑舌鰨性逆轉(zhuǎn)性狀進行GWAS 分析，并與logit 回歸分析結(jié)果作比較，最后檢驗新方法的統(tǒng)計性質(zhì)是否正常和對QTN的檢測效力是否有所提高，同時對新檢測到的位點進行分析。這些研究結(jié)果不但為半滑舌鰨性逆轉(zhuǎn)性狀的研究提供理論依據(jù)，同時也為魚類二元性狀分析提供了新的模型方法。