甜玉米乳熟期果皮厚度相關(guān)基因的加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析

裴 虎,熊才運(yùn),張亞輝,任文闖,李小琴,黃 君

(廣東省植物分子育種重點(diǎn)實(shí)驗(yàn)室,華南農(nóng)業(yè)大學(xué) 農(nóng)學(xué)院,廣東 廣州 510642)

甜玉米(ZeamaysL.sacharatasturt)是玉米屬(ZeamaysL.)甜質(zhì)型玉米亞種,是由su1、sh1、bt1、sh2等隱性突變基因引起的胚乳缺陷型玉米[1]。果皮是玉米籽粒的重要組成部分,具有保護(hù)胚和胚乳不受機(jī)械損傷或病蟲害侵害的作用[2]。與飼用玉米不同,食味品質(zhì)是決定甜玉米品種優(yōu)劣和經(jīng)濟(jì)價(jià)值的重要因素,果皮厚度與甜玉米的柔嫩性和皮渣感有關(guān),是影響其食味品質(zhì)最為重要的關(guān)鍵指標(biāo)之一[3]。目前,改良甜玉米的食味品質(zhì),尤其是果皮柔嫩度,選育果皮較薄的甜玉米品種,提升甜玉米口感,是育種家一直追求的重要育種目標(biāo)。因此,深入挖掘控制甜玉米果皮厚度關(guān)鍵基因,能夠?yàn)樘鹩衩灼焚|(zhì)改良提供基因資源,具有重要的研究意義。

截至目前,針對(duì)甜玉米果皮厚度的研究較少,主要集中在QTL(Quantitative trait locus)定位和遺傳方式分析。玉米果皮厚度的遺傳方式涉及顯性效應(yīng)、加性效應(yīng)和上位性效應(yīng),其中加性效應(yīng)和上位性效應(yīng)起到主要作用[4-6]。前人研究表明,在玉米的10條染色體中除第7號(hào)染色體外,其余9條染色體均定位到玉米果皮厚度相關(guān)QTLs[5,7-8],但多數(shù)QTL難以在不同研究中被重復(fù)檢測到。2020年,華南農(nóng)業(yè)大學(xué)農(nóng)學(xué)院甜玉米課題組利用甜玉米高代導(dǎo)入系群體對(duì)甜玉米果皮厚度進(jìn)行QTL定位和轉(zhuǎn)錄組分析,在第10染色體上鑒定到一個(gè)控制果皮厚度的主效QTL-qPT10-5,能夠解釋表型變異的7.78%～35.38%,結(jié)合轉(zhuǎn)錄組測序和基因注釋分析,認(rèn)為AUX/IAA 轉(zhuǎn)錄因子,ZIM 轉(zhuǎn)錄因子和FATTY ACID EXPORT是控制甜玉米果皮厚度的關(guān)鍵候選基因[9]。

目前,通過正向遺傳學(xué)方法已鑒定出多個(gè)與甜玉米果皮厚度相關(guān)的QTLs,但至今未有相關(guān)基因被克隆的報(bào)道。加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析(Weighted gene co-expression network analysis,WGCNA)是以轉(zhuǎn)錄組、基因芯片等獲得的基因表達(dá)量數(shù)據(jù)為基礎(chǔ),將表達(dá)趨勢相近的基因模塊分類并得到具有高度生物學(xué)意義的模塊和基因[10]。近年來,WGCNA方法已經(jīng)廣泛應(yīng)用于玉米、水稻、高粱等農(nóng)作物的遺傳改良研究[11-14]。然而目前尚未有利用WGCNA對(duì)甜玉米果皮厚度相關(guān)的基因共表達(dá)模塊被鑒定。

本研究以薄果皮甜玉米自交系M03和厚果皮甜玉米自交系M08為材料,在3個(gè)時(shí)間梯度(授粉后15,19,23 d)分析其果皮厚度變化,結(jié)合轉(zhuǎn)錄組和WGCNA分析構(gòu)建果皮厚度特異性模塊并挖掘出調(diào)控玉米果皮厚度的核心基因,為后續(xù)甜玉米果皮厚度分子遺傳機(jī)制的進(jìn)一步研究提供新的線索。

1 材料和方法

1.1 玉米材料及轉(zhuǎn)錄組測序數(shù)據(jù)

