金线莲应答高温胁迫的转录组学特征分析

梅 瑜，王继华，蔡时可

（1广东省农作物遗传改良重点实验室/广东省农业科学院作物研究所，广州 510640；2广东省道地南药资源保护与利用工程技术研究中心，广州 510640）

0 引言

金线莲Anoectochilus sp.(Wall.)Lindl，别名金线兰、金线虎头蕉、鸟人参和金不换等，以全草入药，性味甘、平，入肝、脾、肾三经，具有清热凉血、消炎解毒和强心利尿等功效，属于珍稀名贵的中药材，享有“药中之王”美称[1-2]。金线莲为兰科开唇兰属草本植物，主要分布于中国、日本、斯里兰卡、印度和尼泊尔等东南亚国家，在中国的台湾、福建和广东等地区有人工标准化种植[3-4]，金线莲的药用价值逐渐得到人们的认可，市场需求日益增加[5]。金线莲喜阴凉湿润环境，适宜生长温度为18-25oC，对生长条件要求苛刻，我国华南地区高温早且时间长，人工种植的金线莲因高温难以越夏，高温严重制约着金线莲产业的发展。

目前金线莲在耐热机制的研究还较少，主要集中在栽培和种苗培育方面。高温胁迫影响植物生长、品质和产量，植物对高温胁迫的响应一直以来也是研究的热点[6]。高温胁迫对植物的伤害机制研究较多，主要集中在生理和生化方面的研究，其中高温对叶绿体膜的伤害较为直接[6-7]。近年来高通量测序技术的发展，利用转录技术解析植物对高温胁迫的响应日益增多，如茶树、坛紫菜和玉米等[8-10]。并且和代谢组、蛋白组组学结合，通过转录组及代谢组分析‘郑单958’和‘先玉335’两个主要玉米品种在响应花粒期高温胁迫下基因表达及代谢物差异，挖掘玉米响应高温胁迫的关键基因及代谢物联合解析其耐热机制[10]。本研究所用的金线莲(NYJ2)是前期筛选出来的耐热资源，比较高温胁迫前后的转录组变化有利于解析其耐热机制[11]。但是，金线莲的遗传背景还不清晰，分子生物方面的研究比较薄弱，基因组和转录组也还未报道。因此，本研究对高温胁迫前后的金线莲叶片进行高通量转录组测序，采用Denovo对获得数据进行组装和比对分析。通过对组装出的unigene进行表达量注释和功能注释(GO、KEGG、KOG)，并对高温胁迫前后的差异表达基因进行GO功能分类、KEGG通路富集分析，来揭示其功能差异，这有利于高温胁迫下差异表达基因调控机制的研究和耐热基因的挖掘。

1 材料与方法

1.1 材料

实验材料为广东省农业科学院作物研究所筛选出的耐热金线莲(NYJ2)。

1.2 方法

1.2.1 材料处理方法 驯化后的NYJ2组培苗定植于装有泥炭土:蛭石=1:1的4×4 cm培养钵内，放于人工气候培养箱中培养，控制培养条件为昼温/夜温(12h/12h)25℃/18℃，光照强度2000 Lx，相对空气湿度(80%+5%)，每天上午定时浇水一次，每株20 mL。培养15天后，挑选生长一致的健壮植株进行45℃高温胁迫处理1 h，以处理前样品为对照，每处理2个重复。取相同部位的叶片，用液氮处理15 min，待叶片冷固后迅速放入-80℃冰箱保存备用。

1.2.2 叶片总RNA提取 使用TRNzol Reagent试剂盒(Invitrogen™)提取高温胁迫前后金线莲4个样品的叶片总RNA，通过1%电泳凝胶检测提取RNA的完整性。采用Invitrogen Qubit® 2.0荧光计及试剂盒(Fluorometer Life Tech Invitrogen，Q32886)对总 RNA进行定量。

1.2.3 转录组测序与数据处理 与广州基迪奥生物科技有限公司合作完成，利用Illumina HiSeqTM2000平台测序。对得到的原始文件clean reads再进行更严格的过滤，去除adaptor序列、去除N的比例大于10%的序列、去除低质量序列，保留下是High quality clean reads。通过短reads组装软件Trinity[12]的denovo，对High quality clean reads从头组装，将具有一定长度overlap的reads连成更长的、不含N的组装片段(Contig)，对Contig进行比对、连接得到两端不能再延长的序列，即为组装出来的unigene。

