林 淼 張秋芝 史利玉 劉 蓓 王紅武 潘金豹

（1北京農(nóng)學(xué)院，102206，北京；2中玉金標(biāo)記（北京）生物技術(shù)股份有限公司，102206，北京）

青貯玉米是草食家畜特別是奶牛最重要的高能量飼料來源。優(yōu)質(zhì)青貯玉米具有較高的淀粉含量、低纖維含量和高消化率等特征。玉米秸稈占全株玉米生物量的一半，因此秸稈消化率的高低是影響青貯玉米消化率關(guān)鍵。粗飼料體外干物質(zhì)消化率（in vitro dry matter digestibility，IVDMD）測(cè)定指標(biāo)有30h和48h兩種，48h體外干物質(zhì)消化率（48h IVDMD）是青貯飼料在草食家畜的消化系統(tǒng)中消化并吸收48h的最大量與進(jìn)食量的比值。IVDMD體現(xiàn)青貯玉米秸稈利用和轉(zhuǎn)化的效率，主要受遺傳因素影響，48h IVDMD是青貯玉米品質(zhì)育種首要指標(biāo)[1-4]。

目前，國(guó)內(nèi)外利用不同遺傳背景的群體，采用連鎖分析的方法，在1～10染色體上均發(fā)現(xiàn)了控制玉米秸稈IVDMD的QTL[5-8]。全基因組關(guān)聯(lián)分析（genome-wide association study，GWAS）具有高分辨率且省時(shí)等特點(diǎn)，是檢測(cè)復(fù)雜性狀變異強(qiáng)大的遺傳方法，特別是玉米作為異交作物，重組較多，導(dǎo)致較短的連鎖不平衡（linkage disequilibrium，LD），且表型和遺傳多樣性豐富，分辨率較高[9-11]，但不同群體LD存在差異，如農(nóng)家種1kb、廣泛變異群體1.5kb、優(yōu)良自交系100kb，都利于關(guān)聯(lián)定位的精度[10]。目前關(guān)于青貯玉米秸稈消化率性狀的全基因組關(guān)聯(lián)分析的研究極少。Wang等[9]以559 285個(gè)高質(zhì)量SNPs對(duì)368份玉米自交系采用GWAS方法鑒定了82個(gè)秸稈IVDMD顯著相關(guān)的位點(diǎn)，并發(fā)現(xiàn)了直接參與細(xì)胞壁生物合成的基因ZmC3H2。但Wang等[9]利用關(guān)聯(lián)分析對(duì)秸稈IVDMD進(jìn)行研究時(shí)，群體中沒有青貯玉米品種的親本自交系，從而不利于挖掘更多且直接的優(yōu)良的青貯玉米消化率性狀遺傳位點(diǎn)。

本研究以341份玉米自交系為材料，構(gòu)建關(guān)聯(lián)作圖群體，基于全基因組測(cè)序的高質(zhì)量SNPs進(jìn)行秸稈48h IVDMD的全基因組關(guān)聯(lián)分析，挖掘優(yōu)良等位變異，篩選候選基因，為優(yōu)質(zhì)青貯玉米選育提供參考。

1 材料與方法

1.1 供試材料及田間設(shè)計(jì)

供試材料為194份美國(guó)玉米自交系和147份國(guó)內(nèi)玉米自交系（包括14份國(guó)內(nèi)常用玉米自交系和133份課題組自選玉米自交系）構(gòu)成的341份關(guān)聯(lián)群體，其中6份為青貯玉米品種的親本自交系。

2018年5月初將供試材料分別種植于遼寧省沈陽(yáng)市東亞國(guó)際種業(yè)公司和內(nèi)蒙古自治區(qū)通遼市農(nóng)業(yè)科學(xué)研究院科學(xué)研究基地。行長(zhǎng)5.0m，行距0.6m，株距0.25m，雙行區(qū)，自然散粉，田間管理同一般玉米自交系試驗(yàn)田。

1.2 品質(zhì)測(cè)定方法

