不同叶龄红小豆叶片的转录组分析

董伟欣, 张 磊, 李东晓, 张月辰

(1.河北农业大学农学院/华北作物改良与调控国家重点实验室,河北 保定 071000;2.河北开放大学,河北 石家庄 050080)

红小豆(Vignaangularis)起源于中国[1],是重要的杂粮作物之一,红小豆不仅具有丰富的营养物质,还具有较高的药用和保健价值[2],因此在人类食谱中占有重要地位.然而,红小豆作为一种小杂粮作物,被重视程度不够,品种混杂退化,加之小豆花荚脱落严重,产量低且不稳定成为生产中的难题[3].利用红小豆苗期日照时间调控开花时间对于后期产量和品质的提升具有显著影响,因此,通过转录组学对不同叶龄红小豆叶片进行研究,这对于深入了解红小豆生长发育的代谢过程具有重要意义.

随着测序技术不断发展和成本逐渐降低,转录组学已经广泛应用到不同物种.大豆转录组测序分析表明,lncRNA77580是响应逆境和生长发育过程中有重要功能的长链非编码RNA,另外,通过转录组测序技术挖掘控制大豆百粒重的基因以及糖代谢过程中抗旱基因的表达情况,为大豆抗逆基因的深入研究以及产量提升奠定了基础[4-6].在其他作物上,对野生型粳稻品种‘中花11’进行转录组学研究,鉴定出响应低温胁迫基因31个,其中25个基因上调表达,6个基因下调表达,这为揭示水稻响应低温的分子机制提供了依据[7].玉米转录组测序分析表明,在低磷胁迫下玉米叶片的表型变化与相关基因存在调控规律[8].

转录组测序技术已经被应用到植物育种、逆境胁迫、遗传多样性分析、种质资源挖掘等方面.关于红小豆转录组方面的研究,仅在锈菌侵染[9]和矮秆窄叶[10]等方面有报道,在其生长发育过程中相关基因表达却未见报道.基于此,本试验选择晚熟品种‘唐山红小豆’ 3叶龄和6叶龄叶片作为试验材料,对其进行转录组测序分析,探究红小豆不同叶龄时期生长发育过程中的差异基因数目及表达情况,为深入挖掘红小豆苗期生长过程中显著差异基因的动态变化提供依据.

1 材料与方法

1.1 试验设计

试验于2021年7月在河北农业大学教学实验基地(东经:115°47′,北纬:38°87′)进行,小区长5 m,宽1 m,种植2行,株行距为15 cm×40 cm,随机区组排列,共10个小区,选取3叶龄(真叶展平后5 d,简称:5 d)和6叶龄(真叶展平后15 d,简称:15 d)无病害、长势良好的‘唐山红小豆’(种子由河北省农林科学院粮油作物研究所提供)三出复叶的中间叶片,8～10株混合在一起,每个样品3次重复,共6个样品.以自然光照时间(natural photoperiod,NPD)作为标记,分别为NPD-5 d-1、NPD-5 d-2、NPD-5 d-3、NPD-15 d-1、NPD-15 d-2、NPD-15 d-3,采样后立刻放入液氮带回实验室,于-80 ℃冰箱中保存.

1.2 RNA的提取与检测

以红小豆三出复叶的中间叶片为材料提取RNA,RNA采用试剂盒(mirVanaTMmiRNA ISOlation Kit, Ambion-1561)提取,参照说明书步骤进行测定,用Agilent 2100 Bioanalyzer检测RNA样品的纯度、浓度和完整性.

1.3 测序与分析

试验采用双端测序,以红小豆(网址:https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/001/190/045/)的基因组作为参考基因组[11].用带有Oligo(dT)的磁珠富集mRNA,加入打断试剂将mRNA打断成短片段,以打断后的mRNA为模板,用6碱基随机引物合成一链cDNA,配制二链合成反应体系合成二链cDNA,并使用试剂盒纯化双链cDNA,纯化的双链cDNA再进行末端修复,加A尾并连接测序接头,然后进行片段选择,最后进行PCR扩增,构建的文库利用Illumina HiSeqTM 2500测序仪进行测序.

1.4 差异基因筛选

FPKM[12]和count使用bowtie2[13]和eXpress[14]软件分析得到.通过eXpress软件获取各个样本中转录本(protein coding)的reads数目,使用DESeq R package的estimateSizeFactors函数对数据进行标准化,并使用nbinomTest函数计算差异比较的P-value和fold change[15].

1.5 差异基因富集分析

挑选出P<0.05、差异倍数>2的差异转录本,使用Blast2 GO软件对得到的所有差异基因进行GO注释,用WEGO软件对所有的差异基因进行GO功能分类统计.GO注释包括3个主要功能组别,即:生物学过程(biological process)、分子功能(molecular function)和细胞组分(cellular component).通过BLAST软件将基因序列比对到京都基因和基因组百科全书(Kyoto Encyclopedia of Genes and Genomes,KEGG)数据库,进行基因生化途径注释,对差异基因调控的生物体代谢网络进行鉴定[16].

