馬 娟 劉京寶 朱衛(wèi)紅 黃 璐 宇 婷 喬江方

（河南省農(nóng)業(yè)科學(xué)院糧食作物研究所，河南 鄭州 450002）

玉米（Zea maysL.）是雜種優(yōu)勢(shì)利用最成功的作物之一，主要體現(xiàn)在單交種的生產(chǎn)和利用上［1］。在玉米雜交育種中，選育高配合力的親本材料是培育高產(chǎn)雜交種的先決條件［2］。一般配合力（general combining ability，GCA）是指一個(gè)自交系與其他自交系組配后代在某個(gè)性狀的平均表現(xiàn)，是評(píng)價(jià)親本自交系利用價(jià)值的一個(gè)重要指標(biāo)［1，3］。因此解析GCA 的遺傳機(jī)理對(duì)玉米雜交種產(chǎn)量的提高具有重要意義。

GCA 和農(nóng)藝性狀均是多基因控制的數(shù)量性狀，可以利用連鎖定位和關(guān)聯(lián)分析方法開(kāi)展研究。Qi等［4］利用75 份玉米導(dǎo)入系材料和4 個(gè)測(cè)驗(yàn)種定位到75 個(gè)控制產(chǎn)量及相關(guān)性狀GCA 的數(shù)量性狀位點(diǎn)（quantitative trait loci，QTL），其中穗行數(shù)GCA的顯著QTL有6個(gè)。Lu等［5］以328 個(gè)玉米重組自交系和2 個(gè)測(cè)驗(yàn)種組配的656個(gè)雜交種為材料，利用復(fù)合區(qū)間作圖方法挖掘64 個(gè)控制產(chǎn)量等性狀GCA 的關(guān)鍵位點(diǎn)。Liu 等［6］利用1 個(gè)重組自交系群體和4個(gè)測(cè)驗(yàn)種組配的雜交群體鑒定41個(gè)控制玉米穗行數(shù)、行粒數(shù)等性狀GCA 的關(guān)鍵位點(diǎn)和28個(gè)熱點(diǎn)區(qū)域。李婷等［3］利用EMMAX（efficient mixedmodel association expedited）算法對(duì)246 份玉米雜交組合的籽粒相關(guān)性狀配合力進(jìn)行了關(guān)聯(lián)分析，檢測(cè)到9、7 和5 個(gè)單核苷酸多態(tài)性（singe nucleotide polymorphism，SNPs）分別與粒長(zhǎng)、粒寬和粒厚GCA 顯著關(guān)聯(lián)，單個(gè)位點(diǎn)解釋GCA表型變異率為0.02%～25.52%。

盡管已有研究者挖掘到一些控制GCA 的關(guān)鍵位點(diǎn)，但是這些位點(diǎn)的遺傳利用尚缺乏相關(guān)的分子改良策略?；蚪M選擇（genomic selection，GS）是利用訓(xùn)練群體覆蓋全基因組的分子標(biāo)記和表型數(shù)據(jù)建立數(shù)學(xué)模型，從而對(duì)基因型已知的候選群體進(jìn)行表型預(yù)測(cè)的一種分子育種技術(shù)。GS已廣泛應(yīng)用于預(yù)測(cè)玉米、水稻和小麥等作物親本和雜交種產(chǎn)量等表現(xiàn)，為GS 輔助育種技術(shù)的應(yīng)用提供了重要參考。王欣等［7］基于不完全雙列雜交設(shè)計(jì)（North Carolina II，NCII）的水稻數(shù)據(jù)集，研究了單株產(chǎn)量等性狀GCA的預(yù)測(cè)力，發(fā)現(xiàn)單株產(chǎn)量等8個(gè)性狀的預(yù)測(cè)準(zhǔn)確性介于0.39～0.74。Zhang等［8］利用32個(gè)玉米自交系和9 個(gè)測(cè)驗(yàn)種組配的3 組測(cè)交群體為材料對(duì)產(chǎn)量GCA 開(kāi)展預(yù)測(cè)研究，發(fā)現(xiàn)僅利用自交系的信息對(duì)產(chǎn)量GCA 的預(yù)測(cè)精度較低，為-0.14～0.13，當(dāng)加入測(cè)驗(yàn)種的信息時(shí)，產(chǎn)量GCA 的預(yù)測(cè)準(zhǔn)確性可提高到0.51～0.65。

