花生异型雄蕊转录组分析

王海霞，王铭伦，丁雨龙

(1. 南京林业大学： a． 南方现代林业协同创新中心，b． 生物与环境学院，c． 竹类研究所，江苏 南京 210037；2. 青岛农业大学，山东 青岛 266109)

花生(ArachishypogaeaLinn.)为严格的闭花受精植物，每朵花内有10枚雄蕊，其中，具4枚长椭球状花药的雄蕊为对萼长药雄蕊(antesepal stamen，Sts)，药室较大，花丝较长；具4枚圆球状花药的雄蕊为对瓣圆药雄蕊(antepetal stamen，Stp)，药室较小，花丝较短；另2枚雄蕊为退化雄蕊(staminode，Sta)，仅存花丝，未见花药。花生具有的这类雄蕊为异型雄蕊[1]。异型雄蕊现象从18世纪开始就引起了科学家们的关注[2-3]，一些研究者针对异型雄蕊的种类分布[4-7]以及花粉活性、花粉组织化学成分和结实率等[8-11]方面开展了广泛的研究，然而，异型雄蕊的遗传基础和分子调控机制目前仍未知。

花生的基因组较大、遗传基础狭窄，且其花器官的性状稳定、遗传多态性极低，因而，花生的花器官突变体很少，限制了花生的花发育分子生物学研究。雄蕊发育在被子植物的花器官进化中起关键作用，但是与其他花器官(如花瓣)相比，其分子调控机制尚处于探索阶段，对花生异型雄蕊发育的分子机制进行探讨可为雄蕊分化的分子机制研究提供思路。

作者采用RNA-seq技术对花生异型雄蕊进行转录组分析，挖掘和探究与花生雄蕊分化相关的候选基因和生物学途径，以期解释花生异型雄蕊的形成机制，为花生新品种选育提供理论基础。

1 材料和方法

1.1 材料

供试花生品种为青岛农业大学花生研究中心王铭伦教授提供的品种‘Florunner’，该品种由品种‘Early Runner’与‘Florispan’杂交获得，为普通型蔓生品种，花多且花型较大[12]。于2017年4月末将品种‘Florunner’的种子覆膜种植于青岛农业大学试验田中，并于同年6月中旬(盛花期)选取100株发育良好的植株，每天约17：00时于各植株同一位置采集发育成熟、待开放的花蕾(即花瓣已撑破萼片、微露黄色花瓣)，共采集约400枚花蕾。

将采集的花蕾置于冰盒内，在解剖镜下将雄蕊分为对萼长药雄蕊、对瓣圆药雄蕊和退化雄蕊3组，每组约800枚雄蕊(每组设置3个生物学重复)；将雄蕊分别置于液氮中速冻，-80 ℃保存，送至上海伯豪生物技术有限公司进行RNA提取及测序分析。

1.2 方法

1.2.1 总RNA的提取和检测 采用RNAqueousTMphenol-free total RNA提取试剂盒(美国Ambion公司)分别提取3组雄蕊的总RNA，获得的总RNA经Agilent Bioanalyzer 2100生物分析仪(美国Agilent Technologies公司)电泳检测合格，使用RNAClean XP Kit试剂盒(美国Beckman Coulter公司)和RNase-Free DNase Set试剂盒(德国Qiagen公司)进行纯化。

1.2.2 文库制备和测序 对纯化的总RNA进行mRNA的富集和片段化，采用SuperScript Ⅱ Reverse Transcriptase试剂盒(美国Invitrogen公司)合成cDNA第1链，然后合成cDNA第2链，随之进行末端修复，在3′末端加腺苷酸并连接测序接头，富集、纯化后构建cDNA文库；采用Illumina HiSeqTM2000高通量测序平台(美国Illumina公司)对cDNA文库进行双端测序。

1.2.3 原始数据统计和功能注释 原始测序文件经碱基识别，应用Seqtk软件进行误差过滤；将3组雄蕊样本的测序数据合并形成pool reads，然后应用Trinity软件[13-14]进行denovo拼接；将拼接得到的final unigenes序列与NCBI官方非冗余蛋白质序列数据库(Nr)和UniProt数据库进行BLASTx比对(E≤1×10-5)[15]。

