多花黄精种子微根茎基因表达特征分析

刘保财 陈菁瑛 张武君 黄颖桢 赵云青 刘剑超 危智诚

（1. 福建省农业科学院农业生物资源研究所，福州 350003；2. 福建省农业科学院药用植物研究中心，福州 350003；3. 福建省南平市邵武市农业农村局，邵武 354000；4. 福建和平古镇农业开发有限公司，邵武 354000）

种子是植物贮藏营养和繁育后代的重要器官，也是人及家禽畜牧动物的主要食物来源［1]，种子萌发过程因胚发育程度不同而存在差异，具有成熟胚的种子，遇水吸胀后胚根首先突破种皮并向下生长形成根，胚芽向上生长露出地面形成茎，从而形成完整的植株［2-3]。然而部分植物种子不具有成熟胚，萌发过程存在特殊的形态结构，如黄精（Polygonatum sibiricum）［4]、多花黄精（Polygonatum cyrtonema）［5]、滇重楼（Paris polyphylla）［6]等百合科具有多年生根茎植物的种子，这些种子萌发过程中存在第一节微根茎形成的特殊萌发现象［7-9]。

多花黄精（P. cyrtonema）系百合科（Liliaceae）黄精属（Polygonatum）多年生草本［10]，《中华人民共和国药典》（2020）规定的中药饮品黄精的基原物种之一，具有补气养阴、健脾、润肺、益肾的功效，药食同源物种，甚至一些地方将其列为杂粮作物，近年逐步开始大规模的人工栽培［11-12]。据报道［8]及结合笔者研究，多花黄精种子萌发过程如下，胚根首先突破种皮（图1-A），然后生长一段时间（图1-B），在中胚轴的地方，逐步膨大（图1-C），形成最初的根茎（微根茎，图1-D），进一步在其上端分化形成芽生长点，生长第一片真叶，同时下端分化形成侧根/不定根，从而开始每年在根茎上形成芽和根的生长模式。在此萌发过程中，从胚根突破种皮到微根茎的形成，基本全部利用种子储藏的营养物质，微根茎形成后，逐步变为绿色，开始光合作用，原种子逐步变为空瘪状态，由“异养”变为“自养”（图1-E），进一步生长发育，形成小苗（图1-F）。

图1 多花黄精种子萌发过程Fig. 1 Germination process of P. cyrtonema seeds

当前，为解决多花黄精生产上种苗短缺问题，对多花黄精种子催芽、栽培等进行了研究并报道［8,11-12]，然而关于其萌发过程中调控微根茎发育的关键基因及参与的生理生化变化报道较少，制约了多花黄精种子生理、生产上微根茎的快速生长及其开发利用等研究，为此开展了多花黄精种子萌发过程中微根茎发育研究。本文以胚根突破种皮（图1-A）、微根茎（图1-D）和根状茎变绿（图1-E）3个萌发阶段为实验材料，基于二代高通量测序技术，进一步借助生信分析相关方法，对微根茎发育过程中的表达差异基因进行分析，以期揭示多花黄精种子微根茎发育过程中基因表达特征及其参与的通路变化，为多花黄精的种子生理、形态结构、育种与栽培等方面的研究奠定基础。

1 材料与方法

1.1 材料

试验材料采集于福建省泰宁县野生多花黄精，经福建省农业科学院农业生物资源研究所陈菁瑛研究员鉴定为Polygonatum cyrtonema Hua，种植于邵武市和平镇和平村林下。2019年11月采集10株果皮发黑的成熟果实，经8 d堆置发酵后，用手捏碎变软的果实，漂去果皮等杂质并用水清洗，干净的种子作为试验材料备用。

1.2 方法

1.2.1 材料处理 将洗干净的种子用无菌水冲洗10遍，移到超净工作台上，57%的酒精浸泡1 min，然后用无菌水冲洗5遍，放入经过121℃灭菌30 min的河沙内，再倒入少量的无菌水，保持沙层湿润，进行沙藏，每周定期检查并补充水分；150 d后，挑选胚刚突破种皮的种子（编号A，图1-A）和形成的微根状茎（编号D，图1-D），种植20 d后，挑选变绿的微根状茎（编号E，图1-E），这些种子和微根状茎用清水冲洗干净，吸取表面的水分，立即放入液氮中冻存10 min，然后移入-80℃冰箱中保存，直到RNA提取。

