基于转录组测序对紫花苜蓿细胞质雄性不育系相关代谢通路的鉴定

陈晶晶，王英哲，郭 强，李海静，徐 博

(1吉林农业大学 动物科学技术学院，吉林 长春 130118；2吉林省农业科学院 农业生物技术研究所，吉林 长春 130118)

细胞质雄性不育(CMS)是高等植物常见的生物学特性，是指CMS植物不能产生正常的花药、花粉或雄配子[1-3]。由于细胞质雄性不育系在杂交种子生产中具有重要价值，前人已经鉴定出许多作物品种的细胞质雄性不育系，如棉花(GossypiumhirsutumL.)、大豆(GlycinemaxL.)、胡萝卜(DaucuscarotaL.)、玉米(ZeamaysL.)、洋葱(AlliumcepaL.)、甘蓝型油菜(BrassicanapusL.)、水稻(OryzasativaL.)、向日葵(HelianthusannuusL.)和小麦(TriticumaestivumL.)[4-5]。紫花苜蓿(MedicagosativeL. )是世界上最重要的豆科牧草之一，因其具有优良的营养品质以及较高的饲草产量而被广泛种植[6-7]。吴永敷[8]在1978年从草原1号中发现了6株紫花苜蓿雄性不育系,于洪柱等[9]在2008年发现了紫花苜蓿细胞质雄性不育系材料-MS-GN。

转录组测序是近年来新开发的一种高性能分析方法，运用该技术可以在一个试验中对数千个基因进行比较和分析，有助于发现新的基因和转录因子[10-13]。因此，众多学者采用转录组测序技术对CMS作物的不育机理进行研究。Suzuki等[14]对棉花细胞质雄性不育系和恢复系进行了转录组分析，结果表明与CMS相关的基因主要参与了细胞壁的扩增。Li等[15]采用转录组测序分析大豆细胞质雄性不育系NJCMS1A及其保持系NJCMS1B，认为不育系NJCMS1A的雄性不育可能与一些关键差异表达基因功能紊乱和代谢途径异常相关，如碳水化合物能量代谢、转录因子、花粉发育、活性氧(ROS)清除体系、细胞信号转导和细胞程序性死亡(PCD) 等。然而，截至目前，尚未发现紫花苜蓿细胞质雄性不育系转录组测序的相关研究。

在前人研究的基础上，笔者发现紫花苜蓿细胞质雄性不育系的花粉发育不良及育性丧失可能是由基因突变引起的，但其分子机理尚不清楚。因此，本研究以紫花苜蓿细胞质雄性不育系及其保持系的花药为材料，进行转录组测序，得到与细胞质雄性不育相关的通路及基因，进一步筛选差异基因，为后续研究其不育机理及基因表达分析奠定基础。

1 材料与方法

1.1 试验材料

紫花苜蓿细胞质雄性不育系MSGN1A(以公农-3号为供体母本，保持系MSGN1B作为轮回亲本，回交四代所得)及其保持系MSGN1B，均由吉林省农业科学院提供。2016年4月，将MSGN1A和MSGN1B种植于吉林省长春市吉林农业大学草业科学试验田，选取现蕾期24 h(阶段Ⅰ)、48 h(阶段Ⅱ)和60 h(阶段Ⅲ)的花蕾，去除花瓣及托叶，保留花药部分，每个时期各取样品50 mg，等量混合保存于液氮罐中，备用。试验设置3个生物学重复。

1.2 试验方法

1.2.1 MSGN1A和MSGN1B花药总RNA的提取及质量控制 参照Garg等[16]的方法提取MSGN1A和MSGN1B花药总 RNA，再用紫外分光光度计分别在230，260和280 nm处测定吸光值(OD230、OD260和OD280)，要求OD260/OD280和OD260/OD230均≥1.8。运用Agilent 2100生物分析仪检测总RNA的完整性，总RNA完整性值在8.0～10.0的样品可用于后续试验。

1.2.2 cDNA文库的制备与检测 cDNA文库的制备与检测由北京百迈客生物技术公司完成。总RNA经过DNase消化、mRNA富集、mRNA打断后，反转录合成cDNA第一链。以cDNA第一链为模板，反转录合成cDNA第二链。对双链cDNA进行纯化，然后进行末端修复，最后在3′末端加上A碱基，并连接测序接头(3′端接头为AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC；5′端接头为AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT)。连接产物经胶回收后进行PCR反应及产物回收，cDNA文库构建完成。分别使用Qubit 2.0和Agilent 2100对cDNA文库进行检测，使用q-PCR方法对文库的有效浓度进行准确定量，以保证文库质量。

