基于转录组测序的花生籽粒不同发育时期油脂合成相关基因差异表达分析

陈玉梅，李璐璐，陈锦玲，徐 媛，李惠敏，秦新民

(广西师范大学 生命科学学院，广西 桂林 541004)

花生是我国重要的油料作物和经济作物[1]，其含油量高达40%～60%，蛋白质含量20%～30%[2]，具有很高的经济价值和营养价值。目前，花生在我国农作物中的种植面积居第7位，但年产值居第4位[3]。我国不仅是花生生产和消费大国，同时也是出口大国，出口量占世界花生市场贸易总量的40%。据报道，我国花生消费中46%～48%用于榨油[4]。相关研究表明，榨油原料含油量每提高1个百分点，油脂加工的纯利润可提高7%[5]。长期以来，我国花生育种目标主要是高产、早熟和抗逆性等[6]。但对于一个优良的花生品种来说，高含油量是一个重要指标，提高花生含油量可以大大增加花生的经济价值。谷建中等[7]通过对23个高油酸花生品种(系)进行遗传多样性分析，发现开选016为适应性广、配合力高的优异高油酸花生种质,具有较高的育种价值。通过分子生物学手段在转录水平上深入研究花生油脂合成相关基因，对提高花生油脂含量、改善花生品质等具有重要意义。

参与植物油脂合成的3个关键酶分别为甘油-3-磷酸酰基转移酶(GPAT)、溶血磷脂酸酰基转移酶(LPAAT)和二酰甘油酰基转移酶(DGAT)[8]。其中，DGAT是催化三酰甘油(TAG)合成的最后一步反应的关键酶和限速酶，在油脂合成过程中起着重要作用[9]。目前，在微藻[10]、油茶[11-12]、油棕[13]、油桐[14]、核桃[15]和白檀[16]等多种油料植物研究中已经证实，参与油脂合成的相关基因和酶对于种子油脂合成具有重要的调节作用。近年来，关于花生油脂合成方面的研究已经取得一定的进展。华方静[17]发现，同一花生品种的不同发育时期，Ah-GPAT9基因的表达量差异极显著，种子含油量与该基因表达量呈正相关，并且表达峰值高或峰值持续时间较长均有利于油脂合成。陈四龙等[18]通过构建花生cDNA文库和EST测序发现，AhLPAT基因的表达量与油脂的积累速率变化一致，进一步通过转基因试验验证，成功转化该基因的拟南芥植株,其种子含油量显著增加，说明该基因对种子的油脂合成起到正向调控作用。荧光定量PCR检测结果表明，GPAT、LPAAT、DGAT等3类酰基转移酶基因的表达总量与花生和甘蓝型油菜的籽仁含油量呈显著或极显著正相关[19-20]。并且，DGAT基因能够使酿酒酵母TAG合成缺陷株恢复TAG的合成与积累[21]。和小燕等[22]通过转录组测序对高油和低油2个花生品种进行分析发现，有四大类代谢的484条基因与油脂代谢相关。目前，关于花生油脂合成相关基因的研究大多是对关键基因的单独研究，对花生油脂合成的整个调控网络研究尚不深入，鉴于此，对4个不同发育时期的花生籽粒进行转录组测序，筛选与油脂合成相关的差异表达基因以及转录因子等，为进一步了解花生油脂合成相关基因的调控模式及调控网络奠定基础，同时也为高油花生品种的培育提供依据。

1 材料和方法

1.1 试验材料

供试材料为花生(ArachishypogaeaL.)品种桂花37，由广西壮族自治区农业科学院经济作物研究所提供。

1.2 试验方法

1.2.1 样品采集 在花生盛花期，每天对当天盛开的花朵进行挂牌标记。选取长势一致、无病虫害的花生植株，分别取花后20、35、50、60 d (分别为籽粒发育的前期、中期、成熟期、成熟后期)4个不同发育时期的花生荚果，剥去花生壳后迅速将籽粒(分别以HS20、HS35、HS50、HS60表示)放入液氮中速冻20 min，再转至-80 ℃超低温冰箱保存，作为转录组测序的材料。

1.2.2 花生籽粒总RNA的提取、建库以及测序 提取4个不同发育时期的花生籽粒总RNA，对检测合格的RNA样品进行建库和高通量测序，建库和测序方法参考文献[23]。

