郭宇航，鄭軍，李潤(rùn)植*

（1.山西農(nóng)業(yè)大學(xué)農(nóng)學(xué)院，山西太谷030801；2.中國(guó)農(nóng)業(yè)科學(xué)院作物科學(xué)研究所，北京100081）

玉米是重要的糧食和飼料作物，目前已成為我國(guó)第一大作物。穗發(fā)芽（Preharvest sprouting），又稱(chēng)胎萌（vivipary），是指種子收獲前在母體上發(fā)芽的現(xiàn)象。穗發(fā)芽的種子部分營(yíng)養(yǎng)和貯藏物質(zhì)被分解消耗，導(dǎo)致其品質(zhì)和播種質(zhì)量下降［1］。我國(guó)玉米、水稻、小麥在收獲期經(jīng)常遭遇陰雨連綿的天氣，直接導(dǎo)致穗發(fā)芽，生產(chǎn)造成重大損失。因此，研究穗發(fā)芽相關(guān)基因及調(diào)控網(wǎng)絡(luò)對(duì)培育耐穗發(fā)芽作物品種，保障糧食生產(chǎn)具有重要意義。

目前，在玉米中已經(jīng)發(fā)現(xiàn)了十多個(gè)穗發(fā)芽突變體。這些突變體根據(jù)Robertson［2］提出的標(biāo)準(zhǔn)分為兩類(lèi)：第一類(lèi)突變體包括vp1/vp4、vp6、vp8、vp10/vp13、vp14和vp15，這些突變體的胚乳和幼苗顏色不受影響；第二類(lèi)突變體包括vp2、vp5、vp7/ps1、vp9、vp12/lw2、y9、w3和rea1，這些突變體的類(lèi)胡蘿卜素和葉綠素合成受到影響。目前已知，這些穗發(fā)芽突變體中有9 個(gè)已克隆到基因，包括vp1［3］、vp8［4］、vp10/vp13［5］、vp14［6］、vp15［7］、vp5［8］、vp7/ps1［9］、y9［10］、vp9［11］。VP1 對(duì)玉米種子成熟至關(guān)重要［12］。vp1突變體對(duì)脫落酸（ABA）敏感性降低，不能正常成熟，并引起種子發(fā)芽。玉米穗發(fā)芽基因VP1和擬南芥的ABI3是同源基因，其編碼蛋白是具有多結(jié)構(gòu)域的轉(zhuǎn)錄因子，行使激活或者抑制功能取決于啟動(dòng)子的上下游環(huán)境［12～15］。VP1 的 3 種蛋白結(jié)構(gòu)域 B1，B2，B3 在不同植物中是高度保守的［12，16～21］。VP1 通過(guò) B3 結(jié)構(gòu)域結(jié)合到Sph/RY 元件激活花青素調(diào)節(jié)基因C1基因的表達(dá)［22］。玉米vp1突變體無(wú)法激活C1，進(jìn)而影響糊粉層著色［3］。VP1 還能抑制糊粉層細(xì)胞中發(fā)芽相關(guān)基因如淀粉酶的表達(dá)，表明VP1 也能行使轉(zhuǎn)錄抑制的功能［14，23，24］。

第二代測(cè)序技術(shù)（Next-generation sequencing，NGS），又稱(chēng)高通量測(cè)序技術(shù)（High-throughput sequencing），具有高通量、高準(zhǔn)確性、低成本、耗時(shí)短的特點(diǎn)，廣泛應(yīng)用于各物種的基因組測(cè)序與轉(zhuǎn)錄組分析。例如利用轉(zhuǎn)錄組測(cè)序（RNA-Seq）數(shù)據(jù)，在全基因組的水平上發(fā)掘?qū)τ衩鬃蚜Ｎ镔|(zhì)積累過(guò)程具有重要影響的功能基因［25］。

本研究利用實(shí)驗(yàn)室新發(fā)現(xiàn)的Vp1的等位突變體vp-like8［26］的雜合 體與自交 系 Mo17 和 Zheng58構(gòu)建的F2分離群體，利用生物信息學(xué)比較分析穗發(fā)芽和野生型的轉(zhuǎn)錄組，尋找差異表達(dá)基因（differentially expressed genes，DEGs），通過(guò)基因本體（Gene Ontology，GO）分析和 KEGG 代謝通路富集分析，為揭示玉米vp-like8突變體穗發(fā)芽相關(guān)基因的表達(dá)網(wǎng)絡(luò)提供理論依據(jù)。

1 材料和方法

1.1 試驗(yàn)材料

