謝曉宇，王凱鴻，秦曉曉，王彩香?，史春輝，寧新柱，楊永林，秦江鴻，李朝周，馬麒?，宿俊吉?

1甘肅農(nóng)業(yè)大學(xué)生命科學(xué)技術(shù)學(xué)院/省部共建干旱生境作物學(xué)國家重點(diǎn)實(shí)驗(yàn)室，蘭州 730070；2新疆農(nóng)墾科學(xué)院棉花研究所，新疆石河子 832000；3石河子農(nóng)業(yè)科學(xué)研究院，新疆石河子 832000

0 引言

【研究意義】棉花(Gossypium spp.)是一種生產(chǎn)天然纖維的重要經(jīng)濟(jì)作物，在中國國民經(jīng)濟(jì)中占有重要地位。西北內(nèi)陸是中國棉花的主產(chǎn)區(qū)，但該地區(qū)有效積溫相對(duì)較低、無霜期相對(duì)較短，更適于早熟棉花品種的種植。因此，在西北內(nèi)陸棉區(qū)，棉花的早熟性顯得至關(guān)重要。吐絮率是指吐絮周期內(nèi)某個(gè)特定時(shí)期吐絮棉鈴數(shù)占總棉鈴數(shù)的比率，是衡量棉花吐絮集中性的重要指標(biāo)，也是棉花早熟性的一個(gè)重要性狀。因此，深入解析陸地棉吐絮率的遺傳基礎(chǔ)，鑒定其QTL等位基因，以及挖掘與目標(biāo)性狀相關(guān)的候選基因，對(duì)于培育早熟和適宜機(jī)采棉花新品種等方面具有重要意義。【前人研究進(jìn)展】棉花早熟性是一個(gè)綜合性狀，其相關(guān)性狀主要包括霜前花率、果枝始節(jié)位、果枝始節(jié)高、現(xiàn)蕾期、開花期、花鈴期和全生育期等[1-6]。棉花早熟性是受多基因控制的數(shù)量性狀，存在加性和顯性效應(yīng)[7]。相對(duì)于產(chǎn)量和纖維品質(zhì)性狀，有關(guān)棉花早熟性狀的研究相對(duì)較少。人們利用連鎖分析的方法，采用簡單重復(fù)序列（simple sequence repeats，SSR）標(biāo)記和單核苷酸多態(tài)性（single nucleotide polymorphism，SNP）標(biāo)記定位到棉花全生育期、開花期等早熟相關(guān)性狀的數(shù)量性狀位點(diǎn)（quantitative trait locus，QTL）[8-11]。隨著高通量測序技術(shù)的快速發(fā)展，全基因組關(guān)聯(lián)分析（genome-wide association studies，GWAS）已成為檢測QTL的另一種重要方法，以遺傳變異更豐富的自然群體為研究材料，基于連鎖不平衡的關(guān)聯(lián)作圖方法，大大地提高了QTL的檢測效率[12]。GWAS已被廣泛用于水稻[13]、玉米[14]、油菜[15]、大豆[16]和棉花[17-18]等作物育種目標(biāo)性狀定位的研究中。目前，已報(bào)道棉花GWAS研究主要集中在纖維品質(zhì)和產(chǎn)量等性狀上[19-23]，關(guān)于早熟性GWAS的研究相對(duì)較少[18,24-25]。SU等[24]利用簡化基因組測序開發(fā)的SNP標(biāo)記，開展了6個(gè)早熟相關(guān)性狀的GWAS研究，鑒定出8個(gè)SNP位點(diǎn)與5個(gè)早熟相關(guān)性狀之間存在顯著關(guān)聯(lián)，在與多個(gè)早熟相關(guān)性狀同時(shí)顯著關(guān)聯(lián)的1個(gè)SNP標(biāo)記附近發(fā)現(xiàn)一個(gè)調(diào)控開花的候選基因 CotAD_01947。LI等[25]利用 80K SNP基因芯片對(duì)169份陸地棉的早熟性狀進(jìn)行了大量的研究，獲得29個(gè)顯著的SNP位點(diǎn)，確定了2個(gè)潛在的改良棉花早熟性的候選基因 Gh_D01G0340和Gh_D01G0341。LI等[18]利用基于全基因組重測序GWAS鑒定到了與7個(gè)早熟相關(guān)性狀顯著關(guān)聯(lián)的307個(gè)SNP位點(diǎn)及1個(gè)候選基因Ghir_D03G011310。近年來，隨著勞動(dòng)力資源的緊缺和棉花人工采摘成本的大幅上漲，機(jī)械采收已成為西北內(nèi)陸棉花收獲的主要方式，而機(jī)械采收前需要對(duì)棉花進(jìn)行脫葉催熟處理。然而，脫葉催熟技術(shù)對(duì)棉花產(chǎn)量和纖維品質(zhì)有一定的影響。為了減小化學(xué)脫葉對(duì)棉花產(chǎn)量和纖維品質(zhì)的影響，西北內(nèi)陸棉區(qū)生產(chǎn)上一般要求9月中旬棉田整體吐絮率達(dá)到50%以上，才能噴施脫葉催熟劑[26-27]，因此，9月中旬棉花吐絮率已成為評(píng)價(jià)品種是否適宜機(jī)械采收的關(guān)鍵性狀之一。然而，目前對(duì)陸地棉吐絮率的遺傳解析及QTL定位的研究很少。李俊文等[28]基于連鎖作圖方法，利用陸地棉重組近交系（recombinant inbred lines，RIL）和永久 F2群體，對(duì)控制陸地棉吐絮性狀進(jìn)行 QTL檢測，共檢測到了 11個(gè)與吐絮率相關(guān)的QTL。目前，GWAS方法已被廣泛應(yīng)用于農(nóng)作物育種目標(biāo)性狀的遺傳解析，但是有關(guān)利用GWAS挖掘陸地棉吐絮率QTL的研究鮮有報(bào)道。而且GWAS的混合線性模型（mixed linear model，MLM）不能檢測自然群體中廣泛存在的復(fù)等位變異。同時(shí)進(jìn)行多重測驗(yàn)矯正時(shí)，通常利用提高顯著水平使得該模型僅能檢測到少數(shù)關(guān)聯(lián)位點(diǎn)，無法較為全面地檢測到QTL位點(diǎn)。針對(duì)目前常規(guī)MLM-GWAS方法的局限性，HE等[29]研究提出了一種限制性兩階段多位點(diǎn)全基因組關(guān)聯(lián)分析（restricted two-stage multi-locus GWAS，RTM-GWAS）方法，該方法是一種新發(fā)現(xiàn)的可以全面檢測自然群體和雙（多）親衍生群體中具有不同復(fù)等位變異QTL的分析方法，它解決了常規(guī)GWAS方法不能估計(jì)復(fù)等位變異、不能充分檢測出全基因組QTL等問題，大大提高了QTL的檢測效率[30]。RTM-GWAS方法已成功應(yīng)用于大豆百粒重[31-32]、油酸含量[33]、蛋白質(zhì)含量[34]和異黃酮含量等[35]，棉花產(chǎn)量[36]、纖維品質(zhì)[23]等數(shù)量性狀遺傳基礎(chǔ)的解析?！颈狙芯壳腥朦c(diǎn)】目前，西北棉區(qū)噴施脫葉劑是棉花機(jī)械采收的重要環(huán)節(jié)之一，且吐絮率已成為噴施脫葉劑之前一個(gè)很重要的評(píng)價(jià)指標(biāo)。因此，解析吐絮率性狀的遺傳基礎(chǔ)，對(duì)于快速培育機(jī)采棉專用品種方面具有重要意義。然而，目前有關(guān)陸地棉吐絮率性狀遺傳基礎(chǔ)解析的研究很少，使得吐絮率的研究落后于其他早熟相關(guān)性狀?！緮M解決的關(guān)鍵問題】本研究采用RTM-GWAS方法對(duì)315份陸地棉材料的吐絮率性狀展開研究，挖掘控制吐絮率性狀的優(yōu)異等位變異位點(diǎn)，并進(jìn)行候選基因預(yù)測，以期解析陸地棉吐絮率的遺傳基礎(chǔ)，為陸地棉早熟和適宜機(jī)采相關(guān)性狀的分子育種奠定理論基礎(chǔ)。

