共表达网络分析揭示棕榈藤纤鞭生长动态调控模块

王炯亮 赵韩生 高志民

(国际竹藤中心 北京 100102)

棕榈藤属于棕榈科，多年生常绿攀援类植物，世界上共有600多个种，分别属于13个属，其中27个商业品种在热带地区被广泛种植。棕榈藤每年可产生大约70亿美元的收益，超过500万人在经济上依赖于棕榈藤。棕榈藤是热带森林生态系统的组成部分，其纤鞭带有爪钩，利于在森林环境中攀爬和生长；然而，这也使得棕榈藤种植、管理和采收等经营环节变得困难，导致棕榈藤产品成本较高。因此，纤鞭生长发育的遗传调控机制备受关注。作为竹藤基因组图谱(GABR)项目的一部分，虽然国际竹藤中心在2018年发布了单叶省藤和黄藤2个代表性棕榈藤的全基因组序列，但棕榈藤转录组水平上的基因功能分析仍是空白。

基因协同功能网络在构建整体通路模型、细化基因注释，以及在体内模拟重要调控机制等方面发挥着关键作用。作为其中主要代表的共表达网络(CEN)，它是基于基因转录表达水平以表征相关功能基因的强大工具，可在多种条件下批量注释、聚类和探索海量基因。CEN分析和相应的数据库在植物中得到广泛应用，ATTED-II是一个提供了拟南芥、芥菜、大豆等9个物种的CEN数据集的数据库，CcNET提供了基于CEN的二倍体和多倍体棉花功能基因比较分析，AraNet v2构建了提供28个植物物种间同源关系的CEN，RiceNet v2是一个升级后的水稻CEN分析数据库，Bamboo-NET提供了毛竹的CEN和优化后的基因功能注释。但是，棕榈藤中的CEN研究仍处于空白。

一般CEN的构建遵循以下3个步骤：计算基因间的相似性值；选定阈值提取基因对形成网络；在CEN中寻找功能模块。每一步都有不同的方法可供选择，所以构建一个CEN可以有许多不同的组合方法。在本研究中，通过使用单叶省藤和黄藤的基因组和转录组数据，选择成熟方法构建CEN，最大化了基因的覆盖率和模块的数量，优化对基因的功能注释，在CEN中鉴定了覆盖纤鞭不同发育阶段的功能模块，并使用基因家族分类、同源基因注释、顺式作用元件分析和基因本体分析等工具对优化的基因功能注释进行了评估。此外，构建了棕榈藤的CEN数据库(Rattan-NET)，提供了潜在的功能模块和基因功能注释，方便研究人员使用。

1 材料与方法

1.1 样品收集与测序

不同发育阶段单叶省藤和黄藤纤鞭组织采自中国林业科学研究院广州热带林业研究所(原文附表S1)。按照说明书，使用TRIzol试剂分离出RNA，并使用NanoDrop 2000分光光度计测定浓度和纯度。所使用RNase-free DNase I试剂在37 ℃下处理30 min，以去除所提取RNA中残余的DNA，之后将RNA在反转录系统下反转录成cDNA，并建库测序。使用BGISEQ-500平台对富集得到的文库进行长度为100 bp的双末端测序。同时，还收集了此前测序获得的4个黄藤纤鞭组织的RNA测序数据，共使用了处于不同发育阶段的24个单叶省藤和36个黄藤纤鞭转录组数据集进行分析。

1.2 数据处理和基因表达分析

使用FastQC v0.11.6，在默认参数下进行转录组数据的质量分析。使用Trimmomatic v0.36在LEADING:3、TRAILING:3、SLIDINGWINDOWS 4:15、MINLEN:50、TOPPHRED 64等参数下对接头序列和低质量序列进行修剪和过滤。经质量控制后的序列被输入到HISAT2 v2.1.0中进行转录组序列的回帖，参数如下：-min-intronlen 20、-max-intronlen 4000、-rna-strandness RF。HISAT2输出的数据都被保留下来用于后续分析(原文附表S2)。使用Cufflinks v2.2.1进行FPKM值的计算，使用Cuffdiff v2.2.1评估差异表达基因。此外，此研究使用3σ准则(threshold = average(5% value) + 3*SD)计算FPKM阈值，用于筛除FPKM过低的数据(原文附图S1)。

1.3 CEN的构建