穗行數(shù)是玉米產(chǎn)量的重要構(gòu)成因子之一，籽粒性狀例如粒長(zhǎng)、粒寬與玉米品質(zhì)及產(chǎn)量密切相關(guān)，解析其GCA 遺傳機(jī)理對(duì)玉米雜交種產(chǎn)量和品質(zhì)的遺傳改良具有重要意義。雖然穗行數(shù)和籽粒性狀以及其分子遺傳研究基礎(chǔ)已有相關(guān)報(bào)道，但目前關(guān)于穗行數(shù)和籽粒性狀GCA 全基因組關(guān)聯(lián)分析（genome-wide association study，GWAS）的研究較少，而且其基因組預(yù)測(cè)分析的研究仍鮮見(jiàn)報(bào)道。鑒于此，本研究選用123份玉米自交系和8個(gè)測(cè)驗(yàn)種按照NCII 遺傳交配設(shè)計(jì)獲得537 份雜交組合為材料，采用7 種多位點(diǎn)GWAS（multi-locus GWAS，MGWAS）方法挖掘穗行數(shù)、粒長(zhǎng)、粒寬GCA顯著關(guān)聯(lián)位點(diǎn)，并研究不同顯著關(guān)聯(lián)位點(diǎn)作為固定效應(yīng)對(duì)提高3個(gè)性狀GCA預(yù)測(cè)精度的影響，旨在為基因功能驗(yàn)證以及關(guān)鍵位點(diǎn)的基因組選擇輔助育種提供重要的基因信息和技術(shù)指導(dǎo)。

1 材料與方法

1.1 試驗(yàn)材料、試驗(yàn)設(shè)計(jì)和配合力分析

以123份玉米自選系和8個(gè)測(cè)驗(yàn)種（農(nóng)系531、昌7-2、20H1419、M189、M119、S110T、L119A 和PH4CV）采用NCII遺傳交配設(shè)計(jì)組配537個(gè)F1雜交種。537個(gè)雜交種采用隨機(jī)完全區(qū)組試驗(yàn)設(shè)計(jì)，于2021 年種植在河南省新鄉(xiāng)和周口試驗(yàn)田，2次重復(fù)。新鄉(xiāng)和周口環(huán)境行長(zhǎng)分別是4.0和3.3 m，株距和行距均分別是0.22和0.60 m。收獲后，每份材料選取2 個(gè)代表性果穗人工測(cè)量穗行數(shù)。利用玉米果穗/籽?？挤N流水線儀器（國(guó)家農(nóng)業(yè)信息化工程技術(shù)研究中心）測(cè)量約20 粒籽粒的粒長(zhǎng)和粒寬，獲得單粒粒長(zhǎng)和粒寬用于分析。利用R 語(yǔ)言lme4包的lmer 函數(shù)計(jì)算2 個(gè)環(huán)境的GCA 效應(yīng)值（公式1），并將2 個(gè)環(huán)境數(shù)據(jù)聯(lián)合計(jì)算綜合環(huán)境的GCA 效應(yīng)值（公式2），均用于后續(xù)分析：

其中，y、μ、B、E、GCAL、GCAT、SCA以及ε分別表示雜交種表型、均值、區(qū)組、環(huán)境、自交系GCA、測(cè)驗(yàn)種GCA、特殊配合力以及誤差效應(yīng)。GCAL×E、GCAT×E和SCA×E分別表示自交系GCA 與環(huán)境互作、測(cè)驗(yàn)種GCA 與環(huán)境互作和特殊配合力與環(huán)境互作效應(yīng)。

根據(jù)公式2，利用R 語(yǔ)言stats 包aov 函數(shù)將新鄉(xiāng)和周口環(huán)境數(shù)據(jù)聯(lián)合進(jìn)行多環(huán)境方差分析。不同環(huán)境3個(gè)性狀GCA 相關(guān)性分析（Pearson）和顯著性檢驗(yàn)利用R語(yǔ)言GGally包計(jì)算和展示。

GCA遺傳力的計(jì)算公式如下：

1.2 基因型檢測(cè)和分析

131份玉米自交系（123份自交系和8個(gè)測(cè)驗(yàn)種）葉片的DNA 采用CTAB 方法提取。利用玉米5.5K 液相育種芯片（河南省農(nóng)業(yè)科學(xué)院糧食作物研究所）根據(jù)液相探針雜交原理進(jìn)行靶向測(cè)序基因型鑒定，測(cè)序平臺(tái)為NovaSeq 6000 （Illumina，美國(guó)）。利用BWA（v0.7.17）軟件將過(guò)濾的reads 與玉米B73 第四版參考基因組（http：//www.gramene.org/）比對(duì)。利用GATK （v4.1.2.0）軟件檢測(cè)原始變異位點(diǎn)33 971個(gè)。經(jīng)過(guò)最小等位基因頻率＜0.05，缺失率＞10%和雜合率＞1%的過(guò)濾，獲得11 734個(gè)高質(zhì)量單核苷酸多態(tài)性（single nucleotide polymor hism，SNP）用于后續(xù)分析。

1.3 多位點(diǎn)全基因組關(guān)聯(lián)分析和候選基因挖掘

