整合GWAS和WGCNA篩選鑒定甘藍型油菜生物產量候選基因

王艷花 劉景森 李加納

整合GWAS和WGCNA篩選鑒定甘藍型油菜生物產量候選基因

王艷花1,2,**劉景森1,2,**李加納1,2,*

1西南大學農學與生物科技學院/ 油菜工程研究中心, 重慶 400715;2西南大學現代農業(yè)科學研究院, 重慶 400715

生物產量是作物獲得高產的重要基礎, 對于甘藍型油菜(L.)尤其重要。本研究利用588份甘藍型油菜材料構成的自然群體2年生物產量表型數據的全基因組關聯(lián)分析, 再結合高生物產量材料‘CQ45’和低生物產量材料‘CQ46’的轉錄組測序(RNA-seq)結果, 整合了6個甘藍型油菜材料6個部位(莖稈、葉片、花后30 d主軸與側枝種子、花后30 d主軸與側枝角果皮)的轉錄組數據構建的加權共表達網絡分析(WGCNA), 篩選出與生物產量相關的候選基因。通過相關分析發(fā)現, 2年間甘藍型油菜自然群體中生物產量對大多數產量相關性狀都具有正向效應; 自然群體2年生物產量分析的最佳模型均為K+PCA模型, 共檢測到9個顯著位點(< 1/385691或< 0.05/385691); 根據CQ45和CQ46共36組轉錄組數據, 選擇MAD值為前5%的基因共計5052個用于構建WGCNA, 通過篩選合并共得到了15個模塊, 其中5個基因共表達模塊分別與葉片、莖稈和花后30 d種子顯著性相關; 整合了WGCNA中關鍵模塊的hub gene、GWAS分析得到的顯著SNP位點和極端表型差異基因確定候選基因, 它們的擬南芥同源基因為、、、, 這些基因在光合作用的卡爾文循環(huán)、碳同化、物質積累等方面發(fā)揮重要作用。

GWAS; WGCNA; RNA Seq; 甘藍型油菜

油菜作為四大油料作物(大豆、油菜、向日葵、花生)之一, 由于其較強的適應性, 在世界各地廣泛種植[1-3]。隨著人類對油菜的需求不斷增加, 油菜總產量增速不斷加快[4], 世界油料作物中油菜的種植面積及總產量已位居第二, 僅次于大豆[5-6]。同時油菜也是我國主要農作物, 在作物種植面積中排名第五, 僅次于水稻、玉米、小麥、大豆[7-12]。近年來, 油菜多功能開發(fā)利用的持續(xù)發(fā)展, 形成了油菜的油用、菜用、花用、蜜用、飼用等“一菜多用”的模式, 大幅度的提高了油菜的種植效益[13]。在實際生產中, 我國油菜產量仍較低, 生產成本卻很高, 國際競爭力嚴重不足[14-15]。目前, 我國油菜超過60%依賴進口, 并且種植面積和單位面積產量都停滯不前, 這對我國糧食安全戰(zhàn)略有著嚴重威脅[16]。因此, 提高油菜產量是我國油菜育種不斷追求的目標[17-24]。

生物產量指單位面積土地上獲得的不包括根系的作物干物質的總量, 它是描述莖、葉制造并積累有機物能力強弱的量[25]。生物產量體現作物總的生產力, 或者說作物光合作用碳同化的能力。多項報道已經表明, 生物產量與經濟產量密切相關, 作物提高產量的方法通常有2種: 即在生物產量不變的情況下, 提高收獲指數; 或收獲指數不變的情況下, 提高生物產量[26-30]。近年來, 高通量測序技術的不斷突破、成本不斷降低, 多樣本的轉錄組測序逐漸被應用于系統(tǒng)研究生命科學問題[31], 傳統(tǒng)的少量樣本比較分析已經無法有效處理海量的生物信息數據。因此, 生物信息分析工具和方法應運而生[32]。其中, 加權基因共表達網絡分析(weighted gene co-expression network analysis, WGCNA)能夠特異篩選出與目標性狀高度關聯(lián)的基因, 并進行模塊(module)化分類, 得到具有高度生物學意義的共表達模塊(co-expression module), 能夠有效篩選到核心基因(hub genes)[33]。WGCNA算法作為一種精準、高效的生物信息學及生物數據挖掘方法, 被廣泛應用于生物學各個領域, 利用WGCNA對馬鈴薯C16和C119品種構建抗逆生理性狀高度關聯(lián)的權重網絡, 為深入探究馬鈴薯抗逆分子機制提供新的數據資源[34-35]。本研究以生物產量性狀為主要研究對象, 通過GWAS分析, 鑒定出與生物產量性狀緊密相關的SNP位點; 然后提取SNP位點前后共500 kb的序列片段作為候選區(qū)段, 提取候選區(qū)段內所有基因并結合轉錄組分析結果篩選獲得差異候選基因; 同時, 結合轉錄組數據構建WGCNA共表達網絡, 篩選鑒定出甘藍型油菜生物產量關鍵候選基因。

1 材料與方法

1.1 供試材料

以來自世界各地的588份油菜自交系作為材料, 其中455份來自中國(大部分來自重慶、湖南、湖北等地), 11份來自亞洲其他地區(qū), 102份來自歐洲, 13份來自北美洲, 7份來自澳大利亞。所有材料均由重慶市油菜工程技術研究中心提供。

1.2 材料種植及農藝性狀考察

