基于高通量测序的半夏变异珠芽转录组分析

马琛，汪雷，叶德，徐涛

(浙江理工大学 生命科学学院,浙江 杭州 310018)

半夏Pinelliaternata(Thunb.)Breit，多年生草本植物，天南星科半夏属，在我国，其块茎是高利用的中药之一；据药典记录，其具有燥湿化痰，降逆止呕，消痞散结功效，常用于痰多咳喘，风痰眩晕，痰厥头痛，呕吐反胃，梅核气等。作为重要的传统中药，半夏常以干燥块茎入药，最早记载于《神农本草经》，列为下品，其性温、味辛、略有毒，归脾、胃、肺；具有外用止痛消肿及抗肿瘤潜力[1-4]。珠芽繁殖是半夏繁殖的重要繁殖方式，其具有快速繁殖，产量高等优点，在自然界和农业生产中能增加半夏的存活率。因此，珠芽品质的优劣是决定半夏质量与产量的重要因子，也极大影响了半夏的种植。由于野生的半夏资源已越来越稀少，而半夏的需求却日渐增大，因此，为了培育优质高产的半夏品种，利用基因工程手段就成为缓解半夏资源紧缺的解决手段之一。目前的半夏研究主要集中在优质品种的筛选和高产量、药理药效及临床等方面，对半夏珠芽功能基因的研究较少、进展缓慢。

筛选差异表达基因(differential expressed gene, DEG)可针对性地分析功能基因，利用这些基因对物种进行定向遗传育种改良，提高抗胁迫能力及其分子机制的深入研究[5-6]。基因差异表达是指植株在特定环境下，发生生理应激反应引起功能蛋白的差异表达，导致生长表型出现差异的现象。转录组(transcriptome)是指一定状态下，细胞内所有转录产物的集合，转录组测序能够全面快速地获取某一物种特定器官或组织在某一状态下的几乎所有转录体，并反映其表达水平[7-9]。利用转录组分析能够批量地筛选发掘无参考基因组物种的功能基因，以及探索在不同样品中基因的表达差异。本研究借助高通量测序技术[10]，比较分析半夏叶基部长珠芽和不长珠芽的DEG，并筛选出影响半夏珠芽激素生物合成和信号传导的DEG，为深入研究半夏珠芽相关功能基因的克隆及功能分析提供基础。

1 材料与方法

1.1 材料

三叶半夏种子取自甘肃天水，种于浙江理工大学6号楼317实验室内；选取发育正常、长势优良且无病虫害侵袭的植株进行采样，用解剖刀分别刮取半夏叶基部相同部位长珠芽和不长珠芽的组织，用蒸馏水将其清洗干净，立即用液氮冷冻，放入-80 ℃冰箱保存备用。样品与分组信息：半夏叶基部长珠芽为样品BY，不长珠芽的为样品ZC。本实验利用AxyPrep总RNA小量制备试剂盒，按步骤提取总RNA，并进行质量检测。经质检合格后，使用连接有Oligo(dT)的磁珠富集真核生物mRNA。片段化试剂(fragmentation buffer)将mRNA随机打断成短片段，以片段化的mRNA为模板，利用随机引物(random hexamers)合成一链cDNA，随后合成二链cDNA。AMPure XP beads纯化双链产物，利用T4 DNA聚合酶和Klenow DNA聚合酶将DNA的黏性末端修复为平末端，3’末端加碱基A并加接头，AMPureXP beads进行片段选择，最后进行PCR扩增获得最终测序文库。文库质检合格后用Illumina HiseqTM 4000进行测序，开展转录组研究，共获得单基因(unigene)28 395条。本研究在此基础上分析BY和ZC样品的DEGs。

1.2 DEG的筛选条件

通过RPKM(reads per kilo base of exon model per million mapped reads)法来计算样品BY和样品ZC的unigene表达的丰度值，表示在每百万测序碱基中每千个转录子测序碱基中所包含的测序reads数。计算公式：

RPKM=total exon reads÷mapped reads (million)÷exonlenth(KB)。

本研究根据公式计算变异珠芽和正常珠芽之间表达差异的倍数，用log2fold_change表示样本基因表达的相对差异倍数，计算方法为：log2(Sample_BY FPKM/Sample_ZC FPKM)，即样本BY除以ZC的基因表达水平。将RPKM的log2倍数绝对值>1且FDR(false discovery rate，错误发现率)≤0.05的unigene筛选出，即为DEG。

1.3 DEG的功能注释

将筛选到的unigene序列分别映射到公共数据库Swiss-Prot、NR、Pfam、KEGG、COG、GO中，通过BLASTX得到与指定unigene具有序列最高相似性(取阈值e<=1e-10)的蛋白(表1)，进一步得到相应unigene所注释的基因及其蛋白功能信息。比对到公共数据库NR的结果显示，高达71.22%的unigenes序列找到了相似序列。而近缘物种在相似性序列的匹配中，海枣所占比例最高(44.2%)，其他依次是葡萄(37.5%)、可可树(8.2%)、无油樟(3.2%)、大桉(2.6%)、桐油树(2.2%)和其他物种(2.1%)。

