閆嘉穎，王久紅，趙永鋒，賈曉艷，郭晉杰，祝麗英

（河北農(nóng)業(yè)大學(xué) 農(nóng)學(xué)院/國家玉米改良中心河北分中心/河北省作物種質(zhì)資源實(shí)驗(yàn)室，河北 保定 071001）

植物雄性不育是指植物在有性生殖過程中，雄性器官發(fā)育異常，不能產(chǎn)生正常功能的花藥、花粉或雄配子的現(xiàn)象［1］。在棉花、小麥、大白菜等多種作物中，發(fā)現(xiàn)雄性不育系花藥短小干癟，成熟期花藥不開裂，藥囊內(nèi)無花粉或花粉沒有活力［2］。作物育性的形成依賴于花藥正常發(fā)育，而花藥發(fā)育過程受大量基因的精密調(diào)控［3］，任何發(fā)育時(shí)期的功能異常都會(huì)導(dǎo)致花粉敗育。

轉(zhuǎn)錄組測序技術(shù)是對特定細(xì)胞在某一狀態(tài)下轉(zhuǎn)錄出來的所有RNA 進(jìn)行測序分析，能夠反映細(xì)胞中基因的表達(dá)情況及其調(diào)控規(guī)律，在差異表達(dá)基因（Differentially expressed genes，DEGs）篩選、新基因挖掘、分子標(biāo)記開發(fā)、代謝網(wǎng)絡(luò)等方面的研究中廣泛應(yīng)用［4］。目前，多種作物雄性不育系的轉(zhuǎn)錄組測序工作已相繼開展。曹曉聰?shù)龋?］在轉(zhuǎn)錄組水平上，分析了影響棉花CMS 的基因表達(dá)模式，探討與CMS 相關(guān)的一些潛在基因和代謝途徑。閆東等［6］以紫花苜蓿雄性不育系MS-JN1A 及保持系MS-JN1B 為材料，對5 個(gè)發(fā)育時(shí)期的花藥進(jìn)行轉(zhuǎn)錄組測序，篩選DEGs 進(jìn)行富集分析和功能注釋，探討與紫花苜蓿雄性不育相關(guān)的代謝通路。Li 等［7］比較了CMS-C 型不育系C48-2 及其保持系N48-2在花粉母細(xì)胞和單核期的轉(zhuǎn)錄結(jié)果，檢測到2 069個(gè)與小孢子敗育過程有關(guān)的DEGs，討論與玉米CMS-C 型不育相關(guān)的基因和代謝途徑。

玉米是雜種優(yōu)勢利用最早也是最成功的作物之一，雄性不育在玉米雜交制種和雜種優(yōu)勢利用中發(fā)揮了重要作用［8］。CMS-S 型不育系是玉米雄性不育中資源最豐富的一類，雖然目前對其已有一些細(xì)胞形態(tài)和生理、生化機(jī)制方面的研究，但至今尚未完全闡明其不育機(jī)制，仍需從分子層面進(jìn)行深入研究。鄭58cms-Q1261 是河北農(nóng)業(yè)大學(xué)玉米育種課題組選育的CMS-S 型雄性不育系，其不育特性穩(wěn)定，通過細(xì)胞學(xué)觀察發(fā)現(xiàn)不育系小孢子敗育的關(guān)鍵時(shí)期是單核晚期［9］。本研究以不育系鄭58cms-Q1261和保持系鄭58 單核晚期的花藥為研究對象，利用轉(zhuǎn)錄組測序和生物信息學(xué)分析方法鑒定DEGs，通過GO 和KEGG 富集分析DEGs 參與的生物學(xué)過程和代謝途徑，為CMS-S 型雄性不育分子機(jī)理的深入研究提供參考。

1 材料與方法

1.1 試驗(yàn)材料

選用河北農(nóng)業(yè)大學(xué)玉米育種課題組以玉米自交系鄭58 為背景構(gòu)建的CMS-S 型不育系鄭58cms-Q1261 和保持系鄭58。