利用本課題組選育的2個(gè)甜玉米自交系M03(果皮厚度71.1～135.9 μm)和M08(果皮厚度111.7～171.7 μm)為研究材料。于2017年種植在華南農(nóng)業(yè)大學(xué)增城教學(xué)實(shí)驗(yàn)基地,行距70 cm、株距25 cm、種植密度為52 500株/hm2。為避免材料出現(xiàn)其他品種玉米植株造成的花粉直感現(xiàn)象,采取人工套袋授粉的方式進(jìn)行授粉,分別在授粉后15,19,23 d對(duì)2個(gè)自交系進(jìn)行取樣。每個(gè)自交系隨機(jī)選取3個(gè)玉米果穗作為3個(gè)生物學(xué)重復(fù),取10～15粒果穗中部的籽粒,用雙刃刀片沿胚側(cè)輕輕劃開籽粒,用鑷子剝?nèi)」げ⒘⒓捶湃胍旱鋬?最后放入-80 ℃冰箱保存,用于后續(xù)表型鑒定及轉(zhuǎn)錄組測序。每個(gè)果穗取6粒籽粒進(jìn)行掃描電鏡觀察測量果皮厚度。使用SPSS(Version 19)進(jìn)行表型數(shù)據(jù)分析。

植物總 RNA 利用RNA提取試劑盒(Vazyme,中國南京)進(jìn)行提取。提取的 RNA使用Nanodrop 2000分光光度計(jì)(Thermofisher,美國)、Qubit 2.0熒光計(jì)(Life Technologies)和Aglient 2100生物分析儀系統(tǒng)(Agilent Technologies,CA,美國)來評(píng)估純度(OD260/280≥1.8;OD260/230≥1.0),濃縮RNA樣品(總RNA≥250 ng/μL)和完整性(RIN≥8.0,28S/18S≥0.5)。轉(zhuǎn)錄組測序文庫構(gòu)建主要使用寡核苷酸(dT)磁珠(Invitrogen,Carlsbadcity,CA,美國)從總RNA中分離出每個(gè)樣品1.5 μg mRNA,片段化并反向轉(zhuǎn)錄成cDNA,并通過 PCR進(jìn)行擴(kuò)增。使用Illumina HiSeq 2000平臺(tái)(San Diego,CA,美國)進(jìn)行了轉(zhuǎn)錄組測序。使用 FastQC軟件對(duì)原始測序數(shù)據(jù)進(jìn)行質(zhì)量評(píng)估,去除低質(zhì)量和未知堿基數(shù)目超過15%的 reads。利用 Hisat 2 軟件將Clean reads 比對(duì)到玉米B73參考基因組(ftp://ftp.gramene.org/pub/gramene/release63/fasta/zea_mays/Zea_mays.B73_RefGen_v4.dna.top_level.fa.gz),利用Cufflinks v2.2.1軟件計(jì)算樣品中基因的表達(dá)量,每千個(gè)堿基外顯子百萬片段數(shù)(Fragments per kilobase of exon per million fragments mapped,FPKM)值衡量基因的表達(dá)水平。DESeq2 軟件(Version 1.10.1)用于差異表達(dá)基因的計(jì)算[15]。只有當(dāng)不同組間|Log2Fold Change|>1 且P值<0.05的基因被認(rèn)定為差異表達(dá)基因。

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

通過TBtools v1.09[16]軟件中的WGCNA shinny plugin對(duì)共表達(dá)模塊進(jìn)行劃分。去除掉FPKM值小于1的基因,本研究共得到14 126個(gè)基因用于樣品聚類。通過SFT and Power Select 選項(xiàng)進(jìn)行軟閾值(Soft Thresholding Power)的篩選,計(jì)算得到閾值為14(圖1)。利用Module-net中的Cluster選項(xiàng)將min Module Size調(diào)整為30,Module Cuttree Height調(diào)整為0.25,選擇Max Block Size為15 000對(duì)數(shù)據(jù)進(jìn)行模塊聚類和剪切。利用Module-trait選項(xiàng)計(jì)算基因在模塊內(nèi)的連接度(Eigengene connectivity,KME),取KME值排名靠前且差異表達(dá)的基因作為該模塊內(nèi)核心基因(Hub gene)進(jìn)行后續(xù)分析,利用Cytoscape v3.6.0[17]對(duì)核心基因網(wǎng)絡(luò)進(jìn)行可視化。

B縱坐標(biāo)代表每一個(gè)軟閾值對(duì)應(yīng)的網(wǎng)絡(luò)平均連接程度;C、D.在power值為14時(shí)拓?fù)渚W(wǎng)絡(luò)符合程度的檢測。The vertical coordinate of Fig.B represents the mean connectivity corresponding to each soft threshold;C and D .The detection of topological network compliance at a power value of 14.圖1 基因共表達(dá)網(wǎng)絡(luò)軟閾值的確定Fig.1 Soft threshold determination of gene co-expression network

1.3 差異表達(dá)基因及特異性模塊功能富集分析

將篩選出的差異表達(dá)基因和特異性模塊內(nèi)所包含的基因利用R程序中的 clusterProfiler 軟件包[18]進(jìn)行GO(Gene Ontology)和 KEGG(Kyoto Encyclopedia of Genes and Genomes)富集分析。當(dāng)P值(P-value)小于0.05時(shí)認(rèn)為此 GO功能和 KEGG pathway 功能存在顯著富集情況。使用 Tbtools 軟件和Omicshare平臺(tái)制作熱圖對(duì)基因表達(dá)量進(jìn)行可視化。

