玉米花藥不同發(fā)育時(shí)期的轉(zhuǎn)錄組分析及其差異表達(dá)基因的鑒定

陳慧， 孫曉靖， 郭瑤晴， 連玉杰， 張雪海， 湯繼華， 陳曉陽(yáng)

(河南農(nóng)業(yè)大學(xué)農(nóng)學(xué)院，河南 鄭州 450046)

雄性育性影響玉米種子形成和生存繁衍，雄性不育系在玉米雜交制種和雜種優(yōu)勢(shì)利用中發(fā)揮了重要作用[1-2]。雄性育性的形成依賴于花藥正常發(fā)育，其發(fā)育過(guò)程受一系列遺傳因子精密調(diào)控。依據(jù)花藥發(fā)育的細(xì)胞學(xué)特征，可將玉米花藥發(fā)育劃分為14個(gè)時(shí)期，主要分為減數(shù)分裂前期(S1～S6)，減數(shù)分裂期(S7～S8)和減數(shù)分裂后期(S9～S14)3個(gè)階段[3]。其中，S6期花藥內(nèi)形成小孢子母細(xì)胞，開始進(jìn)入減數(shù)分裂階段；S8b期花藥完成減數(shù)分裂過(guò)程，形成四分體細(xì)胞；S13期花藥已形成成熟花粉粒，包裹花藥的內(nèi)外稃張開。

花藥角質(zhì)層和花粉壁是雄性育性形成過(guò)程中重要屏障?；ㄋ幗琴|(zhì)層覆蓋于花藥表面，防止各種生物脅迫和非生物脅迫破壞花藥發(fā)育[4]?；ㄋ幗琴|(zhì)層主要由角質(zhì)和蠟質(zhì)組成，其中角質(zhì)包含C16和C18脂肪酸、羥基脂肪酸、環(huán)氧脂肪酸、二元脂肪酸等物質(zhì)[5]，蠟質(zhì)由脂肪酸、烴類、醇類、酯類等物質(zhì)構(gòu)成[6]?；ǚ郾谖挥谛℃咦?花粉)表面，保護(hù)小孢子正常發(fā)育和成熟花粉粒形成[7]。利用雄性不育突變體已鑒定到部分參與玉米花藥角質(zhì)層和花粉壁發(fā)育的關(guān)鍵基因。例如，脂肪酸羥化酶基因APV1和MS26[8-9]，羥基脂肪酸脫氫酶基因IPE1/ZmMs20[10-11]，脂肪酰還原酶基因MS6021/ZmMs25/ZmFAR1[12-14]，GDSL脂肪酶基因ZmMs30和IPE2[15-17]，甘油-3-磷酸?；D(zhuǎn)移酶基因ZmMs33[18]，脂類物質(zhì)轉(zhuǎn)蛋白基因MS2/ZmABCG26[14,19]，轉(zhuǎn)錄調(diào)控基因ZmMs7，ZmbHLH51，ZmbHLH122，ZmMYB84，ZmMYB33-1/2和ZmPHD11[20-22]。

