張繼峰，劉華東，王敬國，劉化龍，孫健，楊洛淼，賈琰，吳文申，鄭洪亮，鄒德堂

（東北農(nóng)業(yè)大學(xué)寒地糧食作物種質(zhì)創(chuàng)新與生理生態(tài)教育部重點(diǎn)實(shí)驗(yàn)室，哈爾濱 150030）

0 引言

【研究意義】水稻是世界上主要糧食作物之一，隨著全球人口的增長和耕地面積的日益減少，水稻每年需要達(dá)到 1%的增產(chǎn)才能滿足人類對(duì)食物的需求。在經(jīng)歷矮化育種和雜種優(yōu)勢(shì)2次突破性飛躍后，水稻增產(chǎn)達(dá)到了瓶頸。長期栽培實(shí)踐證明，雜種優(yōu)勢(shì)與理想株型相互利用是實(shí)現(xiàn)水稻進(jìn)一步增產(chǎn)的重要途徑。分蘗數(shù)是水稻理想株型的重要組成部分之一，也是水稻產(chǎn)量形成的重要性狀。分蘗數(shù)決定有效穗數(shù)，適量的分蘗數(shù)會(huì)有助于高產(chǎn)，分蘗數(shù)過多會(huì)導(dǎo)致營養(yǎng)物質(zhì)的浪費(fèi)，分蘗數(shù)過少會(huì)導(dǎo)致過少的有效穗數(shù)，因此，適量的分蘗數(shù)是決定水稻高產(chǎn)的重要條件?！厩叭搜芯窟M(jìn)展】水稻分蘗由分蘗節(jié)上各葉腋下處的分生組織發(fā)育而來，水稻分蘗形成包括2個(gè)部分，分別是分蘗芽的形成和向外生長過程。植物遺傳學(xué)家現(xiàn)在已經(jīng)克隆了MOC1、MOC2、MOC3、LAX1、LAX2和DTE1等與分蘗芽形成相關(guān)的基因。MOC1編碼一個(gè)GRAS家族的核蛋白，主要在腋芽和葉腋原基表達(dá)，在營養(yǎng)生長和生殖生長階段控制葉腋分生組織形成，moc1突變體表現(xiàn)為少蘗，同時(shí)花序軸和小穗減少[1-2]。MOC2/FBP1編碼果糖-1,6二磷酸酶，該基因缺失表現(xiàn)為單孽、矮化、葉色變淡、每穗實(shí)粒數(shù)減少等，moc2突變體與moc1突變體不同，可以生長分蘗芽，但是分蘗芽的生長受到阻礙，不能正常發(fā)育成分蘗，果糖-1,6二磷酸酶是一個(gè)胞質(zhì)FBP酶，該酶失活會(huì)導(dǎo)致植株蔗糖供應(yīng)不足，從而導(dǎo)致moc2缺失突變體分蘗芽受到抑制[3]。MOC3/OsWUS/TAB1編碼一個(gè)WOX蛋白家族成員，是一個(gè)轉(zhuǎn)錄抑制因子，主要在地上部愈傷組織中表達(dá)[4]，該基因正向調(diào)控腋芽發(fā)育，可能是通過促進(jìn)水稻同源異型基因OSH1的表達(dá)，進(jìn)而在莖端分生組織中發(fā)揮重要作用[5]，moc3-1突變體和moc2突變體一樣，在分蘗芽發(fā)育過程中遭受破壞[3-4]，Tab1/wus突變體表現(xiàn)為不能形成腋芽，導(dǎo)致沒有分蘗[5]。LAX1編碼一個(gè)植物特有的bHLH轉(zhuǎn)錄因子，該轉(zhuǎn)錄子含有一個(gè)堿性螺旋-環(huán)-螺旋（bHLH）結(jié)構(gòu)域，該基因在地上部分頂端分生組織與新形成的分生組織之間的邊界區(qū)域表達(dá)，LAX1蛋白在葉腋分生組織中積累，被 mRMA和隨后的蛋白轉(zhuǎn)運(yùn)進(jìn)行嚴(yán)格調(diào)控，該步驟對(duì)葉腋分生組織建立發(fā)揮重要作用[6-7]。LAX2/Gnp4編碼一個(gè)核蛋白，和LAX1蛋白互作共同參與調(diào)控水稻腋生分生組織形成，進(jìn)而調(diào)控水稻分蘗形態(tài)[8]，Lax2與lax1突變體在營養(yǎng)階段葉腋分生組織數(shù)目減少且大多數(shù)側(cè)生小穗腋分生組織減少[8]。DTE1編碼一個(gè)硼酸通道蛋白，在根的外層皮以及成熟葉片節(jié)區(qū)域表達(dá)，在硼饑餓處理下，dtel突變體表現(xiàn)出分蘗數(shù)增加、生長緩慢和育性下降等生長缺陷[9]。研究表明，大部分基因均通過雙親分離群體（RIL、NIL和BIL等）進(jìn)行遺傳位點(diǎn)挖掘，但鑒定的分蘗數(shù)遺傳位點(diǎn)受到雙親遺傳差異或分子標(biāo)記密度的限制，所能挖掘鑒定的位點(diǎn)相對(duì)有限。全基因組關(guān)聯(lián)分析以連鎖不平衡（linkage disequilibrium，LD）為基礎(chǔ)，通過群體基因型數(shù)據(jù)與表型數(shù)據(jù)檢測(cè)有效關(guān)聯(lián)位點(diǎn)，可以實(shí)現(xiàn)對(duì)同一復(fù)雜性狀多數(shù)基因座以及一個(gè)基因座上多個(gè)復(fù)等位基因的檢測(cè)[10-11]。同時(shí)，全基因組關(guān)聯(lián)分析基于自然群體，不需要構(gòu)建專門的作圖群體，并且自然群體經(jīng)過天然雜交多代，重組更加充分，定位精度更高。在人類和動(dòng)物研究的基因挖掘中，關(guān)聯(lián)分析是重要的手段，并且已在植物基因定位中得到廣泛應(yīng)用。2005年，Science雜志首次報(bào)道了人類年齡相關(guān)的黃斑性病癥全基因組關(guān)聯(lián)分析的研究[12]。受到人類遺傳疾病研究的啟示，植物遺傳學(xué)家也將關(guān)聯(lián)分析應(yīng)用于作物遺傳育種中，并取得了重要研究進(jìn)展，為挖掘水稻有效等位基因提供了新方法。隨著基因分型與測(cè)序技術(shù)的日益進(jìn)步，目前，對(duì)于水稻的粒型、株高、抽穗期和劍葉形態(tài)等很多農(nóng)藝性狀開展了GWAS分析，同時(shí)也對(duì)病蟲害抗性等生物脅迫和鹽堿耐受性等非生物脅迫進(jìn)行了全基因組關(guān)聯(lián)分析，檢測(cè)出許多重要的顯著性SNP位點(diǎn)及候選基因，解析了許多重要性狀的遺傳基礎(chǔ)[13-22]?！颈狙芯壳腥朦c(diǎn)】盡管已有很多關(guān)于水稻分蘗數(shù)QTL/基因的研究報(bào)道，但大多是通過突變體基因定位或反向遺傳學(xué)的方法進(jìn)行研究，所挖掘到的QTL/基因有限?！緮M解決的關(guān)鍵問題】本研究以來自不同國家和地區(qū)的295份溫帶粳稻品種為試驗(yàn)材料，結(jié)合高通量重測(cè)序獲得的788 396個(gè)高質(zhì)量多態(tài)性SNP，對(duì)粳稻分蘗數(shù)進(jìn)行全基因組關(guān)聯(lián)分析，鑒定影響粳稻分蘗數(shù)穩(wěn)定存在的主效QTL，并對(duì)QTL區(qū)間內(nèi)所有基因進(jìn)行單倍型分析，根據(jù)基因注釋篩選候選基因，為分子輔助育種提高產(chǎn)量提供理論基礎(chǔ)。

