茅苍术二倍体和同源四倍体叶片转录组比较

张子璇， 谭明谱， 王 将， 张成才， 罗丽娜， 向增旭，①

(南京农业大学： a. 园艺学院， b. 生命科学学院， 江苏 南京 210095)

多倍体植物具有基因剂量效应，其体内重复基因的表达发生变化，导致植物体在生长发育、形态、生理及环境胁迫耐性等方面表现出优势性状[1]。多倍体育种是目前药用植物选育新品种的重要手段，很多药用植物已成功诱导出多倍体且获得的多倍体具有活性成分含量高、抗逆性强等特点。例如：穿心莲〔Andrographispaniculata(Burm. f.) Nees〕同源四倍体叶和茎中穿心莲内酯含量以及茎中脱水穿心莲内脂含量均高于二倍体[2]；甜叶菊(SteviarebaudianaBert.)四倍体的逆境适应能力较二倍体更强，能有效利用弱光和低浓度CO2[3]。与二倍体相比，同源四倍体的叶片发生明显变化，如PlatycodongrandiflormA. De Candolle同源四倍体的叶片长度和宽度、气孔长度和宽度、超氧化物歧化酶和过氧化物酶活性、可溶性糖和蛋白质含量及总叶绿素含量均升高[4]。

茅苍术〔Atractylodeslancea(Thunb.) DC.〕隶属于菊科(Compositae)苍术属(AtractylodesDC.)，为多年生草本植物，其干燥根茎可入药，具有燥湿健脾、祛风散寒和明目的功效[5]，主要产地为江苏、湖北及安徽等省份，道地产区为江苏镇江茅山[6]。近年来，由于掠夺性采挖和生境变化等原因，茅苍术野生资源濒临枯竭，开展茅苍术育种和品质改良研究迫在眉睫。

转录组测序技术能为植物基因表达研究提供数据资料，是研究未完成基因组测序物种的基因及基因功能挖掘的重要手段[7-8]。利用转录组测序技术系统分析茅苍术二倍体和同源四倍体叶片的转录组信息，能进一步明确茅苍术染色体加倍后一系列性状变异的分子机制，为挖掘茅苍术的功能基因提供基础资料。本项目组已经初步研究了不同倍性茅苍术组培苗间的表型差异及生理特征和转录组差异，结果表明：与二倍体相比，茅苍术同源四倍体叶片发生明显的表型变化，包括叶面积指数、叶长、叶宽、叶厚、叶绿素含量、气孔的横径和纵径以及保卫细胞中的叶绿体数显著增加，叶片的可溶性糖、可溶性蛋白质及可溶性淀粉含量升高[9-10]。

鉴于此，以茅苍术二倍体和同源四倍体试管苗为研究对象，采用转录组测序技术对其叶片进行无参转录组测序分析，并在前期研究[10]基础上筛选叶片中的高可信度差异表达unigenes，以揭示茅苍术同源四倍体叶片硕大和抗逆性强的分子机制，并挖掘叶片内直接参与同源四倍体叶片性状变异和抗逆性的相关候选基因。

1 材料和方法

1.1 材料

供试材料为南京农业大学园艺学院中药生物技术实验室前期培养的茅苍术二倍体和同源四倍体试管苗。选取实验室内继代60 d且生长健壮的茅苍术二倍体和同源四倍体试管苗，在生根培养基(添加0.5 mg·L-1NAA的1/2MS培养基)中培养20 d后，分别选取长势基本一致的二倍体和同源四倍体试管苗各5株，每株采集基生叶片2枚，相同倍性的试管苗叶片混匀。按照上述方法重复取样3次。采集的叶片经液氮速冻后置于-80 ℃冰箱中保存、备用。

1.2 方法

1.2.1 总RNA的提取和纯化 采用TRIzol®Reagent(美国Invitrogen公司)分别提取茅苍术二倍体和同源四倍体叶片的总RNA；使用NanoDrop 2000超微量分光光度计(美国Thermo公司)检测总RNA的浓度和纯度；采用DNase Ⅰ〔天根生化科技(北京)有限公司〕纯化总RNA，去除总RNA中残留的基因组DNA；获得的总RNA经Agilent 2100生物分析仪(美国Agilent公司)检测合格后用于后续分析。选取质量合格的总RNA样品(质量在5 μg及以上，质量浓度在200 ng·μL-1及以上，OD260/OD280为1.8～2.2)构建cDNA文库。