1.2.4 Unigene基本功能注释 通过blastx将组装出来的unigene序列与Swiss Prot、NR蛋白数据库进行比对，与基因产物进行直系同源分类的数据库COG/KOG进行比对，与基因产物在细胞中的代谢途径及这些基因产物功能的数据库KEGG进行比对（E-value＜0.00001），获得与unigene具有最高序列相似性的蛋白，从而得到该unigene的蛋白功能注释信息。通过Blast2GO[13]软件，在NR注释信息的基础上得到每个unigene的GO注释信息，然后用WEGO软件[14]对所有unigene做GO功能分类统计，从宏观上认识该物种的基因功能分布特征；根据KEGG的功能注释信息获得unigene的Pathway注释信息并进行代谢通路分析。此外，用软件ESTScan对以上数据库都比对不上的unigene预测其编码区及序列方向，得到其编码区的核酸序列（序列方向5'-＞3'）和氨基酸序列，获得unigene的功能域及蛋白家族注释信息，即编码蛋白框(CDS)预测[15]。

1.2.5 Unigene表达量分析 采用RPKM（Reads Per kb per Million reads）计算Unigene的表达量[16]，得到的基因表达量可用于不同样品间的Unigene的差异比较，其计算见公式(1)：

RPKM表示单个unigene的表达量；C表示比对到这个unigene的reads数；N表示比对到所有unigene的总reads数；L表示这个unigene的碱基数。以热胁迫前为对照，进行比较分析。

2 结果与分析

2.1 金线莲的转录组特征

通过Illunima HiSeqTM2000平台进行高通量测序，获得了金线莲叶片应答高温胁迫前后的转录组信息，4个金线莲的转录组样本共计产生38,341,121,400 bp数据，获得255,607,476条reads。对数据进行过滤，去除无效数据，最终4个金线莲样本获得的Clean reads分别为 52049424、55223682、84465826、54304734，Q20（测序准确率＞=99%）百分比在96.14%～96.91%之间，GC百分比在47.73%～49.73%之间，说明测序质量较好。通过Trinity组装软件将各样本的clean reads从头组装，通过比对软件Bowtie[10]将clean reads与unigene进行比对，得到与4个金线莲样本的参考基因比对上的reads为51641412～83706916，能够唯一比对上参考基因的 reads分别为 43737178、46446403、71481627、45406067，分别占各样本总数84.69%、84.84%、85.40%、84.30%；能够比对上参考基因的reads分别为1252595、1524775、2169327、1317954，占各样本总数的2.43%、2.79%、2.59%、2.45%；不能比对上的reads分别为6651639、6777118、10055962、7139527，占各样本总数的12.88%、12.38%、12.01%和13.25%。继续将这些reads组装，4个样本的转录组库共获得了75688个unigene，GC百分比为40.81%，unigene总长度为64225402 nt，最长为20118nt，最小长度为201 nt，平均长度为848 nt，N50为1448 nt，表明组装结果比较理想；其中，长度在201～999 nt之间的unigene有57198个，占总基因数的75.57%，长度在1000～1999 nt之间的unigene有10828个，占总基因数的14.31%，长度在2000～2999 nt之间的unigene有4620个，占总基因数的6.10%，长度超过3000nt的unigene有3042个，占总基因数的4.02%。

2.2 金线莲转录组的注释

通过BLASTx软件，将获得的unigene依次与NCBI中已公布的Nr、Swiss-Prot、KEGG和KOG蛋白数据库得基因信息进行比对(E-value＜0.00001)，得到与给定unigene具有最高序列相似性的蛋白，从而得到该unigene的蛋白功能注释信息。在75688条unigene中，共有28323个unigene被注释，占所有unigene的37.42%，另有47365个unigene未被注释，占总unigene的62.58%。在被注释的unigene中，与Nr数据库的已知编码蛋白具有同源性的基因有28038个，占被注释基因的98.99%，被注释到Swissprot protein数据库的unigene有19149个，占被注释基因的67.61%，被注释到KOG数据库的基因有17532个，占被注释基因的61.90%，10108个unigene被注释到KEGG数据库，占所有被注释基因的35.69%，被注释到GO数据库的基因有17485个，占被注释基因的62.36%。unigene被Nr、Swissprot、KOG、KEGG4个数据库注释，同时被4个数据库注释的unigene有7854个，被数据库单独注释的unigene数目分别为5870个、104个、56个、42个。金线莲转录组中比对上unigene数目最多的前10个物种分别是油棕(Elaeis guineensis)、海枣(Phoenix dactylifera)、野 芭 蕉 (Musa acuminate)、水 稻 (Oryza sativa)、可 可 (Theobroma cacao)、苜 蓿 (Medicago truncatula)、荷 花 (Nelumbo nucifera)、油 菜 (Brassica napus)、葡 萄 (Vitis vinifera)、棉 花 (Gossypium arboretum)，分别占被注释总数的18.75%、16.00%、8.41%、7.14%、3.91%、2.74%、2.55%、2.71%、1.88%、2.50%，注释到其他物种的unigene逐渐减少。