在玉米收獲時(shí)，玉米穗發(fā)芽的籽粒已干枯致死，因此穗發(fā)芽突變體（viviparous）vp-like8只能以雜合體的形式進(jìn)行保存。利用vp-like8雜合體分別雜交自交系Mo17 和Zheng58，然后再自交，獲得具有穗發(fā)芽表型分離的F2果穗作為后續(xù)轉(zhuǎn)錄組測(cè)序的樣品。

1.2 總RNA 提取及測(cè)序

在授粉后30 d，可在果穗上看到明顯的穗發(fā)芽表型。此時(shí)，隨機(jī)挑選有穗發(fā)芽表型分離的F2果穗，發(fā)芽籽粒和正常籽粒各取40 粒，剝?nèi)シN皮，分別混池Mo17 背景的突變體（Mo17_mut）和野生型（Mo17_wt） ， Zheng58 背 景 的 突 變 體（Zheng58_mut）和野生型（Zheng58_wt）等 4 份實(shí)驗(yàn)樣品，由于利用了不同背景的材料，故每種材料不再單獨(dú)做生物學(xué)重復(fù)。液氮研磨樣品，利用植物總RNA 提取試劑盒（北京天根生化科技有限公司，Cat#DP432）提取籽?？俁NA，在北京貝瑞和康生物技術(shù)有限公司利用HiSeq2000 儀器完成轉(zhuǎn)錄組測(cè)序。

1.3 差異表達(dá)基因的鑒定

測(cè)序得到的原始數(shù)據(jù)，經(jīng)過(guò)去除接頭序列和低質(zhì)量測(cè)序片段之后利用軟件Tophat2 v2.1.0［27］比對(duì)到玉米B73 參考基因組上（B73.AGPv3.27，https：//www.maizegdb.org），所 用參數(shù)“-i 5 -I 60000 --no-novel-juncs --microexon-search -r 100 --mate-std-dev 75”，其它參數(shù)保持默認(rèn)，得到mapped data。 利 用 samtools［28］對(duì) bam 文件進(jìn)行過(guò)濾，設(shè)置參數(shù)“-q 50”。之后利用 HTSeq v0.9.1［29］的htseq-count 腳本計(jì)算出唯一比對(duì)到參考基因組上的 reads 數(shù)，參數(shù)為“-f bam -q -a 10 -m union --nonunique none -s no -r pos”。隨機(jī)挑選基因，使用可視化軟件 Integrative Genomics Viewer（IGV）［30］對(duì)HTSeq-count 統(tǒng)計(jì)計(jì)算的read counts 進(jìn)行驗(yàn)證。最后使用R 語(yǔ)言DESeq2［31］包進(jìn)行差異表達(dá)分析。計(jì)算每個(gè)基因在突變體與野生型中表達(dá)差異倍數(shù)（Fold Change，F(xiàn)C），以（|log2FC|＞1）并且差異顯著性padj（p 值經(jīng)過(guò)多重檢驗(yàn)校正后的值）＜0.01作為判斷基因表達(dá)差異是否顯著的閾值。

1.4 qRT-PCR 驗(yàn)證

隨機(jī)挑選7 個(gè)差異表達(dá)基因（G R M Z M 2 G 133398；GRMZM2G382534；GRMZM2G082487；GRMZM2G123107；GRMZM2G126261；GRMZM 2G017290；GRMZM5G835903）進(jìn)行驗(yàn)證。將 1.2中保留下的RNA 樣品，利用反轉(zhuǎn)錄試劑盒（北京天根生化科技有限公司，Code# AH311-02）反轉(zhuǎn)錄為 cDNA。利用 Primer3［32］設(shè)計(jì)特異性引物，以GAPDH 為內(nèi)參（表 1），利用 ABI 7300 實(shí)時(shí)定量PCR 儀進(jìn)行測(cè)定。qRT-PCR 體系為Dye Ⅰ 0.4 μL ，Q-mix 10 μL，ddH2O 7.8 μL，F(xiàn)/R 引物各 0.4 μL，cDNA 1 μL，qRT-PCR 反應(yīng)條件為 95 ℃預(yù)變性 2 min，95 ℃變性 10 s，58 ℃退火 20 s，72 ℃延伸30 s，設(shè)置 40 個(gè)循環(huán)，在 72 ℃延伸 30 s 階段采集熒光，并添加溶解曲線(xiàn)，結(jié)果用2-ΔΔCt計(jì)算穗發(fā)芽籽粒與正常籽粒中基因表達(dá)差異。試驗(yàn)包含2 次生物學(xué)重復(fù)，每個(gè)生物學(xué)重復(fù)進(jìn)行3 次技術(shù)重復(fù)。