1 材料與方法

1.1 試驗(yàn)材料與田間試驗(yàn)

甘肅農(nóng)業(yè)大學(xué)棉花分子育種實(shí)驗(yàn)室前期構(gòu)建了由315份陸地棉種質(zhì)材料構(gòu)成的一個(gè)自然群體，包含290份來自中國不同時(shí)期、不同地域選育的品種（系），25份從國外引進(jìn)的品種（Foreign Region，F(xiàn)R）[36]。290份國內(nèi)材料按照地域來源進(jìn)行分類，125份來源與黃河流域棉區(qū)（Yellow River Region，YRR）、48份來源于長江流域棉區(qū)（Yangtze River Region，YZRR）、97份來源于西北內(nèi)陸棉區(qū)（Northwest Inland Region，NIR）和20份來源于北部特早熟棉區(qū)（Northern Super Early-maturing Region，NSER）（電子附表 1）。2014—2015年，將該陸地棉自然群體的315份材料種植于中國農(nóng)業(yè)科學(xué)院棉花研究所安陽老所部試驗(yàn)田（河南安陽，2個(gè)環(huán)境分別命名為AY-14和AY-15）；2014年種植于新疆農(nóng)墾科學(xué)院棉花研究所試驗(yàn)田（新疆石河子，該環(huán)境命名為SHZ-14）。每個(gè)環(huán)境設(shè)置3個(gè)重復(fù)，采取隨機(jī)區(qū)組試驗(yàn)設(shè)計(jì)方案，田間管理參照當(dāng)?shù)爻Ｒ?guī)方法。河南安陽試驗(yàn)點(diǎn)：每個(gè)材料單行種植，行長5.00 m，行距0.80 m，株距0.20 m，每個(gè)材料約23個(gè)單株；新疆石河子試驗(yàn)點(diǎn)：每個(gè)材料分2行種植，行長2.00 m，行距0.45 m，株距0.10 m，每個(gè)材料約40個(gè)單株。AY-2014和AY-2015均于當(dāng)年的4月19日播種，SHZ-2014于4月25日播種。

1.2 表型性狀鑒定與數(shù)據(jù)分析

在各試驗(yàn)點(diǎn)的材料正常吐絮后，對(duì)每個(gè)材料的每個(gè)重復(fù)，挑選長勢整齊、均勻一致的10株單株進(jìn)行系紅線標(biāo)記。3個(gè)環(huán)境（AY-14、AY-15和SHZ-14）的吐絮率調(diào)查時(shí)間分別為9月18日、9月10日、9月28日，調(diào)查上述被標(biāo)記10個(gè)單株上的總棉鈴數(shù)及吐絮棉鈴數(shù)（單株中直徑大于2.0 cm的棉鈴被定為有效棉鈴）。計(jì)算每個(gè)種質(zhì)材料的吐絮率，其計(jì)算公式為：吐絮率=（吐絮棉鈴數(shù)/總棉鈴數(shù)）×100%。然后利用SPSS 26.0軟件對(duì)獲得吐絮率的表型數(shù)據(jù)進(jìn)行統(tǒng)計(jì)描述和方差分析，并計(jì)算其廣義遺傳力。