17532个Unigenes被KOG数据库注释的unigene被归为25个蛋白质家族中。其中，7830个unigene注释到一般性功能预测类，这是最大的类群；其次是翻译后修饰、蛋白质转换、伴侣蛋白的unigene有2794个；信号传导机制的unigene有2322个；RNA的加工修饰的unigene有1537个；转录的unigene有1443个；蛋白翻译的unigene有1384个；细胞内运输、分泌和囊泡转运的unigene有1219个；25个蛋白家族中有17个蛋白质家族的unigene数目超过500个；而被归类到细胞流动性、细胞外结构、细胞核结构的unigene较少，分别为15、60和90。

10108个unigenes被KEGG数据库注释，5727个unigene被注释到128个KEGG Pathway中，其中，有512个unigene被注释到核糖体通路中，是数目最多的；其次是被注释到碳代谢通路，有396个unigene；346个unigene注释到氨基酸的生物合成通路；262个unigene注释到淀粉和蔗糖的代谢通路；其他通路被注释的基因富集数目逐渐减少。Gene Ontology(GO)是一个国际标准化的基因功能分类体系，从宏观上表述基因产物的功能，其功能体系涵盖生物学的三个方面：分子功能、细胞组分、生物学过程。根据Nr注释信息，利用Blast2GO[13]软件对unigene进行GO注释信息，在此基础上用WEGO软件[14]对unigene进行GO功能分类统计，结果有17485个unigene被注释到47个GO terms中。其中，被注释到生物学过程的unigene最多，为37170个，尤其是代谢过程(10139 unigenes)、细胞过程(9284 unigenes)和单有机体过程(6412 unigenes)；其次是被注释到细胞组分途径的unigene有32535个，以细胞(8115 unigenes)、细胞部分(8112 unigenes)和细胞器(6724 unigenes)为主；而注释到分子功能过程的基因较少(19889 unigenes)，主要是结合捆绑条目(9180 unigenes)和催化活性条目(8730 unigenes)。

另外，通过MISA软件对所有unigene进行搜索，在75688个序列中，共鉴定出7026个简单重复序列（SSR），其中975个序列至少含有1个SSR位点，如表1。

表1 SSR类型统计表

2.3 金线莲应答高温胁迫转录组差异表达基因分析

利用FPKM值检测样品的重复性，比较样品间表达量的相关性。热处理前后样品的皮尔森相关系数(r)为0.9117和0.9373，表明同一处理的两次平行之间的相关性良好。通过比较高温胁迫处理样本的叶片unigene的表达水平（FDR＜0.05且差异倍数1.5倍以上）。对照样本和高温胁迫处理样本共有777个unigene被差异表达，其中有362个unigene在高温胁迫后表现为上调表达，有415个unigene表现为下调表达。

通过对基因注释结果的进一步分析，发现了15个与光系统Ⅰ(PSⅠ)和光系统Ⅱ(PSⅡ)相关的基因和3个叶绿体基因rbcL(Unigene0012214、Unigene0023 0681、Unigene0010382)、CAT 编码基因(Unigene002 1709)、GAPDH基因(Unigene0055305)、2个PLD基因(Unigene0026596、Unigene0036040)、CYP 编码基因(Unigene0030805)等，其中光系统Ⅱ的3个基因和CYP编码基因在高温胁迫后下调，其余基因均显著上调。此外，MBL基因(Unigene0007875)、Sporamin B基因(Unigene0018863)、PI基因(Unigene0005251)、ZFP 编码基因的表达在高温胁迫后显著上调；液泡铁转运蛋白基因、吲哚乙酸氨基化合成酶基因基因、膜相关激酶调节因子基因均下调表达。

