杜仲雄花芽2个发育时期转录组分析

朱利利 杜庆鑫 何 凤 庆 军 杜红岩*

(1.中国林业科学研究院经济林研究开发中心，郑州 450003； 2.国家林业局杜仲工程技术研究中心，郑州 450003)

杜仲(Eucommiaulmoides)为杜仲科(Eucommiaceae)的单科单属植物，是第四季冰川侵袭后仅留存于我国的孑遗植物,具有重要的经济价值、药用价值和生态价值，广泛分布于我国27个省(市、区)[1]。杜仲为严格的雌雄异株植物，雄株所占比例接近50%，雄花结构简单，为8～12枚针状雄蕊组成一个雄蕊群簇生长于苞片腋间，产量高、易采集[2～3]。杜仲雄花富含京尼平苷酸、绿原酸、黄酮等多种活性成分，同时含有多种人体必需氨基酸，具有较高的药用和营养价值[4～6]。目前，杜仲雄花已作为多种饮品原料进行开发，如杜仲雄花芽、杜仲雄花酒、杜仲雄花功能饮料等，具有广阔的市场开发前景[7]。

植物花芽分化过程是内外因素共同作用的结果，分化数量和质量直接决定花和果实的数量和品质，影响产品开发和商品价值[8]。因此，杜仲雄花芽分化结果特别是雄蕊原基分化对雄花的产量和质量起决定性作用。目前，关于杜仲雄花芽分化过程和时间已进行了观察，但关于调控杜仲雄花芽分化的基因研究较为匮乏。利用高通量测序技术进行转录组学研究，可了解物种在特定环境条件下基因功能表达、生物过程和分子运行机制。本研究前期观察了杜仲雄花芽分化过程，可分为花序原基分化期、苞片原基分化期、雄蕊原基分化初期、雄蕊形态建成期4个时期，与前人研究结果较一致[9]。根据前期结果，本研究以苞片原基分化期和雄蕊原基分化期的雄花芽为材料，采用Illumina Hiseq X-ten平台对杜仲雄花芽进行转录组测序，并对测序的数据进行过滤，过滤数据以杜仲基因组为参考进行序列比对，利用StringTie进行组装。根据基因组注释信息，查询未被注释的转录区，挖掘新基因，同时根据常用数据库信息对新基因进行注释。根据转录组信息，了解杜仲雄花芽在苞叶原基发育期和雄蕊原基分化初期的基因表达情况，挖掘调控雄蕊原基发育的基因，为促进雄花发育及分子育种提供重要的数据参考。

1 材料与方法

1.1 试验材料

试验材料采集于中国林业科学研究院经济林研究中心原阳试验基地(34°55′18″～34°56′27″N，113°46′14″～113°47′35″E)栽培的杜仲良种雄株“华仲11号”(“Huazhong No.11”)。选择3株长势良好、无病虫害、长势一致的植株，根据杜仲花芽分化结果分别于5月30日(苞叶原基分化中期)，6月14日(雄蕊原基分化初期)采集花芽(图1)，每株树采集30个花芽作为一个生物学重复，共3个生物学重复。花芽采集后立即放入液氮中，然后于-80℃超低温冰箱保存备用。样品命名方式：M11-1-1表示苞叶原基分化期雄花芽，M11-2-1表示雄蕊原基分化初期雄花芽，最后一个数字表示生物学重复。

1.2 cDNA文库构建和转录组测序

杜仲雄花芽在液氮中充分研磨后，根据TRIzol试剂盒(Invitrogen)操作步骤提取RNA，其中RNA的浓度、纯度和完整性分别用NanoDrop和Agilent 2100检测。质检合格的RNA，根据cDNA文库构建试剂盒(NEB)说明书构建cDNA文库。利用lllumina Hiseq X-10平台对cDNA文库进行双端测序，测序读长为150 bp，获得原始数据。

1.3 数据组装和新转录本预测

根据转录组分析标准，测序得到的Raw Data去除质量不合格序列后得到高质量的clean reads，用于后续分析。使用杜仲基因组作为参考和分析，利用软件HISAT2将clean reads与基因组进行比对[10]。比对分析完成后利用StringTie将比对上的Reads进行组装，且与杜仲基因组注释信息进行比较(NCBI登录号PRJNA357336、SRP095726及SRS2666014)， 查询未被注释的转录区，挖掘杜仲新转录本，补充完善基因组注释信息[11]。使用BLAST[12]软件将挖掘的新基因与NR[13]，Swiss-Prot[14]，GO[15]，COG[16]，KOG[17]，Pfam[18]和KEGG[19]数据库进行序列比对，使用KOBAS2.0[20]得到新转录本的KEGG Orthology结果，完成新转录本的氨基酸序列预测后使用HMMER[21]软件与Pfam数据库比对，获得新基因的注释信息。

