高粱千粒重全基因組關(guān)聯(lián)分析和候選基因預(yù)測

鄒桂花 丁延慶 徐建霞 曹 寧 陳合云 劉合芹 鄭學(xué)強(qiáng) 張立異,*

(1 浙江省農(nóng)業(yè)科學(xué)院作物與核技術(shù)利用研究所/浙江省數(shù)字旱糧重點實驗室，浙江 杭州 310021；2 貴州省農(nóng)業(yè)科學(xué)院貴州省旱糧研究所，貴州 貴陽 550006)

高粱[Sorghumbicolor(L.)Moench]是我國最早栽培的禾谷類作物之一，因其具有很強(qiáng)的抗逆性和適應(yīng)性，在全國各地均有種植。歷史資料表明，高粱在解決糧食問題中曾發(fā)揮過巨大作用[1]。隨著水稻、小麥、玉米等糧食作物育種水平的提高以及人民生活條件的改善，高粱的利用價值發(fā)生了根本性轉(zhuǎn)變，由曾經(jīng)的口糧變成了如今的工業(yè)原料，用于釀酒、飼用、制醋等，繼續(xù)在我國國民經(jīng)濟(jì)發(fā)展和鄉(xiāng)村振興中發(fā)揮重要作用[2]。

實踐證明，高粱籽粒大小的遺傳改良對單產(chǎn)水平的提高具有重要意義。自上世紀(jì)50年代以來，我國高粱主栽品種經(jīng)歷了從地方品種到新品種到優(yōu)良雜交種的重大更替。每一次品種的更替不僅使得單位面積產(chǎn)量大幅提高，而且導(dǎo)致選育的骨干親本籽粒也在不斷變大[3]。籽粒變大不僅能夠增加產(chǎn)量，而且對種子活力、播后出苗、早期生長和品質(zhì)性狀均具有重要影響[4-6]。因此，高粱籽粒大小的研究一直是國內(nèi)外高粱育種家和遺傳學(xué)家關(guān)注的熱點。早期研究主要從形態(tài)和生理學(xué)機(jī)制上揭示高粱籽粒大小、粒重和產(chǎn)量三者之間的關(guān)系，但上述研究均未與粒重性狀的遺傳機(jī)理結(jié)合，無法確定控制粒重的數(shù)量性狀位點(quantitative trait locus，QTLs)或基因，從而限制了高粱粒重性狀遺傳研究進(jìn)展。隨著高粱分子標(biāo)記技術(shù)的發(fā)展和全基因組測序的完成[7]，利用連鎖分析方法檢測高粱粒重主效QTLs成為熱門課題。大量研究報道了高粱粒重的QTLs定位[8-13]，但由于受作圖群體親本遺傳背景、群體類型和群體大小、標(biāo)記類型差異等的影響，檢測到影響粒重QTLs的數(shù)量和位置存在很大差異，且多數(shù)位點對粒重表型的貢獻(xiàn)率較低，通過QTL克隆的方法克隆這些影響粒重性狀的基因，并深入研究上述基因的功能還存在很大困難。到目前為止，通過QTLs圖位克隆的方法對高粱粒重相關(guān)基因的克隆和功能研究報道較少，遠(yuǎn)遠(yuǎn)落后于水稻、擬南芥和玉米粒重性狀的分子遺傳研究。

基于高通量測序技術(shù)的全基因組關(guān)聯(lián)分析(genome-wide association study，GWAS)方法的發(fā)展，為有效獲取目標(biāo)性狀的候選基因提供了新途徑。與傳統(tǒng)連鎖作圖方法相比，GWAS方法可同時對自然群體內(nèi)遺傳多樣性豐富的每個個體進(jìn)行高通量測序，在獲得數(shù)以百萬計的單核苷酸多態(tài)性(single-nucleotide polymorphism, SNP)分子標(biāo)記的基礎(chǔ)上，運用不斷發(fā)展的統(tǒng)計方法將目標(biāo)性狀的表型與基因型進(jìn)行關(guān)聯(lián)分析，從而快速準(zhǔn)確地定位到影響目標(biāo)性狀的染色體區(qū)段或基因。該方法因研究材料來源廣泛，可有效捕獲個體間豐富的變異，節(jié)省構(gòu)建遺傳群體所需時間，具有更高的分辨率[14]等特點，迅速成為傳統(tǒng)雙親種群作圖的替代方法和研究復(fù)雜性狀遺傳的有力工具。目前GWAS已廣泛應(yīng)用于動植物的功能基因挖掘，并取得了顯著的成績[15]。盡管近年來利用GWAS方法解析農(nóng)作物重要農(nóng)藝性狀的遺傳機(jī)理發(fā)展迅猛，但有關(guān)高粱千粒重全基因組關(guān)聯(lián)分析和候選基因的研究報道較少，僅少數(shù)研究組報道了高粱粒重性狀的GWAS分析[16-20]。這些研究所采用的方法部分為基于單位點分析，或不超過兩種模型，所得結(jié)果難以相互比較?；诖?，本研究以242份籽粒高粱品種(系)組成的自然群體為供試群體，連續(xù)3年在杭州、貴陽、樂東和陵水4個點7個環(huán)境下測定千粒重，利用多位點GWAS分析軟件中的6種模型進(jìn)行全基因組關(guān)聯(lián)分析，鑒定影響千粒重性狀的顯著性位點，挖掘、預(yù)測相關(guān)候選基因，旨在為進(jìn)一步開展高粱千粒重基因的克隆及其遺傳改良提供依據(jù)。

