基于轉錄組測序揭示玉米抗倒伏相關基因和代謝通路

劉松濤 田再民 劉子剛 高志佳 張靜 賀東剛 黃智鴻 蘭鑫

（1 河北北方學院/河北省農產品食品質量安全分析檢測重點實驗室，075000，河北張家口；2 河北巡天農業科技有限公司，075000，河北張家口；3 河北兆育種業集團有限公司，050000，河北石家莊）

玉米是我國重要的糧食作物之一，具有高產、適應性廣和營養物質豐富等特點，因此玉米高產穩產在農業生產中占有重要地位[1]。玉米生長發育過程中經常發生倒伏現象，造成產量大幅下降，據統計[2-3]，倒伏可造成玉米減產15%～20%，導致我國每年玉米產量損失約100 萬t。玉米倒伏主要包括根倒、莖倒和莖折斷3 種形式。導致玉米倒伏的因素有氣候因素、種植密度、土壤質地以及玉米株高、莖稈柔韌性和根系發達程度等自身遺傳因素[4-6]。

莖稈力學特征是影響玉米莖稈倒伏的重要因素之一，代表玉米抗倒伏能力的莖稈力學特征有莖稈拉力、抗推力、穿刺強度和壓碎強度，且與玉米抗倒伏能力呈顯著正相關[7]。莖稈顯微結構是影響玉米莖稈抗倒伏性能的另一個重要因素，維管束數目越多、細胞表皮越厚的玉米莖稈強度越好，抗倒性越強。馮素偉等[8]對莖稈顯微結構進行分析發現，莖稈顯微結構中大維管束數量越多，玉米的抗倒能力越強。前人對玉米倒伏的研究主要集中在莖稈力學特征、形態學特征以及相應的生理生化指標方面，在分子水平對玉米抗倒伏相關基因的挖掘的研究較少。因此，本研究以高抗倒的京農科728、中抗倒的金農738 和低抗倒的先玉335 為材料，在玉米抽雄期取第3 節間莖稈進行轉錄學測序，挖掘與玉米抗倒性相關的基因，為培育高產和高抗倒伏的玉米品種提供理論依據。

1 材料與方法

1.1 試驗設計

試驗于2019 年在河北省張家口市沙嶺子鎮農業科學院（115°05′E，40°6′N）進行，土壤質地為壤質土。供試玉米品種為京農科728（JNK728，高抗倒性，H）、金農738（JN738，中抗倒性，M）和先玉335（XY335，低抗倒性，L），試驗采取完全隨機區組設計，3 次重復，行距55cm，株距30cm，種植密度67 500 株/hm2，每個品種12 行，行長10m。在玉米抽雄期，每個小區取3 株長勢一致的玉米，取第3 節間莖稈液氮速凍后于-80℃保存待測。

1.2 玉米莖稈樣本總RNA 的提取、文庫構建及Illumina 測序

使用康為全能型植物R N A 提取試劑盒（OminiPlant RNA Kit）提取9 個玉米莖稈樣本的總RNA。使用Qubit 2.0 熒光計測定RNA 的濃度，用Agilent 2100 生物分析儀和1%的瓊脂糖凝膠電泳檢測RNA 的完整性。使用Illumina 公司的TruSeqTM RNA 試劑盒（San Diego）構建測序所需cDNA 文庫，用Agilent 2100 生物分析儀檢測文庫質量。由武漢邁特維爾生物科技公司完成建庫測序。

1.3 數據處理、比對及表達量評估

使用Illumina HiSeq 高通量測序平臺對cDNA文庫進行測序，對生成的原始數據進行處理，除去低質量讀段、含接頭或引物的讀段，利用HISAT2將高質量讀段與參考基因組（B73 RefGen_v4）進行序列比對，獲取在參考基因組或基因上的位置信息以及測序樣品特有的序列特征信息。所有下游分析均基于與參考數據完美匹配或只有一處不匹配的讀段。采用FPKM（fragments per kilobase per million mapped reads）計算轉錄本或基因的表達水平。

1.4 差異表達基因的鑒定與功能注釋

為了鑒定與玉米莖折相關基因，用DESeq2 分別對3 個品種組間的差異表達進行分析，基因表達差異倍數在2 倍以上、且FDR＜0.05 確定為差異表達基因（differentially expressed genes，DEGs）。為了進一步鑒定DEGs 的功能，使用Blast2GO（ https://www.blast2go.com/ ） 和 KOBAS 2.1.1(http://kobas.cbi.pku.edu.cn/download.php)進行GO和代謝通路富集分析（FDR ＜0.05）；使用Venny2.1.0 構建Venn 圖。

1.5 實時熒光定量PCR（qRT-PCR）驗證

