利用WGCNA鑒定玉米株高和穗位高基因共表達(dá)模塊

馬 娟 曹言勇 王利鋒 李晶晶 王 浩 范艷萍 李會(huì)勇

利用WGCNA鑒定玉米株高和穗位高基因共表達(dá)模塊

馬 娟 曹言勇 王利鋒 李晶晶 王 浩 范艷萍 李會(huì)勇*

河南省農(nóng)業(yè)科學(xué)院糧食作物研究所, 河南鄭州 450002

株高和穗位高是玉米株型的重要影響因子, 與產(chǎn)量性狀緊密相關(guān)。加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析(weighted gene co-expression network analysis, WGCNA)是探索基因網(wǎng)絡(luò)與特定性狀間關(guān)聯(lián)關(guān)系的重要方法, 為株高和穗位高相關(guān)基因的挖掘提供新途徑。本研究利用鄭58、掖478、昌7-2和黃早四及其組配的雜交種鄭單958、安玉5號(hào)、鄭58/黃早四和掖478/黃早四, 結(jié)合其在45,000株 hm–2和67,500株 hm–2條件下的轉(zhuǎn)錄組數(shù)據(jù), 采用WGCNA構(gòu)建了2種密度條件下的共表達(dá)網(wǎng)絡(luò), 分別得到24個(gè)和21個(gè)共表達(dá)模塊, 并鑒定到與株高和穗位高顯著且高度相關(guān)(相關(guān)系數(shù)的絕對(duì)值>0.50)的共表達(dá)模塊15個(gè), 其中兩性狀相同的模塊共6個(gè)?；蚬δ芨患治鼋Y(jié)果表明, 株高和穗位高共表達(dá)模塊主要參與生長(zhǎng)發(fā)育、光合作用、響應(yīng)光刺激、植物激素、碳水化合物合成/代謝等重要活動(dòng)。根據(jù)模塊內(nèi)基因的連接度, 發(fā)現(xiàn)乙烯響應(yīng)因子、硫胺素酶、磷酸甘油酸激酶、谷胱甘肽轉(zhuǎn)移酶和琥珀酸脫氫酶等是模塊內(nèi)的核心基因。通過構(gòu)建其局部網(wǎng)絡(luò), 發(fā)現(xiàn)與已報(bào)道株高基因、和以及C3H轉(zhuǎn)錄因子、C2C2-GATA轉(zhuǎn)錄因子和乙烯受體同源子等存在關(guān)聯(lián)關(guān)系。此外, 已報(bào)道株高基因和也存在于共表達(dá)模塊中。以這5個(gè)已報(bào)道株高基因?yàn)楹诵? 構(gòu)建其基因網(wǎng)路, 發(fā)現(xiàn)生長(zhǎng)素轉(zhuǎn)錄因子、、、光合系統(tǒng)II氧進(jìn)化多肽和光合系統(tǒng)I N亞基等與其存在關(guān)聯(lián)。15個(gè)共表達(dá)模塊和核心基因的挖掘以及基因生物學(xué)功能和互作網(wǎng)絡(luò)的解析有助于揭示玉米株高和穗位高的遺傳基礎(chǔ)。

玉米; 加權(quán)基因共表達(dá)網(wǎng)絡(luò); 轉(zhuǎn)錄組; 株高; 穗位高

株高和穗位高是玉米株型的重要構(gòu)成因子, 與產(chǎn)量密切相關(guān)。因此解析玉米株高和穗位高的遺傳機(jī)制, 對(duì)玉米株型育種和產(chǎn)量的提高具有重要意義。前人利用數(shù)量性狀定位和全基因組關(guān)聯(lián)分析方法對(duì)玉米株高和穗位高的遺傳構(gòu)成開展了大量研究, 挖掘了許多與之相關(guān)的QTL/QTN[1-6], 也克隆了株高相關(guān)基因如()、、、、、()等[7-11]。盡管株高和穗位高研究取得了重要進(jìn)展, 但從轉(zhuǎn)錄組水平解析其遺傳構(gòu)成的研究較少。

轉(zhuǎn)錄組研究是基因功能及結(jié)構(gòu)研究的基礎(chǔ)和出發(fā)點(diǎn)。加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析(weighted gene co-expression network analysis, WGCNA)是通過網(wǎng)絡(luò)研究基因功能的重要方法[12]。將轉(zhuǎn)錄組數(shù)據(jù)與WGCNA方法結(jié)合可以挖掘出與特定性狀相關(guān)的模塊和基因。目前, 前人已經(jīng)利用WGCNA方法研究了玉米籽粒大小和產(chǎn)量雜種優(yōu)勢(shì)等性狀和組織特異性模塊。Zhang等[13]發(fā)現(xiàn)13個(gè)與玉米粒重顯著相關(guān)的模塊, 其相關(guān)性為?0.87~0.89。Ma等[14]在67,500株 hm–2條件下挖掘出7個(gè)模塊與玉米產(chǎn)量等性狀雜種優(yōu)勢(shì)指數(shù)顯著相關(guān), 通過模塊功能注釋分析, 發(fā)現(xiàn)1461個(gè)基因可能與玉米雜種優(yōu)勢(shì)相關(guān)。Zhan等[15]發(fā)現(xiàn)胚乳隔間、糊粉層和胚乳等組織特異性模塊10個(gè)。楊宇昕等[16]對(duì)B73不同發(fā)育階段的轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行了WGCNA分析, 鑒定到14個(gè)組織特異性模塊, 其中8個(gè)組織特異性模塊內(nèi)基因富集到與開花調(diào)控相關(guān)的代謝通路。此外, WGCNA也可用來挖掘非生物脅迫如鎘脅迫[17]和淹水脅迫的相關(guān)基因[18]、氮響應(yīng)的長(zhǎng)鏈非編碼RNA[19]和肌醇磷酸鹽代謝涉及的新基因和調(diào)控機(jī)制[20]等。

