濒危树种香木莲转录组分析

苗艺明，石 松，杨 梅，刘世男

(1.广西大学林学院，广西 南宁 530004；2.广西壮族自治区都安瑶族自治县自然资源局，广西 都安 530799)

香木莲(Manglietiaaromatica)是木兰科木莲属常绿乔木，属国家Ⅱ级重点保护植物.该树种生长较快，高可达35 m，胸径可达1.9 m，主要分布于广西、云南和贵州海拔400～1 600 m地区.香木莲具有较强的适应性，在石灰岩发育的钙质土，以及土壤贫瘠、生境恶劣的溶洞边缘或石漠化严重的山坡能形成小群落或单株生长.香木莲全株都具香气且均可提取香油，用于调配名贵香料[1]；树体形态优美，花大香艳美丽[2]；木材较轻软、纹理细、抗虫蛀、耐腐蚀、不裂、不变形[3].香木莲是重要的香料，是园林绿化观赏及优良的用材树种.长期以来深受当地群众喜爱，其野生资源几乎被砍伐殆尽，现存野生种群十分稀少，且多为老树，幼树、幼苗少，自然更新能力差.PAN等[4]揭示了香木莲的濒危原因，结果显示，其花粉萌发能力低且传粉条件较差，大孢子、雌配子退化和败育，这些都大大降低了香木莲的结实率.种苗繁育研究对濒危物种扩大种群数量具有重要意义.有性繁殖研究发现，香木莲种皮含油质，致使种子发芽率和保存率低，从而限制了种子苗来源[5]；无性繁殖主要是扦插繁育，ABT生根粉能提高扦插生根率[6].研究幼苗在低温和干旱胁迫下的生理特性发现，香木莲具有一定的抗寒能力，但耐干旱能力较弱[7-8].目前，有关香木莲种苗繁育和濒危原因的研究只有少量报道，而分子生物学方面的研究未见报道.本研究在高通量测序的基础上，利用生物信息学方法分析香木莲转录组序列信息，为后续的香木莲遗传多样性、种质资源保存利用和功能基因挖掘等提供基础数据.

1 材料与方法

1.1 试验材料

本次研究的试验材料采自广西壮族自治区百色市都安县三只羊乡上远村.收集香木莲嫩叶、嫩枝和花，用锡箔纸包裹迅速放入液氮中保存，带回实验室，置于-80 ℃冰箱备用.

1.2 文库构建及测序

使用试剂盒提取嫩叶、嫩枝和花样本的总RNA，质量检验合格后，等量混合不同组织的RNA样品.香木莲cDNA文库构建及测序由杭州科睿迪有限责任公司完成.

1.3 转录组序列组装、注释及分析

RNA-seq测序完成后，过滤带有接头、低质量、冗余的Raw reads得到高质量clean reads数据，统计数量、长度、N%、Q20及GC%.利用Trinity software进行denovo组装，从长到短排序，依次累加不小于总长50%(N50)的拼接转录本长度，统计各文库subreads.将获得的Unigenes分别与NR、GO、KOG、KEGG、SwissProt等数据库进行比对，分析相应的功能注释.

2 结果与分析

2.1 香木莲RNA-seq与组装

通过转录组测序，共得到36 737 304条原始序列，过滤处理后获得35 321 846条有效序列，占总数的96.15%，Q20、Q30序列分别占总数的98.29%和95.07%，GC含量占总数的48.27%，碱基错误率为0.02%.以上结果说明，通过高通量测序平台获得的香木莲序列数量和质量较高，可以用于后续的相关生物信息学分析.

利用Trinity软件组装处理后的片段，共获得48 123条Unigene，全部碱基数达46 188 480 bp.组装后的Unigene长度分布见图1.结果可知：Unigene平均长度为960 nt，N50为1 331 nt.序列长度为200～500 nt的有15 932条，占33.1%；500～1 000 nt的有15 089条，占31.4%；1 000～1 500 nt的有8 098条，占16.85%；1 500～2 000 nt的有4 553条，占9.5%；大于等于2 000 nt的有4 451条，占9.2%.

2.2 香木莲转录组Unigene功能注释