利用7種MGWAS方法，即BLINK（bayesian information and linkage-disequilibrium iteratively nested keyway）［9］、FarmCPU（fixed and random model circulating probability unification）［10］、mrMLM（multi-locus random-SNP-effect mixed linear model）［11］、FASTmrMLM［12］、FASTmrEMMA（FAST multi-locus random-SNP-effect EMMA）［13］、pLARmEB（polygene-background-control-based least angle regression plus empirical Bayes）［14］和ISIS EM-BLASSO（iterative sure independence screening EM-Bayesian least absolute shrinkage and selection operator）［15］對(duì)3個(gè)環(huán)境穗行數(shù)和籽粒性狀GCA進(jìn)行關(guān)聯(lián)分析。7種方法均考慮群體結(jié)構(gòu)Q值和親緣關(guān)系K值。群體結(jié)構(gòu)Q值由Structure v2.3.4［16］計(jì)算，其中亞群數(shù)為1～10，length of burn-in period為5 000，蒙特卡羅重復(fù)個(gè)數(shù)為50 000，每個(gè)亞群數(shù)迭代次數(shù)為3。利用Structure Harvester（v0.6.94）確定最佳的亞群數(shù)為6。3次重復(fù)的Q值由CLUMPP（v1.1.26）軟件的FullSearch算法 確 定［17］。親 緣 關(guān) 系K 值 由TASSEL v5.0 軟 件 的Centered_IBS 方法計(jì)算［18］。BLINK 和FarmCPU 利用R語(yǔ)言GAPIT 包計(jì)算［19］，其他5 種方法均利用R 語(yǔ)言mrMLM包計(jì)算［20］。顯著臨界值設(shè)置為P=0.01/11 734=8.52×10-7。利用SnpEff［21］根據(jù)顯著SNPs 的位置挖掘候選基因，參數(shù)按默認(rèn)設(shè)置。

1.4 基因組預(yù)測(cè)分析

采用貝葉斯A（Bayes A）［22］、Bayes C［23］、Bayesian LASSO［24］、基因組最佳線性無(wú)偏預(yù)測(cè)（genomic best linear unbiased prediction，GBLUP）［25］和再生核希爾伯特空間（reproducing kernel Hilbert space，RKHS）［26］對(duì)3個(gè)環(huán)境穗行數(shù)、粒長(zhǎng)、粒寬GCA 進(jìn)行預(yù)測(cè)分析。5種方法中分子標(biāo)記的效應(yīng)均設(shè)置為隨機(jī)效應(yīng)即隨機(jī)效應(yīng)模型。根據(jù)不同MGWAS 方法挖掘的關(guān)聯(lián)位點(diǎn)信息，將顯著SNPs 作為固定效應(yīng)，其余分子標(biāo)記作為隨機(jī)效應(yīng)（固定效應(yīng)模型），用來(lái)研究MGWAS 先驗(yàn)信息對(duì)GCA 預(yù)測(cè)精度的影響。固定效應(yīng)模型僅考慮GBLUP和RKHS。交叉驗(yàn)證均采用10 倍交叉驗(yàn)證方式，重復(fù)100 次。利用100 次驗(yàn)證群體中基因組估計(jì)育種值與表型值的相關(guān)系數(shù)均值作為預(yù)測(cè)準(zhǔn)確性。所有模型均在R 語(yǔ)言BGLR 包中實(shí)現(xiàn)［27］，其中RKHS 采用核平均模型，寬帶參數(shù)h為0.1、0.5 和2.5。所有模型和方案迭代次數(shù)為12 000 次，預(yù)燒為3 000 次，其他參數(shù)按默認(rèn)設(shè)置。

2 結(jié)果與分析

2.1 一般配合力效應(yīng)表現(xiàn)統(tǒng)計(jì)結(jié)果

新鄉(xiāng)、周口和綜合環(huán)境穗行數(shù)和籽粒性狀GCA 效應(yīng)值的描述性統(tǒng)計(jì)見(jiàn)圖1。以綜合環(huán)境為例，穗行數(shù)、粒長(zhǎng)和粒寬GCA效應(yīng)值分別介于-2.85～2.62、-0.086～0.081和-0.062～0.069（圖1-A、B）。8個(gè)測(cè)驗(yàn)種在3個(gè)性狀中表現(xiàn)不一：對(duì)于穗行數(shù)，測(cè)驗(yàn)種中正向GCA 最高的是20H1419，為1.95～2.44，其次是昌7-2（0.77～1.03），M189 和PH4CV 的配合力負(fù)向最高，為-1.39～-1.25（圖1-C）；粒長(zhǎng)正向和負(fù)向GCA表現(xiàn)最好的測(cè)驗(yàn)種分別是農(nóng)系531 和20H1419，而粒寬正向GCA 效應(yīng)值最高的是M119，負(fù)向最高的是20H1419（圖1-D）。

圖1 不同環(huán)境穗行數(shù)、粒長(zhǎng)和粒寬一般配合力效應(yīng)值統(tǒng)計(jì)Fig.1 Statistics of general combining ability effect values of kernel row number，kernel length，kernel width in different environments

穗行數(shù)與粒長(zhǎng)之間的相關(guān)性較小，而粒寬與穗行數(shù)表現(xiàn)為高度負(fù)相關(guān)，相關(guān)系數(shù)介于-0.80～-0.61。3個(gè)性狀中，2 個(gè)種植環(huán)境GCA 的相關(guān)系數(shù)較高且達(dá)到極顯著為水平，相關(guān)系數(shù)為0.59～0.78，且均與綜合環(huán)境高度極顯著正相關(guān)（r=0.87～0.95）（圖2）。方差分析表明，除了區(qū)組外，環(huán)境、GCA、特殊配合力以及配合力與環(huán)境互作效應(yīng)均達(dá)到顯著或極顯著水平（表1）。3 個(gè)性狀GCA 的遺傳力較高（0.88～0.95），其中穗行數(shù)GCA 的遺傳力最高，為0.95，表明3 個(gè)性狀GCA 主要受遺傳因素影響。

