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

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

中圖法分類(lèi)號(hào)： S917.4 文獻(xiàn)標(biāo)識(shí)碼： A 文章編號(hào)： 1000-2324（2024）04-0531-09

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

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

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

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