辛建攀 李燕 赵楚 田如男

（南京林业大学风景园林学院，南京 210037）

工农业的快速发展、生活污水的不合理排放及化石燃料燃烧等导致我国城乡水体被重金属污染，妨碍了水资源的利用。镉（Cd）是污染水体中常见的重金属污染物之一，分布广泛、难降解，具有极强的生物毒性，能够通过食物链积累与生物放大作用威胁人类及其他水生生物健康［1］。与物理、化学方法相比，植物修复技术因环境干扰小、二次污染易控制及成本低等优势而备受关注，被认为是当前重金属污染环境修复领域的研究热点之一。耐性植物种类的选择是Cd污染水体进行植物绿色修复的基础环节之一。因此，深入研究植物防御Cd胁迫的分子机制可以为植物修复Cd污染提供理论参考。

梭鱼草（Pontederia cordata）属雨久花科（Pontederiaceae）梭鱼草属（Pontederia）多年生观赏湿地植物，其根系发达、生物量大，顶生穗状花序，小花蓝紫色，具有很高的园林观赏价值，是构建人工湿地植被的优良材料。该植物不但可以有效去除水体氮、磷［2-4］，并抑制藻类生长［5-6］，而且还能够有效去除水体Cd、铜及铅等重金属［7-10］。研究表明，根系富集、细胞抗氧化还原系统的调控活性、非蛋白巯基肽含量等因素与梭鱼草Cd胁迫耐性密切相关［10-12］。本课题组研究发现，Cd能够在转录水平上明显改变梭鱼草叶片中黄酮、类黄酮及花青素等次生代谢物生物合成相关基因的表达模式，同时一些差异表达基因显著富集到苯丙烷代谢途径，相似的结果也被发现于重金属胁迫下的杞柳（Salix integra）［13］和旱柳（S. matsudana）［14］等多种植物，表明苯丙烷代谢途径参与梭鱼草Cd胁迫响应。

苯丙烷代谢途径是植物次生代谢途径中的重要分支之一，是多数防御性次生代谢物（如黄酮类、木质素及水杨酸等）的主要来源，能够有效减缓细胞中活性氧的积累，避免细胞形态结构与生理代谢功能被破坏，因此该途径与植物非生物胁迫耐性密切相关［15-18］。鉴于此，本研究利用Illumina Hiseq测序平台对Cd胁迫下梭鱼草叶片进行转录组测序，并进一步挖掘苯丙烷代谢途径中相关基因，为今后深入研究苯丙烷代谢途径参与梭鱼草防御重金属植物毒性的分子机制奠定理论基础，同时也为梭鱼草在修复重金属污染水体中的实际应用提供理论依据。

1 材料与方法

1.1 材料

梭鱼草根状茎的萌蘖生长约2个月后，于5月份选取长势一致、健康的植株（株高45-50 cm）并冲洗干净，然后移植于盛有2 L 1/2 Hoagland营养液的塑料桶中进行适应性培养。每桶1株，共25株，培养周期为20 d，期间适时补充营养液至初始体积。适应性培养结束后，利用0.44 mmol/L Cd2+对植株胁迫处理，于0 h、1 h、3 h、6 h、12 h、24 h和48 h时采取植株顶端成熟叶片用于转录组测序，每个时间点设置3次生物学重复。

1.2 方法

1.2.1 文库构建与Illumina HiSeq测序 采用Trizol试剂盒提取各时间点叶片总RNA，利用NanodropND 2000检测RNA纯度，同时使用Agilent Bioanalyzer 2100系统的RNA nano 6000分析试剂盒精确检测RNA完整性。叶片RNA样品检测合格后，用带有Oligo（dT）的磁珠富集真核生物mRNA。随后加入fragmentation buffer将mRNA打断成短片段，以mRNA为模板，用六碱基随机引物合成一链cDNA，然后加入缓冲液、dNTPs和DNA polymerase I和RNase H合成二链cDNA，再用AMPure XP beads纯化双链cDNA。库检合格后，将cDA文库按照有效浓度及目标下机数据量的需求pooling后进行Illumina HiSeq测序。

