甜玉米乳熟期果皮厚度相关基因的加权基因共表达网络分析

裴 虎,熊才运,张亚辉,任文闯,李小琴,黄 君

(广东省植物分子育种重点实验室,华南农业大学 农学院,广东 广州 510642)

甜玉米(ZeamaysL.sacharatasturt)是玉米属(ZeamaysL.)甜质型玉米亚种,是由su1、sh1、bt1、sh2等隐性突变基因引起的胚乳缺陷型玉米[1]。果皮是玉米籽粒的重要组成部分,具有保护胚和胚乳不受机械损伤或病虫害侵害的作用[2]。与饲用玉米不同,食味品质是决定甜玉米品种优劣和经济价值的重要因素,果皮厚度与甜玉米的柔嫩性和皮渣感有关,是影响其食味品质最为重要的关键指标之一[3]。目前,改良甜玉米的食味品质,尤其是果皮柔嫩度,选育果皮较薄的甜玉米品种,提升甜玉米口感,是育种家一直追求的重要育种目标。因此,深入挖掘控制甜玉米果皮厚度关键基因,能够为甜玉米品质改良提供基因资源,具有重要的研究意义。

截至目前,针对甜玉米果皮厚度的研究较少,主要集中在QTL(Quantitative trait locus)定位和遗传方式分析。玉米果皮厚度的遗传方式涉及显性效应、加性效应和上位性效应,其中加性效应和上位性效应起到主要作用[4-6]。前人研究表明,在玉米的10条染色体中除第7号染色体外,其余9条染色体均定位到玉米果皮厚度相关QTLs[5,7-8],但多数QTL难以在不同研究中被重复检测到。2020年,华南农业大学农学院甜玉米课题组利用甜玉米高代导入系群体对甜玉米果皮厚度进行QTL定位和转录组分析,在第10染色体上鉴定到一个控制果皮厚度的主效QTL-qPT10-5,能够解释表型变异的7.78%～35.38%,结合转录组测序和基因注释分析,认为AUX/IAA 转录因子,ZIM 转录因子和FATTY ACID EXPORT是控制甜玉米果皮厚度的关键候选基因[9]。

目前,通过正向遗传学方法已鉴定出多个与甜玉米果皮厚度相关的QTLs,但至今未有相关基因被克隆的报道。加权基因共表达网络分析(Weighted gene co-expression network analysis,WGCNA)是以转录组、基因芯片等获得的基因表达量数据为基础,将表达趋势相近的基因模块分类并得到具有高度生物学意义的模块和基因[10]。近年来,WGCNA方法已经广泛应用于玉米、水稻、高粱等农作物的遗传改良研究[11-14]。然而目前尚未有利用WGCNA对甜玉米果皮厚度相关的基因共表达模块被鉴定。

本研究以薄果皮甜玉米自交系M03和厚果皮甜玉米自交系M08为材料,在3个时间梯度(授粉后15,19,23 d)分析其果皮厚度变化,结合转录组和WGCNA分析构建果皮厚度特异性模块并挖掘出调控玉米果皮厚度的核心基因,为后续甜玉米果皮厚度分子遗传机制的进一步研究提供新的线索。

1 材料和方法

1.1 玉米材料及转录组测序数据

利用本课题组选育的2个甜玉米自交系M03(果皮厚度71.1～135.9 μm)和M08(果皮厚度111.7～171.7 μm)为研究材料。于2017年种植在华南农业大学增城教学实验基地,行距70 cm、株距25 cm、种植密度为52 500株/hm2。为避免材料出现其他品种玉米植株造成的花粉直感现象,采取人工套袋授粉的方式进行授粉,分别在授粉后15,19,23 d对2个自交系进行取样。每个自交系随机选取3个玉米果穗作为3个生物学重复,取10～15粒果穗中部的籽粒,用双刃刀片沿胚侧轻轻划开籽粒,用镊子剥取果皮并立即放入液氮冷冻,最后放入-80 ℃冰箱保存,用于后续表型鉴定及转录组测序。每个果穗取6粒籽粒进行扫描电镜观察测量果皮厚度。使用SPSS(Version 19)进行表型数据分析。