調(diào)查供試材料生育時(shí)期（主要包括抽雄期和吐絲期）。在玉米籽粒乳線位置1/2時(shí)，選取長(zhǎng)勢(shì)一致的10株玉米，從地上部20cm處刈割，去掉果穗，收獲秸稈。用風(fēng)送式動(dòng)力鍘草機(jī)粉碎后混合均勻，平行取兩份1 000g裝入布袋，立刻放入105℃烘箱中殺青2h，散熱散水氣后調(diào)至65℃烘干至恒重，用錘片式粉碎機(jī)磨碎，過40目篩，裝入牛皮紙袋，于65℃烘箱烘至恒重，取出樣品冷卻至室溫，待用。

試驗(yàn)采用FOSS近紅外光譜儀（NIRS DS2500）[12]對(duì)供試樣品48h IVDMD（%）進(jìn)行測(cè)定，該儀器采用前分光單色儀技術(shù)，探測(cè)器是硅（400～1 100nm）和硫化鉛（1 100～2 500nm），光譜分辨率是0.5nm-1。將每份樣品混勻，取適量于直徑為50mm的石英柱形樣品杯中，波長(zhǎng)為400～2 500nm近紅外光以漫反射形式從樣杯底部掃描，每份樣品掃描2次，每次重復(fù)重新取樣品，結(jié)果取均值（所測(cè)定含量是干基）。

1.3 表型統(tǒng)計(jì)方法

采用SPSS 20.0進(jìn)行描述性統(tǒng)計(jì)分析和相關(guān)性分析，利用R 3.5.1軟件對(duì)表型頻率分布作圖，利用Microsoft Excel 2010進(jìn)行雙因素方差分析估算廣義遺傳力（HB2），HB2=VG/VP×100%=VG（/VG+VE）×100%，VP、VG和VE分別是表型變異方差、遺傳變異方差和環(huán)境變異方差。

1.4 基因型測(cè)定及分析

采用磁珠法[13]提取341份供試材料的基因組DNA，利用Qubit熒光光度計(jì)、NanoDrop 2000和1%瓊脂糖凝膠電泳分別對(duì)樣品DNA進(jìn)行濃度、純度和完整性檢測(cè)。使用Illumina PE150測(cè)序儀對(duì)合格樣品DNA進(jìn)行全基因組重測(cè)序，平均測(cè)序深度是7×，共得到112 626 032個(gè)SNPs標(biāo)記。采用BWA軟件[14]將測(cè)序后短序列和玉米B73參考基因組組裝V4比對(duì)，平均覆蓋度是83%，均勻覆蓋參考基因組。從112 626 032個(gè)SNPs標(biāo)記中去掉缺失率大于30%、雜合率大于10%、最小等位基因頻率小于0.05的標(biāo)記，剩余6 276 612個(gè)高質(zhì)量SNPs用于全基因組關(guān)聯(lián)分析。

1.5 群體結(jié)構(gòu)和親緣關(guān)系

采用fastStructure 1.0軟件[15]，以6 276 612個(gè)MAF＞0.05且均勻分布染色體的SNPs，對(duì)341份供試材料進(jìn)行貝葉斯聚類，計(jì)算不同層次的亞群數(shù)（K=2～10），再根據(jù)最大似然數(shù)LnP（D），確定本群體理論上由5大亞群組成。

1.6 秸稈48h IVDMD全基因組關(guān)聯(lián)分析

在考慮群體結(jié)構(gòu)和親緣關(guān)系的情況下，運(yùn)用MVP軟件的FarmCPU模型[16]將秸稈48h IVDMD表型和本群體的基因型進(jìn)行全基因組關(guān)聯(lián)。Bonferroni校正法即0.05/n（n是本研究的標(biāo)記數(shù)目）確定閾值（P＜1.0×10-8），由于此種方法稍嚴(yán)格，故通過文獻(xiàn)采用閾值（P＜1.0×10-6）[17]。

1.7 秸稈48h IVDMD定位結(jié)果比較分析

借助前人用不同群體構(gòu)建遺傳連鎖圖譜進(jìn)行的QTL定位結(jié)果，若本試驗(yàn)關(guān)聯(lián)定位結(jié)果落入前人QTL定位區(qū)間內(nèi)或相差200kb以內(nèi)說明本試驗(yàn)關(guān)聯(lián)定位效果良好。