1.2.2 组装和功能注释 为了确保转录组测序结果及生物信息学的分析质量，对测序得到的原始序列（raw reads）进行处理，（1）去除含有接头（adapter）的reads；（2）去除N（无法确定碱基信息）比例大于10%的reads；（3）去除低质量reads。采用Trinity对过滤获得clean reads进行拼接，从而获取后续生物信息学分析的参考序列［19］。

基因功能注释采用七大数据库：NR（NCBI non-redundant protein sequences）、NT（NCBI non-redundant nucleotide sequences）、PFAM（portein family）、KOG（eukaryotic ortholog groups）、Swiss Prot（A manually annotated and reviewed protein sequence database）、KEGG（Kyoto Encyclopedia of Genes and Genomes）和GO（gene ontology）。

1.2.3 基因表达量分析 根据Li等［20］的方法对reads和unigenes的比对结果进行统计，获得每个样品比对到每个基因上的readcount数目，并对其进行FPKM（expected number of fragments per kilobase of transcript sequence per millions base pairs sequenced）转换。采用DEGseq软件进行差异基因筛选，差异基因的筛选条件为：表达差异倍数|log2fold change|＞1［fold change=SDt/CK，其中SDt为取样时间点，CK为对照组（0 h）］，显著性P-value＜0.05。同时，结合KEGG代谢通路分析，本研究对0 h、6 h、48 h时梭鱼草叶片苯丙烷代谢途径相关差异基因进行挖掘，从转录水平上判定梭鱼草叶片响应Cd胁迫的动态变化过程。

1.2.4 差异表达基因的趋势化分析 根据Ernst&Bar-Joseph的方法［21］，利用STEM（short time-series expression miner）软件（设置P≤0.05为差异显著的表达模式）对Cd胁迫下各时间点的差异表达基因进行趋势分析。

1.2.5 编码区序列（coding sequence，CDS）预测 按照Swiss-Prot和Nr蛋白库的优先级顺序，将梭鱼草叶片转录组unigenes进行Blast比对。从比对上的序列中提取出transcript开放阅读框（open reading frames，ORF），并将CDS序列翻译成氨基酸序列。未比对上或比对上未预测出结果的序列则采用Estscan（3.0.3）预测其ORF，根据CDS序列翻译获得氨基酸序列。

2 结果

2.1 转录组组装

本研究共建立了21个cDNA文库，得到了970 961 618个raw reads。经过严格的数据过滤和质量评估，共得到926 875 418个clean reads。所有错误率均低于0.10%，且GC值介于53.54%-54.41%。同时，不低于97.00%的clean reads在Q20水平上有较高分数；不低于93.00%的clean reads在Q30水平上有较高分数。这些结果表明转录组测序质量较高，可用于后续生物信息分析。最后，共得到221 392个高质量的unigenes，其总长度为388 100 176 bp，平均长度为1 753 bp，最大长度为25 638 bp，最小长度为201 bp。转录本和unigenes的长度分布见 （表1）。

表1 拼接转录本和unigenes的长度分布Table 1 Length distribution of assembled transcripts and unigenes

2.2 基因注释和功能分类

为了对这些组装的unigenes进行验证和功能注释，所有unigenes以临界值小于1e-5比对到7个数据库（表2）。在221 392个unigenes中，158 229个 （占总数的71.47%）unigenes和121 487个（占总数的54.87%）unigenes分别比对到NR和SwissProt，这两个数据库具有最多的unigenes注释。

表2 梭鱼草叶片unigene功能注释Table 2 Unigene annotation of P. cordatas leaves

通过与NR比对注释，可以获取梭鱼草基因序列与近缘物种基因序列的相似性。如图1所示，梭鱼草与小果野芭蕉（Musa acuminata）基因序列相似度最高（31.50%），其次是棕榈（Elaeis guineensis）（19.70%）、海藻（Phoenix dactylifera）（15.80%）、菠 萝（Ananas comosus）（6.60%）及 芦 笋（Asparagus officinalis）（2.40%），它们均为单子叶植物。

