薯蓣茎生珠芽发生相关的转录组和PEBP家族基因表达特性分析

孙 玲， 王 月， 陈 闽， 孙小芹， 杭悦宇， 张艳梅

〔江苏省中国科学院植物研究所(南京中山植物园)， 江苏 南京 210014〕

薯蓣属(DioscoreaLinn.)为单子叶攀缘植物，其叶腋处常形成一种特殊的侧枝变态器官——茎生珠芽，也称为零余子，是薯蓣属植物重要的无性繁殖器官。薯蓣(D.polystachyaTurcz.)又称山药，是传统的中药和补益食品[1-2]，目前主要采用块茎繁殖。在实际生产中，由于薯蓣块茎繁殖过程中往往存在严重的病虫害感染，会导致品种退化，因而，其茎生珠芽常被用于薯蓣品种的复壮[3-6]。茎生珠芽与卷须、枝刺等其他侧枝变态类似，均由侧生分生组织发育而来，因此，探讨薯蓣茎生珠芽发生的遗传机制，对薯蓣品种繁育及植物侧枝发育机制探讨均具有重要意义。

在自然界中，除薯蓣外，具珠芽的植物还有龙舌兰(AgaveamericanaLinn.)[7]、台闽苣苔〔Titanotrichumoldhamii(Hemsl.) Soler.〕[8]、半夏〔Pinelliaternata(Thunb.) Breit.〕[9]、淡黄花百合(LiliumsulphureumBaker apud Hook. f.)[10]和白花虎眼万年青(OrnithogalumarabicumLinn.)[11]等。目前关于植物珠芽发生机制的研究主要集中在3个方面：一是光照和温度等环境因子的影响效应[12-13]；二是生长素和细胞分裂素等激素的影响效应[14-15]；三是遗传基因的调控机制[7-8，11]。例如：在龙舌兰花序珠芽的不同发育期，KNOXⅠ基因在花序珠芽起始处，特别是在分生组织发育的特定部位中有表达[7]；台闽苣苔的GFLO基因在其花序和花序珠芽原基中均有特异性表达，且表达量在花序珠芽发生过程中明显下降[8]；姜福星等[11]筛选出白花虎眼万年青叶面珠芽形成过程中差异表达的关键基因EMB和FKBP等。上述这些种类的花序珠芽或叶面珠芽的发生机制与茎生珠芽存在一定的差异。

He等[16]利用转录组测序技术，发现在卷丹(LiliumtigrinumKer Gawler)茎生珠芽的发育过程中，细胞分裂素合成基因IPT1/5和细胞分裂素激活基因LOG1/3/5/7的表达量增加，而细胞分裂素降解基因CKX4的表达量减少，导致其内源性iPA型细胞分裂素的含量增加，并促进了细胞分裂素信号通路基因(AHK2/3/4、AHP1和ARR1/2/12)的表达。Wu等[17]通过转录组测序，对参薯(DioscoreaalataLinn.)不同发育阶段的茎生珠芽进行了研究，认为生长素、细胞分裂素和脱落酸合成相关基因(WRKY、bHLH和MYB等)以及转录调节因子等均可能参与了参薯茎生珠芽的发育调控。但表达量发生变化，只能说明上述基因或转录调节因子可能参与了珠芽发育过程的调控，并不能直接佐证他们与珠芽的发生有关，因而，目前关于珠芽发生的分子机制仍不明确。

磷脂酰乙醇胺结合蛋白(PEBP)家族的TFL1-like基因，主要表达于次级分生组织，包括腋生分生组织、维管形成层和中间分生组织，可抑制茎顶端分生组织形成花原基，目前已知其具有延迟开花、维持茎尖分生组织营养生长、促进侧枝产生等功能[18]。例如：过表达RCN基因能导致水稻(OryzasativaLinn.)开花延迟，并产生更多分蘖和多穗的圆锥花序[19]；过表达CmTFL1c基因能导致杭白菊(Chrysanthemum×morifoliumRamat．)的侧枝明显增多、开花延迟，且花形态发生异常[20]。qRT-PCR研究结果[21]表明：在半夏的茎生珠芽不同发育期，PteTFL1基因的表达量存在显著性差异，且随其茎生珠芽的发育，PteTFL1基因表达水平先升高后降低，在茎生珠芽刚发生时的表达量最大，暗示PteTFL1基因在半夏茎生珠芽的生长发育中可能发挥重要作用。据此，在薯蓣茎生珠芽的发生过程中TFL1-like基因是否也具有同样的作用？有待深入探讨。