表1 3個(gè)性狀多環(huán)境方差分析及其遺傳力Table 1 Multi-environment analysis of variance and heritability for three traits

圖2 不同環(huán)境3個(gè)性狀相關(guān)性及顯著性Fig.2 Correlations among three traits in different enviornments and significant levels

2.2 3個(gè)性狀一般配合力顯著關(guān)聯(lián)位點(diǎn)和候選基因

通過(guò)7種MGWAS方法對(duì)3個(gè)環(huán)境穗行數(shù)、粒長(zhǎng)、粒寬GCA進(jìn)行全基因組位點(diǎn)掃描發(fā)現(xiàn)，穗行數(shù)、粒長(zhǎng)、粒寬分別檢測(cè)到23、8和15個(gè)顯著關(guān)聯(lián)SNPs（P＜8.52×10-7）。穗行數(shù)利用mrMLM檢測(cè)的位點(diǎn)最多，粒長(zhǎng)和粒寬均利用mrMLM 和pLARmEB 檢測(cè)的位點(diǎn)最多（電子附表1、圖3）。共有10個(gè)位點(diǎn)被2～5種MGWAS方法同時(shí)檢測(cè)到，其中穗行數(shù)、粒長(zhǎng)和粒寬分別有5、1和4個(gè)（圖3-A～C）。共有8個(gè)位點(diǎn)被至少2個(gè)環(huán)境重復(fù)檢測(cè)到（圖3-D～F）。對(duì)于穗行數(shù)GCA，5_219239213被3個(gè)環(huán)境均檢測(cè)到，5個(gè)為新鄉(xiāng)和綜合環(huán)境共定位的位點(diǎn)。2_197716855 和2_71037706分別為粒長(zhǎng)和粒寬GCA 環(huán)境穩(wěn)定的顯著位點(diǎn)。本研究還發(fā)現(xiàn)6個(gè)SNPs即1_43440622、2_69742504、2_71037706、2_197716855、5_219239213 和8_134634317為環(huán)境穩(wěn)定和MGWAS 方法穩(wěn)定重疊的位點(diǎn)（電子附表1），這些SNPs 可能是控制穗行數(shù)和籽粒性狀GCA效應(yīng)的重要位點(diǎn)。盡管3個(gè)性狀GCA之間表現(xiàn)出一定的相關(guān)性，但在基因組水平?jīng)]有發(fā)現(xiàn)一因多效位點(diǎn)。

圖3 MGWAS方法穩(wěn)定和環(huán)境穩(wěn)定的位點(diǎn)Fig.3 MGWAS method-stable and environment-stable loci

利用SnpEff 軟件在46 個(gè)顯著位點(diǎn)內(nèi)共挖掘候選基因59個(gè)（電子附表1）。6個(gè)環(huán)境穩(wěn)定和MGWAS方法穩(wěn)定重疊位點(diǎn)挖掘的候選基因見(jiàn)表2，其中胚胎敗育基因EMB12（embryo defective12）和類三角形四肽重復(fù)超家族蛋白基因（Zm00001d010946）等是穗行數(shù)GCA重要的候選基因，Zm00001d006084（類尿苷激酶蛋白2，UKL2）和Zm00001d003945分別是粒長(zhǎng)和粒寬GCA重要的候選基因。

表2 MGWAS方法穩(wěn)定和環(huán)境穩(wěn)定重疊位點(diǎn)挖掘的候選基因Table 2 Candidate genes for MGWAS method-stable and environment-stable overlapped loci

2.3 不同隨機(jī)效應(yīng)模型的預(yù)測(cè)準(zhǔn)確性比較

通過(guò)對(duì)穗行數(shù)和籽粒性狀GCA效應(yīng)值進(jìn)行基因組預(yù)測(cè)分析發(fā)現(xiàn)，5種隨機(jī)效應(yīng)模型對(duì)穗行數(shù)和粒寬GCA取得較高的預(yù)測(cè)準(zhǔn)確性，分別為0.62～0.74 和0.63～0.73，而粒長(zhǎng)GCA 基因組預(yù)測(cè)精度較低，為0.28～0.45（圖4）。3個(gè)性狀均利用綜合環(huán)境GCA效應(yīng)值獲得最高的預(yù)測(cè)準(zhǔn)確性。對(duì)于穗行數(shù)和粒寬GCA效應(yīng)值，多數(shù)情況下，RKHS相比其他4種方法略具優(yōu)勢(shì)。對(duì)于粒長(zhǎng)GCA，RKHS在周口和綜合環(huán)境表現(xiàn)出較大優(yōu)勢(shì)，相比其他4種參數(shù)模型，將預(yù)測(cè)準(zhǔn)確性提高了16.63%～31.34%。