图1 Unigene的NR比对同源性分析Fig.1 Sequence homology of the unigenes against NR database

在所有的unigenes中，共有98 391个unigenes（44.44%）被注释到GO。GO可以分为3个本体，包括生物过程（biological process，BP）、细胞组分（cellular component，CC）与 分 子 功 能（molecular function，MF）。其中，共有270 742个unigenes被分到25个BP功能组，175 020个unigenes被分到20个CC功能组，125 388个unigenes被分到10个MF功能组（图2）。在BP中，前3个功能组分 别 为“cellular process”“metabolic process”“single-organism process”，占该组unigenes总数的比例分别为22.06%、20.08%和16.16%。在CC中，共有19.17%、19.17%、12.90%、12.56%、11.10%、10.24%、7.00%的unigenes被分别归类为“cell”“cell part”“organelle”“macromolecular complex”“membrane”“membrane part”“organelle part”。在MF中，最前面的两个亚类分别为“binding”（45.87%）和“catalytic activity”（36.66%），其比例明显高于其他亚类，如“antioxidant activity”（0.39%）和“metallochaperone activity”（0.03%）。

图2 Unigenes在GO 分类下一层级的分布Fig.2 Distribution of unigenes to GO level 2

与KOG分类相关的注释序列可用来评价转录文库的完整性和注释过程的有效性。共有49 597个unigenes归属于KOG分类（图3）。在这25个KOG类别中，最大组为“posttranslational moification，protein turnover，chaperones”，其unigenes数量为6 827个（13.76%），其 次 是“general function prediction only”（13.00%）、“translation，ribosomal structure and biogenesis”（9.22%）和“signal transduction mechanis- ms”（7.96%）。然而，“cell motility”和“extracellular structures”是最小功能组，其unigenes数量占比分别为0.06%和0.09%。

图3 Unigenes的KOG功能注释分布情况Fig.3 Distribution of unigenes for KOG functional annotation

以padj＜0.05且|log2foldchange|＞1为筛选阈值，利用DESeq2将所有unigenes与KEGG数据库进行比对分析。在221 392个unigenes中，共有54 087个unigenes明显匹配到该数据库中，并被分为19个KEGG途径。这些途径包括5个亚组，依次是“metabolism”“genetic information processing”“cellular processes”“environmental information processing”“organismal systems”，其所含的unigenes数量分别为31 186个（57.66%）、14 268个（26.38%）、3 621个（6.70%）、2 582个（4.77%）、2 430个（4.49%）（图4）。

图4 Unigenes的KEGG功能注释的分布情况Fig.4 Distribution of unigenes for KEGG functional annotation

2.3 差异表达基因的趋势分析

如图5所示，Cd胁迫期间梭鱼草叶片差异基因（20 025个）的表达模式可以分为6种。其中，丰度最大的表达模式是Cluster 3，共有5 832个差异表达基因表现出一致的表达趋势，在Cd胁迫3 h、12 h和24 h时，这些基因均上调；在1 h、6 h和48 h时维持在初始水平。Cluster 5的丰度次之，该模式中共有3 733个差异表达基因富集；与0 h相比，Cd处理1-48 h时差异基因均表现为上调趋势，并具有较高的表达量，其可能在梭鱼草应对Cd胁迫中具有重要作用。Cluster 1由3 153个差异表达基因组成，这些基因（除48 h外）的表达量均低于初始水平。Cluster 2共由2 642个差异表达基因组成，这些基因在Cd处理6 h、12 h和24 h时维持在初始水平，而在其他时间点的表达均下调。Cluster 0共由2 631个差异表达基因，其随着各处理时间的延长呈下调趋势。Cluster 4模式中差异表达基因数量最少，这些基因在12 h和24 h时维持在初始水平，而在其他时间点均上调表达。