1.3 全基因組關(guān)聯(lián)分析（GWAS）

陸地棉自然群體的 SNPLDB標(biāo)記構(gòu)建來源于甘肅農(nóng)業(yè)大學(xué)棉花分子育種實(shí)驗(yàn)室前期的研究基礎(chǔ)[36]。利用SLAF-seq（specific-locus simplified fragment sequencing）方法對(duì)自然群體進(jìn)行SNP位點(diǎn)基因分型，參考陸地棉TM-1 v2.1基因組[37]對(duì)測序所獲得的序列進(jìn)行比對(duì)，然后根據(jù)最大缺失率＞0.10，最小等位基因頻率（minor allele frequency，MAF）＜0.05的標(biāo)準(zhǔn)過濾挖掘，獲得13 391個(gè)高質(zhì)量SNP位點(diǎn)。將分布在同一連鎖不平衡（linkage disequilibrium，LD）區(qū)塊內(nèi)的多個(gè)緊密連鎖SNP標(biāo)記，記為同一個(gè)SNP連鎖不平衡區(qū)段（SNP linkage disequilibrium block，SNPLDB），獲得9 244個(gè)SNPLDB標(biāo)記[36]。然后利用這些SNPLDB標(biāo)記和3個(gè)種植環(huán)境下吐絮率的表型值，對(duì) AY-14、AY-15、SHZ-14 3個(gè)單環(huán)境表型和多環(huán)境聯(lián)合分析獲得的表型分別進(jìn)行RTM-GWAS分析。利用SAS軟件PROC GLM 程序中方差分析的線性模型，按照郝曉帥等[38]提出的方法，估算獲得多環(huán)境聯(lián)合分析的表型值。然后按照HE等[29]提出的RTM-GWAS程序進(jìn)行關(guān)聯(lián)分析，該程序主要分為2個(gè)分析階段：（1）基于單位點(diǎn)簡單線性模型對(duì)SNPLDB標(biāo)記進(jìn)行初步篩選；（2）使用多位點(diǎn)逐步回歸模型進(jìn)行篩選，同時(shí)估計(jì)每個(gè)顯著位點(diǎn)的等位基因效應(yīng)值。根據(jù)等位基因效應(yīng)值，利用 R 3.6.0軟件建立目標(biāo)性狀的QTL-allele矩陣。在RTMGWAS 2個(gè)階段的分析程序中，顯著水平（P）均設(shè)置為0.05。另外，利用R 3.6.0軟件，以PCs和親緣關(guān)系系數(shù)K為協(xié)變量進(jìn)行MLM-GWAS分析，鑒定吐絮率性狀的顯著SNP位點(diǎn)，利用R 3.6.0軟件繪制出關(guān)聯(lián)位點(diǎn)的曼哈頓圖，以-lg(P)≥4.43的SNP為顯著位點(diǎn)。

1.4 穩(wěn)定關(guān)聯(lián)位點(diǎn)優(yōu)異等位變異的鑒定

能在多環(huán)境聯(lián)合分析、單環(huán)境2個(gè)分析模式中同時(shí)檢測到的SNPLDB位點(diǎn)，被認(rèn)為與吐絮率顯著關(guān)聯(lián)的穩(wěn)定關(guān)聯(lián)位點(diǎn)。利用 SPSS 26.0軟件計(jì)算穩(wěn)定SNPLDB位點(diǎn)各種等位變異類型表型的平均值，并進(jìn)行雙尾T檢驗(yàn)，確定優(yōu)異等位變異類型，使用Origin 2018繪制相對(duì)應(yīng)的統(tǒng)計(jì)圖。

1.5 優(yōu)異等位變異類型頻率的比較分析

在315份陸地棉種質(zhì)材料中，分別選取擁有最高和最低吐絮率的品種系各50個(gè)，分析比較穩(wěn)定關(guān)聯(lián)位點(diǎn)優(yōu)異等位變異類型的頻率；然后比較YRR、YZRR、NIR和NSER 4個(gè)不同地理來源種質(zhì)材料間優(yōu)異等位變異類型的頻率差異。

1.6 候選基因的挖掘

根據(jù)SNPLDB位點(diǎn)在陸地棉TM-1 v2.1參考基因組[38]上的物理位置，結(jié)合 LD的衰減距離，將穩(wěn)定顯著關(guān)聯(lián)SNPLDB位點(diǎn)的兩端各擴(kuò)展500 kb，形成約1 Mb基因組區(qū)間。利用棉花功能基因組學(xué)數(shù)據(jù)庫（https://cottonfgd.org/），列出分布在上述顯著關(guān)聯(lián)的穩(wěn)定SNPLDB位點(diǎn)兩側(cè)1 Mb區(qū)間內(nèi)的全部基因。然后分別利用南京農(nóng)業(yè)大學(xué)（Nanjing Agricultural University，NAU）和中國農(nóng)業(yè)科學(xué)院棉花研究所（Cotton Research Institute of Chinese Academy of Agricultural Sciences，CRI）2組不同組織的RNA-Seq數(shù)據(jù)，計(jì)算它們的FPKM（fragments per kilobase of transcript per million fragments mapped）值，篩選出在棉鈴胚珠和纖維組織中高表達(dá)的基因作為候選基因，并對(duì)候選基因進(jìn)行功能注釋。候選基因表達(dá)方式的熱圖（heatmap）和維恩圖（venn diagram）分析使用OmicShare在線工具（http://www.omicshare.com/tools）繪制完成。最后利用棉花功能基因組學(xué)數(shù)據(jù)庫對(duì)候選基因進(jìn)行GO富集分析和KEGG代謝途徑分析。