1 材料與方法

1.1 試驗(yàn)材料以及性狀考察

試驗(yàn)材料由國內(nèi)外廣泛收集的295份粳稻品種構(gòu)成，其中，國內(nèi)材料主要來自黑龍江省、吉林省和遼寧省等，國外材料主要來自日本、韓國和俄羅斯等（電子附表1）。所有材料于2018—2019年種植在東北農(nóng)業(yè)大學(xué)阿城水稻試驗(yàn)基地，每份材料種植8行，每行20株，行、株間距為30 cm×13.3 cm，單株插秧。田間管理依照常規(guī)大田生產(chǎn)的基本方法。為了減少邊際效應(yīng)對(duì)表型的影響，于水稻分蘗盛期，在每個(gè)品種的第4行，選取表型相似的連續(xù)5株，考察每株水稻分蘗數(shù)，并計(jì)算平均值。

1.2 全基因組重測(cè)序

使用植物基因組DNA提取試劑盒提取295份粳稻品種的DNA，將檢測(cè)合格的樣品DNA采用Illumina HiSeq XTen進(jìn)行高通量測(cè)序。利用GATK軟件工具包進(jìn)行SNP檢測(cè)，使用Plink軟件從初步獲得的群體SNP中篩選出最小等位基因頻率（minimum allele frequency，MAF）大于 5%，缺失率（missing rate）小于20%的788 396個(gè)SNP用于后續(xù)分析[23]。