薯蓣品种‘铁棍山药’(Dioscoreapolystachya‘Tiegun’)和‘花籽山药’(D.polystachya‘Huazi’)的遗传背景极为一致，但前者茎中上部叶腋通常具茎生珠芽，基部叶腋则不具茎生珠芽，属于具茎生珠芽品种；后者茎上不具有珠芽，属于无茎生珠芽品种。作者以上述2个薯蓣品种为实验材料，通过转录组测序，挖掘可能参与调控薯蓣茎生珠芽发生的相关基因，并探究与侧枝发生密切相关的TFL1-like基因及其在薯蓣茎生珠芽发生过程中的作用，以期为薯蓣品种繁育及植物侧枝发育机制研究奠定基础。

1 材料和方法

1.1 材料

于2015年6月在徐州市农业科学院实验苗圃选取长势良好的‘铁棍山药’和‘花籽山药’植株，在‘铁棍山药’的茎上采集具茎生珠芽的叶腋部位(S1)、茎基部无茎生珠芽的叶腋部位(S2)和茎上部无茎生珠芽的叶腋部位(S3)3组样品，在‘花籽山药’的茎上采集叶腋部位(S4)1组样品，用于转录组测序。

在‘铁棍山药’茎上采集不同发育期茎生珠芽的叶腋部位，包括茎生珠芽的未发育期(叶腋部位未见珠芽，T0期)、发育初期(叶腋部位冒出白色小珠芽，T1期)、发育中期(叶腋部位珠芽由白色刚变为棕色，T2期)和发育成熟期(叶腋部位珠芽全变为棕色，T3期)的4组样品，用于DpTFL1基因表达量分析。

取样时，使用刀片在叶腋部位割取面积约4 mm2的表皮，包含珠芽但不包含维管组织。样品经液氮速冻后于-70 ℃保存、备用。每个部位或每个时期各采集约10株样株，混合后作为1份样品。

1.2 方法

1.2.1 总RNA的提取及纯化 将300 mg样品置于液氮中研磨成粉末，用RNA isolater Total RNA Extraction Reagent试剂(南京诺唯赞生物科技股份有限公司)提取总RNA；获得的RNA用不含RNase的DNase Ⅰ(美国Invitrogen公司)于37 ℃下消化处理30 min，去除DNA污染，然后用Agilent 2100生物分析仪(美国Agilent公司)检测RNA样品的完整性和纯度；最后使用Dynabeads mRNA Purification Kit试剂盒(美国Invitrogen公司)对mRNA进行分离纯化。

1.2.2 cDNA文库构建 用Fragmentation Buffer试剂(美国Invitrogen公司)将mRNA随机切割成150～200 nt的片段；以片段化的mRNA为模板，用6碱基随机引物和逆转录酶将mRNA片段反转录成cDNA的第1链，并在RNase H和DNA聚合酶Ⅰ的作用下合成双链cDNA；对双链cDNA片段进行末端修复及3′末端加‘A’碱基，并在片段两端连接上特定的测序接头；经质量体积分数1.2%琼脂糖凝胶电泳(120 V，20 min)后回收长度约200 bp的cDNA片段，通过PCR扩增富集测序样本。扩增体系总体积25.0 μL，包括20.0 μL金牌Mix (green)(北京擎科新业生物科技有限公司)、10 μmol·L-1正向和反向检测引物各2.0 μL以及0.05 ng·μL-1cDNA 1.0 μL。扩增程序：98 ℃预变性2 min；98 ℃变性10 s、58 ℃退火10 s、72 ℃延伸10 s，共30个循环；72 ℃延伸1 min。并用QIAquick Gel Purification Kit试剂盒(德国Qiagen公司)对PCR产物进行纯化。