图5 Cd胁迫下梭鱼草叶片差异基因的表达趋势Fig.5 Tendency graphs of differentially expressed genes in P. cordata’s leaves with cadmium exposure

2.4 梭鱼草响应Cd胁迫过程中叶片次生代谢物合成途径相关酶的鉴定

如表3所示，Cd胁迫下梭鱼草叶片unigenes参与苯丙素、萜类化合物骨架、类固醇、类黄酮、花青 素、玉米素及油菜素内酯等生物合成相关的31个次生代谢通路。其中，苯丙素生物合成（ko00940）代谢通路中unigenes数量最多，为696条；缬氨酸、亮氨酸和异亮氨酸降解（ko00280）、萜类化合物骨架生物合成（ko00900）、类胡萝卜素生物合成（ko00906）代谢通路中的unigenes数量依次次之。参与油菜素内酯生物合成（ko00905）、二萜生物合成（ko00904）、咖啡因代谢（ko00232）、花青素生物合成（ko00942）、甜菜红碱生物合成（ko00965）、异黄酮生物合成（ko00943）及吲哚生物碱生物合成（ko00901）等代谢通路的unigenes数量较少。

表3 Cd胁迫下梭鱼草叶片转录组unigenes次生代谢KEGG通路注释Table 3 Unigenes for KEGG pathway annotation involved in the secondary metabolites in the leaves of P. cordata with cadmium exposure

2.4.1 木质素生物合成相关的差异表达基因 木质素单体生物合成代谢通路首先是莽草酸生物合成代谢通路，其次是经苯丙氨酸经脱氨基、甲基化、羟基活化等过程的苯丙烷代谢途径，然后是对羟基肉桂酸及其辅酶A脂类化合物氧化还原为木质素单体，最后强化脱氢聚合生成木质素，此过程涉及到苯丙氨酸解氨酶（PAL）、肉桂酸-4-羟基化酶（C4H）、4-香豆酸辅酶A连接酶（4CL）、查耳酮合酶（CHS）、过氧化物酶（POD）、肉桂酸脱氢酶（CAD）及咖啡酰莽草酸酯酶（CSE）等多种酶的参与。在6_0 h组中，编码对-羟基苯基丙烷此类型木质素合成的关键酶共有3个，包括1条4CL序列（Log2fold change=-9.48）、POD序列（Log2fold change=2.71）及CAD序列（Log2fold change=-5.87）（图6）。在48_0 h组中，梭鱼草叶片中与愈创木基木质素生物合成相关的关键酶共有2个，包括3条C4H序列（Log2fold change=5.22-7.89）和1条CSE序列（Log2fold change=2.26）。在48_6 h组中，梭鱼草叶片中编码愈创木基木质素生物合成途径的关键酶共有5个，包括2条PAL序列（Log2fold change=5.22；6.17），1条咖啡酰辅酶A-O-甲基转移酶（CCoAOMT）序列（Log2fold change=3.93），6条C4H序列（Log2fold change=6.19-7.84），1条CSE序列（Log2fold change=2.26）及4条4CL序列（Log2fold change=2.23-9.69）（表4）。

图6 Cd胁迫下梭鱼草叶片木质素与黄酮类化合物生物合成途径Fig.6 Lignins and flavonoids biosynthesis pathways in the leaves of P. cordata with cadmium exposure

2.4.2 黄酮类生物合成相关的差异表达基因 黄酮类化合物是一类含有C6-C3-C6基本结构骨架的化合物，起源于苯丙烷代谢途径，涉及到查耳酮合酶（CHS）、查尔酮异构酶（CHI）、二氢黄酮3-羟化酶（F3H）等多种酶。在6_0 h组中，梭鱼草叶片中编码黄酮类生物合成相关的关键酶仅有1个，为二氢黄酮3-羟化酶（F3H）（Log2fold change=-2.59）。在48_0 h组中，梭鱼草叶片中编码黄酮类化合物合成相关的关键酶有3个，包括3条C4H序列（Log2fold change=5.22-7.89）、1条 黄 酮 醇3-O-甲 基 转移酶（F3OMT）序列（Log2fold change=4.39）和1条（类黄酮-3，5′羟化酶）F3′ 5′H序列（Log2fold change=-2.29）。在48_6 h组中，梭鱼草叶片中编码在黄酮类化合物合成相关的关键酶有5个，包括6条C4H序 列（Log2fold change=6.19-7.84），1条咖啡酰辅酶A-O-甲基转移酶（CCoAOMT）序列（Log2fold change=3.93），1条花青素合酶（ANS）序列（Log2fold change=8.73），1条IFS序列（Log2fold change=-3.67）及3条花青素-3-O-葡糖基转移酶序列（3GT）（Log2fold change=3.13-10.90）（表4）。