1.2.3 差异表达基因的筛选 对4个不同发育时期的花生籽粒进行转录组测序，将得到的数据按发育时期的先后顺序进行两两比对，即分为HS20-vs-HS35、HS35-vs-HS50和HS50-vs-HS60 3个比较组，旨在分析出在4个时期中都有的差异表达基因，再进一步对差异表达基因进行功能注释。对表达量的原始数据进行标准化处理，并进行泊松分布计算，对差异检验的P值作多重假设检验校正，减少假阳性。差异表达基因设定为FDR≤0.001且倍数差异在2倍以上，同时校正P值<0.05。

1.2.4 差异表达基因的GO功能注释分类以及富集分析 根据差异表达基因筛选结果进行GO(Gene ontology)功能分类，主要包括分子功能、细胞组分和生物过程三大功能类。使用R软件的Phyper函数进行富集分析，计算P值，并进行FDR校正，FDR≤0.01视为显著富集。

1.2.5 差异表达基因的KEGG Pathway分类和富集分析 对差异基因进行生物通路分类，使用R软件中Phyer函数进行富集分析，得到显著富集的通路，对P值进行FDR校正，FDR≤0.01视为显著富集。

1.2.6 转录因子家族分析 转录因子是一类能与5′端上游特定序列专一性结合，从而保证目的基因以特定强度在特定时间与空间表达的蛋白质分子。通过getorf检测Unigene的ORF，然后通过hmmsearch检索，比对出转录因子的蛋白质结构域，根据植物转录因子数据库(PlantTDB)描述的转录因子家族特征对Unigene进行功能鉴定。

1.2.7 差异表达基因的RT-PCR验证 通过Tri-zol法提取花生籽粒4个时期的总RNA，采用Aidlab公司反转录试剂盒(TUREscript 1st Strand cDNA Synthesis Kit)进行反转录操作，利用Analytikjena-qTOWER 2.2型荧光定量PCR仪进行荧光定量PCR反应，各个样品中目的基因相对表达量通过qPCRsoft 3.2软件的Pfaffl method公式计算。

2 结果与分析

2.1 不同发育时期花生籽粒总RNA的提取与质量检测

对不同发育时期的花生籽粒样品HS20、HS35、HS50、HS60提取总RNA并进行质量检测。分光光度计检测结果表明，各样品的总RNA样品OD260/280在2.04～2.90，RIN/RQN值在9.1～9.4，其完整性高，符合后续试验的要求。

2.2 不同发育时期花生籽粒转录组文库测序结果质量分析

对花生籽粒样品HS20、HS35、HS50、HS60进行转录组测序，得到的原始数据(Raw reads)经数据过滤，除去低质量、接头污染以及未知碱基N含量大于5%的Reads，得到过滤后数据(Clean reads)。各个时期所得到的Clean reads占原始数据的94%以上，其中Q20高达97%以上，Q30高达90%以上，说明转录组测序数据质量较高(表1)，可进行后续生物信息学分析。

表1 测序数据统计Tab.1 Sequencing date statistics

2.3 不同发育时期花生籽粒差异表达基因的筛选

根据不同发育时期花生籽粒样品基因表达水平，检测显著差异表达基因，结果如图1所示。在HS20-vs-HS35中，差异表达基因总数为25 769个，其中上调的差异表达基因有3 235个，下调的差异表达基因有22 534个；HS35-vs-HS50中，差异表达基因总数为18 403个，其中上调的差异表达基因有11 579个，下调的差异表达基因有6 824个；HS50-vs-HS60中，差异表达基因总数为12 308个，其中上调的差异表达基因有3 235个，下调的差异表达基因有9 073个。

图1 不同发育时期花生籽粒差异表达基因火山图Fig.1 The volcano plot of differential expression genes of the different developmental stages in peanut seed

2.4 不同发育时期花生籽粒差异表达基因的GO功能注释分析

2.4.1 GO功能分类 HS20-vs-HS35、HS35-vs-HS50、HS50-vs-HS60差异表达基因的GO功能分类结果如图2—4所示，在HS20-vs-HS35的差异表达基因中，参与分子功能(Molecular function)、细胞组分(Cellular component)、生物过程(Biological process)的基因分别有15 649、22 489、18 913个；在HS35-vs-HS50中，分别有11 445、16 336、13 383个；在HS50-vs-HS60中，分别有7 782、10 741、9 216个。

图2 HS20-vs-H35差异表达基因GO分类图Fig.2 GO classification map of differential expression genes of HS20-vs-H35