1.2.2 建库测序 采用DP441-RNAprep pure多糖多酚植物总RNA提取试剂盒（天根生化科技（北京）有限公司）提取总RNA，NanoPhotometer spectrophotometer检测RNA纯度（OD260/280及OD260/230比值）和Agilent 2100 bioanalyzer检测RNA完整性。进一步使用试剂盒Illumina 的 NEBNext®UltraTMRNA Library Prep Kit，参照说明书构建文库。文库构建完成后，用Qubit2.0 Fluorometer进行初步定量，用Agilent 2100 bioanalyzer对文库插入片段大小检测，用RT-qPCR对文库有效浓度进行准确定量。库检合格后，由北京诺禾致源科技股份有限公司应用Illumina 测序完成序列的测定。

1.2.3 基因表达水平分析 采用RSEM软件［13]，调用Bowtie2对比并对结果进行统计，以Trinity拼接得到的转录组作为参考序列（Ref），将每个样品的clean reads匹配到Ref，获得单个样品的匹配率，统计每个样品中基因比对到的reads数目，并进行FPKM（expected number of fragments per kilobase of transcript sequence per millions base pairs sequenced）转换，获得单个样品中基因和转录本的表达水平。

1.2.4 差异基因表达分析 使用基于负二项分布的模型确定数字基因表达数据中的差异表达地DESeq2 R包［14]（1.10.1），进行两个条件/组的差异表达分析，使用Benjamini和Hochberg的方法校正得到的P值（Padj）以控制错误发现率，Padj＜0.05的基因认定为差异表达。基于 GOseq［15]所述方法对筛选到的DEGs进行GO富集分析。基于KEGG注释结果，使用 KOBAS，Padj＜0.05的基因认定为差异表达，进行 Pathway富集分析。

1.2.5 RT-qPCR 验证 随机选取通路上的7个差异基因，设计引物和合成引物序列见表1，以种子萌发的3个不同阶段的RNA 为模板，按照 TransScript®All-in-One First-Strand cDNA Synthesis SuperMix for qPCR（One-Step gDNA Removal）（AT341）试剂盒进行反转录合成cDNA。参照 PerfectStart® Green qPCR SuperMix（+Dye II）试剂盒说明书，以cDNA为模板，在ABI QuantStudio 3仪上完成 RT-qPCR。用2-ΔΔCt方法［16]计算每个基因在不同样品中的相对表达量，Actin为内参基因。

表1 RT-qPCR引物序列Table 1 Gene-specific primers used in RT-qPCR

2 结果

2.1 微根茎发育过程基因表达水平

以组装后的序列为模板，将胚根突破种皮（A1、A2、A3）、微根状茎形成（D1、D 2、D3）和微根状茎变为绿色（E1、E2、E3）的clean reads与之比对，各样本的匹配reads数目和匹配率见表2，胚根突破种皮时匹配率为71.63%-74.23%，微根状茎形成时的匹配率为73.60%-75.09%，绿色微根状茎的匹配率为74.81%-76.03%。对各样品中表达的Cluster数目进行统计并转换为FPKM，再根据FPKM大小划分区间，各区间的表达量及所占比例见表3，各样品的FPKM＞15的Cluster数量没有较大的差异，即各样品中高表达的基因数目差异较少。

表2 各样品reads 与组装转录本比对Table 2 Mapped results of sample reads and assembly transcripts

表3 样品表达水平FPKM区间数量统计Table 3 FPKM interval statistics of sample expressed levels

2.2 微根茎形成与胚根突破种皮的差异基因分析

