整合GWAS和WGCNA筛选鉴定甘蓝型油菜生物产量候选基因

2021-06-09王艳花刘景森李加纳

作物学报 2021年8期

王艳花刘景森李加纳

王艳花1,2,**刘景森1,2,**李加纳1,2,*

1西南大学农学与生物科技学院/ 油菜工程研究中心, 重庆 400715;2西南大学现代农业科学研究院, 重庆 400715

生物产量是作物获得高产的重要基础, 对于甘蓝型油菜(L.)尤其重要。本研究利用588份甘蓝型油菜材料构成的自然群体2年生物产量表型数据的全基因组关联分析, 再结合高生物产量材料‘CQ45’和低生物产量材料‘CQ46’的转录组测序(RNA-seq)结果, 整合了6个甘蓝型油菜材料6个部位(茎秆、叶片、花后30 d主轴与侧枝种子、花后30 d主轴与侧枝角果皮)的转录组数据构建的加权共表达网络分析(WGCNA), 筛选出与生物产量相关的候选基因。通过相关分析发现, 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位点和极端表型差异基因确定候选基因, 它们的拟南芥同源基因为、、、, 这些基因在光合作用的卡尔文循环、碳同化、物质积累等方面发挥重要作用。

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

油菜作为四大油料作物(大豆、油菜、向日葵、花生)之一, 由于其较强的适应性, 在世界各地广泛种植[1-3]。随着人类对油菜的需求不断增加, 油菜总产量增速不断加快[4], 世界油料作物中油菜的种植面积及总产量已位居第二, 仅次于大豆[5-6]。同时油菜也是我国主要农作物, 在作物种植面积中排名第五, 仅次于水稻、玉米、小麦、大豆[7-12]。近年来, 油菜多功能开发利用的持续发展, 形成了油菜的油用、菜用、花用、蜜用、饲用等“一菜多用”的模式, 大幅度的提高了油菜的种植效益[13]。在实际生产中, 我国油菜产量仍较低, 生产成本却很高, 国际竞争力严重不足[14-15]。目前, 我国油菜超过60%依赖进口, 并且种植面积和单位面积产量都停滞不前, 这对我国粮食安全战略有着严重威胁[16]。因此, 提高油菜产量是我国油菜育种不断追求的目标[17-24]。

生物产量指单位面积土地上获得的不包括根系的作物干物质的总量, 它是描述茎、叶制造并积累有机物能力强弱的量[25]。生物产量体现作物总的生产力, 或者说作物光合作用碳同化的能力。多项报道已经表明, 生物产量与经济产量密切相关, 作物提高产量的方法通常有2种: 即在生物产量不变的情况下, 提高收获指数; 或收获指数不变的情况下, 提高生物产量[26-30]。近年来, 高通量测序技术的不断突破、成本不断降低, 多样本的转录组测序逐渐被应用于系统研究生命科学问题[31], 传统的少量样本比较分析已经无法有效处理海量的生物信息数据。因此, 生物信息分析工具和方法应运而生[32]。其中, 加权基因共表达网络分析(weighted gene co-expression network analysis, WGCNA)能够特异筛选出与目标性状高度关联的基因, 并进行模块(module)化分类, 得到具有高度生物学意义的共表达模块(co-expression module), 能够有效筛选到核心基因(hub genes)[33]。WGCNA算法作为一种精准、高效的生物信息学及生物数据挖掘方法, 被广泛应用于生物学各个领域, 利用WGCNA对马铃薯C16和C119品种构建抗逆生理性状高度关联的权重网络, 为深入探究马铃薯抗逆分子机制提供新的数据资源[34-35]。本研究以生物产量性状为主要研究对象, 通过GWAS分析, 鉴定出与生物产量性状紧密相关的SNP位点; 然后提取SNP位点前后共500 kb的序列片段作为候选区段, 提取候选区段内所有基因并结合转录组分析结果筛选获得差异候选基因; 同时, 结合转录组数据构建WGCNA共表达网络, 筛选鉴定出甘蓝型油菜生物产量关键候选基因。

1 材料与方法

1.1 供试材料