轉(zhuǎn)錄組測(cè)序技術(shù)(RNA-Seq)是利用高通量測(cè)序技術(shù)對(duì)生物進(jìn)行轉(zhuǎn)錄組分析,是一種檢測(cè)全基因組基因表達(dá)譜的有效方法[23]。通過(guò)比較不同發(fā)育時(shí)期基因表達(dá)譜，可以挖掘到影響植物生長(zhǎng)發(fā)育的關(guān)鍵基因。WAN等[24]利用玉米自交系Oh43和W23花藥轉(zhuǎn)錄組數(shù)據(jù)，鑒定到125個(gè)參與雄性育性形成的脂類物質(zhì)代謝相關(guān)基因。JIANG等[22]分析3個(gè)自交系S5～S12時(shí)期花藥轉(zhuǎn)錄譜，發(fā)現(xiàn)1 100個(gè)轉(zhuǎn)錄因子在3個(gè)自交系中均表現(xiàn)出時(shí)期差異表達(dá)，進(jìn)一步通過(guò)CRISPR/Cas9技術(shù)驗(yàn)證了12個(gè)轉(zhuǎn)錄因子的功能。然而，目前尚未利用轉(zhuǎn)錄組測(cè)序技術(shù)挖掘影響玉米花藥角質(zhì)層和花粉壁發(fā)育的關(guān)鍵基因。本研究對(duì)玉米自交系B73的3個(gè)發(fā)育時(shí)期花藥和苗期葉片進(jìn)行轉(zhuǎn)錄組測(cè)序，通過(guò)GO和KEGG等生物學(xué)信息分析鑒定影響花藥角質(zhì)層和花粉壁形成的關(guān)鍵基因及其參與的代謝途徑，為玉米雄性育性調(diào)控機(jī)制的解析奠定基礎(chǔ)，同時(shí)為雄性不育系的創(chuàng)制提供潛在靶標(biāo)位點(diǎn)。

1 材料與方法

1.1 試驗(yàn)材料

供試的玉米自交系為B73，2020年4月種植于河南省新鄉(xiāng)市原陽(yáng)縣河南農(nóng)業(yè)大學(xué)原陽(yáng)試驗(yàn)基地(35.1°N，113.6°E)。依據(jù)花藥長(zhǎng)度和碘-碘化鉀染色判斷花藥發(fā)育時(shí)期，選取S6期(記為B1)，S8b期(記為B2)，S13期(記為B3)花藥，以及苗期葉片(記為B4)于-80 ℃保存。每組樣品包含4個(gè)生物學(xué)重復(fù)。

1.2 RNA提取、文庫(kù)構(gòu)建和測(cè)序

利用天根試劑盒提取樣品總RNA。Nanodrop檢測(cè)RNA純度，瓊脂糖凝膠電泳分析RNA降解程度及是否有污染，Qubit2.0對(duì)RNA濃度進(jìn)行精確定量，Agilent 2100精確檢測(cè)RNA完整性。樣品檢測(cè)合格后，按如下流程構(gòu)建文庫(kù)：(1)利用帶有Oligo(dT)磁珠富集mRNA。(2)加入fragmentation buffer將mRNA打斷成短片段。(3)以mRNA為模板，用六堿基隨機(jī)引物合成第1條cDNA鏈，加入緩沖液、dNTPs、RNase H和DNA polymerase I合成第2條cDNA鏈。(4)雙鏈cDNA末端修復(fù)且在3’端加A。(5)QiaQuick PCR試劑盒純化并加EB緩沖液洗脫，末端修復(fù)，連接測(cè)序接頭。(6)瓊脂糖凝膠電泳選擇片段大小，PCR擴(kuò)增富集cDNA。文庫(kù)質(zhì)檢合格后，通過(guò)Illumina Novaseq平臺(tái)進(jìn)行PE150測(cè)序。

1.3 測(cè)序數(shù)據(jù)過(guò)濾、參考基因組比對(duì)和基因表達(dá)定量

測(cè)序得到原始序列raw reads，通過(guò)數(shù)據(jù)過(guò)濾得到干凈、有效、高質(zhì)量的clean reads。數(shù)據(jù)過(guò)濾步驟：(1)當(dāng)單端測(cè)序reads中N>3時(shí)，去除此對(duì)paired reads；(2)當(dāng)單端測(cè)序read中質(zhì)量值低于5的堿基比例≥20%時(shí)，去除此對(duì)paired reads；(3)去除adapter序列時(shí)，adapter序列至少要匹配8 bp。采用Hisat2軟件把過(guò)濾后的reads與玉米自交系B73參考基因組(V4版本)比對(duì)[25]。利用featureCounts軟件對(duì)每個(gè)樣品進(jìn)行基因水平的定量，合并得到所有樣本的基因表達(dá)量矩陣[26]。