1 材料與方法

1.1 試驗材料

供試材料為242份籽粒高粱種質(zhì)資源，主要為我國的地方品種、育成品種及高世代育種材料，來自中國農(nóng)業(yè)科學(xué)院作物科學(xué)研究所作物種質(zhì)資源中心、浙江省農(nóng)業(yè)科學(xué)院作物與核技術(shù)利用研究所以及貴州省旱糧研究所。

1.2 田間試驗與表型調(diào)查

2018—2020年連續(xù)3年分別在貴州貴陽(18GY、19GY和20GY)、浙江杭州(18HZ和19HZ)、海南樂東(20LD)和海南陵水(20LS)共進(jìn)行7次試驗。小區(qū)長2.5 m，行距60 cm，株距25 cm，人工點播，適時進(jìn)行間苗定苗、追肥、除草和病蟲害防治等管理工作。在抽穗期，每小區(qū)選擇具有代表性的3個穗子用硫酸鈉紙袋套袋自交，在整株完成授粉時，將紙袋換成網(wǎng)袋，直至種子成熟。收獲后自然曬干，再脫粒、精選和考種分析。千粒重具體測量方法：單株全株脫粒后，選擇飽滿無損傷的籽粒(>100粒)，利用SC-G型自動考種分析及千粒重儀系統(tǒng)(杭州萬深檢測科技有限公司)進(jìn)行計數(shù)和測定，重復(fù)取樣3次，保證3次稱重結(jié)果差異小于0.1 g(5%以內(nèi))為有效的3個數(shù)據(jù)，取平均值。利用SPSS 21.0軟件對各個環(huán)境的千粒重表型進(jìn)行平均值、標(biāo)準(zhǔn)差、頻率分布等描述性統(tǒng)計分析和相關(guān)性分析，利用R 4.03軟件繪圖。

1.3 基因型檢測