圖4 5種隨機(jī)效應(yīng)模型對(duì)3個(gè)性狀一般配合力的預(yù)測(cè)準(zhǔn)確性Fig.4 Prediction accuracy of five random effect models for general combining ability of three traits

2.4 多位點(diǎn)全基因組關(guān)聯(lián)分析先驗(yàn)信息對(duì)預(yù)測(cè)準(zhǔn)確性的影響

率為27.58%～49.81%。在周口環(huán)境，多數(shù)情況下MGWAS先驗(yàn)信息加入能夠提高3 個(gè)性狀GCA 效應(yīng)值的預(yù)測(cè)力，粒長(zhǎng)的提高率較高，為17.65%～83.05%，而穗行數(shù)和粒寬的提高率為0.66%～12.53%。相比隨機(jī)效應(yīng)模型，綜合環(huán)境中固定效應(yīng)模型均可進(jìn)一步提高穗行數(shù)和粒寬GCA 效應(yīng)值的基因組預(yù)測(cè)力，提高率為2.24%～12.46%。但對(duì)于粒長(zhǎng)，將綜合環(huán)境FASTmrMLM挖掘的1 個(gè)顯著位點(diǎn)作為固定效應(yīng)會(huì)提高預(yù)測(cè)精度（9.26%～10.88%），其他2 種MGWAS 方法先驗(yàn)信息與預(yù)測(cè)方法的整合不會(huì)改變或略降低預(yù)測(cè)準(zhǔn)確性。總體而言，將MGWAS 挖掘的顯著關(guān)聯(lián)位點(diǎn)作為固定效應(yīng)，并加入GS 方法，是一種提高穗行數(shù)和籽粒性狀GCA預(yù)測(cè)精度的有效途徑。

3 討論

3.1 MGWAS解析GCA變異位點(diǎn)的優(yōu)勢(shì)

GWAS 能夠檢測(cè)覆蓋全基因組的變異位點(diǎn)，是解析數(shù)量性狀遺傳結(jié)構(gòu)的一種統(tǒng)計(jì)分析方法。2006 年，Yu 等［28］首次提出整合親緣關(guān)系K 矩陣和群體結(jié)構(gòu)Q矩陣的混合線性模型方法，用來(lái)擬合基因型和表型的關(guān)系。為了解決混合線性模型方法存在的混雜問(wèn)題和統(tǒng)計(jì)功效低等問(wèn)題，多種單位點(diǎn)模型和MGWAS 模型如BLINK、FarmCPU、mrMLM 等相繼提出，并已廣泛應(yīng)用于玉米等作物重要性狀的遺傳研究中。本研究利用7 種MGWAS 方法對(duì)穗行數(shù)和2 個(gè)籽粒性狀GCA 效應(yīng)進(jìn)行了關(guān)聯(lián)分析，共鑒定得到46個(gè)顯著關(guān)聯(lián)SNPs，其中mrMLM方法的檢測(cè)位點(diǎn)最多。An等［29］利用6種MGWAS對(duì)穗行數(shù)進(jìn)行定位研究發(fā)現(xiàn)，mrMLM 在4 個(gè)環(huán)境中的檢測(cè)功效均最高。玉米雄穗分枝數(shù)的MGWAS 研究也得到了相同的結(jié)論［30］。本研究中，僅10個(gè)位點(diǎn)被多種MGWAS方法重復(fù)檢測(cè)到，多數(shù)位點(diǎn)表現(xiàn)為方法特異性。但對(duì)于穗行數(shù)本身，An 等［29］發(fā)現(xiàn)40%～62%的顯著位點(diǎn)至少被2 種MGWAS 方法重復(fù)檢測(cè)到。Zhou 等［31］對(duì)成熟期籽粒含水量進(jìn)行GWAS 分析時(shí)，也發(fā)現(xiàn)方法相同的位點(diǎn)具有較高的占比（約45%）。本研究中，8 個(gè)位點(diǎn)被至少2 個(gè)環(huán)境重復(fù)檢測(cè)到，是控制穗行數(shù)和籽粒性狀GCA 環(huán)境穩(wěn)定的SNPs。其中6 個(gè)SNPs 不僅被多種MGWAS 方法重復(fù)檢測(cè)到，還是環(huán)境穩(wěn)定的位點(diǎn)。因此，利用多種MGWAS 方法有助于挖掘GCA 穩(wěn)定的變異位點(diǎn)。

3.2 不同定位研究結(jié)果比較

