海松发育胚胎转录组的综合组装与分析

背景

有被子植物和裸子植物之间在胚胎发育明显的差异。大多数植物胚胎发育过程中基因表达与调控当前所知的已经从研究中得出的被子植物物种，特别是从模型厂拟南芥．少数发表的关于针叶树胚胎发生的转录本谱的研究表明，存在许多假定的胚胎特异性转录本，但没有指定的功能。为了扩展对针叶树胚胎发生中转录组表达的认识，我们对几个发育阶段的合子胚胎的转录组进行了测序，覆盖了大部分发育阶段松果体松树胚胎发生(海上松)。

结果

从五个合子胚的发育阶段收集的总RNA样品与Illumina公司的技术进行测序。从头转录组组装，因为没有基因组序列尚未公布松果体松树．海松合子胚胎发生时期的参考转录组包含67429个转录本，可能编码58527个蛋白质。注释显示，31%的预测蛋白只存在于松树胚胎发生中。差异表达转录本的功能分类和富集分析表明，在早期胚胎阶段碳水化合物运输和代谢过度，许多推定的糖苷水解酶的鉴定突出了这一点，可能与细胞壁修饰和碳水化合物运输转录本有关。此外，染色质重塑事件的优势在胚胎发生的早期到中期被检测到，与组蛋白及其与转录增加相关的翻译后修饰物的活跃合成，以及转座子的沉默有关。

结论

我们的结果延长了针叶树术期间对基因表达和调节的理解，并是一种有价值的资源，以支持针叶树种类植被繁殖的体细胞胚胎发生的进一步改善。与碳水化合物代谢，单糖传输和表观遗传调节相关的特异性转录物似乎在杉木早期胚胎发生中发挥着重要作用，并且可能是早期胚胎发生的可靠分子标记的来源。

在更高的植物中，胚胎发生从受精卵形成开始，并理解整个发育过程，导致由种子组织包围的完全成熟和休眠胚胎[1］．目前关于植物胚胎发生的大部分知识来源于被子植物物种的研究，特别是模式生物拟南芥(综述[2])。然而，裸子植物和被子植物谱系估计在3亿年前就分离了[3.]和它们的差异，特别是在胚胎发生阶段，是公知的（通过审查[4])。在裸子植物中胚胎发生的分子研究，特别是针叶树，在过去几年中获得了兴趣（[5，6])。这已经由如何在被子植物和裸子植物之间胚胎发育的特性差异在分子水平上建立更好的理解，和他们进化的影响带动。此外，体细胞胚胎发生的进一步改善，一个有吸引力的技术，经济上重要针叶树大规模无性繁殖，在很大程度上依赖于对基本流程控制胚胎发育更多的知识。