1.5 差異表達(dá)基因的注釋

基因本體（Gene Ontology，GO）分析是根據(jù)GO 功能顯著性富集來(lái)分析差異表達(dá)基因的功能，本研究利用AgriGo v2.0 在線(xiàn)分析網(wǎng)站（http：//systemsbiology.cau.edu.cn/agriGOv2/index.php）進(jìn)行分析，每條條目至少有3 個(gè)基因，并且FDR＜0.05 作為篩選標(biāo)準(zhǔn)。

京都基因與基因組百科全書(shū)（Kyoto Encyclopedia of Genes and Genomes，KEGG）數(shù)據(jù)庫(kù)可以系統(tǒng)分析基因的代謝通路，本研究利用KOBAS3.0［33，34］在線(xiàn)分析網(wǎng)站（http：//kobas.cbi.pku.edu.cn/index.php）進(jìn)行分析，P 值經(jīng) Benjamini and Hochberg 方法校正后，篩選Corrected P-Value＜0.01的通路作為顯著富集的通路。

表1 qRT-PCR 所用引物Table 1 Primers used in qRT-PCR

2 結(jié)果與分析

2.1 序列比對(duì)結(jié)果

在本研究中，Mo17 背景的突變體（Mo17_mut）和野生型（Mo17_wt），Zheng58 背景的突變體（Zheng58_mut）和野生型（Zheng58_wt）等4 份樣品的測(cè)序原始數(shù)據(jù)，經(jīng)過(guò)去除接頭序列和低質(zhì)量測(cè)序片段，之后通過(guò)Tophat2 比對(duì)到玉米B73 參考基因組（B73.AGPv3.27），比對(duì)率均達(dá)80%以上（表2），可以用于進(jìn)一步分析。對(duì)于每份數(shù)據(jù)的比對(duì)結(jié)果，通過(guò)HTSeq-count 對(duì)其片段計(jì)數(shù)（read counts），以基因組注釋為基礎(chǔ)，計(jì)算每個(gè)基因區(qū)域中比對(duì)上的片段數(shù)，以此作為每個(gè)基因的原始表達(dá)數(shù)據(jù)。以GRMZM2G000743 基因?yàn)槔?，利用IGV 對(duì)HTSeq-count 統(tǒng)計(jì)計(jì)算的read counts 進(jìn)行驗(yàn)證，在 Mo17_mut 材料中的 read counts 為 12，IGV 結(jié)果顯示有 12 個(gè)成對(duì)的 reads 比對(duì)到該基因區(qū)域（圖1），與計(jì)算得到的結(jié)果一致，驗(yàn)證了HTSeq-count 統(tǒng)計(jì)計(jì)算read counts 結(jié)果的準(zhǔn)確性。

表2 Tophat 比對(duì)結(jié)果Table 2 The tophat mapped results

2.2 差異表達(dá)基因分析

分別以Mo17 背景的野生型和Zheng58 背景的野生型數(shù)據(jù)為對(duì)照，計(jì)算每個(gè)基因在突變體與野生型中表達(dá)量的比值，即表達(dá)差異倍數(shù)。根據(jù)DESeq2 分 析結(jié)果，以|log2FC|＞1 并且 padj＜0.01 為條件進(jìn)行篩選（圖2A），共篩選出27 50 個(gè)DEGs。其中，1 812 個(gè)基因表達(dá)上調(diào)，占總DEGs 的65.89%左右，938 個(gè)基因表達(dá)下調(diào)，占總DEGs 的34.11%左右（圖2B）。Vp1基因（GeneID：GRMZM2G1333 98）確實(shí)存在于差異基因中，并且在突變體中表達(dá)下調(diào)（表3）。

2.3 qRT-PCR 驗(yàn)證

本研究從發(fā)現(xiàn)的DEGs 中隨機(jī)選取7 個(gè)基因進(jìn)行熒光定量PCR 驗(yàn)證，通過(guò)與差異表達(dá)分析結(jié)果比較，雖然差異基因在上調(diào)和下調(diào)倍數(shù)上存在差異，但是總體變化趨勢(shì)是一致的，驗(yàn)證了RNA-seq 分析結(jié)果的可靠性（圖3）。

2.4 差異表達(dá)基因GO 富集分析

圖1 利用IGV 驗(yàn)證read counts 結(jié)果Fig.1 Using IGV to validate results of read counts

圖2 差異表達(dá)基因火山圖和數(shù)量Fig.2 Vocalno plot and number of DEGs