1.8 秸稈48h IVDMD候選基因預(yù)測(cè)

采用Haploview 4.1 jar[18]分析本群體連鎖不平衡（LD）衰減距離，本群體LD衰減到0.1時(shí)的距離約是2.7kb，在此情況下，為了防止盲目基因或區(qū)段擴(kuò)展，故基于玉米B73基因組組裝V4序列信息和GWAS檢測(cè)到的所有顯著SNPs位點(diǎn)（P＜1.0×10-6），選取基因內(nèi)部及其上下游1kb潛在調(diào)控區(qū)域作為候選基因區(qū)段[19-20]。

1.9 群體變異位點(diǎn)分析

利用GATK軟件[21]采用Haplotype Caller模式進(jìn)行變異位點(diǎn)檢測(cè)，借助ANNOVAR軟件[17]對(duì)基于候選基因的SNPs進(jìn)行注釋。

2 結(jié)果與分析

2.1 秸稈48h IVDMD表型分析

沈陽(yáng)和通遼兩個(gè)環(huán)境下秸稈48h IVDMD頻率分布圖和基本描述統(tǒng)計(jì)見圖1和表1，結(jié)果說明遺傳變異范圍非常廣泛，均值分別為80.96%和67.94%，變異系數(shù)分別為4.85%和9.50%，說明同一環(huán)境的不同自交系間的秸稈48h IVDMD存在差異，不同環(huán)境下群體的整體表現(xiàn)也有差異，但整體符合正態(tài)分布。Person相關(guān)性結(jié)果表明，兩個(gè)環(huán)境下秸稈48h IVDMD顯著相關(guān)。秸稈48h IVDMD的廣義遺傳力為71.98%。這些結(jié)果都說明秸稈48h IVDMD性狀主要受到遺傳效應(yīng)影響。

2.2 秸稈48h IVDMD全基因組關(guān)聯(lián)分析

圖1 秸稈48h IVDMD頻率分布Fig.1 Frequency distribution of 48h IVDMD

表1 兩個(gè)環(huán)境秸稈48h IVDMD統(tǒng)計(jì)分析Table 1 Statistics analysis of 48h IVDMD in two environments

利用FarmCPU模型對(duì)341份玉米自交系的48h IVDMD進(jìn)行全基因組關(guān)聯(lián)分析。繪制曼哈頓圖（Manhattan）（圖2），超過水平線（兩條水平線代表全基因組顯著水平閾值6.0和8.0）的SNPs即為與本試驗(yàn)?zāi)繕?biāo)性狀顯著相關(guān)的位點(diǎn)。兩個(gè)環(huán)境中FarmCPU模型的QQ圖表示群體結(jié)構(gòu)較為良好地被控制（圖3）。在沈陽(yáng)和通遼的兩個(gè)環(huán)境中共檢測(cè)到153個(gè)顯著SNPs位點(diǎn)（P＜1.0×10-6），分布在玉米第1～10染色體，其中，沈陽(yáng)環(huán)境中，在顯著水平P＜1.0×10-8上檢測(cè)到4個(gè)顯著SNPs位點(diǎn)，分布于第2、3、6、8染色體上（圖4）。在第1～3、6～8、10染色體上，沈陽(yáng)和通遼兩個(gè)環(huán)境中都檢測(cè)出了顯著SNPs位點(diǎn)（P＜1.0×10-6），但未發(fā)現(xiàn)相同的顯著SNPs位點(diǎn)，位于10號(hào)染色體的chr10_131199866標(biāo)記（沈陽(yáng)）和 chr10_130000984、chr10_130001456、chr10_130001612、chr10_130001742、chr10_130001760 5個(gè)連鎖標(biāo)記（通遼）最近，距離是1.20Mb。

圖2 不同環(huán)境（沈陽(yáng)和通遼）秸稈48h IVDMD全基因組關(guān)聯(lián)曼哈頓圖Fig.2 Manhattan plot of genome-wide association under the environments of Shenyang and Tongliao