1.3 群體結(jié)構(gòu)與親緣關(guān)系分析

利用 ADMIXTURE軟件計(jì)算每個(gè)品種的遺傳成分，將K值設(shè)置為1—10，運(yùn)行結(jié)束后，提取每個(gè)K值的交叉驗(yàn)證（cross validation，CV）值，最小CV值所對(duì)應(yīng)的K值為最優(yōu)亞群分組數(shù)，所用群體被劃分為3個(gè)亞群[23]。另外使用GCTA軟件[24]進(jìn)行主成分分析（principal component analysis，PCA），利用R語言繪制三維 PCA圖。利用 MEGA 7[25-26]的鄰接法（neighbor-joining）構(gòu)建系統(tǒng)進(jìn)化樹（NJ樹）。使用TASSEl 5.0軟件[27]對(duì)群體材料間的親緣關(guān)系進(jìn)行評(píng)估[23]。

1.4 粳稻分蘗數(shù)全基因組關(guān)聯(lián)分析

利用TASSEl 5.0軟件的混合線性模型（mixed linear model，MLM），結(jié)合群體SNP基因型數(shù)據(jù)，進(jìn)行粳稻分蘗數(shù)的全基因組關(guān)聯(lián)分析，對(duì)GEC軟件計(jì)算出的有效獨(dú)立SNP數(shù)目進(jìn)行閾值的確定，最終將P＜5.46×10-6作為閾值，判定SNP標(biāo)記與目標(biāo)性狀關(guān)聯(lián)的顯著性。利用R語言的“CMplot”包繪制曼哈頓圖和 Q-Q圖。依據(jù)每條染色體 LD衰減結(jié)果[23]，以GWAS檢測(cè)到的峰值SNP所在染色體上下游LD衰減距離區(qū)域界定為該QTL的染色體區(qū)間。

1.5 候選基因單倍型分析

