刘嘉智,刘长宁
(1.中国科学院西双版纳热带植物园,云南勐腊 666303;2.中国科学院大学,北京 100049)
橡胶籽是天然橡胶种植生产中的重要副产品,橡胶籽的含油量为35%~50%,年产量可达2 060 kg/hm2[1],与其 他 非 食用 植 物 油相 比,橡胶籽油燃烧的CO2和NO2排放率更低,是理想的生物柴油来源[2-3]。然而,目前对橡胶树的利用主要集中于获取天然橡胶,忽略了橡胶籽作为生物柴油的巨大潜力,与之相关的研究也并不充足。
通常,具有同源序列的蛋白编码基因具有相似的功能,因而基于同源分析可以在近缘物种中预测相关功能基因。本研究基于共线性分析和同源分析,利用目前研究中较为深入的蓖麻全基因组关联分析数据,将蓖麻野生型品种中种子农艺性状相关的基因映射到橡胶树reyan7-33-97 品种基因组上,使用生物信息学手段预测橡胶树中种子性状相关基因并进行基因功能分析,为橡胶籽的合理开发利用提供理论参考。
蓖麻野生型品种的基因组文件以及基因组注释文件下载于http://oilplants.iflora.cn/Download/castor_download.html,橡胶树reyan7-33-97 品种的基因组文件以及基因组注释文件从https://www.ncbi.nlm.nih.gov/data-hub/genome/GCF_001654055.1/下载。从上述文件分别提取基因ID、序列信息、基因组位置信息和功能描述等内容。从NCBI SRA 数据库中收集具有文献来源且为双端测序的蓖麻和橡胶树根、茎、种子3 种组织的转录组数据,样本编号分别为SRR3136156、SRR6127585、SRR3136168、SRR17218208、SRR17218209、SRR17218216。
1.2.1 橡胶树种子性状相关基因预测
根据文献调研[4],在蓖麻野生型品种中收集基于全基因组关联分析得到的种子特征性状相关的SNP 位点。根据基因组注释信息,在上述SNP 位点前后50 kb 范围内提取基因,并认为这些基因可能与种子特征有关。
橡胶树与蓖麻同属大戟科,具有较近的亲缘关系,因此,首先判断蓖麻种子性状相关基因所在区块在橡胶树中是否存在共线区块。由于橡胶树reyan7-33-97 品种基因组只组装到scaffold 级别,使用BLAST 进行共线性分析后认为存在3 个及以上基因,且基因之间相互间隔不超过50 kb的片段具有共线性。然后,使用Orthofinder 软件[5]对这两个物种进行同源组鉴定,得到的结果再根据双向BLAST 最优策略进一步确定同源基因,并根据基因结构验证预测的同源基因是否可靠,最终将收集到的蓖麻野生型种子性状相关基因比对到橡胶树基因组上,来挖掘橡胶树中与种子特征性状相关的基因。
1.2.2 表达谱分析和差异表达分析
下载蓖麻和橡胶树根、茎、种子3 种组织对应的转录组原始测序文件,使用Trim Galore 进行质控检测去除低质量碱基和接头。利用HISAT2 软件[6]对基因组序列构建索引后进行序列比对,使用Stringtie 工具[7]进行定量分析来获取基因水平上的Raw counts、TPM 值以及FPKM 值。将TPM 值进行归一化处理,使用R 包pheatmap 绘制表达量热图,并对上述两个物种得到的表达量结果进行比较。使用GFOLD 1.1.4 软件[8]对橡胶树基因进行差异表达分析,进一步挖掘橡胶树中与种子性状相关基因在种子与根、茎之间的表达差异。
1.2.3 miRNA 靶位点预测
根据文献调研[9]以及pmiren 数据库[10]中获取橡胶树的成熟miRNA 序列,使用gffread 软件提取转录本序列,使用psRNATarget 在线工具[11]默认参数进行miRNA 靶基因预测,从预测结果中筛选期望值大于3.5 的相互作用的miRNA 和靶基因,使用Cytoscape 软 件[12]将miRNA 与 靶 基因 间的 网络关系进行可视化展示。
1.2.4 蛋白互作网络预测
基于同源组划分,以拟南芥蛋白互作网络作为参考映射到橡胶树中,获得橡胶树reyan7-33-97 的蛋白互作网络。从String 数据库[13]中下载拟南芥蛋白互作物理子网,使用Orthofinder 软件预测拟南芥与橡胶树reyan7-33-97 之间的蛋白同源群组。使用Cytoscape 软件对橡胶树种子性状相关候选基因进行网络可视化。
在野生型蓖麻中共发现18 个与种子特征性状相关的SNP 位点,这些SNP 位点分别与种子面积(SA)、种子长度(SL)、种子厚度(ST)、种子宽度(SW)以及单粒重(SSW)等5 个特征形状相关。根据蓖麻基因组注释文件,在上述SNP 位点上下50 kb 范围内共有67 个基因(图1),分别位于蓖麻的1、2、3、5、6 以及9 号染色体上(表1)。其中Rc01G001599、Rc01G001600、Rc01G001601、Rc01 G001602、Rc01G001603、Rc01G001604、Rc01G001605、Rc01G001606、Rc01G001607、Rc01G001608、Rc01G 001609 对上述5 个性状都产生影响。
表1 蓖麻野生型种子特征相关基因
图1 蓖麻野生型种子性状相关基因统计
基于共线性分析,在蓖麻野生型与橡胶树reyan7-33-97 中发现4 个共线区块,分别位于蓖麻的2、5、6 以及9 号染色体上,对应橡胶树的NW_018746954.1、NW_018745856.1、NW_018746236.1 以及NW_018746078.1 这4 条scaffold,共 覆 盖19 个基因(表2)。根据Orthofinder 软件,在蓖麻野生型与橡胶树reyan7-33-97 中共鉴定到16 926 个同源组,根据双向BLAST 最优确定二者之间序列相似程度,上述67 个蓖麻基因最终在橡胶树中比对得到49 个基因。对上述橡胶树基因以及同源蓖麻基因根据基因结构作图进行验证(图2),发现基因结构与最终得到的基因序列预期结果相符。
表2 蓖麻野生型与橡胶树reyan7-33-97共线区块基因统计
图2 蓖麻野生型与橡胶树reyan7-33-97同源基因结构图
为进一步探究蓖麻和橡胶树种子性状相关基因的功能,进行了组织表达谱的分析,图3 是蓖麻和橡胶树种子性状相关基因的组织表达谱。从中发现,有些同源基因在3 个组织中的表达量类似,但有些差别很大。对蓖麻和橡胶对应的49 个基因在种子中表达结果进行统计,发现在两个物种中有5 个基因均为正表达,25 个基因均为负表达,19 个基因表达相反,这可能是由于物种在分化过程中部分基因功能改变导致的。
图3 蓖麻和橡胶树种子性状相关基因组织表达谱
根据差异表达分析,选出16 个在种子中表达量远高于根茎的基因、17 个在种子中表达量远低于根茎的基因以及5 个在三者之中都没有差异表达的基因,其余基因在3 个组织中无交集(表3)。
表3 种子与根茎差异表达基因
基因的miRNA 靶位点预测也是研究基因调控的重要途径。本研究共收集到橡胶树的217 条成熟miRNA 序列,对橡胶树种子性状相关基因进行miRNA 靶位点预测,结果发现LOC110671483、LOC110673189、 LOC110645416、 LOC110651200、LOC110672934、 LOC110641727、 LOC110653311、LOC110659520 受到不同miRNA 调控(图4)。由于miRNA 在物种进化中相当保守,大部分miRNA 在特定的组织和发育阶段表达,我们预测得到的miRNA 许多功能都与种子性状相关。例如之前研究发现miR156 参与调控拟南芥种子成熟过程[14]、miR408 在水稻种子发育中具有重要作用[15]等,因此推测上述8 个基因可能受到miRNA 调控影响种子发育,从而影响了种子特征性状。
图4 橡胶树种子性状相关基因miRNA靶位点预测
基于同源组的划分,橡胶树reyan7-33-97 品种中共有22 815 个基因映射到19 762 个拟南芥基因。根据预测得到的橡胶树蛋白互作网络,我们将橡胶树reyan7-33-97 品种中与种子性状相关的49 个基因放入蛋白互作网络中,发现除LOC110654127 和LOC110668244 外,其 它 基 因 之 间不存在边。这可能是影响种子性状的橡胶树候选基因与其他基因也存在相互作用导致的。为此,提取这49 个基因以及最近的相邻基因重新进行网络构建,结果发现其中16 个橡胶树种子性状候选基因通过一级相互连接构成网络(图5),因此认为这16 个基因存在间接蛋白相互作用从而影 响 种 子 性 状,LOC110654127 和LOC110668244 同时存在直接蛋白相互作用影响种子性状,而其它33 个基因可能单独影响种子性状或者存在不同于拟南芥的蛋白互作成员。
图5 橡胶树种子性状候选基因蛋白互作网络
本研究使用全基因组关联分析、基因组共线性分析和基因序列同源分析等多种技术手段,预测了49 个与橡胶树种子性状相关的基因并进行了多种组学分析,研究结果为进一步开发利用橡胶籽提供了一定的工作基础。然而,橡胶籽种子性状基因与种子产油之间的关系是一个复杂的过程,植物激素的生物合成和信号传递对种子大小和形态的控制至关重要,例如赤霉素可以促进种子的扩大和发育[16],而乙烯则可以抑制种子的生长[17]。因此,后续研究可深入挖掘橡胶籽种子性状基因和植物激素生物合成和信号传递通路之间的关系,通过对橡胶树中这些激素合成和信号传递通路进行深入分析,揭示橡胶籽种子性状基因在种子性状形成过程中的作用机制,以便更好地理解种子大小和形态的调控机制,为橡胶树的育种和生产提供指导。