1.4 差異表達(dá)基因鑒定和功能富集分析

利用edge R軟件分析基因在各樣品中的差異表達(dá)，計(jì)算差異表達(dá)的Pvalue和Padj值，Padj是校正后的Pvalue[27]。差異表達(dá)基因篩選條件為：Padj<0.05，|log2FoldChange(FC)|>1。應(yīng)用topGO軟件進(jìn)行差異表達(dá)基因GO(gene ontoloy)富集分析。應(yīng)用KOBAS軟件進(jìn)行差異表達(dá)基因KEGG(kyoto encyclopedia of genes and genomes)通路富集分析[28]。

1.5 qRT-PCR驗(yàn)證

通過(guò)NCBI(https://www.ncbi.nlm.nih.gov/)設(shè)計(jì)qRT-PCR引物，引物序列見表1。利用2×SYBR Green qPCR Premix(Universal)試劑盒在伯樂(lè)C1000 TouchTMThermal Cycler儀器進(jìn)行qRT-PCR分析。以ZmACTIN為內(nèi)參，2-ΔΔCt法計(jì)算相對(duì)表達(dá)量。

表1 qRT-PCR試驗(yàn)所用引物Table 1 Primers used in the qRT-PCR

2 結(jié)果與分析

2.1 玉米花藥發(fā)育不同時(shí)期轉(zhuǎn)錄組測(cè)序數(shù)據(jù)的質(zhì)量分析

通過(guò)Illumina Novaseq平臺(tái)對(duì)16個(gè)RNA樣品文庫(kù)進(jìn)行測(cè)序，測(cè)序共產(chǎn)生456 100 906條clean reads，質(zhì)量值≥30的堿基所占比例(Q30)均大于90%，說(shuō)明測(cè)序數(shù)據(jù)質(zhì)量較好，可用于進(jìn)一步分析。將獲得的有效測(cè)序數(shù)據(jù)比對(duì)到玉米自交系B73參考基因組(V4版本)上，比對(duì)率為92.46%～96.49%，平均比對(duì)率為95.56%(表2)。

表2 玉米花藥發(fā)育不同時(shí)期RNA-Seq數(shù)據(jù)的統(tǒng)計(jì)分析

為驗(yàn)證測(cè)序數(shù)據(jù)的可靠性，分析了8個(gè)玉米雄性育性基因在測(cè)序數(shù)據(jù)中的表達(dá)模式。MS32在玉米花藥發(fā)育早期表達(dá)，調(diào)控絨氈層分化過(guò)程[29]。在測(cè)序數(shù)據(jù)中，MS32基因在S6期花藥中特異性表達(dá)。IPE1/ZmMs20，APV1，IPE2，MS2/ZmABCG26，ZmMs7，ZmMYB84調(diào)控花藥角質(zhì)層和花粉外壁發(fā)育，主要在S8b和S9期花藥表達(dá)[8，10-11，14，17，19-22]。RNA-Seq結(jié)果表明，6個(gè)基因在S8b期花藥中表達(dá)量最高。LBD10被報(bào)道影響花粉后期發(fā)育，通過(guò)檢測(cè)發(fā)現(xiàn)，該基因在S13時(shí)期花藥表達(dá)量最高[22]。上述結(jié)果表明，8個(gè)基因在測(cè)序數(shù)據(jù)中的表達(dá)模式與已報(bào)道基因表達(dá)模式一致。同時(shí)，qRT-PCR分析顯示，8個(gè)基因RNA-Seq和qRT-PCR結(jié)果吻合，說(shuō)明測(cè)序數(shù)據(jù)真實(shí)可靠(圖1)。

FPKM:每千個(gè)堿基的轉(zhuǎn)錄每百萬(wàn)映射讀取的碎片。FPKM: Fragments per kilobase of exon model per million mapped fragments.

2.2 玉米花藥發(fā)育不同時(shí)期差異表達(dá)基因分析

以Padj<0.05和|log2FC|>1作為差異表達(dá)基因的篩選標(biāo)準(zhǔn)，在B1 vs B4，B2 vs B4，B3 vs B4，B2 vs B1，B3 vs B1，B2 vs B3 6組之間分別鑒定到15 900，14 785，17 141，8 115，17 858，16 497個(gè)差異表達(dá)基因(表3)。不同發(fā)育時(shí)期花藥與苗期葉片相比，S13期花藥與苗期葉片差異表達(dá)基因數(shù)最多(17 141)，其次是S6期(15 900)，S8b期最少(14 785)。不同發(fā)育時(shí)期花藥相比，S6與S13期花藥差異表達(dá)基因最多(17 858)，其次是S8b與S13(16 497)，S6與S8b最少(8 115)，說(shuō)明S6與S8b期花藥轉(zhuǎn)錄譜相似。此外，S6與S13期花藥差異表達(dá)基因多于任一時(shí)期花藥與苗期葉片差異表達(dá)基因，表明不同發(fā)育時(shí)期花藥基因表達(dá)差異較大。

表3 3個(gè)時(shí)期花藥和苗期葉片差異表達(dá)基因數(shù)目

2.3 玉米花藥角質(zhì)層發(fā)育相關(guān)基因的鑒定

為挖掘影響玉米花藥角質(zhì)層發(fā)育的候選基因，通過(guò)KEGG方法分析不同發(fā)育時(shí)期花藥差異表達(dá)基因參與的生化代謝途徑。B2 vs B1和B2 vs B3上調(diào)差異表達(dá)基因均富集到角質(zhì)和蠟質(zhì)合成代謝途徑，而兩組下調(diào)差異表達(dá)基因及B1 vs B3差異表達(dá)基因均不能富集到角質(zhì)和蠟質(zhì)合成代謝途徑(圖2)。這表明S8b時(shí)期是玉米花藥角質(zhì)層發(fā)育的關(guān)鍵時(shí)期，該時(shí)期相關(guān)基因的表達(dá)對(duì)角質(zhì)層形成至關(guān)重要。

圖2 B2 vs B1(A)和B2 vs B3(B)上調(diào)差異表達(dá)基因KEGG富集分析

B2 vs B1中8個(gè)差異表達(dá)基因富集到角質(zhì)和蠟質(zhì)合成代謝途徑。Zm00001d049950和Zm00001d048337編碼脂肪酰還原酶，Zm00001d012274和Zm00001d042723編碼ω羥基棕櫚酸阿魏酸?；D(zhuǎn)移酶，Zm00001d025833和Zm00001d017953編碼過(guò)加氧酶，Zm00001d042814和Zm00001d012326編碼細(xì)胞色素單加氧酶CYP86A1(表4)。

表4 B2 vs B1上調(diào)差異表達(dá)基因中角質(zhì)和蠟質(zhì)合成相關(guān)基因

B2 vs B3中15個(gè)差異表達(dá)基因富集到角質(zhì)和蠟質(zhì)合成代謝途徑，其中Zm00001d049950，Zm00001d048337，Zm00001d012274，Zm00001d042723，Zm00001d025833，Zm00001d042814和Zm00001d012326在B2 vs B1中被鑒定到。Zm00001d014055，Zm00001d051932，Zm00001d017251編碼醛脫羧酶，Zm00001d032284和Zm00001d020238編碼ω羥基脂肪酸脫氫酶，Zm00001d051800編碼過(guò)加氧酶，Zm00001d017528和Zm00001d002687分別編碼細(xì)胞色素單加氧酶CYP86A35和CYP86A7(表5)。

表5 B2 vs B3上調(diào)差異表達(dá)基因中角質(zhì)和蠟質(zhì)合成相關(guān)基因

2.4 玉米花粉壁發(fā)育相關(guān)基因的鑒定

花粉壁位于小孢子外表面，其形成起始于四分體時(shí)期的初生外壁。為鑒定影響花粉壁發(fā)育的關(guān)鍵基因，分析了四分體時(shí)期花藥與其他組樣品之間(B2 vs B1，B2 vs B3，B2 vs B4)的差異表達(dá)基因，共鑒定到3 555個(gè)共同差異表達(dá)基因(圖3)。其中，1 348個(gè)基因在S8b時(shí)期花藥中均上調(diào)表達(dá)，344個(gè)基因均下調(diào)表達(dá)。

圖3 B2 vs B1，B2 vs B3，B2 vs B4差異表達(dá)基因韋恩圖

GO分析表明，1 348個(gè)共同上調(diào)差異表達(dá)基因在花粉壁發(fā)育相關(guān)生物過(guò)程顯著富集。11個(gè)基因參與該過(guò)程，其中7個(gè)基因?yàn)橐褕?bào)道的玉米雄性育性基因，分別為脂肪酰還原酶基因MS6021/ZmMs25/ZmFAR1，ABCG轉(zhuǎn)運(yùn)蛋白基因MS2/ZmABCG26，細(xì)胞色素單加氧酶基因APV1和MS26，PHD-finger型轉(zhuǎn)錄因子基因ZmMs7，α-吡喃酮還原酶基因ZmDFR1和ZmDFR2。Zm00001d002342編碼多聚半乳糖醛酸酶，Zm00001d047858編碼異胡豆苷合酶，Zm00001d019478編碼查爾酮合成酶，Zm00001d003702編碼4-香豆酰輔酶A連接酶(表6)。