1.2.3 Illumina测序和RNA-seq分析 构建完成的cDNA文库用Agilent 2100生物分析仪和StepOnePlus Real-Time PCR System(美国ABI公司)质检合格后，使用Illumina HiSeqTM2000测序平台(美国Illumina公司)进行测序，并将所得的图像数据转化为相应的核苷酸序列数据，获得原始测序数据。

在原始测序数据中去除接头序列和低质量reads(Q<20)，获得clean reads；使用Trinity分析软件(https：∥trinityrnaseq.sourceforge.net/)对获得的高质量序列进行转录组denovo拼接；使用TGICL软件对拼接的转录组进行聚类去冗余得到unigenes序列，并对各样品的unigenes再次进行聚类去冗余得到最终的unigenes，命名为“All-unigenes”，用于后续分析；对组装的contigs和unigenes进行长度分布统计。使用DESeq2软件对各样品进行差异表达水平分析，采用FPKM(每千个碱基的转录每百万映射读取的碎片)方法[22]进行定量分析，其中差异unigenes的筛选使用表达差异倍数(fold change，FC)的方法，判断条件为︱log2FC︱≥1(P≤0.05)[23]。

使用Blast2GO软件(https：∥www.blast2go.com)将所得的unigenes序列与NR数据库(https：∥www.ncbi.nlm.nih.gov/refseq/)进行BLASTx批量比对(E≤10-5)，获取最佳功能注释；根据NR注释信息进行GO功能注释，在获得每个unigene的GO注释后，使用WEGO软件(http：∥wego.genomics.org.cn/cgi-bin/wego/index.pl)对所有unigenes进行GO功能分类统计；将BLASTx所得unigenes与COG数据库(https：∥www.ncbi.nlm.nih.gov/COG)进行比对分析(E<10-5)，获得基因的COG功能注释和COG功能分类；通过KEGG数据库(https：∥www.genome.jp/kegg)的注释信息得到unigenes的Pathway注释；按NR、NT、GO、COG、KEGG、SwissProt和InterPro的顺序将unigenes序列与以上数据库进行BLASTx比对，翻译后得到该unigene编码区的核酸序列(序列方向为5′→3′)和氨基酸序列并预测编码蛋白框(CDS)，没有得到比对结果的unigenes则使用ESTScan软件预测其编码区。

1.2.4 PEBP家族基因的鉴定和表达量分析 从Phytozome数据库(https：∥phytozome.jgi.doe.gov/pz/portal.html)下载拟南芥〔Arabidopsisthaliana(Linn.) Heynh.〕和水稻基因组序列和注释信息，从Ensemble 数据库(ftp：∥ftp.ensemblgenomes.org/pub/plants/release-49/fasta/dioscorea_rotundata/dna/)下载几内亚薯蓣(DioscorearotundataPoir.)基因组序列和注释信息。以PEBP家族基因中保守的乙醇胺结合域(PF01161)为问询序列，采用BLAST和HMM(隐马尔科夫模型)方法，对获得的‘铁棍山药’和‘花籽山药’转录组数据库以及拟南芥、水稻和几内亚薯蓣基因组蛋白序列进行同源性比对(E值设为10-4)，使用MEGA 7.0软件中的ClustalW对同源序列进行比对，并使用PhyML软件对同源序列进行系统演化树构建[24-25]。以FPKM值表示转录本丰度，对具茎生珠芽和无茎生珠芽样本的转录组中全部PEBP家族基因的表达量进行分析和比较。