将微根茎形成时与胚根突破种皮比较（图2-A），共有显著性差异的Unigenes 17907条，其中上调的有9797条，下调的为8110条。进一步将上调或下调差异基因在GO数据库中富集分析，选取最显著的40个Terms绘制柱状图（图2-B），微根茎形成时与胚根突破种皮相比，所有的Terms中上调的Unigenes数量大于下调的Unigenes数量。就差异Unigenes数量而言，无论是上调还是下调，主要集中富集在生物过程（biological process, BP）和分子功能（molecular function, MF），而细胞组成（cellular component, CC）的Unigenes数量相对较少。在BP中参与代谢过程（metabolic process）的差异基因有4952条，其中上调的有3135条（如Cluster-68615.36386、Cluster-68615.82941等），下调的有1817条（如Cluster-68615.49491、Cluster-68615.48139等）；参与蛋白质磷酸化（protein phosphorylation）、磷酸化（phosphorylation）、含磷化合物代谢过程（phosphate-containing compound metabolic process）、磷代谢过程（phosphorus metabolic process）等与磷相关的Terms均有较高的差异表达。在MF中主要参与的有催化活性（catalytic activity）、转移酶活性（transferase activity）、水解酶活性（hydrolase activity）、离子结合（ion binding）等Terms，其中参与催化活性的差异基因有4196条，上调的有2814条（如Cluster-68615.86087、Cluster-68615.91384等），下调的有1382条（如Cluster-68615.14912、Cluster-68615.64310等）。CC与MF和BP相比，参与的差异基因较少，显著差异的前40个Terms中，Terms数量也较少。

图2 DvsA差异表达基因Fig. 2 DvsA different expressed genes

通过Pathway显著性富集确定差异基因间相互协调及其生物学功能，将差异表达基因参与的最主要生化代谢途径和信号转导途径富集后，挑选了富集最显著的20条pathway条目进行展示（图3），差异基因主要富集于植物激素信号转导（plant hormone signal transduction）、淀粉和蔗糖代谢（starch and sucrose metabolism）、光合作用（photosynthesis）、光合作用-天线蛋白（photosynthesis - antenna proteins）、苯丙烷生物合成（phenylpropanoid biosynthesis）、黄酮类生物合成（flavonoid biosynthesis）等通路中，而这些通路与种子萌发过程中胚形态结构发育和代谢有直接的关联。富集到差异最显著苯丙烷生物合成通路中基因有240条，有显著差异的有108条（如Cluster-32281.0、Cluster-68615.50869等）。富集到植物激素信号转导通路中的差异基因有354条，有显著差异的有133条（如Cluster-68615.83466、Cluster-68615.49580等）。富集到淀粉和蔗糖的代谢通路的基因有397条，有显著差异的有132条（如Cluster-68615.91668、Cluster-68615.7364等）。

图3 DvsA差异基因KEGG通路富集Fig. 3 KEGG pathway enrichment of DvsA differential gene

为了进一步解析差异基因中上调和下调基因参与的生理与代谢通路，分别将上调的差异基因的前20条通路（图4-A）和下调的前20条通路进行了KEGG通路富集（图4-B）。在上调的差异表达基因中，主要富集在了芪类、二芳基庚烷和姜醇生物合成（stilbenoid, diarylheptanoid and gingerol biosynthesis）、淀粉和蔗糖代谢（starch and sucrose metabolism）、植物激素信号转导（plant hormone signal transduction）、光合作用-天线蛋白（photosynthesis-antenna proteins）、苯丙烷生物合成（phenylpropanoid biosynthesis）等通路中，其中富集到淀粉和蔗糖代谢、植物激素信号转导等通路的差异基因显著上升。下降的差异基因主要富集在了缬氨酸、亮氨酸和异亮氨酸的生物合成（valine, leucine and isoleucine biosynthesis）、辅酶Q及其他萜-醌合成（ubiquinone and other terpenoid-quinone biosynthesis）、植物激素信号转导（plant hormone signal transduction）、剪接体（splicesome）、RNA转运（RNA transport）等通路中，其中植物激素信号转导、剪接体、RNA转运等通路中富集的差异基因数量明显下降。而无论是表达上调还是表达下调的差异基因中，富集到植物激素信号转导通路中的基因均有大量表达，即植物激素在微根茎的形成中发挥了重要作用。

图4 TOP 20 差异基因KEGG富集结果Fig. 4 KEGG enrichment results of TOP 20 differential genes