本研究結(jié)合轉(zhuǎn)錄組測(cè)序(RNA-Seq)數(shù)據(jù)和WGCNA方法, 構(gòu)建玉米加權(quán)基因共表達(dá)網(wǎng)絡(luò)并劃分模塊, 挖掘株高和穗位高特異性模塊, 通過對(duì)特異性模塊分析, 挖掘株高和穗位高顯著相關(guān)核心基因及其互作網(wǎng)絡(luò), 為進(jìn)一步研究玉米株高和穗位高的分子機(jī)制提供新的思路。

1 材料與方法

1.1 材料和轉(zhuǎn)錄組測(cè)序數(shù)據(jù)

玉米自交系鄭58、掖478、昌7-2、黃早四及其組配的雜交種鄭單958 (鄭58/昌7-2)、安玉5號(hào)(掖478/昌7-2)、鄭58/黃早四和掖478/黃早四于2016年種植在河南省農(nóng)業(yè)科學(xué)院新鄉(xiāng)原陽試驗(yàn)田。采用隨機(jī)區(qū)組試驗(yàn)設(shè)計(jì), 3次重復(fù)。種植密度為45,000株 hm–2和67,500株 hm–2。授粉后15 d, 采集8個(gè)材料的穗位葉提取RNA, 用于轉(zhuǎn)錄組測(cè)序。采用Illumina X Ten對(duì)2種密度條件8個(gè)材料3個(gè)生物學(xué)重復(fù)進(jìn)行雙端測(cè)序。48個(gè)樣品的轉(zhuǎn)錄組測(cè)序數(shù)據(jù)見NCBI (National Center for Biotechnology Information)的Sequence Read Archive數(shù)據(jù)庫(登錄號(hào)為 SRP136913)。利用Tophat v2.0.10將clean reads比對(duì)到玉米B73參考基因組(ftp://ftp.ensemblgenomes. org/pub/release-29/plants/fasta/zea_mays/dna/Zea_mays.AGPv3.29.dna.toplevel.fa.gz)。利用Cufflinks v.2.1.1計(jì)算基因在兩種密度條件下8個(gè)材料3個(gè)生物學(xué)重復(fù)轉(zhuǎn)錄本的表達(dá)量。利用每千堿基外顯子百萬片段數(shù)(fragments per kilobase of exon per million fragments mapped, FPKM)值衡量基因的表達(dá)水平。在兩種密度條件下, 分別選取24個(gè)樣品FPKM值均大于1的基因, 用于構(gòu)建基因共表達(dá)網(wǎng)絡(luò)分析。授粉后45 d, 每個(gè)材料每個(gè)重復(fù)測(cè)量5個(gè)植株的株高和穗位高。利用R語言計(jì)算方差分析和Duncan’s多重比較。

1.2 加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析

利用R軟件WGCNA包構(gòu)建基因共表達(dá)網(wǎng)絡(luò)分析(https://horvath.genetics.ucla.edu/html/Coexpression Network/Rpackages/WGCNA/index.html)。為了確保符合無尺度網(wǎng)絡(luò)分布, WGCNA需要選擇合適的加權(quán)系數(shù)β。β利用WGCNA包中的pickSoftThreshold函數(shù)計(jì)算。首先設(shè)置β=1–20, 分別計(jì)算其對(duì)應(yīng)的相關(guān)系數(shù)和基因連接度均值。β的選取標(biāo)準(zhǔn)即滿足相關(guān)系數(shù)的平方接近0.8, 同時(shí)還需要保證一定的基因連接度。根據(jù)圖1結(jié)果, 在45,000株 hm–2和67,500株 hm–2條件下, β值均取9來構(gòu)建共表達(dá)網(wǎng)絡(luò)。采用動(dòng)態(tài)切割法(dynamic tree cut)識(shí)別共表達(dá)模塊。利用自動(dòng)網(wǎng)絡(luò)構(gòu)建函數(shù)blockwiseModules構(gòu)建網(wǎng)絡(luò), 其中最小模塊大小(minModuleSize)為50, mergeCutHeight=0.25合并相似性為0.75的模塊, 其他參數(shù)按照默認(rèn)設(shè)置。利用softConnectivity函數(shù)計(jì)算基因的連接度, 篩選出模塊內(nèi)連接度前10的基因, 定義這些基因?yàn)槟K的核心基因。根據(jù)拓?fù)渲丿B矩陣計(jì)算出模塊內(nèi)不同基因間的權(quán)重值。權(quán)重值越高表示基因間的關(guān)聯(lián)程度越高。根據(jù)權(quán)重值, 選取與核心基因和已報(bào)道株高基因存在互作網(wǎng)絡(luò)前10%的基因, 利用Cytoscape v3.6.1進(jìn)行可視化展示。

圖1 45,000株hm–2(A) 和67,500株 hm–2(B)條件下軟閾值β的確定

A和B: 左圖縱坐標(biāo)是無尺度網(wǎng)絡(luò)模型指數(shù); 右圖縱坐標(biāo)每一個(gè)軟閾值對(duì)應(yīng)的平均連接度; 橫坐標(biāo)均代表軟閾值β。

A and B: The ordinate represents the index of scale free network model in left figure. The ordinate represents the average link degree of each soft threshold in right figure. The abscissa represents the soft threshold β.

1.3 株高和穗位高顯著關(guān)聯(lián)模塊的鑒定和基因功能分析

為了找到與株高和穗位高相關(guān)性較高的基因, 本研究定義相關(guān)系數(shù)絕對(duì)值在0.5以上且顯著性達(dá)到0.05水平的模塊為顯著模塊。利用R語言的clusterProfiler包對(duì)顯著模塊內(nèi)基因進(jìn)行GO (gene ontology)分析。以FDR=0.05為閾值來篩選顯著生物學(xué)過程(biological process, BP)。

2 結(jié)果與分析

2.1 兩種密度條件下8個(gè)材料株高和穗位高的表現(xiàn)

4個(gè)自交系和4個(gè)雜交種在67,500株 hm–2的株高和穗位高均高于45,000株 hm–2(圖2)。只有掖478/黃早四的株高在兩種密度條件下存在顯著差異(雙尾測(cè)驗(yàn),=0.011)(圖2-A)。方差分析結(jié)果顯示同一密度條件下8個(gè)材料間存在顯著差異(附表1和附表2)。Duncan’s多重比較結(jié)果表明同一密度條件下, 4個(gè)雜交種的株高和穗位高均顯著高于4個(gè)自交系(圖2)。兩種密度條件下, 安玉5號(hào)的株高和穗位高均顯著高于鄭單958。2種密度條件下, 4個(gè)自交系的株高存在顯著差異, 以昌7-2的株高最高、黃早四次之、掖478和鄭58最低。兩種密度條件下, 昌7-2和黃早四的穗位高沒有顯著差異, 掖478和鄭58的穗位高也沒有顯著差異, 但昌7-2和黃早四的穗位高均顯著高于掖478和鄭58。

圖2 45,000株 hm–2和67,500株hm–2條件下8個(gè)材料株高和穗位高

圖中數(shù)字表示均值±標(biāo)準(zhǔn)誤。柱中字母為Duncan’s多重比較結(jié)果, 不同字母表示材料間在< 0.05水平差異顯著。*表示0.05顯著水平。

Number represents mean±SE. Bars with different letters are significantly different at< 0.05 as determined by Duncan’s multiple comparison.*denotes significant at the 0.05 probability level.

2.2 基因共表達(dá)網(wǎng)絡(luò)的構(gòu)建

經(jīng)過濾FPKM值<1的基因, 24個(gè)樣品在45,000株 hm–2和67,500株 hm–2密度下, 分別有14,251個(gè)和14,190個(gè)基因用于構(gòu)建加權(quán)基因共表達(dá)網(wǎng)絡(luò)。首先根據(jù)基因的表達(dá)量進(jìn)行聚類分析, 聚類程度高的基因?qū)⒈环值揭粋€(gè)模塊。然后利用動(dòng)態(tài)切割法識(shí)別共表達(dá)模塊。不同的模塊用不同的顏色表示。其中, Grey 模塊表示未分配到任何模塊的基因。低密度條件下共獲得24個(gè)共表達(dá)模塊(圖3-A), 不同模塊間包含的基因個(gè)數(shù)差異較大。其中, Turquoise模塊的基因個(gè)數(shù)最多, 為3777, Darkturquoise模塊的基因個(gè)數(shù)最少, 為57個(gè)(圖4)。高密度條件下共獲得21個(gè)模塊(圖3-B)。其中, Turquoise模塊的基因個(gè)數(shù)最多(3283), Royalblue基因個(gè)數(shù)最少(71)(圖5)。

2.3 株高和穗位顯著相關(guān)的共表達(dá)模塊及其功能注釋

在45,000株hm–2條件下, 24個(gè)模塊中與株高和穗位高相關(guān)性絕對(duì)值大于0.50, 且顯著的模塊分別有7個(gè)和4個(gè), 兩性狀共同的模塊有2個(gè)(圖4)。其中, Turquoise模塊與株高和穗位高均具有較高的相關(guān)性, 相關(guān)系數(shù)分別為?0.62 (=0.0013)和?0.80 (=2.33E-06)。Lightgreen模塊與株高和穗位高均表現(xiàn)出顯著負(fù)相關(guān), 相關(guān)系數(shù)介于?0.51~ ?0.64。Turquoise模塊富集到較多生物學(xué)過程, 共有1667個(gè)BP條目(附圖1)。其中215個(gè)條目與發(fā)育相關(guān), 121個(gè)條目與轉(zhuǎn)運(yùn)相關(guān), 67個(gè)條目與定位相關(guān)。Turquoise模塊內(nèi)基因還顯著富集在信號(hào)途徑(58)、信號(hào)傳導(dǎo)(29)、碳水化合物代謝合成(44)、泛素相關(guān)的過程(23)、光合作用和光響應(yīng)(24)以及激素相關(guān)的生物學(xué)過程(22)等。Lightgreen模塊顯著富集到255個(gè)BP條目, 其中43.92%為發(fā)育相關(guān)的過程(附圖2)。

圖3 45,000株hm–2(A)和67,500株hm–2(B)條件下基因聚類樹和模塊構(gòu)建

圖4 45,000株hm–2條件下株高(PH)和穗位高(EH)共表達(dá)模塊及其基因個(gè)數(shù)和相關(guān)系數(shù)

在67,500株hm–2條件下, 與株高和穗位顯著相關(guān)的共表達(dá)模塊分別有4個(gè)和6個(gè)(圖5)。其中, 與兩性狀均相關(guān)的基因模塊有4個(gè)。Midnightblue和Pink模塊與株高和穗位高均存在高度相關(guān)性, 相關(guān)系數(shù)介于?0.70~ ?0.80。Purple和Tan模塊與兩性狀也均有較高的相關(guān)性(= ?0.64 ~ ?0.72)。Midnightblue模塊顯著富集到362個(gè)BP條目, 主要包括發(fā)育相關(guān)過程(80)、碳水化合物代謝合成(18)、光合作用相關(guān)和光響應(yīng)(15)、信號(hào)傳導(dǎo)(15)、泛素相關(guān)的途徑(11)和植物激素(3)等(附圖3)。Pink模塊基因也主要顯著富集在發(fā)育相關(guān)(29)、定位(41)、光合作用(6)和光響應(yīng)(3)等的生物學(xué)過程(附圖4)。Purple模塊基因主要富集在生物脅迫相關(guān)的過程(附圖5)。Tan模塊顯著富集40個(gè)BP條目, 其中對(duì)光、熱等刺激的響應(yīng)有19條, 發(fā)育相關(guān)的條目有11條, 光合作用有2條(附圖6)。

圖5 67,500株hm–2條件下株高(PH)和穗位高(EH)共表達(dá)模塊及其基因個(gè)數(shù)和相關(guān)系數(shù)

2.4 顯著共表達(dá)模塊的核心基因和已報(bào)道株高基因的基因網(wǎng)絡(luò)

選取顯著模塊內(nèi)連接度最高的前10個(gè)基因作為該模塊的核心基因, 并對(duì)核心基因的互作網(wǎng)絡(luò)利用Cytoscape軟件進(jìn)行了可視化。低密度條件下, Turquoise模塊內(nèi)連接度最高的前10個(gè)基因中, 有注釋信息的基因?yàn)楹土虬匪孛?。是乙烯響?yīng)的轉(zhuǎn)錄因子, 可能是株高和穗位高重要的候選基因。由于與互作的基因較多, 共有2750個(gè), 圖6僅對(duì)權(quán)重值前150個(gè)有關(guān)聯(lián)的基因進(jìn)行可視化。與具有較高互作網(wǎng)絡(luò)關(guān)系的基因?yàn)镃3H轉(zhuǎn)錄因子、谷胱甘肽轉(zhuǎn)移酶、已報(bào)道株高基因[21]、乙烯受體同源子和C2C2-GATA轉(zhuǎn)錄因子等。已知株高基因[10]和[22]與也存在關(guān)聯(lián), 但其權(quán)重值不在前10%中。的關(guān)聯(lián)基因有、Alfin類轉(zhuǎn)錄因子、、肌動(dòng)蛋白相關(guān)蛋白、F-box/kelch 重復(fù)蛋白SKIP11、FAR1相關(guān)序列、順烏頭酸、二硫鍵異構(gòu)酶、核結(jié)合蛋白和亞油酸酯13S脂氧合酶等(圖6)。在Lightgreen模塊中, 核心基因有注釋信息的是磷酸甘油酸激酶, 其關(guān)聯(lián)的基因有SAR同源子和蛋白質(zhì)二硫鍵異構(gòu)酶等。67,500株 hm–2條件下, Pink模塊核心基因是和琥珀酸脫氫酶。圖6顯示這2個(gè)基因間存在關(guān)聯(lián)。與這2個(gè)基因均有網(wǎng)絡(luò)關(guān)系的基因?yàn)?、磷酸烯醇丙酮酸羧化酶、鋁誘導(dǎo)蛋白同源子和葉綠體外膜的轉(zhuǎn)位子等。