表6 S8b時(shí)期花藥共同上調(diào)差異表達(dá)基因中花粉壁發(fā)育相關(guān)基因

1 348個(gè)共同上調(diào)差異表達(dá)基因在轉(zhuǎn)錄調(diào)控等分子功能顯著富集，其中191個(gè)基因被注釋為轉(zhuǎn)錄因子和轉(zhuǎn)錄調(diào)節(jié)因子，占比14.17%(191/1 348)。MYB，AP2/ERF-ERF和NAC家族轉(zhuǎn)錄因子數(shù)最多，分別為21，19，17(圖4)。4個(gè)轉(zhuǎn)錄因子基因(ZmMYB84，ZmMs7，MS9，ZmPHD11)為已報(bào)道的玉米雄性育性基因。

圖4 S8b時(shí)期花藥共同上調(diào)差異表達(dá)基因中轉(zhuǎn)錄因子分析

3 結(jié)論與討論

花藥角質(zhì)層形成過(guò)程包含多類物質(zhì)代謝途徑。目前，通過(guò)雄性不育突變體代謝組學(xué)和蛋白生化功能分析，部分代謝途徑被證明參與玉米花藥角質(zhì)層發(fā)育。例如脂肪酸ω羥基化途徑，羥基脂肪酸氧化途徑，脂酰ACP還原途徑等。