利用Blast软件将香木莲48 123条Unigene与各数据库进行比对，结果显示：比对到NR、GO、SwissProt、KEGG和KOG数据库的Unigene分别有37 877、32 125、27 988、30 143和37 199条，占比依次为78.7%、66.8%、58.2%、62.6%和77.3%.

2.2.1 香木莲转录组Unigene的NR功能注释

利用Blast软件将香木莲全部Unigene与NR数据库进行比对，结果见图2.由图2可见：在匹配的近缘物种中，香木莲与莲花(Nelumbonucifera)同源序列最多，为10 427条，占总数的27.5%；其次为博落回(Macleayacordata)，有6 084条，占总数的16.1%；葡萄(Vitisvinifera) 2 287条、棕榈(Elaeisguineensis) 1 993条、海枣(Phoenixdactylifera)1 622条、洛矶山耧斗菜(Aquilegiacoerulea)1 465条、无油樟(Amborellatrichopoda)854条、核桃(Juglansregia) 553条、栓皮槠(Quercussuber)521条、小果野芭蕉(Musaacuminatasubsp.malaccensis)516条、芦笋(Asparagusofficinalis) 501条、菠萝(Ananascomosus) 435条、橡胶树(Heveabrsiliensis)427条、甜橙(Citrussinensis) 346条，分别占总数的6.0%、5.3%、4.3%、3.9%、2.3%、1.5%、1.4%、1.4%、1.3%、1.2%、1.1%和0.9%；其余9 846条分布于其他物种中，占总数的26%.

图1香木莲转录组组装序列长度分布Fig.1Length distribution of transcriptome Unigenes for Manglietia aromatica

图2香木莲转录组Unigene注释匹配的物种分布Fig.2Species distribution of Manglietia aromatic with Unigenes annotation

2.2.2 香木莲转录组Unigene的GO功能注释

香木莲转录组Unigene的GO功能注释见图3.结果表明：32 125条Unigene共获得246 974个GO功能注释，分成生物学过程、分子功能以及细胞组分3大类.其中，生物学过程获得了最多的注释，有107 579个，占43.6%；其次是细胞组分，获得98 271个注释，占39.8%；第3为分子功能，获得41 124个注释，占16.6%.3大功能进一步又分成58个亚类，其中，生物学过程包含28个亚类，细胞过程、代谢过程亚类获得的注释偏多，分别占该类型的20.0%和17.0%，碳利用率、细胞失活、行为以及硫利用率亚类注释最少，均占该类型的0.01%及以下；细胞过程包含18个亚类，其中，细胞、细胞组分亚类获得的注释偏多，均占该类型的22.2%；分子功能包含12个亚类，其中，结合、催化活性亚类获得的注释偏多，分别占该类型的47.1%和39.3%，翻译调节器活性和蛋白标亚类获得的注释较少，分别占该类型的0.02%和0.01%.

图3香木莲转录组Unigene的GO功能分类Fig.3GO functional annotation of transcriptome Unigenes for Manglietia aromatica

2.2.3 香木莲转录组Unigene的KOG功能注释

香木莲转录组Unigene序列KOG蛋白数据库分类注释见图4.结果显示：有37 199条Unigenes能够匹配在KOG数据库中，共获得25 525条注释，可分为25个功能大类.其中，一般功能预测基因最多，有4 542条，占总数的17.8%；其次是信号传导机制，有2 999条，占总数的11.8%；再次是翻译后修饰、蛋白质转换、伴侣有2 492条，占9.8%；转录功能有1 412条，占5.5%；碳水化合物运输和代谢有1 314条，占5.2%；细胞内、分泌、囊泡运输有1 293条，占5.1%；胞外结构和细胞运动所占比例最少，分别占0.3%和0.05%.

图4香木莲转录组Unigene的KOG注释分类Fig.4KOG functional classification of transcriptome Unigenes for Manglietia aromatica

2.2.4 香木莲转录组Unigene的KEGG代谢通路分析