应用rpstblastn程序(http：∥www.biomedcentral.com/content/supplementary/1471-2105-13-42-s1/Cloud- BioLinux-Package-Documentation/docs/rpstblastn.html)将unigenes序列与真核生物直系同源序列数据库(KOG)进行比对(E≤1×10-5)，取排名前5的比对结果进行KOG功能分析[16-17]。将unigenes序列与基因功能分类体系数据库(GO)进行BLASTx比对(E≤1×10-5)，取最佳比对结果进行GO功能分析[18]。

应用KEGG KAAS在线pathway比对分析工具(https：∥www.genome.jp/kaas-bin/kaas_main)对unigenes进行KEGG代谢途径分析[19]。

2 结果和分析

2.1 测序结果和数据组装结果分析

测序结果显示：测序样品的Q20值均大于96%，数据量约31.46 G，clean reads所占比例均在96.0%以上。说明转录组测序质量较高，可以用于denovo拼接。

对clean reads进行转录组装和拼接，得到197 193个contigs，片段总长度为240 463 900 bp，最长片段为15 810 bp，非冗余unigenes片段的平均长度为1 219 bp，GC含量为39.1%，片段长度N50为1 966 bp。统计结果显示：长度大于等于2 000 bp的unigenes有38 740个，占unigenes总数的19.6%；长度小于2 000 bp但大于等于1 000 bp 的unigenes有46 857个，占unigenes总数的23.8%；长度小于1 000 bp但大于等于200 bp 的unigenes有111 596个，占unigenes总数的56.6%。总体上长度大于等于1 000 bp的unigenes占unigenes总数的43.4%，表明组装完整性较高，可用于注释分析。

2.2 unigenes功能注释分析