對(duì)2年環(huán)境下共同檢測(cè)到的QTL內(nèi)所有基因進(jìn)行單倍型分析。具體操作過程如下：從3KRGP的Rice SNP-Seek Database網(wǎng)站[28]中提取位于QTL區(qū)間內(nèi)所有基因外顯子區(qū)非同義突變SNP。利用之前研究中篩選得到的788 396個(gè)高質(zhì)量SNP在位于該QTL區(qū)間內(nèi)的所有SNP與3KRGP中提取的所有非同義突變SNP進(jìn)行比對(duì)，確定其是否為非同義突變。若基因存在非同義突變，則先利用DnaSP軟件對(duì)該QTL區(qū)間內(nèi)所有存在非同義突變的基因進(jìn)行單倍型分析，對(duì)不同單倍型（≥15份材料）的表型值進(jìn)行多重比較，篩選不同單倍型表型值之間存在顯著性差異的基因。對(duì)于不同單倍型表型值之間差異不顯著的基因，再用啟動(dòng)子區(qū)（ATG前1.5 kb）的SNP進(jìn)行單倍型分析。結(jié)合基因注釋預(yù)測(cè)可能的候選基因。

2 結(jié)果

2.1 295份粳稻材料分蘗數(shù)表現(xiàn)

2018和2019年295份粳稻品種的分蘗數(shù)均具有較大的表型分布，分別為7.00—49.00個(gè)和6.00—33.80個(gè)（表1和圖1），且2年總趨勢(shì)基本一致，分蘗數(shù)平均值分別為18.35和15.75個(gè)。

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

利用TASSLE 5.0軟件的MLM模型分析，2018年在水稻第 8、9和 10染色體共檢測(cè)到 3個(gè) QTL（qTiller8、qTillerq9和qTiller10），2019年在水稻第9染色體檢測(cè)到1個(gè)QTL（qTiller9），其中，qTiller9在2年中均被檢測(cè)到（表2）。依據(jù)第9染色體LD衰減距離（52 kb）[23]，選取峰值 SNP所在染色體的上下游LD衰減距離區(qū)域，界定為該QTL的區(qū)間（表2），結(jié)果顯示，qTiller9被定位于水稻第9染色體15.01—15.12 Mb，與已定位影響水稻分蘗數(shù)QTL（qNOT9-1）[29]位于相同區(qū)間。

圖1 2年環(huán)境條件下分蘗數(shù)頻次分布直方圖Fig.1 Histogram of tillering number frequency distribution under 2-year environmental conditions

表1 2年環(huán)境下295份粳稻材料的分蘗數(shù)表型值Table 1 The phenotypic values of tillering number among 295 japonica rice varieties in two years

表2 粳稻分蘗數(shù)顯著關(guān)聯(lián)位點(diǎn)Table 2 Significant loci associated with tillering number in japonica rice

圖2 分蘗數(shù)GWAS結(jié)果的曼哈頓圖和QQ散點(diǎn)圖Fig.2 Manhattan plots and quantile-quantile (Q-Q) plots of genome-wide association studies for the tillering number

2.3 候選基因分析

2018—2019年共檢測(cè)到 3個(gè)影響粳稻分蘗數(shù)的QTL，其中，2年中能夠穩(wěn)定表達(dá)的QTL——qTiller9位于已定位影響水稻分蘗數(shù)QTL（qNOT9-1）區(qū)間內(nèi)。因此，對(duì)qTiller9區(qū)間內(nèi)全部15個(gè)基因的外顯子區(qū)非同義突變SNP和啟動(dòng)子區(qū)SNP進(jìn)行單倍型分析。結(jié)果顯示，共有 6個(gè)基因（LOC_Os09g25090、LOC_Os09g25100、LOC_Os09g25150、LOC_Os09g25190、LOC_Os09g25200和LOC_Os09g25220）的不同單倍型分蘗數(shù)之間均檢測(cè)到極顯著差異水平（圖 3）。LOC_Os09g25090被啟動(dòng)子區(qū)SNP分為2種單倍型，hap2（TAA）分蘗數(shù)極其顯著大于 hap1（AGG）；LOC_Os09g25100被非同義突變SNP分為2種單倍型，hap2（GAGA）分蘗數(shù)極其顯著大于hap1（AGCC）；LOC_Os09g25150被非同義突變SNP分為2種單倍型，hap2（ATG）分蘗數(shù)極其顯著大于 hap1（GCC）；LOC_Os09g25190被啟動(dòng)子區(qū)SNP分為2種單倍型，hap2（GCATCGCATCGACGCCGA）分蘗數(shù)極其顯著大于 hap1（ATGCTGATGAAGTCATCC）；LOC_Os09g25200被非同義突變SNP分為2種單倍型，hap2（TAG）分蘗數(shù)極其顯著大于 hap1（AGA）；LOC_Os09g25220被非同義突變SNP分為2種單倍型，hap1（GG）分蘗數(shù)極其顯著大于 hap2（AA）（表 3和圖3）。其中，LOC_Os09g25090和LOC_Os09g25100編碼鈣調(diào)蛋白依賴性蛋白激酶，它是脫落酸（ABA）表達(dá)所必需Ca2+的傳感器。因此，推測(cè)LOC_Os09g25090和LOC_Os09g25100可能是影響水稻分蘗數(shù)的候選基因。