試驗(yàn)材料于2021 年5 月播種于河北農(nóng)業(yè)大學(xué)農(nóng)業(yè)科技園，播種后50 d 起，取不育系和保持系的花藥，利用涂片法鏡檢，確認(rèn)小孢子的發(fā)育時(shí)期。待發(fā)育至單核晚期進(jìn)行取樣，每個(gè)樣品由多株花藥組成，3次重復(fù)。保持系鄭58 和不育系鄭58cms-Q1261 分別命名為A 和B，3 次重復(fù)分別用A1、A2、A3 和B1、B2、B3 表示，液氮速凍后-80 ℃保存?zhèn)溆谩?/p>

1.2 掃描電鏡觀察

取單核晚期的花藥置于2.5%的戊二醛固定液中固定。分別用50%、75%、85%、90%、95% 乙醇逐級脫水，每個(gè)梯度30 min，最后用無水酒精脫水重復(fù)2～3 次。在預(yù)冷的醋酸異戊酯置換20 min，臨界點(diǎn)干燥。將干燥樣品固定在樣品臺(tái)上，用離子濺射法對材料進(jìn)行噴金鍍膜。最后用掃描電鏡（日立SU8100）進(jìn)行觀察。

1.3 花藥RNA 提取

使用植物總RNA 提取試劑盒（華越洋）提取花藥樣品的總RNA，使用NanoDrop 2000 (Thermo Fisher Scientific,Wilmington,DE) 對RNA 純度和濃度進(jìn)行檢測。利用1%的瓊脂糖凝膠電泳（SDSPAGE）和Agilent 生物分析儀2100 系統(tǒng)的RNA Nano 6000 檢測試劑盒對RNA 完整性進(jìn)行評估。

1.4 轉(zhuǎn)錄組測序及DEGs 分析

使用Illumina 測序平臺(tái)對樣品測序并分析，原始數(shù)據(jù)經(jīng)過處理后，得到clean reads 用于后續(xù)分析。利用HiSat2 軟件將過濾后的數(shù)據(jù)比對到玉米參考基因組（ZmB73_RefGen_v4），利用stringtie 軟件對轉(zhuǎn)錄本拼接并估算基因表達(dá)量。對樣本FPKM 值進(jìn)行統(tǒng)計(jì)，利用R 軟件包DESeq2 對基因表達(dá)值進(jìn)行DEGs 分析。以基因相對表達(dá)量的差異倍數(shù)|log2(fold change)| ＞1.5 且錯(cuò)誤發(fā)現(xiàn)率（False discovery rate，F(xiàn)DR）＜ 0.05 為DEGs 篩選的閾值標(biāo)準(zhǔn)。

1.5 qRT-PCR 驗(yàn)證

選取12 個(gè)DEGs 進(jìn)行實(shí)時(shí)定量PCR 驗(yàn)證。采用TRIZOL 法提取總RNA，合成相應(yīng)的cDNA。利用Premier5.0 software 設(shè)計(jì)熒光實(shí)時(shí)定量PCR 引物，引物序列如表1 所示。以玉米肌動(dòng)蛋白基因作為內(nèi)參基因，計(jì)算相對表達(dá)量（2-ΔΔCt法），保持系A(chǔ)中各目標(biāo)基因中的表達(dá)量為1，不育系B 目標(biāo)基因中的表達(dá)量與之進(jìn)行比較。

表1 qRT-PCR 中所用差異表達(dá)基因及其引物序列Table 1 Lists the DEGs and their primer sequences used in qRT-PCR

2 結(jié)果與分析

2.1 雄穗和花藥特征

觀察散粉期保持系A(chǔ) 和不育系B 的雄穗和花藥，如圖1 所示，發(fā)現(xiàn)保持系雄穗發(fā)達(dá)，穎殼開裂，花藥正常散粉，花粉量大；不育系雄穗外形和保持系無明顯差異，但穎殼不開裂，花藥不外露。進(jìn)一步對保持系和不育系的小花進(jìn)行觀察，保持系花藥飽滿，而不育系花藥干癟。