為了驗證RNA 測序結果的準確性，選擇10個DEGs 進行qRT-PCR 驗證。采用Primer 5 對DEGs 進行熒光定量引物設計。內參基因為玉米GAPDH，引物信息見表1。使用Bio-Rad iQ5 熒光定量PCR 儀（Bio-RAD，美國），選用2×Fast Super EvaGreen ?qPCR Mastermix（EverbrightInc.，美國）進行熒光定量分析。反應體系為cDNA 模板1μL、前后引物各1μL、ddH2O 7μL、2×AugeGreen TMMaster Mix 10μL。用2-ΔΔCT法[9]計算基因的相對表達量。

表1 熒光定量引物序列Table 1 The primer sequence of RT-qPCR

2 結果與分析

2.1 轉錄組測序結果統計與質量評估

用HiSeq6000 平臺對9 個玉米莖稈樣本進行雙端測序后共獲得59.26Gb 高質量讀段，與玉米參考基因組比對后發現，84.84%～87.28%的高質量讀段比對到參考基因的唯一位置，3.04%～3.58%的有效數據比對到參考基因的多個位置，Q30 和GC 含量分別高于93.90%和53.16%（表2）。其次，用主成分分析（PCA）評估樣本之間相關性，結果（圖1）表明，生物學重復之間高度相關，每個材料間明確分離。以上結果證實轉錄組測序數據是可靠的，可用于進一步分析。

圖1 轉錄組測序樣本PCA 分析Fig.1 Principal component analysis of samples used for transcriptome sequencing

表2 轉錄組測序數據統計Table 2 The summary of RNA-Seq data

2.2 差異表達基因鑒定

為了挖掘玉米抗倒伏相關基因，對3 個抗倒性不同的玉米材料進行品種間基因差異表達分析。由表3 可知，先玉335 與京農科728 比較（L-vs-H）鑒定到差異基因最多，其中4411 個DEGs 上調表達，3368 個DEGs 下調表達。先玉335 與京農738 比較（L-vs-M）鑒定到5373 個DEGs，包括2896 個DEGs 下調表達，2477 個DEGs上調表達。京農738 與京農科728 比較（M-vs-H）鑒定到4905 個DEGs，2831 個DEGs 下調表達，2074 個DEGs 上調表達。3 個比較組共鑒定到10 093 個DEGs，其中2095、769、606 個DEGs分別在L-vs-H、L-vs-M、M-vs-H 分組特異表達，1341 個DEGs 在3 個分組中共同表達（圖2）。

圖2 分組鑒定到的DEGs Venn 圖分析Fig.2 Venn diagram analysis of DEGs observed in comparison groups

表3 差異表達基因統計Table 3 Statistics of differentially expressed genes

2.3 DEGs GO 功能富集分析

為了探究抗倒性不同玉米品種差異表達基因的生物學功能，對鑒定到的DEGs 進行GO 功能富集分析。由圖3 可知，生物學過程分類中，L-vs-H分組顯著富集到植物型細胞壁的生物發生、光合作用、木質素合成過程，L-vs-M 分組顯著富集到植物型次生細胞壁的生物發生，M-vs-H 分組顯著富集到次生代謝過程、細胞壁發生、苯丙素的代謝過程。分子功能分類中，L-vs-H、M-vs-H 分組顯著富集到血紅素結合、單加氧酶活性、碳水化合物結合；L-vs-M 分組顯著富集到葉綠素結合、四吡咯結合、糖基轉移酶活性。細胞成分分類中，光合體系、光合體系Ⅰ、光合體系Ⅱ在L-vs-H、L-vs-M 分組中顯著富集，非原質體在M-vs-H 分組中顯著富集。因此，3 個玉米品種抗倒性不同可能是因為富集到相同GO 條目的DEGs 數量不同。

圖3 DEGs 的GO 富集分析Fig.3 Gene ontology(GO)enrichment analysis of DEGs

2.4 DEGs 代謝通路富集分析

為了近一步了解DEGs 的功能，對3 個分組鑒定到的DEGs 進行KEGG 代謝通路富集分析。結果（圖4）表明，L-vs-H 和M-vs-H 分組均富集到6 條代謝通路，L-vs-M 分組富集到5 條代謝通路。其中，苯丙素生物合成、次生代謝產物的生物合成、類黃酮生物合成在L-vs-H 和M-vs-H 分組顯著富集。L-vs-M 分組中，光合作用―天線蛋白、植物―病原體相互作用、次生代謝產物的生物合成顯著富集。

圖4 DEGs 代謝通路富集分析Fig.4 Metabolic pathway enrichment analysis of DEGs

2.5 3 個玉米品種莖稈的顯微結構差異

3 個玉米品種的顯微鏡可視范圍內莖稈維管束數目、單個維管束面積、莖稈表皮細胞厚度如圖5所示。先玉335、金農738 和京農科728 顯微鏡可視范圍內莖稈維管束數目分別為32、39、44 個（圖5a 和圖6a）；先玉335 單個維管束面積最小，京農科728 單個維管束面積最大（圖5b 和圖6b）；莖稈表皮細胞厚度的順序依次為京農科728＞金農738＞先玉335（圖5c 和圖6c）。由此可知，京農科728 的高抗倒性可能與較少的維管束數目、較大的單個維管束面積以及較厚的莖稈表皮細胞厚度有關。

