基于RNA-Seq的向日葵脂肪酸合成基因筛选与分析

郭树春,张艳芳,李素萍,邵 盈,邢慧军,于海峰*

(1 内蒙古自治区农牧业科学院 作物科学研究所,呼和浩特 010031;2 内蒙古农业大学 园艺与植物保护学院,呼和浩特 010010;3 四子王旗林业和草原局,内蒙古乌兰察布 011800)

植物油脂是由脂肪酸和甘油化合而成的天然高分子化合物。它是人体必需营养成分,可为人体提供必要的能量、维持特定的功能,同时它还是脂溶性维生素的载体[1]。脂肪酸中的单不饱和油酸、多不饱和亚油酸都是有益脂肪酸,其含量及组成是决定其营养价值和用途的关键[2]。向日葵(HelianthusannuusL.)是世界四大油料作物之一。向日葵籽油(又称葵花籽油)是日常消费、食品加工中的一种优质植物油,其脂肪酸(fatty acid,FA)主要由油酸、亚油酸、棕榈酸和硬脂酸组成[3]。向日葵籽油饱和脂肪酸(棕榈酸和硬脂酸)相对含量很低,介于4%～13%之间,不饱和脂肪酸(油酸、亚油酸)含量高,在87%～94%之间[4-5],不饱和脂肪酸中油酸和亚油酸含量呈负相关关系。早在1999年联合国粮农组织根据向日葵油酸含量的不同,将其基因型分为三大类:①低油酸向日葵(传统向日葵),油酸含量14%～39%;②中油酸向日葵,油酸含量42%～72%;③高油酸向日葵,油酸含量75%～91%。现在市场上的向日葵油原料主要是传统向日葵,其不饱和脂肪酸中亚油酸占55%～65%,油酸占20%～30%[6]。油酸因其分子结构中不饱和烯键数目少于亚油酸而更加稳定,高油酸向日葵油烹制的食品货架期会更长,更加健康[7-8],因此市场对高油酸向日葵油更加青睐。现已证实,向日葵籽油中油酸和亚油酸的占比是决定其品质的关键,其脂肪酸含量和组分由基因型决定,还受到一些环境因子如温度、光照和湿度等影响[9-10],且这些环境因子中对其影响最大的是温度[11-12]。进一步研究表明,随着温度的升高,低油酸向日葵籽中的油酸含量明显升高,亚油酸含量随之降低,而温度越低,则相反。而高油酸向日葵中的油酸含量对环境变化的敏感性较低,其脂肪酸组成变化也较小[10]。但是,从转录组层面在不同生长条件下对向日葵脂肪酸代谢的分子机制研究较少。

随着中国农业的高质量发展,不同脂肪酸组分表型的向日葵需求将不断增加。为了迎合社会发展的需要,弥补相关研究的不足,本研究以高油酸自交系J9、低油酸自交系P50为试验材料,采用错期播种方法,通过转录组分析,解析向日葵脂肪酸合成及组分变化的分子调控机制,促进高油酸向日葵分子育种稳步发展,为向日葵产业高质量发展提供理论依据。

1 材料和方法

1.1 试验材料

以自育高油酸向日葵自交系J9(油酸相对含量大于85.0%)和低油酸自交系P50(油酸含量小于35.0%)为试验材料,于2022年4月20日开始分期播种,每间隔10 d设1个播期,共设6个不同的播种期,每个播期处理设3次重复。试验在呼和浩特市(N40.49°,E111.41°,海拔106 m)内蒙古自治区农牧业科学院试验田进行,试材采取相同的栽培管理措施。

1.2 样品采集