通過(guò)GO 功能分類(lèi)將差異表達(dá)基因進(jìn)行功能分類(lèi)注釋?zhuān)? 791 個(gè)差異基因得到注釋?zhuān)植加?99 個(gè)GO 分類(lèi)條目。在生物學(xué)過(guò)程（biological process）、細(xì)胞組分（cellular component）和分子功能（molecular function）中所占的比例分別為 49.6%（347）、11.9%（83）、38.5%（269）（圖4）。進(jìn)一步對(duì)差異基因的生物進(jìn)程類(lèi)別富集分析發(fā)現(xiàn)（圖5），碳水化合物代謝過(guò)程（carbohydrates metabolic process）、細(xì)胞氮化合物代謝過(guò)程（cellular nitrogen compound metabolic process）、脂質(zhì)代謝過(guò)程（lipid metabolic process）、光合作用（photosynthesis）和對(duì)刺激響應(yīng)（response to stress）顯著富集。穗發(fā)芽需要大量能量，涉及多種代謝過(guò)程，我們發(fā)現(xiàn)在碳水化合物代謝過(guò)程中，α-淀粉水解酶基因和其他糖基水解酶基因在突變體中表達(dá)上調(diào)；在細(xì)胞氮化合物代謝過(guò)程中，多種與氮代謝相關(guān)的基因表達(dá)上調(diào)；在脂質(zhì)代謝過(guò)程中，酰基輔酶A 氧化酶基因表達(dá)上調(diào)，并且與ABA 合成相關(guān)的β-胡蘿卜素羥化酶基因表達(dá)下調(diào)；在光合作用中，多種光合作用相關(guān)的基因表達(dá)上調(diào)；在脅迫響應(yīng)過(guò)程中，多個(gè)過(guò)氧化物酶基因表達(dá)上調(diào)（表3）。這些顯著富集的GO 類(lèi)別的相關(guān)基因共同作用，為玉米籽粒發(fā)芽提供能量和保護(hù)，與玉米穗發(fā)芽是緊密相關(guān)的。

圖3 qRT-PCR 驗(yàn)證隨機(jī)挑選的DEGsFig.3 Validation of the random selected DEGs by qRT-PCR

表3 部分差異表達(dá)基因信息Table 3 Part of differentially expressed genes’informatation

續(xù)表

2.5 差異基因代謝途徑分析

圖4 差異表達(dá)基因的GO 分類(lèi)Fig.4 Go annotation of differentially expressed genes

圖5 差異表達(dá)基因的生物進(jìn)程（有向無(wú)環(huán)圖）Fig.5 Biological process of differentially expressed genes（Directed acyclic graph）

利用KEGG 數(shù)據(jù)庫(kù)對(duì)差異表達(dá)基因所參與的代謝途徑進(jìn)行顯著性富集分析，結(jié)果顯示，這些差異基因分布于110 個(gè)代謝通路中，分別為代謝通路（Metabolic pathways）、次生代謝產(chǎn)物的生物合成（Biosynthesis of secondary metabolites）、苯丙素生物合成（Phenylpropanoid biosynthesis）、光合作用（Photosynthesis）、植物激素信號(hào)轉(zhuǎn)導(dǎo)（Plant hormone signal transduction）、碳代謝（Carbon metabolism）、淀粉和蔗糖代謝（Starch and sucrose metabolism）等。其中代謝通路和次級(jí)代謝產(chǎn)物生物合成通路富集到的基因是最多的，分別有314 和211 個(gè)基因被富集（表4）。植物內(nèi)源激素與植物生長(zhǎng)發(fā)育密切相關(guān)，有48 個(gè)基因富集到植物激素信號(hào)傳導(dǎo)通路（圖6），圖中紅框?yàn)椴町惐磉_(dá)基因在通路中位置，涉及生長(zhǎng)素（Auxin）、細(xì)胞分裂素（Cytokinin，CTK）、赤霉素（Gibberellin，GA）、脫落酸（Abscisic acid，ABA）、乙烯（Ethylene，ETH）等，其中ABA 直接影響種子休眠，我們發(fā)現(xiàn)蛋白磷酸酶2C基因（PP2C）、絲氨酸/蘇氨酸蛋白激酶基因（SnRK2）、ABA 響應(yīng)元件結(jié)合因子（ABF）在該通路中受到影響。這些受影響的代謝通路中的基因?yàn)檠芯坑衩姿氚l(fā)芽相關(guān)表達(dá)網(wǎng)絡(luò)提供了豐富的基因資源。

表4 顯著富集的KEGG 代謝通路Table 4 Highly enriched KEGG terms

3 討論與結(jié)論