2.3.1 高温胁迫前后差异表达基因分析 为了进一步研究分析高温胁迫下组间样本的差异基因的生物学功能，从中寻找与高温胁迫相关的基因，将两个比较组中所有的差异基因与GO数据库的各个term进行比对（图1）。

图1 差异基因的GO分类及表达

结果生物学过程((Biological Process)的GO term富集的最多，有620个terms。按富集基因数量依次列出较多的10个terms为：代谢过程(181 unigenes)、细胞过程(158 unigenes)、单一有机体过程(142 unigenes)、细胞代谢过程(126 unigenes)、有机物质代谢过程(120 unigenes)、主要代谢过程(98 unigenes)、单个有机体代谢过程(95 unigenes)、单个有机体的细胞过程(93 unigenes)、高分子代谢过程(77 unigenes)、细胞高分子代谢过程(57 unigenes)。其次是分子功能(Molecular Function，216 terms)富集，基因数量较多的10个term为：催化活性(169 unigenes)、结合绑定(155 unigenes)、有机环化合物结合绑定(84 unigenes)、杂环化合物结合(80 unigenes)、离子结合(64 unigenes)、阳离子绑定 55 unigenes)、水解酶活性(53 unigenes)、氧化还原酶活性(49 unigenes)、转移酶活性(45 unigenes)、小分子结合绑定(41 unigenes)等。

富集term最少的是细胞组分(Cellular Component，85 terms)，基因数量较多的10个term：细胞部分(153 unigenes)、细胞(153 unigenes)、细胞内(145 unigenes)、细胞内部分(142 unigenes)、细胞内膜细胞器(125 unigenes)、膜旁细胞器(125 unigenes)、细胞内细胞器(125 unigenes)、细胞器(125 unigenes)、膜(83 unigenes)、胞质部分(75 unigenes)。

在生物学过程中，代谢过程显著上调的基因有105个，显著下调的有76个；细胞过程显著上调94个unigenes，显著下调64个unigenes；单一有机体过程显著上调基因79个，显著下调基因63个，其余term逐渐减少。分子功能中，差异基因在催化活性途径显著上调86个，显著下调83个；在结合绑定term中显著上调87个，显著下调68个；在膜中显著上调49个，显著下调34个。在细胞组分中，差异基因在细胞和细胞部分两个term中显著上调的数目均为92个，显著下调的为61个；在细胞器term有68个unigenes显著上调，57个unigenes显著下调；49个基因在膜中显著下调，34个下调。

2.3.2 高温胁迫前后差异表达基因的KEGG分析 自然界生物体内的不同基因相互协调的，基于代谢途径的分析有助于对基因的生物学功能进一步了解，通过KEGG数据库对差异基因进行富集分析。139个差异基因被注释到75个代谢途径中，其中光合作用途径中被注释的基因最多（27个），有16个基因显著上调，11个基因显著下调；其次是碳代谢途径，有23个基因与碳代谢有关。从富集前20的代谢途径中可以看出（图2），前4个途径富集结果较好，且上调的基因与光合作用中的生物固碳有关。CYP81E1/E7、PLD、CHS、CYP、GH3、FabF、FAR、Lhca、LHC、pepA、CAlM、BKI1、JAZ等基因在KEGG途径的特定阶段出现。

图2 差异基因的KEGG富集结果

3 讨论与结论

高温胁迫下，植物的光合作用途径、碳代谢途径、光合固碳作用、代谢途径、氨基酸代谢途径、黄酮类代谢途径、酯类代谢途径、萜类合成途径和脂肪酸合成代谢等途径相关基因表达量受影响而变化[17-19]。本研究借助高通量测序技术较早的对高温胁迫前后的金线莲进行转录组测序和解析，共获得75688个unigene，其中28323个在数据库中得到功能注释。17532个unigene被KOG数据库归类到25个家族，与128个Pathway有关；其中，17485个unigene被注释到GO三大功能分类下的47个条目中。对比高温胁迫前后共有777个基因表达差异明显，大部分为光合途径和抗氧化途径相关，KEGG注释结果显示光合作用的差异基因最多。其中与PSⅠ、PSⅡ相关的基因15个、3个叶绿体基因rbcL、2个PLD基因，CAT编码基因、GAPDH基因、CYP编码基因各1个，除了光系统Ⅱ的3个基因和CYP编码基因在高温胁迫后下调外，其余基因均显著上调。