1.4 轉(zhuǎn)錄組測序數(shù)據(jù)的qRT-PCR驗(yàn)證

利用實(shí)時(shí)熒光定量PCR(Quantitative real-time PCR,qRT-PCR)驗(yàn)證RNA-seq結(jié)果的可靠性。使用HiScriptⅢ RT SuperMix for qPCR(+gDNA wiper)試劑盒將果皮RNA反轉(zhuǎn)錄成cDNA作為qRT-PCR模板。以u(píng)biquitin(登錄號(hào):NM_001138130)為內(nèi)參基因,使用CFX96(Bio-Rad)實(shí)時(shí)熒光定量 PCR 系統(tǒng),選用2×T5 Fast qPCR Mix(SYBR GreenI)(TsingKe)試劑盒,熒光染料為SYBR Green。所有反應(yīng)均在20 μL的體系下進(jìn)行,其中包含200 ng cDNA、0.8 μmol/L核心基因特異性引物以及10 μL的2×T5 Fast qPCR Mix。qPCR反應(yīng)程序:95 ℃預(yù)變性1 min;95 ℃ 10 s,55 ℃ 5 s,40個(gè)循環(huán);溶解曲線95 ℃ 10 s,55 ℃ 10 s,72 ℃ 15 s。使用2-ΔΔCT法分析基因相對(duì)表達(dá)量[19],所用基因引物見表1。

2 結(jié)果與分析

2.1 不同材料果皮厚度變化分析

在相同的栽培條件下,分別取自交系M03和M08授粉后15,19,23 d的果皮進(jìn)行掃描電鏡觀察(Scanning electron microscope,SEM)(圖2-A—F)。研究結(jié)果顯示,隨著授粉時(shí)間的增加,甜玉米籽粒果皮厚度整體呈現(xiàn)逐漸降低的趨勢,授粉后15～19 d的果皮厚度降低的速度較為緩慢,從授粉后19 d開始果皮厚度呈現(xiàn)出快速變薄的趨勢。其中,M08的果皮厚度相比M03更厚,不同自交系間果皮厚度的差異隨時(shí)間的推移逐漸增加,在乳熟期(授粉后19 d)開始達(dá)到顯著水平(圖2-G)。

A—C.M03授粉后15,19,23 d的果皮掃描電鏡結(jié)果;D—F.M08授粉后15,19,23 d的掃描電鏡結(jié)果;G.M03和M08的果皮厚度變化,柱中字母為 Duncan′s多重比較結(jié)果,不同字母表示材料間在P<0.05水平差異顯著。A—C.SEM results of 15,19 and 23 days after M03 pollination;D—F.SEM results of 15,19 and 23 days after M08 pollination;G.The change of pericarp thickness of M03 and M08.The letters in the column are Duncan′s multiple comparison results,different letters indicate significant difference between materials at the level of P<0.05.圖2 M08和M03在3個(gè)不同時(shí)期的籽粒果皮掃描電鏡圖片以及果皮厚度變化Fig.2 SEM of seed pericarp and changes in pericarp thickness of M08 and M03 at three different periods

2.2 差異表達(dá)基因篩選及GO和KEEG富集分析

2.2.1 轉(zhuǎn)錄組數(shù)據(jù)分析 去除含有接頭污染和低質(zhì)量的reads,共得到480 628 079條高質(zhì)量reads,與玉米B73參考基因組進(jìn)行序列比對(duì),其比對(duì)率為78.82%～ 84.50%,外顯子比對(duì)率為87.38%～91.97%(表2)。

表2 參試樣品轉(zhuǎn)錄組數(shù)據(jù)質(zhì)量Tab.2 Quality of the transcriptome data of each sample

通過M03和M08自交系內(nèi)部不同發(fā)育時(shí)期以及不同發(fā)育時(shí)期材料間的比較,共篩選出4 748個(gè)DEGs。其中,M03材料內(nèi)部比較(M03-DAP15 vs M03-DAP23)共鑒定到1 234個(gè)DEGs(上調(diào)641個(gè),下調(diào)593個(gè));M08材料在不同發(fā)育時(shí)期的果皮厚度差異相對(duì)于M03較少,其內(nèi)部比較(M08-DAP15 vs M08-DAP23)獲得了較少的差異基因,為471個(gè)DEGs(上調(diào)255個(gè),下調(diào)216個(gè))。相同發(fā)育時(shí)期的材料間比較(M03-DAP15 vs M08-DAP15;M03-DAP19 vs M08-DAP19;M03-DAP23 vs M08-DAP23)獲得了更多的差異基因,分別為2 334,2 238和2 039個(gè)DEGs(圖3-A)。材料內(nèi)部與材料間比較共有的差異基因?yàn)?1個(gè)(圖3-B)。