基于上述分析，淀粉和蔗糖的代谢与植物激素代谢在根茎初步形成过程中具有关键作用，因此进一步将DvsA的差异基因映射到淀粉和蔗糖代谢（ko00500）通路，结果（图5）表明，参与编码磷酸蔗糖合酶（sucrose-phosphate synthase, 2.4.1.14）、己糖激酶（hexokinase, 2.7.1.1）和α-海藻糖酶（alphatrehalase, 3.2.1.28）的差异转录本均全部下调，参与编码α-淀粉酶（alpha-amylase, 3.2.1.1）、糖原磷酸化酶（glycogen phosphorylase, 2.4.1.1）、4-α-葡聚糖转移酶（4-alpha-glucanotransferase, 2.4.1.25）、麦芽糖酶-葡糖淀粉酶（3.2.1.20maltase-glucoamylase）、β-D-呋喃果糖苷水解酶（beta-fructofuranosidase, 3.2.1.26）、葡萄糖-1-磷酸腺苷酰转移酶（glucose-1-phosphate adenylyltransferase, 2.7.7.27）等的差异转录本均全部上调，而参与编码果糖激酶（fructokinase, 2.7.1.4）、β-葡萄糖苷酶（beta-glucosidase, 3.2.1.21）等的差异转录本则既有上调也有下调。

图5 参与淀粉和蔗糖代谢的DvsA差异表达基因Fig. 5 DvsA differentially expressed genes involved in starch and sucrose metabolism

为了进一步剖析植物激素相关基因对微根茎初步形成的影响，将DvsA的差异基因映射到植物激素信号转导（ko04075）通路上，结果（图6）表明，参与编码生长素（包括AUX/IAA、ARF、GH3、SAUR）、细胞分裂素（B-ARR、A-ARR）和赤霉素（TF）各关键酶的基因均有上调和下调，然而这些基因中，除 ARF和GH3外，所有表达上调的基因数多于表达下调的基因数，如AUX/IAA中，仅有1条基因表达下调，11条表达上调。参与编码油菜素内酯各关键酶基因表达均上调，参与编码脱落酸关键酶（PYR/PYLH、ABF）基因表达全部下调，编码关键酶（PP2C、SnPK2）基因的下调转录本多于上调转录本，再次说明参与编码植物激素的基因在微根茎形成中具有重要作用，而油菜素内酯可能是根茎初步形成的一类关键激素。

图6 参与植物激素信号转导的DvsA差异表达基因Fig. 6 DvsA differentially expressed genes involved in plant hormone signal transduction

2.3 微根茎变绿前后差异基因分析

比较微根茎变绿前后差异基因表达（图7-A），共有显著性差异的Unigenes 26833条，其中上调的有14567条，下调的为12266条。进一步将上调或下调差异基因在GO数据库中富集，选取最显著的40个Terms绘制柱状图（图7-B），微根茎变绿前后相比，除了与过氧化酶反应有关的3个Terms外，所有的同一Terms中上调的Unigenes数量大于下调的Unigenes数量。就差异Unigenes数量而言，无论是上调还是下调，主要富集在生物过程（bological process, BP），而分子功能（molecular function, MF）和细胞组成（cellular component, CC）的Unigenes数量相对较少。在BP中参与代谢过程（metabolic process）的差异基因有6896条，其中上调的有3758条（如Cluster-68615.21356、Cluster-68615.32031等），下调的有3138条（如Cluster-68615.92980、Cluster-68615.3727等）；参与肽生物合成过程（peptide biosynthetic process）、酰胺生物合成过程（amide biosynthetic process）、有机氮化合物生物合成过程（organonitrogen compound biosynthetic process）等生物合成通路中的基因均得到强烈富集，即这些基因参与了微根茎变绿前后的形态结构或生理活性。而在MF中主要参与的有催化活性（catalytic activity）、氧化还原酶活性（oxidoreductase activity）等Terms，其中参与催化活性的差异基因有5518条，上调的有2893条（如Cluster-15152.0、Cluster-68615.16027等），下调的有2625条（如Cluster-68615.69338、Cluster-68615.31773等）。在CC过程中，差异基因主要参与细胞质（cytoplasm）、大分子复合物差异（macromolecular complex）等Terms中。

图7 EvsD差异表达基因Fig. 7 EvsD differently expressed genes