以来自世界各地的588份油菜自交系作为材料, 其中455份来自中国(大部分来自重庆、湖南、湖北等地), 11份来自亚洲其他地区, 102份来自欧洲, 13份来自北美洲, 7份来自澳大利亚。所有材料均由重庆市油菜工程技术研究中心提供。

1.2 材料种植及农艺性状考察

588份材料于2016、2018年种植于重庆北碚歇马油菜基地(29º45′39.99″, 106º22′38.47″, 海拔238.57 m)。每年播种时间为9月20日左右, 按完全随机区组设计, 2个重复同时种植, 每个材料种植2行, 每行15株, 行距40 cm, 株距20 cm。田间管理选择常规生产方式, 2个重复按照相同的管理措施和施肥水平进行, 收获时间为次年5月5日左右。于成熟期, 在每个小区中选择5株长势一致的植株。将植株分成上部分枝和下部茎秆两部分, 并将上部分枝装进网袋中。上部分枝和下部茎秆分别彻底晒干后, 使用天平秤进行称量, 将两部分重量相加得到植株的总干重。每份材料的5个样本干重的平均值为该材料的生物产量(biomass yield, BY)。

1.3 表型数据分析

利用Microsoft Excel软件初步整理588份自然群体的表型数据, 并计算单个环境及差值条下生物产量的平均值、标准差、变异系数等, 并用SPSS绘制生物产量的正态分布图; 利用SPSS统计分析软件对自然群体进行各性状进行相关分析; 利用R软件对多年环境下自然群体的生物产量表型数据进行最佳线性无偏预测(best linear unbiased prediction, BLUP)处理。

1.4 全基因组关联分析

利用本实验室已有588份甘蓝型油菜重测序的数据, 通过基因型分析, 最终获得385,691个可利用的SNP标记数据[36]。本研究利用这些标记结合重庆2017、2019年和BLUP环境的生物产量(BY)表型数据行全基因组关联分析。本研究使用6种统计模型在Tassel 5.2.1[37]软件中执行对甘蓝型油菜生物产量的关联分析, 即: 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图确定最佳模型, 并选择最佳模型下生物产量性状的全基因组关联分析结果作为进一步分析对象。本试验使用的SNP数据为385,691个, 因此规定值小于阈值(1/385,691=2.593E-06和0.05/385,691=1.296E-07)的位点为显著关联位点。利用Haploview[38]软件制作Manhattan图, 将通过关联分析检测到的染色体组上与各性状显著相关的标记位点可视化。

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层次聚类法对基因进行聚类, 按照混合动态剪切树的标准, 并设置每个基因网络模块最少的基因数为30。在使用动态剪切法在确定基因模块后, 我们依次计算每个模块的特征向量值(eigengenes), 然后对模块进行聚类分析, 将距离较近的模块合并成新的模块, 设置height = 0.25、deepSplit = 3、minModuleSize = 30。计算所得到模块的特征向量与样本来源组织的相关性, 选择相关系数>0.6的模块作为与组织部位显著相关模块。通过特征向量基因分析确定目标基因模块, 对所有的模块中具有代表性的基因-特征向量基因(module eigengene, ME)进行聚类分析, 进一步进行两两模块之间的ME相关性分析, ME之间的相关性越高, 其所在的模块的相关性也就越高; 将模块中的所有基因和ME基因在所有样本中分别进行表达水平分析, GO和KEGG分析。最后使用Cytoscape[39]将所关注模块的ME基因的共表达网络进行可视化。

2 结果与分析

2.1 表型数据统计及各性状的相关性分析

对2017年和2019年588份自然群体的生物产量的表型变异情况进行统计分析发现, 在不同环境下自群体的生物产量呈连续变异, 近似正态分布, 说明生物产量是受多基因控制的数量性状, 符合GWAS分析的要求。各产量相关性状的相关分析结果表明, 每角果粒数在2017年和2019年与茎秆干重、经济产量、收获指数和生物产量均呈正相关, 而茎秆干重仅在2017年呈现出显著正相关; 茎秆干重与上部生物产量、经济产量、生物产量在2年间均呈现显著正相关; 上部生物产量与经济产量、收获指数、生物产量在2年间均呈现显著正相关; 而经济产量和收获指数与生物产量在2年间同样呈现正相关。多数性状间的相关性均达到显著或极显著水平, 表明产量相关的各性状间存在显著相关作用(表1)。