已報(bào)道株高相關(guān)基因如、和在45,000株hm–2條件下的Turquoise模塊中。以這3個(gè)基因?yàn)楹诵臉?gòu)建局部網(wǎng)絡(luò), 發(fā)現(xiàn)Turquoise模塊的核心基因與、和均存在互作(圖7)。與、、半胱氨酸蛋白酶、半胱氨酸蛋白酶抑制劑、mTERF蛋白結(jié)構(gòu)域和等存在關(guān)聯(lián)。與、氯化鈉脅迫蛋白和等表現(xiàn)出互作網(wǎng)絡(luò)關(guān)系。與丙酮酸脫羧化酶和MAP激酶等存在互作關(guān)系。已報(bào)道株高基因[23]和[24]分別在Brown和Pink模塊中。生長(zhǎng)素轉(zhuǎn)錄因子和與關(guān)聯(lián)。與光合系統(tǒng)II氧進(jìn)化多肽、光合系統(tǒng)I N亞基、苯丙氨酸氨裂解酶和等存在互作網(wǎng)絡(luò)關(guān)系。

圖6 顯著共表達(dá)模塊內(nèi)核心基因的基因網(wǎng)絡(luò)

圖中青綠色、淺綠色和粉紅色分別表示Turquoise、Lightgreen和Pink模塊。