1.2.5 qRT-PCR分析 按前述方法提取‘铁棍山药’不同发育期茎生珠芽叶腋样品的RNA，并对DpTFL1基因表达进行qRT-PCR分析。使用Primer 6.0软件设计引物，DpTFL1基因正向引物qDpTFL1-F和反向引物qDpTFL1-R的序列分别为5′-GAGGGAGGTGGTGGAGTATGAA-3′和5′-ATTGAAGAAGACAGCGGCG-3′；参照文献[17]，根据薯蓣Actin基因设计引物，正向引物DpActin-F和反向引物DpActin-R的序列分别为5′-TGACGAGGATATTCAACCCCT-3′和5′-GATACCCCTCTTGGATTGAGC-3′。

采用ABI StepOnePlus荧光定量PCR仪(美国ABI公司)进行qRT-PCR分析。扩增体系总体积20.0 μL，包括SYBR Green Master Mix(南京诺唯赞生物科技股份有限公司)10.0 μL、10 μmol·L-1正向和反向引物各0.8 μL以及0.05 ng·μL-1cDNA 2.0 μL，用ddH2O补足至20.0 μL。扩增程序：95 ℃预变性5 min；95 ℃变性10 s、60 ℃退火30 s，共40个循环。扩增反应结束后进行融解曲线分析，分析程序：95 ℃ 15 s、60 ℃ 1 min、95 ℃ 15 s。采用2-ΔΔCT法[26]进行基因相对表达量分析。

1.3 数据处理和统计分析

使用SPSS 20.0软件对‘铁棍山药’不同发育期叶腋样品DpTFL1基因的相对表达量进行显著性检验(t-检验，α=0.05)。

2 结果和分析

2.1 转录组测序和de novo组装结果

‘铁棍山药’和‘花籽山药’不同部位样品中转录组unigenes的质量指标见表1，2个品种叶腋部位转录组unigenes的长度分布图见图1。

表1 薯蓣品种‘铁棍山药’和‘花籽山药’不同部位样品中转录组unigenes的质量指标

图1 薯蓣品种‘铁棍山药’和‘花籽山药’叶腋部位转录组unigenes的长度分布图

结果显示：共获得77 981个All-unigenes，长度主要集中于300～3 000 bp，占总量的99.44%。All-unigenes的平均长度为1 050 bp，N50值为1 764 bp，GC含量为42.66%。4组样品的unigenes总数、总长度、平均长度、N50值、N90值和GC含量均有不同程度的差异。

2.2 功能注释分析

2.2.1 功能注释结果 将获得的unigenes序列比对到7大功能数据库，结果显示：NR数据库注释了50 030个unigenes，NT数据库注释了44 565个unigenes，GO数据库注释了5 392个unigenes，COG数据库注释了22 961个unigenes，KEGG数据库注释了38 681个unigenes，SwissProt数据库注释了35 872个unigenes，InterPro数据库注释了38 807个unigenes。其中，NR数据库提供了比其他数据库更多的unigenes注释，大约占All-unigenes数的64.2%。

将‘铁棍山药’和‘花籽山药’不同部位样品的转录组unigenes与NR数据库进行同源序列比对，结果显示：供试2个品种的转录组unigenes与4个种类具有较高的同源性，其中，与油棕(ElaeisguineensisJacq.)的同源性最高，有17 580个unigenes相似，占注释unigenes总数的35.1%；与海枣(PhoenixdactyliferaLinn.)有12 768个unigenes相似，占注释unigenes总数的25.5%；与小果野蕉(MusaacuminataColla)的1个亚种Musaacuminatasubsp.malaccensis(Ridl.) N. W. Simmonds有4 658个unigenes相似，占注释unigenes总数的9.3%；与莲(NelumbonuciferaGaertn.)有1 556个unigenes相似，占注释unigenes总数的3.1%；与其他种类的unigenes相似性较低，占注释unigenes总数的26.9%。从种类的属性看，薯蓣为单子叶植物，与薯蓣转录组注释结果同源性较高的前3个种类也同为单子叶植物。

2.2.2 GO注释分析 在GO数据库对具茎生珠芽的叶腋样品(S1)与无茎生珠芽的叶腋样品(S2、S3和S4)的转录组unigenes功能进行两两比较，结果见图2。结果显示：分别有6 098、8 840和9 141个差异unigenes注释在细胞组分、分子功能和生物过程3大类43个亚类中。