2 結(jié)果

2.1 陸地棉吐絮率的表型變異特征

315個(gè)陸地棉品種（系）吐絮率在3個(gè)環(huán)境中的變異范圍為 39.27%—91.32%，平均值為 69.14%，變異系數(shù)為14.25%。AY-14環(huán)境下的吐絮率分布范圍相對(duì)較小，變幅為46.73%—90.00%，平均值為80.62%，變異系數(shù)為9.86%。AY-15環(huán)境下的吐絮率分布范圍相對(duì)較大，變幅為 10.42%—100.00%，平均值為55.99%，變異系數(shù)為34.60%。SHZ-14環(huán)境下的吐絮率分布范圍也較小，變幅為37.78%—98.05%，平均值為70.80%，變異系數(shù)為14.42%（表1）。3個(gè)環(huán)境下吐絮率的方差分析結(jié)果顯示，基因型、環(huán)境以及基因型×環(huán)境互作均呈現(xiàn)極顯著差異（P＜0.0001），說明吐絮率的基因型與環(huán)境之間存在顯著的互作效應(yīng)（表2）。利用3個(gè)環(huán)境計(jì)算的吐絮率廣義遺傳力相對(duì)較大，為67.03%（表2）。以上結(jié)果表明，陸地棉吐絮率性狀主要受遺傳因素影響，也受環(huán)境因素影響。

表1 陸地棉吐絮率的次數(shù)分布和描述統(tǒng)計(jì)Table 1 Frequency distribution and descriptive statistics of boll opening rate in upland cotton

表2 陸地棉吐絮率的方差分析Table 2 Variance analysis of boll opening rate in upland cotton

2.2 全基因組關(guān)聯(lián)分析（GWAS）

甘肅農(nóng)業(yè)大學(xué)棉花分子育種實(shí)驗(yàn)室前期已利用SNPLDB標(biāo)記對(duì)陸地棉自然群體進(jìn)行了聚類分析和主成分分析（principal component analysis，PCA），把315個(gè)品種（系）劃分為2個(gè)類群；LD分析發(fā)現(xiàn)該群體在全基因組水平上的平均LD值約為500 kb[36]。通過RTM-GWAS方法的多環(huán)境分析模型，在9 244個(gè)SNPLDB中檢測到52個(gè)SNPLDB位點(diǎn)與吐絮率顯著關(guān)聯(lián)，這些顯著關(guān)聯(lián)位點(diǎn)分布在陸地棉全部26條染色體上，每條染色體上分布1—5個(gè)位點(diǎn)（表3）。第2、4、6、7、14、18、20、24和26染色體上最少，僅檢測到1個(gè)SNPLDB標(biāo)記；第15染色體上最多，檢測到5個(gè)SNPLDB標(biāo)記。在52個(gè)顯著關(guān)聯(lián)的SNPLDB位點(diǎn)中，10個(gè)關(guān)聯(lián)位點(diǎn)的-lg(P)值大于10.00，其中，LDB_16_37952328上具有最高的-lg(P)值（46.14）。進(jìn)一步分析發(fā)現(xiàn)，在多環(huán)境聯(lián)合分析中檢測到的 52個(gè)顯著關(guān)聯(lián)位點(diǎn)中，8個(gè)位點(diǎn)同時(shí)也至少在1個(gè)單環(huán)境中也被檢測到，例如，位點(diǎn)LDB_16_37952328在2個(gè)單環(huán)境（AY-15和SHZ-14）中被檢測到；LDB_5_96395565、LDB_16_49503485、LDB_10_6908012、LDB_15_9697358和LDB_4_81118668在AY-15單個(gè)環(huán)境中被檢測到。根據(jù)顯著關(guān)聯(lián)位點(diǎn)的-lg(P)≥10，以及在單環(huán)境和多環(huán)境聯(lián)合分析中同時(shí)出現(xiàn)的穩(wěn)定性，最終篩選出6個(gè)SNPLDB可能是與吐絮率穩(wěn)定顯著關(guān)聯(lián)的主效位點(diǎn)，它們分別為 LDB_16_37952328、LDB_5_96395565、LDB_16_49503485、LDB_10_6908012、LDB_15_9697358和LDB_4_81118668（圖1-A）。這6個(gè)位點(diǎn)均為單個(gè)SNP構(gòu)成SNPLDB位點(diǎn)，共可解釋7.04%的表型變異（表3）。本研究利用MLM-GWAS方法僅發(fā)現(xiàn)1個(gè)位點(diǎn)SNP_16_37952328與吐絮率顯著關(guān)聯(lián)（圖1-B）。

圖1 利用RTM-GWAS和MLM-GWAS方法對(duì)陸地棉吐絮率的遺傳解析Fig.1 Genetic analysis of boll opening rate in upland cotton by MLM-GWAS and RTM-GWAS procedure

表3 與陸地棉吐絮率顯著關(guān)聯(lián)的SNPLDB位點(diǎn)Table 3 SNPLDB loci significantly associated with boll opening rate in upland cotton

續(xù)表3 Continued table 3