图1 杜仲雄花芽2个花芽分化时期形态特征 Br.苞片；SAM.茎尖分生组织；Sq.鳞片；Sta.雄蕊Fig.1 Morphological characteristics of male floral bud at two differentiated stages Br.Bract; SAM.Stem apical meristem; Sq.Squama; Sta.Stamen

1.4 基因表达量定量和差异基因筛选

转录本组装完成采用FPKM(Fragments Per Kilobase of transcript per Million fragments mapped)计算基因表达水平[22]。使用Ballgown对不同分化阶段的花芽差异表达基因进行计算和统计[22]。差异基因筛选条件设置为|log2foldchange|>1且校正P-value(FDR)<0.01。

2 结果与分析

2.1 杜仲雄花芽转录组测序数据统计

杜仲6个雄花芽转录组测序(表1)，原始数据经过滤后共获得40.48 Gb clean reads,各样品clean reads均达到6.09 Gb及以上，Q30碱基百分比在90.75%及以上，GC平均含量为46.86%，数据质量较高可用于后续分析。分别将各样品的clean reads与杜仲基因组进行序列比对，比对效率为91.59%～92.05%。

表1 测序数据统计

图2 新转录本功能注释结果统计Fig.2 Summary of annotated new transcripts in seven databases

图3 杜仲雄花芽不同发育时期表达基因和差异表达基因 a.2个发育时期表达基因韦恩图；b.2个发育时期差异表达基因火山图Fig.3 Identified and differentially expressed genes at two differentiated stages of E.ulmoides a. Expressed genes at two differentiated stages; b. Differentially expressed genes at two differentiated stages

图4 差异基因GO功能富集 生物学过程(A.代谢过程；B.细胞过程；C.单个有机体过程；D.生物调节；E.定位；F.刺激应答；G.细胞成分组织或生物发生；H.发育过程；I.信号；J.多细胞生物的过程；K.生殖；L.生殖过程；M.多有机体过程；N.解毒；O.增长；P.免疫系统过程；Q.运动；R.细胞杀伤；S.生物节奏的过程；T.生物附着)；细胞组分(A.细胞；B.细胞部分；C.膜；D.细胞器；E.膜部分；F.细胞器部分；G.大分子复合体；H.胞外区；I.膜密封腔；J.细胞连接；K.超分子复合物；L.细胞外区域部分；M.病毒体；N.病毒部分；O.类核)；分子功能(A.催化活性；B.绑定；C.载体活性；D.结构分子活性；E.核酸结合转录因子活性；F.分子功能调节；G.信号传感器活性；H.电子载体活性；I.抗氧化活性；J.分子传感器活性；K.营养库活性；L.转录因子活性、蛋白质绑定；M.金属伴侣活性；N.蛋白质标记)Fig.4 GO analyses of differentially expressed genesBiological process(A.Metabolic process; B.Cellular process; C.Single-organism process; D.Biological regulation; E.Localization; F.Response to stimulus; G.Cellular component organization or biogenesis; H.Developmental process; I.Signaling; J.Multicellular organismal process; K.Reproduction; L.Reproductive process; M.Multi-organism process; N.Detoxification; O.Growth; P.Immune system process; Q.Locomotion; R.Cell killing; S.Rhythmic process; T.Biological adhesion); Cellular component(A.Cell; B.Cell part; C.Membrane; D.Organelle; E.Membrane part; F.Organelle part; G.Macromolecular complex; H.Extracellular region; I.Membrane-enclosed lumen; J.Cell junction; K.Supramolecular complex; L.Extracellular region part; M.Vvirion; N.Virion part; O.Nucleoid; P.Catalytic activity); Molecular function(A.Catalytic activity; B.Binding; C.Transporter activity; D.Structural molecule activity; E.Nucleic acid binding transcription factor activity; F.Molecular function regulator; G.Signal transducer activity; H.Electron carrier activity; I.Antioxidant activity; J.Molecular transducer activity; K.Nutrient reservoir activity; L.Transcription factor activity,protein binding; M.Metallochaperone activity; N.Protein tag)