： 样品S1与样品S2的差异unigenes富集通路Enrichment pathway of differential unigenes between sample S1 and sample S2； ： 样品S1与样品S3的差异unigenes富集通路Enrichment pathway of differential unigenes between sample S2 and sample S3； ： 样品S1与样品S4的差异unigenes富集通路Enrichment pathway of differential unigenes between sample S1 and sample S4. S1，S2，S3. 具茎生珠芽的薯蓣品种‘铁棍山药’的叶腋样品Leaf axil samples from Dioscorea polystachya ‘Tiegun’with stem bulbils： S1. 具茎生珠芽的叶腋样品Leaf axil sample with stem bulbils； S2. 茎基部无茎生珠芽的叶腋样品Leaf axil sample without stem bulbil on base part of stem； S3. 茎上部无茎生珠芽的叶腋样品Leaf axil sample without stem bulbil on upper part of stem. S4. 无茎生珠芽的薯蓣品种‘花籽山药’的叶腋样品Leaf axil sample from D. polystachya ‘Huazi’ without stem bulbil.

在生物过程的19个亚类中，3组差异unigenes主要注释在细胞过程、代谢过程和单一生物过程3个亚类中。其中，注释在细胞过程的差异unigenes分别有569、844和864个，各占对应注释unigenes总数的9.3%、9.5%和9.5%；注释在代谢过程的差异unigenes分别有647、926和978个，各占对应注释unigenes总数的10.6%、10.5%和10.7%；注释在单一生物过程的差异unigenes分别有417、575和619个，各占对应注释unigenes总数的6.8%、6.5%和6.8%。

在细胞组分的14个亚类中，3组差异unigenes主要注释在细胞、细胞部分和细胞器3个亚类中。其中，注释在细胞的差异unigenes分别有588、838和825个，各占对应注释unigenes总数的9.6%、9.5%和9.0%；注释在细胞部分的差异unigenes分别有588、838和825个，各占对应注释unigenes总数的9.6%、9.5%和9.0%；注释在细胞器的差异unigenes分别有470、681和653个，各占对应注释unigenes总数的7.7%、7.7%和9.3%。

在分子功能的10个亚类中，3组差异unigenes主要注释在2个亚类中：一是结合亚类，分别有538、809和786个差异unigenes，各占对应注释unigenes总数的8.8%、9.2%和8.6%；二是催化活性亚类，分别有628、849和960个unigenes，各占对应注释unigenes总数的10.3%、9.6%和10.5%。

2.2.3 KEGG代谢通路分析 在KEGG数据库中对具茎生珠芽的叶腋样品与无茎生珠芽的叶腋样品的转录组unigenes功能进行两两比较，结果见图3。结果显示：3组转录组分别有11 182、13 552和15 918个差异unigenes注释在KEGG代谢通路，共5大类18亚类，且主要集中在代谢和遗传信息处理。

： 样品S1与样品S2的差异unigenes富集通路Enrichment pathway of differential unigenes between sample S1 and sample S2； ： 样品S1与样品S3的差异unigenes富集通路Enrichment pathway of differential unigenes between sample S2 and sample S3； ： 样品S1与样品S4的差异unigenes富集通路Enrichment pathway of differential unigenes between sample S1 and sample S4. S1，S2，S3. 具茎生珠芽的薯蓣品种‘铁棍山药’的叶腋样品Leaf axil samples from Dioscorea polystachya ‘Tiegun’with stem bulbils： S1. 具茎生珠芽的叶腋样品Leaf axil sample with stem bulbils； S2. 茎基部无茎生珠芽的叶腋样品Leaf axil sample without stem bulbil on base part of stem； S3. 茎上部无茎生珠芽的叶腋样品Leaf axil sample without stem bulbil on upper part of stem. S4. 无茎生珠芽的薯蓣品种‘花籽山药’的叶腋样品Leaf axil sample from D. polystachya ‘Huazi’ without stem bulbil.