表4 Cd胁迫下梭鱼草叶片中编码木质素与黄酮类生物合成相关酶的unigenesTable 4 Unigenes enconding enzyme involved in lignin and flavonoid biosynthesis in the leaves of P. cordata with cadmium exposure

2.5 CDS序列分析

如图7和图8所示，随着序列长度的增加，unigene数量呈减少趋势。梭鱼草叶片unigenes通过Blast比对得到168 355条CDS序列。其中，其中长度介于101-1 000 nt的序列有124 583条，占总数的74.00%；长度介于1 001-2 000 nt的序列有32 375条，占总数的19.23%。同时，未比对上unigenes则通过Estscan软件预测其ORF，共得到84 673条序列，其长度介于50-13 530 nt。如图7所示，CDS序列长度主要集中在50-900 nt之间，共有80 634条，占比为84.32%。

图7 Cd胁迫下梭鱼草叶片所有unigenes的Blast分析CDS长度分布Fig.7 Blast CDS length distribution of unigenes in the leaves of P. cordata with cadmium exposure

图8 Cd胁迫下梭鱼草叶片所有unigenes的Estscan分析CDS长度分布Fig.8 Estscan CDS length distribution for transcriptome in the leaves of P. cordata with cadmium exposure

3 讨论

转录组测序为研究重金属胁迫下植物细胞生理代谢途径和挖掘抗逆基因等提供了丰富信息。本研究通过对Cd胁迫下梭鱼草叶片进行转录组测序，获得了221 392条高质量的unigenes，其中共有170 175个unigenes被成功注释到数据库中，但仍有23.13%的unigenes没有获得比对结果，这可能与序列片段过短、注释信息缺乏及物种特异性等因素有关［22］。N50值越大且长度不低于800 bp，表明转录组序列完整性较好；Q30不低于80.00%时，即可认为转录组测序质量可靠［23-24］。梭鱼草叶片unigenes平均长度为1 753 bp，N50为2 609 bp，明显高于赶黄草［25］（Penthorum chinense）（1 119 bp；1 601 bp）、张溪香芋（Colacasia escuenla ‘Zhang Xi’）［26］（875.11 bp；1 658 bp）、 百 香 果（Passiflora edulia）［27］（1145 bp；1 991 bp）及 滇 黄 精［28］（Polygonatum kingianum）（770.38 bp；1 350 bp）等多种植物。同时，本研究所得到的Q30均大于91.00%，表明梭鱼草叶片转录组测序和拼接质量较高，可以确保unigenes功能注释能获得较为准确、可靠的比对结果，因此可利用转录组数据进行后续分析。此外，本研究发现，共有49 597条unigenes归属于25个KOG分类，其中注释到“secondary metabolites biosynthesis，transport and catabolism”（次生代谢物生物合成、运输和分解代谢）类群的unigenes共有785条。KEGG通路注释结果表明，Cd胁迫下梭鱼草叶片unigenes参与黄酮类、萜类及生物碱等次生代谢物生物合成相关的31个次生代谢通路，相似的结果也被发现于苹果［29］（Malus domestica）、多裂骆驼蓬［30］（Peganum mulstecutm）及夏枯草（Prunella vulgaris）［31］等植物。Unigenes注释信息的完成为今后进一步考察次生代谢途径在梭鱼草防御重金属胁迫中的作用及其机制提供了理论基础。