圖3 不同材料不同時(shí)期DEGs數(shù)量統(tǒng)計(jì)及韋恩圖Fig.3 Number of DEGs in different datasets and venn diagram

2.2.2 差異基因富集分析 對(duì)同一材料不同發(fā)育時(shí)期及不同材料相同發(fā)育時(shí)期進(jìn)行差異表達(dá)基因鑒定,對(duì)獲得的差異表達(dá)基因進(jìn)行GO功能注釋,GO富集條目主要分為生物過程、分子功能及細(xì)胞組分3 類。結(jié)果顯示,差異表達(dá)基因GO功能富集主要富集到細(xì)胞組分組織(GO:0071840)、對(duì)刺激的反應(yīng)(GO:0050896)、生物過程調(diào)節(jié)(GO:0050789)等生物過程;結(jié)構(gòu)分子活性(GO:0005198)、核酸結(jié)合轉(zhuǎn)錄物(GO:0001071)、轉(zhuǎn)運(yùn)活性(GO:0005215)等分子功能;膜(GO:0016020)、細(xì)胞組分(GO:0044464)、細(xì)胞(GO:0005623)等細(xì)胞組分(圖4)。

KEGG代謝通路富集分析結(jié)果顯示差異表達(dá)基因主要參與丙氨酸、天冬氨酸和谷氨酸代謝(Alanine,aspartate and glutamate metabolism)、萜類骨架的生物合成(Terpenoid backbone biosynthesis)和植物MAPK信號(hào)通路(MAPK signaling pathway-plant);M03與M08材料內(nèi)部的不同發(fā)育時(shí)期比較(M03-DAP15 vs M03-DAP23與M08-DAP15 vs M08-DAP23)結(jié)果顯示,差異表達(dá)基因還可共同富集到氨基酸生物合成(Biosynthesis of amino acids)、淀粉和蔗糖代謝(Starch and sucrose metabolism);相同時(shí)期不同材料之間比較得出差異表達(dá)基因主要參與β-丙氨酸代謝(Beta-alanine metabolism)、精氨酸和脯氨酸代謝(Arginine and proline metabolism)、半胱氨酸和蛋氨酸代謝(Cysteine and methionine metabolism)途徑(圖4)。

2.3 共表達(dá)基因聚類和模塊構(gòu)建

通過對(duì)FPKM值小于1的基因進(jìn)行過濾,最終獲得14 126個(gè)基因用于構(gòu)建基因聚類樹和共表達(dá)模塊分析。根據(jù)基因的FPKM值將表達(dá)模式相似的基因進(jìn)行聚類,利用動(dòng)態(tài)切割法(Dynamic tree cut)將得到的基因聚類樹切割成不同模塊(圖5-A)。將切割得到的18個(gè)不同模塊分別用不同的顏色進(jìn)行區(qū)分,其中,基因數(shù)目最多的模塊是Turquoise模塊,最少的是Grey60模塊(圖5-B)。

A.基因聚類樹,每個(gè)聚類樹對(duì)應(yīng)一個(gè)基因模塊,聚類等級(jí)相同的模塊用同樣的顏色表示;B.聚類出的18個(gè)模塊以及每個(gè)模塊所包含基因個(gè)數(shù)。A.A gene cluster tree,each cluster tree corresponds to a gene module,and with the same clustering level are represented by the same color;B.The 18 modules clustered and the number of genes contained in each module.圖5 基因共表達(dá)模塊以及各模塊所包含的基因數(shù)Fig.5 Gene co-expression module and the number of genes contained in each module

2.4 特異性模塊篩選及目標(biāo)模塊KEGG富集分析