注释在代谢相关通路的差异unigenes数最多，分别有7 100、8 265和9 882个，各占对应注释unigenes总数的63.5%、61.0%和62.1%，涉及碳水化合物代谢、脂质代谢、氨基酸代谢和其他次生代谢产物的生物合成等11个亚类。

注释在遗传信息处理相关通路的差异unigenes分别有1 969、2 925和3 220个，各占对应注释unigenes总数的17.6%、21.6%和20.2%，涉及翻译、转录、复制和修复以及折叠、分类和降解4个亚类。

在环境信息处理相关通路中，分别有459、521和612个差异unigenes注释在信号转导亚类，各占对应注释unigenes总数的4.1%、3.8%和3.8%。在细胞过程相关通路中，分别有691、829和966个差异unigenes注释在运输和分解代谢亚类，各占对应注释unigenes总数的6.2%、6.1%和6.1%。在生物系统相关通路中，分别有757、816和951个差异unigenes注释在环境适应性亚类，各占对应注释unigenes总数的6.8%、6.0%和6.0%。

2.3 基因表达分析

2.3.1 unigenes的表达量比较 对‘铁棍山药’具茎生珠芽的叶腋样品(S1)与其他3个无茎生珠芽的叶腋样品(S2、S3和S4)的转录组unigenes的表达量进行比较，结果显示：与‘花籽山药’的叶腋样品相比，表达量上调和下调的unigenes分别有12 258和7 597个；与茎上部无茎生珠芽的叶腋样品相比，表达量上调和下调的unigenes分别有10 679和5 610个；与茎基部无茎生珠芽的叶腋样品相比，表达量上调和下调的unigenes分别有5 938和7 386个。其中，前两组样品中表达量上调的unigenes明显多于表达量下调的unigenes数量，前者是后者的1.6～1.9倍；而后一组样品中表达量上调的unigenes则少于表达量下调的unigenes。

将‘铁棍山药’具茎生珠芽的叶腋样品与其他3组无茎生珠芽的叶腋样品进行总体比较，结果显示：表达量均上调的unigenes有3 531个，表达量均下调的unigenes有1 084个。在具茎生珠芽的叶腋样品中特异性表达的unigenes有235个，而在3组无茎生珠芽的叶腋样品中均特异性表达的unigene仅1个。

2.3.2 PEBP家族基因的鉴定和表达量分析 在薯蓣品种‘铁棍山药’和‘花籽山药’不同部位样品的转录组数据中共鉴定出11个PEBP家族成员，将这些PEBP家族成员与拟南芥、水稻和几内亚薯蓣的基因组数据进行同源性比对，依次比对出同源基因6、19和11个。

薯蓣转录组与拟南芥、水稻和几内亚薯蓣基因组中PEBP家庭成员的系统演化树见图4。由图4可以看出：供试2个薯蓣品种的PEBP家族基因主要分布在TFL1-like和FT-like2个分支上，其中编号为Unigene40804_All的基因(命名为DpTFL1)与拟南芥中已知功能的TFL1基因(基因编号为At5G03840.1)以较高的支持率(88%)聚为一支。

Os： 水稻基因组Oryza sativa Linn. genome； AT： 拟南芥基因组Arabidopsis thaliana (Linn.) Heynh. genome； Dr： 几内亚薯蓣基因组Dioscorea rotundata Poir. genome； CL，Unigene： 薯蓣转录组Dioscorea polystachya Turcz. transcriptome.

‘铁棍山药’和‘花籽山药’不同部位样品中PEBP家族成员的表达丰度见表2。由表2可以看出：分布在TFL1-like分支的5个基因(基因编号分别为CL1415.Contig6_All、CL1415.Contig3_All、Unigene40804_All、Unigene8060_All以及CL1415.Contig7_All)均在‘铁棍山药’具茎生珠芽的叶腋样品中表达丰度最高，在其他3个无茎生珠芽的叶腋样品中表达丰度均较低。另外，Unigene40804_All(即DpTFL1基因)只在‘铁棍山药’具茎生珠芽的叶腋样品中特异性表达，推测该基因可能与茎生珠芽的发生有关。