Colors in figure represent Turquoise, Lightgreen, and Pink modules.

圖7 已報(bào)道株高基因D8、DWF1、ZmGRF10、An1和GA20ox3的基因網(wǎng)絡(luò)

圖中青綠色、棕色和粉紅色分別表示Turquoise、Brown和Pink模塊。

Colors in figure represent Turquoise, Brown, and Pink modules.

3 討論

隨著測(cè)序技術(shù)的發(fā)展, 利用共表達(dá)網(wǎng)絡(luò)對(duì)基因組、轉(zhuǎn)錄組和代謝組等海量組學(xué)數(shù)據(jù)進(jìn)行系統(tǒng)研究, 是目前比較流行的大數(shù)據(jù)處理方法。WGCNA是共表達(dá)網(wǎng)絡(luò)分析有效的分析方法, 能夠特異地篩選出與目標(biāo)性狀具有高度生物學(xué)意義的共表達(dá)模塊, 在玉米等植物中已經(jīng)被證明是一種高效的數(shù)據(jù)挖掘方法[25]。

本研究利用WGCNA方法對(duì)8個(gè)玉米材料轉(zhuǎn)錄組測(cè)序數(shù)據(jù)構(gòu)建了加權(quán)基因共表達(dá)網(wǎng)絡(luò), 兩種密度條件下挖掘出株高和穗位高高度相關(guān)的模塊共15個(gè), 其中株高和穗位高相同的模塊6個(gè)。通過對(duì)模塊內(nèi)基因的注釋分析, 發(fā)現(xiàn)共表達(dá)基因主要與發(fā)育、光合作用、響應(yīng)光刺激、植物激素和碳水化合物合成/代謝等重要活動(dòng)相關(guān)。Wang等[26]對(duì)3個(gè)雜交種在拔節(jié)期、大喇叭口期和抽雄期的莖進(jìn)行了轉(zhuǎn)錄組測(cè)序分析, 利用WGCNA挖掘了株高相關(guān)的共表達(dá)模塊, 發(fā)現(xiàn)模塊內(nèi)基因主要參與植物激素相關(guān)代謝途徑。已經(jīng)克隆的株高基因如、和參與赤霉素合成和信號(hào)轉(zhuǎn)導(dǎo)途徑[8,10],基因參與調(diào)控生長(zhǎng)素極性運(yùn)輸[11]。Wang等[26]挖掘的株高基因與本研究高低密度株高基因相同的分別有364個(gè)和1833個(gè), 與高低密度穗位高基因相同的分別有725個(gè)和1722個(gè)。除了植物激素相關(guān)的途徑外, Wang等也發(fā)現(xiàn)抽雄期3個(gè)雜交種間的差異基因顯著富集在糖代謝和分生組織發(fā)育相關(guān)過程[26]。因此, 這些重要的生物學(xué)活動(dòng)可能與株高和穗位高的形成有關(guān)。