通过Pathway显著性富集确定差异基因间相互协调及其生物学功能，将差异表达基因参与的最主要生化代谢途径和信号转导途径富集后，挑选了富集最显著的20条pathway条目进行展示（图8），差异基因主要富集于核糖体（ribosome）、淀粉和蔗糖代谢（starch and sucrose metabolism）、光合作用（photosynthesis）、光合作用-天线蛋白（photosynthesis-antenna proteins）、吞噬体（phagosome）、卟啉和叶绿素代谢（porphyrin and chlorophyll metabolism）等通路中，而这些通路直接体现了种子由“异养”向“自养”的转变，如富集到差异最显著核糖体通路中基因的有1403条，有显著差异的有640条；富集到光合作用、光合作用-天线蛋白的差异基因分别有97条和41条，有显著差异的有63条和33条。

图8 TOP 20差异基因KEGG富集结果Fig. 8 KEGG enrichment results of TOP 20 differential genes

为了进一步解析差异基因中上调和下调基因参与的生理与代谢通路，分别将上调的差异基因的前20条通路和下调的前20条通路进行了KEGG通路富集（图9-A和9-B）。整体而言，上调的差异基因数多于下调的差异基因数。在上调的差异基因中，主要富集在了光合作用（photosynthesis）、核糖体（ribosome）、光合作用-天线蛋白（photosynthesisantenna proteins）等通路中，其中注释到光合作用的有61条、光合作用-天线蛋白的有33条、光合生物中的碳固定的60条。差异基因中表达量降低的基因主要富集在了苯丙烷生物合成（phenylpropanoid biosynthesis）、淀粉和蔗糖代谢（starch and sucrose metabolism）、氨基糖和核苷酸糖代谢（amino sugar and nucleotide sugar metabolism）、吞噬体（phagosome）等通路中，其中富集于淀粉和蔗糖代谢通路中的差异基因有114条。值得注意的是植物-病原体相互作用通路中富集了74条。基于通路分析表明，微根茎变绿前，仍然是以淀粉、蔗糖代谢为主，变绿后，光合相关通路的基因进行了表达与上调，同时与植物病原互作的基因开始高表达。

微根茎变绿前后，光合通路中相关基因具有关键作用，因此进一步将EvsD的差异基因比对到光合作用（ko00195）通路上（图10），参与编码光合系统I（PsaA、PsaB、PsaD、PsaE、PsaF、PsaG、PsaH、PsaK、PsaL、PsaN、PsaO）、光合系统Ⅱ（PsbA、PsbB、PsbC、PsbK、PsbO、PsbP、PsbQ、PsbR、PsbS、PsbY、PsbZ、Psb27、Psb28）等关键酶的基因表达均为上调，仅有少量的如编码三磷酸腺苷酶ε基因表达下调。

2.4 差异基因RT-qPCR分析

为了验证转录组获得的数据，随机选择具有显著差异的7个基因，采用RT-qPCR分析（图11），各基因在各样本的表达量与转录组表达的趋势基本一致。

图11 7个差异基因的RT-qPCR验证Fig. 11 Verification of seven selected DEGs by RT-qPCR

3 讨论

3.1 微根茎发育前后的基因表达特征：淀粉和蔗糖的代谢与植物激素信号转导