玉米穗發(fā)芽是一個(gè)復(fù)雜的生物學(xué)過(guò)程，影響因素較多，往往各因素相互影響。在本研究中，借助本實(shí)驗(yàn)室新發(fā)現(xiàn)的玉米穗發(fā)芽vp-like8突變體［26］進(jìn)行研究。趙燕［35］對(duì)玉米穗發(fā)芽差異基因分析發(fā)現(xiàn)光合作用通路最顯著富集，我們發(fā)現(xiàn)與光合作用相關(guān)的基因大部分都表達(dá)上調(diào)（表3），說(shuō)明光合作用與穗發(fā)芽確實(shí)存在緊密聯(lián)系?；ㄇ嗨卣{(diào)節(jié)基因C1與糊粉層著色相關(guān)，本研究中，該基因表達(dá)下調(diào)（表 3），與前人研究結(jié)果一致［3］。VP1 能夠抑制糊粉層中和發(fā)芽相關(guān)的α-淀粉酶的表達(dá)［14，23］，本研究中我們發(fā)現(xiàn)在淀粉和蔗糖代謝途徑中，α-淀粉酶基因表達(dá)上調(diào)（表3）。這些結(jié)果印證了我們分析的可靠性，同時(shí)為利用轉(zhuǎn)錄組測(cè)序技術(shù)研究玉米穗發(fā)芽機(jī)制提供了支持。

研究表明，種子內(nèi)源激素 ABA［36，37］對(duì)種子發(fā)芽有重要作用。Suzuki 等人研究發(fā)現(xiàn)［38］，VP1 能影響ABA 信號(hào)傳導(dǎo)中的bZIP基因和蛋白磷酸酶2C基因（PP2C）的表達(dá)。本研究中，同樣發(fā)現(xiàn)這些基因在ABA 信號(hào)傳導(dǎo)通路中是差異表達(dá)的，并且還發(fā)現(xiàn)絲氨酸/蘇氨酸蛋白激酶基因（SnRK2）也在ABA 信號(hào)傳導(dǎo)中受到影響（表3）。同時(shí)，還發(fā)現(xiàn)生長(zhǎng)素轉(zhuǎn)運(yùn)載體，生長(zhǎng)素響應(yīng)蛋白相關(guān)基因的表達(dá)也有顯著差異（表3）。綜上，這些差異基因表達(dá)水平的變化為穗發(fā)芽相關(guān)調(diào)控網(wǎng)絡(luò)的研究提供了參考。

圖6 植物內(nèi)源激素信號(hào)傳導(dǎo)Fig.6 Plant hormone signal transduction

其它作物穗發(fā)芽相關(guān)研究也較為豐富。廖泳祥等［39］通過(guò)比較不同的水稻材料種子，發(fā)現(xiàn)在不耐穗發(fā)芽的水稻材料的種子中α-淀粉酶含量明顯高于耐穗發(fā)芽的水稻種子，類(lèi)似結(jié)果在小麥中也見(jiàn)報(bào)道［40］。孫果忠等［41］人對(duì)不同生理狀態(tài)的小麥離體胚萌發(fā)過(guò)程研究發(fā)現(xiàn)，ABA 可以抑制未成熟胚與休眠成熟胚中α-淀粉酶-1 的表達(dá)，但在不休眠成熟胚中則誘導(dǎo)α-淀粉酶-1 表達(dá)。此外，也有部分學(xué)者認(rèn)為，種皮顏色、子粒和穗部其它物理因素與穗發(fā)芽存在關(guān)聯(lián)。

本研究利用轉(zhuǎn)錄組測(cè)序技術(shù)系統(tǒng)分析了玉米vp-like8突變體和野生型的籽粒轉(zhuǎn)錄組。通過(guò)比較突變體和野生型間的表達(dá)水平，發(fā)現(xiàn)了2 750 個(gè)差異表達(dá)基因，其中1 812 個(gè)基因在突變體中上調(diào)表達(dá)，938 個(gè)下調(diào)表達(dá)。通過(guò)功能注釋分析發(fā)現(xiàn)，這些差異基因主要涉及代謝通路、光合作用、植物激素信號(hào)傳導(dǎo)等途徑，表明這些途徑在穗發(fā)芽過(guò)程中具有一定作用。本文通過(guò)轉(zhuǎn)錄組測(cè)序的方法初步研究了vp-like8基因的表達(dá)差異，為深入了解Vp1基因功能和玉米穗發(fā)芽的機(jī)制及基因調(diào)控網(wǎng)絡(luò)提供了數(shù)據(jù)信息。