应用于mRNA发现和分析(RNA-seq)的下一代测序(NGS)技术已被证明对研究植物基因调控非常有用，特别是对于仍然缺少参考基因组的非模式物种(由[7，8])。单独的1KP项目等大型RNA测序项目达到了超过1000种不同植物物种的转录组测序，这代表了在系统发育和土地植物演化研究中具有高影响力的巨大努力[http://www.onekp.com；［9]]。RNA-seq数据极大地促进了我们对植物胁迫反应调控的认识[10，植物发育[11，12，13]、商业/生物技术相关植物产品的合成[14甚至是特定基因的进化[15］．几种植物物种，包括水稻的胚胎基因转录（栽培稻) ［16),玉米玉米) ［17，18)、菜籽油(芸苔栗鸟) ［19]和萝卜（Raphanus SativusL.）[20.，都是通过RNA-seq技术生成的。在针叶树中，苏格兰松(Pinus sylvestris.) ［21]主要使用高通量测序来分析主导胚胎。此外，其他研究使用NGS寻址裸子植物胚胎发生（通过[22)，包括Yakovlev et al.(2014)关于胚胎转录组变化的报告挪威云杉在不同温度条件下[23在胚胎发生和非胚胎发生组织的转录组上云杉balfouriana［24关于早期体细胞胚胎形成和种子发育的转录组对比分析南洋杉angustifolia［25的一项全面的转录组调查松果体lambertiana组织类型包括胚胎[26的体细胞胚转录组分析挪威云杉和[27]上的具有相关联碳水化合物介导的应答的识别南洋杉angustifolia胚胎形成。

以前，时间课程转录组研究松果体松树指出胚胎发育过程中表观遗传调控因子与特异性转录因子的相关性[28］．在这种研究中，遵循跨物种微阵列杂交方法，限制了候选转录物的鉴定到来自的阵列探针Pinus Taeda.根和针组织。在该研究的基础上，我们通过使用高通量测序方法扩展了转录组分析的范围，该方法比微阵列具有已知的好处(在[29])，包括检索新的和/或低表达转录本的能力，或可能被微阵列分析遗漏的替代剪接变体[12］．在这项工作中，我们产生了旨在具有在胚胎发生期间表达的胚胎发霉转录物的显着扩展胚胎目录的特异性转录组，靶向胚胎发育期间表达差异较高的转录物。通过遵循这种方法，在胚胎发育的特定阶段以及所涉及的特定转录物中，我们已经确定了过度代表的过程，即碳水化合物代谢和表观遗传调节。海洋松木中的这种独特资源进一步有助于深化针叶草中胚胎发生期间转录活动的知识。

全面的转录组组装

根据先前报道的分期系统，我们对5个胚胎发生阶段(Day0, Day5, Day11, Day15和Day25)进行了rna测序[30.，覆盖合子胚的整个发育时期p .松树，直至成熟阶段。为了捕捉胚胎的转录组景观和整个发育过程中基因表达的主要差异，利用Illumina短读技术对RNA-seq文库进行测序。在没有公布参考基因组的情况下p .松树，一种全面的转录组组装方法([31];综述了在32])，它使用Trinity(版本2.0.6)将重新组装的读本组合在一起[33以及通过将读取的文本和重新组装的文本映射到一起的向导汇编p . taeda基因组(版本1.01)[34］．共获得约319 M读对，平均每个样本有63.8 M读对，每个样本的读对从44.9 M(14.1%)到75.4 M(23.6%)不等。5.3%的读取通过过滤步骤被删除。由此产生的约302个读对随后被用于从头转录组组装和映射到p . taeda基因组。在获得完整的转录组组装后，在正确的距离和方向上共绘制了183.4 M读对(正确配对)，用于表达式分析，这代表57.5%的原始读取(表1)．

Maritimime Pine Zygotic胚胎发生的最终参考转录组含有67,429转录物（在加入GGEX01000000下沉积在DDBJ / ENA / Genbank，https://www.ebi.ac.uk/ena/data/view/GGEX01000000)，长度从148 bp到12,752 bp不等，平均长度为999 bp(表2)．

转录组注释

所有的ORF可能性都是从新组装的转录组中产生的，每个转录组只翻译一个(多个可能情况下最长的一个)，生成最终包含58,527个蛋白的蛋白质组(附加文件)1)．使用两个信息来源进行组装的发育胚胎转录组的注释：最佳的互易次数（BRH）到蛋白质Omp . taeda，p . lambertiana和答:芥，和NCBI数据库中蛋白质的同源性(附加文件2)．从注释的BRHs的功能和基因本体论(GO)术语被分配回p .松树蛋白质，以注释结果14,211p .松树蛋白质。从转录组的BLASTX比对到NCBI非冗余蛋白数据库(e值< 10)开始，使用Blast2GO(3.1版)注释转录本⁻³额外的文件3.)．近30000序列与至少一个GO项注解，并围绕7500成绩单曾在NCBI数据库中的同源性，但无法清楚地关联到GO项。有28780，26585和24241的转录与至少一个F：GO，P：GO或C：GO术语，分别。有16,056成绩单从每个三类至少一个GO项。超过20,000p .松树蛋白质序列没有与数据库中的任何蛋白质对齐(附加文件4)．同源蛋白质呈现为71％的平均相似性，而54％的爆炸点击率为70％以上的序列相似性（附加文件5)．分析得分最高的相应序列p .松树记录显示,西加云杉的近亲p .松树是迄今为止最具代表性的物种，有16000多种爆炸顶部命中．另一种同源属于来自不同植物分类群的物种，包括裸子植物，植物植物和苔藓。四个其他松树物种被代表，即p . taeda，p .放射虫纲，P.樟子松和p . monticola(附加文件6)．与ebi的蛋白质序列和函数的eproproxto数据库（IPS）的比较显示约37,500p .松树具有IPS结果的序列;IPS分析贡献了GO条款到超过20,000p .松树序列(额外的文件7)．映射的结果表明，和的UniProtKB TAIR是GO术语的两个主要来源数据库相关联，以p .松树序列(额外的文件8)．

胚胎发育过程中的功能调节

的蛋白质松树、针叶松、蓝柏根据EMBL的功能注释蛋白eggNOG数据库中各自最佳同源序列的eggNOG组聚在一起(附加文件9)．在比较三种物种的群体数量时，p .松树有最多的独家团体（4355）。仍然，大多数群体总共有5698种来自三种物种的蛋白质（图。1)．具有更多蛋白质成员的群体（附加文件10)是三种针叶树共有的，并被注释为包含“pentatricopeptide - ppr - repeats”，或涉及复制，重组和修复或未知函数“终端逆重复-TIR-”，“富含亮氨酸的重复序列”或“NB-ARC结构域”，它是在牵涉的分子开关信号转导机制．另一方面，这些团体独家p .松树其中至少有10个蛋白被注释为“锌指蛋白”(19个成员)，“姐妹染色单体聚合蛋白PDS5”(13个成员)，“gpd - l -半乳糖磷酸化酶”和“锌离子结合蛋白”(各11个成员)，还有几个簇被注释为“逆转录转座子蛋白”。然而，大多数团体都不包括p .松树只包含一种蛋白质。

因为每个蛋白质组都被分类为一个功能性的eggNOG注释[35]，它是可以将功能信息冷凝为24个类别。大多数预测的蛋白质和相关eggNOGs有未知函数(共有14134个蛋酒，其中5684个来自p .松树)(额外的文件10)．对于剩下的，整个功能类别的热线图p .松树胚胎发育(图2）显示在早期（Day0），中间（Day5和Day11）或晚期（Day15和Day25）胚胎发生阶段的三个主要表现功能类别。在早期胚胎阶段，与之相关的功能细胞骨架，能源生产与转换，碳水化合物运输和新陈代谢，氨基酸转运和代谢和细胞内运输、分泌和囊泡运输中普遍存在，与前三只被这些阶段上调。相比之下，复制，重组和修复,细胞周期控制，细胞分裂，染色体分裂,染色质结构和动力学在后期胚胎阶段遵循相反的趋势，并上调。

关于发育胚胎的整体代谢活动，碳水化合物运输和新陈代谢在早期胚胎阶段上调，在成熟胚胎阶段降低，氨基酸转运和代谢在早期胚至早期子叶胚阶段占优势，而次级代谢产物生物合成、运输和分解代谢在成熟胚胎阶段达到高峰。核苷酸输送和新陈代谢和辅酶运输和新陈代谢在成熟胚期表达量最高，而在成熟胚期表达量最高脂质运输与代谢在先子宫内和成熟的胚胎阶段都是相关的。

沿胚胎发育的差异表达转录物

从67,429个组装的转录本中，64,766个在至少一个发展阶段明确表达，39,838个在所有五个发展阶段都表达(附加文件)11)．每个阶段表达的转录本数量相等，每个阶段超过10,000个表达的转录本。共有4953个转录本(7.3%)仅在胚胎发育的一个阶段表达，其中848个、338个、282个、1703个、1782个转录本分别在Day0、Day5、Day11、Day15和Day25特异表达。