588份材料于2016、2018年種植于重慶北碚歇馬油菜基地(29o45′39.99″, 106o22′38.47″, 海拔238.57 m)。每年播種時間為9月20日左右, 按完全隨機區(qū)組設計, 2個重復同時種植, 每個材料種植2行, 每行15株, 行距40 cm, 株距20 cm。田間管理選擇常規(guī)生產方式, 2個重復按照相同的管理措施和施肥水平進行, 收獲時間為次年5月5日左右。于成熟期, 在每個小區(qū)中選擇5株長勢一致的植株。將植株分成上部分枝和下部莖稈兩部分, 并將上部分枝裝進網袋中。上部分枝和下部莖稈分別徹底曬干后, 使用天平秤進行稱量, 將兩部分重量相加得到植株的總干重。每份材料的5個樣本干重的平均值為該材料的生物產量(biomass yield, BY)。

1.3 表型數據分析

利用Microsoft Excel軟件初步整理588份自然群體的表型數據, 并計算單個環(huán)境及差值條下生物產量的平均值、標準差、變異系數等, 并用SPSS繪制生物產量的正態(tài)分布圖; 利用SPSS統(tǒng)計分析軟件對自然群體進行各性狀進行相關分析; 利用R軟件對多年環(huán)境下自然群體的生物產量表型數據進行最佳線性無偏預測(best linear unbiased prediction, BLUP)處理。

1.4 全基因組關聯(lián)分析

利用本實驗室已有588份甘藍型油菜重測序的數據, 通過基因型分析, 最終獲得385,691個可利用的SNP標記數據[36]。本研究利用這些標記結合重慶2017、2019年和BLUP環(huán)境的生物產量(BY)表型數據行全基因組關聯(lián)分析。本研究使用6種統(tǒng)計模型在Tassel 5.2.1[37]軟件中執(zhí)行對甘藍型油菜生物產量的關聯(lián)分析, 即: Q模型、naive模型、K模型、PCA模型、Q+K模型和PCA+K模型。Q代表群體結構, PCA代表主成分, K代表親緣關系; naive、Q和PCA模型在一般線性模型(general linear model, GLM)下運行, K、Q+K和PCA+K模型在混合線性模型(mixed linear model, MLM)下運行。利用SAS軟件對這6種模型的運算結果中的-lg()的觀測值和-lg()期望值繪制QQ (Quantile-quantile)圖。通過比較QQ圖確定最佳模型, 并選擇最佳模型下生物產量性狀的全基因組關聯(lián)分析結果作為進一步分析對象。本試驗使用的SNP數據為385,691個, 因此規(guī)定值小于閾值(1/385,691=2.593E-06和0.05/385,691=1.296E-07)的位點為顯著關聯(lián)位點。利用Haploview[38]軟件制作Manhattan圖, 將通過關聯(lián)分析檢測到的染色體組上與各性狀顯著相關的標記位點可視化。

1.5 轉錄組測序分析