1.6 实时荧光定量PCR(qRT-PCR)

采用3叶龄和6叶龄红小豆叶片的总RNA作为模板,分别在昼夜节律、天线蛋白和植物激素信号转导代谢途径中各选择差异性显著的2个基因进行验证,使用定量试剂盒(QuantiFast®SYBR®Green PCR Kit)制备qRT-PCR反应体系,在荧光定量PCR仪(C1000型PCR仪)上进行测定,反应程序为95 ℃ 10 min,95 ℃ 10 s,60 ℃ 30 s,40个循环,2-ΔΔCt计算处理,引物设计使用primer 5.0(表1).

表1 转录组数据qRT-PCR验证引物序列

2 结果与分析

2.1 叶片测序数据质量分析

2个样本经过测序后共得到94.1×106个原始reads(图1A),过滤后得到92.31×106个clean reads(图1B);2个样本的Q30所占比例分别为94.29和94.05%,均值为94.17%(图1C);GC含量所占比例分别为42.94%和44.41%,均值为43.68%(图1D);定位到基因组上的序列数量所占比例分别为97.45%和98.09%,均值为97.77%(图1E);在参考序列上有唯一比对位置的序列数量所占比例分别为92.95%和93.19%,均值为93.07%(图1F).测序结果表明,2个红小豆样本测序质量较好,可以进行后续分析.

A:原始reads数目;B:过滤后clean reads数目;C:Q30占比;D:GC占比;E:定位到基因组上的序列占比;F:参考序列上唯一比对序列占比.不同小写字母表示差异显著(P<0.05).

2.2 差异表达基因的筛选

2个样本共筛选出6 314个差异表达基因,其中,3 059个基因上调表达,所占比例为48.45%,3 255个基因下调表达,所占比例为51.55%(图2A).红小豆2个发育时期上调和下调基因的数目和表达相差较小,上调基因略多于下调基因(图2B).

A:差异基因数目;B:基因表达;灰色表示非显著差异表达基因,红色表示上调显著差异表达基因,绿色表示下调显著差异表达基因.

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

为进一步验证转录组数据的准确性,从昼夜节律(LOC108319708、LOC108322606)、天线蛋白(LOC108328079、LOC108330684)、植物激素信号转导代谢通路(LOC108319999、LOC108320001)中各选取2个差异显著的基因进行qRT-PCR验证,结果与FPKM值趋势相一致,3叶龄的LOC108319708(图3A)、LOC108322606(图3D)、LOC108328079(图3E)、LOC108330684(图3F)的4个基因表达与6叶龄相比显著升高,而LOC108319999(图3B)和LOC108320001(图3C)显著降低.此结果表明红小豆转录组测序数据准确,可靠性较高.

A:LOC108319708基因不同叶龄期;B:LOC108319999基因不同叶龄期;C:LOC108320001基因不同叶龄期;D:LOC108322606基因不同叶龄期;E:LOC108328079基因不同叶龄期;F:LOC108330684基因不同叶龄期.

2.4 差异表达基因的GO分析

GO功能富集分为生物过程、细胞组分和分子功能.通过转录组测序共获得2 933个差异表达基因,其中生物过程中包含585个差异基因,占比为19.95%,富集最多的是防御反应,其次是水分缺失响应、几丁质反应、氧化还原过程等;细胞组分中包含1 461个差异表达基因,占比为49.81%,显著差异基因主要富集在质膜、叶绿体、叶绿体类囊体膜和类囊体功能上;分子功能中包含887个差异表达基因,占比为30.24%,显著差异基因主要富集在DNA结合转录因子活性、血红素结合、铁离子结合和单加氧酶活性等功能上(表2).

表2 差异表达基因GO功能分类

从以上结果可以看出,细胞组分所占比例最大,其中质膜和叶绿体功能上显著差异表达基因最多,分别为678个(其中350个基因上调表达,328个基因下调表达)和467个(其中127个基因上调表达,340个基因下调表达),说明随着红小豆叶龄的增加,细胞组分发挥了重要作用.红小豆叶片的主要功能是进行光合作用,通过基因的上调或下调表达调控细胞质膜对营养物质的吸收,从而促进红小豆的生长发育.

2.5 差异表达基因的KEGG分析

KEGG功能富集分析表明,前20条代谢通路分别是淀粉和蔗糖代谢,植物的昼夜节律,植物激素信号转导,缬氨酸、亮氨酸和异亮氨酸的降解,亚麻酸代谢,脂肪酸降解,色氨酸代谢,丙氨酸代谢,植物MAPK信号通路,柠檬烯和针烯降解,苯丙素生物合成,类胡萝卜素生物合成,半乳糖代谢,DNA复制,油菜素类固醇生物合成,过氧物酶体,光合作用—天线蛋白,丙酮酸代谢,糖酵解/糖异生和亚油酸代谢(图4).

图4 差异表达基因的KEGG富集分析