1.2.2 cDNA文库构建和转录组测序 采用磁珠法[11]分离总RNA中的mRNA；采用Fragmentation buffer(美国Invitrogen公司)将获得的mRNA随机断裂成长度约200 bp的小片段；以mRNA为模板，采用SuperScript Ⅱ Reverse Transcriptase试剂盒(美国Agilent公司)反转录合成cDNA第1链，并使用10×DNA Polymerase buffer和DNA Polymerase Ⅰ试剂(美国Promega公司)合成cDNA第2链。使用End Repair Mix(美国Enzymatics公司)将cDNA的粘性末端补成平末端，随后在3′末端加上1个A碱基并连接测序接头；采用Illumina HiSeq 2000高通量测序平台(美国Illumina公司)测序，获得原始测序数据。

1.2.3 测序数据过滤及组装 使用Cutadapt v1.16软件过滤原始测序数据，除去污染、含N比高于10%和低质量的reads，获得高质量测序数据(clean data)。由于茅苍术无参考基因组序列，使用Trinity v2.6.6软件对所有clean data进行从头组装，去冗余后得到unigenes。

1.2.4 基因功能注释和差异表达基因分析 使用BLASTx v2.2.25软件将得到的unigenes比对到NR、String、SwissProt和KEGG数据库，分别进行基因功能注释；使用Blast2GO软件(http：∥www.blast2go.com)进行GO注释；使用HMMER v3.2.1软件进行Pfam注释。

使用edgeR v3.24软件分析茅苍术二倍体和同源四倍体叶片基因的表达量差异，将log2FC|≥1且FDR≤0.05作为筛选差异表达基因的标准，其中，FC为差异倍数，FDR为误判率。对筛选出的unigenes进行GO功能分类统计和KEGG富集分析；比较项目组前期对茅苍术二倍体和同源四倍体试管苗全株[10]及本研究中茅苍术二倍体和同源四倍体试管苗叶片的基因差异表达情况，据此筛选出叶片中的高可信度差异表达unigenes。

2 结果和分析

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

茅苍术二倍体和同源四倍体叶片转录组的测序分析结果(表1)表明：原始reads有448 453 428条，clean reads有447 891 390条，所有样本的Q20均大于98%，Q30均大于94%，GC含量在46%以上。

表1 茅苍术二倍体和同源四倍体叶片转录组测序分析

对获得的clean reads进行组装和拼接，共获得750 700个转录本数据，总碱基数364 718 682 bp，最长39 024 bp，最短201 bp，平均长度485.8 bp，N50和N90分别为585和239 bp；去冗余后获得495 797个unigenes，这些unigenes的总碱基数为199 717 026 bp，GC含量46.54%，平均长度402.8 bp，N50和N90分别为413和225 bp(表2)。

表2 茅苍术叶片转录本和unigene分析

2.2 基因功能注释分析

将茅苍术转录组去冗余后获得的unigenes分别与GO、NR、String、SwissProt、KEGG和Pfam数据库进行比对，结果(表3)显示：共注释到259 841个unigenes，占unigenes总数的52.41%。其中，注释到String数据库的unigenes数最多(180 360)，占unigenes总数的36.38%；注释到Pfam数据库的unigenes数最少(66 311)，占unigenes总数的13.37%。

表3 茅苍术转录组unigenes与6个数据库的比对结果

NR数据库中茅苍术与其他植物转录组unigenes序列的比对结果(图1)显示：茅苍术与Cynaracardunculusvar.scolymus(Linn.) Fiori的unigenes序列同源性最高，共有28 614个unigenes序列相似，占该数据库注释unigenes数的20.77%；与向日葵(HelianthusannuusLinn.)的unigenes序列同源性较高，有21 191个unigenes序列相似，占该数据库注释unigenes数的15.38%；与大麦(Hordeumvulgaresubsp.vulgareLinn.)、克里藻〔Klebsormidiumnitens(Kützing) Lokhorst〕、黄石斛(DendrobiumcatenatumLindl.)、旋蒴苣苔〔Boeahygrometrica(Bunge) R. Br.〕、衣藻(Chlamydomonaseustigma)、小立碗藓〔Physcomitrellapatens(Hedw.) Bruch et Schimp.〕、蓖麻(RicinuscommunisLinn.)和胶球藻C-169(CoccomyxasubellipsoideaC-169)的unigenes序列同源性较低，有2 012～5 475个unigenes序列相似，均占该数据库注释unigenes数的4%以下。另外，茅苍术与白梨(PyrusbretschneideriRehd.)、甜菜(Betavulgarissubsp.vulgarisLinn.)和绿藻(Ostreococcustauri)等植物的unigenes序列同源性更低，均占该数据库注释unigenes数的1.8%以下。

