玉米CMS-S 型不育系花藥的轉錄組分析

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

（河北農業大學 農學院/國家玉米改良中心河北分中心/河北省作物種質資源實驗室，河北 保定 071001）

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

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

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

1 材料與方法

1.1 試驗材料

選用河北農業大學玉米育種課題組以玉米自交系鄭58 為背景構建的CMS-S 型不育系鄭58cms-Q1261 和保持系鄭58。

試驗材料于2021 年5 月播種于河北農業大學農業科技園，播種后50 d 起，取不育系和保持系的花藥，利用涂片法鏡檢，確認小孢子的發育時期。待發育至單核晚期進行取樣，每個樣品由多株花藥組成，3次重復。保持系鄭58 和不育系鄭58cms-Q1261 分別命名為A 和B，3 次重復分別用A1、A2、A3 和B1、B2、B3 表示，液氮速凍后-80 ℃保存備用。

1.2 掃描電鏡觀察

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

1.3 花藥RNA 提取

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

1.4 轉錄組測序及DEGs 分析

使用Illumina 測序平臺對樣品測序并分析，原始數據經過處理后，得到clean reads 用于后續分析。利用HiSat2 軟件將過濾后的數據比對到玉米參考基因組（ZmB73_RefGen_v4），利用stringtie 軟件對轉錄本拼接并估算基因表達量。對樣本FPKM 值進行統計，利用R 軟件包DESeq2 對基因表達值進行DEGs 分析。以基因相對表達量的差異倍數|log2(fold change)| ＞1.5 且錯誤發現率（False discovery rate，FDR）＜ 0.05 為DEGs 篩選的閾值標準。

1.5 qRT-PCR 驗證

選取12 個DEGs 進行實時定量PCR 驗證。采用TRIZOL 法提取總RNA，合成相應的cDNA。利用Premier5.0 software 設計熒光實時定量PCR 引物，引物序列如表1 所示。以玉米肌動蛋白基因作為內參基因，計算相對表達量（2-ΔΔCt法），保持系A中各目標基因中的表達量為1，不育系B 目標基因中的表達量與之進行比較。

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

2 結果與分析

2.1 雄穗和花藥特征

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

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

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

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

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

2.2 轉錄組學分析

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

表2 轉錄組測序質量分析Table 2 Transcriptome sequencing quality analysis

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

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

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

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

2.3 差異表達基因qRT-PCR 驗證

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

圖6 12 個差異基因的RT-qPCR 驗證Fig. 6 RT-qPCR validation of 12 DEGs

3 討論

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

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

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

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

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