植物总 RNA 利用RNA提取试剂盒(Vazyme,中国南京)进行提取。提取的 RNA使用Nanodrop 2000分光光度计(Thermofisher,美国)、Qubit 2.0荧光计(Life Technologies)和Aglient 2100生物分析仪系统(Agilent Technologies,CA,美国)来评估纯度(OD260/280≥1.8;OD260/230≥1.0),浓缩RNA样品(总RNA≥250 ng/μL)和完整性(RIN≥8.0,28S/18S≥0.5)。转录组测序文库构建主要使用寡核苷酸(dT)磁珠(Invitrogen,Carlsbadcity,CA,美国)从总RNA中分离出每个样品1.5 μg mRNA,片段化并反向转录成cDNA,并通过 PCR进行扩增。使用Illumina HiSeq 2000平台(San Diego,CA,美国)进行了转录组测序。使用 FastQC软件对原始测序数据进行质量评估,去除低质量和未知碱基数目超过15%的 reads。利用 Hisat 2 软件将Clean reads 比对到玉米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软件计算样品中基因的表达量,每千个碱基外显子百万片段数(Fragments per kilobase of exon per million fragments mapped,FPKM)值衡量基因的表达水平。DESeq2 软件(Version 1.10.1)用于差异表达基因的计算[15]。只有当不同组间|Log2Fold Change|>1 且P值<0.05的基因被认定为差异表达基因。

1.2 加权基因共表达网络分析

通过TBtools v1.09[16]软件中的WGCNA shinny plugin对共表达模块进行划分。去除掉FPKM值小于1的基因,本研究共得到14 126个基因用于样品聚类。通过SFT and Power Select 选项进行软阈值(Soft Thresholding Power)的筛选,计算得到阈值为14(图1)。利用Module-net中的Cluster选项将min Module Size调整为30,Module Cuttree Height调整为0.25,选择Max Block Size为15 000对数据进行模块聚类和剪切。利用Module-trait选项计算基因在模块内的连接度(Eigengene connectivity,KME),取KME值排名靠前且差异表达的基因作为该模块内核心基因(Hub gene)进行后续分析,利用Cytoscape v3.6.0[17]对核心基因网络进行可视化。

B纵坐标代表每一个软阈值对应的网络平均连接程度;C、D.在power值为14时拓扑网络符合程度的检测。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 基因共表达网络软阈值的确定Fig.1 Soft threshold determination of gene co-expression network

1.3 差异表达基因及特异性模块功能富集分析

将筛选出的差异表达基因和特异性模块内所包含的基因利用R程序中的 clusterProfiler 软件包[18]进行GO(Gene Ontology)和 KEGG(Kyoto Encyclopedia of Genes and Genomes)富集分析。当P值(P-value)小于0.05时认为此 GO功能和 KEGG pathway 功能存在显著富集情况。使用 Tbtools 软件和Omicshare平台制作热图对基因表达量进行可视化。

1.4 转录组测序数据的qRT-PCR验证

利用实时荧光定量PCR(Quantitative real-time PCR,qRT-PCR)验证RNA-seq结果的可靠性。使用HiScriptⅢ RT SuperMix for qPCR(+gDNA wiper)试剂盒将果皮RNA反转录成cDNA作为qRT-PCR模板。以ubiquitin(登录号:NM_001138130)为内参基因,使用CFX96(Bio-Rad)实时荧光定量 PCR 系统,选用2×T5 Fast qPCR Mix(SYBR GreenI)(TsingKe)试剂盒,荧光染料为SYBR Green。所有反应均在20 μL的体系下进行,其中包含200 ng cDNA、0.8 μmol/L核心基因特异性引物以及10 μL的2×T5 Fast qPCR Mix。qPCR反应程序:95 ℃预变性1 min;95 ℃ 10 s,55 ℃ 5 s,40个循环;溶解曲线95 ℃ 10 s,55 ℃ 10 s,72 ℃ 15 s。使用2-ΔΔCT法分析基因相对表达量[19],所用基因引物见表1。