本研究發(fā)現(xiàn)，15個(gè)穗行數(shù)GCA顯著位點(diǎn)與已報(bào)道的穗行數(shù)GCA、性狀本身、產(chǎn)量QTL 和Meta-QTL（MQTL）存在重疊（電子附表1）。其中，位點(diǎn)4_203093326位于3個(gè)控制穗行數(shù)GCA QTL［6］以及穗行數(shù)和產(chǎn)量MQTL 區(qū)間內(nèi)（MQTL69）［32］。位于2號(hào)染色體上的4個(gè)顯著位點(diǎn)（2_153730555、2_159709723、2_159709745和2_225337449）被發(fā)現(xiàn)與穗行數(shù)QTL以及產(chǎn)量、穗部等性狀MQTL區(qū)間存在重疊［32-35］。其中2_153730555 是MGWAS 方法穩(wěn)定的一個(gè)控制穗行數(shù)GCA的主效位點(diǎn)。4個(gè)環(huán)境穩(wěn)定和MGWAS方法穩(wěn)定的SNPs均發(fā)現(xiàn)位于前人定位的產(chǎn)量QTL或MQTL區(qū)間內(nèi)（電子附表1），表明這4個(gè)SNPs可能是控制穗行數(shù)GCA效應(yīng)的重要位點(diǎn)。

圖5 GBLUP和RKHS整合MGWAS顯著位點(diǎn)的預(yù)測(cè)準(zhǔn)確性Fig.5 Prediction accuracy of the integration of GBLUP and RKHS with significant loci from MGWAS

19個(gè)籽粒性狀GCA位于前人定位的籽粒性狀GCA、百粒重和產(chǎn)量等性狀QTL或MQTL區(qū)間內(nèi)。粒寬GCA顯著位點(diǎn)6_144964691 不僅位于一個(gè)粒寬GCA QTL（qKW6-4）的置信區(qū)間內(nèi)［5］，還與控制百粒重和單穗粒重QTL區(qū)間存在重疊［34，36］。3_19172706也是粒寬GCA顯著關(guān)聯(lián)位點(diǎn)，被發(fā)現(xiàn)位于Lu 等［5］挖掘的一個(gè)對(duì)粒寬和百粒重GCA具有一因多效的QTL區(qū)間內(nèi)。1_268584002是一個(gè)控制粒長(zhǎng)GCA的主效SNP，與前人鑒定的粒長(zhǎng)、粒寬及其GCA、百粒重以及單穗粒重QTL 位置存在重疊［5，34，36-37］。粒長(zhǎng)關(guān)聯(lián)位點(diǎn)4_186588788 不僅與Lu 等［5］和Liu 等［36］定位的粒寬及粒寬GCA QTL 區(qū)間重疊，還與一個(gè)控制籽粒性狀、產(chǎn)量及穗部性狀MQTL 存在重疊［33］。2_197716855 和2_71037706 分別是控制粒長(zhǎng)和粒寬GCA 的主效位點(diǎn)，不僅被多個(gè)MGWAS 重復(fù)檢測(cè)到，還被2個(gè)環(huán)境同時(shí)檢測(cè)到。位點(diǎn)2_197716855位于Lu 等［5］和Liu 等［36］鑒定的粒寬及產(chǎn)量QTL 的置信區(qū)間內(nèi)，而2_71037706 位于Chen 等［33］利用元分析整合鑒定的一個(gè)控制籽粒性狀MQTL區(qū)間內(nèi)。粒寬GCA位點(diǎn)2_209342026 和6_90473152 至 少 被2 種MGWAS 方 法同時(shí)檢測(cè)到，發(fā)現(xiàn)與Liu 等［36］和Zhang 等［37］定位的2個(gè)影響粒寬的數(shù)量性狀位點(diǎn)存在重疊。

3.3 一般配合力候選基因功能挖掘

MGWAS方法穩(wěn)定和環(huán)境穩(wěn)定重疊位點(diǎn)5_219239213位于一個(gè)控制玉米產(chǎn)量QTL的區(qū)間內(nèi)［5］（電子附表1），該位點(diǎn)對(duì)穗行數(shù)GCA效應(yīng)的平均變異解釋率為12.57%，是一個(gè)主效位點(diǎn)，其預(yù)測(cè)的基因?yàn)镋MB12。該基因編碼質(zhì)體翻譯起始因子3 在玉米胚胎發(fā)生發(fā)育過(guò)程中起著關(guān)鍵作用［38］。8_134634317 也是一個(gè)MGWAS 方法穩(wěn)定和環(huán)境穩(wěn)定重疊的位點(diǎn)，其挖掘的候選基因Zm00001d010946是一類四肽重復(fù)序列（tetratricopeptide repeat domain，TPR）蛋白（表2）。TPR 蛋白在擬南芥根發(fā)育、葉綠體早期發(fā)育和生長(zhǎng)素極性運(yùn)輸中起著重要作用［39-40］。Flo2（FLOURY ENDOSPERM2）是一個(gè)具有TPR基序且保守的基因，在調(diào)控水稻和小麥粒重和淀粉品質(zhì)中發(fā)揮重要作用［41-42］。此外，三角狀五肽重復(fù)蛋白基因PPR（pentatricopeptide repeat-containing protein）和抗壞血酸鹽過(guò)氧化物酶基因APX2基因（ascorbate peroxidase2）也是可能控制穗行數(shù)GCA 的重要基因（附表1），其中APX2對(duì)應(yīng)的SNP與Zhang等［35］定位的一個(gè)穗行數(shù)QTL 位置重疊。水稻育性恢復(fù)基因Rf4、玉米雄穗分枝數(shù)的一個(gè)主效QTL以及玉米種子發(fā)育關(guān)鍵基因Dek39均編碼PPR 基因［43-45］。OsAPX2的功能缺失會(huì)影響水稻（Oryza sativaL.）幼苗的生長(zhǎng)發(fā)育，導(dǎo)致水稻幼苗半矮化、葉片黃綠色和種子不育［46］。但這些基因調(diào)控穗行數(shù)GCA 的遺傳機(jī)理尚不清楚，其調(diào)控機(jī)理有待于進(jìn)一步挖掘和驗(yàn)證。