2.4.2 GO功能富集分析 对不同发育时期的花生籽粒的差异表达基因进行GO功能富集分析，HS20-vs-HS35的差异表达基因被富集到5 790个GO 条目中，其中，参与生物过程、细胞组分、分子功能的差异表达基因含义不同分别占56.98%、12.31%、30.71%；HS35-vs-HS50的差异表达基因被富集到5 301个GO 条目中，参与生物过程、细胞组分、分子功能的差异表达基因含义不同分别占56.97%、11.94%、31.09%；HS50-vs-HS60的差异表达基因被富集到4 826个GO 条目中，参与生物过程、细胞组分、分子功能的差异表达基因分别占57.44%、12.20%、30.36%。通过分析各组分占比可以发现，3个比较组中大多数差异表达基因参与生物过程，其次是分子功能，最后是细胞组分。3个比较组中差异表达基因最显著的20个GO富集功能如图5—7所示。

图3 HS35-vs-HS50差异表达基因GO分类图Fig.3 GO classification map of differential expression genes of HS35-vs-HS50

2.5 不同发育时期花生籽粒差异表达基因的KEGG Pathway富集分析

对不同发育时期花生籽粒差异表达基因进行KEGG Pathway富集分析，将所有的差异表达基因分为135类代谢通路，HS20-vs-HS35、HS35-vs-HS50、HS50-vs-HS60的代谢通路以及差异基因数量如表2所示。结果表明，在代谢通路中，参与花生油脂代谢的通路主要有脂肪酸生物合成(ko00061)、不饱和脂肪酸的生物合成(ko01040)、脂肪酸延长(ko00062)、脂肪酸代谢(ko01212)、花生四烯酸代谢(ko00590)等。其中，在HS20-vs-HS35、HS35-vs-HS50和HS50-vs-HS60 3个比较组中，参与脂肪酸生物合成的差异表达基因分别有63、63、43个；参与不饱和脂肪酸生物合成的分别有45、38、24个；参与脂肪酸延长的分别有78、54、31个；参与脂肪酸代谢的分别有113、98、67个；参与花生四烯酸代谢的分别有57、56、34个。

2.6 基因转录因子家族分类

对具有编码转录因子能力的基因进行预测，并对基因所属的转录因子家族进行分类，结果如表3所示。由表3可见，可能参与油脂合成的转录因子有ABI3VP1家族、AP2-EREBP家族以及C2C2-Dof家族的成员，其中,属于ABI3VP1家族的转录因子有154个、属于AP2-EREBP家族的转录因子有236个， 属于C2C2-Dof家族的转录因子有66个。

图4 HS50-vs-HS60差异表达基因GO分类图Fig.4 GO classification map of differential expression genes of HS50-vs-HS60

图5 HS20-vs-HS35差异表达基因GO富集图 Fig.5 GO enrichment map of differential expression genes of HS20-vs-HS35

图6 HS35-vs-HS50差异表达基因GO富集图 Fig.6 GO enrichment map of differential expression genes of HS35-vs-HS50

图7 HS50-vs-HS60差异表达基因GO富集图 Fig.7 GO enrichment map of differential expression genes of HS50-vs-HS60

表2 不同发育时期花生籽粒代谢通路及其差异表达基因数量Tab.2 The metabolic pathways and the number of differential expression genes of the different developmental stages in peanut seed

续表2 不同发育时期花生籽粒代谢通路及其差异表达基因数量Tab.2(Continued) The metabolic pathways and the number of differential expression genes of the different developmental stages in peanut seed

续表2 不同发育时期花生籽粒代谢通路及其差异表达基因数量Tab.2(Continued) The metabolic pathways and the number of differential expression genes of the different developmental stages in peanut seed

续表2 不同发育时期花生籽粒代谢通路及其差异表达基因数量Tab.2(Continued) The metabolic pathways and the number of differential expression genes of the different developmental stages in peanut seed

表3 基因所属转录因子家族的分类Tab.3 Classification of transcription factor families of genes

2.7 不同发育时期花生籽粒差异表达基因的RT-PCR验证

为验证转录组测序结果中基因表达量的准确性，以18SRNA为内参基因，选取转录组测序中的4个差异表达基因进行RT-PCR验证，4个差异表达基因和内参基因的引物如表4所示。由图8可见，RT-PCR验证的基因表达量变化趋势与转录组测序结果(表5)中基因表达量的变化趋势基本一致，表明转录组测序所得到的数据具有较高的可靠性。

表4 RT-PCR验证的差异表达基因及其引物Tab.4 The differential expression genes and primers of RT-PCR validation

图8 差异表达基因的RT-PCR验证结果Fig.8 The RT-PCR validation results of the differential expression genes