根據(jù)基因的連接度, 我們篩選了株高和穗位高顯著模塊內(nèi)的核心基因, 并對(duì)其基因網(wǎng)絡(luò)進(jìn)行了研究。其中, 乙烯響應(yīng)因子是低密度條件下株高和穗位高顯著相關(guān)模塊Turquoise的核心基因。乙烯響應(yīng)因子家族在激素、非生物學(xué)脅迫和生長(zhǎng)發(fā)育等方面具有重要作用[27-29]。本研究中也顯著富集在生長(zhǎng)素響應(yīng)和其刺激響應(yīng)、生長(zhǎng)素介導(dǎo)的信號(hào)途徑、發(fā)育過程、非生物學(xué)脅迫和碳水化合物代謝等過程。該基因與已報(bào)道株高基因具有較高的關(guān)聯(lián)關(guān)系。其他株高基因和也與存在關(guān)聯(lián)。這些表明可能是影響株高的重要基因, 但其作用機(jī)制有待進(jìn)一步探究。

已報(bào)道的株高相關(guān)基因如、、、和均存在于本研究株高和穗位高共表達(dá)模塊中。其中、和在45,000株hm–2條件下株高和穗位高顯著相關(guān)的Turquoise模塊中。是赤霉素響應(yīng)突變體基因, 其突變體在DELLA區(qū)、VHYNP區(qū)及同時(shí)這2個(gè)區(qū)發(fā)生氨基酸不等數(shù)目的缺失導(dǎo)致赤霉素信號(hào)傳導(dǎo)受阻而產(chǎn)生矮稈或半矮稈表型[10,30]。是擬南芥油菜素類固醇合成基因的同源基因, 其轉(zhuǎn)基因玉米植株表現(xiàn)出不同程度的矮稈情況[21]。是生長(zhǎng)調(diào)節(jié)轉(zhuǎn)錄因子, 其過表達(dá)會(huì)通過降低細(xì)胞增殖從而導(dǎo)致玉米植株的降低和葉片的變小, 但產(chǎn)量相關(guān)性狀并沒有發(fā)生改變[22]。除了前面提到的Turquoise模塊的核心基因外, 我們發(fā)現(xiàn)很多轉(zhuǎn)錄因子與這3個(gè)株高基因存在網(wǎng)絡(luò)關(guān)系。例如C2C2-CO-like轉(zhuǎn)錄因子與這3個(gè)株高基因均存在關(guān)聯(lián)關(guān)系。bZIP轉(zhuǎn)錄因子、、、和與和均存在互作關(guān)系。和與和存在互作網(wǎng)絡(luò)關(guān)系。同源框轉(zhuǎn)錄因子和分別與和具有關(guān)聯(lián)關(guān)系。盡管這些轉(zhuǎn)錄因子與已報(bào)道株高基因的權(quán)重值不在其局部網(wǎng)絡(luò)的前10%內(nèi), 但其與已知株高基因網(wǎng)絡(luò)關(guān)系的發(fā)掘有助于揭示株高的遺傳機(jī)制。