圖3 不同環(huán)境秸稈48h IVDMD全基因組關(guān)聯(lián)QQ圖Fig.3 Quantile-quantile plot of genome-wide association under the environments of Shenyang and Tongliao

圖4 顯著SNPs在染色體上分布Fig.4 Significant SNPs were distributed on chromosomes

采用混合線性模型（MLM）對(duì)341份玉米自交系的48h IVDMD進(jìn)行全基因組關(guān)聯(lián)，共發(fā)現(xiàn)分布在第2、4～9染色體的10個(gè)SNPs顯著位點(diǎn)（P＜1.0×10-6），結(jié)果仍沒有發(fā)現(xiàn)兩個(gè)環(huán)境相同的SNPs顯著位點(diǎn)（P＜1.0×10-6）（圖 5）。利用 MLM不僅沒有關(guān)聯(lián)到兩個(gè)環(huán)境中和兩個(gè)性狀顯著相關(guān)的相同的位點(diǎn)，而且兩個(gè)環(huán)境中在共同染色體上也沒有檢測(cè)到顯著位點(diǎn)。從兩個(gè)模型的曼哈頓圖和QQ圖可以看出，F(xiàn)armCPU模型關(guān)聯(lián)的顯著位點(diǎn)優(yōu)于MLM的，說明FarmCPU模型關(guān)聯(lián)效率要高于MLM，候選基因易從FarmCPU模型關(guān)聯(lián)的顯著位點(diǎn)中預(yù)測(cè)。

圖5 兩個(gè)環(huán)境秸稈48h IVDMD全基因組關(guān)聯(lián)曼哈頓圖及其QQ圖（MLM）Fig.5 Manhattan plots and its QQ plots of genome-wide association in two environments (MLM)

2.3 候選基因分析

基于FarmCPU模型關(guān)聯(lián)的153個(gè)NDF的SNPs顯著位點(diǎn)（P＜1.0×10-6）預(yù)測(cè)候選基因，結(jié)合玉米B73參考基因組組裝V4注釋信息，共預(yù)測(cè)到38個(gè)候選基因。在38個(gè)候選基因中，有3個(gè)候選基因編碼未知功能蛋白，有1個(gè)候選基因編碼MYB轉(zhuǎn)錄因子，其余大多數(shù)候選基因參與生長(zhǎng)發(fā)育和防御反應(yīng)的生物過程，少數(shù)候選基因涉及信號(hào)轉(zhuǎn)導(dǎo)途徑、細(xì)胞代謝途徑和RNA相關(guān)途徑。最重要的是，本次關(guān)聯(lián)分析發(fā)現(xiàn)了Zm00001d009009基因和Zm00001d028874基因可能與玉米秸稈纖維素和半纖維素合成有關(guān)[19-20]。

群體變異位點(diǎn)分析表明，64個(gè)SNPs顯著位點(diǎn)（P＜1.0×10-6）位于38個(gè)候選基因中的34.38%內(nèi)含子區(qū)、28.13%外顯子區(qū)、23.44% 5'/3'非翻譯區(qū)、12.50%基因上下游區(qū)、及1.56%非編碼RNA區(qū)（圖6）。其中，在28.13%（18個(gè)SNPs）外顯子變異中，有13個(gè)非同義SNPs位于9個(gè)候選基因中，這些基因功能涉及RNA/DNA結(jié)合蛋白、開花調(diào)控、ABA信號(hào)轉(zhuǎn)導(dǎo)負(fù)調(diào)控因子、葉和花形態(tài)建成、耐熱、十八烷類防御通路[21-25]。

圖6 變異基因的位置Fig.6 The location of the variant gene

3 討論

3.1 GWAS定位結(jié)果分析