表5 4个转录组测序基因的表达量Tab.5 The expression of four genes sequenced by transcriptome

3 结论与讨论

植物种子中的油脂主要以三酰甘油的形式储存，其积累呈慢—快—慢的变化规律[24-25]。脂肪酸是油脂的重要组成部分，本研究中发现，花生不同发育时期差异表达基因涉及脂肪酸合成的代谢途径主要有脂肪酸生物合成、不饱和脂肪酸的生物合成、脂肪酸延长、脂肪酸代谢、花生四烯酸代谢等几种代谢通路。HS20-vs-HS35、HS35-vs-HS50、HS50-vs-HS60参与脂肪酸生物合成的差异表达基因分别有63、63、43个，参与不饱和脂肪酸生物合成的分别有45、38、24个，参与脂肪酸延长的分别有78、54、31个，参与脂肪酸代谢的分别有113、98、67个，参与花生四烯酸代谢的分别有57、56、34个。这些基因的表达呈现出慢—快—慢的趋势。如参与脂肪酸合成代谢的基因fabG在花后20、35、50、60 d时的表达量(Reads per kb per million reads，FPKM)分别为3.67、551.79、1 960.28、2 100.08。参与亚油酸代谢途径的基因LOX1_5的表达量分别为3.05、53.18、1 096.15、1 268.16，脂肪酸代谢中的FAD2基因表达量分别为23.68、153.80、226.88、268.06。另外，参与油脂合成的相关基因和酶中，催化三酰甘油合成的最后一步反应的DGAT基因有79个。其他参与油脂合成的酶，如酰基ACP 硫酯酶(Acyl-ACP thioesterase)基因有11个，磷酸烯醇式丙酮酸羧化酶(Phosphoenolpyruvate carboxylase)基因有23个，β-酮脂酰-ACP 合酶(β-ketoacyl-ACP synthase)基因有3个，β-酮脂酰-ACP 还原酶(β-ketoacyl-ACP reductase)基因有2个。

相关研究表明，多种转录因子在种子油脂合成过程中起调控作用[26]。ABI3、FUS3和LEC2转录因子在油脂合成和种子成熟过程中相互调节，其中，ABI3主要参与蛋白质含量的调节，FUS3主要调节脂质的合成。另外，LEC2和FUS3在种子胚的不同发育阶段发挥作用[27]。WRI1(WRINKLED1)转录因子是AP2/EREBP家族的一员。拟南芥中过表达WRI1可以上调脂肪酸合成的一组基因，包括丙酮酸激酶、乙酰辅酶A羧化酶、酰基载体蛋白和酮脂酰ACP合成酶等基因[28]。玉米中过表达ZmLEC1转录因子可以提高含油量[29]。Dof(DNA binding with one finger)是植物特有的转录因子，N末端含有一个高度保守的C2-C2单锌指结构，能够特异识别启动子序列中的AAAG/CTTT元件，从而调节植物基因的表达[30]。相关研究表明，Dof转录因子调控脂肪酸的合成[31]。大豆中已有28个Dof基因被鉴定参与油脂合成调控，过表达GmDof4和GmDof11基因可以提高乙酰辅酶A羧化酶和长链乙酰辅酶A合成酶基因的表达，同时GmDof4和GmDof11基因转化拟南芥植株的种子中，总脂肪酸和脂质含量都有所增加[32]。本研究中，与油脂合成相关的转录因子有ABI3VP1、AP2-EREBP、C2C2-Dof家族的成员，分别有154、236、66个。

本研究中，HS20-vs-HS35、HS35-vs-HS50和HS50-vs-HS60的差异表达基因总数分别是25 769个、18 403个和12 308个，其中上调表达基因分别占12.55%、62.92%和26.28%。 从上述统计数据可以发现，在花生籽粒油脂合成前期，下调差异表达基因数量(22 534)所占比例大，籽粒油脂积累缓慢；随着籽粒的不断生长，上调差异表达基因数量(11 579)所占比例上升，多于下调差异表达基因，有利于促进籽粒中油脂的迅速积累；而到了油脂合成后期，下调差异表达基因数量(9 073)所占比例较大，又多于上调差异表达基因，导致了籽粒中油脂合成速度变慢。以上数据表明，花生油脂的积累是个非常复杂的过程，受到许多基因和转录因子的共同调节，构成了一个非常复杂的调节网络，其主要功能基因及其在油脂合成与积累过程的相互关系尚有待深入研究。