一种differential expression analysis (FDR < 0.05) between each pair of consecutive stages identified 1738 transcripts (2.6%) differentially expressed in at least one transition (Additional file12)．在第1 (Day0 ~ Day5)、第2 (Day5 ~ Day11)、第3 (Day11 ~ Day15)和第4 (Day15 ~ Day25)转换中，差异表达的转录本分别为798 383、591和568个。差异表达转录本的上调在第一次转化中占主导地位，1185个上调的差异表达转录本中有594个上调(图1)。3.）对于第一次转变，而在最后过渡的差异表达转录物的下调更丰富（图。4)，约占全部下调转录本的38%。每一个发育转变与连续的转变共享很少的上调或下调转录本，支持本研究中选择的发育阶段在持续的转录活性方面有明显的差异。

通过对差异表达的转录本进行基因富集分析，获得了与第一阶段和最后阶段过渡相关的过度表达的GO术语(见附加文件)13，14和15)．这两个过渡共享相关的下调成绩单的富集细胞周期的调节和细胞分裂．此外，第一过渡似乎是由参与转录的下调来表征单糖转运和植物型细胞壁纤维素代谢．最后一个转变的特点是与表观遗传学相关的转录下调和注释DNA依赖性DNA复制术语。至于上调转录本，在第一次转变中发现的转录本与萜类化合物分解代谢和叔醇代谢，而不仅丰富GO术语最后过渡发现是分子功能营养物储存活动(数据没有显示)。总的来说，术语碳水化合物运输和代谢(表3.)和表观遗传学相关术语(表4）在差异表达的转录物列表中高度代表。

在差异表达的转录集群

通过K -means聚类分析，1738个在胚胎发育过程中表达差异的转录本可以被分为12个不同的表达模式相似的簇(K)(图)。5)．总体而言，所述簇可进一步在四组根据胚胎阶段划分（A至d），其中表达的峰：（A）早期胚胎阶段（簇1，2，3，9），（B）早期子叶和/或子叶阶段（簇6，8，10，12），（C）预子叶胚胎阶段（簇4,5），和（d）成熟胚阶段（簇7，11）。以下在每个集群执行的富集分析，几个过程和功能被发现在（图胚胎发育的不同阶段过表达。5)．A组的转录本丰富了几个GO术语，包括碳水化合物新陈代谢，单糖转运，固醇代谢，细胞壁组织或生物发生，DNA依赖性DNA复制，细胞周期，基因表达调控，表观遗传(图)5)．在集群2中，与碳水化合物新陈代谢和单糖转运和细胞壁的高度代表(附加文件16)．在集群3中，其与簇2不同，主要是因为表达在中胚阶段的中胚阶段保持相对较高，然后朝着成熟的胚胎阶段降低，DNA依赖性DNA复制是否与贡献成绩单的数量(89份)特别相关染色质的组织，DNA修饰，DNA甲基化，组蛋白赖氨酸甲基化，DNA包装，组蛋白甲基化，基因表达调控，表观遗传(附加文件17)．B组的转录本富含氧化石墨烯，例如类固醇脱氢酶活性，倍半萜类分解代谢，种子休眠释放和叔醇代谢．在组d，最具代表性的条件是营养物储存活动，碳氧裂解酶的活动，的α-红没药烯合酶活性和萜烯合酶活性．

验证通过qPCR

根据表达谱和碳水化合物代谢的推测，我们选择了八个转录本的子集3.）或表观遗传调控（表4)，采用RT-qPCR独立验证RNA-seq结果(引物信息见附加文件)18)．其中5个基因在胚胎发育过程中有差异表达(FDR < 0.05)，并包含在不同的转录谱中(簇2、3、5、10)。RT-qPCR和RNA-seq的表达结果基本一致，Pearson相关值介于0.52 (Pp34678)和0.93 (Pp38781)之间，证实了特定发育阶段的表达峰值(图3)。6)．对于未被预测为差异表达的转录本(Pp29536, Pp34388和Pp34678)， RT-qPCR和RNA-seq也有很好的一致性。