52個(gè)顯著關(guān)聯(lián)的SNPLDB位點(diǎn)中，38個(gè)為單個(gè)SNP的等位變異位點(diǎn)，14個(gè)為含有多個(gè)SNP的復(fù)等位變異位點(diǎn)。這些顯著關(guān)聯(lián)SNPLDB位點(diǎn)等位變異類型的分布范圍為3—10個(gè)，共179個(gè)，其中90個(gè)為增效等位變異/單倍型，89個(gè)為減效等位變異/單倍型。增效等位變異類型的效應(yīng)值分布在0.014—19.43，減效等位變異類型的效應(yīng)值分布在-21.49—-0.039。顯著關(guān)聯(lián)位點(diǎn)中，復(fù)合位點(diǎn)LDB_19_52309050_52309284的等位變異效應(yīng)值變異范圍較大（-21.49—19.43），其他等位變異/單倍型的效應(yīng)值變異范圍較為集中，約96%的等位變異/單倍型效應(yīng)值分布在-6—5（圖1-C）。為了更加明晰地展示陸地棉種質(zhì)群體的遺傳基礎(chǔ)，利用52個(gè)顯著SNPLDB位點(diǎn)的等位變異效應(yīng)值，構(gòu)建了一個(gè) 52×315（位點(diǎn)×品種系）的QTL等位變異矩陣，該矩陣能夠直觀地反映了315個(gè)陸地棉種質(zhì)系中控制吐絮率的遺傳變異構(gòu)成（圖1-D）。

2.3 穩(wěn)定顯著關(guān)聯(lián)位點(diǎn)優(yōu)異等位變異的鑒定

為了鑒定上述6個(gè)穩(wěn)定顯著關(guān)聯(lián)位點(diǎn)的優(yōu)異等位變異類型，對(duì)其每種等位變異所對(duì)應(yīng)的吐絮率進(jìn)行了顯著性比較。關(guān)聯(lián)位點(diǎn)LDB_16_37952328存在有CC（n=217）、CT（n=23）和TT（n=67）3種等位變異類型，攜帶有TT和CT等位變異種質(zhì)的吐絮率（分別為80.08%和77.17%）極顯著高于攜帶有CC的種質(zhì)系（67.07%）（圖2-A）。LDB_5_96395565存在AA（n=212）、AG（n=17）和 GG（n=53）3種等位變異類型，AA型的吐絮率（71.87%）極顯著高于 GG型（67.26%）（圖2-B）。LDB_16_49503485存在有CC（n=36）、CT（n=30）和TT（n=249）3種等位變異類型，TT型種質(zhì)的吐絮率平均值（71.12%）顯著高于含有CC型（67.04%）和CT型（67.25%）（圖2-C）。LDB_4_81118668存在 CC（n=76）、CT（n=37）和TT（n=182)3種等位變異類型，CT型（72.84%）和 TT型（71.82%）種質(zhì)系的吐絮率極顯著高于 CC型（67.83%）（圖 2-D）。其余 2個(gè)穩(wěn)定關(guān)聯(lián)位點(diǎn)LDB_10_6908012和LDB_15_9697358的3種等位變異類型所對(duì)應(yīng)種質(zhì)的吐絮率無顯著性差異（圖2-E和圖 2-F）。按照不同等位變異類型對(duì)應(yīng)種質(zhì)的吐絮率有顯著性差異，吐絮率高的等位變異類型為優(yōu)異等位變異原則，4個(gè)穩(wěn)定顯著關(guān)聯(lián)位點(diǎn)及其優(yōu)異等位變異被鑒定，其優(yōu)異等位變異分別為 LDB_16_37952328（TT）、LDB_5_96395565（AA）、LDB_16_49503485（TT）、和LDB_4_81118668（TT）。

圖2 6個(gè)穩(wěn)定顯著關(guān)聯(lián)SNPLDB位點(diǎn)不同等位變異的吐絮率比較Fig.2 Comparison of boll opening rate with different allelic variation at six stable and significant association SNPLDB loci

2.4 優(yōu)異等位變異分布頻率的比較

為了研究優(yōu)異等位變異分布頻率與吐絮率的關(guān)系，分別選取2組高、低吐絮率的極端材料各50份。高吐絮率種質(zhì)的表型分布在81.88%—91.32%，低吐絮率種質(zhì)的表型分布在 39.27%—60.59%。比較結(jié)果發(fā)現(xiàn)，4個(gè)穩(wěn)定關(guān)聯(lián)SNPLDB位點(diǎn)在高吐絮率品種（系）中的優(yōu)異等位變異頻率均高于低吐絮率品種（系）（圖3-A）。

為了比較不同地理來源的材料間吐絮率及其優(yōu)異等位變異頻率的差異性，選取來自YRR、YZRR、NSER和NIR4個(gè)不同地理來源區(qū)域的125、48、20和95個(gè)品種（系）。不同地理來源品種（系）的吐絮率平均值大小順序?yàn)镹SER（75.91%）＞YRR（73.95%）＞NIR（69.96%）＞YZRR（64.72%）。方差分析結(jié)果表明，NSER品種（系）的吐絮率極顯著高于 NIR和YZRR，且YRR和NIR品種（系）的吐絮率極顯著高于YZRR（P＜0.01，圖3-B）。此外，對(duì)4個(gè)穩(wěn)定關(guān)聯(lián)位點(diǎn)的優(yōu)異等位基因頻率進(jìn)行比較，發(fā)現(xiàn)4個(gè)位點(diǎn)的優(yōu)異等位基因頻率在中國 4個(gè)不同地理來源區(qū)域的品種（系）間存在差異。進(jìn)一步分析發(fā)現(xiàn)，位點(diǎn)LDB_16_49503485在YRR、YZRR和NIR有較高的優(yōu)異等位基因頻率，均在 70%以上（圖3-C）。另外3個(gè)穩(wěn)定關(guān)聯(lián)位點(diǎn)（LDB_16_37952328、LDB_5_96395565和LDB_4_81118668），在YRR和NSER品種（系）中的優(yōu)異等位變異位點(diǎn)頻率要高于來源于YZRR和NIR的品種（系）（圖3-D—F），優(yōu)異等位變異頻率的變化與不同地理來源品種（系）的吐絮率表型性狀變化趨勢一致。