圖1 鄭58 和鄭58cms-Q1261 的雄穗和花藥Fig. 1 Tassel and anthers of Zheng 58 and Zheng 58cms Q1261

利用掃描電鏡（SEM）對保持系A(chǔ) 和不育系B單核晚期的花藥進(jìn)行觀察（圖2）?？梢姳３窒祷ㄋ庯枬M（圖A-1），花藥外壁的表皮細(xì)胞排列相對疏松（圖A-2），花藥內(nèi)壁表面纖維層存在突起，其表面存在烏氏小體（圖A-3），藥室內(nèi)有大量圓球狀飽滿的花粉粒（圖A-4，A-5）。不育系花藥干癟（圖B-1），花藥外壁表皮細(xì)胞排列相對緊密（圖B-2），內(nèi)壁較為平滑，其表面烏氏小體顆粒數(shù)量、大小、分布與保持系相似（圖B-3），藥室內(nèi)花粉粒皺縮，形狀不規(guī)則（圖B-4，B-5）。

圖2 鄭58 和鄭58cms-Q1261 花藥的掃描電鏡觀察Fig. 2 Sem observation of anther of Zheng58 and Zheng58 cms-Q1261

圖3 差異表達(dá)基因火山圖Fig. 3 Volcano map of DEGs

2.2 轉(zhuǎn)錄組學(xué)分析

2.2.1 測序質(zhì)量分析 利用Illumina HiSeq 2500 測序平臺(tái)對保持系A(chǔ) 和不育系B 單核晚期的花藥3 次重復(fù)樣品進(jìn)行轉(zhuǎn)錄組測序。所測數(shù)據(jù)經(jīng)過嚴(yán)格的質(zhì)控后，共得到41.67 Gb 的序列數(shù)據(jù)，各樣品Clean Data 均達(dá)到6.35 Gb。如表2 所示，獲得序列的GC含量在55.75%～56.56%，各樣品Q20 堿基百分比均不低于98.00%，Q30 堿基百分比在94.32%以上。表明樣品測序質(zhì)量較高，能夠用于后續(xù)的分析。

表2 轉(zhuǎn)錄組測序質(zhì)量分析Table 2 Transcriptome sequencing quality analysis

2.2.2 差異表達(dá)基因分析 基于基因的表達(dá)量，以|log2(fold change)| ＞1.5 且FDR＜ 0.05 為標(biāo)準(zhǔn)篩選差異表達(dá)基因。通過對保持系和不育系單核晚期花藥的轉(zhuǎn)錄本的比較分析，共鑒定到14 565 個(gè)DEGs，其中在不育系鄭58 cms-Q1261 中7 551 個(gè)基因上調(diào)表達(dá)，7 014 個(gè)基因下調(diào)表達(dá)（圖 3）。

2.2.3 差異表達(dá)基因GO 分析 對上述篩選到的DEGs 進(jìn)行GO 富集分析，共有11 232 個(gè)DEGs 注釋到生物學(xué)過程、細(xì)胞組分和分子功能3 個(gè)部分(圖4)。有8 440 個(gè)DEGs 注釋到生物學(xué)過程，涉及生物過程中的22 個(gè)功能條目，其中代謝進(jìn)程（上調(diào)DEGs 有3 079 個(gè)，下調(diào)有DEDGs 2 634 個(gè)）；細(xì)胞進(jìn)程（上調(diào)DEGs 有2 944 個(gè)，下調(diào)有DEDGs 2 643 個(gè)）和單組織過程（上調(diào)DEGs 有2 350 個(gè)，下調(diào)有DEDGs 1 930 個(gè)）等二級功能的基因富集比例較大。有8 539 個(gè)DEGs 注釋到細(xì)胞組分，涉及16 個(gè)功能條目，細(xì)胞（上調(diào)DEGs 有3 845 個(gè)，下調(diào)有DEDGs 3 499 個(gè)）；細(xì)胞組成（上調(diào)DEGs 有3 850 個(gè)，下調(diào)有DEDGs 3 504 個(gè)）和細(xì)胞器官（上調(diào)DEGs 有3 073 個(gè)，下調(diào)有DEDGs 2 995 個(gè)）這3 個(gè)二級功能占比較大。有8 568 個(gè)DEGs 注釋到分子功能的15 個(gè)條目，其中結(jié)合功能（上調(diào)DEGs 有3 292 個(gè)，下調(diào)有DEDGs 2 364 個(gè)）和催化活性（上調(diào)DEGs 有2 898 個(gè)，下調(diào)有DEDGs 2 272 個(gè)）所占比例較高。