Li等[5]以鄭58和HD568組成的220份RIL群體為材料，用6個(gè)消化性狀定位了47個(gè)QTLs，其中位于2、7～10號(hào)染色體的9個(gè)QTLs控制IVDMD，解釋單個(gè)QTL 5.3%的表型變異，其中定位在9.03染色體上QTL和Chr9_25266410標(biāo)記最小僅距88kb。Wang等[7]以CE03005高油玉米自交系和B73組成的211份材料對(duì)玉米消化性狀進(jìn)行QTL研究，其中，在4～6和10號(hào)染色體上發(fā)現(xiàn)有6個(gè)位點(diǎn)與IVDMD顯著相關(guān)，解釋單個(gè)位點(diǎn)10.6%的表型變異。本試驗(yàn)和48h IVDMD顯著相關(guān)的標(biāo)記 chr4_59965077、chr4_63943611、chr4_69161939、chr4_74067130、chr4_74067534和 chr4_74243482均定位在了bin4.05染色體區(qū)段內(nèi)，本試驗(yàn)和48h IVDMD顯著相關(guān)的標(biāo)記chr6_31376191、chr6_8473014、chr6_8473017、chr6_54601996和 chr6_54602105均定位在了bin6.01染色體區(qū)段內(nèi)，本試驗(yàn)和48h IVDMD顯著相關(guān)的標(biāo)記chr10_46816004、ch10_63652102、chr10_65363792、chr10_15312554、chr10_17307934和chr10_48100033均定位在了bin10.03染色體區(qū)段內(nèi)，與本研究結(jié)果基本一致。王琪[6]用B73和By804作親本構(gòu)建RIL群體，在3號(hào)染色體鑒定了消化率的遺傳位點(diǎn)，本試驗(yàn)在這個(gè)區(qū)域也發(fā)現(xiàn)了相關(guān)QTL。

本試驗(yàn)并非所有顯著SNP位點(diǎn)都在已知的QTL定位區(qū)間，可能原因在于：QTL作圖受親本遺傳背景或檢測(cè)微效QTL功效低的原因；雙親群體中的功能位點(diǎn)可能在關(guān)聯(lián)分析群體中是稀有等位變異，因而導(dǎo)致在關(guān)聯(lián)分析群體中不能檢測(cè)；關(guān)聯(lián)分析群體存在雙親群體中可能不存在的遺傳變異；關(guān)于IVDMD的QTL研究非常少見。

3.2 有潛力的候選基因分析

鑒定玉米秸稈消化率的相關(guān)基因能更好地理解細(xì)胞壁合成的分子遺傳機(jī)制。在本試驗(yàn)38個(gè)候選基因中，基因Zm00001d009009編碼產(chǎn)物類糖基轉(zhuǎn)移酶KOBITO1是高度保守的、具有植物特異性的新型質(zhì)膜結(jié)合蛋白，該蛋白是擬南芥細(xì)胞擴(kuò)張過程中合成纖維素所必需的蛋白。在所有光和黑暗環(huán)境中，KOB1基因突變導(dǎo)致擬南芥下胚軸變短，與野生型植株相比，該突變體中纖維素的相對(duì)含量下降了33%，證明該蛋白通過調(diào)節(jié)纖維素合成，促進(jìn)植株下胚軸生長(zhǎng)[26-27]?；騔m00001d028874編碼產(chǎn)物是木葡聚糖6-木糖基轉(zhuǎn)移酶2，木葡聚糖是禾本科植物初生細(xì)胞壁半纖維素合成的主要組分，木糖基轉(zhuǎn)移酶2（XXT2）將核苷酸二磷酸糖中的木糖殘基轉(zhuǎn)移到木聚糖主鏈中的葡聚糖O-6位置，研究表明，XXT2突變使木聚糖含量下降近30%，另外，木葡聚糖6-木糖基轉(zhuǎn)移酶突變導(dǎo)致水稻根毛發(fā)育異常；有研究表明，編碼木葡聚糖6-木糖基轉(zhuǎn)移酶基因OsXXT1在維持水稻細(xì)胞壁結(jié)構(gòu)和抗拉強(qiáng)度方面很重要[28-30]。