2.2.1 序列相似性分析 花生异型雄蕊转录组unigenes序列与Nr数据库的比对结果(图1)显示：在花生异型雄蕊转录组中共有129 223个unigenes注释到Nr数据库，占unigenes总数的65.5%；其中，与Nr数据库中的同源序列相比，E值较低(1×10-30

图1 花生异型雄蕊转录组unigenes序列的E值分布图

在Nr数据库中对花生异型雄蕊转录组unigenes与其他植物种类转录组的注释unigenes进行同源序列匹配，结果见图2。结果显示：在花生异型雄蕊转录组中共有106 900个unigenes注释到已探明的同源序列中，占注释unigenes总数的82.7%。其中，与大豆〔Glycinemax(Linn.) Merr.〕的32 466个unigenes序列相似，占注释unigenes总数的25.1%，同源性最高；与蒺藜苜蓿(MedicagotruncatulaGaertn.)的17 343个unigenes序列相似，占注释unigenes总数的13.4%；与菜豆(PhaseolusvulgarisLinn.)的12 710个unigenes序列相似，占注释unigenes总数的9.8%；与小豆〔Vignaangularisvar.angularis(Willd.) Ohwi et H. Ohashi〕的9 804个unigenes序列相似，占注释unigenes总数的7.6%；与赤豆〔Phaseolusangularis(Willd.) Ohwi et H. Ohashi〕的3 958个unigenes序列相似，占注释unigenes总数的3.1%；与花生(Nr数据库)的3 771个unigenes序列相似，占注释unigenes总数的2.9%；与百脉根(LotuscorniculatusLinn.)的3 273个unigenes序列相似，占注释unigenes总数的2.5%；与葡萄(VitisviniferaLinn.)、桃(AmygdaluspersicaLinn.)和可可(TheobromacacaoLinn.)等13种植物的unigenes序列同源性均较低，在已注释unigenes总数中所占比例均在2.0%以下。

图2 花生异型雄蕊转录组unigenes与其他植物转录组unigenes的同源序列数量分布图

2.2.2 KOG功能分析 在KOG数据库中对花生异型雄蕊转录组unigenes序列进行功能比对，结果见图3。结果显示：在花生异型雄蕊转录组unigenes序列中，共有133 387个unigenes注释到25个KOG功能类别上，其中有8 776个unigenes注释到未知功能(function unknown)，占注释unigenes总数的6.6%，这些unigenes的功能仍不明确。

图3 花生异型雄蕊转录组unigenes的KOG功能分类

功能明确的unigenes均与配子体发育功能有关。其中，有17 135个unigenes注释到信号转导机制(signal transduction mechanism)功能，占注释unigenes总数的12.8%；有16 407个unigenes注释到一般功能预测(general function prediction only)，占注释unigenes总数的12.3%；有9 893个unigenes注释到翻译后修饰、蛋白质折叠和分子伴侣(posttranslational modification，protein turnover and chaperone)功能，占注释unigenes总数的7.4%；有8 870个unigenes注释到转录(transcription)功能，占注释unigenes总数的6.6%；有5 138个unigenes注释到细胞周期调控、细胞分裂和染色体分离(cell cycle control，cell division and chromosome partitioning)功能，占注释unigenes总数的3.9%；有2 710个unigenes注释到染色体结构和活力(chromatin structure and dynamic)功能，占注释unigenes总数的2.0%。

2.2.3 GO功能分析 在GO数据库中对花生异型雄蕊转录组unigenes序列进行功能比对，结果见表1。结果显示：在花生异型雄蕊转录组unigenes序列中，注释到生物过程、细胞组分和分子功能3大类55个亚类的unigenes共有48 497个。

在生物过程的23个亚类中，分别有31 817、25 479、16 828、6 327、6 192和5 866个unigenes注释到代谢过程(metabolic process)、细胞过程(cellular process)、单一生物过程(single-organism process)、生物调节(biological regulation)、应激反应(response to stimulus)和生物过程调节(regulation of biological process)功能，数量较多，各占注释unigenes总数的65.6%、52.5%、34.7%、13.0%、12.8%和12.1%。其中，在核酸代谢过程(GO： 0090304，8334个)、磷酸化(GO： 0016310，4428个)和转录调节(GO：0006355，2322个)等与代谢过程相关的GO功能类别中unigenes高度富集；而注释到节律过程(rhythmic process)、运动(locomotion)、排毒(detoxification)、生物附着(biological adhesion)和行为(behavior)功能的unigenes数量则较少。

在细胞组分的18个亚类中，分别有13 334、13 334、9 427和8 905个unigenes注释到细胞(cell)、细胞成分(cell part)、细胞器(organelle)和细胞膜(cell membrane)功能，数量较多，各占注释unigenes总数的27.5%、27.5%、19.4%和18.4%。其中，在与核(GO：0005634，4 069个)和核糖体(GO：0005840，639个)等细胞器相关的GO功能类别中unigenes高度富集；而注释到类核(nucleoid)、病毒体部分(virion part)、病毒(virion)、膜外区域部分(extracellular region part)、胞外基质(extracellular matrix)和胞外基质组分(extracellular matrix component)功能的unigenes数量则较少。

表1 花生异型雄蕊转录组unigenes的GO功能分类

在分子功能的14个亚类中，分别有28 066和26 091个unigenes注释到结合(binding)和催化活性(catalytic activity)功能，数量较多，各占注释unigenes总数的57.9%和53.8%。其中，在ATP结合(GO：0005524，7 229个)、转录因子(GO：0003700，1 161个)和钙离子结合(GO：0005509，728个)等与结合相关的GO功能类别中unigenes高度富集；而注释到营养储存活性(nutrient reservoir activity)、转录调控因子活性(translation regulator activity)和金属伴侣蛋白活性(metallochaperone activity)功能的unigenes数量则较少。

2.2.4 KEGG代谢通路分析 以KEGG代谢库为参考，对花生异型雄蕊转录组unigenes的代谢通路进行统计和分类，结果见图4。结果显示：在花生异型雄蕊转录组中共有30 428个unigenes注释到2 708个酶功能，映射了354个代谢通路，可分为遗传信息处理(genetic information processing)、细胞过程(cellular process)、生物系统(organismal system)、环境信息处理(environmental information processing)和代谢 (metabolism) 5大类32个亚类，其中与遗传信息处理、细胞过程和环境信息处理相关的亚类较少，分别仅有4、4和3个。

在遗传信息处理代谢通路中，有3 307个unigenes注释到翻译(transcription)亚类，有3 015个unigenes注释到折叠、分类和降解(folding，sorting and degradation)亚类。在细胞过程代谢通路中，有2 713个unigenes注释到运输和分解代谢(transport and catabolism)亚类。在环境信息处理代谢通路中，有4 911个unigenes注释到信号转导(signal transduction)亚类，其中，有1 055个unigenes注释到植物激素信号转导通路(ko04075)，有285个unigenes注释到钙信号通路(ko01100)。

： 遗传信息处理Genetic information processing； ： 细胞过程Cellular process； ： 生物系统Organismal system； ： 环境信息处理Environmental information processing； ： 代谢Metabolism.

共有5 475个unigenes注释到与生物系统相关的代谢通路中，占注释unigenes总数的18.0%，分别注释到9个亚类。其中，有1 913个unigenes注释到内分泌系统(endocrine system)亚类，有1 453个unigenes注释到免疫系统(immune system)亚类。

注释到与代谢相关通路的unigenes数量最多，有14 231个，占注释unigenes总数的46.8%，分别注释到12个亚类。其中，有6 529个unigenes注释到全球和总览图(global and overview maps)亚类，有4 635个unigenes注释到糖类代谢(carbohydrate metabolism)亚类，有2 802个unigenes注释到氨基酸代谢(amino acid metabolism)亚类，有2 538个unigenes注释到脂质代谢(lipid metabolism)亚类，有2 337个unigenes注释到能量代谢(energy metabolism)亚类，有1 278个unigenes注释到核苷酸代谢(nucleotide metabolism)亚类。

3 讨论和结论

花生隶属于豆科(Fabaceae)落花生属(ArachisLinn.)，从同源序列的匹配结果看，其异型雄蕊转录组unigenes与同科的大豆和蒺藜苜蓿的同源匹配率较高，与大豆的32 466个unigenes序列相似，与蒺藜苜蓿的17 343个unigenes序列相似，这一结果佐证了这些种类亲缘关系。而供试花生品种的异型雄蕊转录组中仅有3 771个unigenes与Nr数据库中花生的同源unigenes相匹配，一方面可能是由于本实验是以花生异型雄蕊为实验材料，而Nr数据库中的花生基因组注释信息是以果针为研究材料有关；另一方面，花生栽培种是由2个二倍体野生种A.duranensis和A.ipa⊇nsis自然杂交后经过染色体加倍形成的异源四倍体[20-21]，导致花生栽培品种的遗传关系复杂。

有花植物的雄蕊和花粉发育是一个复杂的过程，与花被器官相比，生殖器官的特有基因更多[22]。对花生异型雄蕊转录组测序获得数据量约31.46 G，其中，有129 223个unigenes在Nr数据库比对出同源序列信息，且许多unigenes功能与雄蕊产生雄配子的功能有关。许多转录因子决定着花器官中萼片、花瓣、雄蕊和心皮[23]的发育时序与数量，进而参与控制被子植物花发育。在花生异型雄蕊转录组中，有8 870个unigenes注释到转录功能，其中5 360个unigenes具有转录因子功能，说明转录因子对植物基因表达的调控有至关重要的作用[24]。值得注意的是，有6.6%的unigenes基因功能仍不明确，这与花生的基因组大、重复序列比例高等原因有关。

雄蕊和花粉的发育涉及一系列精密的细胞质和核基因相互作用[25]；Ca2+可作为细胞信号转导的信使，通过钙调蛋白(CaM)将胞外信号转换成胞内的生理生化反应，包括细胞分裂、分化和凋亡，对花粉发育和萌发起着非常重要的作用[26]。在GO数据库中对花生异型雄蕊转录组unigenes序列进行功能比对，结果显示其unigenes在与核酸代谢过程、磷酸化、转录调节、核、核糖体、ATP结合、转录因子和钙离子结合相关的功能类别中高度富集；而KEGG代谢通路分析结果表明：在花生异型雄蕊转录组中注释到代谢相关通路的unigenes数量最多，并以注释到糖类代谢通路、氨基酸代谢通路、脂质代谢通路、能量代谢通路和核苷酸代谢亚类的unigenes数量较多，说明有大量涉及孢粉素合成、氨基酸合成、核糖体组装、内质网蛋白质加工和脂质转移等功能的unigenes被注释；此外，与钙调蛋白相关的unigenes也被注释到无机离子运输与代谢以及信号转导机制这2个类别上。由此可知，在花生异型雄蕊转录组中，一部分基因可参与不同的代谢途径，同时花药发育的精密调控需要多个代谢途径参与，其中，糖类代谢不仅为花药发育提供能量，还可作为信号物质影响花药和花粉的发育[27]；而在花粉发育过程中通常会观察到高呼吸率和高能量需求[28]，主要代谢途径的功能障碍将对花粉粒发育产生不利影响[29]，雄性不育系中与能量供应相关的ATP含量显著降低[30-31]。因而，对于花生异型雄蕊，特别是退化雄蕊与可育雄蕊进行差异表达基因功能分析，可为花生雄蕊发育的分子调控提供更多启示。

花生异型雄蕊的发育调控是一个复杂的生物过程，应结合花生基因组测序成果开展大量的基因组分析，以期挖掘更多有价值的遗传信息，从而更全面地了解花生异型雄蕊分化的分子调控机制。