本研究采用的策略是将皮尔森相关系数(PCC)和互相排序(MR)整合在一起用于CEN的构建(原文附图S2)：1)使用FPKM计算基因对之间的PCC值，并过滤掉弱相关的基因对，只保留强相关的基因对(原文附图S3)；2)使用PCC计算基因对之间的MR值，并提取高置信度的基因对，剔除不可靠的基因对；3)进行Receiver operating characteristic分析，用于筛选最优的MR参数值。最后，本研究选定了单向MR<3+双向MR<30的基因对作为CEN。

1.4 功能模块的鉴定和基因注释优化

使用基于集团渗透算法的CFinder v2.0.6测定每个节点的密度，进而在网络中鉴定模块。为使基因的覆盖度和模块数量最大化，分别在单叶省藤和黄藤中选定k=6和k=5作为CFinder的参数(原文附图S4)。此外，使用基因集富集分析(GSEA)工具对模块进行功能注释，不可靠的模块在Fisher检验和多重假设检验过程中被剔除。然后，使用模块和CEN数据对2个棕榈藤物种的基因功能注释进行优化，即：用模块的功能注释对属于该模块的基因进行注释优化；对每个基因的共表达基因集进行GSEA，FDR<0.05的词条也做为是该基因优化后注释。

1.5 棕榈藤和拟南芥之间的直系同源基因鉴定

使用BLAST相互最佳匹配(RBH)方法在棕榈藤和拟南芥之中鉴定直系同源基因，排位前三的RBH被认为是最优的直系同源对，所有BLAST结果的E-value值分布的峰值作为阈值，小于该阈值的基因对被认定为次优的直系同源基因对。

1.6 顺式作用元件的显著性分析

参照前人研究，使用基于Z-score和P-value的筛选策略用于顺式作用元件的显著性测试，在扫描棕榈藤基因3 kb启动子区域时，选择P-value<0.05的元件模体作为显著富集的调控元件。

1.7 捕光复合物基因家族的鉴定与系统发育分析