通過(guò)對(duì)籽粒顯著關(guān)聯(lián)位點(diǎn)候選基因預(yù)測(cè)，4 個(gè)轉(zhuǎn)錄因子基因即C3H7（C3H-transcription factor 317）、HB81（homeobox-transcription factor 81）、E2F11（E2FDP-transcription factor 211）和E2F17以及鉀離子高親和性轉(zhuǎn)運(yùn)蛋白基因HAK25（potassium high-affinity transporter25）、核糖體蛋白基因RPL23（60S ribosomal protein L23-like）和 質(zhì) 膜ATP 酶 基 因PM4（plasma membrane ATPase4）等可能是調(diào)控籽粒性狀GCA 的關(guān)鍵基因（附表1）。DiZF-C3H1是一個(gè)鋅指結(jié)構(gòu)轉(zhuǎn)錄因子，其過(guò)表達(dá)會(huì)導(dǎo)致擬南芥轉(zhuǎn)基因植株根系發(fā)育遲緩、花期延遲、花器官異常和角果變形［47］。轉(zhuǎn)錄因子VvHB58通過(guò)多種激素信號(hào)通路影響葡萄種子和果實(shí)的發(fā)育，其在轉(zhuǎn)基因番茄中的過(guò)表達(dá)會(huì)導(dǎo)致頂端優(yōu)勢(shì)喪失、果實(shí)大小和種子數(shù)量減少以及胚乳細(xì)胞變大［48］。E2F 轉(zhuǎn)錄因子家族在玉米種子萌發(fā)發(fā)育過(guò)程中起著重要作用［49］。敲除OsHAK26或質(zhì)膜定位的H+-K+轉(zhuǎn)運(yùn)體基因OsHAK1或OsHAK5，均可引起水稻花粉花藥畸變、花粉粒數(shù)減少以及花粉粒萌發(fā)率降低，而且體外試驗(yàn)結(jié)果表明，外源K+濃度對(duì)oshak5花粉萌發(fā)率有顯著影響［50］。此外，粒寬GCA 候選基因Zm00001d046560 被注釋為質(zhì)膜ATP酶PM4。該基因家族通過(guò)參與生長(zhǎng)素介導(dǎo)的細(xì)胞延長(zhǎng)過(guò)程進(jìn)而影響小麥胚乳發(fā)育［51］。核糖體蛋白在玉米籽粒發(fā)育中起著重要作用。Wang 等［52］通過(guò)對(duì)轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行共表達(dá)網(wǎng)絡(luò)分析，發(fā)現(xiàn)38個(gè)核糖體蛋白基因?yàn)樽蚜Ｐ誀铒@著關(guān)聯(lián)模塊的核心基因。玉米Dek44編碼線粒體核糖體蛋白R(shí)PL9，其突變體dek44籽粒發(fā)育遲緩、具有小粒表型［53］。

3.4 顯著位點(diǎn)作為固定效應(yīng)有助于提高一般配合力基因組預(yù)測(cè)精度

由于穗行數(shù)GCA效應(yīng)的遺傳力較高，即使利用5種隨機(jī)效應(yīng)模型也獲得了較高的預(yù)測(cè)準(zhǔn)確性（0.62～0.74）。多數(shù)研究表明，玉米自交系和雜交種穗行數(shù)性狀本身具有較高的預(yù)測(cè)力，適合用來(lái)開(kāi)展基因組預(yù)測(cè)［29，54-55］。粒寬和穗行數(shù)利用隨機(jī)效應(yīng)模型也獲得了較高的預(yù)測(cè)準(zhǔn)確性（0.63～0.73）。盡管粒長(zhǎng)GCA 的遺傳力較高，但隨機(jī)模型對(duì)其預(yù)測(cè)的準(zhǔn)確性較低（0.28～0.45）。由此可見(jiàn)，穗行數(shù)和粒寬GCA 預(yù)測(cè)精度主要由其遺傳力決定，而粒長(zhǎng)GCA 除了受遺傳力影響外，還可能主要受性狀遺傳結(jié)構(gòu)等因素的影響。