和也是株高相關(guān)的基因, 分別在Brown和Pink模塊中。通過影響油菜素類固醇合成而導(dǎo)致玉米莖節(jié)間長(zhǎng)度縮短[23]。生長(zhǎng)素轉(zhuǎn)錄因子和與具有互作關(guān)系。在擬南芥中,受到遠(yuǎn)紅外光的刺激會(huì)導(dǎo)致葉柄伸長(zhǎng)[24], 進(jìn)而導(dǎo)致植株高度變化。已克隆株高基因是調(diào)控玉米株高的主效QTL的候選基因[31]。赤霉素累積會(huì)引起該主效QTL染色體片段代換系的節(jié)間增長(zhǎng)[31]。Pink模塊中,與光合系統(tǒng)II氧進(jìn)化多肽、光合系統(tǒng)I N亞基、苯丙氨酸氨裂解酶和等存在互作關(guān)系。

此外, QTL定位和全基因組關(guān)聯(lián)分析方法鑒定的株高和穗位高候選基因也在本研究的共表達(dá)模塊內(nèi)。例如, 鄭雷等在元分析整合的株高一致性QTL置信區(qū)間內(nèi)發(fā)掘的候選基因(編碼棕色葉中脈)[32]存在于45,000株hm–2下顯著模塊Lightyellow中。45,000株hm–2下發(fā)現(xiàn)的與株高顯著相關(guān)的基因是轉(zhuǎn)錄因子, 也是穗位高顯著關(guān)聯(lián)的SNP (SYN32392)的候選基因[4]。45,000株hm–2下發(fā)現(xiàn)的株高顯著相關(guān)基因(線粒體底物載體家族蛋白)和(亮氨酸重復(fù)蛋白激酶家族蛋白)分別是株高顯著關(guān)聯(lián)SNP 位點(diǎn)S3_18231963和S6_ 164148781的候選基因[5]。而且, 這2個(gè)SNP解釋株高表型變異率比較高, 均為9.3%[5]。45,000株hm–2條件下的株高和穗位高共同基因是穗位高顯著關(guān)聯(lián)位點(diǎn)S2_25387853的候選基因[5]。67,500株hm–2條件下的株高和穗位高共同基因(轉(zhuǎn)導(dǎo)蛋白/WD40重復(fù)超家族蛋白)是Li等[5]研究中穗位高關(guān)聯(lián)位點(diǎn)S5_69308290的候選基因。

上述前人研究結(jié)果與本研究株高和穗位高共表達(dá)模塊基因存在諸多重合, 說明利用WGCNA從轉(zhuǎn)錄組水平鑒定株高和穗位高共表達(dá)模塊及其生物學(xué)意義是可行的。因此, 利用WGCNA方法解析株高和穗位高等重要農(nóng)藝性狀, 分析其共表達(dá)模塊的生物學(xué)功能, 挖掘顯著關(guān)聯(lián)模塊內(nèi)的關(guān)鍵基因, 可以為解析復(fù)雜農(nóng)藝性狀的分子機(jī)制提供重要途徑。

4 結(jié)論

利用WGCNA方法構(gòu)建了不同密度條件下株高和穗位高的加權(quán)共表達(dá)網(wǎng)絡(luò), 得到高度關(guān)聯(lián)模塊15個(gè)。揭示了株高和穗位高關(guān)聯(lián)模塊的生物學(xué)功能, 鑒定了關(guān)聯(lián)模塊內(nèi)核心基因, 解析了其基因網(wǎng)絡(luò)關(guān)系。

[1] 李清超, 李永祥, 楊釗釗, 劉成, 劉志齋, 李春輝, 彭勃, 張巖, 王迪, 譚巍巍, 孫寶成, 石云素, 宋燕春, 張志明, 潘光堂, 黎裕, 王天宇. 基于多重相關(guān)RIL群體的玉米株高和穗位高QTL定位. 作物學(xué)報(bào), 2013, 39: 1521–1529. Li Q C, Li Y X, Yang Z Z, Liu C, Liu Z Z, Li C H, Peng B, Zhang Y, Wang D, Tan W W, Sun B C, Shi S Y, Song C Y, Zhang Z M, Pan G T, Li Y, Wang T Y. QTL mapping for plant height and ear height by using multiple related RIL populations in maize., 2013, 39: 1521–1529 (in Chinese with English abstract).

[2] 何坤輝, 常立國(guó), 崔婷婷, 渠建洲, 郭東偉, 徐淑兔, 張興華, 張仁和, 薛吉全, 劉建超. 多環(huán)境下玉米株高和穗位高的QTL定位. 中國(guó)農(nóng)業(yè)科學(xué), 2016, 49: 1443–1452. He K H, Chang L G, Cui T T, Qu J Z, Guo D W, Xu S T, Zhang X H, Zhang R H, Xue J Q, Liu J C. Mapping QTL for plant height and ear height in maize under multi-environments., 2016, 49: 1443–1452 (in Chinese with English abstract).