本研究通過(guò)不同發(fā)育時(shí)期花藥差異表達(dá)基因KEGG分析，鑒定到調(diào)控玉米花藥角質(zhì)層發(fā)育的新代謝途徑。Zm00001d012274和Zm00001d042723編碼一種乙酰轉(zhuǎn)移酶，可能介導(dǎo)阿魏酰輔酶A和ω羥基棕櫚酸酯化途徑調(diào)控花藥角質(zhì)層發(fā)育。水稻乙酰轉(zhuǎn)移酶DPW2催化阿魏酰輔酶A和ω羥基棕櫚酸形成阿魏酰棕櫚酸酯，參與花藥角質(zhì)形成[30]。環(huán)氧脂肪酸是花藥角質(zhì)層角質(zhì)組成成分，然而其遺傳調(diào)控機(jī)制尚不清楚。Zm00001d025833，Zm00001d017953和Zm00001d051800編碼過(guò)加氧酶，可能通過(guò)脂肪酸環(huán)氧化途徑影響花藥角質(zhì)層發(fā)育。烴類物質(zhì)是花藥角質(zhì)層蠟質(zhì)組成成分，水稻OsCER1已被報(bào)道通過(guò)調(diào)節(jié)超長(zhǎng)鏈烴形成而參與角質(zhì)層蠟質(zhì)合成[31]。Zm00001d014055，Zm00001d051932和Zm00001d017251是水稻OsCER1同源基因，可能通過(guò)影響超長(zhǎng)鏈烴合成途徑而調(diào)控玉米花藥角質(zhì)層發(fā)育。