香木莲转录组Unigene的KEGG分类见图5.结果显示：共有30 143条Unigenes获得注释并涉及142个代谢通路.进一步将这142个代谢通路划分为5大类，包括代谢、遗传信息处理、环境信息处理、细胞过程以及生物系统相关通路.该5大类又分为19个亚类，其中，代谢相关通路中11个亚类，以全局和概述地图居多，占总数的33.3%，第2是碳水化合物代谢相关通路，占该通路的17.5%，萜类和聚酮代谢、多糖合成和代谢以及核苷酸代谢相关通路相对较少，占比均在3.0%以下；遗传信息处理相关通路有4个亚类，其中，翻译相关通路所占比例最大，占总数的38.7%，折叠、分选和降解处理相关通路占32.3%，复制和修复相关通路最少，占比仅为13.0%；环境处理相关代谢通路中包含2个亚类，信号传导通路明显居多，高达82.3%；细胞过程和生物系统相关通路都仅包括1个亚类.

图5香木莲转录组Unigene的KEGG分类Fig.5KEGG classification of transcriptome Unigenes for Manglietia aromatica

3 小结与讨论

转录组测序技术在植物分子标记开发、功能基因挖掘及鉴定研究中起着重要作用[9-10]，已在木兰科树种中广泛应用.已有研究利用RNA-Seq技术对多种植物进行了测序，结果显示：红花玉兰(Magnoliawufengesis)共获得94 805条Unigene，平均长度为695 nt，N50为1 038 nt[11]；景宁玉兰(M.sinostellata)获得52 441条Unigene，平均长度为648 nt，N50为1 126 nt[12]；乐东拟单性木兰(Parakmerialotugensis)获得273 252条Unigene，平均长度为590 nt，N50为752 nt[13]；鹅掌楸(LiriodendronchinenseSarg)获得162 092条Unigene，平均长度为547 nt，N50为719 nt[14].本研究利用IIIumina HiSeqTM4000对香木莲转录组进行了测序，共获得48 123条Unigene，平均长度为960 nt，N50为1 331 nt.香木莲的N50和平均长度都高于上述木兰科树种.N50是评价组装序列完整性的重要指标，长度越长，代表组装的完整性越好[15].研究结果表明，香木莲测序获得的序列质量高且拼接完整性较好，有利于后续开展基因组方面的研究.

通过几种木兰科树种转录组Unigene与NR数据库比对发现，红花玉兰Unigene获得注释数较高的前3个物种为莲花、葡萄和可可树(Thenobromacacao)[11]；华木莲获得注释较高的前3个物种为莲花、油棕和海枣[16]；乐东拟单性木兰和景宁木兰获得注释数较高的前3个物种均为葡萄、海枣和可可树[12-13].而香木莲获得匹配率较高的前3个物种为莲花、博落回和葡萄，且博落回是以上树种都没有匹配到的.由此表明，香木莲与同科的树种存在相似功能基因，也可能存在特异功能基因.由香木莲GO功能注释分析结果可知，被注释的32 125条Unigene在功能上划分为3大类58个亚类，其中，注释到生物学过程大类的Unigene数量最多，主要是细胞过程和代谢过程；通过KEGG pathway分析可知，有30 143条Unigene获得注释，涉及5大类19个亚类共142个代谢通路，其中，以代谢相关通路和碳水化合物代谢相关通路为主.这一结果与华木莲、红花玉兰的Unigene GO功能注释和KEGG pathway分析结果类似[11，16]，都以生物学过程中的细胞过程和代谢过程为主.可见，香木莲在细胞和代谢活动的基因表达量相对较高，具有较强的代谢能力和丰富的生物过程.此外，分析玉兰的KEGG pathway发现，通路“类黄酮合成”富集44个Unigene，通路“花青素合成”富集9个Unigene，通路“黄酮和黄酮醇”富集6个Unigene[11]，这些通路可能参与花青素苷合成，进一步影响红花玉兰花色.本研究中的KEGG pathway分析表明，通路“类黄酮合成” “花青素合成” “黄酮和黄酮醇”分别富集140、4和19个Unigene.前期的野外调查发现，不同居群香木莲的花被片颜色有淡红色和白色两种，因此，通过测序分析香木莲转录组，有利于挖掘花青素苷合成的相关通路及关键基因，可为阐明不同花色形成机理提供重要的研究基础.