圖5 3 個玉米品種莖稈的顯微結構Fig.5 Microstructure of the stems of the three maize varieties

圖6 3 個玉米品種莖稈的顯微結構Fig.6 Microstructures of the stalk of the three maize varieties

2.6 DEGs 的熒光定量驗證

為了驗證轉錄組測序的準確性，通過qRT-PCR檢測了10 個已報道與莖稈抗倒相關或注釋到與莖稈抗倒相關途徑的DEGs 表達量，結果（圖7）表明，所有10 個基因的表達水平與RNA-seq 測序的結果一致，且相關系數達到0.9171，因此qRT-PCR的分析結果證實了轉錄組測序結果的可靠性。

圖7 10 個DEGs 的RNA-seq 測序結果的qRT-PCR 驗證Fig.7 The qRT-PCR validation of the RNA-seq data for the ten DEGs

3 討論

玉米在我國糧食安全戰略中占有重要地位，隨著玉米種植密度增加，莖稈增高，莖稈倒伏問題日益嚴重，成為限制玉米產量和品質的因素之一。作物倒伏指由各種外界因素引發的植株莖稈由自然直立狀態到永久錯位的現象，是農作物生產中存在的一個普遍性問題。導致作物倒伏的原因主要有遺傳、栽培管理和自然環境3 個方面，其中由作物基因型決定的遺傳因素是內因，也是決定作物自身是否抗倒的最直接、最根本的因素[10]。植物生長處于風、雨等自然條件的影響下，當作物的莖稈彈性不足以使彎曲的植株恢復時容易發生大面積倒伏。玉米莖稈單位面積內維管束數目和表皮細胞厚度決定了玉米莖稈的抗壓強度，本研究通過比較分析不同玉米品種的維管束數目和表皮細胞厚度，結果表明維管束數目越多，表皮細胞越厚，玉米莖稈抗折強度越高。楊碩等[11]對玉米莖稈顯微結構研究表明，莖稈維管束數目、大維管束數目和小維管束數目在不同環境下均與莖稈抗倒伏能力呈正相關。同樣，對水稻的研究[12]發現，抗倒伏能力越強的水稻，其莖壁厚，大、小維管束數目多。

莖稈抗倒伏是由多個基因控制的復雜的數量性狀，目前對于玉米莖稈抗倒伏相關基因與代謝通路的研究較少。因此，本研究對3 個抗倒性不同玉米品種進行轉錄學測序，并對DEGs 進行了GO 和KEGG 富集分析，報道了玉米抗倒伏的關鍵差異表達基因和代謝通路。木質素是一類重要的大分子有機物質，由類苯丙酸途徑的單體衍生形成，在維持細胞壁結構完整性方面起著關鍵性作用，莖稈的木質素含量越高作物的抗倒性越強[13]。木質素合成途徑中，CoA 連接酶（4CL）引導光合代謝產物向木質素合成途徑的流動，處于代謝分支點的位置[14]。研究[10]表明，當4CL 活性比原來提升1 倍時，木質素含量也會相應提高1/4。本研究中苯丙素生物合成、木聚糖生物合成過程、細胞壁生物發生和植物型次生細胞壁的生物發生在3 個分組中顯著富集，對大豆莖倒伏的研究[10]表明DEGs 也與這些功能相關。已報道的與莖倒伏相關的苯丙烷類代謝途徑中的關鍵酶苯丙氨酸解氨酶、肉桂酸4-羥化酶和連接植物苯丙烷復合途徑和木質素特異生物合成途徑的4CL 顯著差異表達。推測這些顯著富集的通路和差異基因可能與玉米抗倒性相關，其功能還需進一步驗證。

4 結論

對3 個抗倒性不同的玉米品種抽雄期莖稈進行轉錄組學分析，3 個比較組共鑒定到10 093 個DEGs，其中2095、769、606 個DEGs 分別在L-vs-H、L-vs-M 和M-vs-H 分組特異表達，1341 個DEGs在3 個分組中共同表達。GO 富集分析表明，3 個玉米品種抗倒性不同可能是因為富集到相同GO 條目的DEGs 數量不同。代謝通路富集分析表明，L-vs-H、M-vs-H 分組顯著富集到苯丙素生物合成、次生代謝產物的生物合成、類黃酮生物合成途徑，而L-vs-M 分組顯著富集到光合作用―天線蛋白、植物―病原體相互作用途徑，注釋到這些途徑的基因為與玉米抗倒伏相關的關鍵候選基因。莖稈的顯微結構分析表明，京農科728 的抗倒性強可能與莖稈維管束數目多、單個維管束面積大、莖稈表皮細胞厚有關。此外，qRT-PCR 驗證了轉錄組測序的結果。本研究確定的關鍵基因和代謝途徑可作為未來定向克隆和下游分析研究的寶貴遺傳資源或選擇目標。