2.2.4 差異表達(dá)基因KEGG 分析 為了進(jìn)一步了解DEGs 的生物學(xué)功能，將鑒定到的DEGs 進(jìn)行KEGG 分析。結(jié)果發(fā)現(xiàn)，DEGs 共富集到127 條代謝通路，表明玉米CMS-S 型不育系花粉敗育是多通路協(xié)調(diào)調(diào)控的復(fù)雜過程。分析顯著富集的前20 條通路（圖5），可見DEGs 主要富集到淀粉和蔗糖代謝、碳代謝、氧化磷酸化、糖酵解/糖異生、氨基糖和核苷酸糖代謝和植物-病原互作，而脂肪酸生物合成和糖胺多糖的降解富集較少。

圖5 差異表達(dá)基因 KEGG Pathway 富集氣泡圖Fig. 5 Bubble chart of KEGG Pathway enrichment of DEGs

2.3 差異表達(dá)基因qRT-PCR 驗(yàn)證

根據(jù)基因富集結(jié)果選取12 個(gè)DEGs 進(jìn)行qRTPCR 驗(yàn)證，詳細(xì)信息見表3，其中包括與戊糖和葡萄糖醛酸的相互轉(zhuǎn)化(ko00040)相關(guān)的2 個(gè)果膠裂解酶超家族蛋白；與淀粉和蔗糖代謝(ko00500)相關(guān)的5 個(gè)果膠裂解酶類；與氨基糖和核苷酸糖代謝(ko00520)相關(guān)的2 個(gè)bHLH 轉(zhuǎn)錄因子；與內(nèi)質(zhì)網(wǎng)蛋白加工(ko04141)相關(guān)的3 個(gè)熱休克蛋白。qRTPCR 數(shù)據(jù)通過2-ΔΔCT換算為相對表達(dá)量（Relative expression），并進(jìn)行獨(dú)立t檢驗(yàn)。圖6 結(jié)果顯示，這些 DEGs 在不育系B 和保持系A(chǔ) 中表達(dá)量差異極顯著，qRT-PCR 驗(yàn)證結(jié)果與 RNA-Seq 的差異表達(dá)趨勢一致，說明本研究的 RNA-Seq 測序結(jié)果可靠，推測這些基因可能在花粉發(fā)育過程中發(fā)揮重要的調(diào)控作用。

圖6 12 個(gè)差異基因的RT-qPCR 驗(yàn)證Fig. 6 RT-qPCR validation of 12 DEGs

3 討論