秸稈細(xì)胞壁相關(guān)性狀不僅能調(diào)節(jié)玉米生長(zhǎng)，而且在作物育種上，與抗逆性表現(xiàn)出很強(qiáng)的相關(guān)性[31]，可以增強(qiáng)對(duì)病蟲的免疫機(jī)能，提高對(duì)不利環(huán)境的適應(yīng)能力。在本試驗(yàn)38個(gè)候選基因中，基因Zm00001d031449含有1個(gè)非同義SNP，將333位和678位的異亮氨酸突變成甲硫氨酸，編碼葉綠體脂氧合酶Ⅱ。蟲害的番茄葉片中TomLoxD基因上調(diào)，編碼1個(gè)葉綠體LOX，該LOX可能是十八烷類防御信號(hào)通路的組成部分，該通路在一系列酶的作用下產(chǎn)生的茉莉酸與其受體結(jié)合，活化相應(yīng)基因，產(chǎn)生能削弱或阻斷害蟲消化道內(nèi)蛋白酶的蛋白酶抑制素，使害蟲厭食或消化不良而死[32-33]。Zm00001d044201含有 1個(gè)非同義 SNP，將221位的丙氨酸突變成絲氨酸，編碼ATP依賴Zn2+蛋白酶FTSH11，僅位于葉綠體和線粒體，研究表明，在高于30℃的環(huán)境中，擬南芥FtSH11突變體的生長(zhǎng)發(fā)育受到抑制，而野生型植株不受影響，故該基因產(chǎn)物可能和植物耐熱相關(guān)[34]?；騔m00001d031449和Zm00001d044201可能因自然選擇而發(fā)生非同義突變，但編碼產(chǎn)物對(duì)植物生長(zhǎng)和適應(yīng)性存在正調(diào)節(jié)作用?；騔m00001d043581編碼產(chǎn)物是一類類枯草桿菌蛋白酶，與病原體識(shí)別及互作，同時(shí)被病原體誘導(dǎo)，啟動(dòng)細(xì)胞程序性死亡和超敏反應(yīng)，對(duì)植物產(chǎn)生免疫反應(yīng)[35]。綜上，這些抗性基因的發(fā)現(xiàn)說明秸稈細(xì)胞壁是保護(hù)植株免受生物和非生物傷害的天然屏障。

秸稈消化率性狀在轉(zhuǎn)錄水平上也受到調(diào)控。位于3'非翻譯區(qū)的Zm00001d028777是ZmMs9基因，編碼R2/R3植物特異性MYB轉(zhuǎn)錄因子，在植物次生代謝中起重要作用，與1號(hào)染色體上P1基因有關(guān)，P1基因和調(diào)節(jié)基因C1有較高的同源性，通過調(diào)節(jié)A1和C2基因，參與玉米黃酮醇生物合成過程，可能和防御紫外線相關(guān)[36-37]。

綜上，玉米秸稈細(xì)胞壁是玉米健壯生長(zhǎng)的必需結(jié)構(gòu)，也是草食牲畜主要的纖維和能量來源，理解秸稈細(xì)胞壁合成機(jī)制和玉米生長(zhǎng)適應(yīng)性機(jī)制對(duì)找到有關(guān)消化率性狀的基因有幫助。本文初探了青貯玉米秸稈48h IVDMD的遺傳機(jī)制，為分子改良纖維含量和消化率提供了數(shù)據(jù)支撐，也為優(yōu)質(zhì)青貯玉米選育提供了理論參考。

4 結(jié)論

利用6 276 612個(gè)高質(zhì)量SNPs（MAF＞0.05）對(duì)341份玉米自交系秸稈48h IVDMD進(jìn)行全基因組關(guān)聯(lián)分析，兩個(gè)環(huán)境中采用FarmCPU模型在所有染色體上檢測(cè)到153個(gè)SNPs位點(diǎn)（P＜1.0×10-6）和秸稈48h IVDMD顯著相關(guān)，共挖掘38個(gè)候選基因，主要涉及細(xì)胞生長(zhǎng)發(fā)育、防御反應(yīng)、信號(hào)轉(zhuǎn)導(dǎo)和細(xì)胞代謝等相關(guān)生物學(xué)功能。Zm00001d009009基因和Zm00001d028874基因可能分別與玉米秸稈纖維合成和初生壁合成有關(guān)。