由于同一模塊內(nèi)的基因表達(dá)模式相近,因此,這些基因很可能具有相似的生物學(xué)功能。結(jié)合果皮厚度表型數(shù)據(jù)計(jì)算每個(gè)模塊與果皮厚度之間的相關(guān)性,最終在18個(gè)模塊中得到了4個(gè)相關(guān)性系數(shù)絕對(duì)值大于0.60的模塊,分別是Pink模塊(-0.63)、Yellow模塊(-0.82)、Turquoise模塊(0.66)和Magenta模塊(0.81)(圖6)。通過模塊之間的相關(guān)性熱圖可以看出4個(gè)模塊之間具有較高的相關(guān)性,可將其作為目標(biāo)模塊進(jìn)行后續(xù)分析。為進(jìn)一步探究與甜玉米果皮厚度顯著相關(guān)模塊的生物學(xué)功能,本研究利用OmicShare Tools(https://www.omicshare.com/)對(duì)模塊進(jìn)行KEGG功能富集分析。結(jié)果表明,除Magenta模塊外,其他3個(gè)模塊均富集到氨基酸的生物合成、α-亞麻酸代謝(Alpha-linolenic acid metabolism)、半胱氨酸和蛋氨酸代謝等通路。此外,泛酸和輔酶A生物合成(Terpenoid backbone biosynthesis)在Turquoise模塊被特異性富集;脂肪酸延長(Fatty acid elongation)在Yellow模塊特異富集;泛酸和輔酶A的生物合成(Pantothenate and CoA biosynthesis)在Pink模塊特異富集;植物MAPK信號(hào)通路、植物病原體相互作用等代謝通路只在Magenta模塊特異富集(圖7)。將不同材料及不同發(fā)育階段鑒定到的DEGs與所有所有模塊相關(guān)基因的KEGG富集結(jié)果進(jìn)行對(duì)比,結(jié)果顯示,不同方法得到的基因可以共同富集到丙氨酸、天冬氨酸和谷氨酸代謝、植物MAPK信號(hào)通路、半胱氨酸和蛋氨酸代謝等通路。

A.各模塊與果皮厚度之間的相關(guān)系數(shù)熱圖,紅色框?yàn)樘禺愋阅K;B.特異性模塊與DEGs共有基因數(shù)。A.The heat map of correlation coefficient between each module and peel thickness,and the red box is the specific module; B.The number of genes shared by specific modules and DEGs.圖6 表型-模塊相關(guān)性熱圖和模塊-DEGs共有基因數(shù)Fig.6 Phenotype-module correlation heatmap and number of genes shared by significant modules and DEGs

圖7 特異性模塊KEGG富集分析Fig.7 KEGG enrichment analysis of significant module

2.5 核心基因挖掘及互作網(wǎng)絡(luò)構(gòu)建與分析