[3] 劉坤, 張雪海, 孫高陽, 閆鵬帥, 郭海平, 陳思遠(yuǎn), 薛亞東, 郭戰(zhàn)勇, 謝惠玲, 湯繼華, 李衛(wèi)華. 玉米株型相關(guān)性狀的全基因組關(guān)聯(lián)分析. 中國(guó)農(nóng)業(yè)科學(xué), 2018, 51: 821–834. Liu K, Zhang X H, Sun G Y, Yan P S, Guo H P, Chen S Y, Xue Y D, Guo Z Y, Xie H L, Tang J H, Li W H. Genome-wide association studies of plant type traits in maize., 2018, 51: 821–834 (in Chinese with English abstract).

[4] 李凱, 張曉祥, 管中榮, 沈亞歐, 潘光堂. 玉米株高和穗位高的全基因組關(guān)聯(lián)分析. 玉米科學(xué), 2017, 25(6): 1–7. Li K, Zhang X X, Guan Z R, Shen Y O, Pan G T. Genome-wide association analysis of plant height and ear height in maize., 2017, 25(6): 1–7 (in Chinese with English abstract).

[5] Li X, Zhou Z, Ding J, Wu Y, Zhou B, Wang R, Ma J, Wang S, Zhang X, Xia Z, Chen J, Wu J. Combined linkage and association mapping reveals QTL and candidate genes for plant and ear height in maize., 2016, 7: 833.

[6] Weng J, Xie C, Hao Z, Wang J, Liu C, Li M, Zhang D, Bai L, Zhang S, Li X. Genome-wide association study identifies candidate genes that affect plant height in Chinese elite maize (L.) inbred lines., 2011, 6: e29229.

[7] Fujioka S, Yamane H, Spray C R, Gaskin P, Macmillan J, Phinney B O, Takahashi N. Qualitative and quantitative analyses of gibberellins in vegetative shoots of normal,,,, andseedlings ofL., 1988, 88: 1367–1372.

[8] Winkler R G, Helentjaris T. The maizegene encodes acytochrome P450-mediated early step in gibberellin biosynthesis., 1995, 7: 1307–1317.

[9] Thornsberry J M, Goodman M M, Doebley J, Kresovich S, Nielsen D, Buckler E S.polymorphisms associate with variation in flowering time., 2001, 28: 286–289.

[10] Lawit S J, Wych H M, Xu D, Kundu S, Tomes D T. Maize DELLA proteinsandas modulators of plant development., 2010, 51: 1854–1868.

[11] Multani D S, Briggs S P, Chamberlin M A, Blakeslee J J, Murphy A S, Johal G S. Loss of an MDR transporter in compact stalks of maizeand sorghummutants., 2003, 302: 81–84.

[12] Zhang B, Horvath S. A general framework for weighted gene co-expression network analysis., 2005, 4: Article 17.

[13] Zhang X, Hirsch C N, Sekhon R S, De Leon N, Kaeppler S M. Evidence for maternal control of seed size in maize from phenotypic and transcriptional analysis., 2016, 67: 1907–1917.

[14] Ma J, Zhang D, Cao Y, Wang L, Li J, Lübberstedt T, Wang T, Li Y, Li H. Heterosis-related genes under different planting densities in maize (L.)., 2018, 69: 5077–5087.

[15] Zhan J, Thakare D, Ma C, Lloyd A, Nixon N M, Arakaki A M, Burnett W J, Logan K O, Wang D, Wang X, Drews G N, Yadegari R. RNA sequencing of laser-capture microdissected compartments of the maize kernel identifies regulatory modules asso-ciated with endosperm cell differentiation., 2015, 27: 513–531.

[16] 楊宇昕, 桑志勤, 許誠(chéng), 代文雙, 鄒棖. 利用WGCNA進(jìn)行玉米花期基因共表達(dá)模塊鑒定. 作物學(xué)報(bào), 2019, 45: 161–174. Yang Y X, Sang Z Q, Xu C, Dai W S, Zou C. Identification of maize flowering gene co-expression modules by WGCNA., 2019, 45: 161–174 (in Chinese with English abstract).

[17] Peng H, He X, Gao J, Ma H, Zhang Z, Shen Y, Pan G, Lin H. Transcriptomic changes during maize roots development responsive to Cadmium (Cd) pollution using comparative RNA seq- based approach., 2015, 464: 1040–1047.

[18] Thirunavukkarasu N, Hossain F, Mohan S, Shiriga K, Mittal S, Sharma R, Singh R K, Gupta H S. Genome-wide expression of transcriptomes and their co-expression pattern in subtropical maize (L.) under waterlogging stress., 2013, 8: e70433.

[19] Lyu Y, Liang Z, Ge M, Qi W, Zhang T, Lin F, Peng Z, Zhao H. Genome-wide identification and functional prediction of nitrogen-responsive intergenic and intronic long non-coding RNAs in maize (L.)., 2016, 17: 350.

[20] Zhang S, Yang W, Zhao Q, Zhou X, Jiang L, Ma S, Liu X, Li Ye, Zhang C, Fan Y, Chen R. Analysis of weighted co-regulatory networks in maize provides insights into new genes and regulatory mechanisms related to inositol phosphate metabolism., 2016, 17: 129–146.

[21] Tao Y, Zheng J, Xu Z, Zhang X, Zhang K, Wang G. Functional analysis of, a maize homolog of thebrassinosteroids biosynthetic DWF1/DIM gene., 2004, 167: 741–751.