图1 NR数据库中茅苍术与部分植物转录组unigenes序列的比对结果Fig. 1 Comparison results of unigene sequences in transcriptomes of Atractylodes lancea (Thunb.) DC. and some plants in NR database

2.3 差异表达unigenes分析

比较茅苍术二倍体和同源四倍体叶片unigenes的表达量，筛选出923个差异表达unigenes。与二倍体相比，同源四倍体中表达量上调的unigenes有317个，表达量下调的unigenes有606个。

2.3.1 GO功能分类结果 GO功能分类结果(表4)表明：茅苍术二倍体和同源四倍体叶片差异表达unigenes的功能被分成生物过程、细胞组分和分子功能3个大类61个小类。在生物过程大类中，注释为细胞过程的差异表达unigenes最多，共有342个；注释为代谢过程、应激反应、生物调节和生物过程调节的差异表达unigenes较多，分别有307、249、219和201个。在细胞组分大类中，注释为细胞、细胞部分、细胞器、细胞器部分和细胞膜的差异表达unigenes较多，分别有362、362、337、232和231个。在分子功能大类中，注释为结合和催化活性的差异表达unigenes较多，分别有291和210个。

表4 茅苍术二倍体和同源四倍体转录组差异表达unigenes的GO功能分类

2.3.2 KEGG富集结果 茅苍术二倍体和同源四倍体共有155个差异表达unigenes注释到KEGG通路中，并对其进行KEGG分类和通路富集分析，结果分别见表5和表6。结果显示：这些差异表达unigenes被注释到代谢、生物系统、遗传信息处理、环境信息处理和细胞过程5个大类中，以注释到代谢中的差异表达unigenes最多(75)，占注释到KEGG通路的差异表达unigenes总数的48.39%。这些差异表达unigenes富集在8条代谢通路中，包括次生代谢物生物合成(ko01110)、淀粉和蔗糖代谢(ko00500)、脂肪酸延伸(ko00062)、光合作用-天线蛋白(ko00196)、NF-κB信号通路(ko04064)、植物激素信号转导(ko04075)、植物MAPK信号通路(ko04016)和核糖体(ko03010)。

表5 茅苍术二倍体和同源四倍体转录组差异表达unigenes的KEGG分类

表6 茅苍术二倍体和同源四倍体转录组差异表达unigenes的KEGG通路富集

2.4 高可信度差异表达unigenes分析

筛选出的茅苍术二倍体和同源四倍体叶片和全株的共有高可信度差异表达unigenes有298个，其中，156个共有高可信度差异表达unigenes有功能注释信息，而且其中的58个在KEGG代谢途径中有注释信息，包括热激蛋白基因(DnaK/Hsp70)、ACC氧化酶基因(ACO)、乙烯受体基因(ETR/ERS)、脂肪酸脱氢酶基因(FAD2)、丙二烯氧化物合酶基因(AOS)、脂氧合酶基因(LOX)、长链-3-羟酰基-CoA脱水酶基因(PHS1/PAS2)和14-3-3结构域蛋白基因(YWHAE)。茅苍术二倍体和同源四倍体叶片和全株的这些共有高可信度差异表达unigenes的表达分析见表7。结果表明：DnaK/Hsp70、ACO、FAD2和PHS1/PAS2基因的表达量在同源四倍体中上调，而ETR/ERS、AOS、LOX和YWHAE基因的表达量在同源四倍体中下调。

表7 茅苍术二倍体和同源四倍体叶片和全株部分共有高可信度差异表达unigenes的表达分析

3 讨论和结论

本研究中，茅苍术二倍体和同源四倍体的unigenes在GO、NR、String、SwissProt、KEGG和Pfam数据库的比对结果表明不是所有unigenes都能在已知数据库中比对出功能，这一结果与目前测序的许多植物的相关研究结果相符[11-13]。这主要是因为多数非模式植物尚未完成全基因组测序工作，致使许多植物的基因组信息不全，不可能通过数据库比对得知所有unigenes的功能。值得注意的是，与NR数据库比对结果显示，茅苍术与同为菊科的C.cardunculusvar.scolymus和向日葵的同源性较高，说明同科植物的unigenes序列同源性更高。