表2 薯蓣品种‘铁棍山药’和‘花籽山药’不同部位样品中PEBP家族成员的表达丰度

2.3.3 在茎生珠芽不同发育期DpTFL1基因的表达量变化 扩增结果显示：使用DpTFL1基因的正向引物qDpTFL1-F和反向引物qDpTFL1-R能够扩增出单一条带，获得的目的片段长度为187 bp；使用内参正向引物DpActin-F和反向引物DpActin-R也能够扩增出单一条带，获得的目的片段长度为190 bp。

采用qRT-PCR技术对DpTFL1基因在‘铁棍山药’茎生珠芽4个不同发育期的相对表达量进行分析，结果(图5)显示：DpTFL1基因相对表达量在‘铁棍山药’的茎生珠芽未发育期最低，但随发育过程推进，DpTFL1基因表达水平逐渐升高，DpTFL1基因的相对表达量在茎生珠芽发育初期、中期和成熟期均显著高于茎生珠芽未发育期；DpTFL1基因的相对表达量在茎生珠芽发育成熟期达到最高，为茎生珠芽发生初期的10倍。表明DpTFL1基因是薯蓣茎生珠芽特异性表达基因，且其表达量随珠芽的发育而变化。

T0： 茎生珠芽未发育期Undeveloped stage of stem bulbil； T1： 茎生珠芽发育初期Early development stage of stem bulbil； T2： 茎生珠芽发育中期Middle development stage of stem bulbil； T3： 茎生珠芽发育成熟期Mature development stage of stem bulbil.

3 讨 论

上述研究结果表明：薯蓣品种‘铁棍山药’的具茎生珠芽叶腋部位与其茎上部无茎生珠芽叶腋部位以及无茎生珠芽品种‘花籽上药’的叶腋部位相比，表达量上调的unigenes数明显多于表达量下调的unigenes数，说明大部分差异unigenes在茎生珠芽的起始和形成过程中上调，提示茎生珠芽的启动过程可能需要更多基因的上调表达。前人在白花虎眼万年青和卷丹的叶面珠芽形成过程中也发现上调unigenes数多于下调unigenes数[11，16]。值得注意的是，Weber等[27]和姜福星等[28]的研究结果表明：注释为BTB/POZ基因的unigenes表达量在具叶面茎生珠芽的叶腋部位明显上调，该基因能够调控其下游基因表达，进而促进细胞的增殖、分裂和分化进程。由于在白花虎眼万年青的叶面珠芽发育过程中BTB/POZ基因的表达量也显著上调[28]，因此，该基因可能与启动珠芽发生相关。另外，在‘铁棍山药’的具茎生珠芽叶腋部位显著上调表达的基因还有蔗糖合成酶基因SUS(另文发表)，该基因是淀粉合成通路中的重要基因，而淀粉合成在珠芽的起始过程中起重要作用，多糖能够在茎生珠芽细胞的成熟过程中积累并促使茎生珠芽的形成[6]，而SUS基因在卷丹茎生珠芽发育过程中也上调表达[16]。因此，推测与淀粉合成相关的SUS等基因与薯蓣茎生珠芽发生相关。

在‘铁棍山药’的具茎生珠芽叶腋部位，还发现多个与生长素、细胞分裂素和赤霉素等与激素相关的基因表达量明显上调，如生长素生物合成基因YUCCA1、生长素转运蛋白相关基因ABCB19和LAX2以及细胞分裂素转运基因ENT3等(另文发表)；这些基因的表达量在参薯的茎生珠芽发育期也显著上调[17]；在卷丹的茎生珠芽形成过程中，与赤霉素合成相关的4个GA2OX基因的表达量不同程度上调[15]。因此，在‘铁棍山药’的具茎生珠芽叶腋部位中，表达量显著上调的与激素相关基因也很可能参与了薯蓣茎生珠芽发育的调控。