2.2 新转录本发掘和注释

基于clean reads在基因组比对结果，进行新转录本的发掘，共发掘新转录本2 929个,将发掘的新转录本与NR，Swiss-Prot，GO，COG，KOG，Pfam，KEGG 7大数据库进行序列比对，共2 166个新转录本得到功能注释，占总得到新转录本的73.95%(图2)。其中，注释到NR蛋白数据库中的新转本数量最多为2 136个占总新转录本数72.93%，然后依次注释到Swiss-Prot数据库数量为1 509个占51.52%、GO数据库1 336个占45.61%，Pfam数据库1 235个占42.16%，KOG数据库为1 174个占40.08%，KEGG数据库865个占29.53%，而在COG数据库注释数量最少为571个占19.49%。

2.3 杜仲雄花芽2个发育时期差异基因表达分析

杜仲雄花芽2个发育时期共鉴定获得20 226个表达基因，2个发育时期共表达基因18 094个，苞叶基分化期特异表达778个，雄蕊原基分化初期特异表达1 354个(图3a)。结果显示约89.5%基因在杜仲雄花芽2个发育阶段均表达，仅约3.8%和6.7%基因在苞叶原基和雄蕊原基发育阶段特定表达。根据差异基因筛选条件，在2个发育时期筛选出583个差异表达基因，其中在雄蕊原基发育初期上调基因315个，下调基因267个(图3b)，可能是这些特异表达基因特别是差异表达基因促进了杜仲雄花芽的雄蕊原基萌发。

2.4 杜仲雄花芽差异基因GO功能分析

为了进一步挖掘参与杜仲雄蕊发育的基因，对雄花芽苞叶原基发育时期和雄蕊原基发育初期的差异基因进行GO功能富集分析，筛选出雄花芽中与成花诱导紧密相关的GO条目(图4)。雄花芽在苞叶原基发育期和雌蕊原基发育初期之间共有359个差异基因富集到个1 388个GO条目(GO terms)，其中显著性富集GO条目(P-value<0.05)有91个，其中有52个GO条目显著富集到生物学过程，11个显著富集到细胞组分，28个显著性富集到分子功能。雄花芽显著性富集的GO条目包括与分生组织发育相关的过程如分生组织从营养生长向生殖生长的转变(vegetative to reproductive phase transition of meristem)、生长调节(regulation of growth)、生长负向调控(negative regulation of growth)，与激素相关的激素生物合成过程如吲哚硫代葡萄糖苷代谢过程(indole glucosinolate metabolic process)、乙烯代谢过程(ethylene metabolic process)和乙烯生物合成过程(ethylene biosynthetic process)，与光周期相关的远红光反应(response to far red light)，光周期和开花(photoperiodism,flowering)，光周期现象(photoperiodism)。同时差异基因被富集到核酸结合转录因子活性(nucleic acid binding transcription factor activity)、发育过程(developmental process)，生殖(reproduction)和生殖过程(reproductive process)等与开花紧密相关的生物过程。说明杜仲雄花芽中参与营养生长向生殖生长的转变、激素合成与代谢、光周期过程、核酸结合转录因子活性、生殖过程等基因的差异表达可能是促进杜仲雄花芽由苞叶原基发育向雌、雄蕊原基发育转变的直接促进因子。

2.5 杜仲雄花芽差异基因KEGG功能分析

杜仲雄花芽2个发育时期有98个差异表达基因被富集到70条KEGG pathway。其中富集比例前5的pathway分别是苯丙素生物合成富集比例11.22%，其次是淀粉和蔗糖代谢富集比例11.22%，碳代谢富集比例10.20%，内质网蛋白质加工富集比例10.20%，氰基氨基酸代谢富集比例8.16%。同时差异基因被富集到植物激素信号转导富集比例为3.06%，以及显著富集(P-value<0.05)到植物生物节律富集比例6.12%(表2)。代谢通路富集结果说明杜仲雄花芽在雄蕊发育过程中需要碳水化合物作为基础，且需要激素和一些次生代谢物调控，同时杜仲雄花芽形态分化与植物节律密切相关。

表2 差异表达基因代谢通路注释

***P<0.001

2.6 开花相关差异基因分析结果