其中差异性最明显的3条代谢通路分别是淀粉蔗糖、昼夜节律和植物激素,但每个通路中的基因数目和表达情况有所差异.植物激素代谢通路中共有108个基因,其中 74个基因上调表达,34个基因下调表达;淀粉和蔗糖代谢通路中共有75个基因,其中35个基因上调表达,40个基因下调表达;昼夜节律代谢通路中共有31个基因,其中9个基因上调表达,22个基因下调表达(表3).从3条代谢通路中差异基因的热谱图可以看出,随着叶龄的增加,植物激素代谢途径中的基因多数上调表达,而淀粉和蔗糖代谢途径中的上调和下调基因相差不多,昼夜节律代谢途径中的基因多数下调表达(图5).随着植株的生长,激素代谢途径中的基因发挥着重要作用,植株体内的碳水化合物不断增加,激素也不断发生变化,最后通过激素代谢途径中的基因上调表达响应红小豆的生长发育.

表3 红小豆叶片叶龄对淀粉和蔗糖代谢、昼夜节律和激素代谢通路的影响

数值为正、负分别表示基因表达上调和下调.A:激素信号转导代谢通路,通路ID是var04075;B:淀粉和蔗糖代谢通路,通路ID是var00500;C:昼夜节律代谢通路,通路ID是var04712.

3 讨论

本试验对光照比较敏感的晚熟品种‘唐山红小豆’3叶龄和6叶龄叶片进行转录组测序,结果表明,Q30>90%,GC含量>40%,说明测序比较成功.进一步分析表明,通过转录组测序技术共获得6 314个差异表达基因,其中,3 059个基因上调表达,3 255个基因下调表达,GO结果表明细胞组分所占的比例最大,其中,质膜和叶绿体功能上差异基因最多,分别为678和467个,说明红小豆幼苗生长发育过程中叶绿体和质膜功能发挥着重要作用.这在其他作物上也有相关报道;黄睿等[17]研究表明,随着吉祥草的生长,光合作用增强,植株体内的碳水化合物含量较高,维持较低的质膜透性,通过提高酶活性抑制脂膜过强的氧化作用;牛祥芬等[18]的研究也表明,随着烟草的生长,叶片内的超微结构变得复杂,叶绿体基粒片层多,基质浓厚,多片层基粒比例较大;艾叶等[19]通过转录组测序技术对建兰‘铁骨水晶’的叶片超微结构进行研究,结果表明,水晶叶片中与光合作用等代谢途径相关的基因表达量显著下降,进一步说明建兰水晶艺叶形成是由于光合作用相关基因的下调表达造成的.以上研究与红小豆具有相似之处,说明红小豆随着叶龄的增加,通过上调或下调质膜和叶绿体基因调控叶绿体基粒片层和质膜透性,进而调控植株的生长发育.

红小豆3叶龄和6叶龄的GO和KEGG通路富集分析表明,前20条通路中差异最显著的分别为淀粉和蔗糖代谢、昼夜节律和激素代谢通路,每个通路中分别有75(35个基因上调表达,40个基因下调表达)、31(9个基因上调表达,22个基因下调表达)和108个(74个基因上调表达,34个基因下调表达)差异表达基因,其中,激素信号转导代谢途径上的差异基因数目最多,说明激素代谢途径在响应红小豆生长发育过程中发挥着重要作用.热谱图结果表明随着红小豆叶龄的增加,基因表达也呈现动态变化.从以上结果可以看出,随着叶龄的增长,差异明显的通路主要与光照和开花相关.前期试验发现,对不同叶龄红小豆进行遮光后,其开花时间提前[20-21].红小豆遮光10 d后与自然光相比较,KEGG中显著富集的通路包括昼夜节律、激素、淀粉和蔗糖代谢等(论文未发表),只是各个通路中的基因数目不同,上下调表达也不同.这在其他作物上也有相关报道.如:王泽涵等[22]以‘黄棪’茶树品种为试验材料,通过对不同时期茶树花进行转录组分析发现,差异基因主要富集在花的生长发育、激素、昼夜节律和天线蛋白等代谢途径上;卢园萍等[23]通过转录组测序技术对不同发育阶段的蘑菇子实体进行KEGG通路分析,差异基因主要富集在碳水化合物代谢和氨基酸代谢等途径上;唐玉娟等[24]通过转录组测序技术对芒果新梢停长期和基部膨大期的顶芽进行分析,结果表明,KEGG富集最多的通路分别为糖代谢、激素、光合作用和次生物质的生物合成积累等,本试验结果与前人研究具有相似之处,基本趋势一致,只是转录组技术在不同的植物上稍有差异,侧重点也不同.以上结果说明,红小豆叶片是感受光照的主要器官,随着红小豆叶龄的生长,促进其物质积累和开花的通路不断显著富集,与短日照诱导相比,富集的基因数目不同,至于每个通路中的主导基因以及基因之间如何相互作用而促进红小豆生长和开花的分子机制还有待进一步探讨.