本研究鑒定到11個(gè)基因在玉米花粉壁發(fā)育過(guò)程中發(fā)揮重要作用，7個(gè)為已報(bào)道的雄性育性相關(guān)基因。其中MS6021/ZmMs25/ZmFAR1，APV1，MS26，ZmMs7基因突變后均表現(xiàn)出花粉外壁發(fā)育缺陷。四聚酮α-吡喃酮還原酶基因ZmDFR1和ZmDFR2單突變體雄性育性下降，其雙突變體完全不育[32]。水稻四聚酮α-吡喃酮還原酶基因OsTKPR1突變破壞花粉外壁發(fā)育過(guò)程[33]。ABCG轉(zhuǎn)運(yùn)蛋白基因MS2/ZmABCG26突變體完全不育，其水稻同源基因OsABCG15參與花粉外壁發(fā)育和雄性育性形成[34]。

多聚半乳糖醛酸酶基因Zm00001d002342催化果膠分子多聚α-(1，4)-聚半乳糖醛酸的裂解。棉花多聚半乳糖醛酸酶GhNSP參與花粉外壁形成過(guò)程[35]。異胡豆苷合酶基因Zm00001d047858調(diào)控萜類吲哚生物堿的合成。玉米MS45和水稻OsSTRL2為異胡豆苷合酶基因，調(diào)控花粉壁形成[36-37]。查爾酮合成酶基因Zm00001d019478以香豆酰輔酶A和丙二酰輔酶A為底物生成查爾酮。擬南芥LAP5和LAP6具有查爾酮合成酶功能，突變體表現(xiàn)為花粉壁和雄性育性異常[38]。4-香豆酰輔酶A連接酶基因Zm00001d003702在香豆酰輔酶A和阿魏酰輔酶A形成中發(fā)揮重要功能，可能通過(guò)芳香族化合物合成途徑調(diào)控花粉壁發(fā)育。

玉米花藥角質(zhì)層和花粉壁發(fā)育受一系列轉(zhuǎn)錄因子調(diào)控。ZmMs7，ZmbHLH51，ZmbHLH122，ZmMYB84，ZmMYB33-1/-2，ZmPHD11已被證明參與花藥角質(zhì)層發(fā)育和雄性育性形成。ZmMs7與核因子Y亞基形成復(fù)合體直接調(diào)節(jié)ZmMT2C表達(dá)和絨氈層細(xì)胞PCD過(guò)程，同時(shí)間接調(diào)節(jié)MS6021和ZmACOS5等基因表達(dá)而參與花藥角質(zhì)層和花粉壁發(fā)育[20-21]。ZmMYB84通過(guò)與ZmMs25啟動(dòng)子結(jié)合調(diào)控花藥角質(zhì)層和花粉外壁形成[13]。本研究鑒定了191個(gè)轉(zhuǎn)錄因子和轉(zhuǎn)錄調(diào)節(jié)因子,這些因子可能參與玉米花藥角質(zhì)層和花粉壁發(fā)育。后續(xù)對(duì)這些基因進(jìn)行功能解析，將完善玉米花藥角質(zhì)層和花粉壁發(fā)育的轉(zhuǎn)錄調(diào)控機(jī)制。