2.2 甘蓝型油菜生物产量性状的GWAS分析

各生物产量性状在6种模型下的关联分析后得到的结果进行Quantile-Quantile散点图(QQ Plot)的绘制, 选择与预测值最接近的一种模型作为最佳模型, 2017年、2019年和BLUP环境下最佳模型为K+PCA模型(图1)。在最佳模型下, 利用本实验室已有的385,691个SNP, 以值小于阈值(1/385,692 = 2.593E-06)确定显著关联SNP位点, 并绘制曼哈顿图(图2)。生物产量2017、2019年和BLUP环境下共检测到9个显著关联标记位点(表2)。其中, 在2017年环境中检测到7个显著关联标记位点, 分别位于A03、A07、C01 (2个)、C04 (2个)、C06染色体, 并且位于同一条染色体的显著关联位点其置信区间高度重合; 2019年检测到2个显著关联位点, 分别位于A07和C09染色体; 而BLUP环境下共检测到4个显著关联位点, 分别位于A07、C03、C06 (2个), 位于C06染色体的2个显著关联位点其置信区间高度重合。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)。侧枝角果皮中分子功能相关性最高的是营养库活性(GO:0045735)等, 参与的生物过程中相关性较高的是对含氧化合物的反应(GO:1901700)等(图3-B); 这些差异基因相关性最高的是能量代谢途径(图3-C)。

花后30 d主轴种子中检测到6948个显著差异表达基因, CQ45对于CQ46有4080个显著上调表达基因, 2868个显著下调基因(图3-A)。分子功能相关性最高的是微管运动活性(GO:0003777)、硫葡萄糖苷酶活性(GO:0019137)等, 生物过程相关性较高的是细胞循环(GO:0007049)、细胞分裂(GO:0051301) (图3-B); 相关性较高的代谢途径是碳水化合物代谢中的淀粉和蔗糖代谢(图3-C)。

花后30 d侧枝种子中检测到17,309个显著差异表达基因, CQ45对于CQ46有11489个显著上调表达基因, 5820个显著下调基因(图3-A)。分子功能相关性最高的是微管运动活性(GO:0003777)、营养库活性(GO:0045735)硫葡萄糖苷酶活性(GO:0019137)等, 参与的生物过程相关性较高的是细胞循环(GO:0007049)等, 其结果与30 d主轴角果种子比较相似(图3-B); 相关性较高的是氨基酸代谢中的丙氨酸、天冬氨酸和谷氨酸代谢(图3-C)。

表2 最佳模型下生物产量显著位点表

(图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)等, 参与较多的生物过程是植物型次生细胞壁生物发生(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),共表达网络符合无尺度网络。使用动态剪切法在确定基因模块后, 依次计算每个模块的特征向量值, 然后对模块进行聚类分析, 将距离较近的模块合并成新的模块后共得到了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基因在所有样本中分别进行表达水平分析发现, 各个模块内的基因其在不同组织部位中的表达情况呈现不同的特性。每个模块内的基因表达水平高度相关, ME表达水平与模块整体表达水平也高度相关(图5-D和附图1), 说明目标模块的ME可充分代表其模块中的整体基因。

2.4.3 目标基因模块中的关键基因为获得与组织部位生长发育和功能高度相关的5个模块中的核心基因, 利用Cytoscape软件对基因互作调控网络进行可视化处理, 筛选出模块中共表达权重>0.1并且连通性排名前10的基因, 并结合NCBI数据库确定模块内的关键基因(图6)。Turquoise模块中确定了6个关键基因, 分别为、、、、、。Green模块中筛选出了4个核心基因, 分别为、、、。Magenta模块中筛选出了4个核心基因, 分别为、、、。Pink模块中筛选出了6个核心基因, 分别为、、、、、。Yellow模块中筛选出了8个核心基因, 分别为、、、、、、、。

2.5 候选基因的筛选

检测到显著性SNP位点上下游500 kb的LD置信区间内寻找甘蓝型油菜基因, 结合转录组差异表达基因, 通过WGCNA筛选到的与组织部位生长发育高度相关的关键模块中的核心基因(hub genes), 初步确定与甘蓝型油菜生物产量相关的基因, 将这些基因的蛋白质序列与拟南芥基因蛋白质序列进行BLAST对比, 结合前人已报道的拟南芥同源基因的功能以及特性, 筛选到生物产量性状中发挥重要作用的候选基因, 在拟南芥中的同源基因分别是、、、、、(表3)。