表1 BLAST注释结果统计情况

将筛选到的半夏unigenes比对到数据库COG里，结果显示，共有17 130条unigenes(占总数60.33%)被注释到25种COG分类中(表2)。被注释的unigenes涵盖大部分生命活动，涉及KOG功能类别也比较全面；其中，最大的分类是“基因功能预测”，包含2 063条(12.04%)unigenes；其次是“信号传导机制”，包含了1 788条(10.44%)unigenes；而“细胞运动”是最少的，仅4条(0.02%)unigenes；在其他类别中，各自的基因表达丰度都不相同。

表2 COG功能分类

2 结果与分析

2.1 DEG筛选

De novo测序中最重要的过程就是基因差异表达分析，通过差异表达分析可以了解不同样本间基因的差异表达情况。根据条件|Log2(样品BY RPKM/样品ZC RPKM)|≥1、且FDR≤0.05，共筛选到13 656条DEG，其中包含6 139条上调表达的DEG和7 517条下调表达的DEG。

2.2 DEG 功能分析

将筛选所得的13 656条DEGs导入数据库GO中，深入探索其功能定义。结果表明，在GO数据库中获得注释的DEG共有4 595条，占DEG总量的33.65%。其中分别将2 324、592和1 679条DEG归类至三大主要功能类别，即生物学途径(biological process)、细胞组件(cellular component)和分子生物学功能(molecular function)。将被注释的DEG功能定义可进一步分在3个主要功能类别的50个亚类中(图1)，在“生物学途径”类别中，“依赖DNA的转录调节(regulation of transcription, DNA-dependent)”是DEG归入最多的亚类；最少的则是“多细胞有机体的发育(multicellular organismal development)”。在“细胞组件”类别中，DEG归入最多的亚类是“膜的完整性(integral to membrane)”；最少的是“叶绿体类囊体膜(chloroplast thylakoid membrane)”。在“分子生物学功能”类别中，DEG归入最多的亚类是“ATP结合(ATP binding)”；最少的是“GTP结合(GTP binding)”。

图1 DEG的GO分类

2.3 DEG KEGG富集分析

生物体内，生物学功能的行使依赖于不同基因的相互协调，通过分析代谢途径可以确定DEG的主要代谢途径及信号通路等。获得的13 656条DEG比对到KEGG代谢途径数据库，获得功能注释的DEG共3 551条，占DEG总数的26.00%，包含在239条KEGG代谢途径中。在KEGG途径的功能定义中，主要包括产物代谢与合成及细胞作用等，表3中是富集DEG最多的10条KEGG代谢途径。

2.4 影响激素合成和信号传导的DEGs分析

将筛选到的DEGs与Nr、GO、KEGG等蛋白公共数据库分别比对，得到DEGs的功能注释，并筛选到参与半夏珠芽激素生物合成及信号传导的关键酶基因17个(表4)。通过分析该基因在半夏变异植株和正常植株中的RPKM值显示，表达上调基因为6个，即在BY中的表达量高于在ZC中的表达量；表达下调基因为11个，即在ZC中的表达量高于BY中的表达量。

表3 半夏DEG富集最多的10条KEGG代谢途径

3 讨论

植物生长中的发育差异及生理结构的变化甚至病变，细胞功能多样性和形态差异的产生等决定性因子是基因的差异性表达[11]。因此，分析基因差异性表达，可有效挖掘影响植物性状的相关功能基因，促进植物遗传育种的改良。20世纪90年代起至今，差减杂交、mRNA差异显示PCR、抑制性消减杂交、基因表达的系列分析(SAGE)、表达序列标签(EST)、cDNA限制性片段长度多态性分析、基因芯片和荧光定量PCR等经典基因表达差异分析方法已经形成[12]。但这些方法都各自存在一定的不足，如差减杂交经过层析，所回收的cDNA量很低，严重影响实验的进一步进行，且只能检测2组mRNA之间的表达差异；mRNA差异显示逆转录PCR需反复操作，进行多次电泳及PCR。本研究通过高通量测序技术挖掘获得相关功能基因，并进行基因新功能的探讨；采用自动化程序，同时分析多组mRNA，检测出存在差异的单个核苷酸。分析DEG的同时，还能获得SSR[13]信息及组装出的完整基因序列。本研究从半夏珠芽转录组数据的28 395条unigenes中，高效批量地筛选出DEGs 13 656条；通过GO归类，该DEGs被包含在50个功能亚类中，显著集中在半夏生长发育过程中的性状表现及次生代谢等。分析KEGG途径显示，有3 551条半夏珠芽差异表达基因富集在生理代谢过程和次生代谢产物合成途径中，占DEGs总量的26.00%。本研究获得的DEGs数据丰富了半夏珠芽的基因资源，对研究半夏珠芽的基因调控机制意义重大。