植物多倍化后的性状差异主要由基因表达变化引起，且基因表达只发生微小改变也可能引起表型变异[14-15]。因此，基因差异表达是茅苍术同源四倍体植株性状变异的基础。本研究发现，注释为细胞过程和代谢过程以及细胞、细胞部分、细胞器、细胞器部分和细胞膜的茅苍术二倍体和同源四倍体叶片的差异表达unigenes较多，说明茅苍术的倍性变化会影响其细胞发育和代谢；还有部分差异表达unigenes注释为信号和信号传感器活性，说明茅苍术同源四倍体叶片在信号转导方面与二倍体也存在差异。在其他植物多倍体研究中，酸枣〔Ziziphusjujubavar.spinosa(Bunge) Hu ex H. F. Chow.〕四倍体的差异表达基因富集在糖类和氨基酸代谢、信号传递以及能量代谢方面，并且这些基因的表达量在淀粉和蔗糖代谢及植物激素传导通路中显著上调[15]，杨树(Populusspp.)三倍体中参与糖类代谢和细胞生长的差异表达基因的表达量也上调[16]。本研究的KEGG富集结果表明：茅苍术二倍体和同源四倍体的差异表达unigenes富集在脂肪酸延伸、淀粉和蔗糖代谢及次生代谢物生物合成等通路。淀粉、蔗糖和脂肪均为植物细胞活动提供能量，并且，在植物体内，淀粉和蔗糖代谢、脂肪酸延伸及次生代谢物生物合成过程均伴随着能量转化和物质积累。分析表明：差异表达基因可参与植物体的能量转化和物质积累过程，为细胞发育提供充足的能量和原料、催化相关酶活性，这也可能是茅苍术同源四倍体叶片较二倍体大且厚的一个重要原因。

植物在适应环境胁迫过程中需要能量、物质、激素和酶等共同发挥作用。在受到环境胁迫时，植物体会产生大量的转导信号，与淀粉和蔗糖等代谢过程相关的基因的表达量明显上调，体内的淀粉和蔗糖等代谢过程加强[17]，使植物体适应逆境条件。在本研究的GO功能分类中，注释到生物过程类的细胞过程、代谢过程、应激反应、生物调节和生物过程调节等分类的差异表达unigenes较多，且在KEGG富集结果中，部分差异表达unigenes富集于植物激素信号转导、植物MAPK信号通路和NF-κB信号通路共3个信号转导相关通路中。其中，植物MAPK信号通路是将细胞膜表面信号放大并传递到细胞内，从而引起细胞内转录及代谢水平变化的信号转导途径之一，与植物生长发育和抗逆信号转导具有密切关系[18]。茅苍术同源四倍体叶片在能量转化、物质积累和信号转导方面都与二倍体叶片存在差异，这些差异能够影响植物体对环境变化的响应，据此推断茅苍术同源四倍体在应对逆境胁迫时可能存在更强的适应能力，以保证整个植株包括其药用部位根状茎的正常生长。后续研究应对不同倍性茅苍术药用部位根茎进行转录组和代谢组联合分析，以揭示二者在外观性状、生长发育及抗逆性等方面的分子生物学差异机制，为进一步定位茅苍术的药用成分合成关键基因及通过基因工程手段进行茅苍术育种等奠定基础。

将茅苍术二倍体和同源四倍体叶片和全株[10]共有的差异表达unigenes进行比较，共筛选出298个共有高可信度差异表达unigenes，其中，叶片高可信度差异表达unigenes(包括DnaK/Hsp70、ACO、ETR/ERS、FAD2、AOS、LOX、PHS1/PAS2和YWHAE等基因)可能与同源四倍体出现叶变大、变厚，叶色深绿，气孔增大等表型变化和植株抗逆性相关[19-34]。本研究中，DnaK/Hsp70、ACO、FAD2和PHS1/PAS2基因的表达量在茅苍术同源四倍体叶片和全株中均上调，而ETR/ERS、AOS、LOX和YWHAE基因的表达量在茅苍术同源四倍体叶片和全株中均下调，这可能是茅苍术同源四倍体植株抗逆性强于二倍体的原因。另外，茅苍术二倍体和同源四倍体叶片和全株[10]非共有的差异表达unigenes也值得研究，这些基因可能调控其植株叶片和地下部(如根状茎)的生长发育。

本研究结果显示：茅苍术同源四倍体叶片在细胞发育、代谢、能量转化、物质积累及信号转导方面与二倍体间存在差异，差异表达unigenes可能参与叶片发育调控。高可信度差异表达基因DnaK/Hsp70、ACO、FAD2和PHS1/PAS2的表达量在茅苍术同源四倍体叶片和全株中均上调，而ETR/ERS、AOS、LOX和YWHAE的表达量均下调，这些基因直接参与茅苍术同源四倍体叶片发育和抗逆性，可作为茅苍术倍性育种的关键候选基因。