高粱千粒重全基因組關聯分析和候選基因預測

鄒桂花 丁延慶 徐建霞 曹 寧 陳合云 劉合芹 鄭學強 張立異,*

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

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

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

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

1 材料與方法

1.1 試驗材料

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

1.2 田間試驗與表型調查

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

1.3 基因型檢測

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

1.4 全基因組關聯分析

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

1.5 候選基因預測

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

1.6 候選基因的轉錄本表達水平差異分析

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

1.7 候選基因的單倍型分析

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

2 結果與分析

2.1 千粒重表型統計分析

242份高粱種質資源千粒重表型在7個環境中均存在廣泛的遺傳變異。單個環境下測得的千粒重最大值和最小值的差值(極值)范圍在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年杭州。同一環境不同年份極值存在明顯差異，2018年杭州與2019年杭州極值相差3.29 g，2019年貴陽與2018年貴陽極值相差3.66 g。單個環境千粒重平均值范圍為21.87～26.89 g，變異系數范圍為0.19～0.25(表1)。

表1 3年7個環境關聯群體千粒重的統計分析

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

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

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

注：** 表示在P<0.01水平相關性顯著。

2.2 千粒重全基因組關聯分析

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

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

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

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

2.3 千粒重候選基因分析

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

2.4 千粒重候選基因的轉錄本表達水平差異分析

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

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

2.5 候選基因單倍型分析

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

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

3 討論

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

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

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

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

4 結論

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