和在拟南芥中的同源基因是, 编码拟南芥叶绿体-1,6-二磷酸酶[40]。豌豆中叶绿体-1,6-二磷酸果糖()干扰载体的转基因植株表现出叶鲜重明显增加、光合作用的测量结果显示具有更高的碳同化率, 表明FBPase作为二氧化碳同化中的关键酶的作用, 并且还可以协调碳和氮的代谢[41]。在拟南芥中的同源基因是, 编码S-腺苷同型半胱氨酸水解酶(SAHH)。SAHH是维持真核生物甲基化稳态的关键酶[42], 会竞争性地抑制甲基转移酶(MT)活性, 使植株生长缓慢, 繁殖力低, 种子发芽减少[43-44]。在拟南芥的同源基因是, 编码在叶绿体基质中发现的一种小肽, 在黑暗中氧化的CP12与甘油醛-3-磷酸脱氢酶和磷酸布洛激酶(碳同化循环的两种酶)形成一种非活性超分子复合物[45],的转录物在决定成熟叶片和光合作用能力方面起着重要作用[46-48], Marri等人发现, CP12蛋白氧化对植物不同环境条件下光合过程的整体平衡至关重要[49]。在拟南芥中的同源基因是一种卡尔文循环酶, 具有光合固碳作用[48]。在拟南芥中的同源基因为, 磷酸化甘油醛-3-P脱氢酶(GAPC-1)是一种高度保守的胞浆酶, 纯合缺失突变株表现出生长延迟、角果形态改变和种子数量低[51]。在拟南芥中的同源基因是, 是拟南芥营养组织中主要肌动蛋白基因, 对植物生长发育有着极其重要的作用[52]。多项研究结果表明, 甘蓝型油菜生物产量相关候选基因的同源基因功能, 在植物氧化还原、能量和碳水化合物代谢、光合作用、物质积累等生长发育过程中发挥作用。

A: 样本聚类与数据矫正; B和C: 基于规模独立性和均值连通性选择软阈值; 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: 基因共表达网络模块与组织部位关联热图; 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)

表3 生物产量候选基因

3 讨论

生物产量体现作物总的生产力, 对作物的经济产量有着密不可分的关系。本研究所使用甘蓝型油菜自然群体中的生物产量变化范围广、变异系数大, 说明生物产量是比较复杂的性状。本研究在2017、2019年和BLUP环境下共检测到9个生物产量显著关联标记位点, 贡献率在5.64%~7.98%。与生物产量显著关联的SNP位点分布在各个染色体, Lu等[53]利用520份来自世界各地的甘蓝型油菜材料作为自然群体, 使用MLM进行全基因组关联分析, 在除了A01、A06、C02、C07以外的15条染色体上共检测到26个生物产量SNP, 贡献率在4.33%~17.03%之间; 60K SNP芯片对520分材料构成的自然群体进行全基因组关联分析, 在除了A01、A02和A10意外的16条染色体共检测到89个生物产量SNP位点, 贡献率在3.47%~8.36%之间[54]; Luo等[55]利用155分甘蓝型油菜材料构成的自然群体并使用60K SNP芯片检测到1个生物产量SNP, 其贡献率为0.76%。本研究中, 生物产量在2017、2019年和BLUP环境下共检测到9个显著关联标记位点, 分别位于A03、A07、C01、C03、C04、C06、C09染色体, 贡献率在5.64%~7.98%。前人检测到与生物产量显著关联的部分SNP位点也定位在这些染色体上, 说明生物产量的构成因素比较复杂, 且由多个SNP位点协同控制。