前人關(guān)于不同作物雄性不育轉(zhuǎn)錄組的研究發(fā)現(xiàn)不育過程常常伴有大量基因參與，不同作物的DEGs富集到的主要基因功能和代謝途徑不盡相同。魏霽桐等［10］對小麥生理型雄性不育系PHYMS-1376花藥5 個(gè)時(shí)期的轉(zhuǎn)錄組分析，顯示有36 058 個(gè)DEGs，主要富集在淀粉和蔗糖代謝、苯丙素生物合成和植物激素信號(hào)傳導(dǎo)等通路。薛亞東等［11］在玉米CMS-C 型“三系”材料減數(shù)分裂前期Ⅰ、中期Ⅰ及四分體時(shí)期花藥的轉(zhuǎn)錄組分析中，篩選出DEGs 5 676 個(gè)，主要富集在氧化磷酸化、碳代謝及糖酵解等能量代謝相關(guān)的途徑。Wei 等［12］關(guān)于棉花雄性不育的轉(zhuǎn)錄組研究結(jié)果顯示，在減數(shù)分裂和四分體階段分別有9 595 個(gè)、10 407 個(gè)DEGs，主要富集在蔗糖和淀粉代謝、激素合成和糖酵解磷酸戊糖途徑合成等通路。本研究對CMS-S 型不育系鄭58cms-Q1261 和保持系鄭58 單核晚期花藥進(jìn)行RNA-seq 測序，結(jié)果顯示有14 565 個(gè)DEGs，顯著富集在淀粉和蔗糖代謝、碳代謝、氧化磷酸化、糖酵解/糖異生、氨基糖和核苷酸糖代謝等途徑中。Xing 等［13］在油菜雄性不育的轉(zhuǎn)錄組分析中，發(fā)現(xiàn)DEGs 主要參與淀粉和蔗糖代謝、苯丙烷生物合成。Lui 等［14］關(guān)于煙草雄性不育研究中，DEGs 主要富集在內(nèi)質(zhì)網(wǎng)蛋白加工和氧化磷酸化。高煥庭等［15］在小麥雄性不育研究中，KEGG 代謝途徑主要有糖酵解途徑和三羧酸循環(huán)。與本研究中富集到的代謝途徑部分一致，表明這些代謝通路可能與花粉敗育相關(guān)。進(jìn)一步說明不同物種雄性不育均受眾多基因影響，是一個(gè)多通路協(xié)同調(diào)控過程，育性調(diào)控機(jī)理復(fù)雜。

淀粉和蔗糖代謝為小孢子發(fā)育提供營養(yǎng)，糖類、淀粉等物質(zhì)的缺失可能造成代謝紊亂，影響小孢子的形成，導(dǎo)致植物雄性不育［16］。Su 等［17］采用RNA-seq 技術(shù)，在玉米CMS-S 型雄性不育中鑒定出2 108 個(gè)DEGs，這些DEGs 在能量產(chǎn)生和轉(zhuǎn)化、碳水化合物代謝和信號(hào)轉(zhuǎn)導(dǎo)等生物過程中發(fā)揮作用，進(jìn)一步研究了淀粉和蔗糖代謝的途徑，找到一些與蔗糖生物合成過程有關(guān)的DEGs。果膠是花粉細(xì)胞壁的主要組分，在細(xì)胞壁結(jié)構(gòu)、信號(hào)轉(zhuǎn)導(dǎo)和植物生長發(fā)育中具有重要作用［18］，果膠降解是淀粉與蔗糖的代謝的一個(gè)重要分支。本研究選取與果膠酶相關(guān)的淀粉蔗糖代謝通路進(jìn)行分析，發(fā)現(xiàn)編碼果膠降解相關(guān)的一系列酶（包括2 個(gè)果膠裂解酶超家族蛋 白Zm00001d046389、Zm00001d053597 和5 個(gè)果膠裂解酶類Zm00001d015820、Zm00001d036820、Zm00001d015821、Zm00001d015825、Zm00001d036824）的DEGs 均顯著下調(diào)，推測在不育系中由于果膠裂解酶基因的缺失導(dǎo)致果膠異常積累，影響花粉的正常成熟進(jìn)而導(dǎo)致雄性不育。與Li 等［19］人在小麥花藥和花粉發(fā)育調(diào)控相關(guān)基因的研究中，果膠降解相關(guān)基因的下調(diào)導(dǎo)致花藥不開裂的結(jié)果一致。