根据GO和KEGG富集分析，营养生长向生殖生长的转变、激素合成与代谢、激素信号转导、光周期过程和植物生物节律，核酸结合转录因子活性等生物过程与杜仲雄花芽发育密切相关，对调节这些生物途径的差异基因进行详细分析。分析结果显示，与激素相关的差异表达基因有16个，其中与生长素(IAA)相关的有6个差异基因，有4个上调差异基因分别是ARF9(EUC06536-RA)、GH3.10(EUC06844-RA)、ABCB19(EUC13060-RA)、PILS5(EUC02516-RA)和下调基因PCNT115(EUC23432-RA和EUC25737-RA)；与乙烯(ETH)代谢和传导途径相关的有9个差异基因，包含上调因子RAP2-7(EUC00521-RA)、AFP3(EUC_newGene_23655)、AHK1(EUC08386-RA)、EREBP-like(EUC00228-RA)，FER(EUC20697-RA)、AFP2(EUC15307-RA)、ACO(EUC22557-RA)、CTR1(EUC03132-RA)和一个下调因子ERF13(EUC_newGene_3252)；而与脱落酸(ABA)相关的仅有下调因子PP2C基因(EUC22646-RA)。

杜仲雄花芽2个发育时期发现了调控光周期开花途径重要转录因子COL差异表达，包括上调转录因子COL4(EUC12564-RA)和COL2(EUC15205-RA)，下调因子COL9(EUC17024-RA)，同时还发现与生物节律相关2个上调基因CCA1(EUC05190-RA)和光敏色素PHYA(EUC15838-RA)以及五个下调基因PRR73(EUC00514-RA)，APRR5(EUC16182-RA)，ADO3(EUC22871-RA)，GI(EUC19998-RA)和ELF4(EUC18826-RA)。此外，结果显示有多种转录因子家族成员差异表达调控杜仲雄花芽从营养生长向生殖生长转变，包括MADS-box中MIKC亚家族基因上调转录因子AGL3(EUC07615-RA)和AGL8(EUC07648-RA)，以及开花整合子FLC(EUC06393-RA)和SOC1(EUC03328-RA)；WRKY家族上调转录因子WRKY41(EUC07015-RA、EUC22632-RA和EUC07417-RA)、WRKY24(EUC00375-RA)；MYB家族上调转录因子MYB61(EUC06946-RA)、CCA1(EUC05190-RA)和RVE8(EUC03052-RA)，以及下调转录因子LUX(EUC07116-RA和EUC11242-RA)；GATA家族上调因子GATA18(EUC11132-RA)和下调因子GATA9(EUC18806-RA)；同时还包括bHLH、TCP、PIF1以及新转录因子家族PLATZ参与雄花芽的雄蕊发育过程。在本研究中，杜仲雄花芽转录组差异分析结果还显示与植物分生组织相关的基因AGO2和AGO7在雄花芽苞叶原基分化期和雄蕊原基分化初期差异表达。

3 讨论

本研究以杜仲雄花良种“华仲11号”为材料，对处于苞叶原基发育时期和雄蕊原基发育时期的花芽进行转录组测序，功能富集结果显示差异表达基因显著富集在淀粉和蔗糖代谢途径和碳代谢途径。研究表明糖类在成花过程中可作为能量物质和结构物质促进花芽分化，而且可以作为成花信号物质直接作用于茎尖调控花芽分化[23]。植物成花过程中除了需要以碳水化合物作为基础外，而且与激素调控密切相关。大量研究证明激素参与植物的开花过程，对植物花芽分化有显著调控作用。生长素极性运输对花的形成具有重要作用，其中拟南芥(Arabidopsisthaliana)pin-1突变体的花序轴的内源生长素和生长素的极性运输能力比起野生型明显降低，其花芽不能正常发育成完整的花[24]。本文结果显示有15个差异基因与生长素、乙烯和ABA的生物代谢和传导途径相关，其中与乙烯和生长素相关的基因除ERF13外其他都在雄蕊原基发育阶段表现为上调，而与ABA相关的基因则表现为下调。其中与生长素响应元件结合的生长素响应因子(ARF9)，具有维持植物生长素动态平衡的GH3基因，生长素转运蛋白基因ABCB19和PILS5在杜仲雄蕊原基发育时期表达上调，而生长素诱导蛋白基因PCNT115则表达下调。乙烯具有刺激开花及花的衰老，其中参与乙烯信号传导的RAP2，AFP3，CTR1等是乙烯信号途径中一个重要的调节因子，在杜仲雄蕊原基发育时期表达上调。乙烯响应因子ERF在玉米中可以调控胚乳发育[25]，本研究显示ERF在杜仲雄蕊原基发育时期表达下调。对ABA调控花芽分化研究显示，低含量的ABA有利于板栗(Castaneamollissima)花簇原基的分化和生长[26]。ABA受体PYL蛋白是ABA信号转导通路中重要成员，可与ABA结合抑制PP2C表达[23]；在本研究中，PP2C在杜仲雄花芽苞叶原基时期表达上调，而在雄蕊原基发育期可能受到抑制表达下调。综合以上结果，生长素、乙烯可能正调控杜仲雄蕊原基的形态发育，而ABA可能负调控杜仲雄蕊原基发育。