1.2.3 转录组测序及测序数据组装 转录组测序及测序数据组装均由北京百迈客生物技术公司完成。应用Illumina Hiseq 2000平台进行转录组文库测序，测序条件为：PE150，双端序列150 bp。对原始序列进行数据过滤，去除其中的接头序列及低质量序列，获得高质量的去冗余序列。采用Trinity软件进行序列组装，获得紫花苜蓿的功能基因。

1.2.4 功能基因差异表达分析及功能注释 根据功能基因在不育系MSGN1A和保持系MSGN1B样品中的表达量进行差异表达分析(以FDR<0.01和FC>2作为筛选MSGN1A和MSGN1B之间差异表达基因(DEGs)的阈值)及差异表达基因的功能注释。功能注释由北京百迈客生物技术公司进行，功能注释的数据库为GO (Gene Ontology，http://www.geneontology.org/)、KEGG (Kyoto Encyclopedia of Genes and Genomes，http://www.genome.jp/kegg/)、COG (Clusters of Orthologous Groups of proteins/orthologous groups of genes)。

2 结果与分析

2.1 紫花苜蓿MSGN1A和MSGN1B总RNA质量分析

对提取的紫花苜蓿MSGN1A和MSGN1B花药总RNA的质量进行检测，检测结果如表1所示。由表1可知，紫花苜蓿不育系MSGN1A和保持系MSGN1B花药总RNA完整性值均大于8.0，且OD260/OD280均大于1.9，OD260/OD230值均大于2.0，总RNA质量浓度分别为158和163 ng/μL，说明样品总RNA的完整性与纯度合格达到测序要求。

表1 紫花苜蓿MSGN1A和MSGN1B总RNA质量检测结果Table 1 Result of quality detection of total RNA in alfalfa MSGN1A and MSGN1B

2.2 紫花苜蓿转录组测序及测序数据分析

利用Illumina Hiseq 2000平台对紫花苜蓿细胞质雄性不育系MSGN1A及其保持系MSGN1B的花药进行转录组测序分析，结果如表2所示。由表2可知，紫花苜蓿转录组测序共获得去冗余序列54.43 G，其中不育系去冗余序列为27.43 G，G+C含量为41.98%，质量值(Q)大于等于30的碱基占91.21%；保持系去冗余序列为27.00G，G+C含量为41.90%，质量值大于等于30的碱基占90.99%。

采用Trinity软件对去冗余后的高质量序列进行组装，获得172 483个转录本(Transcripts)和95 679个功能基因(Unigenes)，不同长度转录本、功能基因的条数及其比例见表3。转录本和功能基因的总长度分别达到160 772 822和73 481 347 bp；转录本平均长度为932.11 bp,N50为1 437 bp；功能基因平均长度为768 bp，N50为1 209 bp。

表2 紫花苜蓿转录组测序数据统计结果Table 2 Results of illumina transcription sequencing for alfalfa

表3 紫花苜蓿转录本和功能基因长度分布和数目统计Table 3 Number and length range of transcripts and unigenes of alfalfa genome

2.3 紫花苜蓿功能基因差异表达分析及其功能注释

2.3.1 DEGs分析 与保持系MSGN1B相比，紫花苜蓿不育系MSGN1A有187个DEGs，其中169个表达下调，18个表达上调。

2.3.2 GO功能注释 在序列同源性(≥80%)的基础上，对187个DEGs进行GO功能注释，结果(图1)有3个上调和84个下调的DEGs被标记为30个功能类别，其中包括14个生物学过程，8个细胞组分和8个分子功能。生物学过程功能类别中，代谢过程是主要的功能类别，其次是细胞过程、信号生物过程、细胞成分组织或生物合成、生物调节等。在细胞组分功能类别中，细胞和细胞分离是主要的功能类别，其次是生物膜、膜分离等。分子功能类别中，催化活性是最主要的功能类别，紧随其后的是合成、转运活性、分子功能调控等功能类别。这证明了催化活性、代谢过程、细胞过程、合成、信号生物过程、生物膜和膜分离等功能可能与紫花苜蓿细胞质雄性不育相关。

1.生物调节;2.细胞成分组织或生物合成；3.细胞过程；4.发育过程；5.生长；6.定位；7.代谢过程；8.多细胞生物过程；9.多生物过程；10.繁殖；11.繁殖过程；12.刺激反应；13.信号；14.信号生物过程；15.细胞；16. 细胞分离；17.胞外区；18.高分子配合物；19.生物膜；20.膜分离；21.细胞器；22.细胞器分离；23.抗氧化活性；24.合成；25.催化活性；26.电子传递体活性；27.分子功能调控；28.核酸合成转录因子活性；29.结构分子活性；30.转运活性1.Biological regulation;2.Cellular component organization or biogenesis;3.Cellular process;4.Developmental process;5.Growth;6.Localization;7.Metabolic process;8.Multicellular organismal process;9.Multi-organism process;10.Reproduction;11.Reproductive process;12.Response to stimulus;13.Signaling;14.Single-organism process;15.Cell;16.Cell part;17.Extracellular region;18.Macromolecular complex;19.Membrane;20.Membrane part;21.Organelle;22.Organelle part;23.Antioxidant activity;24.Binding;25.Catalytic activity;26.Electron carrier activity;27.Molecularfunction regulator;28.Nucleic acid binding transcription factor activity;29.Structural molecule activity;30.Transporter activity图1 紫花苜蓿细胞质雄性不育系MSGN1A差异表达基因的GO功能注释Fig.1 Histogram showing GO function annotation of DEGs of alfalfa CMS line MSGN1A

2.3.3 KEGG通路分析 为了识别DEGs参与的代谢途径，通过KEGG通路数据库对187个DEGs进行分析，结果如图2所示。

1.类黄酮生物合成；2.苯丙素的生物合成；3.缬氨酸，亮氨酸和异亮氨酸的生物合成；4.泛素介导的蛋白水解作用；5.戊糖、葡萄糖醛酸转换；6.泛酸酯和CoA生物合成；7.内质网蛋白加工；8.淀粉和蔗糖代谢；9.糖胺聚糖降解；10.缬氨酸、亮氨酸和异亮氨酸降解；11.自噬调节；12.2-氧代羧酸代谢；13.氨基糖和核苷酸糖代谢；14.植物病原菌互作；15.氨基酸的生物合成；16.植物激素信号转导；17.RNA转运；18.半胱氨酸和蛋氨酸代谢；19.氨酰生物合成；20.囊泡运输中的SNARE蛋白相互作用;21.鞘脂类代谢;22.核糖体;23.氧化磷酸化; 24.色氨基酸代谢1.Flavonoid biosynthesis;2.Phenylpropanoid biosynthesis;3.Valine,leucine and isoleucine biosynthesis;4.Ubiquitin mediated proteolysis;5.Pentose and glucuronate interconversions;6.Pantothenate and CoA biosynthesis;7.Protein processing in endoplasmic reticulum;8.Starch and sucrose metabolism;9.Glycosaminoglycan degradation;10.Valine,leucine and isoleucine degradation;11.Regulation of autophagy;12.2-Oxocarboxylic acid metabolism;13.Amino sugar and nucleotide sugar metabolism;14.Plant-pathogen interaction;15.Biosynthesis of amino acids;16.Plant hormone signal transduction;17.RNA transport;18.Cysteine and methionine metabolism;19.Aminoacyl-tRNA biosynthesis;20.SNARE interactions in vesicular transport;21.Sphingolipid metabolism;22.Ribosome;23.Oxidative phosphorylation;24.Cyan amino acid metabolism图2 紫花苜蓿细胞质雄性不育系MSGN1A差异表达基因的KEGG功能分类Fig.2 Histogram showing KEGG function classification of DEGs of alfalfa CMS line MSGN1A

由图2可知，在187个DEGs中有4个上调和33个下调DEGs被注释到24个KEGG通路，其中淀粉和蔗糖代谢是最具代表性的途径(KEGG通路编号为ko00500，共14个DEGs，下同)，其次是戊糖、葡萄糖醛酸转换(ko00040，12),氧化磷酸化(ko00190，3)，植物激素信号转导(ko04626，2)，植物-病原菌互作(ko04626，2)，RNA转运(ko04130，2)，苯丙素的生物合成(ko00940，2)。少数几个基因参加了自噬调节(ko04140，1)，半胱氨酸和蛋氨酸代谢(ko00270，1)，泛酸酯和CoA生物合成(ko00770,1)等。试验结果表明，与紫花苜蓿细胞质雄性不育相关的主要代谢通路包括碳水化合物代谢通路，戊糖、葡萄糖醛酸转换通路，氧化磷酸化通路等。

2.3.4 COG功能分类 对187个DEGs进行COG功能分类，共有50个DEGs(原始结果为109个DEGs，但其中包括59个重复DEGs，剔除重复后实际为50个DEGs)注释到16个COG功能类别(图3)，其中3个为上调DEGs，47个为下调DEGs。在COG功能分类中，一般功能预测(有23个DEGs，占比为21.1%(23/109)，下文同)是最具代表性的分类；其次是信号传导机制(18，16.51%)，转录(17，15.6%)，复制、重组和修复(17，15.6%)；再次是碳水化合物运输和代谢(8，8.34%)，细胞壁/膜的生物合成(4，3.67%)，次级代谢产物生物合成、运输和分解代谢(4，3.6%)，无机离子的转运和代谢(4，3.67%)；其他功能分类的DEGs及其占比较低。结果表明，信号传导机制，转录、复制、重组和修复，碳水化合物运输和代谢，细胞壁/膜/包膜生物发生可能涉及到紫花苜蓿的细胞质雄性不育。

1.能量产生与转换;2.细胞周期控制，细胞分裂，染色体分裂；3.氨基酸运输和代谢；4.碳水化合物运输和代谢；5.辅酶的运输和代谢；6.翻译、核糖体结构和生物起源；7.转录；8.复制、重组和修复；9.细胞壁/膜/包膜生物发生；10.无机离子的转运和代谢；11.次生代谢物生物合成、转运和分解代谢；12.一般功能预测；13.功能未知；14.信号传导机制；15.胞内运输、分泌和囊泡运输；16.细胞骨架1.Energy production and conversion;2.Cell cycle control,cell division,chromosome partitioning;3.Amino acid transport and metabolism;4.Carbohydrate transport and metabolism;5.Coenzyme transport and metabolism;6.Translation,ribosomal structure and biogenesis;7.Transcription;8.Replication,recombination and repair;9.Cell wall/membrance/envelope biogenesis;10.Inorganic ion transport and metabolism;11.Secondary metabolites biosynthesis,transport and catabolism;12.General function prediction only;13.Function unknown;14.Signal transduction mechanisms;15.Intracellular trafficking, secretion,and vesicular transport;16.Cytoskeleton图3 紫花苜蓿细胞质雄性不育系MSGN1A差异表达基因的COG功能分类Fig.3 Histogram showing COG function classification of DEGs of alfalfa CMS line MSGN1A

3 结论与讨论

在花粉发育过程中，酶可以催化葡萄糖、淀粉和其他糖类物质的合成[17]。而糖类物质是植物体内重要的能源物质[18]，因此，糖类物质代谢与花粉发育息息相关。由KEGG通路注释结果可知，本研究筛选出了14个涉及淀粉、蔗糖代谢的DEGs和12个涉及戊糖和葡萄糖酸盐互换的DEGs，这些基因全部下调表达，表明紫花苜蓿细胞质雄性不育系MSGN1A中这17个基因(剔除重复的9个基因)的表达水平低于其保持系MSGN1B，导致不育系糖类物质合成与代谢异常，从而致使不育系花粉育性丧失。另外，从COG分析结果可知，8个DEGs与碳水化合物运输和代谢过程相关。这与Dong等[19]在大白菜CMS中发现碳水化合物代谢相关DEGs下调表达的结果类似。

Cao等[20]研究发现，在绒毡层细胞中Ca2+-ATPase的异常分布会影响正常细胞程序性死亡(PCD)，且不能正常降解绒毡层细胞，小孢子发育所需的能量及营养不能及时供应，导致花粉败育。本研究KEGG数据库中注释到，2个能量产生和代谢的DEGs下调表达， 3个DEGs与氧化磷酸化有关且下调表达，而氧化磷酸化是合成三磷酸腺苷(ATP)的重要途径之一，也是能量转换的重要方式[21]，因此认为MSGN1A中与能量代谢和氧化磷酸化相关的基因表达异常，导致MSGN1A花粉中ATP合成受阻，最终使得花粉因能量不足而败育。

Rato等[22]研究发现，花粉管的生长依赖于多种信号通路，其中钙调蛋白是钙离子信号通路中的关键因子。Yang等[23]研究认为，在大豆CMS中，12个DEGs的差异表达可能导致钙信号通路被破坏、花粉发育异常以及细胞质雄性不育系NJCMS1A不育。而在本研究中，由COG对差异表达基因的功能分类可知，有18个DEGs参与了信号传导，其中17个DEGs在MSGN1A中表达下调，表明这些DEGs影响了信号传导机制，也充分证明了信号传导机制对紫花苜蓿的雄性不育性至关重要。

本研究运用转录组测序技术，对紫花苜蓿细胞质雄性不育系MSGN1A及其保持系MSGN1B进行转录组分析和差异表达分析，通过不同数据库对其功能基因进行注释，筛选出与细胞质雄性不育性相关的代谢通路，初步探讨了紫花苜蓿不育系的分子机理，为紫花苜蓿利用不育系的育种工作提供了理论基础。