圖3 4個(gè)穩(wěn)定關(guān)聯(lián)位點(diǎn)優(yōu)異等位變異的頻率分布Fig.3 Frequency distribution of superior alleles of four stable association loci

2.5 候選基因的預(yù)測

根據(jù)陸地棉TM-1v2.1參考基因組[38]，在4個(gè)穩(wěn)定關(guān)聯(lián)位點(diǎn)兩側(cè)1 Mb擴(kuò)展區(qū)域中存在178個(gè)基因。利用NAU測序獲得的RNA-seq計(jì)算基因FPKM值，發(fā)現(xiàn)其中61個(gè)基因，至少在陸地棉10個(gè)組織或器官中一個(gè)中有較高表達(dá)量。依據(jù)基因在不同組織中的表達(dá)情況，可把61個(gè)基因分為4組（用1—4表示）。第1組（17個(gè)）主要在第5、10和20 DPA（days post anthesis，開花后生長天數(shù)）的纖維中高表達(dá)；第2組（9個(gè)）主要在花瓣和雄蕊中高表達(dá)；第3組（17個(gè)）主要在胚珠的8個(gè)發(fā)育階段中高表達(dá)；第4組（18個(gè)）主要在根、莖、葉組織和花器官中較高表達(dá)（圖4-A）。利用另一組CRI測序獲得的RNA-seq數(shù)據(jù)，發(fā)現(xiàn)59個(gè)基因至少在 12個(gè)組織或器官的一個(gè)中有較高表達(dá)量，將59個(gè)基因分為4組（用I—Ⅳ表示）。第Ⅰ組（16個(gè)）在花瓣、花托、花萼、苞片、花藥和花絲中高表達(dá)；第Ⅱ組（13個(gè)）在第10、15、20和25 DPA的纖維中高表達(dá)；第Ⅲ組（9個(gè)）在根、莖和葉中高表達(dá)；第Ⅳ組（21個(gè)）在第10、15、20和25 DPA的胚珠中高表達(dá)（圖4-B）。

通過NAU和CRI的2組RNA-seq數(shù)據(jù)，各發(fā)現(xiàn)34個(gè)基因在胚珠和纖維中高表達(dá)。它們中的23個(gè)基因在兩組數(shù)據(jù)中的表達(dá)模式一致，推測它們可能是與陸地棉吐絮率相關(guān)的候選基因（圖 4-C）。將預(yù)測的23個(gè)候選基因可分為3類基因，a類基因（9個(gè)）主要在0—5 DPA的胚珠中表達(dá)，b類基因（4個(gè)）主要在10—35 DPA的胚珠中表達(dá)，c類基因（10個(gè)）主要在5—25 DPA的纖維中表達(dá)。進(jìn)一步分析發(fā)現(xiàn)，20個(gè)候選基因具有生物學(xué)功能注釋，其中過氧化物酶體(S)-2-羥基氧化酶基因Gh_D03G1603和鈣結(jié)合線粒體載體蛋白基因 Gh_D03G1058可能是促進(jìn)棉鈴成熟或早衰的調(diào)控基因（表4）。

表4 與陸地棉吐絮率相關(guān)的候選基因Table 4 Candidate genes related to BOR in upland cotton

圖4 與陸地棉吐絮率相關(guān)候選基因的預(yù)測Fig.4 Prediction of candidate genes related to the BOR of upland cotton

對(duì)上述預(yù)測的23個(gè)候選基因進(jìn)行GO富集分析和KEGG代謝途徑分析。GO結(jié)果顯示，在20個(gè)基因中總共發(fā)現(xiàn)了 55個(gè)調(diào)控途徑，其中 Gh_D03G1608、Gh_D03G1610、Gh_D03G1616、Gh_D03G1629和Gh_D03G1067與蛋白質(zhì)氨基酸糖蛋白結(jié)合相關(guān)（GO：0005515）。KEGG結(jié)果顯示，在11個(gè)基因中發(fā)現(xiàn)了32個(gè)調(diào)控途徑，其中 Gh_A05G3660、Gh_D03G1063和Gh_D03G1067有共同的核糖體代謝途徑（ko01100）（電子附表2）。

3 討論

3.1 吐絮率在機(jī)采棉生產(chǎn)中的重要性

2020年全國棉花播種面積已達(dá)316.99萬hm2，總產(chǎn)量為591.00萬t；西北內(nèi)陸區(qū)棉花種植面積為251.85萬hm2，總產(chǎn)量為519.10萬t，其產(chǎn)量占全國棉花總產(chǎn)量的 87.83%（http://www.stats.gov.cn/tjsj/zxfb/202012/t20201218_1810113.html）。棉花已成為西北內(nèi)陸最主要的經(jīng)濟(jì)作物之一。然而該地區(qū)春季溫度回升慢、秋季下降快、初霜來臨早，一旦發(fā)生霜凍，會(huì)造成棉花減產(chǎn)，因此，培育早熟棉品種顯得尤為重要[39]。早熟性也是適宜機(jī)采棉品種最為關(guān)鍵的性狀之一。機(jī)械采棉是目前降低中國棉花生產(chǎn)成本，實(shí)現(xiàn)棉花生產(chǎn)可持續(xù)發(fā)展，提升棉花國際競爭力的重要措施之一。在西北內(nèi)陸棉區(qū)，為了提高棉花機(jī)收效率而不降低其產(chǎn)量和纖維品質(zhì)，所以在9月中下旬棉田整體吐絮率達(dá)到50%以上時(shí)，才噴施脫葉催熟劑。如果棉田吐絮率低于50%時(shí)，噴施脫葉催熟劑會(huì)顯著降低中上部棉鈴的單鈴重及纖維品質(zhì)，從而造成棉花減產(chǎn)降質(zhì)，因此，吐絮率已成為西北機(jī)采棉育種中不可忽視的一個(gè)重要性狀。