2 结果与分析

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

在相同的栽培条件下,分别取自交系M03和M08授粉后15,19,23 d的果皮进行扫描电镜观察(Scanning electron microscope,SEM)(图2-A—F)。研究结果显示,随着授粉时间的增加,甜玉米籽粒果皮厚度整体呈现逐渐降低的趋势,授粉后15～19 d的果皮厚度降低的速度较为缓慢,从授粉后19 d开始果皮厚度呈现出快速变薄的趋势。其中,M08的果皮厚度相比M03更厚,不同自交系间果皮厚度的差异随时间的推移逐渐增加,在乳熟期(授粉后19 d)开始达到显著水平(图2-G)。

A—C.M03授粉后15,19,23 d的果皮扫描电镜结果;D—F.M08授粉后15,19,23 d的扫描电镜结果;G.M03和M08的果皮厚度变化,柱中字母为 Duncan′s多重比较结果,不同字母表示材料间在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个不同时期的籽粒果皮扫描电镜图片以及果皮厚度变化Fig.2 SEM of seed pericarp and changes in pericarp thickness of M08 and M03 at three different periods

2.2 差异表达基因筛选及GO和KEEG富集分析

2.2.1 转录组数据分析 去除含有接头污染和低质量的reads,共得到480 628 079条高质量reads,与玉米B73参考基因组进行序列比对,其比对率为78.82%～ 84.50%,外显子比对率为87.38%～91.97%(表2)。

表2 参试样品转录组数据质量Tab.2 Quality of the transcriptome data of each sample

通过M03和M08自交系内部不同发育时期以及不同发育时期材料间的比较,共筛选出4 748个DEGs。其中,M03材料内部比较(M03-DAP15 vs M03-DAP23)共鉴定到1 234个DEGs(上调641个,下调593个);M08材料在不同发育时期的果皮厚度差异相对于M03较少,其内部比较(M08-DAP15 vs M08-DAP23)获得了较少的差异基因,为471个DEGs(上调255个,下调216个)。相同发育时期的材料间比较(M03-DAP15 vs M08-DAP15;M03-DAP19 vs M08-DAP19;M03-DAP23 vs M08-DAP23)获得了更多的差异基因,分别为2 334,2 238和2 039个DEGs(图3-A)。材料内部与材料间比较共有的差异基因为11个(图3-B)。

图3 不同材料不同时期DEGs数量统计及韦恩图Fig.3 Number of DEGs in different datasets and venn diagram

2.2.2 差异基因富集分析 对同一材料不同发育时期及不同材料相同发育时期进行差异表达基因鉴定,对获得的差异表达基因进行GO功能注释,GO富集条目主要分为生物过程、分子功能及细胞组分3 类。结果显示,差异表达基因GO功能富集主要富集到细胞组分组织(GO:0071840)、对刺激的反应(GO:0050896)、生物过程调节(GO:0050789)等生物过程;结构分子活性(GO:0005198)、核酸结合转录物(GO:0001071)、转运活性(GO:0005215)等分子功能;膜(GO:0016020)、细胞组分(GO:0044464)、细胞(GO:0005623)等细胞组分(图4)。

KEGG代谢通路富集分析结果显示差异表达基因主要参与丙氨酸、天冬氨酸和谷氨酸代谢(Alanine,aspartate and glutamate metabolism)、萜类骨架的生物合成(Terpenoid backbone biosynthesis)和植物MAPK信号通路(MAPK signaling pathway-plant);M03与M08材料内部的不同发育时期比较(M03-DAP15 vs M03-DAP23与M08-DAP15 vs M08-DAP23)结果显示,差异表达基因还可共同富集到氨基酸生物合成(Biosynthesis of amino acids)、淀粉和蔗糖代谢(Starch and sucrose metabolism);相同时期不同材料之间比较得出差异表达基因主要参与β-丙氨酸代谢(Beta-alanine metabolism)、精氨酸和脯氨酸代谢(Arginine and proline metabolism)、半胱氨酸和蛋氨酸代谢(Cysteine and methionine metabolism)途径(图4)。