在开花授粉后(days after flowering,DAF)20 d,取整齐一致单株的向日葵籽仁为试验样品,累计采集12份样品,液氮速冻。编号为:J9_S1、J9_S2、J9_S3、J9_S4、J9_S5、J9_S6和P50_S1、P50_S2、P50_S3、P50_S4、P50_S5、P50_S6,每份样品采集6份,3份(3次重复)用于转录组测序,3份(3次重复)用于脂肪酸含量与组分测定。采样期气象数据来源于中国气象数据网(http://data.cma.cn/)。

1.3 脂肪酸组分与含量测定

利用液相色谱法测定各样品(累计36个)主要脂肪酸相对含量,如:油酸、亚油酸、硬脂酸、亚麻酸、棕榈酸等组分参照GB 5009.168—2016第三法测定其含量。

1.4 转录组测序

利用TaKaRa公司的RNAisoPlus提取样本总RNA[13],并文库构建,构建好的文库用Illumina HiSeqTM 4000进行转录组测序。

1.5 测序数据过滤与组装

使用Trinity软件对转录组的De novo进行从头组装[14]。原始序列数据经杂质去除后得到clean reads,然后将每个样品clean reads用短reads比对软件TMAP比对到参考基因序列上(DDBJ/ENA/GenBank,MNCJ00000000[15]),组装得到的Unigene[16],得到最终的Unigene信息,具体信息分析过程及分析流程见图1。

图1 数据处理流程

1.6 基因表达量计算和功能注释

使用RPKM法(reads per kb per million reads)[17]计算基因的表达量,两样本间的差异表达基因筛选参考Audic[18]基于测序的差异基因检测方法,此方法发表于GenomeResearch,公式为:MRPKM=(1000000×C)/(N×L/1000)。通过控制FDR(fa1se discovery rate)≤0.001且︱log2R︱≥1(倍数差异在2倍以上)的值校正P阈值。

RPKM法公式各符号含义:A基因的表达量用RPKM(A)表示,则C表示能唯一比对到基因A区域的Reads数,N为唯一比对到参考基因的总reads数,L为基因A的碱基数。具体计算过程:首先将测序得到的reads匹配到参考基因上,将能匹配到参考基因上的reads数标准化到106,以此为基础来换算某一基因上匹配上的reads数,再除以该基因长度就得到该基因的表达量(RPKM值)。RPKM不仅对测序深度作了归一化,而且对基因长度也作了归一化,使得不同长度的基因在不同测序深度下得到的基因表达水平估计值有了可比性,是目前最可靠的基因表达估计方法,计算得到的基因表达量可直接用于比较不同样品间的基因表达差异分析。

每个样品中的每个Unigene表达量确定后,为了便于后续对基因的功能研究,对每个基因进行功能注释,将每个基因序列首先跟NCBI的NR库用(http://www.geneontology.org/)BLAST软件进行比对,再将每个基因注释到NT、Swiss-Prot、KEGG、COG、GO功能库。例如:用Blast2 GO软件注释到GO的各级条目下。基因序列通过BLAST比对到KEGG数据库进行KEGG注释,Q≤0.05的Pathway定义为在差异表达基因中显著富集的Pathway,通过Pathway显著性富集能确定差异表达基因参与的主要生物代谢途径和信号转导通路[19]。

2 结果与分析

2.1 油酸含量测定结果

错期播种同期采样试验条件下,测定自交系J9和P50在6个播期的油酸相对含量,并收集采样期(从开始授粉到20 DAF)平均温度。结果(图2)显示:高油酸品系J9油酸相对含量受日平均温度变化影响不显著,其油酸相对含量基本稳定在85.0%左右;而低油酸品系P50油酸相对含量随日平均温度的降低而显著降低,在S1播种期油酸相对含量最高达42.0%,而S6播种期油酸相对含量仅为20.49%。由此可见,向日葵脂肪酸组分和含量易受温度影响,且高油酸向日葵品系脂肪酸组分和含量受温度影响小,而低油酸品系易受温度影响。因此,笔者设定了错期播种同期采样试验,在授粉期不同温度下筛选向日葵脂肪酸合成基因,并研究其调控关系。

*表示两品系在同一时期的显著水平(P<0.05)。下同。

2.2 RNA-Seq测序质量评估结果

通对36个(12个样品,3个重复)样品使用Hiseq 4000测序平台进行RNA-Seq测序,2个品系共测得265.79 G数据,平均每个样品数据约7.38 G,单个样品数据最少的为J9_S2_2,共测得6.12 G数据;最多的为P50_S2_2,共测得9.06 G数据。样品平均Q20为96.59%,平均Q30为92.93%。碱基错误率分布符合一般规律。GC含量平均值为46.50%,且不存在明显分离现象,结果见表1。上述结果显示,RNA-Seq测序数据质量达到预期,可用于后续分析。

表1 测序数据质量评估

2.3 RNA-Seq测序重复性检验

根据每个样本基因表达量的数值,计算皮尔逊相关系数的平方(R2),得出重复样本之间的相关性。结果(图3)为:在J9品系中,相关系数最高的为J9_S5_1和J9_S5_3之间,相关系数高达0.935;相关系数最低的为J9_S1_2和J9_S1_3之间,相关系数为0.754。在P50品系中,相关系数最高的为P50_S4_1和P50_S4_3之间,相关系数高达0.918;相关系数最低的为P50_S3_2和P50_S3_1之间,相关系数为0.769。平均相关系数高于0.8。因此,本研究中各个生物学重复间相关性较好,样品的重复性良好。

图中横向、纵向数字1,2,3表示3个重复。

2.4 Unigene基因表达量分析和功能注释

利用测序得到的clean reads,采用RPKM算法确定基因的表达量。首先将某一样品测序得到的reads匹配到参考基因上,并将reads数标准化为106,以此为基础来换算该基因上匹配上的reads数,再除以该基因长度得到该样品该基因的基因表达量,每个样本每个基因的基因表达量用RPKM方法表示。按RPKM值的大小,将基因分为高表达基因(RPKM≥10)、中表达基因(1

图4 不同向日葵品系基因表达量的叠加直方图

对上述所有Unigene进行功能注释,这些基因被注释到了GO功能库下的1 283个GO条目下。在3个GO本体(ontology)中,生物学过程(biological process,BP)中条目最多,包含665个条目;其次分子功能(molecular function,MF)和细胞组分(cellular component,CC),分别包含451,166个条目。在图5中列出了位于3个GO本体中前10位的GO条目,KEGG库注释的Unigene共计3 286个,其中位于前20位中的角质层、木栓质和蜡生物合成(cutin, suberine and wax biosynthesis)属于脂质代谢。

图5 所有Unigene GO和Unigene KEGG分类结果

2.5 差异表达基因筛选与分析结果

首先对品种内不同处理样品差异表达分析,将同一品种不同播期品与其对照(S1)间的比对所获得的差异表达基因总和定义为S-CK数据集,S-CK数据集(S2-CK、S3-CK、S4-CK、S5-CK、S6-CK,累计5个数据集,表2)反映了J9和P50在脂肪酸积累关键时期,在不同播期下条件下,品种内响应环境变化的DEGs;把J9和P50在相同播期的测序结果经比对所获得的差异表达基因之和定义为S数据集(累计6个数据集,表2),该数据集反映J9和P50在脂肪酸积累各关键时期,品种间响应环境变化的DEGs。

表2 差异表达基因筛选数据集命名

对S-CK数据集进行品种内差异分析发现,在S-CK数据集中,高油酸品系J9在S5-CK中的差异基因数(different expression genes,DEGs)最多,为1 022个;低油酸品系P50在S3-vs-S1中的差异基因数最多,为1 712个。上述结果表明,高油酸品系在S5响应环境变换的差异表达基因较多,而低油酸品系在S3期响应环境变换的差异表达基因较多;总体上高油酸品系响应环境变换的差异表达基因较少,这和2.1节向日葵籽粒中油酸含量和温度变化关系研究结果相吻合。

在基因表达调控类型研究中,S-CK数据集的5种比较(S2 vs S1、S3 vs S1、S4 vs S1、S5 vs S1和S6 vs S1)中,高油酸品系J9中含有258,382,100,311,207个上调DEGs,有551,467,346,711,1 427个下调DEGs;低油酸品系P50中含有12,930,34,305,873个上调DEGs,有12,782,42,355,980个下调DEGs(图6,A和图6,B);上述与油酸合成相关DEGs的获得,为油酸合成的共表达分析以及通路分析提供了基础。

A.品种内上调差异分析结果;B.品种内下调差异分析结果;C.品种内特异表达分析。粉色填充表上调;绿色填充表下调;黄色线表J9品系;蓝色线表P50品系。

上述结果经试验内特异表达分析显示,在S-CK数据集的5种比较(S2 vs S1、S3 vs S1、S4 vs S1、S5 vs S1和S6 vs S1)策略中,高油酸品系J9中含有372,279,115,322,147个特异上调DEGs;有664,386,718,660,1 249个特异下调DEGs(图6,C),分析发现J9和P50的品种特异性表达基因数远超其品种共表达数量,该结论证明了本研究材料选择的可靠性,提取特异表达基因功能条目,为后续研究提供信息。同时,为了进一步确定影响高油酸和低油酸品系的关键差异基因,对S数据集(来源于RNA-Seq数据)进行品种间特异表达分析(见表2),结果如图7。由图7可知,在S数据集中,无论是上调基因还是下调基因,各时期品种间特异表达基因相当。最后,品种特异性表达基因和品种间特异表达基因去除重复后,在J9和P50在-CK和S数据集得到15 885个与向日葵脂肪酸代谢相关的品种间特异表达DEGs用于后续功能分析。

图7 J9 vs P50在相同播种期差异表达基因分析

2.6 特异表达基因功能分析

将2.5节分析得到的15 885个特异表达DEGs映射到GO注释中,经GO功能富集,得到与脂质代谢相关的GO条目:如,脂质生物合成过程(GO:0008610)、脂肪酸代谢过程(GO:0006631)、类固醇生物合成过程(GO:0006694)和类固醇代谢过程(GO:0008202)等。J9在中特异表达DEGs的GO注释中,累计获得96与脂质代谢相关DEGs(表3),P50中获得43个(表4);J9和P50在品种间特异表达GO条目映射研究中,获得278个DEGs(表5)。

表3 J9特异表达DEGs GO term映射结果

表4 P50特异表达DEGs GO term映射结果

表5 J9和P50特异表达DEGs GO term映射结果

总之,在品种内特异表达分析和品种间特异表达分析中,累计403个DEGs被注释到与脂质代谢相关的GO条目下。在品种内特异表达分析结果中,将通过GO功能富集获得的与脂质代谢直接相关的403个DEGs映射到与脂质代谢相关的KEGG代谢通路,累计29个DEGs直接参与了19条相关KEGG通路(图8,A)。其中,显著富集于糖基磷脂酰肌醇(GPI)生物合成KEGG通路的DEGs 5个、脂肪酸延伸的4个、脂肪酸代谢的4个、类固醇生物合成的3个、脂肪酸生物合成的3个、萜类骨架生物合成的3个、不饱和脂肪酸生物合成的2个、α-亚麻酸代谢的2个和甘油磷脂代谢通路的3个(图8)。

图8 特异表达DEGs的KEGG功能富集

3 讨 论

将在品种间特异表达分析结果与品种内特异表达分析得到DEGs和去除重复DEGs,共得到42个DEGs,经分析其隶属于22种酶。以此为基础绘制了向日葵脂肪酸代谢相关DEGs通路图,向日葵脂肪酸代谢过程论述如下。

首先,向日葵通过光合作用合成蔗糖(suc),然后被蔗糖转用酶(suc transporter,SUC)运送到向日葵种子细胞中,经过磷酸化、糖酵解、氧化生成脂肪酸合成的碳源乙酰辅酶A(acetyl-CoA),在质体(plastid)中开始其脂肪酸的从头合成(图9)。首先,乙酰-CoA在乙酰-CoA羧化酶(ACCase)催化下,与酰基载体蛋白(acyl carrier protein,ACP)结合,生成向日葵脂肪酸合成的直接原料丙二酰-CoA(malonyl-CoA);然后,丙二酰-CoA在脂肪酸合成酶复合体(fatty acid synthase,FAS)或酰基载体蛋白转移酶(MCAT)催化下[20]生成脂肪酸碳链延长供体丙二酰-ACP(malonyl-ACP);同时,乙酰-CoA,ACP转酰基酶将1个丙二酰基团由乙酰-CoA转移至ACP的丝氨酸残基上,生成乙酰-ACP,并释放CoA(图10)。本试验中,基因110889119和110925594参与了脂肪酸合成的限速酶乙酰-CoA羧化酶(ACCase)的调节,基因110889119在J9和P50的脂肪酸积累前期相对表达量偏高,其在J9中更高;基因110925594P50却在P50中的相对表达量更高,说明该基因更容易受到环境因素的影响。

图9 脂肪酸生物合成通路中的表达模式

其次,由β-酮酯酰-ACP还原酶(β-ketoacyl-ACP reductase,KAR)催化乙酰-ACP生成β-ketoacyl-ACP,然后在NADPH等酶的参与下,生成β-羟基丁酰-ACP(β-hydroxyacyl-ACP),此步是脂肪酸生物合成的最初还原步骤。本试验中,有累计6个基因参与了上述过程的调节,其中调节β-酮酯酰-ACP还原酶(KAR)的2个基因在A-试验中表达特征明显,基因110864403在积累前期低表达,后期高表达,而基因110893771则表现相反。4个基因参与了NADPH等酶的调节,基因1109382561和110868833表达特征明显,基因110868833在整个播期在J9中表达量较高,而1109382561则相反,它们在A-试验中也表现出同样的特征;值得注意的是,虽然基因110868833在J9整个脂肪酸积累期表达量相对较高,但是其在油酸积累关键期(20 DAFs)表达量最高。

再次,β-羟基丁酰-ACP由β-羟基丁酰-ACP脱水酶(β-hydroxacyl-ACP dehydratase,HAD)催化[21],在其α位和β位碳原子之间脱水生成2,3-反式丁烯酰(巴豆酰)-ACP(enoyl-ACP)。然后,反式丁烯酰-ACP被多功能蛋白(multifunctional protein,MFP)催化[22],转化生成相应的丁酰-ACP(butyryl-ACP)。基因110893751参与了此过程中MFP的激活,该基因在J9中脂肪酸合成28 DAFs和35 DAFs表达量最高。

最后,饱和丁酰-ACP再经转酰基作用与新的丙二酸单酰-ACP依次重复上述4个过程,经数次循环后,脂肪酸合成通过终止反应来终止合成,可生成C16∶0和C18∶1饱和脂肪酸;然后,在脂酰基-ACP硫脂酶(acyl-ACP thioesterase,Fat)作用下,将ACP上的酰基部分(C16∶0和C18∶1饱和脂肪酸)通过水解反应释放出来,成为游离脂肪酸,并释放1个巯基-ACP,使反应终止[23]。碳链聚合终止后,产生的不同碳链长度游离脂肪酸,由酰基-CoA合成酶(acyl-CoA sythetase,ACS)负责催化合成脂酰-CoA(C16∶0-CoA和C18∶1-CoA),并采用某种未知机制从质体转运到内质网或胞质中,完成脂肪酸链的延长。此过程中,基因110890752和110871748参与了C16∶0和C18∶1饱和脂肪酸的合成,而基因110941185和110893751参与了脂酰-CoA的合成与转运。

4 结 论

在向日葵籽粒形成过程中,虽然其脂肪酸组分和含量随着积累时间的推移而变化,但是向日葵脂肪酸主要以不饱和脂肪酸油酸和亚油酸为主,其合成的关键转折时期是20 DAF;通过播种期调控可以明显影响向日葵脂肪酸组分和含量,向日葵脂肪酸组分和含量易受温度影响,且高油酸向日葵材料脂肪酸组分和含量受温度影响小,而低油酸材料易受温度影响。在错期播种试验条件下,对36份样品进行转录组测序,2个品系共测得265.79 G数据。

试验中,J9在S5-vs-S1中的差异基因数最多,为1 022个;P50在S3-vs-S1中的差异基因数最多,为1 712个。经基因差异表达分析,得到的15 885个特异表达DEGs,经GO功能富集,得到与脂质代谢相关的GO条目累计获得96与脂质代谢相关DEGs;最后,J9和P50品种间特异表达GO条目映射研究中,获得278个DEGs。

为进一步揭示向日葵种子发育过程中与油脂积累及油酸生物合成相关的关键基因,将通过GO功能富集获得的与脂质代谢直接相关的403个DEGs映射到与脂质代谢相关的KEGG代谢通路,累计29个DEGs直接参与了19条相关KEGG通路;去除重复DEGs,共得到43个DEGs,隶属于22种酶。其中KAS、ACX、DGAT、MGD和OLE是向日葵不饱和脂肪酸生物合成和油脂积累的关键基因。