棉花從播種出苗到群體一半的第一個(gè)棉鈴?fù)滦跛璧奶鞌?shù)稱為全生育期。該性狀全生育期因品種、氣候及栽培條件的不同而異，通常陸地棉品種全生育期約為100—140 d。然而全生育期并不包含吐絮期（第一棉鈴?fù)滦醯阶詈笠粋€(gè)棉鈴?fù)滦醯臅r(shí)間段）的生育階段。以前的研究中，把霜前花率作為評(píng)價(jià)吐絮期吐絮進(jìn)程的一項(xiàng)重要指標(biāo)[1-2]。霜前花率是指早霜以前及早霜后5 d內(nèi)所采摘的棉花產(chǎn)量占總產(chǎn)量的百分率。然而該指標(biāo)在目前的機(jī)采棉生產(chǎn)模式下已失去了實(shí)際意義，無法指導(dǎo)當(dāng)前的機(jī)采棉生產(chǎn)。如果一個(gè)棉花品種的吐絮期過長，造成吐絮不夠集中，會(huì)影響機(jī)械化采收進(jìn)程，故吐絮集中已成為早熟的一個(gè)關(guān)鍵性狀。吐絮率不僅僅是一個(gè)反映棉花早熟性的重要指標(biāo)之一，而且是反映吐絮集中程度的重要指標(biāo)。吐絮越早越集中，吐絮率越高，反之則越低。然而目前有關(guān)棉花吐絮率的研究主要集中在脫葉劑方面[40]，缺乏從遺傳基礎(chǔ)方面的解析和研究，對(duì)控制棉花吐絮率QTL定位和候選基因挖掘等方面的研究較少。本研究中全面解析陸地棉吐絮率性狀遺傳基礎(chǔ)，挖掘調(diào)控棉花吐絮率QTL及候選基因，對(duì)于今后開展早熟機(jī)采棉分子育種具有重要意義。

3.2 與吐絮率穩(wěn)定關(guān)聯(lián)位點(diǎn)及其優(yōu)異等位變異的發(fā)掘

準(zhǔn)確的表型鑒定是利用GWAS挖掘農(nóng)作物目標(biāo)性狀優(yōu)異等位變異和候選基因的關(guān)鍵。本研究發(fā)現(xiàn)315個(gè)陸地棉品種（系）吐絮率存在廣泛的表型差異，且3個(gè)環(huán)境間表型值存在極顯著差異。其原因是吐絮率大小易受溫度、光照、有效積溫等外界環(huán)境因素的影響。在不同環(huán)境種植陸地棉種質(zhì)材料，由于鑒定吐絮率的外界環(huán)境因素在不同地點(diǎn)、不同年份之間會(huì)存在不小的差異，因此獲得表型性狀也會(huì)存在較大的差異。如AY-2014和AY-2015，同一地點(diǎn)不同年份之間其表型數(shù)據(jù)也存在顯著差異。故分別利用單個(gè)環(huán)境的表型值和多環(huán)境聯(lián)合估算的表型值進(jìn)行了 RTM-GWAS分析，鑒定能夠與目標(biāo)性狀穩(wěn)定關(guān)聯(lián)SNPLDB位點(diǎn)。本研究中利用多環(huán)境聯(lián)合估算的表型值進(jìn)行了 RTM-GWAS分析，共獲得了52個(gè)SNPLDB位點(diǎn)與吐絮率顯著關(guān)聯(lián)。將該52個(gè)與吐絮率顯著關(guān)聯(lián)的SNPLDB位點(diǎn)與前人發(fā)現(xiàn)的其他早熟性狀相關(guān) QTL進(jìn)行比較，發(fā)現(xiàn)位點(diǎn)LDB_16_37952328、LDB_17_51465593、LDB_5_23314512和 LDB_21_22036171與前人研究報(bào)道的QTL相重疊[10,25,41-42]。其中分布在D03染色體上位點(diǎn)LDB_16_37952328在多篇文章中被發(fā)現(xiàn)[23-24,28]，這些結(jié)果證實(shí)了本研究結(jié)果的可靠性，也證明該位點(diǎn)是控制陸地棉早熟相關(guān)性狀的一個(gè)關(guān)鍵位點(diǎn)。它不僅與開花期和全生育期顯著相關(guān)，而且與吐絮率顯著相關(guān)，值得在后期的研究中重點(diǎn)關(guān)注。此外，其余48個(gè)位點(diǎn)應(yīng)該是新發(fā)現(xiàn)與吐絮率相關(guān)的位點(diǎn)?？傊?，目前與其他早熟相關(guān)性狀 QTL或關(guān)聯(lián)位點(diǎn)的研究相對(duì)較多[3,5,22,28]，但有關(guān)與棉花吐絮率穩(wěn)定關(guān)聯(lián)的位點(diǎn)或候選基因的研究很少。本研究利用兩組前人報(bào)道的 RNA-seq數(shù)據(jù)，發(fā)現(xiàn) 23個(gè)基因在胚珠和纖維中高表達(dá)，推測它們可能是與陸地棉吐絮率相關(guān)的候選基因。更為重要的是，推測Gh_D03G1603和Gh_D03G1058可能通過調(diào)節(jié)棉鈴成熟或早衰，促進(jìn)棉花集中早絮。Gh_D03G1603是過氧化物酶體(S)-2-羥基氧化酶基因（GLO1），水稻中該類基因通過改變H2O2含量來調(diào)節(jié)生長素含量，進(jìn)而促進(jìn)水稻早熟[43]。由此我們推測 Gh_D03G1603通過調(diào)控棉花子房中的生長素含量提高，促進(jìn)子房及其周圍組織的膨大，加速了棉鈴發(fā)育；當(dāng)生長素達(dá)到一定含量時(shí)還會(huì)誘導(dǎo)乙烯的形成，進(jìn)而促進(jìn)棉花集中早絮成熟。Gh_D03G1058是鈣結(jié)合線粒體載體蛋白基因（slc25a24），擬南芥中調(diào)節(jié)激活Ca2+運(yùn)輸通道介導(dǎo)脫落酸合成，導(dǎo)致氣孔關(guān)閉，促進(jìn)擬南芥早衰[44]。推測Gh_D03G1058可能調(diào)節(jié)棉鈴中Ca2+運(yùn)輸，激活大量脫落酸合成，促進(jìn)棉鈴的發(fā)育和成熟。