表3 候選基因單倍型分組及每種單倍型的SNP組成Table 3 Candidate gene haplotype group and the composition of each haplotype SNP

圖3 候選基因不同單倍型之間分蘗數(shù)的箱線圖Fig.3 Boxplots for the tillering number based on the haplotypes(hap) for candidate gene

3 討論

水稻分蘗數(shù)是由多基因控制的數(shù)量性狀，挖掘水稻分蘗數(shù)相關(guān)QTL和基因?qū)λ驹霎a(chǎn)具有重要意義。連鎖分析已經(jīng)被證明可以有效地挖掘水稻分蘗數(shù)QTL，但是由于親本數(shù)目較少，后代群體重組機(jī)會(huì)有限，因此，當(dāng)?shù)任换蛟陔p親中無差異時(shí)，QTL無法通過連鎖分析檢測(cè)出來。與連鎖分析相比，全基因組關(guān)聯(lián)分析以數(shù)百個(gè)品種組成的自然群體為研究材料，自然群體已經(jīng)過天然雜交多代，重組更加充分，定位精度更高，并且能夠同時(shí)檢測(cè)同一基因座上的多個(gè)等位基因。本研究利用全基因組關(guān)聯(lián)分析方法，于2018和2019年在水稻第8、9和10染色體共檢測(cè)到3個(gè) QTL，分別為qTiller8、qTiller9和qTiller10。HEMAMALINI等[29]利用IR64和Azucena雜交得到的56個(gè)隨機(jī)DH系在水稻第9染色體11.80—15.55 Mb區(qū)間內(nèi)定位到一個(gè)影響水稻分蘗數(shù)的 QTL（qNOT9-1），與本研究qTiller9位于相同區(qū)間，這驗(yàn)證了本研究結(jié)果可靠性的同時(shí)進(jìn)一步縮小了 QTL區(qū)間，為后續(xù)候選基因分析提供了堅(jiān)實(shí)的理論依據(jù)；qTiller8和qTiller10在前人研究中未見報(bào)道，且2個(gè)QTL僅在2018年檢測(cè)到，需要通過遺傳群體進(jìn)一步驗(yàn)證。本研究進(jìn)一步通過QTL所在的物理位置與前人已克隆基因進(jìn)行比較，未發(fā)現(xiàn)QTL區(qū)間內(nèi)有已知基因的報(bào)道。

表4 候選基因的基因注釋Table 4 Candidate gene of gene annotation