多花黄精种子从胚根突破种皮到微根茎形成，发生一系列表观形态特征变化，陈怡等［8]对多花黄精的萌发过程进行了形态解剖观察，明确了成熟的多花黄精种子，其胚未分化出明显的子叶、胚芽和胚根，随着种子的萌发伴随着胚乳的降解而分化形成成熟胚。从胚根突破种皮到微根茎形成的差异基因经GO数据库中富集分析，主要参与代谢过程、催化活性、单生物代谢过程、氧化还原酶活性、核糖体等，KEGG显著性富集表明差异基因主要富集于核糖体、植物激素信号转导、淀粉和蔗糖代谢等通路中，上调的差异基因主要富集在了淀粉和蔗糖代谢、植物激素信号转导等通路中，下调的差异基因主要富集在了植物激素信号转导、剪接体、RNA转运等通路中，与前人报道多花黄精种子含有淀粉、粗脂肪、可溶性糖、游离氨基酸和植物甾醇含量较高，具有较高的营养价值和抗氧化酶活性结果一致［17-19]。进一步淀粉和蔗糖的代谢通路比对分析，参与编码磷酸蔗糖合酶、己糖激酶和α-海藻糖酶的差异转录本均全部下调，参与编码α-淀粉酶、糖原磷酸化酶、4-α-葡聚糖转移酶、麦芽糖酶-葡糖淀粉酶、β-D-呋喃果糖苷水解酶、葡萄糖-1-磷酸腺苷酰转移酶等的差异转录本均全部上调，而参与编码果糖激酶、β-葡萄糖苷酶等的差异转录本则既有上调也有下调，淀粉水解成葡萄糖，可能为微根茎形成提供必要的能量［19]。激素种类和含量的变化对多花黄精种子萌发具有显著影响，多花黄精种子萌发过程中出现IAA、GA3和TZR含量呈上升-下降-上升趋势，ABA含量降低并维持在较低水平［5,20]，这与本文所述的编码生长素、细胞分裂素和赤霉素各关键酶的基因均有上调和下调的结果基本一致。而参与编码脱落酸关键酶（PYR/PYLH、ABF）基因表达全部下调，编码油菜素内酯各关键酶均为上调，可见油菜素内酯可能是根茎初步形成的一类关键激素，与报道的油菜素内酯调控马铃薯块茎的形成结果一致［21]，其能否作为外源激素用以促进微根茎的膨大，有待进一步研究与分析。

3.2 外界环境通过内在成分影响种子萌发与微根茎的发育

低温可以有效打破和缩短黄精种子休眠，促进其发育不全的胚成熟、突破种皮和发芽，主要原因也是低温诱导了激素相关基因的表达，与本文研究的参与植物激素信号转导的基因可能调控多花黄精微根茎发育结果一致［22-24]。层积处理是打破多花黄精种子休眠的有效方法之一，种子沙藏后糖分种类、含量等均发生了变化，提高了种子的发芽率和发芽势［20,25]。外源激素的处理也有利于打破多花黄精种子的休眠，6-苄氨基嘌呤可打破上胚轴休眠，从而促进多花黄精种子的发芽［26]。有些微生物的分泌物也具有促进多花黄精种子萌发的效果［27]，这些分泌物是否促进微根茎发育，有待进一步研究。上述可知，多花黄精种子胚根突破种皮到第一节根茎的形成，涉及许多生理变化与形态结构发育，因此受到内部和外界多种因素综合影响。

3.3 微根茎变绿后即可进行光合作用

微根茎变绿前，以种子内贮藏的营养物质消耗为主，变绿后，种子内贮藏的营养物质基本消耗完毕，开始逐步依靠自身光合作用进行生长发育。因此微根茎变绿前后显著性差异基因GO富集分析，主要参与了代谢过程、肽生物合成过程、催化活性等；KEGG显著性富集表明，差异基因主要富集于核糖体、淀粉和蔗糖代谢、光合作用等通路中，这与前人报道的多花黄精种子萌发过程中差异基因主要参与碳代谢、次生代谢产物的生物合成和多糖的代谢结果一致［20,28]。微根茎变绿前后的差异基因在光合通路上比对后，除编码三磷酸腺苷酶ε等少量基因表达下调外，均为上调，表明微根茎已逐步开始光合作用，与光照影响多花黄精种子的萌发一致［22]，但光合作用的强弱尚待进一步研究。

3.4 微根茎营养价值亟需评价

近年对多花黄精根状茎的多糖［29-30]、甾体皂苷生物合成［31-32]等关键基因进行了研究，并解析了根茎多糖的药理［33-34]，然而多花黄精微根茎发育过程中所产生的多糖是否具有类似的功效，有待进一步研究。所产生的微根茎能否作为类似“豆芽”的芽菜类开发，也亟需进一步对其营养、药理、毒理等开展评价。

综上，通过本文阐述，明确了微根茎基因表达特征，为多花黄精种子萌发的生理、生产上促进微根茎快速膨大、微根茎的开发利用等相关研究奠定了基础。

4 结论

本研究获得了多花黄精种子萌发过程中不同萌发阶段的转录组数据，参与淀粉和蔗糖代谢和植物激素信号转导通路的基因在微根茎形成中发挥了重要作用，尤其是编码油菜素内酯通路关键酶的基因。微根茎变绿后，光合相关基因得到大量表达，已具有光合作用能力。