苯丙烷代谢途径是植物体内多种防御性次生代谢产物（如黄酮类、木质素及水杨酸等）的主要来源。本研究中，梭鱼草叶片苯丙烷代谢途径是受Cd胁迫影响最为显著的次生代谢途径，相似的结果也被发现于逆境中的大豆［32］、棉花［33］及牛皮杜鹃（Rhododendron chrysanthum）［34］等多种植物，表明苯丙烷代谢途径在植物抵抗非生物胁迫中可能发挥了重要防御作用。PAL是苯丙烷代谢途径的关键酶，可为木质素和黄酮类化合物的生物合成提供代谢底物［16］。在48_6 h组中，梭鱼草叶片中有2条上调表达的PAL序列被鉴定出来，这可能有利于增强PAL酶活性，从而促进下游次生代谢物的合成与积累［35］。高等植物PAL家族中这种仅有少量基因表达而其余同源基因沉默的现象［35］也被发现于赶黄草［25］、夏枯草［31］及八爪金龙［36］（Ardisia crispa）等多种植物。

C4H、4CL也是木质素与黄酮类化合物生物合成途径中上游阶段的关键酶，可催化肉桂酸形成对-香豆酸辅酶A。在48_0 h和48_6 h组中，梭鱼草叶片中共有6条差异表达的C4H序列被鉴定出来，其中主要表达的序列有1条，这可能有利于增强整条代谢通路的通量。作为一类复杂的酚类聚合物，木质素在植物应对重金属胁迫中发挥了重要的防御作用［37-38］。木质素生物合成起始于苯丙氨酸，进入苯丙酸途径后经过一系列的脱氨基、羟基化、甲基化与还原等步骤生成木质素单体，然后经POD、漆酶和酚酶的催化作用聚合为高分子木质素［39］。我们发现，对-羟基苯基木质素生物合成过程中，差异表达的4CL、CAD和POD基因均下调表达；作为单子叶植物体内主要的木质素类型，梭鱼草叶片中与愈创木基木质素生物合成相关的CSE、4CL及CCoAMOT基因上调表达，推测该植物可能通过调控愈创木基木质素的生物合成以吸附叶片细胞中过量的Cd2+［40］。黄酮类化合物作为一类重要的抗逆性物质，在植物抵御重金属胁迫过程中发挥了抗氧化作用［17，41］。CHS、F3H、查尔酮异构酶（CHI）、二氢黄酮醇4-还原酶（FDR）及无色花青素双加氧酶（LDOX）等多种酶参与植物黄酮类化合物的生物合成。本研究通过进一步筛选从各处理组梭鱼草叶片中共获得了15个与黄酮类化合物生物合成的差异表达基因，共编码8个关键酶。其中，包括1个苯丙烷代谢公共途径所需酶（C4H）和7个黄酮类化合物生物合成所需酶（ANS、3GT、IFS、F3H、F3′ 5′H、F3OMT、CCoAMOT），其中PAL、C4H、ANS、3GT、CCoAMOT、F3OMT及POD基因上调表达，这可能与叶片中Cd积累与解毒过程有关［17，42-43］，可以作为今后梭鱼草防御Cd植物毒性机制研究的方向之一。本研究对Cd胁迫下梭鱼草叶片进行了转录组测序，构建了梭鱼草叶片转录组数据库，为深入研究次生代谢途径在梭鱼草防御重金属植物毒性中的作用机制提供了基础数据。

4 结论

在Cd胁迫下，梭鱼草叶片unigenes参与萜类化合物、类黄酮类、生物碱等化合物生物合成相关的31个次生代谢通路，同时叶片差异表达基因具有3个显著的表达模式。Cd胁迫下梭鱼草叶片苯丙烷代谢途径中共有26个差异表达基因分别编码3个苯丙烷公共代谢途径关键酶和9个黄酮类化合物生物合成关键酶，其中CSE、CCoAMOT及3GT等基因在防御Cd胁迫中可能具有重要作用。