以往的全基因組關(guān)聯(lián)分析挖掘相關(guān)性狀候選基因，大多依據(jù)顯著性峰值SNP位點(diǎn)上下游固定物理距離確定候選基因區(qū)域，然后通過基因功能注釋進(jìn)行候選基因的預(yù)測(cè)，然而水稻每條染色體的LD衰減距離各不相同，僅僅依靠固定距離確定候選基因區(qū)域，降低了該區(qū)域的精準(zhǔn)性。本研究的檢測(cè)到的QTL區(qū)域是依據(jù)每條染色體LD衰減距離及顯著性峰值SNP位點(diǎn)確定的，這大大增加候選基因區(qū)域的精準(zhǔn)性，同時(shí)本研究對(duì)共定位QTL區(qū)域內(nèi)所有基因進(jìn)行單倍型分析，再通過基因注釋篩選單倍型表型值之間存在顯著性差異的基因，這大大增加了候選基因篩選的效率及可靠性。

本研究通過對(duì)qTiller9區(qū)間內(nèi)全部15個(gè)基因進(jìn)行單倍型分析，發(fā)現(xiàn)有6個(gè)基因不同單倍型的分蘗數(shù)存在極顯著差異。進(jìn)一步通過基因功能注釋篩選候選基因，發(fā)現(xiàn)6個(gè)候選基因分別屬于鈣調(diào)蛋白依賴性蛋白激酶、肉桂酰輔酶 A還原酶和蛋白結(jié)合蛋白。其中LOC_Os09g2509和LOC_Os09g25100均預(yù)測(cè)編碼鈣調(diào)蛋白依賴性蛋白激酶，其結(jié)構(gòu)包括CaM結(jié)合域、絲氨酸/蘇氨酸激酶域和3個(gè)EF手基序，是多個(gè)Ca2+傳感器中的一員，可以檢測(cè)到Ca2+的瞬時(shí)濃度[30-31]，同時(shí)，Ca2+對(duì)脫落酸（ABA）的表達(dá)是必需的[37-39]。ZHU等[32]研究發(fā)現(xiàn)，9順式-環(huán)氧類胡蘿卜素雙加氧酶OsNCED2是控制ABA合成代謝的關(guān)鍵酶，過表達(dá)OsNCED2的水稻轉(zhuǎn)基因品系中，其分蘗數(shù)減少[33]；另外，擬南芥基因MAX3和MAX4均編碼類胡蘿卜素裂解雙加氧酶，是 ABA生物合成中的關(guān)鍵酶[34]，擬南芥突變體max4存在一種側(cè)芽抑制信號(hào)，該信號(hào)由根部向上運(yùn)輸，并與生長素互作進(jìn)而抑制植株側(cè)芽發(fā)生，這種信號(hào)的前體是類胡蘿卜素，而植物中 ABA的前體同樣為類胡蘿卜素[35]；STEVEN等[36]進(jìn)一步研究發(fā)現(xiàn)，擬南芥中的max3和max4突變體是由MAX3和MAX4損傷引起的，導(dǎo)致2種突變體側(cè)向分支急劇增加。以上研究表明，ABA在植物側(cè)芽發(fā)生中發(fā)揮重要作用。本研究在水稻表達(dá)數(shù)據(jù)庫（Rice Expression Profile Database）中查詢LOC_Os09g25090和LOC_Os09g25100在田間整個(gè)生長生育時(shí)期各組織的時(shí)空表達(dá)量，發(fā)現(xiàn)2個(gè)候選基因在水稻莖部均有較高強(qiáng)度的表達(dá)。因此，推測(cè)2個(gè)候選基因可能在水稻分蘗處表達(dá)，通過調(diào)節(jié)Ca2+濃度介導(dǎo)ABA合成代謝進(jìn)而影響水稻分蘗數(shù)。

4 結(jié)論

2018和2019年共檢測(cè)到3個(gè)與粳稻分蘗數(shù)相關(guān)QTL（qTiller8、qTiller9和qTiller10），其中，qTiller9在 2年中均穩(wěn)定表達(dá)，并與已定位影響水稻分蘗數(shù)QTL（qNOT9-1）重合。LOC_Os09g25090和LOC_Os09g25100可能為影響qTiller9的候選基因。