从TAIR网站(https://www.arabidopsis.org)下载拟南芥的捕光复合物(LHC)基因序列，作为查询序列，在单叶省藤、黄藤和水稻基因序列组成的数据库中进行相似性搜索。搜索工具为BLAST+，参数为E-value<1e-5。相似性搜索的结果进一步使用pfamscan.pl脚本进行检验，提取该基因家族的编码序列，用于后续分析。

在系统发育分析中，使用MUSCLE v3.8.31在默认参数下进行多序列比对。使用Gblocks server提取保守区序列。使用jModeltest v2.1.6寻找最优的替换模型，参数如下：-f、-g 4、-i、-s 203、-S BEST。最后，使用PhyML v20120412进行系统发育树的构建，其中模型参数由jModeltest得到，自检参数设为1 000。

1.8 查询与可视化平台构建

Rattan-NET的构建基于Linux系统、Apache软件、MySQL软件和PHP语言。涉及计算相关的分析工具是基于Python、Perl和R脚本语言开发的。Cytoscape.js是Linux版本的Cytoscape，它用于模块与CEN的动态展示。GBrowse和SequenceServer用于提供序列扫描和BLAST+服务。

2 结果和分析

2.1 CEN的构建和功能模块的鉴定

使用处于不同发育阶段的24个单叶省藤和36个黄藤的纤鞭样本进行CEN的构建。分别在单叶省藤和黄藤确定了符合单向MR<3+双向MR≤30条件的630 081和670 502个基因对作为CEN的边(原文表1)。单叶省藤和黄藤的CEN分别包含了31 847和36 769个基因，它们各自覆盖了62.16%和68.93%的基因(原文表1)，2个棕榈藤CEN平均边的数量相近。CEN基因所相连边的数量密度图表明只有很少一部分基因有着高连接度，这符合生物网络的一般特征(原文附图S5A)。

使用基于集团渗透算法的CFinder软件和GSEA工具(阈值为FDR<0.05)在2个CEN中鉴定潜在的功能模块，分别在单叶省藤和黄藤中鉴定得到了3 504和3 027个潜在的功能模块(原文表1)。它们的分布图表明随着每个模块包含基因数量的增大，该规模的模块数量下降(原文附图S5B)。单叶省藤和黄藤功能模块平均包含的基因数为8.1和6.1(原文表1)，这些功能模块使用包含了PO、GFam、GO和KEGG等基因集的GSEA工具进行功能注释。注释结果表明，这些功能模块覆盖了光合作用、植物次级细胞壁生物合成、木质素生物合成、类黄酮生物合成和苯丙素生物合成等在发育中重要的生物过程。进一步分析模块中显著富集的GO词条(Z-score>4)(原文附表S3)，发现了2个棕榈藤物种之间保守的词条，都包括泛素依赖蛋白代谢过程、质子跨膜转运、囊泡介导转运、细胞分裂素响应等。此外，Rattan-NET网站的Download页面提供了由GSEA产生的功能注释文件。共有54.6%的单叶省藤和62.4%黄藤的基因在模块和CEN分析中得到注释优化。

2.2 棕榈藤中类黄酮和木质素生物合成通路的CEN分析

使用RBH方法在2个棕榈藤物种中鉴定参与苯丙素生物合成通路的PAL、C4H和4CL基因(原文图1)，结果表明，单叶省藤中有8个PAL基因、12个C4H基因、9个4CL基因，黄藤中有13个PAL基因、12个C4H基因和15个4CL基因(原文附图S4)。通过检测功能注释，在这些基因中调查与类黄酮生物合成和木质素生物合成相关的共表达基因(原文图1)。结果表明，在单叶省藤和黄藤的上述3类基因家族中分别有12个和14个类黄酮生物合成相关的基因(原文图1B与附表S5)，与发现的类黄酮生物合成相关基因相比，在单叶省藤和黄藤中分别只有2个和5个与木质素生物合成相关的PAL、C4H和4CL基因(原文图1)。由此表明，在棕榈藤的PAL、C4H和4CL这3个基因家族中，类黄酮相关的基因比木质素生物合成相关的基因更多。

2.3 棕榈藤中次级细胞壁生物合成的CEN分析

选择参与木质素生物合成和次级细胞壁(SCW)转录调控的代表性基因，构建它们的CEN(原文图2A、B)。这些基因包括4CL1(Calsi_gene34733和Daeje_Gene51484)、CCR1(Calsi_gene01533和Daeje_Gene15213)和MYB103(Calsi_gene15542和Daeje_Gene04728)。通过分析，在网络中鉴定到参与SCW生物合成的调控因子，在单叶省藤中有MYB20、MYB54和MYB52，而在在黄藤中有MYB55和MYB4，它们可以被回贴到SCW转录调控网络的简化图中(原文图2A、B、C)。因为木质素是SCW中最重要的成分之一，本研究也在网络中鉴定了参与木质素生物合成的相关基因，并将它们回贴到简化图中(原文图2A、B、C)。此外，作为纤维素生物合成酶复合物的核心成分的纤维素合酶基因也在网络中被鉴定出来，包括IRX3(Calsi_gene01341和Calsi_gene30420)、IRX1(Calsi_gene32935)、CESA4(Calsi_gene08952)、IRX6(Daeje_Gene45742)、IRX3(Daeje_Gene59357)、CESA4(Daeje_Gene59357)(原文图2A、B)。随后通过在CCR1和4CL1基因上游3 kb区域的顺式作用分析鉴定得到显著富集(P-value<0.05)的“MYB61”元件(原文附表S6)。

使用包含GO、KEGG和GFam等数据集的GSEA工具，对2个CEN进行分析以检测潜在的生物过程(原文图2D)。SCW生物合成(GO:0009834)、木质素生物合成(GO:0009809)和类黄酮生物合成(map00941)等GO词条同时在单叶省藤和黄藤中被富集，这表明这些生物学过程在2个棕榈藤物种之间具有保守性。此外，在单叶省藤中富集到了纤维素生物合成过程(GO:0030244)、黄酮和黄烷醇生物生物合成(map00944)等词条，而黄藤中富集到了木质素分解代谢过程(GO:0046274)、苯丙素生物合成(map00940)和苯丙素生物合成基因家族等词条。

另外，还鉴定出与SCW生物合成相关的潜在功能模块(原文附图S6)。在GSEA中，SCW合成相关的词条得到显著富集(FDR<0.05)，包括纤维素生物合成过程(GO:0030244)、植物细胞壁生物合成家族、纤维素合成酶类似物，木质素生物合成过程(GO:0009809)、细胞版合成等。

2.4 棕榈藤中光合作用相关基因的分析

在棕榈藤中鉴定得到了LHC基因，并重构了系统发育树以研究棕榈藤中LHC基因家族的进化(原文图3A和附表S7)。系统发育树表明，LHC基因家族可分成13个亚家族(原文图3B)。在单叶省藤和黄藤中分别有21和23个LHC基因家族成员，这与拟南芥中LHC基因的数量(21)类似，比水稻中(15)的多。虽然棕榈藤与拟南芥的LHCB1亚家族成员数量相近，但是其它亚家族成员数量的分布却并不相似。比如，在棕榈藤LHCB6、LHCA2、LHCB5和LHCB3基因亚家族中各有着双倍于拟南芥的成员数量。

通过GSEA分析LHC基因家族的共表达基因(原文附图S7)，发现光合作用相关的GO词条，如光合作用(GO:0015979)、叶绿体类囊体膜(GO:0009535)和光刺激响应(GO:0009416)等显著富集(FDR<0.05)，这与前人研究LHC蛋白的主要功能是通过捕获光能供给光化学反应的结果相符，也表明了2个棕榈藤CEN具有可靠性。同时，通过对棕榈藤LHCA1基因的PCC top300 CEN和拟南芥LHCA1基因的PCC top300 CEN(下载自ATTED-II和AraNet)之间的比较分析发现，2个网络存在高度的一致性，而且还发现一些LHC直系同源基因对，这同样表明网络的可靠性。本研究试着通过分析黄藤中不同LHC基因之间的共表达关系去研究它们之间可能的合作关系(原文图3C)，发现大多数LHC基因关系较近，但属于LHCB1亚家族的4个基因(Daeje_Gene30151、Daeje_Gene53180、Daeje_Gene67755和Dae_Gene18901)游离在其它LHC基因之外。

2.5 Rattan-NET：一个在线的棕榈藤CEN数据库

为提供方便直观的棕榈藤CEN分析功能，构建了Rattan-NET数据库。该数据库整合了CEN分析、顺式作用元件分析、GSEA、GBrowse和Sequenceserver等工具，有助于研究人员在转录水平上对棕榈藤基因功能注释进行精炼。其中，Cytoscape网页版工具提供了CEN和模块的交互式展示功能。在Search页面，用户可以通过输入单个或多个基因的ID进行CEN搜索。Rattan-NET还提供了GO、KEGG和GFam等功能注释的基因ID和关键词搜索功能，在基因详情页提供了每个基因在不同组织下的FPKM值，GBrowse提供了基因组浏览和序列提取功能。总之，Rattan-NET为用户提供了全面的数据和便捷的工具。此外，它会随着分析工具的扩充、棕榈藤物种数据的增加将持续得到维护和优化。

3 讨论

国际竹藤中心的研究团队于2018年发布了2个棕榈藤物种(单叶省藤和黄藤)的基因组序列和它们在全基因组水平上的基因功能注释，这为棕榈藤分子研究提供了重要的基础数据。本研究通过使用上文描述的策略构建CEN，将转录组层面的信息引入到基因功能注释。2个CEN符合无标度网络的特征，表明了它们的可靠性(原文附图S5A)。在单叶省藤和黄藤中分别有3 504个和3 027个，覆盖了多个在纤鞭生长发育过程中重要的生物学过程的模块得到了GSEA鉴定。泛素依赖蛋白代谢过程、质子跨膜转运、囊泡介导转运、细胞分裂素响应等词条是2个棕榈藤物种之间功能模块保守的GO注释。然而，只有62.16%的单叶省藤和68.93%的黄藤基因被CEN覆盖，这低于在毛竹中>90%的CEN覆盖率，这可能是棕榈藤的样本数量过少和单叶省藤与黄藤的基因组scaffold过短导致的。随着对竹藤基因组的维护和GABR计划的推进，将得到棕榈藤中更高质量的基因组、基因注释和源自更多组织的转录组数据集，这将为未来CEN的构建和分析提供全面且高质量的数据。同时，其他棕榈藤物种的CEN将会被构建并进行棕榈藤物种之间的CEN比较分析。

4 结论

本研究利用单叶省藤和黄藤基因组数据和不同发育阶段纤鞭的转录组数据，构建了CEN并优化了基因功能注释，并通过功能富集分析、基因家族分类、顺式作用元件分析和GO分析等证实了基因功能注释和预测的可靠性。通过数据挖掘系统，丰富了对单叶省藤和黄藤重要农艺性状的遗传基础的理解，并对其他重要性状研究提供了新的思路。构建的在线的CEN数据库(Rattan-NET)，将有助于棕榈藤CEN的应用，对促进棕榈藤分子生物学发展，揭示棕榈藤重要性状分子调控机制具有重要科学价值。

原文出处

Wang J, Ma X, Yang J, Hui Y, She J, Tian T, Li Z, Xu W, Gao Z, Su Z, Zhao H. Coexpression analysis reveals dynamic modules regulating the growth and development of cirri in the rattans (CalamussimplicifoliusandDaemonoropsjenkinsiana).Frontiers in Genetics, 2020,11:378. DOI: 10.3389/fgene.2020.00378. eCollection 2020.