將不同模塊內(nèi)的基因與不同材料、不同發(fā)育時(shí)期鑒定到的差異表達(dá)基因進(jìn)行比對(duì),選擇二者的交集進(jìn)行基因互作網(wǎng)絡(luò)構(gòu)建。通過篩選,Turquoise模塊中有1 351個(gè)基因、Yellow模塊中有543個(gè)基因、Pink模塊中有10個(gè)基因、Magenta模塊中有34個(gè)基因用于互作網(wǎng)絡(luò)構(gòu)建。將每個(gè)模塊中與表型相關(guān)性較低(|GS|<0.5)的基因進(jìn)行過濾,選取模塊內(nèi)連通性(KME)最高的前20個(gè)基因作為核心基因的候選基因。為獲取候選基因的功能注釋信息,通過玉米參考基因組數(shù)據(jù)庫MaizeGDB(https://www.maizegdb.org/)和擬南芥同源注釋對(duì)相關(guān)基因的功能進(jìn)行查找,最終篩選出與甜玉米果皮厚度相關(guān)的核心基因(Hub gene)。在Turquoise模塊共篩選出3個(gè)核心基因,Yellow模塊篩共選出4個(gè)核心基因,Pink模塊共篩選出2個(gè)核心基因,Magenta模塊共篩選出4個(gè)核心基因(表3)。將核心基因作為互作網(wǎng)絡(luò)的核心節(jié)點(diǎn)構(gòu)建互作網(wǎng)絡(luò),篩選基因網(wǎng)絡(luò)中的轉(zhuǎn)錄因子(Transcript factor,TF)和與果皮厚度相關(guān)的其他基因,并選擇權(quán)重值(Weight)最高的前50個(gè)基因繪圖,利用Cytoscape v3.6.0軟件進(jìn)行可視化(圖8)。

表3 模塊中核心基因在擬南芥中的同源基因功能注釋Tab.3 Functional notes of homologous genes of hub genes in module in Arabidopsis

A～D分別為Turquoise模塊、Yellow模塊、Pink模塊、Magenta模塊的基因網(wǎng)絡(luò);紅色大點(diǎn)為核心基因,灰色矩形為模塊內(nèi)與表型相關(guān)基因的標(biāo)注。A—D are the gene networks of Turquoise module,Yellow module,Pink module and Magenta module respectively; the large red dot is the hub gene,and the gray rectangle is the annotation of phenotype related genes in the module.圖8 特異性模塊的基因共表達(dá)網(wǎng)絡(luò)及其核心基因Fig.8 Gene co-expression network and hub genes of significant modules

通過基因功能注釋發(fā)現(xiàn)核心基因與已報(bào)道的可能與果皮厚度發(fā)育相關(guān)的基因具有互作網(wǎng)絡(luò)關(guān)系。例如:Turquoise模塊中的核心基因與細(xì)胞分裂素N-葡萄糖基轉(zhuǎn)移酶ZmCNGT(Zm00001d019265),GA2氧化酶GA2ox9(Zm00001d008909)和果膠甲酯酶ZmPME(Zm00001d050082)等具有互作網(wǎng)絡(luò)關(guān)系;Yellow模塊中的核心基因與聚半乳糖醛酸酶PGL23(Zm00001d013032)具有互作網(wǎng)絡(luò)關(guān)系。這些核心基因互作網(wǎng)絡(luò)中還包含ARF、AUX/IAA、bZIP、MYB、ERF、WRKY等可能與細(xì)胞分裂與果皮厚度發(fā)育相關(guān)的轉(zhuǎn)錄因子。

2.6 核心基因的實(shí)時(shí)熒光定量PCR驗(yàn)證

為驗(yàn)證轉(zhuǎn)錄組測序結(jié)果中基因表達(dá)量的可靠性,本研究挑選了9個(gè)核心基因進(jìn)行qRT-PCR表達(dá)量驗(yàn)證。結(jié)果表明,這些基因在2個(gè)不同的甜玉米自交系授粉后15,19,23 d的果皮組織中的表達(dá)量水平與轉(zhuǎn)錄組測序結(jié)果一致(圖9),進(jìn)一步說明了轉(zhuǎn)錄組測序結(jié)果的可靠性以及數(shù)據(jù)的可重復(fù)性。

3 結(jié)論與討論

甜玉米作為鮮食玉米的一種,果皮厚度是決定其口感優(yōu)劣的重要因素[20]。果皮角質(zhì)化程度、細(xì)胞大小、細(xì)胞數(shù)量細(xì)胞壁厚度以及細(xì)胞的排布等因素都可能會(huì)影響果皮的厚度[21]。這些因素除受基因型影響外還與栽培環(huán)境、氣候變化、采收時(shí)間等有關(guān)。因此,本研究在同等的栽培條件下,嚴(yán)格控制了授粉時(shí)間和采收時(shí)間,最大程度避免了上述因素對(duì)甜玉米果皮厚度表型的影響。為解析甜玉米果皮厚度的遺傳基因,鑒定控制甜玉米果皮厚度的關(guān)鍵候選基因,本研究選擇果皮厚度有明顯差異的2個(gè)甜玉米自交系(M03和M08)進(jìn)行轉(zhuǎn)錄組測序,通過基因差異表達(dá)分析和加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析,選取二者的交集進(jìn)行后續(xù)研究。WGCNA作為識(shí)別基因與樣本間相關(guān)性模式的常用方法,可以特異性的篩選出與目標(biāo)性狀相關(guān)的基因并聚類成共表達(dá)模塊,是相關(guān)性狀協(xié)同表達(dá)調(diào)控控制解析的重要工具[22]。

3.1 差異表達(dá)基因和特異性模塊內(nèi)基因與碳水化合物代謝、淀粉和蔗糖代謝及MAPK信號(hào)通路有關(guān)

本研究利用乳熟期果皮較薄的甜玉米自交系M03和果皮較厚自交系M08的乳熟期果皮材料的RNA-seq數(shù)據(jù)為基礎(chǔ)進(jìn)行了差異表達(dá)分析,在相同材料的不同發(fā)育時(shí)期及不同材料的相同發(fā)育時(shí)期共鑒定了4 748個(gè)差異表達(dá)基因,對(duì)相關(guān)基因進(jìn)行GO和KEGG功能富集分析發(fā)現(xiàn)“碳水化合物代謝過程”、 “淀粉和蔗糖代謝”等與甜玉米果皮厚度發(fā)育有關(guān)。在禾本科植物研究中,如小麥、高粱等作物,發(fā)現(xiàn)其籽粒果皮中含有較多的淀粉顆粒,果皮中淀粉的含量會(huì)隨著籽粒的發(fā)育而變化[23-24]。在小麥果皮厚度發(fā)育研究中發(fā)現(xiàn)小麥植株受精后珠被開始發(fā)育為果皮,淀粉隨著果皮的發(fā)育不斷積累,在DAP(Day after pollination)5左右果皮中淀粉含量到達(dá)峰值,隨后逐漸分解并在DAP22左右消失殆盡,在果皮細(xì)胞淀粉含量處于峰值時(shí)可以明顯觀察到細(xì)胞呈現(xiàn)出輕微的膨脹變形,說明果皮細(xì)胞中的淀粉含量可以通過影響細(xì)胞大小進(jìn)而影響果皮厚度,不同的淀粉含量可能會(huì)造成果皮厚度的差異,與本研究結(jié)果一致,甜玉米果皮厚度發(fā)育的關(guān)鍵基因均與碳水化合物代謝、淀粉和蔗糖代謝相關(guān)[24]。

通過WGCNA分析,本研究共鑒定到了4個(gè)具有高度生物學(xué)意義并與果皮厚度顯著相關(guān)的模塊,通過富集分析揭示了這些模塊的生物學(xué)功能。將所有模塊的KEGG富集結(jié)果與差異表達(dá)基因的功能富集結(jié)果進(jìn)行比較,鑒定到天冬氨酸、丙氨酸與谷氨酸代謝、植物MAPK信號(hào)通路、半胱氨酸和蛋氨酸代謝等通路與甜玉米果皮厚度的發(fā)育有關(guān)。其中MAPK信號(hào)通路涉及包含細(xì)胞防御、抗逆、細(xì)胞增殖和細(xì)胞凋亡在內(nèi)的多種生命活動(dòng),推測該通路在調(diào)節(jié)甜玉米果皮厚度發(fā)育過程中具有十分重要的作用。

3.2 MYB轉(zhuǎn)錄因子、脫落酸、細(xì)胞壁pH值及細(xì)胞周期蛋白能夠影響甜玉米果皮厚度

在4個(gè)特異性模塊中選擇與差異表達(dá)基因的交集進(jìn)行后續(xù)分析,共篩選到13個(gè)核心基因,這些核心基因可能在參與調(diào)控果皮厚度發(fā)育過程中具有重要作用。Turquoise模塊的Zm00001d023282和Yellow模塊的Zm00001d022227分別為MYB69和MYB70轉(zhuǎn)錄因子,編碼MYB結(jié)構(gòu)域蛋白。前人在玉米MYB轉(zhuǎn)錄因子的研究中發(fā)現(xiàn),ZmFDL1/MYB94功能的缺失導(dǎo)致ω-羥基脂肪酸和多羥基脂肪酸的合成顯著降低,而角質(zhì)層又是由這些超長鏈脂肪酸(VLCFA)化合物組成,該基因的功能缺失會(huì)造成玉米發(fā)育前期包括籽粒、幼苗葉片等部位的角質(zhì)層形成受阻,同時(shí)該研究還指出干旱脅迫介導(dǎo)ABA的產(chǎn)生會(huì)抑制ZmFDL1的表達(dá),從而增加角質(zhì)層的厚度防止水分過快蒸發(fā)[25]。在擬南芥研究中發(fā)現(xiàn)幾個(gè)與表皮蠟質(zhì)生物合成的基因受MYB轉(zhuǎn)錄因子家族成員調(diào)控,包括AtMYB30、AtMYB94和AtMYB96[26-28]。本研究鑒定到的核心基因Zm00001d005884為bZIP轉(zhuǎn)錄因子基因,在擬南芥的同源基因?yàn)锳T1G49720,該基因與脫落酸反應(yīng)元件結(jié)合,介導(dǎo)ABA反應(yīng)的轉(zhuǎn)錄調(diào)控。核心基因Zm00001d014291在擬南芥的同源基因?yàn)锳T4G33500,該基因編碼蛋白磷酸酶2C家族蛋白,在番茄中的研究表明,ABA與受體結(jié)合會(huì)影響編碼蛋白磷酸酶蛋白的基因SIPP2C5的表達(dá),該基因表達(dá)量的變化會(huì)影響果實(shí)的各種性狀,包括果形、果皮厚度、可溶性固形物含量等[29]。核心基因Zm00001d012460編碼BAG家族蛋白,該蛋白與細(xì)胞凋亡調(diào)節(jié)有關(guān),Earp等[23]通過對(duì)禾本科作物的研究發(fā)現(xiàn),除了高粱以外其他大部分禾本科作物隨著籽粒的成熟,果皮細(xì)胞將經(jīng)歷PCD過程,使果皮呈現(xiàn)出厚度逐漸降低的趨勢。Magenta模塊中的核心基因Zm00001d041214編碼種質(zhì)膜H+-ATP酶,該基因的突變會(huì)導(dǎo)致擬南芥植株發(fā)育受限,生長素誘導(dǎo)的信號(hào)級(jí)聯(lián)反應(yīng)中,最終作用到質(zhì)膜H+-ATP 酶,被生長素信號(hào)激活增強(qiáng)的質(zhì)子泵將H+送進(jìn)細(xì)胞壁降低pH值,激活細(xì)胞壁內(nèi)對(duì)pH值敏感的酶和蛋白質(zhì),造成細(xì)胞壁松動(dòng)、延伸和生長[30]。核心基因Zm00001d048143編碼一種細(xì)胞周期蛋白,細(xì)胞周期蛋白依賴性激酶(CDK)是一個(gè)重要的絲氨酸/蘇氨酸蛋白激酶家族,在細(xì)胞有絲分裂中起著重要的作用[31-32]。前人在番茄果皮發(fā)育的研究中發(fā)現(xiàn),過表達(dá)或抑制CDK蛋白相關(guān)基因的表達(dá)會(huì)顯著影響番茄果皮的細(xì)胞層數(shù)和厚度[33]。