2.3 共表达基因聚类和模块构建

通过对FPKM值小于1的基因进行过滤,最终获得14 126个基因用于构建基因聚类树和共表达模块分析。根据基因的FPKM值将表达模式相似的基因进行聚类,利用动态切割法(Dynamic tree cut)将得到的基因聚类树切割成不同模块(图5-A)。将切割得到的18个不同模块分别用不同的颜色进行区分,其中,基因数目最多的模块是Turquoise模块,最少的是Grey60模块(图5-B)。

A.基因聚类树,每个聚类树对应一个基因模块,聚类等级相同的模块用同样的颜色表示;B.聚类出的18个模块以及每个模块所包含基因个数。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 基因共表达模块以及各模块所包含的基因数Fig.5 Gene co-expression module and the number of genes contained in each module

2.4 特异性模块筛选及目标模块KEGG富集分析

由于同一模块内的基因表达模式相近,因此,这些基因很可能具有相似的生物学功能。结合果皮厚度表型数据计算每个模块与果皮厚度之间的相关性,最终在18个模块中得到了4个相关性系数绝对值大于0.60的模块,分别是Pink模块(-0.63)、Yellow模块(-0.82)、Turquoise模块(0.66)和Magenta模块(0.81)(图6)。通过模块之间的相关性热图可以看出4个模块之间具有较高的相关性,可将其作为目标模块进行后续分析。为进一步探究与甜玉米果皮厚度显著相关模块的生物学功能,本研究利用OmicShare Tools(https://www.omicshare.com/)对模块进行KEGG功能富集分析。结果表明,除Magenta模块外,其他3个模块均富集到氨基酸的生物合成、α-亚麻酸代谢(Alpha-linolenic acid metabolism)、半胱氨酸和蛋氨酸代谢等通路。此外,泛酸和辅酶A生物合成(Terpenoid backbone biosynthesis)在Turquoise模块被特异性富集;脂肪酸延长(Fatty acid elongation)在Yellow模块特异富集;泛酸和辅酶A的生物合成(Pantothenate and CoA biosynthesis)在Pink模块特异富集;植物MAPK信号通路、植物病原体相互作用等代谢通路只在Magenta模块特异富集(图7)。将不同材料及不同发育阶段鉴定到的DEGs与所有所有模块相关基因的KEGG富集结果进行对比,结果显示,不同方法得到的基因可以共同富集到丙氨酸、天冬氨酸和谷氨酸代谢、植物MAPK信号通路、半胱氨酸和蛋氨酸代谢等通路。

A.各模块与果皮厚度之间的相关系数热图,红色框为特异性模块;B.特异性模块与DEGs共有基因数。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 表型-模块相关性热图和模块-DEGs共有基因数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 核心基因挖掘及互作网络构建与分析

将不同模块内的基因与不同材料、不同发育时期鉴定到的差异表达基因进行比对,选择二者的交集进行基因互作网络构建。通过筛选,Turquoise模块中有1 351个基因、Yellow模块中有543个基因、Pink模块中有10个基因、Magenta模块中有34个基因用于互作网络构建。将每个模块中与表型相关性较低(|GS|<0.5)的基因进行过滤,选取模块内连通性(KME)最高的前20个基因作为核心基因的候选基因。为获取候选基因的功能注释信息,通过玉米参考基因组数据库MaizeGDB(https://www.maizegdb.org/)和拟南芥同源注释对相关基因的功能进行查找,最终筛选出与甜玉米果皮厚度相关的核心基因(Hub gene)。在Turquoise模块共筛选出3个核心基因,Yellow模块筛共选出4个核心基因,Pink模块共筛选出2个核心基因,Magenta模块共筛选出4个核心基因(表3)。将核心基因作为互作网络的核心节点构建互作网络,筛选基因网络中的转录因子(Transcript factor,TF)和与果皮厚度相关的其他基因,并选择权重值(Weight)最高的前50个基因绘图,利用Cytoscape v3.6.0软件进行可视化(图8)。

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

A～D分别为Turquoise模块、Yellow模块、Pink模块、Magenta模块的基因网络;红色大点为核心基因,灰色矩形为模块内与表型相关基因的标注。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 特异性模块的基因共表达网络及其核心基因Fig.8 Gene co-expression network and hub genes of significant modules

通过基因功能注释发现核心基因与已报道的可能与果皮厚度发育相关的基因具有互作网络关系。例如:Turquoise模块中的核心基因与细胞分裂素N-葡萄糖基转移酶ZmCNGT(Zm00001d019265),GA2氧化酶GA2ox9(Zm00001d008909)和果胶甲酯酶ZmPME(Zm00001d050082)等具有互作网络关系;Yellow模块中的核心基因与聚半乳糖醛酸酶PGL23(Zm00001d013032)具有互作网络关系。这些核心基因互作网络中还包含ARF、AUX/IAA、bZIP、MYB、ERF、WRKY等可能与细胞分裂与果皮厚度发育相关的转录因子。

2.6 核心基因的实时荧光定量PCR验证

为验证转录组测序结果中基因表达量的可靠性,本研究挑选了9个核心基因进行qRT-PCR表达量验证。结果表明,这些基因在2个不同的甜玉米自交系授粉后15,19,23 d的果皮组织中的表达量水平与转录组测序结果一致(图9),进一步说明了转录组测序结果的可靠性以及数据的可重复性。

3 结论与讨论

甜玉米作为鲜食玉米的一种,果皮厚度是决定其口感优劣的重要因素[20]。果皮角质化程度、细胞大小、细胞数量细胞壁厚度以及细胞的排布等因素都可能会影响果皮的厚度[21]。这些因素除受基因型影响外还与栽培环境、气候变化、采收时间等有关。因此,本研究在同等的栽培条件下,严格控制了授粉时间和采收时间,最大程度避免了上述因素对甜玉米果皮厚度表型的影响。为解析甜玉米果皮厚度的遗传基因,鉴定控制甜玉米果皮厚度的关键候选基因,本研究选择果皮厚度有明显差异的2个甜玉米自交系(M03和M08)进行转录组测序,通过基因差异表达分析和加权基因共表达网络分析,选取二者的交集进行后续研究。WGCNA作为识别基因与样本间相关性模式的常用方法,可以特异性的筛选出与目标性状相关的基因并聚类成共表达模块,是相关性状协同表达调控控制解析的重要工具[22]。

3.1 差异表达基因和特异性模块内基因与碳水化合物代谢、淀粉和蔗糖代谢及MAPK信号通路有关

本研究利用乳熟期果皮较薄的甜玉米自交系M03和果皮较厚自交系M08的乳熟期果皮材料的RNA-seq数据为基础进行了差异表达分析,在相同材料的不同发育时期及不同材料的相同发育时期共鉴定了4 748个差异表达基因,对相关基因进行GO和KEGG功能富集分析发现“碳水化合物代谢过程”、 “淀粉和蔗糖代谢”等与甜玉米果皮厚度发育有关。在禾本科植物研究中,如小麦、高粱等作物,发现其籽粒果皮中含有较多的淀粉颗粒,果皮中淀粉的含量会随着籽粒的发育而变化[23-24]。在小麦果皮厚度发育研究中发现小麦植株受精后珠被开始发育为果皮,淀粉随着果皮的发育不断积累,在DAP(Day after pollination)5左右果皮中淀粉含量到达峰值,随后逐渐分解并在DAP22左右消失殆尽,在果皮细胞淀粉含量处于峰值时可以明显观察到细胞呈现出轻微的膨胀变形,说明果皮细胞中的淀粉含量可以通过影响细胞大小进而影响果皮厚度,不同的淀粉含量可能会造成果皮厚度的差异,与本研究结果一致,甜玉米果皮厚度发育的关键基因均与碳水化合物代谢、淀粉和蔗糖代谢相关[24]。

通过WGCNA分析,本研究共鉴定到了4个具有高度生物学意义并与果皮厚度显著相关的模块,通过富集分析揭示了这些模块的生物学功能。将所有模块的KEGG富集结果与差异表达基因的功能富集结果进行比较,鉴定到天冬氨酸、丙氨酸与谷氨酸代谢、植物MAPK信号通路、半胱氨酸和蛋氨酸代谢等通路与甜玉米果皮厚度的发育有关。其中MAPK信号通路涉及包含细胞防御、抗逆、细胞增殖和细胞凋亡在内的多种生命活动,推测该通路在调节甜玉米果皮厚度发育过程中具有十分重要的作用。

3.2 MYB转录因子、脱落酸、细胞壁pH值及细胞周期蛋白能够影响甜玉米果皮厚度

在4个特异性模块中选择与差异表达基因的交集进行后续分析,共筛选到13个核心基因,这些核心基因可能在参与调控果皮厚度发育过程中具有重要作用。Turquoise模块的Zm00001d023282和Yellow模块的Zm00001d022227分别为MYB69和MYB70转录因子,编码MYB结构域蛋白。前人在玉米MYB转录因子的研究中发现,ZmFDL1/MYB94功能的缺失导致ω-羟基脂肪酸和多羟基脂肪酸的合成显著降低,而角质层又是由这些超长链脂肪酸(VLCFA)化合物组成,该基因的功能缺失会造成玉米发育前期包括籽粒、幼苗叶片等部位的角质层形成受阻,同时该研究还指出干旱胁迫介导ABA的产生会抑制ZmFDL1的表达,从而增加角质层的厚度防止水分过快蒸发[25]。在拟南芥研究中发现几个与表皮蜡质生物合成的基因受MYB转录因子家族成员调控,包括AtMYB30、AtMYB94和AtMYB96[26-28]。本研究鉴定到的核心基因Zm00001d005884为bZIP转录因子基因,在拟南芥的同源基因为AT1G49720,该基因与脱落酸反应元件结合,介导ABA反应的转录调控。核心基因Zm00001d014291在拟南芥的同源基因为AT4G33500,该基因编码蛋白磷酸酶2C家族蛋白,在番茄中的研究表明,ABA与受体结合会影响编码蛋白磷酸酶蛋白的基因SIPP2C5的表达,该基因表达量的变化会影响果实的各种性状,包括果形、果皮厚度、可溶性固形物含量等[29]。核心基因Zm00001d012460编码BAG家族蛋白,该蛋白与细胞凋亡调节有关,Earp等[23]通过对禾本科作物的研究发现,除了高粱以外其他大部分禾本科作物随着籽粒的成熟,果皮细胞将经历PCD过程,使果皮呈现出厚度逐渐降低的趋势。Magenta模块中的核心基因Zm00001d041214编码种质膜H+-ATP酶,该基因的突变会导致拟南芥植株发育受限,生长素诱导的信号级联反应中,最终作用到质膜H+-ATP 酶,被生长素信号激活增强的质子泵将H+送进细胞壁降低pH值,激活细胞壁内对pH值敏感的酶和蛋白质,造成细胞壁松动、延伸和生长[30]。核心基因Zm00001d048143编码一种细胞周期蛋白,细胞周期蛋白依赖性激酶(CDK)是一个重要的丝氨酸/苏氨酸蛋白激酶家族,在细胞有丝分裂中起着重要的作用[31-32]。前人在番茄果皮发育的研究中发现,过表达或抑制CDK蛋白相关基因的表达会显著影响番茄果皮的细胞层数和厚度[33]。

3.3 特异性模块内与果皮厚度相关基因挖掘

甜玉米乳熟期果皮厚度容易受环境因素影响,目前尚未有控制甜玉米果皮厚度发育相关基因被克隆的报道。本研究通过基因共表达网络构建与基因功能注释信息分析,鉴定到与甜玉米乳熟期果皮厚度相关的重要候选基因。在Turquoise模块中包含编码细胞分裂素N-葡萄糖基转移酶的基因ZmCNGT(Zm00001d019265),研究表明,该基因控制玉米籽粒中细胞分裂素的含量,并且该基因还可以触发玉米生殖器官细胞的程序性死亡[34]。Fahima等[35]对荔枝施用外源细胞分裂素,发现荔枝果皮的细胞分裂受到显著影响,荔枝果皮厚度变厚。在Turquoise模块中鉴定到编码赤霉素氧化酶的GA2ox9(Zm00001d008909),该基因影响玉米组织器官中的赤霉素(Gibberellin,GA)含量,通过对玉米叶片不同部位的激素含量测定,发现在细胞分裂区(Division zone)的赤霉素含量显著高于其他部位,GA含量对细胞分裂的速度以及组织伸长具有重要的影响[36]。果胶是植物细胞壁的主要组成成分,它被合成后作为高度甲基酯化的聚合物分泌到细胞壁,Turquoise模块中的ZmPME(Zm00001d050082)编码果胶甲基酯酶,该酶可以将高度甲基酯化的聚合物去甲酯化,它的活性影响细胞与细胞之间的聚合力,进而影响细胞壁松弛或紧密程度[37-39]。细胞壁的松弛程度会直接影响到细胞壁的厚度以及果皮细胞之间的排布紧密程度,而这2个因素是影响果皮厚度关键因素[40]。Yellow模块中的PGL17(Zm00001d052103)控制聚半乳糖醛酸酶的合成,该酶是果胶降解的关键酶,参与各种发育过程,如果实成熟、细胞扩增和器官脱落[41],可能是参与甜玉米果皮厚度调控的重要候选基因。

本研究鉴定的4个与果皮厚度显著相关的模块均与生长素转录因子构成网络关系,如ARF1、ARF7、ARF22、IAA2、IAA25等。其他作物研究发现,生长素相关基因与果皮的厚度有关,例如Liu等[42]利用番茄的gib-3突变体对番茄果皮细胞层数进行研究发现,施加IAA的突变体番茄果皮的平均细胞层数为24.5,而没有施加IAA的番茄果皮的平均细胞层数只有13.3。除生长素转录因子外共表达网络分析还鉴定到bZIP、ERF和MYB等转录因子与甜玉米果皮厚度具有共表达关系。bZIP转录因子可以调节多种生物过程,如种子成熟、生殖器官的发育、蔗糖信号转导和病原体防御等[43]。通过前人的研究证实生长素、细胞分裂素、赤霉素、脱落酸等植物激素在植物的果皮发育过程中可能起着关键的作用,对果皮厚度和层数的变化都有着十分重要的影响。

本研究利用乳熟期不同果皮厚度甜玉米自交系M03和M08的籽粒果皮为材料,结合果皮厚度表型数据和转录组测序数据,联合差异表达基因分析和WGCNA分析共鉴定了4个与果皮厚度显著相关的共表达模块和13个核心基因。通过基因功能注释发现这些模块及其核心基因与果皮厚度调控紧密相关。同时,在对模块内基因网络关系的分析过程中发现了可能参与调控甜玉米乳熟期果皮厚度的其他基因和转录因子。结果表明,淀粉含量、细胞分裂和细胞凋亡等相关基因在甜玉米乳熟期果皮厚度的调控中发挥着重要的作用。本研究结果为甜玉米果皮厚度发育的分子机制研究提供了线索,也为培育优质口感的甜玉米品种提供理论支持。