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

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

玉米是重要的糧食和飼料作物，目前已成為我國(guó)第一大作物。穗發(fā)芽（Preharvest sprouting），又稱胎萌（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)分為兩類：第一類突變體包括vp1/vp4、vp6、vp8、vp10/vp13、vp14和vp15，這些突變體的胚乳和幼苗顏色不受影響；第二類突變體包括vp2、vp5、vp7/ps1、vp9、vp12/lw2、y9、w3和rea1，這些突變體的類胡蘿卜素和葉綠素合成受到影響。目前已知，這些穗發(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），又稱高通量測(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）提取籽粒總RNA，在北京貝瑞和康生物技術(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 階段采集熒光，并添加溶解曲線，結(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 在線分析網(wǎng)站（http：//systemsbiology.cau.edu.cn/agriGOv2/index.php）進(jìn)行分析，每條條目至少有3 個(gè)基因，并且FDR＜0.05 作為篩選標(biāo)準(zhǔn)。

京都基因與基因組百科全書（Kyoto Encyclopedia of Genes and Genomes，KEGG）數(shù)據(jù)庫(kù)可以系統(tǒng)分析基因的代謝通路，本研究利用KOBAS3.0［33，34］在線分析網(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 功能分類將差異表達(dá)基因進(jìn)行功能分類注釋，1 791 個(gè)差異基因得到注釋，分布于699 個(gè)GO 分類條目。在生物學(xué)過(guò)程（biological process）、細(xì)胞組分（cellular component）和分子功能（molecular function）中所占的比例分別為 49.6%（347）、11.9%（83）、38.5%（269）（圖4）。進(jìn)一步對(duì)差異基因的生物進(jìn)程類別富集分析發(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ò)程中，?；o酶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 類別的相關(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 分類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ā)芽的水稻種子，類似結(jié)果在小麥中也見報(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ù)信息。