采用CTAB法[21]提取242份高粱資源葉片DNA，質(zhì)檢合格后利用HiSeq 4000高通量測序儀(Illumina，美國)進(jìn)行測序。利用BWA軟件將過濾后的測序數(shù)據(jù)比對到高粱參考基因組BTx623的序列上(https://phytozome-next.jgi.doe.gov/info/Sbicolor_v3_1_1)。利用GATK軟件進(jìn)行群體SNP檢測，用Plink軟件進(jìn)行SNP過濾，篩選條件：SNP最小等位基因頻率(minimum allele frequency，MAF)>5%，缺失率(missing rate, MR)<20%。最終獲得2 015 850個高質(zhì)量SNPs用于后續(xù)全基因組關(guān)聯(lián)分析和單倍型分析。

1.4 全基因組關(guān)聯(lián)分析

對2018—2020年貴陽、杭州、樂東和陵水7個環(huán)境獲得的千粒重性狀進(jìn)行全基因組關(guān)聯(lián)分析。為了提高關(guān)聯(lián)分析的準(zhǔn)確性和可靠性，本研究利用mrMLM 4.0軟件[22]進(jìn)行性狀多位點關(guān)聯(lián)分析(multi-locus GWAS methods)，即同一染色體上物理距離≤20 kb 的所有顯著關(guān)聯(lián)SNP位點總結(jié)為1個。該軟件包含mrMLM[23]、FASTmrMLM[24]、FASTmrEMMA[25]、ISIS EM-BLASSO[26]、pLARmEB[27]、pKWmEB[28]6種關(guān)聯(lián)模型。K矩陣(親緣關(guān)系矩陣)由mrMLM軟件進(jìn)行計算。GWAS分析時所有參數(shù)均設(shè)置為默認(rèn)值。顯著關(guān)聯(lián)的SNPs閾值為LOD(logarithm of odds)≥3.0[26]。當(dāng)顯著SNP位點(物理距離不超過50 kb)滿足至少在3個模型中檢測到，或者在2個以上試驗中檢測到，即認(rèn)為該SNP位點可靠，將其作為1個數(shù)量性狀核苷酸位點(quantitative trait nucleotides，QTN)[29]。QTN的命名參照Linge等[29]的方法，如qtnTGW1.1表示在第1染色體上控制千粒重的第一個QTN，利用RectChr 1.31軟件繪圖。

1.5 候選基因預(yù)測

根據(jù)已克隆的水稻(Oryzasativa)和玉米(ZeamaysL.)粒重相關(guān)基因的蛋白序列，在Phytozome 13數(shù)據(jù)庫中搜索高粱的同源基因，將搜索出來的同源基因的物理位置與本研究中檢測到的千粒重顯著關(guān)聯(lián)的QTNs所在的物理位置進(jìn)行比較，以QTNs峰值所在位點上下游50 kb范圍為界，若同源基因落于該范圍內(nèi)，則視為可能的候選基因。

1.6 候選基因的轉(zhuǎn)錄本表達(dá)水平差異分析

對小粒品種654和大粒品種LTR108在穗發(fā)育的3個階段進(jìn)行取樣，選取<5 cm幼穗、10～15 cm幼穗以及抽穗后5 d三個時期的穗樣品，提取總RNA[13]。按照NEBNext?UltraTMRNA文庫構(gòu)建試劑盒(New England Biolabs,美國)操作說明構(gòu)建轉(zhuǎn)錄組文庫，質(zhì)檢合格后在NovaSeq 6000平臺(Illumina，美國)進(jìn)行測序。測序完成后，得到每個樣品的原始數(shù)據(jù)(raw data)，經(jīng)過濾處理，通過軟件Filter-FQ去除接頭序列、未知序列超過10%的序列以及低質(zhì)量的序列后獲得清除數(shù)據(jù)(clean data)，利用Hisat2軟件將clean data錨定到高粱參考基因組BTx623上?；蜣D(zhuǎn)錄組表達(dá)水平是通過FPKM(fragments per kilobase of transcript per million mapped reads)方法進(jìn)行歸一化處理。每時期的樣品均取3次生物學(xué)重復(fù)。

1.7 候選基因的單倍型分析

基于242份高粱種質(zhì)資源的重測序數(shù)據(jù)，挑選出千粒重候選基因(Sobic.002G367300、Sobic.002G367600、Sobic.003G411900、Sobic.009G024600和Sobic.010G210100)的多態(tài)性SNP位點，利用CandiHap 1.30軟件對候選基因進(jìn)行單倍型分析，每個候選基因選擇前3種主要單倍型(≥10份材料)進(jìn)行千粒重差異顯著性分析。

2 結(jié)果與分析

2.1 千粒重表型統(tǒng)計分析

242份高粱種質(zhì)資源千粒重表型在7個環(huán)境中均存在廣泛的遺傳變異。單個環(huán)境下測得的千粒重最大值和最小值的差值(極值)范圍在23.64～31.96 g之間，其中2020年陵水的差異最小，然后依次是2018年貴陽(24.32 g)、2020年樂東(26.71 g)、2020年貴陽(27.07 g)、2019年貴陽(27.98 g)、2018年杭州(28.67 g)，差異最大為2019年杭州。同一環(huán)境不同年份極值存在明顯差異，2018年杭州與2019年杭州極值相差3.29 g，2019年貴陽與2018年貴陽極值相差3.66 g。單個環(huán)境千粒重平均值范圍為21.87～26.89 g，變異系數(shù)范圍為0.19～0.25(表1)。

表1 3年7個環(huán)境關(guān)聯(lián)群體千粒重的統(tǒng)計分析

從頻率分布圖(圖1)可以看出，242份高粱種質(zhì)資源的千粒重在不同年份呈現(xiàn)廣泛的連續(xù)變異。2018年貴陽和2021年陵水的千粒重分布更趨于平展，2020年貴陽的千粒重分布具有右側(cè)長尾。不同年份千粒重的峰度和偏度的絕對值均小于1(表1)，表明千粒重值在不同年份均服從正態(tài)分布，呈現(xiàn)出典型的數(shù)量性狀特點，受多基因控制。

注：A：2018年杭州；B：2018年貴陽；C：2019年杭州；D：2019年貴陽；E：2020年樂東；F：2020年陵水；G：2020年貴陽。

相關(guān)系數(shù)可以反映同一品種在不同環(huán)境中表現(xiàn)的一致性和可重復(fù)性。2018年杭州與2019年杭州兩個環(huán)境間的相關(guān)系數(shù)雖達(dá)到極顯著水平，但低于0.50。2018年貴陽與2019年貴陽以及2019年貴陽與2020年貴陽間的相關(guān)系均為0.71，而2018年貴陽與2020年貴陽的相關(guān)系數(shù)為0.60；其余環(huán)境間的相關(guān)系數(shù)均界于0.52～0.69之間(圖2)。表明相同地點，因年際間氣象因素的差異，同一品種的表現(xiàn)會出現(xiàn)差異。

注：** 表示在P<0.01水平相關(guān)性顯著。

2.2 千粒重全基因組關(guān)聯(lián)分析

利用mrMLM 4.0軟件包中的mrMLM、FASTmrMLM、FASTmrEMMA、ISIS EM-BLASSO、pLARmEB、pKWmEB 6種多位點關(guān)聯(lián)模型對7個環(huán)境的千粒重進(jìn)行全基因組關(guān)聯(lián)分析，共檢測到323個與千粒重顯著關(guān)聯(lián)的SNPs，分布在10條染色體上。LOD值范圍是3.0～35.8, 可解釋0.4%～26.6%的表型變異。比較不同模型檢測到的顯著性位點數(shù)目發(fā)現(xiàn)，F(xiàn)ASTmrMLM模型檢測到的最多(74個)，然后依次是pKWmEB模型(66個)、pLARmEB模型(62個)、mrMLM模型(52個)、ISIS EM-BLASSO模型(42個)、FASTmrEMMA模型(27個)。有4種模型(FASTmrMLM、ISIS EM-BLASSO、pLARmEB和pKWmEB)均在2018年杭州環(huán)境點檢測到最多的顯著性位點，另外2種模型(mrMLM和FASTmrEMMA)則在2019年貴陽環(huán)境點檢測到最多的顯著性位點。有4種模型(mrMLM、FASTmrMLM、FASTmrEMMA和pKWmEB)在2019年杭州環(huán)境點檢測到的顯著性位點數(shù)最少，另有ISIS EM-BLASSO和pLARmEB模型分別在2020年陵水和2020年貴陽環(huán)境點檢測到的顯著性位點數(shù)最少。同一模型檢測到的顯著性位點在各條染色體上的分布不均衡。mrMLM模型在7個環(huán)境中檢測到的顯著性位點主要分布在第2和第5染色體上，第8染色體上的最少，僅在2020年樂東環(huán)境點檢測到1個顯著性位點。FASTmrMLM模型在第3染色體上檢測到的顯著性位點最多，第6染色體上最少，僅在2018年杭州環(huán)境點檢測到1個顯著性位點。FASTmrEMMA模型在第2和第6染色體上各檢測到6個顯著性位點，在第1、第3和第9染色體上沒有檢測到顯著性位點。ISIS EM-BLASSO模型在第10染色體上檢測到的顯著性位點最多，其他染色體上分布著2～7個不等的顯著性位點。pLARmEB模型在第3染色體上檢測到的顯著性位點最多，其他染色體上分布著1～11個不等的顯著性位點。pKWmEB模型在第3和第5染色體上各檢測到13個顯著性位點，在第6染色體上顯著性位點數(shù)最少，僅1個。

分析上述顯著性位點，發(fā)現(xiàn)不同模型間存在重復(fù)和交疊結(jié)果，可以合并成同一個QTN，合并條件為：至少在3種模型中檢測到，或者在2個以上環(huán)境中檢測到相同的SNP，即視為一個可信的QTN。結(jié)果發(fā)現(xiàn)，共有96個顯著性QTNs位點與千粒重相關(guān)聯(lián)，這些位點不均勻地分布在10條染色體上，依次為6個、16個、15個、7個、16個、2個、10個、3個、7個和14個(表2、圖3)。同一種模型在不同環(huán)境中均檢測到的QTNs有7個，占總數(shù)的7.4%。在同一環(huán)境中，3種以上模型均檢測到的QTNs有48個，占51.1%。有6個QTNs在3個環(huán)境下被檢測到，單個QTN可解釋1.0%～7.1%的表型變異。qtnTGW3.2和qtnTGW9.1均在2019年貴陽、2020年貴陽和2020年陵水環(huán)境點被檢測到，qtnTGW3.13在2018年杭州、2020年貴陽和2020年陵水環(huán)境點被檢測到。qtnTGW4.7和qtnTGW5.15在2019年貴陽、2020年貴陽和2020年樂東環(huán)境點被檢測到。qtnTGW10.2在2019年杭州、2020年貴陽和2020年陵水環(huán)境點被檢測到。在2個環(huán)境下檢測到的QTNs有35個，占36.5%。值得注意的是，未發(fā)現(xiàn)一個模型在7個環(huán)境下均檢測到與千粒重顯著關(guān)聯(lián)的共同QTNs，也未發(fā)現(xiàn)6個模型均檢測到與千粒重顯著關(guān)聯(lián)的共同QTNs。

表2 利用mrMLM4.0軟件中三個以上模型或兩個環(huán)境中檢測到與千粒重顯著性關(guān)聯(lián)的QTNs

圖3 影響千粒重的QTNs在染色體上的分布

2.3 千粒重候選基因分析

根據(jù)已克隆的水稻和玉米粒重相關(guān)基因的蛋白質(zhì)序列，在高粱中共搜索到111個同源基因(https://www.ricedata.cn/；https://www.maizegdb.org/；https://phytozome-next.jgi.doe.gov/)。這些基因不均勻地分布在10條染色體上，其中第1染色體上數(shù)量最多，達(dá)23個，第5染色體上最少，僅1個[19]。將這些基因的物理位置與本研究中所檢測到的QTNs所在區(qū)域進(jìn)行比較，發(fā)現(xiàn)有5個與水稻粒重相關(guān)基因同源的高粱基因位于4個QTNs區(qū)間內(nèi)，其他QTNs區(qū)間內(nèi)均不包含已克隆的水稻粒重相關(guān)基因的同源基因。在2018年杭州和2020年陵水環(huán)境點同時檢測到的千粒重qtnTGW2.15，所在區(qū)間內(nèi)包含Sobic.002G367300和Sobic.002G367600基因：前者與水稻GW7/GL7同源；生物學(xué)功能顯示GL7編碼擬南芥LONGIFOLIA蛋白的同源蛋白，調(diào)節(jié)細(xì)胞的縱向伸長[30]；GW7表達(dá)量上調(diào)，能增加谷粒的縱向細(xì)胞分裂，并減少橫向細(xì)胞分裂，導(dǎo)致谷粒變得細(xì)長。后者與水稻BG2同源；與野生型相比，bg2-D突變體增高，同時粒長、粒寬、粒厚以及千粒重均顯著增加[31]。在2018年貴陽環(huán)境點利用mrMLM、FASTmrMLM、pLARmEB、pKWmEB和ISIS EM-BLASSO 5種模型均檢測到的千粒重qtnTGW3.13，所在區(qū)間內(nèi)包含Sobic.003G411900基因，與水稻OsARF4同源；生物學(xué)功能顯示OsARF4為生長素應(yīng)答因子，與OsGSK5/OsSK41互作負(fù)向調(diào)控籽粒大小和粒重[32]。在2019年貴陽、2020年貴陽和2020年陵水環(huán)境點同時檢測到的千粒重qtnTGW9.1，所在區(qū)間內(nèi)包含Sobic.009G024600基因，與水稻RSR1同源；生物學(xué)功能顯示RSR1是調(diào)控淀粉合成的AP2/EREBP家族的轉(zhuǎn)錄因子；淀粉作為水稻胚乳的主要組成成分，其含量和結(jié)構(gòu)直接影響著水稻種子的品質(zhì)[33]。在2018年貴陽環(huán)境點利用FASTmrMLM、pLARmEB和pKWmEB 3種模型均檢測到千粒重qtnTGW10.11，所在區(qū)間內(nèi)包含Sobic.010G210100基因，與水稻TGW6同源；生物學(xué)功能顯示TGW6不僅直接控制胚乳的長度，而且間接參與源到庫的碳水化合物的運輸。在庫器官(籽粒)中，日本晴的tgw6等位基因能通過控制IAA供應(yīng)影響合子體到細(xì)胞化階段的轉(zhuǎn)變，從而限制細(xì)胞數(shù)目的增加和籽粒長度，而Kasalath中tgw6等位基因功能喪失會通過影響源器官的生長增加粒重，從而導(dǎo)致水稻產(chǎn)量顯著增加[34]。

2.4 千粒重候選基因的轉(zhuǎn)錄本表達(dá)水平差異分析

圖4為高粱小粒品種654和大粒品種LTR108穗發(fā)育時期5個候選基因的轉(zhuǎn)錄組表達(dá)量差異分析，可以看出5個候選基因的表達(dá)模式存在明顯差異。其中Sobic.002G367300基因在小粒品種654三個穗發(fā)育時期表達(dá)量緩慢升高，而大粒品種LTR108在前兩個時期表達(dá)量幾乎沒有變化，僅在抽穗后5 d表達(dá)量明顯下降。在前兩個時期小粒品種654的表達(dá)量低于大粒品種LTR108，而抽穗后5 d的表達(dá)量顯著高于大粒品種LTR108。Sobic.002G367600基因在小粒品種654的表達(dá)量隨著穗發(fā)育的完成逐漸降低，而大粒品種LTR108的表達(dá)量呈現(xiàn)高-低-高的趨勢，且在<5 cm幼穗時期的表達(dá)量最高。大粒品種LTR108的表達(dá)量在三個時期均高于小粒品種654，且在<5 cm幼穗時期和抽穗后5 d的表達(dá)量是小粒品種654的3倍以上。Sobic.003G411900基因在兩個品種中的表達(dá)量趨勢呈現(xiàn)一致，均為隨著穗發(fā)育的完成，表達(dá)量逐漸降低。三個時期小粒品種654的表達(dá)量均高于大粒品種LTR108。Sobic.009G024600基因在兩個品種中表達(dá)量的趨勢也呈現(xiàn)一致，均為前兩個時期低，在抽穗后5 d表達(dá)量明顯升高。Sobic.010G210100基因在兩個品種中的表達(dá)量均較低(FPKM<1)，該基因在小粒品種654中三個時期的表達(dá)趨勢為高-低-高，而在大粒品種LTR108中的表達(dá)趨勢是低-高-低。

注：S1: < 5cm幼穗；S2: 10～15cm幼穗；S3: 抽穗后5 d的穗子。

2.5 候選基因單倍型分析

為了調(diào)查候選基因中變異SNP位點對表型的影響，本試驗開展了Sobic.002G367300、Sobic.002G367600、Sobic.003G411900、Sobic.009G024600和Sobic.010G210100等5個基因的包含10個及以上品種的單倍型分析，如圖5所示。候選基因Sobic.002G367300在第1個外顯子(包含3個非同義突變和1個同義突變)、第1個內(nèi)含子(包含13個多態(tài)性SNPs)、第2個外顯子(包含1個同義突變)和第5個外顯子(包含2個非同義突變)共存在20個多態(tài)性SNPs，這些多態(tài)性SNPs將110個籽粒高粱劃分為Hap1、Hap2和Hap3，千粒重分別為24.63、23.35和25.93 g。在該候選基因中，Hap3千粒重最大，與Hap2差異達(dá)顯著水平(P<0.05)，為優(yōu)勢單倍型。候選基因Sobic.002G367600在其5′-UTR非編碼區(qū)(包含1個多態(tài)性SNP)和第1個外顯子(包含2個非同義突變)共存在3個多態(tài)性SNPs，這些多態(tài)性SNPs將188個籽粒高粱劃分為Hap1、Hap2和Hap3三種主要單倍型，千粒重分別為23.77、24.77和25.14 g，3種主要單倍型間的差異未達(dá)到顯著水平。候選基因Sobic.003G411900在第2、第6、第7、第14個內(nèi)含子，第8個外顯子和3′-UTR共存在6個多態(tài)性SNPs，這些多態(tài)性SNPs將98個籽粒高粱劃分為Hap1、Hap2和Hap3三種主要單倍型，千粒重分別為23.61、25.44和24.75 g，3種主要單倍型間的差異也未達(dá)到顯著或極顯著水平。候選基因Sobic.010G210100在其5′-UTR非編碼區(qū)(包含1個多態(tài)性SNP)、第1個內(nèi)含子(包含2個多態(tài)性SNPs)、第3個外顯子(包含1個多態(tài)性SNP)和3′-UTR(包含1個多態(tài)性SNP)共存在5個多態(tài)性SNPs，這些多態(tài)性SNPs將171個籽粒高粱劃分為Hap1、Hap2和Hap3三種主要單倍型，千粒重分別為23.46、25.53和22.49 g，Hap1和Hap2之間千粒重的差異達(dá)極顯著水平(P<0.01)，Hap1和Hap3之間千粒重的差異不顯著，Hap2和Hap3之間千粒重的差異達(dá)顯著水平。候選基因Sobic.009G024600的內(nèi)含子和外顯子未發(fā)現(xiàn)差異SNP位點，故未進(jìn)行單倍型分類。

圖5 候選基因單倍型分組及每種單倍型SNP 組成

3 討論

作物粒重是典型的由微效多基因控制的數(shù)量性狀，其遺傳基礎(chǔ)可以通過基于連鎖作圖的QTL定位和基于連鎖不平衡的GWAS分析進(jìn)行闡明。為提高QTL定位和GWAS分析的準(zhǔn)確性和后續(xù)開展基因克隆的成功率，在不同遺傳背景和多個環(huán)境下調(diào)查粒重性狀表型，開展QTL定位和GWAS分析，鑒定穩(wěn)定的主效位點是前提。為闡明高粱粒重遺傳基礎(chǔ)，本研究利用同一套高粱自然群體在4個試驗點不同年份調(diào)查千粒重表型數(shù)據(jù)并開展GWAS分析。7個環(huán)境千粒重的最大值、最小值、變幅、頻率分布圖均存在明顯差異。同一個地點不同年份之間、不同地點相同年份之間，以及在不同氣候類型(亞熱帶濕潤溫和型氣候的貴陽、亞熱帶季風(fēng)氣候的杭州以及熱帶季風(fēng)氣候的陵水和樂東)條件下的表現(xiàn)都不完全一致。GWAS分析結(jié)果同樣反映了這種差異：7個環(huán)境下未檢測到相同的QTN，6個QTNs在3個環(huán)境下檢測到，2個環(huán)境下檢測到的QTNs有35個，在單一環(huán)境下檢測到的QTNs占多數(shù)，表明高粱粒重的遺傳基礎(chǔ)非常復(fù)雜。本研究檢測到的影響千粒重的QTNs，為探究高粱千粒重的遺傳提供了有用信息。

GWAS是基于連鎖不平衡來識別分子標(biāo)記之間或候選基因與性狀之間關(guān)系的方法。如果兩個基因座位間的不平衡是來源于群體結(jié)構(gòu)或親緣關(guān)系，則會影響GWAS結(jié)果的準(zhǔn)確性。此外，GWAS研究中僅能檢測到一小部分變異，對于由稀有變異、結(jié)構(gòu)變異、上位性互作及基因與環(huán)境互作等因素造成的遺傳變異難以檢測到[35]。為了提高GWAS分析的準(zhǔn)確性，需進(jìn)一步完善統(tǒng)計分析方法。基于混合線性模型(mixed linear model，MLM)的GWAS因兼顧固定效應(yīng)和隨機(jī)效應(yīng)，顯著提高了數(shù)量性狀位點檢測的功效，被視為標(biāo)準(zhǔn)方法[36]。然而，MLM是單位點模型的分析方法，一次只能進(jìn)行單個位點的檢測，與受多基因控制的復(fù)雜性狀的真實遺傳模型不符[37]。隨著GWAS模型的不斷發(fā)展，同時兼顧所有位點信息的多位點GWAS方法得以開發(fā)[23]。與單位點全基因組關(guān)聯(lián)分析模型相比，多位點SNP隨機(jī)效應(yīng)混合線性模型更加接近真實的動植物遺傳模型[38]，即使應(yīng)用了較為寬松的關(guān)聯(lián)閾值，也具有更高的統(tǒng)計功效和更低的假陽性率[39]。本研究在高粱中首次利用包含mrMLM、FASTmrMLM、FASTmrEMMA、pLARmEB、pKWmEB、ISIS EM-BLASSO 6種多位點模型的mrMLM 4.0軟件檢測影響千粒重的QTN，共檢測到323個顯著性位點。不同模型檢測到的顯著性位點存在差異，本研究中FASTmrEMMA模型檢測到顯著關(guān)聯(lián)位點最少，F(xiàn)ASTmrMLM模型檢測到的最多，在開展GWAS分析時，選擇不同的統(tǒng)計分析方法會影響檢測效果，因此，進(jìn)行多種統(tǒng)計分析方法的比較和甄選，將有助于提高GWAS分析的準(zhǔn)確性。本研究中最終合并和整合到96個影響千粒重的QTNs，遠(yuǎn)遠(yuǎn)多于以前單個研究報道的高粱粒重GWAS分析的數(shù)量，檢測效率大大提高。

到目前為止，高粱中僅少數(shù)文獻(xiàn)報道了粒重GWAS分析結(jié)果。其中Upadhyaya等[16]利用39個微衛(wèi)星標(biāo)記(simple sequence repeat，SSR)，在第10染色體的50 197 002-51 527 907 bp區(qū)間內(nèi)檢測到一個影響粒重的主效QTL；Sakhi等[17]利用98個SSR標(biāo)記，對107個高粱種質(zhì)采用MLM中的條件型Q+K模型在第8染色體的57 579 092-58 285 925 bp檢測到一個影響百粒重的主效QTL。兩個研究因標(biāo)記數(shù)量少，基因組上的覆蓋度低，檢測到的QTL數(shù)量極其有限。Boyles等[18]利用268 830個SNPs、390個品種組成的關(guān)聯(lián)群體以及混合線性模型(MLM)和貝葉斯稀疏線性混合模型(Bayesian sparse linear mixed model，BSLMM)對千粒重性狀開展GWAS分析，兩年環(huán)境共檢測到19個顯著關(guān)聯(lián)的SNPs。Tao等[19-20]利用同一套材料對千粒重開展GWAS分析，在2020年的報道中利用FarmCPU模型在自然群體中檢測到58個影響千粒重的顯著關(guān)聯(lián)位點，巢式關(guān)聯(lián)群體(backcross-nested association mapping，BC-NAM)群體中檢測到59個影響千粒重的顯著性關(guān)聯(lián)位點，有42個位點是兩個群體中同時被檢測到的；2021年的報道中利用始花時去除半個穗子處理的方法對千粒重開展GWAS分析，在自然群體中檢測到79個影響千粒重的顯著關(guān)聯(lián)位點，BC-NAM群體中檢測到69個影響千粒重的顯著性關(guān)聯(lián)位點，有57個位點是兩個群體中同時檢測到的。去掉半個穗子和保留完整穗子檢測到的共同顯著性位點有18個。將本研究與上述5個研究結(jié)果進(jìn)行比較，發(fā)現(xiàn)所檢測到的眾多顯著性位點中，僅qtnTGW2.12、qtnTGW2.13、qtnTGW3.13、qtnTGW4.6分別與Tao等[19-20]檢測到的qGS2.5、qGS2.8、qGS3.7和qFHGS4.11重疊。推測可能由于GWAS采用的分析群體不同，本研究分析群體主要由中國高粱種質(zhì)組成，而上述關(guān)聯(lián)群體所涉及的品種主要來源于非洲和美國的高粱。

水稻、玉米和高粱同屬于禾本科二倍體作物。在漫長的馴化過程中，一些馴化性狀(如果實增大、種子不易落粒、休眠被打破、生育期同步等)被人們有意或無意選擇并保存下來，而涉及這些相似性狀的基因因存在同源性和共線性，也被很好地保存下來，因而可能存在功能保守的基因組區(qū)域。水稻中克隆的粒重GS3在玉米和高粱中的同源基因均被證實影響著籽粒大小和粒重[13,40]。本研究檢測到4個千粒重QTNs，其所在區(qū)間內(nèi)含5個水稻中已經(jīng)克隆的粒重同源基因，這些基因在小粒品種654和大粒品種LTR108穗發(fā)育的三個時期表達(dá)量存在差異，有可能調(diào)控高粱千粒重，但需要進(jìn)一步驗證相關(guān)基因的功能。

4 結(jié)論

本研究利用多位點全基因組關(guān)聯(lián)分析方法中的6種模型對高粱千粒重性狀進(jìn)行了全基因組關(guān)聯(lián)分析，在7個環(huán)境下共檢測到96個與千粒重顯著關(guān)聯(lián)的QTNs位點，其中6個位點在3個環(huán)境中被同時檢測到，另外35個在2個環(huán)境中被同時檢測到。有5個與水稻粒重相關(guān)基因同源的高粱基因位于在4個QTNs區(qū)間內(nèi)，這些基因可能是調(diào)控籽粒發(fā)育的候選基因。