3.3 特異性模塊內(nèi)與果皮厚度相關(guān)基因挖掘

甜玉米乳熟期果皮厚度容易受環(huán)境因素影響,目前尚未有控制甜玉米果皮厚度發(fā)育相關(guān)基因被克隆的報(bào)道。本研究通過基因共表達(dá)網(wǎng)絡(luò)構(gòu)建與基因功能注釋信息分析,鑒定到與甜玉米乳熟期果皮厚度相關(guān)的重要候選基因。在Turquoise模塊中包含編碼細(xì)胞分裂素N-葡萄糖基轉(zhuǎn)移酶的基因ZmCNGT(Zm00001d019265),研究表明,該基因控制玉米籽粒中細(xì)胞分裂素的含量,并且該基因還可以觸發(fā)玉米生殖器官細(xì)胞的程序性死亡[34]。Fahima等[35]對(duì)荔枝施用外源細(xì)胞分裂素,發(fā)現(xiàn)荔枝果皮的細(xì)胞分裂受到顯著影響,荔枝果皮厚度變厚。在Turquoise模塊中鑒定到編碼赤霉素氧化酶的GA2ox9(Zm00001d008909),該基因影響玉米組織器官中的赤霉素(Gibberellin,GA)含量,通過對(duì)玉米葉片不同部位的激素含量測定,發(fā)現(xiàn)在細(xì)胞分裂區(qū)(Division zone)的赤霉素含量顯著高于其他部位,GA含量對(duì)細(xì)胞分裂的速度以及組織伸長具有重要的影響[36]。果膠是植物細(xì)胞壁的主要組成成分,它被合成后作為高度甲基酯化的聚合物分泌到細(xì)胞壁,Turquoise模塊中的ZmPME(Zm00001d050082)編碼果膠甲基酯酶,該酶可以將高度甲基酯化的聚合物去甲酯化,它的活性影響細(xì)胞與細(xì)胞之間的聚合力,進(jìn)而影響細(xì)胞壁松弛或緊密程度[37-39]。細(xì)胞壁的松弛程度會(huì)直接影響到細(xì)胞壁的厚度以及果皮細(xì)胞之間的排布緊密程度,而這2個(gè)因素是影響果皮厚度關(guān)鍵因素[40]。Yellow模塊中的PGL17(Zm00001d052103)控制聚半乳糖醛酸酶的合成,該酶是果膠降解的關(guān)鍵酶,參與各種發(fā)育過程,如果實(shí)成熟、細(xì)胞擴(kuò)增和器官脫落[41],可能是參與甜玉米果皮厚度調(diào)控的重要候選基因。