[22] Wu L, Zhang D, Xue M, Qian J, He Y, Wang S. Overexpression of the maize, an endogenous truncated growth regulating factor protein, leads to reduction in leaf size and plant height., 2014, 56: 1053–1063.

[23] Hartwig T, Chuck G S, Fujioka S, Klempien A, Weizbauer R, Potluri D P, Choe S, Johal G S, Schulz B. Brassinosteroid control of sex determination in maize., 2011, 108: 19814–19819.

[24] Tamotsu H, Rod W K, Chris A H, Masaji K. The involvement of gibberellin 20-oxidase genes in phytochrome-regulated petiole elongation of., 2005, 138: 1106–1116.

[25] Zhao W, Langfelder P, Fuller T, Dong J, Li A, Hovarth S. Weighted gene coexpression network analysis: state of the art., 2010, 20: 281–300.

[26] Wang H, Gu L, Zhang X, Liu M, Jiang H, Cai R, Zhao Y, Cheng B. Global transcriptome and weighted gene co-expression network analyses reveal hybrid-specific modules and candidate genes related to plant height development in maize., 2018, 98: 187–203.

[27] Mizoi J, Shinozaki K, Yamaguchi-Shinozaki K. AP2/ERF family transcription factors in plant abiotic stress responses., 2012, 1819: 86–96.

[28] Hinz M, Wilson I W, Yang J, Buerstenbinder K, Llewellyn D, Dennis E S, Sauter M, Dolferus R.RAP2: 2. An ethylene response transcription factor that is important for hypoxia survival., 2010, 153: 757–772.

[29] Licausi F, Ohme Takagi M, Perata P. APETALA2/ethylene responsive factor (AP2/ERF) transcription factors: mediators of stress responses and developmental programs., 2013, 199: 639–649.

[30] Cassani E, Bertolini E, Cerino B F, Landoni M, Gavina D, Sirizzotti A, Pilu R. Characterization of the first dominant dwarf maize mutant carrying a single amino acid insertion in the VHYNP domain of thegene., 2009, 24: 375–385.

[31] Teng F, Zhai L, Liu R, Bai W, Wang L, Huo D, Tao Y, Zheng Y, Zhang Z., a candidate gene for a major QTL,, for plant height in maize., 2013, 73: 405–416.

[32] 鄭雷, 周羽, 曾興, 邸宏, 翁建峰, 李新海, 王振華. 玉米株高QTL定位研究進(jìn)展. 作物雜志, 2016, (2): 8–13.Zheng L, Zhou Y, Zeng X, Di H, Weng J F, Li X H, Wang Z H. QTL Mapping of plant height in maize., 2016, (2): 8–13 (in Chinese with English abstract).

Identification of gene co-expression modules of maize plant height and ear height by WGCNA

MA Juan, CAO Yan-Yong, WANG Li-Feng, LI Jing-Jing, WANG Hao, FAN Yan-Ping, and LI Hui-Yong*

Institute of Cereal Crops, Henan Academy of Agricultural Sciences, Zhengzhou 450002, Henan, China

Plant height (PH) and ear height (EH) are important factors for maize plant type and grain yield. Weighted gene co-expression network analysis (WGCNA) is an important method to explain the relationships between gene network and complicated traits and identify the PH and EH associated genes. In this study, we used Zheng 58, Ye 478, Chang 7-2, Huangzaosi and its combinations Zhengdan 958, Anyu 5, Zheng 58/Huangzaosi, and Ye 478/Huangzaosi as materials and utilized transcriptome data under the planting densities of 45,000 plants hm–2and 67,500 plants hm–2to construct a co-expression network by WGCNA, getting 24 and 21 co-expression modules, respectively. Among them, a total of 15 co-expression modules were significantly correlated with PH and EH, with the absolute correlation coefficients higher than 0.50. Six modules were overlapped between PH and EH. By gene function analysis, these overlapped modules were significantly enriched in development, photosynthesis, response to light stimulus, plant hormone, and carbohydrate biosynthesis/metabolism related activities. According to connectivity of genes in modules, AP2-EREBP transcription factor, thiaminase, phosphoglyceric kinase, glutathione transferase, and succinate dehydrogenasewere considered as hub genes. From gene networks,was connected with three known PH genes,,, and(C3H transcription factor),(C2C2-GATA transcription factor), and ethylene homology. Reported PH genesandwere also found in our co-expression modules. From the networks of the five known PH genes, ARF-transcription factor 7 (),,, photosystem II oxygen evolving polypeptide, and photosystem I N subunithad connections with these known PH genes. The identification of 15 co-expression modules and their hub genes, and analysis of their gene function and gene networks of key genes will be helpful for revealing the genetic basis of PH and EH.

maize; weighted gene co-expression network; transcriptome; plant height; ear height

2019-04-02;

2019-09-26;

2019-10-10.

10.3724/SP.J.1006.2020.93021

李會(huì)勇, E-mail: lihuiyong1977@126.com

E-mail: majuanjuan85@126.com

本研究由國(guó)家重點(diǎn)研發(fā)計(jì)劃項(xiàng)目(2016YFD100103)和河南省科技攻關(guān)項(xiàng)目(192102110008)資助。

This study was supported by the National Key Research and Development Program of China (2016YFD100103) and the Science and Technology Project of Henan province (192102110008).

URL:http://kns.cnki.net/kcms/detail/11.1809.S.20191010.1139.002.html