通过RNA-Seq对极端表型材料的差异表达基因分析表明, 叶片和角果皮主要富集到光合作用、淀粉和蔗糖代谢等生物途径, 这些代谢途径可能在油菜植株的物质积累过程中发挥重要作用, 因此, 这些差异基因的不同表达模式可能是造成2个材料生物产量有显著差异的主要原因。茎秆中的差异表达基因主要富集在细胞壁合成等途径, 油菜茎秆次生细胞壁的发育对其抗倒伏有着决定性影响, 茎秆的发育状况良好对油菜碳水化合物的运输积累有着积极作用[55]。因此, 本研究极端材料茎秆中的差异基因有可能在这些方面影响其生物产量。种子是油菜主要的储存和收获器官, 其发育状况直接影响着最终产量[56], 因此也直接影响生物产量。本研究极端生物产量材料种子中的差异基因主要富集在营养库活性、蔗糖淀粉代谢等功能和途径。

WGCNA是共表达网络分析有效的分析方法, 能够特异地筛选出与目标性状具有高度生物学意义的共表达模块, 在玉米等植物中已经被证明是一种高效的数据挖掘方法。本研究得到15个生物产量加权共表达网络模块, 分析发现一些模块与部分组织有着显著的相关性, 分别是叶片与Turquoise模块, 茎秆与Green、Magenta、Pink模块, 开花后30 d主轴种子与Yellow模块相关。叶片、茎秆和种子均在作物的生物产量的构成中发挥极其重要的作用, 直接影响到农产品最终的产量。结合GWAS和转录组差异基因得到的结果、整合WGCNA挖掘得到关联模块中的核心基因(hub genes), 最终筛选到一批可能与甘蓝型油菜生物产量密切相关的基因, 分别是、、、、、、。这些基因的干扰、敲除或者超量表达均引起了植株的整体重量的变化[40-52], 在拟南芥生物产量中发挥重要作用。

由于生物组学大数据存在复杂、多层次和信息互补的特点, 分析这些数据的一个关键目标是确定可预测表型性状的有效模型, 发现重要性状的关键调节因子并阐明其生物功能。本研究根据GWAS分析得到的与生物产量高度关联SNP标记, 同时结合转录组分析确定的差异基因以及WGCNA得到的基因富集模块进一步筛选出与甘蓝型油菜生物产量相关的基因。通过BLAST对比找到这些基因在拟南芥中的同源基因, 并根据其注释以及已有报道进一步筛选与生物产量高度相关的候选基因, 为油菜中生物产量候选基因的功能研究奠定基础。

4 结论

本研究通过2年间甘蓝型油菜自然群体中生物产量对大多数产量相关性状都具有正向效应, 说明甘蓝型油菜的生物产量是其他产量相关性状的基础和保障。结合GWAS关联基因与转录组差异基因得到178个基因, 根据拟南芥同源基因的功能, 筛选得到与生物产量性状相关的和为关键候选基因, 它们在能量和碳水化合物代谢以及光合作用等方面中发挥作用。此外, 选取WGCNA所得到的重点模块中连通性前10的基因作为hub gene, 通过同源基因功能注释分析, 得到生物产量候选基因、、、和, 在光合作用的卡尔文循环、碳同化、物质积累等方面发挥作用。

[1] 金森森, 罗帅. 全球视角下的南美大豆贸易发展趋势. 中国粮食经济, 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定位. 南京农业大学硕士学位论文, 江苏南京, 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] 张卫. 明确大宗油料作物发展规划提升国内大宗油料作物自给力. 中国食品, 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] 徐露萍. 浅谈油菜种植技术的应用与推广. 南方农业, 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] 李培武, 丁小霞, 张文, 谢立华, 周海燕. 中加双低油菜质量标准体系比较及我国发展对策探讨. 农产品质量与安全, 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] 张哲, 殷艳, 刘芳, 王积军, 傅廷栋. 我国油菜多功能开发利用现状及发展对策. 中国油料作物学报, 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] 姚东和, 杨民胜, 李志辉. 林分密度对巨尾桉生物产量及生产力的影响. 中南林业科技大学学报, 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] 李跃建, 宋荷仙, 朱华忠. 小麦收获指数、生物产量和籽粒产量的稳定性分析. 西南农业学报, 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的马铃薯根系抗旱相关共表达模块鉴定和核心基因发掘. 作物学报, 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 controlling 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] 柳寒. 油菜灌浆期种子和角果皮基因表达差异及相关基因的功能分析. 华中农业大学博士学位论文, 湖北武汉, 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)和高等学校学科创新引智基地项目(“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