光周期途径是植物成花诱导的重要途径，植物在感受适宜的光周期后，通过生物钟调控开花基因表达。本研究显示与光周期和物种生物钟密切相关的基因COL成员(COL4、COL2、COL9)、CCA1、PHYA、PRR73、APRR5、ADO3、GI和ELF4在杜仲雄花芽发育不同时期差异表达。其中，CO同源基因COL是光周期途径成花的关键基因，拟南芥ATCOL2表达受光照调控，在黎明时表达最高，傍晚时表达急剧下降，但在转基因植物中其表达量对开花时间无太大影响[27]；而从小麦(Triticumaestivum)中克隆的TACOL4转化到水稻(Oryzasativa)中可延迟水稻抽穗[28]；在玉米(Zeamays)中COL9可通过调控Ehd1途径延长玉米开花时间[29]。GI是拟南芥核心生物钟基因之一，调控下游CO基因的表达来影响植物开花时间和进程[30]。光敏色素基因PHYA调控CO表达，生物钟蛋白基因CCA1和LHY通过降低SVP蛋白累积加速拟南芥开花进程，此外早花基因ELF4、PRR73、APRR5等基因也是调控生物钟和开花过程的重要因子[31～32]。因此，光周期途径很可能是诱导杜仲成花的主要途径。

转录因子家族在植物花器官发育、果实发育以及逆境胁迫等方面发挥重要作用，其中MADS-box基因是重要的花器官属性决定基因。MADS-box家族成员FLC、SOC1、AGL24和SVP是植物成花诱导途径春化途径、自主途径、赤霉素途径和光周期途径中的关键节点作用基因[33]。本文研究发现了4个TypeⅡMADS-box型MIKC亚家族成员在杜仲雄花芽2个发育时期差异表达，其中AGL3和AGL8在苞叶原基发育期表达下调而在雄蕊原基发育初期表达上调，而SOC1、FLC与AGL基因表达正好相反。FLC是一个开花抑制因子，春化途径与自主途径中发挥作用，在拟南芥中过表达FLC则可推迟开花[34]。开花整合子SOC1受FLC的负调控抑制植物成花，同时也受CO信号在长日照条件下诱导植物成花，而且也能接受赤霉素的调控诱导植物开花[35～36]。在ABC(DE)模型中，SEP类基因被命名为E类，而AGL3基因为第四个SEP基因，影响着四轮花器官和胚珠发育[37]。此外，在拟南芥中花原基中发现了MADS-box基因AGL8，在营养生长阶段不表达，而在花序顶端分生组织、花梗和茎生叶中高表达[38]。所以，在杜仲雄花芽发育过程中SOC1和FLC在苞叶原基发育时期高表达调控花分生组织发育，然后SEP类基因被激活进而调控杜仲雄蕊原基发育，但这4个基因在杜仲雄蕊发育过程中的具体功能还需进一步实验验证。

4 结论

杜仲雄花芽2个发育时期获得583个差异表达基因，雄蕊原基发育初期上调基因315个，下调基因267个。差异基因的GO和KEGG富集结果显示，差异基因富集在营养生长向生殖生长的转变、乙烯激素合成与代谢，核酸结合转录因子活性，光周期现象等生物过程以及碳水化合物代谢、植物激素信号转导和植物生物节律等在花芽分化过程中起重要调控作用的代谢通路。根据差异基因分析结果，发现在杜仲雄花芽形态分化过程中除需要以碳水化合物作为基础外，还受生长素、乙烯和脱落酸等激素调控，同时光周期途径可能是杜仲成花诱导的重要途径。此外，MADS-box家族成员FLC、SOC1、AGL3和AGL8与杜仲雄蕊器官发育密切相关。本研究结果从转录水平上分析了杜仲雄花芽发育的相关代谢途径和基因，可在一定程度上解析杜仲雄花芽形态分化的分子调控模式与机制。