根據(jù)顯著關(guān)聯(lián)位點(diǎn)的-lg(P)值的大小以及在多環(huán)境聯(lián)合分析和單環(huán)境分析中顯著位點(diǎn)同時(shí)出現(xiàn)的穩(wěn)定性，在52個(gè)SNPLDB位點(diǎn)中最終篩選出6個(gè)位點(diǎn)可能與吐絮率穩(wěn)定關(guān)聯(lián)的主效位點(diǎn)。通過不同等位變異類型對(duì)應(yīng)種質(zhì)的吐絮率的顯著差異比較分析，最終確定4個(gè)位點(diǎn)是穩(wěn)定顯著關(guān)聯(lián)的主效位點(diǎn)，并鑒定出該4個(gè)主效位點(diǎn)的優(yōu)異等位變異類型。進(jìn)一步分析發(fā)現(xiàn)，4個(gè)穩(wěn)定關(guān)聯(lián)主效位點(diǎn)的優(yōu)異等位變異在高吐絮率種質(zhì)系中出現(xiàn)頻率高于低吐絮率種質(zhì)系中出現(xiàn)的頻率。例如，特早熟陸地棉優(yōu)良品種新石 K18中，同時(shí)包含了該4個(gè)與吐絮率關(guān)聯(lián)SNPLDB位點(diǎn)的優(yōu)異等位變異類型。此外，在4個(gè)不同地理來源的品種系間，YRR和NSER的平均吐絮率較高，YZRR和NIR的平均吐絮率較低，同時(shí)位點(diǎn)LDB_16_37952328、LDB_5_96395565、和LDB_4_81118668在4個(gè)不同地域的優(yōu)異等位變異頻率與之相對(duì)應(yīng)，進(jìn)一步確定了這些優(yōu)異等位變異的可靠性。因此，可以根據(jù)這些優(yōu)異等位變異篩選出優(yōu)良種質(zhì)材料，為分子標(biāo)記輔助選擇提供優(yōu)異親本資源。

3.3 RTM-GWAS方法的優(yōu)越性

RTM-GWAS方法能夠較全面地解析陸地棉吐絮率性狀QTL及其復(fù)等位變異，分析結(jié)果用于個(gè)別基因挖掘的同時(shí)，還可用于群體遺傳結(jié)構(gòu)分析以及作物育種的優(yōu)化組合設(shè)計(jì)等方面[30]。RTM-GWAS方法基于包含不同單倍型/等位基因的 SNP連鎖不平衡區(qū)塊，這種方法比常規(guī)GWAS更適合QTL鑒定；該方法可以提供具有比較總體遺傳信息的QTL-等位基因矩陣，包括顯著的單核苷酸多態(tài)性位點(diǎn)及其等位基因效應(yīng)、群體等位基因組成和每個(gè)位點(diǎn)等位基因頻率的差異，這是以往傳統(tǒng)的單位點(diǎn)GWAS模型無法做到的[45-46]。本研究利用RTM-GWAS方法對(duì)陸地棉吐絮率性狀進(jìn)行分析，檢測到52個(gè)顯著關(guān)聯(lián)SNPLDB位點(diǎn)，其中4個(gè)位點(diǎn)在多環(huán)境聯(lián)合分析和單環(huán)境分析中同時(shí)被檢測到，被認(rèn)定為穩(wěn)定關(guān)聯(lián)的主效SNPLDB位點(diǎn)；而利用MLM-GWAS進(jìn)行吐絮率的研究時(shí)，僅檢測到1個(gè)顯著 SNP位點(diǎn)。此外我們?cè)谇捌诶?MLM-GWAS進(jìn)行其他6個(gè)早熟相關(guān)性狀的研究中，檢測到8個(gè)SNP位點(diǎn)與5個(gè)早熟性狀之間存在13個(gè)顯著關(guān)聯(lián)，發(fā)現(xiàn)1個(gè)穩(wěn)定SNP位點(diǎn)（rsDt3：13562854）同時(shí)多個(gè)早熟相關(guān)性狀（-lg(P)≥6.21）之間存在顯著關(guān)聯(lián)[24]。這些結(jié)果都充分證明了RTM-GWAS方法在全基因組關(guān)聯(lián)位點(diǎn)檢測方面具有優(yōu)勢。

4 結(jié)論

在陸地棉自然群體中檢測到52個(gè)與吐絮率性狀顯著相關(guān)的SNPLDB位點(diǎn)。其中4個(gè)被鑒定為與吐絮率穩(wěn)定關(guān)聯(lián)的主效SNPLDB位點(diǎn)。在該4個(gè)主效位點(diǎn)的鄰近區(qū)域共注釋了178個(gè)基因，預(yù)測發(fā)現(xiàn)了其中的 23個(gè)基因可能是調(diào)控陸地棉吐絮率的候選基因。