本研究鑒定的4個(gè)與果皮厚度顯著相關(guān)的模塊均與生長素轉(zhuǎn)錄因子構(gòu)成網(wǎng)絡(luò)關(guān)系,如ARF1、ARF7、ARF22、IAA2、IAA25等。其他作物研究發(fā)現(xiàn),生長素相關(guān)基因與果皮的厚度有關(guān),例如Liu等[42]利用番茄的gib-3突變體對(duì)番茄果皮細(xì)胞層數(shù)進(jìn)行研究發(fā)現(xiàn),施加IAA的突變體番茄果皮的平均細(xì)胞層數(shù)為24.5,而沒有施加IAA的番茄果皮的平均細(xì)胞層數(shù)只有13.3。除生長素轉(zhuǎn)錄因子外共表達(dá)網(wǎng)絡(luò)分析還鑒定到bZIP、ERF和MYB等轉(zhuǎn)錄因子與甜玉米果皮厚度具有共表達(dá)關(guān)系。bZIP轉(zhuǎn)錄因子可以調(diào)節(jié)多種生物過程,如種子成熟、生殖器官的發(fā)育、蔗糖信號(hào)轉(zhuǎn)導(dǎo)和病原體防御等[43]。通過前人的研究證實(shí)生長素、細(xì)胞分裂素、赤霉素、脫落酸等植物激素在植物的果皮發(fā)育過程中可能起著關(guān)鍵的作用,對(duì)果皮厚度和層數(shù)的變化都有著十分重要的影響。

本研究利用乳熟期不同果皮厚度甜玉米自交系M03和M08的籽粒果皮為材料,結(jié)合果皮厚度表型數(shù)據(jù)和轉(zhuǎn)錄組測序數(shù)據(jù),聯(lián)合差異表達(dá)基因分析和WGCNA分析共鑒定了4個(gè)與果皮厚度顯著相關(guān)的共表達(dá)模塊和13個(gè)核心基因。通過基因功能注釋發(fā)現(xiàn)這些模塊及其核心基因與果皮厚度調(diào)控緊密相關(guān)。同時(shí),在對(duì)模塊內(nèi)基因網(wǎng)絡(luò)關(guān)系的分析過程中發(fā)現(xiàn)了可能參與調(diào)控甜玉米乳熟期果皮厚度的其他基因和轉(zhuǎn)錄因子。結(jié)果表明,淀粉含量、細(xì)胞分裂和細(xì)胞凋亡等相關(guān)基因在甜玉米乳熟期果皮厚度的調(diào)控中發(fā)揮著重要的作用。本研究結(jié)果為甜玉米果皮厚度發(fā)育的分子機(jī)制研究提供了線索,也為培育優(yōu)質(zhì)口感的甜玉米品種提供理論支持。