在这项工作中，RNA-SEQ用于提供海洋松树发育胚胎的转录组的全面概述。以前，使用与此处分析相同的发育阶段的海洋松树胚胎的第一次转录组分析依赖于使用储质松树（p . taeda）微阵列，其在两个物种[之间是保守的胚胎发生过程中鉴定的转录本36］．虽然微阵列已成功用于交叉物种杂交[37，38，它包含大约25000个独特的cdnap . taeda这限制了我们的分析范围，因为我们主要关注胚胎发育。在目前的研究中，这种限制被克服了，并且在胚胎发生期间的转录组景观被显著地扩展了。利用综合策略，从5个发育阶段的胚胎中获得的300多对阅读对进行组装，生成了海松合子胚胎发生的参考转录组，包含67,429个转录本。这个数字是三倍的表达转录本的数量，已经被鉴定为p . taeda微阵列(36]，对应于引燃蛋白数量的类似增加。尽管如此，西加云杉，葡萄和里纳斯市政府在所有顶级击中的分布中仍然是前三种最代表的物种。

的注释p .松树，p . taeda和p . lambertiana预测蛋白及其最佳同源物的eggNOG数表明，大多数注释蛋白是三个物种共有的(40.3% inp .松树51.1%,p . lambertiana和86.2％，中p . taeda)．p .松树和p . lambertiana预测蛋白质组共同具有大约10倍以上的蛋白比用共享p . taeda而后者与两种物种的蛋白质含量大致相同。此外，独家预测蛋白的百分比较高p .松树(30.8%)，低于3.0%p . taeda．同源组分析突出了三个物种共有的部分eggnog，其中蛋白成员数量较多，包括“PPR repeat”和“Retrotransposon protein”参与复制，重组和修复或用未知函数，几个被描述为“亮氨酸丰富重复(LRR)”和“带有ARC基序的核苷酸结合域(NB-ARC域)”的蛋酒参与信号转导机制和“toll/interleukin-1 receptor (TIR)”未知函数．这些结果与报道的芽转录组注释密切一致孔雀草松其中发现的最大家族包括1794个成员，包括LRR、TIR、NB-ARC、高尔基转运复合物5 (COG5)和痘病毒A32蛋白基元[31］．

早期胚胎发生中的碳水化合物代谢和运输

这项工作聚集了几个证据，无论是在整个胚胎发育中所代表的功能类别的分析，还从分析特定转录物的分析，指向碳水化合物运输和新陈代谢的突出作用p .松树胚胎发生。对蛋酒的功能分类进行分析，结果显示，蛋白的表达上调碳水化合物运输和新陈代谢在胚胎早期和胚胎发生后期的表达减少。同时，对所有差异表达的转录本进行富集分析，发现一组Day0(聚类2谱)达到峰值的高表达转录本富集于碳水化合物新陈代谢， 包含细胞碳水化合物代谢．在同一份档案中，细胞碳水化合物代谢似乎与初级细胞壁成分纤维素、果胶和葡聚糖的代谢有关。另外，在Day0到Day5期间下调的204个差异表达的转录本富集单糖转运和植物型细胞壁纤维素代谢(见附加文件13)．最近的一项研究P.樟子松种子发展也指出了一个过度代表碳水化合物代谢过程和细胞壁修改在早期发育阶段过度表达的差异转录本中的术语[21］．

许多糖苷水解酶编码的转录本在第一胚胎期(聚类2)表达较多，包括AT4g02290的Pp38781同源物(糖基水解酶9家族)，其表达谱已通过RT-qPCR成功验证。特别地，在这个表达谱中存在几种松树同源物的α-木糖苷酶和木葡聚糖内转糖基酶(以前命名为木葡聚糖内转-β-葡聚糖酶)，指出了木葡聚糖动员和/或降解在松树早期胚胎发生中的重要性[39］．Xyloglucan是除草外的初级细胞壁的主要半纤维素组分（由[40]），能够在细胞壁膨胀期间在木葡聚糖的周转或再循环中调整木糖葡聚糖侧链的糖磷酶（在[41])。拟南芥α-木糖苷酶1 /改变的木葡聚糖3 /热侵蚀萌发1（XYL1 / AXY3 / TRG1），其推定的海洋松同源物（PP8434）在早期胚胎发生中也上调，参与木葡聚糖降解到游离单糖中的酶的核查编码[42和唯一对木葡聚糖有活性的α-木糖苷酶[43］．拟南芥的发芽种子xyl1由于原发性细胞壁完整性的变化，函数损失突变体展示细胞壁松动和降低的种子休眠[44］．也是拟南芥推定的同源物(Pp38450)木葡聚糖内转移酶a3 /木葡聚糖内转移酶/水解酶27（EXGT-A3 / XTH27），参与气管元素通过沉积在细胞壁中的Xyloglucan的降解的分化[45，用相同的表达谱（簇2）识别，以及推定的同源物XTH5，在胚胎轴和种子胚胎中表达的嗜酸甘油蛋白（Ga） - 培养基，参与发芽[46］．综上所述，在萌发过程中参与细胞壁修饰的糖苷水解酶对松树胚胎发生的早期控制也有一定的作用。事实上，负责木葡聚糖降解的酶已经被发现与不同物种的体细胞胚胎发生诱导有关(由[47，48]），包括针叶树如p .放射虫纲其中检测到α- d -半乳糖苷酶(SEPR1)上调[49］．在挪威云杉体细胞胚胎发生A.Xyloglucan: xyloglucosyl转移酶和一个糖基转移酶家庭1蛋白发现在早期胚胎从增殖到分化以及从早期胚胎到晚期胚胎发育的转变中差异表达和下调[38］．

在我们的研究结果中发现了糖基水解酶家族的其他成员的差异表达，包括推定的同源物(Pp46170)细胞壁转化酶1（CWINV1 / CWI1)，被认为在植物库组织中起作用，它将蔗糖通过质质体从韧皮部卸下，水解成果糖和葡萄糖(在[50])。棉花和拟南芥种子的研究显示CWINV在整个球形胚胎中都有表达，但似乎仅限于鱼雷胚胎的中央区域[51］．比较了成株茎原基对体细胞胚诱导的响应云杉glauca，上调CWINV1仅在非应答基因型中检测到，并讨论了与生物应激反应的可能联系[52］．也有许多证据表明cwv介导的蔗糖水解与早期胚胎发生中促进细胞分裂的糖信号有关(由[50，53])。

许多推定碳水化合物转运体的转录本编码也在第一胚胎期(Day0)达到高峰，其中大多数为己糖(或单糖)转运体，参与运输和摄入细胞质，可能是在CWIN切割了存在于外质体中的蔗糖后[51］．虽然仍然很少是众所周知的这些己糖转运蛋白在早期植物胚胎发生中的作用（在[54]），对其他库组织以前的工作，像发育中的胚胎，也symplastically分离，已经表明它们对转运蛋白依赖性摄取通过质膜的蔗糖衍生物（通过审查[55])。拟南芥糖转运蛋白14 (ATSTP14)，属于AtSTPs家族，但特异于半乳糖运输，在源组织(绿叶)和库组织(种子胚乳和子叶)中均有表达[56，也被检测到。有几个证据表明AtSTP14在细胞壁循环中发挥作用，即在β-半乳糖苷酶(糖基水解酶)和其他酶降解细胞壁时释放的细胞壁源半乳糖的运输[56］．另一种推测的拟南芥松树同源物(Pp32670)ATSTP13在簇2中发现，可能参与从细胞壁中提取的单糖的运输(在[54])。

来自簇2和3的富集产生，其聚集几种糖基水解酶，碳水化合物转运蛋白和激酶，以及细胞周期相关的转录物，似乎支持碳水化合物代谢和细胞周期调节剂之间的主动沟通p .松树早期胚胎发生。到目前为止，营养状态和遗传控制相互作用对植物细胞周期控制和细胞分裂的影响仅在胚胎后期发育研究中得到阐明[57］．在拟南芥中显示的意见，即碳水化合物可用性由糖信令与特定按键的细胞周期调节剂，例如CYCB1相互作用播放一段G2 / M过渡对细胞决定作用; 1和CDKB1;分生组织的1，这直接影响增殖58］．在五细胞周期蛋白依赖性激酶b在簇2和3之间发现的转录物存在PP42651，拟南芥周细胞周期蛋白依赖性激酶B2; 2（CDKB2; 2），其已被证明是细胞周期进展和SAM组织的调节因素，并参与激素信令[59］．

胚胎发生早期到中期的表观遗传学相关转录本

与表观遗传调控的不同组件相关的许多成绩单被发现沿海松树胚胎发育过程中差异表达。他们中的大多数被包括在第三组简档，在早期胚胎发育，其特征在于表达的峰值。

5个假定的组蛋白亚基同源物在松树胚胎发育过程中有差异表达，并遵循不同的表达谱。Pp36206和Pp38724分别编码核心组蛋白H3和H2A亚基，在胚胎发生的早期表现出较高的丰度，在胚胎发生的后期普遍下降(第二个小峰出现在第11天)。也在挪威云杉体细胞胚胎发生A.组蛋白3发现在胚胎发生后期有差异表达和下调[60］．组蛋白H3是已知的磷酸化靶点，在细胞周期依赖的方式下，所有三种拟南芥极光激酶[61，62］．有趣的是，Aurora-2的松树同系物(Pp32543)也存在于同一簇中。此外，推测由Pp46360、Pp14332和Pp46359松树转录本编码的H1(连接子)组蛋白在D11/D15阶段(聚类6或10)出现表达高峰，已被报道参与DNA甲基化和去甲基化、细胞周期进程和植物发育(由[63])。体细胞胚胎发生材料的微阵列分析挪威云杉揭示了上调的组蛋白H1在胚发生前的质量（的PEMs）植物生长调节剂，当PEM到胚胎转变诱发的戒断后一天，并在过渡从增殖到早期胚胎的分化〔38］．

染色质形成或重塑也是相关的基因横跨松胚胎发育出现差异表达。一种锯齿状，花期早/SWR1复合体亚基6（小子/ SWC6)松树同系物(Pp30270)在整个胚胎发育过程中增加了其表达(聚类7)。小子编码了SWR1染色质重塑复合物的一个亚基，该亚基负责atp依赖的组蛋白H2A被H2A取代。Z变异，并通过开花抑制因子的正调控与拟南芥的开花抑制有关开花轨迹C（方法） 和疯狂的开花4（MAF4) ［64，65，66］．另外两个snf2相关的染色质重塑推定转录本遵循不同的表达谱。Pp14163，公认的松木同源物染色质重塑31（CHR31）在第5天和第15天达到峰，而同源物RAD5.(Pp26994)，在早期胚胎发生中表达过多，表明这些基因的表达依赖于发育阶段。在早期胚胎中，Pp34781也达到了高峰FASCIATA 1（FAS1）编码组蛋白蛋白染色质组件系数-1（CAF-1）的三个亚基之一的松树同源物。其拟南芥对应于后胚胎发育期间需要，用于适当的组织和功能的功能，但胚胎发育期间似乎不需要，至少在拟南芥晚期胚胎发生过程中[67］．

的转录p .松树胚胎发生中有大量的转录本，可能编码组蛋白亚基翻译后修饰的酶，特别是E2泛素结合酶。Pp30887是一种公认的E2泛素结合酶28基因的松树同源物，在胚胎早期和中期有差异表达和上调(聚类1)。已知这些蛋白进行H2B单泛素化，并与转录激活相关(综述于[68])。组蛋白赖氨酸甲基化也在松树胚胎发生中发挥了相关作用，根据对差异表达的转录本Pp44003和Pp33894的鉴定，推测分别编码拟南芥trithorax相关蛋白6 (ATXR6)(聚类3)和SU(VAR)3 - 9相关蛋白5 (SUVR5)(聚类10)。ATXR5和ATXR6参与抑制染色质修饰H3K27me1，有助于保持构成性异染色质状态，在大多数情况下是转座子和其他重复和沉默的元件，并防止在相同的细胞周期中发生重复复制[69，70］．另一方面，SUVR5负责抑制染色质修饰H3K9me2，独立于DNA甲基化的存在[71，其存在通常与转座子沉默和DNA甲基化控制有关，主要存在于中粒/着丝粒区[72］．

甲基化1/ORTHRUS 2 (VIM1/ORTH2)基因变体的推测同源物，一种与DNA甲基转移酶1 (MET1)合作促进CpG甲基化和着丝粒组织的甲基胞嘧啶结合蛋白[73，74]，在我们的数据中也发现(Pp11214)在胚胎早期(聚类3)中表达量最大，并在成熟胚胎中逐渐减少。在挪威云杉，vim1.在两种不同温度(18°C和30°C)条件下生长的早期形态发生期体细胞胚之间存在差异表达，在较高温度下表达上调[23］．然而，Pp3794是该植物的同源物畴重排甲基2（DRM2)，是一种主要的DNA甲基化转移酶基因，负责所有序列环境(CG、CHG和CHH)的DNA甲基化(由[75]），在第5天进行表达曲线峰，但从中间胚胎发生朝向成熟（簇5）增加。DRM2通过Novo甲基化维持ChH甲基化，通常在RNA定向的DNA甲基化（RDDM）途径内（在[76])。虽然在达到成熟前的松针推定DRM2峰的表达，但检测到从中间胚胎阶段稳步增加到成熟胚胎的趋势。该观察结果与拟南芥的研究一致表明成熟胚胎表现出CHH甲基化位点的饱和度，以及RDDM的更高活性和表达DRM2与早期胚胎相比[77，78］．

最后，在松树胚胎转录组中发现了许多与RNA沉默相关的转录本。在这一组中，值得一提的是Pp12441，它是RNA沉默player ARGONAUTE 2的公认松树同源物，与生物应激反应有关[79]，在胚胎发生的早期和中期上调(聚类1)。

这项工作提供了一个额外的资源来帮助理解基因调控，并与针叶树胚胎发育进程相关的重大事件。通过使用RNA-seq的技术接入的特定胚胎发育阶段中表达的基因，我们扩大已用DNA芯片杂交技术获得的沿海松树合子胚的先前公布的转录组分析。这样一个重要的结果是，碳水化合物的运输和代谢的发现清楚地过度表现在早期胚胎阶段。分配给eggNOGs功能类别无论是分析，或差异表达的转录本的富集分析和许多假定糖苷水解酶和糖转运基因，朝松胚胎发育及其相关的作用点的识别。另一个相关的结果提供给以前的研究强有力的支持是占主导地位，在早中期胚胎染色体的几个事件的，由相关联的转录增加，以及转座子的沉默组蛋白及其翻译后修饰的活跃合成重塑证明。

虽然没有公布的基因组p .松树该研究结果对致力于通过体细胞胚胎发生改善针叶树无性繁殖的植物研究领域具有重要意义。

植物材料

未成熟的雌球从收集授粉p .松树河中的小岛。生长在葡萄牙Mata Nacional do Escaroupim(经度8°44'W，纬度39°4'N)无性系种子果园的树木。该种子园是通过半同胞子代试验遗传选择的无性系顶端嫁接而建立的。加树最初是1963/64年在葡萄牙的Mata Nacional de Leiria(葡萄牙)选择的[80］．这些树是在葡萄牙政府的土地研究建立了一个实验种植园的一部分。从INIAV（Oeiras的，葡萄牙）获得的锥体，Ministério da Agricultura, Florestas e Desenvolvimento Rural，并在森林工程师Alexandre Aguiar和Isabel Carrasquinho的许可后提供，符合机构和国家指南。

六月中旬和7月结束之间发生的收集期。种子取出并用于分离胚胎如先前在[描述36］．通过[的分段系统后，快速评估每个胚胎的发育阶段30.]，胚柄被除去，并根据该阶段在液氮中胚胎立即冷冻成不同的池。五个不同的胚胎发育阶段被认为如下：Day0包括在早期胚胎阶段T0，T1和T2;第五天包括预子叶胚胎阶段T3和T4;Day11包括早期子叶胚胎阶段T4B;第15日包括在子叶胚阶段T5;和Day25包括在成熟胚阶段T7。根据胚胎阶段，每个池包含20-65合子胚。Several separate pools were prepared for each stage and samples were stored at − 80 °C until further analysis.

RNA提取和测序

根据制造商的说明，使用RNeasy Plant Mini试剂盒(Qiagen, Valencia CA, USA)从每个胚胎池中提取RNA，使用缓冲液RLC。使用ND-1000分光光度计(NanoDrop, Wilmington DE, USA)测定RNA的产量和纯度，0.8%琼脂糖凝胶电泳和RedSafe™核酸染色液(iNtRON Biotechnology)染色检查完整性。根据制造商的说明，使用RNase-Free DNase I (Qiagen)从DNA污染中清洗RNA样本。来自5个胚胎期池(一个生物复制)的总RNA样本被送到测序服务提供商，Illumina RNA-seq库在那里准备，并使用HiSeq 2000平台进行测序。

RNA-SEQ数据预处理和综合组件

全面的转录组合组装方法[31]在[32使用]生成参考转录组(除非另有说明，软件使用默认参数运行):使用Trimmomatic过滤原始读取(v 0.32;［81.)，使用默认选项来删除对端读取illumina适配器序列，低质量或低复杂度读取(SLIDINGWINDOW:4:5)，从两端5 bp (LEADING:5 TRAILING:5)，读取小于25 bp (MINLEN:25)。所有干净的reads用于与Trinity进行从头组装，但只有两个reads都保留的对用于引导组装或表达分析。三一(v 2.0.6;［33)用于生成带有默认参数的从头组装，加上“——min_glue 4——CuffFly——group_pairs_distance 600——genome_guided_max_intron 10000”。的p . taeda基因组和注释（v 1.01）从加利福尼亚大学的基因组项目下载[34］．来自每个胚胎阶段的干净的读取对独立地与此对齐p . taeda使用GSNAP无基因注释（V 2014年8月4日的基因组; [82.])。We only used “concordant paired” alignments where both reads in a pair align with a minimum length of 25 bp (50 bp per alignment), in the right forward-reverse orientation and insert length distance. The read alignments from all the stages, together with thep . taeda基因注释，用作袖扣中的输入（v.2.2.1; [83.)来重建另一组转录本。最后一步包括结合两组成绩单，从Trinity和Cufflinks:这些成绩单对齐p . taeda基因组与GMAP (v 2014-12-22， [84.]），并且所述对准转录物由基因座PASA（释放20140417群集; [85.])。然而，没有一致的文字记录p . taeda基因组通过使用“-pid1”的CD-ind序列聚集（V 4.6.3; [86.]），后来通过检查每个用TransDecoder（V 2.0.1; [87.]）使用所有的Viridiplantae作为参考，最后连接到之前的转录本(通过位点聚集的)。为了生成最终的综合转录组，我们筛选了602个没有ORF且小于200 bp的短转录本。602份成绩单列在附加文件中2．

功能注释

利用NCBI BLASTX将转录本与NCBI非冗余(nr)和拟南芥TAIR蛋白数据库进行比较，e值为1e-10。结果在Blast2GO中输入[88.注释氧化石墨烯术语、酶蛋白编码和KEGG通路。将转录本中的保守基序和结构与EBI InterPro (http://www.ebi.ac.uk/interpro/interproscan.html)．我们使用植物转录因子数据库(PLNTFDB，http://plntfdb.bio.uni-potsdam.de./)作为在我们的转录组中识别转录因子和其他转录调节因子的参考。该数据库包含了来自不同植物物种的近30,000条实验鉴定的元素的蛋白质序列，并通过HMM方法将它们按照蛋白质结构域进行分类。差异表达的转录本序列使用BLASTX与PLNTFDB进行比对，最小e值为1e-10。我们认为任何结果低于该阈值的转录本都是TF/转录调节因子，并将其在同源基因家族中以较低的e值进行注释。

如前所述，所有的ORF可能性都是使用TransDecoder从新组装的转录组中产生的，但只有每个转录本最长的ORF保留在最终p .松树蛋白质组。我们分类了任何组装的成绩单，其中Transdecoder无法识别ORF作为非编码。其他针叶树的转录组从戴维斯大学的Pinerefseq项目下载。他们每个人的蛋白质组成与类似的方式建造p .松树TransDecoder。通过将所有蛋白与BLASTP进行比对，确定最佳互反命中(BRH)， e值为1e-5。任何来自近亲属的注释都被分配回原注释p .松树成绩单。的蛋白质P. Pinaster，P. Taeda和p . lambertiana根据各自在EBI的功能注释蛋白eggNOG数据库中最佳同源序列的eggNOG组聚在一起。每个蛋白质首先与已经分类的蛋白质数据库进行比对，然后用其最佳同源物的eggNOG编号进行注释。将三种针叶树的蛋白以相同的egggnog编号连接在一起，形成蛋白质簇。蛋白质序列及描述注释数据库可于EMBL下载(http://eggnogdb.embl.de/)。

表达分析、基因富集及聚类分析

每个发育阶段的干净阅读与Bowtie的综合转录组(v 2.2.5;［89.]）和仅保留对在右的距离和方向（--very敏感-X 1000 --no混合--no-不和谐）读取映射。从这些比对，每个转录物的丰度在每个阶段使用快速（V量化1.5.1; [90.]）制作在Edger下游使用的计数表[91.］．少于10个计数的转录本此时被丢弃。由于在发育的早期阶段分离足够数量的合子胚胎的技术困难，我们没有复制的文库。为了将缺少复制的影响降到最低，我们遵循了开发人员在这种情况下推荐的协议[92.:用“x = .]”计算整个数据集的共同离散度dgelist（rowd.delim（'express.eft.eft_count'）），group = c（1,1,1,1,1,1）; x = calcnormfactors（x）; x = emplatedisp（x）;”，然后在FDR < 0.05时，使用edgeR’s exact test“exactTest()”确定每对连续阶段之间的差异表达转录本，该测试允许共同分散和单因素实验(我们的时间)。

通过Mayday中的K-means聚类，将差异表达的转录本按照每个阶段对齐读的归一化数量划分为聚类[93.]基于表达式值之间的欧几里德相关性。BLAST2GO中使用每个群集中的成绩单列表以识别丰富的GO条款。富集分析基于F-FISHER测试（FDR <0.05）。使用Resvik算法选项使用Revigo分配Go项之间的关系[94.]，并使用Treemap库在R内绘制(github.com/mtennekes/treemap.git)．按功能分类构建表达热图，将该分类中属于某一蛋白组的每个转录本的数量相加，转化为Z-scores，聚类，利用Mayday绘制热图[93.］．

通过RT-qPCR表达验证

来自De Novo组装的基因子集p .松树基于差异表达和/或表观遗传相关的注释选择转录组，以通过RT-QPCR验证RNA-SEQ结果。primer3plus（http://primer3plus.com/cgi-bin/dev/primer3plus.cgi)设计RT-qPCR引物，条件为:气相色谱含量为50 ~ 60%，Tm在50 ~ 65℃(以最近邻Tm为准);应该结合转录本的5 '或3 '较少保守区域;放大器大小为75-200 bp(附加文件18)．用PCR引物统计（http://www.bioinformatics.org/sms2/pcr_primer_stats.html.）和oligoanalyzer 3.1（https://eu.idtdna.com/calc/analyzer)．

首先使用RNA BR检测试剂盒(Thermo Fisher Scientific)在Qubit 3.0荧光仪中对RNA样本进行定量。根据定量RT-PCR标准程序中的制造商说明，使用转录子高保真cDNA合成试剂盒(Roche Diagnostics)进行cDNA合成，每20 uL反应混合添加1000 ng总RNA。每个胚胎期池准备3个生物重复。

所有QPCR实验都在用96孔白板（Roche诊断）（Roche诊断）（Roche诊断）中进行，其中每20μlQPCR反应混合物包括1x Sybr Green I Master（Roche诊断），每个引物和2μl的500nm1:20稀释cDNA。为每种生物复制制备三种技术复制。所有基因的扩增程序是相同的，在退火温度下具有轻微差异：95℃，10分钟，在95℃，20s时的40℃，20s在退火温度下和72℃下的8秒（退火温度在附加文件中18)．为了在放大时检查是否存在单峰，我们得到了熔化曲线。为了丢弃任何残余gDNA污染的可能性,控制同一rna样本准备但没有逆转录酶的酶,并用作模板与ELF1A qPCR扩增引物,哪里没有信号证明先前的有效性RNase-free DNase我(试剂盒)治疗。作为额外的对照，非模板对照(NTC)和阳性对照被包括在所有的试验中。此外，在所有的试验中，还包括从5个胚胎期池中分别提取等量总RNA的3个池样本的生物复制。每对引物的效率按[95.，使用ROCHE LC480软件检索的Cp值(附加文件18)．每个感兴趣基因(GOI)的相对表达被有效修正，如[95.]，使用ELF1a，atup和histo3作为参考基因，其在先前的研究中显示出可靠的参考基因p .松树体细胞胚胎发生[96.］．

将RNA-seq实验和RT-qPCR得到的数据进行比较。从RNA-seq中，一个基因在每个发育阶段的计数(来自eXpress)与同一基因在所有发育阶段的平均计数之间的碱基2的对数比率被计算出来。对RT-qPCR获得的数据采用了类似的方法，即对感兴趣基因在每个发育阶段的归一化数量(delta-delta-Ct)与同一基因在所有发育阶段的平均归一化数量进行碱基2的对数比分析。如果35个周期后RT-qPCR无法检测到表达，我们将表达值设为低于之前检测的值。

BRH:

最佳倒数命中

去：

基因本体论

“诱导多能性”:

InterPro数据库

门店:

新一代测序

ORF：

开放阅读框

RNA-SEQ：

RNA序列

来自INIAV的Isabel Carrasquinho和Alexandre Aguiar因提供植物材料而获得认可。Marta Simões因参与了植物材料的准备而受到认可。

资金

This work was supported through projects funded by (1) the European Commission Seventh Framework Programme (FP7, Grant Agreement N°289841-PROCOGEN), and (2) Fundação para a Ciência e a Tecnologia (FCT), through grants GREEN-it (UID/Multi/04551/2013), IF/01168/2013 and the doctoral fellowship SFRH/BD/79779/2011 (to ASR). JDV was supported by an Institute Programme Grant to Earlham Institute (BBS/E/T/000PR9818) from the Biotechnology and Biological Science Research Council (BBSRC), UK. The funding bodies had no role in the design of the study, collection, analysis, and interpretation of data, or in writing the manuscript.

可用性数据和材料

的数据沉积在欧洲核苷酸存档（ENA）的研究PRJEB21602下，与运行种质ERR2238559（阶段T0 / T1 / T2），ERR2238560（阶段T3 / T4），ERR2238561（阶段T4B），ERR2021873（阶段T5）和ERR2021874（阶段T7）[https://www.ebi.ac.uk/ena/data/view/prjeb21602.］．转录组鸟枪组装项目已在DDBJ / ENA / GenBank中加入GGEX00000000下沉积。这里所描述的版本是第一个版本，GGEX01000000 [https://www.ebi.ac.uk/ena/data/view/GGEX01000000］．

作者指出

1. Andreia S. Rodrigues和José J. De Vega对这项工作做出了同样的贡献。

从属关系

1. Instituto de Biologia实验ETecnológica（IBET），Apartado 12,2780-901，Oeiras，葡萄牙

   安德烈S. Rodrigues＆CéliaM.Miguel

2. Instituto de technologia Química e Biológica António Xavier, Universidade Nova de Lisboa (ITQB Nova)， Av. da República, 2780-157, Oeiras，葡萄牙

   安德烈S. Rodrigues＆CéliaM.Miguel

3. 英国诺维奇NR4 7UZ诺维奇研究园Earlham研究所

   JoséJ.Dega

4. Universidade de Lisboa, Faculdade de Ciências, BioISI - Biosystems & Integrative Sciences Institute, Campo Grande, 1749-016, Lisbon, Portugal

   CéliaM.Miguel.

贡献

本研究被CMM构思和引导。ASR进行了实验工作，包括制备用于测序和RT-QPCR的总RNA和基因表达验证。JDV进行了生物信息学分析。asr;JDV和CMM参与了结果分析及其生物解释。ASR，JDV和CMM写了这篇论文。所有作者阅读并认可的终稿。

通讯作者

对应到CéliaM.Miguel.．

伦理批准和同意参与

不适用。

同意出版物

不适用。

竞争利益

两位作者宣称他们没有相互竞争的利益。

出版商的注意

施普林格《自然》杂志对已出版的地图和机构附属机构的管辖权要求保持中立。

开放获取本文遵循知识共享署名4.0国际许可协议(http://creativecommons.org/licenses/by/4.0/)，它允许在任何媒体上无限制地使用、分发和复制，前提是你给予原作者和来源适当的荣誉，提供一个到知识共享许可协议的链接，并指出是否作出了更改。创作共用及公共领域专用豁免书(http://creativecommons.org/publicdomain/zero/1.0/)适用于本条提供的数据，除非另有说明。

重印和权限