從供試群體中選擇1個高生物產量(CQ45)和1個低生物產量(CQ46)材料, 分別選取蕾苔期莖稈(J)和葉片(Le)、采用掛花標記的方法選取開花后30 d主軸種子(30S)、主軸角果皮(30P)、側枝種子(30CS)和側枝角果皮(30CP) 6個時期的樣品, 每個樣品取2個生物重復, 送北京諾禾致源生物信息科技有限公司進行轉錄組測序, 在OmicShare (基迪奧)云平臺完成(http://www.omicshare.com) GO和KEGG富集分析, 篩選差異表達基因。

1.6 加權共表達網絡分析

使用R軟件包preprocessCore中的normalize. quantiles函數對轉錄組數據進行標準化。然后計算每個基因的MAD (median absolute deviation), 選擇的MAD值為前5%的基因作為樣本間變異大的基因群用于WGCNA構建。使用average-linkage層次聚類法對基因進行聚類, 按照混合動態(tài)剪切樹的標準, 并設置每個基因網絡模塊最少的基因數為30。在使用動態(tài)剪切法在確定基因模塊后, 我們依次計算每個模塊的特征向量值(eigengenes), 然后對模塊進行聚類分析, 將距離較近的模塊合并成新的模塊, 設置height = 0.25、deepSplit = 3、minModuleSize = 30。計算所得到模塊的特征向量與樣本來源組織的相關性, 選擇相關系數>0.6的模塊作為與組織部位顯著相關模塊。通過特征向量基因分析確定目標基因模塊, 對所有的模塊中具有代表性的基因-特征向量基因(module eigengene, ME)進行聚類分析, 進一步進行兩兩模塊之間的ME相關性分析, ME之間的相關性越高, 其所在的模塊的相關性也就越高; 將模塊中的所有基因和ME基因在所有樣本中分別進行表達水平分析, GO和KEGG分析。最后使用Cytoscape[39]將所關注模塊的ME基因的共表達網絡進行可視化。

2 結果與分析

2.1 表型數據統(tǒng)計及各性狀的相關性分析

對2017年和2019年588份自然群體的生物產量的表型變異情況進行統(tǒng)計分析發(fā)現, 在不同環(huán)境下自群體的生物產量呈連續(xù)變異, 近似正態(tài)分布, 說明生物產量是受多基因控制的數量性狀, 符合GWAS分析的要求。各產量相關性狀的相關分析結果表明, 每角果粒數在2017年和2019年與莖稈干重、經濟產量、收獲指數和生物產量均呈正相關, 而莖稈干重僅在2017年呈現出顯著正相關; 莖稈干重與上部生物產量、經濟產量、生物產量在2年間均呈現顯著正相關; 上部生物產量與經濟產量、收獲指數、生物產量在2年間均呈現顯著正相關; 而經濟產量和收獲指數與生物產量在2年間同樣呈現正相關。多數性狀間的相關性均達到顯著或極顯著水平, 表明產量相關的各性狀間存在顯著相關作用(表1)。

2.2 甘藍型油菜生物產量性狀的GWAS分析

各生物產量性狀在6種模型下的關聯(lián)分析后得到的結果進行Quantile-Quantile散點圖(QQ Plot)的繪制, 選擇與預測值最接近的一種模型作為最佳模型, 2017年、2019年和BLUP環(huán)境下最佳模型為K+PCA模型(圖1)。在最佳模型下, 利用本實驗室已有的385,691個SNP, 以值小于閾值(1/385,692 = 2.593E-06)確定顯著關聯(lián)SNP位點, 并繪制曼哈頓圖(圖2)。生物產量2017、2019年和BLUP環(huán)境下共檢測到9個顯著關聯(lián)標記位點(表2)。其中, 在2017年環(huán)境中檢測到7個顯著關聯(lián)標記位點, 分別位于A03、A07、C01 (2個)、C04 (2個)、C06染色體, 并且位于同一條染色體的顯著關聯(lián)位點其置信區(qū)間高度重合; 2019年檢測到2個顯著關聯(lián)位點, 分別位于A07和C09染色體; 而BLUP環(huán)境下共檢測到4個顯著關聯(lián)位點, 分別位于A07、C03、C06 (2個), 位于C06染色體的2個顯著關聯(lián)位點其置信區(qū)間高度重合。9個SNP的貢獻率在5.64%~7.98%, 其中2017年檢測到的位于C06染色體的SNP (S16_ 26169577)貢獻率最高。

表1 自然群體2年各產量相關性狀間的相關分析

*、**分別表示在0.05和0.01水平相關性顯著(雙尾)。.

*,**: significant correlation at the 0.05 and 0.01 probability levels, respectively. SWSI: seeds weight per silique index; TSW: 1000-grain weight; SNPS: seeds number per silique; ST: stem dry weight; CBY: canopy biomass yield; SY: seed yield; HI: harvest index; BY: biomass yield.

2.3 極端表型材料RNA-Seq分析

在葉片中檢測到極端材料的差異基因共6820個, CQ45對CQ46顯著上調表達的基因有3516個, 顯著下調表達的基因有3304個(圖3-A)。生物分子功能類中相關性最高的是氧化還原酶活性(GO:0016903)、催化活性(GO:0003824)等, 而參與基因數量較多的生物過程是單體代謝過程(GO:0044710)、光合作用(GO:0015979)等(圖3-B); 葉片中差異基因主要集中在碳水化合物代謝途徑(圖3-C)。

花后30 d主軸角果皮中檢測到17,309個顯著差異表達基因, CQ45對于CQ46有2682個顯著上調表達基因, 3570個顯著下調基因(圖3-A)。差異基因中分子功能相關性較高的是葡萄酶轉移(GO:0046527)等, 生物過程相關性最高的有應激反應(GO:0006950)、單體生物合成過程(GO:0044711)等(圖3-B); 相關性較高的途徑是氨基酸的生物合成以及脂質代謝途徑(圖3-C)。

花后30 d側枝角果皮中檢測到5431個顯著差異表達基因, 其中CQ45對于CQ46有2173個顯著上調表達基因, 3258個顯著下調基因(圖3-A)。側枝角果皮中分子功能相關性最高的是營養(yǎng)庫活性(GO:0045735)等, 參與的生物過程中相關性較高的是對含氧化合物的反應(GO:1901700)等(圖3-B); 這些差異基因相關性最高的是能量代謝途徑(圖3-C)。

花后30 d主軸種子中檢測到6948個顯著差異表達基因, CQ45對于CQ46有4080個顯著上調表達基因, 2868個顯著下調基因(圖3-A)。分子功能相關性最高的是微管運動活性(GO:0003777)、硫葡萄糖苷酶活性(GO:0019137)等, 生物過程相關性較高的是細胞循環(huán)(GO:0007049)、細胞分裂(GO:0051301) (圖3-B); 相關性較高的代謝途徑是碳水化合物代謝中的淀粉和蔗糖代謝(圖3-C)。

花后30 d側枝種子中檢測到17,309個顯著差異表達基因, CQ45對于CQ46有11489個顯著上調表達基因, 5820個顯著下調基因(圖3-A)。分子功能相關性最高的是微管運動活性(GO:0003777)、營養(yǎng)庫活性(GO:0045735)硫葡萄糖苷酶活性(GO:0019137)等, 參與的生物過程相關性較高的是細胞循環(huán)(GO:0007049)等, 其結果與30 d主軸角果種子比較相似(圖3-B); 相關性較高的是氨基酸代謝中的丙氨酸、天冬氨酸和谷氨酸代謝(圖3-C)。

表2 最佳模型下生物產量顯著位點表

(圖3)

(圖3)

(圖3)

A: 極端表型材料各組織差異基因數量; B: 極端表型材料各組織差異基因的GO富集分析; C: 極端表型材料各組織差異基因的KEGG分析。

A: the number of differentially expressed genes in tissues of extreme phenotypic materials; B: GO enrichment analysis of differentially expressed genes in tissues of extreme phenotypic materials; C: KEGG analysis of differentially expressed genes in tissues of extreme phenotype materials.

莖稈中檢測到極端表型材料的差異基因12,867個, CQ45對于CQ46有6452個顯著上調表達基因, 6415個顯著下調基因(圖3-A)。分子功能相關性最高的是轉移酶活性(GO:0016758)、UDP葡萄糖基轉移酶活性(GO:0035251)等, 參與較多的生物過程是植物型次生細胞壁生物發(fā)生(GO:0009834)、細胞壁合成(GO:0042546)等(圖3-B); 參與相關性較高的代謝途徑是到次生代謝產物的生物合成等途徑(圖3-C)。

2.4 加權共表達網絡(WGCNA)分析

2.4.1 基因共表達網絡的模塊生成 根據36個樣本得到的RNA-Seq的結果, 計算各基因的MAD值, 取MAD為前5%共計5052個基因, 進行基因表達水平的層次聚類(圖4-A); 使用R軟件包WGCNA構建權重共表達網絡, 選取擬合曲線第1次接近0.9時的閾值參數β (β=10)篩選共表達模塊(圖4-B, C),共表達網絡符合無尺度網絡。使用動態(tài)剪切法在確定基因模塊后, 依次計算每個模塊的特征向量值, 然后對模塊進行聚類分析, 將距離較近的模塊合并成新的模塊后共得到了15個模塊(圖4-D), 聚類樹中每個枝代表1個模塊, 分別是Black、Cyan、Green、Magenta、Pink、Purple、Red、Salmon、Tan、Turquoise、Yellow、Greenyellow、Brown、Blue模塊, Grey模塊是無法聚集到其他模塊的基因集合(包含11個基因)。在15個模塊中基因數量33~1843不等, 其中包含基因數量最多的是Turquoise模塊(1843個), 最少的是Cyan模塊(33個), 各個基因之間的表達關系如圖4-E。

2.4.2 模塊與組織之間相關性 本研究分別計算這15個模塊的特征向量與樣本來源組織的相關性, 得到5個顯著相關的模塊, 其中與葉片(L)極顯著相關的模塊是Turquoise模塊, 相關系數達到0.77; 與莖稈顯著相關的模塊為Green、Magenta、Pink模塊, 其相關系數分別為0.9、0.62、0.75; 與30 d主軸種子顯著相關的模塊是Yellow模塊, 其相關系數為0.6 (圖5-A)。

通過特征向量分析確定目標基因模塊, 對所有的模塊中具有代表性的基因-特征向量基因(module eigengene, ME)進行聚類分析(圖5-B), ME之間的相關性越高, 其所在的模塊的相關性也就越高, 進一步進行兩兩模塊之間的ME相關性分析(圖5-C); 將模塊中的所有基因和ME基因在所有樣本中分別進行表達水平分析發(fā)現, 各個模塊內的基因其在不同組織部位中的表達情況呈現不同的特性。每個模塊內的基因表達水平高度相關, ME表達水平與模塊整體表達水平也高度相關(圖5-D和附圖1), 說明目標模塊的ME可充分代表其模塊中的整體基因。

2.4.3 目標基因模塊中的關鍵基因 為獲得與組織部位生長發(fā)育和功能高度相關的5個模塊中的核心基因, 利用Cytoscape軟件對基因互作調控網絡進行可視化處理, 篩選出模塊中共表達權重>0.1并且連通性排名前10的基因, 并結合NCBI數據庫確定模塊內的關鍵基因(圖6)。Turquoise模塊中確定了6個關鍵基因, 分別為、、、、、。Green模塊中篩選出了4個核心基因, 分別為、、、。Magenta模塊中篩選出了4個核心基因, 分別為、、、。Pink模塊中篩選出了6個核心基因, 分別為、、、、、。Yellow模塊中篩選出了8個核心基因, 分別為、、、、、、、。

2.5 候選基因的篩選

檢測到顯著性SNP位點上下游500 kb的LD置信區(qū)間內尋找甘藍型油菜基因, 結合轉錄組差異表達基因, 通過WGCNA篩選到的與組織部位生長發(fā)育高度相關的關鍵模塊中的核心基因(hub genes), 初步確定與甘藍型油菜生物產量相關的基因, 將這些基因的蛋白質序列與擬南芥基因蛋白質序列進行BLAST對比, 結合前人已報道的擬南芥同源基因的功能以及特性, 篩選到生物產量性狀中發(fā)揮重要作用的候選基因, 在擬南芥中的同源基因分別是、、、、、(表3)。

和在擬南芥中的同源基因是, 編碼擬南芥葉綠體-1,6-二磷酸酶[40]。豌豆中葉綠體-1,6-二磷酸果糖()干擾載體的轉基因植株表現出葉鮮重明顯增加、光合作用的測量結果顯示具有更高的碳同化率, 表明FBPase作為二氧化碳同化中的關鍵酶的作用, 并且還可以協(xié)調碳和氮的代謝[41]。在擬南芥中的同源基因是, 編碼S-腺苷同型半胱氨酸水解酶(SAHH)。SAHH是維持真核生物甲基化穩(wěn)態(tài)的關鍵酶[42], 會競爭性地抑制甲基轉移酶(MT)活性, 使植株生長緩慢, 繁殖力低, 種子發(fā)芽減少[43-44]。在擬南芥的同源基因是, 編碼在葉綠體基質中發(fā)現的一種小肽, 在黑暗中氧化的CP12與甘油醛-3-磷酸脫氫酶和磷酸布洛激酶(碳同化循環(huán)的兩種酶)形成一種非活性超分子復合物[45],的轉錄物在決定成熟葉片和光合作用能力方面起著重要作用[46-48], Marri等人發(fā)現, CP12蛋白氧化對植物不同環(huán)境條件下光合過程的整體平衡至關重要[49]。在擬南芥中的同源基因是一種卡爾文循環(huán)酶, 具有光合固碳作用[48]。在擬南芥中的同源基因為, 磷酸化甘油醛-3-P脫氫酶(GAPC-1)是一種高度保守的胞漿酶, 純合缺失突變株表現出生長延遲、角果形態(tài)改變和種子數量低[51]。在擬南芥中的同源基因是, 是擬南芥營養(yǎng)組織中主要肌動蛋白基因, 對植物生長發(fā)育有著極其重要的作用[52]。多項研究結果表明, 甘藍型油菜生物產量相關候選基因的同源基因功能, 在植物氧化還原、能量和碳水化合物代謝、光合作用、物質積累等生長發(fā)育過程中發(fā)揮作用。

A: 樣本聚類與數據矯正; B和C: 基于規(guī)模獨立性和均值連通性選擇軟閾值; D: 已識別模塊的樹狀聚類圖; E: 已識別模塊的熱圖。

A: sample clustering to detect outliers; B and C: the selection of the soft threshold based on scale independence and mean connectivity; D: cluster dendrogram of the identified modules; E: the heatmap of identified modules.

A: 基因共表達網絡模塊與組織部位關聯(lián)熱圖; B: 各樣本中Turquoise模塊的所有基因與相應ME的表達水平; C: 不同模塊兩兩之間ME的相關性; D: ME聚類樹。

A: association analysis of gene co-expression network modules with tissues; B: expression levels of all genes and corresponding ME in turquoise module of each sample; C: ME correlation between different modules; D: ME cluster tree.

(圖6)

(圖6)

表3 生物產量候選基因

3 討論

生物產量體現作物總的生產力, 對作物的經濟產量有著密不可分的關系。本研究所使用甘藍型油菜自然群體中的生物產量變化范圍廣、變異系數大, 說明生物產量是比較復雜的性狀。本研究在2017、2019年和BLUP環(huán)境下共檢測到9個生物產量顯著關聯(lián)標記位點, 貢獻率在5.64%~7.98%。與生物產量顯著關聯(lián)的SNP位點分布在各個染色體, Lu等[53]利用520份來自世界各地的甘藍型油菜材料作為自然群體, 使用MLM進行全基因組關聯(lián)分析, 在除了A01、A06、C02、C07以外的15條染色體上共檢測到26個生物產量SNP, 貢獻率在4.33%~17.03%之間; 60K SNP芯片對520分材料構成的自然群體進行全基因組關聯(lián)分析, 在除了A01、A02和A10意外的16條染色體共檢測到89個生物產量SNP位點, 貢獻率在3.47%~8.36%之間[54]; Luo等[55]利用155分甘藍型油菜材料構成的自然群體并使用60K SNP芯片檢測到1個生物產量SNP, 其貢獻率為0.76%。本研究中, 生物產量在2017、2019年和BLUP環(huán)境下共檢測到9個顯著關聯(lián)標記位點, 分別位于A03、A07、C01、C03、C04、C06、C09染色體, 貢獻率在5.64%~7.98%。前人檢測到與生物產量顯著關聯(lián)的部分SNP位點也定位在這些染色體上, 說明生物產量的構成因素比較復雜, 且由多個SNP位點協(xié)同控制。

通過RNA-Seq對極端表型材料的差異表達基因分析表明, 葉片和角果皮主要富集到光合作用、淀粉和蔗糖代謝等生物途徑, 這些代謝途徑可能在油菜植株的物質積累過程中發(fā)揮重要作用, 因此, 這些差異基因的不同表達模式可能是造成2個材料生物產量有顯著差異的主要原因。莖稈中的差異表達基因主要富集在細胞壁合成等途徑, 油菜莖稈次生細胞壁的發(fā)育對其抗倒伏有著決定性影響, 莖稈的發(fā)育狀況良好對油菜碳水化合物的運輸積累有著積極作用[55]。因此, 本研究極端材料莖稈中的差異基因有可能在這些方面影響其生物產量。種子是油菜主要的儲存和收獲器官, 其發(fā)育狀況直接影響著最終產量[56], 因此也直接影響生物產量。本研究極端生物產量材料種子中的差異基因主要富集在營養(yǎng)庫活性、蔗糖淀粉代謝等功能和途徑。

WGCNA是共表達網絡分析有效的分析方法, 能夠特異地篩選出與目標性狀具有高度生物學意義的共表達模塊, 在玉米等植物中已經被證明是一種高效的數據挖掘方法。本研究得到15個生物產量加權共表達網絡模塊, 分析發(fā)現一些模塊與部分組織有著顯著的相關性, 分別是葉片與Turquoise模塊, 莖稈與Green、Magenta、Pink模塊, 開花后30 d主軸種子與Yellow模塊相關。葉片、莖稈和種子均在作物的生物產量的構成中發(fā)揮極其重要的作用, 直接影響到農產品最終的產量。結合GWAS和轉錄組差異基因得到的結果、整合WGCNA挖掘得到關聯(lián)模塊中的核心基因(hub genes), 最終篩選到一批可能與甘藍型油菜生物產量密切相關的基因, 分別是、、、、、、。這些基因的干擾、敲除或者超量表達均引起了植株的整體重量的變化[40-52], 在擬南芥生物產量中發(fā)揮重要作用。

由于生物組學大數據存在復雜、多層次和信息互補的特點, 分析這些數據的一個關鍵目標是確定可預測表型性狀的有效模型, 發(fā)現重要性狀的關鍵調節(jié)因子并闡明其生物功能。本研究根據GWAS分析得到的與生物產量高度關聯(lián)SNP標記, 同時結合轉錄組分析確定的差異基因以及WGCNA得到的基因富集模塊進一步篩選出與甘藍型油菜生物產量相關的基因。通過BLAST對比找到這些基因在擬南芥中的同源基因, 并根據其注釋以及已有報道進一步篩選與生物產量高度相關的候選基因, 為油菜中生物產量候選基因的功能研究奠定基礎。

4 結論

本研究通過2年間甘藍型油菜自然群體中生物產量對大多數產量相關性狀都具有正向效應, 說明甘藍型油菜的生物產量是其他產量相關性狀的基礎和保障。結合GWAS關聯(lián)基因與轉錄組差異基因得到178個基因, 根據擬南芥同源基因的功能, 篩選得到與生物產量性狀相關的和為關鍵候選基因, 它們在能量和碳水化合物代謝以及光合作用等方面中發(fā)揮作用。此外, 選取WGCNA所得到的重點模塊中連通性前10的基因作為hub gene, 通過同源基因功能注釋分析, 得到生物產量候選基因、、、和, 在光合作用的卡爾文循環(huán)、碳同化、物質積累等方面發(fā)揮作用。

[1] 金森森, 羅帥. 全球視角下的南美大豆貿易發(fā)展趨勢. 中國糧食經濟, 2014, (4): 30–32.Jin S S, Luo S. The development trend of South American soybean trade from a global perspective., 2014, (4): 30–32 (in Chinese with English abstract).

[2] 陶伯玉. 甘藍型油菜主要農藝、品質性狀的遺傳分析和QTL定位. 南京農業(yè)大學碩士學位論文, 江蘇南京, 2015. Tao B Y. Genetic Analysis and QTL Mapping of Main Agronomic and Quality Traits in. MS Thesis of Nanjing Agricultural University, Nanjing, Jiangsu, China, 2015 (in Chinese with English abstract).

[3] 張衛(wèi). 明確大宗油料作物發(fā)展規(guī)劃提升國內大宗油料作物自給力. 中國食品, 2016, (19): 82–89. Zhang W. Clarify the development plan of bulk oil crops to enhance the self-sufficiency of domestic bulk oil crops., 2016, (19): 82–89 (in Chinese with English abstract).

[4] 徐露萍. 淺談油菜種植技術的應用與推廣. 南方農業(yè), 2015, 9(3): 17–18. Xu L P. Talking about the application and extension of rape planting technology., 2015, 9(3): 17–18 (in Chinese with English abstract).

[5] 李培武, 丁小霞, 張文, 謝立華, 周海燕. 中加雙低油菜質量標準體系比較及我國發(fā)展對策探討. 農產品質量與安全, 2006, (6): 23–26. Li P W, Ding X X, Zhang W, Xie L H, Zhou H Y. Comparison of quality standard systems for double-low rapeseed between China and Canada and discussion on my country’s development countermeasures., 2006, (6): 23–26 (in Chinese with English abstract).

[6] Mirza H, Karim F M. Rapeseed () cultivation in dry season., 2010, 12: 27–80.

[7] 傅廷棟. 油菜的品種改良. 作物研究, 2007, 21(3): 159–162. Fu T D. Rapeseed variety improvement., 2007, 21(3): 159–162 (in Chinese with English abstract).

[8] Liersch A, Bocianowski J, Henryk W, Laurencja S. Assessment of genetic relationships in breeding lines and cultivars ofand their implications for breeding winter oilseed rape., 2016, 56: 15–40.

[9] Zhao J, Meng J. Genetic analysis of loci associated with partial resistance toin rapeseed (L.)., 2003, 106: 759–764.

[10] Wen Y C, Fu T D, Tu J X, MA C Z, Shen J X, Wen J, Zhang S F. Factors analysis of silique shatter resistance in rapeseed (L.)., 2010, 32: 25–29.

[11] Chen S, Qian T. Analysis on technical efficiency and influencing factors of rapeseed production., 2016, 8: 1–7.

[12] Jiang H D, Zhou Q, Li N, Sun X F. Effect of Cd on the growth and physiological characteristics of rape seedlings., 2006, 28: 39–43.

[13] 張哲, 殷艷, 劉芳, 王積軍, 傅廷棟. 我國油菜多功能開發(fā)利用現狀及發(fā)展對策. 中國油料作物學報, 2018, 40: 618–623. Zhang Z, Yin Y, Liu F, Wang J J, Fu T D. The current status and development countermeasures of multifunctional development and utilization of rapeseed in my country., 2018, 40: 618–623 (in Chinese with English abstract).

[14] Chen H, He H, Zou Y J, Chen W, Yu R B, Liu X, Yang Y, Gao Y M, Xu J L, Fan L M, Li Y, Li Z K, Deng X W. Development and application of a set of breeder-friendly SNP markers for genetic analyses and molecular breeding of rice (L.)., 2011, 123: 869–879.

[15] Ma T Q, Lin L B. Application of molecular marker technologies in rape breeding., 2004, 2: 728–732.

[16] Xiong Q F, Wen J, Li X H, Shen J X. Technological innovation and industrial development of rapeseed in China., 2014, 16: 14–22.

[17] Beecher H G, Dunn B W, Thompson J A, Humphreys E, Mathews S K, Timsina J. Effect of raised beds, irrigation and nitrogen management on growth, water use and yield of rice in south- eastern Australia., 2006, 46: 1363–1372.

[18] Xu Z, Yu Z, Zhao J. Theory and application for the promotion of wheat production in China: past, present and future., 2013, 93: 2339–2350.

[19] Cheng L, Li H P, Qu B, Huang T, Tu J X, Fu T D, Liao Y C. Chloroplast transformation of rapeseed () by particle bombardment of cotyledons., 2010, 29: 371–381.

[20] Wang G C, Liu Z, Yang G S. Development of high oil content lines ofthrough isolated microspore culture., 2007, 29: 382–386.

[21] Liao X, Wang H Z. Present status of rapeseed science and technology innovation system in China and suggestion for further development., 2007, 6: 28–34.

[22] Cheng L, Li H P, Qu B. Chloroplast transformation of rapeseed () by particle bombardment of cotyledons., 2010, 29: 371–381.

[23] Zhang X K, Chen J, Chen L, Wang H Z, Li J N. Imbibition behavior and flooding tolerance of rapeseed seed (L.) with different testa color., 2008, 55: 1175–1184.

[24] Palmieri S, Leoni O. Industrial prospectives of rapeseed oil. Physiologic technological and qualitative aspects., 1992, 14: 63–70.

[25] Weilenmann M E, Lúquez J, Suárez W. Variability of biological yield, economic yield and harvest index in soybean cultivars grown in irrigated and rainfed environments in Argentina., 2000, 21: 39–40.

[26] Li S Y, Zhu R K, Varshney J, Zhan X, Zheng J, Shi X, Wang G, Wang H. A systematic dissection of the mechanisms underlying the natural variation of silique number in rapeseed (L.) germplasm., 2020, 18: 568–580.

[27] Luo X, Ma C Z, Yue Y, Hu K N, Li Y Y, Duan Z Q, Wu M, Tu J X, Shen J X, Yi B, Fu T D. Unravelling the complex trait of harvest index in rapeseed (L.) with association mapping., 2015, 16: 379–388.

[28] 姚東和, 楊民勝, 李志輝. 林分密度對巨尾桉生物產量及生產力的影響. 中南林業(yè)科技大學學報, 2000, 20(3): 20–23. Yao D H, Yang M S, Li Z H. The effect of stand density on the biological yield and productivity of., 2000, 20(3): 20–23(in Chinese with English abstract).

[29] 李躍建, 宋荷仙, 朱華忠. 小麥收獲指數、生物產量和籽粒產量的穩(wěn)定性分析. 西南農業(yè)學報, 1998, (1): 25–30. Li Y J, Song H X, Zhu H Z. Analysis of stability of wheat harvest index, biological yield and grain yield., 1998, (1): 25–30 (in Chinese with English abstract).

[30] Dai X L, Zhao J X, Xiang Y, Zhang T, Ren T B, Cheng G P. Correlation analysis of harvest index and plant traits of hybrid., 2018, 2: 47–58.

[31] Tian J, Ma K, Saaem I. Advancing high-throughput gene synthesis technology., 2009, 5: 714–722.

[32] Dillies M, Rau A, Aubert J, Hennequet C, Jeanmougin M, Servant A, Keime C, Marot G, Castel D, Estelle J, Guernec G, Jagla B, Jouneau L, Lalo? C, Gall L, Scha?ffer B, Le C, Guedj M, Jaffrézic F, Consortium. A comprehensive evaluation of normalization methods for Illumina high-throughput RNA sequencing data analysis., 2013, 14: 671–683.

[33] Wang J X, Zhang X Q, Shi M L, Gao L J, Niu X F, Te R, Chen L, Zhang W W. Metabolomic analysis of the salt-sensitive mutants reveals changes in amino acid and fatty acid composition important to long-term salt stress insp. PCC 6803., 2014, 14: 431–440.

[34] 秦天元, 孫超, 畢真真, 梁文君, 李鵬程, 張俊蓮, 白江平. 基于WGCNA的馬鈴薯根系抗旱相關共表達模塊鑒定和核心基因發(fā)掘. 作物學報, 2020, 46: 1033–1051. Qin T Y, Sun C, Bi Z Z, Liang W J, Li P C, Zhang J L, Bai J P. Identification of drought-related co-expression modules and hub genes in potato roots based on WGCNA., 2020, 46: 1033–1051 (in Chinese with English abstract).

[35] 宋長新, 雷萍, 王婷. 基于WGCNA算法的基因共表達網絡構建理論及其R軟件實現. 基因組學與應用生物學, 2013, 32: 135–141. Song C X, Lei P, Wang T. Gene co-expression network construction theory based on WGCNA algorithm and its R software implementation., 2013, 32: 135–141 (in Chinese with English abstract).

[36] Lu K, Wei L J, Li X L, Wang Y T, Wu J, Liu M, Zhang C, Chen Z Y, Xiao Z C, Jian H J, Cheng F, Zhang K, Du H, Cheng X C, Qu C M, Qian W, Liu L Z, Wang R, Zou Q Y, Ying J M, Xu X F, Mei J Q, Liang Y, Chai Y R, Tang Z L, Wan H F, Ni Y, He Y J, Lin N, Fan Y H, Sun W, Li N N, Zhou G, Zheng H K, Wang X W, Andrew H P, Li J N. Whole-genome resequencing revealsorigin and genetic loci involved in its improvement., 2019, 10: 1154–1165.

[37] Bradbury P J, Zhang Z, Kroon D E, Casstevens T M, Ramdoss Y, Buckler E S. TASSEL: software for association mapping of complex traits in diverse samples., 2007, 23: 2633–2635.

[38] Barrett J C, Fry B, Maller J, Daly M J. Haploview: analysis and visualization of LD and haplotype maps., 2005, 21: 263–265.

[39] Kohl M, Wiese S, Warscheid B. Cytoscape: software for visualization and analysis of biological networks., 2011, 696: 291–303.

[40] Livingston A K, Cruz J A, Kohzuma K, Dhingra A. Anmutant with high cyclic electron flow around photosystem I (hcef) involving the NADPH dehydrogenase complex., 2010, 22: 221–233.

[41] Mariam S. Increased sucrose level and altered nitrogen metabolism intransgenic plants expressing antisense chloroplastic fructose-1,6-bisphosphatase.,2004, 55: 2495–2503.

[42] Mull L, Ebbs M L, Bender J. A histone methylation-dependent DNA methylation pathway is uniquely impaired by deficiency inS-adenosylhomocysteine hydrolase., 2005, 174: 1161–1171.

[43] Ouyang B, Joung J G, Kolenovsky A, Koh C, Nowak J, Caplan A, Keller W A, Cui Y, Cutler A J, Tsang E W. Transcriptome profiling and methyl homeostasis of anmutant deficient in S-adenosylhomocysteine hydrolase1 (SAHH1)., 2012, 79: 315–331.

[44] Godge M R, Kumar D, Kumar P P.HOG1 gene and its petunia homolog PETCBP act as key regulators of yield parameters., 2008, 27: 1497–1507.

[45] Pedro S C F R, Mazhar S, Rosalba M, Mathilde F, Stéphanie B, Rebecca L, Barbara M, Conrad W, Hervé V, Ian F. TheH OMOLOGY-DEPENDENT GENE SILENCING1 gene codes for an-adenosyl--homocysteine hydrolase required for DNA methylation-dependent gene silencing., 2005, 17: 404–417.

[46] López C, Abuzaid O, Lawson T, Raines C A.CP12 mutants have reduced levels of phosphoribulokinase and impaired function of the calvin-benson cycle., 2017, 68: 2285–2298.

[47] Fermani S, Trivelli X, Sparla F, Thumiger A, Calvaresi M, Marri L, Falini G, Zerbetto F, Trost P. Conformational selection and folding-upon-binding of intrinsically disordered protein CP12 regulate photosynthetic enzymes assembly., 2012, 287: 21372–21383.

[48] Singh P, Kaloudas D, Raines C A. Expression analysis of theCP12 gene family suggests novel roles for these proteins in roots and floral tissues.,2008, 59: 3975–3985.

[49] Marri L, Thieulin P G, Lebrun R, Puppo R, Zaffagnini M, Trost P, Gontero B, Sparla F. CP12-mediated protection of Calvin-Benson cycle enzymes from oxidative stress., 2014, 97: 228–237.

[50] Ferro M, Salvi D, Seigneurin B D, Court M, Moyet L, Ramus C, Miras S, Mellal M, Le G S, KiefferJ S, Bruley C, Garin J, Joyard J, Masselon C, Rolland N. AT_CHLORO, a comprehensive chloroplast proteome database with subplastidial localization and curated information on envelope proteins.,2010, 9: 1063–1084.

[51] Henry E, Fung J, Liu G, Drakakak G, Coaker G. Beyond glycolysis: GAPDHs are multi-functional enzymes involved in regulation of ROS, autophagy, and plant immune responses., 2015, 11: e1005199.

[52] Ringli C, Baumberger N, Diet A, Frey B, Keller B. ACTIN2 is essential for bulge site selection and tip growth during root hair development of., 2002, 129: 1464–1472.

[53] Lu K, Liu P, Zhang C, Lu J H, Yang B, Xiao Z C, Liang Y, Xu X F, Qu C M, Zhang K, Liu L Z, Zhu Q L, Fu M L, Yuan X Y, Li J N. Genome-wide association and transcriptome analyses reveal candidate genes underlying yield-determining traits in., 2017, 8: 206–220.

[54] Lu K, Xiao Z C, Jian H J, Liu P, Qu C M, Fu M L, He B, Tie L M, Liang Y, Xu X F, Li J N. A combination of genome-wide association and transcriptome analysis reveals candidate genes contro-lling harvest index-related traits in., 2016, 6: 36452.

[55] Sun Y Y, Liu T T, Yang H Y, Zuo Q S, Zhou G S, Wu J S. The correlation between the growth characteristics of rapeseed stalk and its yield and the lodging resistance., 2014, 53: 30–35.

[56] 柳寒. 油菜灌漿期種子和角果皮基因表達差異及相關基因的功能分析. 華中農業(yè)大學博士學位論文, 湖北武漢, 2015.Liu H. Gene Expression and Functional Analysis in Seeds and Siliques of Rapeseed during the Filling Stage. PhD Dissertation of Huazhong Agricultural University, Wuhan, Hubei, China, 2015 (in Chinese with English abstract).

附圖1 各樣本中關鍵模塊的所有基因與相應ME的表達水平

Fig. S1 Expression levels of all genes and corresponding ME in key modules of each sample

Integrating GWAS and WGCNA to screen and identify candidate genes for biological yield inL.

WANG Yan-Hua1,2,**, LIU Jing-Sen1,2,**, and LI Jia-Na1,2,*

1College of Agronomy and Biotechnology, Southwest University / Chongqing Engineering Research Center for Rapeseed, Chongqing 400715, China;2Academy of Agricultural Sciences, Southwest University, Chongqing 400715, China

Biomass yield is especially important for, as it is the basis for high yields of crops. In this study, the phenotypic data of the natural populations composed of 588 materials were used for genome-wide association analysis (GWAS). We performed the transcriptome sequencing (RNA-seq) of biomass yield using ‘CQ45’ (high biological yield material) and ‘CQ46’ (low biological yield material). A weighted gene co-expression network analysis (WGCNA) network was constructed by integrating transcriptome data of six tissues of the extreme materials, such as stalks, leaves, 30 day after flowering (DAF) seeds of main inflorescence and lateral branch, 30 DAF pod keratin of main branch and lateral branch. We finally screened the candidate genes related to biomass yield. The main results are as follows: Biomass yields inhad positive effects on most yield-related traits; K + PCA model was the best model for biomass analysis of the natural population, and nine significant loci were detected in the best model (< 1/385691 or< 0.05/385691); according to 36 groups of transcriptome data, MAD value of each gene was calculated. A total of 5052 genes with MAD value of the top 5% were selected to construct WGCNA. Fifteen gene modules were obtained, among which, five genes co-expression modules were significantly correlated with leaves, stems, and seeds of 30 DAF. The hub genes of the key modules in WGCNA, the significant SNP loci obtained from GWAS, and the extreme phenotypic differential genes were integrated to identify the candidate genes. Theirhomologous genes were,,, and, which played the important roles in the Calvin cycle, carbon assimilation, and material accumulation of photosynthesis.

GWAS; WGCNA; RNA seq;

10.3724/SP.J.1006.2021.04175

本研究由中國博士后科學基金面上項目(2019M653319), 重慶市自然科學基金博士后科學基金項目(cstc2019jcyj-bshXO116)和高等學校學科創(chuàng)新引智基地項目(“111”項目)(B12006)資助。

This study was supported by the Project of China Postdoctoral Science Foundation (2019M653319), the Project of Chongqing Natural Science Foundation Postdoctoral Science Foundation (cstc2019jcyj-bshXO116), and the Project of Intellectual Base for Discipline Innovation in Colleges and Universities (“111” Project) (B12006).

李加納, E-mail: ljn1950@swu.edu.cn

**同等貢獻(Contributed equally to this work)

E-mail: hawer313@163.com

2020-07-30;

2020-12-01;

2021-01-11.

URL: https://kns.cnki.net/kcms/detail/11.1809.S.20210108.1700.010.html