本研究获得的psaA(Unigene0004500；Unigene 0047840)和psaB(Unigene0004499)基因是PSⅠ的核心亚基，由它们编码的蛋白A1和A2是PSⅠ复合体中心的组成部分，能够束缚一些还原因子从而减缓上游的移动电子到psaC；从PSⅠ到NADP+穿梭电子的铁硫蛋白为铁还原蛋白，其复合体镶嵌在psaC上，具有电子转移载体的作用，由psaD(Unigene0006392)和psaE(Unigene0013764)翻译的蛋白能够协助电子从铁还原到PSⅠ[20]，这两个基因在高温胁迫后上调，而且它们是PSⅠ的受体位点亚基，金线莲可通过提高或激发PSⅠ的功能来抵御高温。另外，PSⅡ是应答外界胁迫最为敏感的光合系统，参与水的裂解、释放氧是它的重要功能。本研究共获得10个参与PSⅡ的unigene，叶绿体基因psbB、psbC、psbM、psbD参与编码叶绿素结合蛋白CP47、反应中心蛋白D2，前3个基因下调表达；有7个基因上调表达，psbO、psbP、psbQ分别编码叶绿体类囊体膜放氧复合体的外周蛋白OEE1、OEE2、OEE3来维持放氧反应[21]；编码放氧蛋白的3个基因在高温胁迫后显著上调，与金线莲的高温反应正相关，这与温度胁迫诱导的甘蔗psbO和psbP基因的表达一致[22]。此外，据报道psbR连接psbP和psbQ蛋白组成放氧复合体[23]；psbY是一种新型的与锰结合的低分子蛋白质，位于PSⅡ外围，靠近细胞色素b559的α和β亚基[24]。

本研究中获得的3个叶绿体基因rbcL在高温胁迫后表达下降，这与学者在研究高温胁迫下龙须菜rbcL转录表达一致，同时Rubisco活性降低，当使用外源激素SA和MJ50时，Rubisco活性逐渐恢复、rbcL继续表达[26]。可能由于rbcL来源于叶绿体DNA，属母性遗传保守性较高，Rubisco在植物光合作用中具有双重作用，既有活化Rubisco的活性，也有ATP水解酶的活性[25]。此外，还获得了3个磷脂酰胆碱磷脂水解酶D(phospholipaseD，PLD)基因，在高温胁迫后均差异表达，据报道PLD的水解产物PA能参与植物的热激反应中细胞信号传导，并能积极调动和激发细胞内的抗性反应来应答逆境[27-29]，YANG等[30]发现高温胁迫下的高山离子芥大量表达CbPLD基因，其产物P在膜脂降解中有比较重要的作用[31]，推测该基因可能参与了金线莲应答高温胁迫的调控。

作为糖酵解反应酶之一，3-磷酸甘油全脱氢酶GAPDH基因缓解外界胁迫导致的氧化胁迫中的作用明显，在植物抗病和非生物胁迫中具有重要的作用。有研究表明，在拟南芥中GAPCs能够与PLDδ作用产生H2O2来抵抗ROS和干旱胁迫，并阐明了膜质的信号传递、能量代谢及生长调控这三个途径的关系[32]。转ERD15基因的烟草和翦股颖的GAPDH基因干旱胁迫下均被诱导表达，并且与其对应的蛋白的表达量也显著增加[8,33-34]。本研究中GAPDH编码基因在高温胁迫后显著上调表达，可见金线莲抵抗高温胁迫与该基因密切相关。热激蛋白(HSPs)家族的基因对高温、低温的非生物胁迫具有明显的应答反应，金线莲中与HSP合成有关的基因(Unigene0029507、Unigene0052917)，在高温胁迫后显著上调，研究表明HSP家族的转录因子参与调控金线莲对高温胁迫的响应[9,35]。

综上所述，转录组分析是鉴定植物转录调控最为有效和低成本的工具之一，也是非模式植物功能基因组学研究的重要手段之一，是本草基因组学重要组成之一。本研究可为兰科植物分子生物学研究提供参考。此外，金线莲资源的收集、保护、利用，以及种质资源创制是今后研究的重点，为金线莲产业可持续发展提供保障。