表4 半夏珠芽激素合成与信号转导相关DEG

植物生长发育过程中，参与调控功能的内源性激素非常重要，其调节生物合成及信号传导的酶基因也成为研究热点。基于半夏珠芽转录组DEGs分析，本研究共筛选获得参与半夏珠芽内源性激素生物合成和信号传导相关过程的功能基因17个，这些基因在半夏叶基部长珠芽和不长珠芽中呈现出明显的表达差异，因此，半夏珠芽生长发育过程与这些基因及其调节的代谢产物密切相关。

生长素即吲哚乙酸，是最早发现的促进植物生长的激素，但对茎、芽、根生长的促进作用因浓度而异，调控细胞分裂和分化及顶端优势等一系列生理过程[14]。本研究表明，在ZC中，半夏生长素调控基因低表达，而在BY中高表达，推断在BY中生长素积累高于ZC，这可能是导致半夏叶基部长珠芽的因素之一。赤霉素最显著的功能是加速细胞的伸长，促进细胞分裂和扩大，且不引起细胞壁的酸化，除此之外，赤霉素还有着抑制成熟，侧芽休眠，衰老，块茎形成的生理作用[15]。本研究从半夏珠芽转录组中筛选到的GA 2-oxidase，是赤霉素氧化酶的编码基因，该基因主要表达于ZC 中，表明半夏珠芽赤霉素合成具有特异性。细胞分裂素(cytokinin kinetin)是一类调节植物细胞生长和发育的植物激素，与植物生长素有协同作用；从玉米或其他植物中分离得到，一般在植物根部产生。其有两种明显的生理作用：一是促进细胞分裂和调控其分化，二是延缓蛋白质和叶绿素的降解，延迟衰老。本研究从半夏转录组中鉴定到的乙酰辅酶A羧化酶/生物素羧化酶[16]，在半夏BY中明显高于半夏ZC中。脱落酸是一种抑制生长的植物激素，广泛分布于高等植物，除促使叶子脱落外，还可使芽进入休眠状态，抑制细胞的分裂；因此可促进叶和果实的衰老和脱落，增强植物抗逆性[17-18]。本研究筛选到9-顺式-环氧类胡萝卜素双加氧酶[19]，是一个ABA生物合成的限速关键酶基因，密切影响ABA的积累，且渗透胁迫诱导合成的关键性调节位点即为该基因催化的步骤。相比BY中，9-顺式-环氧类胡萝卜素双加氧酶在ZC表达更高，显示ABA合成在ZC中更为活跃。乙烯在种子萌发到成熟衰老和植物花、果实的成熟及衰老等一系列生理过程，起重要的调节作用，间接影响植物的性状差异及品质优劣[20]。研究挖掘到调控乙烯生物合成的重要关键酶基因——1-氨基环丙烷-1-羧酸合酶[21]，是催化乙烯生物合成的限速酶，将S-腺苷甲硫氨酸转变成1-氨基环丙烷基-1-羧酸和甲基硫腺苷。其在半夏BY中表达下调，在半夏ZC中表达高于在BY中，说明乙烯主要在ZC中合成明显。油菜素甾醇被称为第六类植物激素，在植物中含量低微，仅有其他五大激素任意一种的千分之一，但其具有的高生理活性、广谱、无毒及独特作用机理，对植物的增产和抗逆性具有重大意义[22]。并且生长素和油菜素甾醇通过共同调控下游基因的转录、二者信号组分间的互作、合成代谢的调控以及调控生长素极性运输等，在多层次相互作用从而对植物生长发育实现精确地调控[23]。植物中油菜素甾醇的信号转导及其与脱落酸信号的相互调节对半夏的生长发育也具有重要影响。茉莉酸在种子的萌发、果实的成熟、根的生长、球茎的形成以及孕育花粉等生长发育相关的生理过程，同时在系统防御生物过程起着关键的调控作用[24]。本研究中茉莉酸生物合成在半夏BY中明显高于半夏ZC，可以推断，半夏BY植株在生长过程中有更多的优势，对于半夏产量和质量的提高都大有裨益。

本研究利用高通量测序技术对三叶半夏叶基部长珠芽和不长珠芽进行转录组和基因组分析,丰富了半夏基因组的信息,并挖掘到大量差异表达基因，为后续半夏功能基因的分析及克隆提供了丰富的序列数据基础。本研究将筛选到的半夏unigenes比对到数据库COG里，结果显示：共有17 130条unigenes被注释到25种KOG分类中，其中721条基因功能未知。希望未来更多物种可以进行转录组研究，在丰富基因公共数据库的同时，发掘已知基因的新功能和新的功能基因。本研究筛选得到13 656条显著性差异表达基因，需要对其进行深入研究，发掘出半夏生长发育中控制性状的关键基因，探索生物学调控机制，为半夏优质高产品种的改良提供科学依据，促进中药半夏的深度开发。