本研究結(jié)果表明，多數(shù)情況下，利用7 種MGWAS挖掘的顯著關(guān)聯(lián)位點(diǎn)作為固定效應(yīng)加入預(yù)測(cè)方法可進(jìn)一步提高穗行數(shù)和2 個(gè)籽粒性狀GCA 的預(yù)測(cè)準(zhǔn)確性。An 等［29］認(rèn)為將多種MGWAS 方法挖掘的tagSNP 作為固定效應(yīng)可以顯著提高穗行數(shù)本身的預(yù)測(cè)力。相比隨機(jī)效應(yīng)模型，Ma 等［55］研究發(fā)現(xiàn)將4～8 個(gè)顯著SNPs 作為固定效應(yīng)整合到GBLUP、RKHS、Bayes A、Bayes B 和Bayes C 中可以提高玉米穗行數(shù)和穗長(zhǎng)基因組估計(jì)育種值的準(zhǔn)確性，但這種策略并不能提高其他性狀例如單株產(chǎn)量、單穗重和百粒重的預(yù)測(cè)力。在玉米巢式關(guān)聯(lián)群體中，整合最顯著關(guān)聯(lián)位點(diǎn)作為固定效應(yīng)加入預(yù)測(cè)模型有助于提高株高、小斑病抗性和灰斑病抗性的預(yù)測(cè)精度［56］。在小麥中，整合共同關(guān)聯(lián)性較強(qiáng)的GWAS 信號(hào)作為固定效應(yīng)加入GBLUP 模型中可提高9%～10%產(chǎn)量的基因組預(yù)測(cè)力，而穩(wěn)定性指數(shù)Pi 沒(méi)有表現(xiàn)出優(yōu)勢(shì)［57］。通過(guò)模擬研究，Rice 等［58］發(fā)現(xiàn)絕大多數(shù)情況下固定效應(yīng)作為協(xié)方差加入預(yù)測(cè)模型并不能提高預(yù)測(cè)準(zhǔn)確性。在雙親群體中，Herter 等［59］指出將顯著的QTL 作為固定效應(yīng)可以提高小麥抽穗期、株高、赤霉病和葉枯病的預(yù)測(cè)準(zhǔn)確性。在玉米BC3F5群體中，相比隨機(jī)效應(yīng)模型，將1～2 個(gè)QTL 作為固定效應(yīng)可以提高穗行數(shù)、穗粒數(shù)、行粒數(shù)和葉夾角基因組估計(jì)育種值的準(zhǔn)確性，而將1～5 個(gè)QTL 作為固定效應(yīng)會(huì)降低穗位高的預(yù)測(cè)力［60］。

除了性狀、群體類型的影響外，顯著位點(diǎn)作為固定效應(yīng)能否提高預(yù)測(cè)能力也取決于不同的預(yù)測(cè)模型和顯著位點(diǎn)的個(gè)數(shù)。例如，玉米BC3F5群體中，Bayes A和RKHS方法整合1～5 個(gè)固定效應(yīng)QTL 能夠提高株高的預(yù)測(cè)力，但整合6 個(gè)固定效應(yīng)QTL 會(huì)降低預(yù)測(cè)力；而B(niǎo)ayes C、Bayesian LASSO、Bayesian ridge regression 和GBLUP 僅整合5個(gè)固定效應(yīng)QTL不會(huì)降低株高的預(yù)測(cè)準(zhǔn)確性［60］。本研究中，GBLUP 和RKHS 利用相同MGWAS 信息在不同環(huán)境表現(xiàn)一致。盡管大部分鑒定的顯著SNPs 作為固定效應(yīng)能夠有效地提高穗行數(shù)和2 個(gè)籽粒性狀GCA 的預(yù)測(cè)準(zhǔn)確性，但是各個(gè)MGWAS 方法在不同環(huán)境表現(xiàn)不一。這是由于不同MGWAS 檢測(cè)功效不同導(dǎo)致了位點(diǎn)信息的差異。綜上所述，利用多種MGWAS方法檢測(cè)顯著位點(diǎn)，不僅可以為后續(xù)基因功能驗(yàn)證提供可靠基因信息，還能進(jìn)一步對(duì)分子標(biāo)記信息進(jìn)行優(yōu)化，提高GCA 基因組估計(jì)育種值的準(zhǔn)確性，為選擇高配合力親本材料提供準(zhǔn)確的預(yù)測(cè)信息。

4 結(jié)論

利用7種MGWAS方法共鑒定46個(gè)控制玉米穗行數(shù)及籽粒性狀GCA顯著位點(diǎn)，其中6個(gè)至少被2種MGWAS方法和2 個(gè)環(huán)境重復(fù)檢測(cè)到。共挖掘到59 個(gè)候選基因，其中EMB12、TPR、PPR和APX2是控制穗行數(shù)GCA的關(guān)鍵基因，而UKL2、E2F11、E2F17、C3H17、HB81、HAK25、RPL23和PM4是粒長(zhǎng)和粒寬GCA 重要的候選基因。5 種隨機(jī)效應(yīng)模型可以有效預(yù)測(cè)玉米穗行數(shù)和粒寬GCA效應(yīng)。不同MGWAS顯著位點(diǎn)作為固定效應(yīng)加入GS預(yù)測(cè)方法有利于獲得更高的預(yù)測(cè)精度，尤其適用于粒長(zhǎng)GCA的預(yù)測(cè)。