轉(zhuǎn)錄因子對植物育性具有一定的調(diào)控作用，在擬南芥［20］、水稻［21］和玉米［22］等作物中能夠控制花藥分裂和雄性育性。其中bHLH 轉(zhuǎn)錄因子在所有轉(zhuǎn)錄因子家族中數(shù)量最多，已有研究表明，bHLH轉(zhuǎn)錄因子與花藥絨氈層和小孢子發(fā)育相關(guān)，其功能異常時(shí)，絨氈層細(xì)胞提前降解，花粉形成受阻導(dǎo)致小孢子發(fā)育受阻，最終表現(xiàn)敗育［23］。Jung 等［24］在水稻中鑒定了1 個(gè)與絨氈層發(fā)育相關(guān)的bHLH 轉(zhuǎn)錄因子UDT1，其逆轉(zhuǎn)錄轉(zhuǎn)座子Tos17 插入形成的不育植株表現(xiàn)雄蕊周緣細(xì)胞分化異常，絨氈層液泡化，花粉粒無法形成。Dukowic 等［25］比較擬南芥和玉米在減數(shù)分裂時(shí)期的轉(zhuǎn)錄組發(fā)現(xiàn)，14 個(gè)擬南芥bHLH基因和3 個(gè)玉米bHLH基因上調(diào)表達(dá)，并且發(fā)現(xiàn)玉米和擬南芥有對應(yīng)的同源基因。Nakata 等［26］在擬南芥中鑒定出參與植物激素信號(hào)傳導(dǎo)的3 個(gè)bHLH 轉(zhuǎn)錄因子（JAM1、JAM2 和JAM3），發(fā)現(xiàn)它們在雄蕊發(fā)育中作為負(fù)調(diào)控因子發(fā)揮作用。本研究也鑒定到2 個(gè)bHLH 轉(zhuǎn)錄因子（Zm00001d021019、Zm00001d038357）在單核晚期的不育系中上調(diào)表達(dá)。

熱休克蛋白是植物在逆境條件下產(chǎn)生的應(yīng)激蛋白，能夠平衡植物體內(nèi)的代謝，在生長發(fā)育中具有重要的作用［27］。有研究表明，在沒有逆境條件下，一些熱休克蛋白與配子體的發(fā)育相關(guān)，TMS1是擬南芥熱敏雄性不育基因，敲除TMS1 基因后發(fā)現(xiàn)在30℃條件下花粉管生長受阻，導(dǎo)致雄性生育能力顯著降低［28］。Mei 等［29］在蘿卜細(xì)胞質(zhì)雄性不育研究中，鑒定出5 個(gè)熱休克轉(zhuǎn)錄因子(HSFs)和19 個(gè)熱休克蛋白(HSPs) 基因在蘿卜CMS 系中下調(diào)表達(dá)，對小孢子發(fā)育產(chǎn)生影響。本研究中有3 個(gè)熱休克蛋白（Zm00001d012420、Zm00001d028557、Zm00001d020898）在不育系中顯著高表達(dá)，推測是其他代謝途徑紊亂導(dǎo)致熱休克蛋白高度表達(dá)，影響蛋白質(zhì)加工或降解，在一定程度上抑制其他種類的蛋白質(zhì)的合成，從而影響小孢子發(fā)育。

本研究對不育系鄭58cms-Q1261 和保持系鄭58 單核晚期花藥進(jìn)行轉(zhuǎn)錄組測序，篩選出14565 個(gè)DEGs，GO 和KEGG 分析DEGs 主要富集在淀粉和蔗糖代謝、碳代謝和氧化磷酸化等代謝通路，發(fā)現(xiàn)果膠裂解酶類、熱休克蛋白、bHLH 轉(zhuǎn)錄因子等與雄性不育有關(guān)，為進(jìn)一步深入解析玉米S 型細(xì)胞質(zhì)雄性不育的分子機(jī)制提供了理論基礎(chǔ)。