在‘铁棍山药’具茎生珠芽叶腋样品中，特异性表达的unigenes数(235)明显多于3组无茎生珠芽的叶腋样品(1)，提示在茎生珠芽的形成过程中可能需要更多基因在特定部位表达。在‘铁棍山药’具茎生珠芽叶腋部位特异表达的基因中，包含一些与细胞壁修饰相关的基因，如木葡聚糖内糖基转移酶基因XET和果胶酶基因PE等，以及18个与糖合成和运输有关的基因，如蔗糖转化酶基因VIF1和海藻糖-磷酸合成酶基因TPS等(另文发表)。在参薯珠芽发育初期和马铃薯(SolanumtuberosumLinn.)块茎发生过程中，这些与细胞壁修饰以及糖合成和运输有关的基因也显著上调表达[17，29]；而在‘铁棍山药’和参薯茎生珠芽的发育过程中可溶性糖含量呈曲线上升[6，17]。说明这些参与糖代谢的相关基因可能也参与了薯蓣茎生珠芽发育的调控。

在‘铁棍山药’具茎生珠芽叶腋部位特异性表达的与侧枝发生相关的unigenes中，Unigenes40804_All被注释为PEBP家族的TFL1-like同源基因(命名为DpTFL1基因)，其表达丰度最高。TFL1-like基因可能与植物的茎生珠芽发生、侧枝形成和块茎生长相关，如：在杭白菊中CmTFL1c基因的过表达可导致其侧枝数明显增多[20]；在半夏茎生珠芽的生长发育过程中PteTFL1基因可能发挥重要作用[21]；而马铃薯的StTFL1基因则可通过调控糖代谢途径诱导促进马铃薯块茎的形成[30]等。据此推测DpTFL1基因可能与薯蓣茎生珠芽发育密切相关，这一结论仍需后续实验验证。此外，在3组无茎生珠芽的叶腋样品中特异性表达的1个unigene的表达丰度并不高(FPKM≤1.09)，且该基因在7大功能数据库中均未被注释为已知功能基因，该基因的功能仍需进一步研究。

研究结果显示：随‘铁棍山药’茎生珠芽发育，DpTFL1基因表达量逐渐升高，并在茎生珠芽发育成熟期达到最高。实际上TFL1-like基因在不同植物的不同分生组织发育过程中的表达模式存在一定差异，例如：半夏的PteTFL1基因在半夏茎生珠芽原基形成期表达量高，而在其茎生珠芽原基分化期、膨大期和成熟期表达量则有所减少[21]；在白梨(PyrusbretschneideriRehd.)的花器官发育过程中，pbTFL1基因在花芽分化初期表达量最高，在花启动期明显下降，在花器官原基开始时降至很低，之后表达量不再升高并保持低水平[31]；而金鱼草(AntirrhinummajusLinn.)的CEN基因(TFL1同源基因)在花序分生组织中表达量呈现逐渐升高的趋势，且仅在花序发育后期表达量较高，在营养生长期则不表达[32]。Goretti等[33]认为，TFL1-like基因编码的蛋白能够通过14-3-3蛋白与bZIP转录因子FD形成复合物，并与特定的染色质结合调控下游基因的表达，如LFY(LEAFY)、AP1(APETALA1)以及FDL1(FD-like1)基因等，从而调控植物的营养生长与生殖生长，并决定分生组织的命运。

此外，TFL1-like基因调控的下游基因也包括与多种激素通路相关的基因，以及与糖合成、运输相关的基因，且与多种激素通路相关联[34]。但不同植物在各类分生组织发育过程中相关激素的合成和转运以及糖的合成、运输和储存机制并不相同，导致在不同的植物中TFL1-like基因表达量的变化规律可能存在异同。本研究仅仅阐明了DpTFL1基因在‘铁棍山药’茎生珠芽发育过程中的表达模式，但该基因是否调控茎生珠芽的发生及其调控机制还有待进一步的实验研究。