随着工业的发展，镉（Cd）成为全球污染最主要的重金属之一，如工业废气会将Cd带到气流层中，并以大气沉降、降雨以及下雪等方式进入土壤中，造成土壤Cd污染；同时，农业上的污水灌溉以及含Cd药剂的使用都会使Cd进入土壤，从而加重土壤的Cd污染。此外土壤中的Cd主要通过根的吸收进入植物体内，然后与根细胞的细胞壁结合，后经木质部运输到茎、叶、果实等器官中积累，抑制植物的生长发育^[1]。Cd可以通过食物链进入动物和人体内，对机体造成严重伤害。过量Cd积累导致肺纹理增多、紊乱而模糊，白细胞和中性粒细胞数量增多，引起肺泡Ⅱ型上皮细胞毛细血管内皮细胞受损，使肺泡间隔增厚，间质纤维增生，也会直接作用于骨骼，使有机体骨骼矿化发生障碍、骨钙溶出的增加以及骨胶原和骨的固化作用异常等^[2-3]。20世纪60年代发生在日本神通川流域的“骨痛病”，原因就是当地居民食用富含重金属Cd的大米造成的。在我国也发生了由于土壤Cd污染日益严重出现了例如“镉大米”，农作物减产等一系列恶性事件^[4]。

考虑到传统的物理、化学修复重金属污染的方法中存在不可避免的缺陷^[5]，植物修复以其成本低、安全可靠、对环境干扰小且改善生态环境等优点而备受关注^[6]，木本植物因为其具有庞大的根系、巨大的生物量、发达的维管组织、坚固的木质组织、以及高效的蒸腾作用，并且可以有吸收多种重金属的特点而成为环境重金属污染治理的重要方向^[7]。旱柳（Saliz matsudana）作为一种常见的木本植物，具有生长快、易繁殖、生物量大、根系发达、再生能力强和对多种重金属耐受的特点^[8]，此外旱柳的叶与根等新陈代谢较快的器官可以大量积累重金属Cd，而其它的主要营养储存器官（果实、籽粒、块茎等）对于重金属积累较少，因此旱柳具有巨大的治理重金属污染的潜力^[9]。

当前对旱柳重金属耐受性的研究主要集中在对重金属的植物体内富集情况的调查、影响因素、生理变化等方面，对重金属富集后植物耐受的分子机理研究较少，尤其在转录水平上研究很少。Claudia等^[10]对柳树耐Cd无性系叶片中Cd的沉积研究发现Cd主要分布于嫩叶的叶尖以及老叶的叶颈，水培时间较短的旱柳主要以嫩叶为主。目前我国土壤一般Cd污染浓度为2.5 mg/L，徐爱春^[11]研究表明Cd浓度为50 mg/L时对旱柳具有较大影响效应，本研究以旱柳叶片为实验材料从转录组水平着手，探究旱柳在Cd胁迫条件下的耐受性反应情况，为旱柳修复土壤Cd污染提供理论指导。

1 材料与方法 1.1 材料

本研究所用的旱柳枝条取样自昆明市北京路延长线的一株旱柳。一年生旱柳枝条被截取长度为25 cm继而用改良后霍格兰德营养液^[12]进行水培。水培后的旱柳随机分成3组（每组5盆，每盆3棵），1组为对照组，其余2组为试验组。Cd（NO₃）₂作为唯一的Cd来源对试验组进行胁迫处理，使营养液中Cd浓度分别为2.5 mg/L和50 mg/L，营养液每3 d换一次。研究材料为3组的第0、1、7与30天并处于同一位置的叶片，采样方法使用随机混采法；所有试验叶片采集完成后液氮急速冷冻后放置于-80℃冰箱内保存备用。对照组用CK表示，2.5 mg/L处理组用L表示，50 mg/L处理组用H表示；样品编号如下：第0天为CK0，第1天分别为CK1、L1、H1，第7天分别为CK7、L7、H7，第30天分别为CK30、L30、H30。

1.2 方法 1.2.1 总RNA提取及样品检测

总RNA提取使用OMEGA公司的E.Z.N.ATM Plant RNA Kit试剂盒进行提取，提取方法严格按照试剂盒说明书进行。并将获得的总RNA用琼脂糖凝胶电泳分析RNA的降解程度及是否存在污染，用液相色谱Agilent 2100和OD_260/280的比值检测RNA的纯度。旱柳总RNA的提取、质量检测及无参考基因组序列的转录组高通量测序分析均由深圳华大基因生物科技有限公司进行。

1.2.2 文库构建及转录组测序

样品检测合格后，加入Fragmentation buffer将mRNA打断成短片段，用六碱基随机引物（Random hexamers）以mRNA为模板合成一链cDNA，然后加入缓冲液、dNTPs和DNA polymerase I合成二链cDNA，双链cDNA再利用AMPure XP beads纯化。纯化后双链cDNA末端修复、加A尾并连接测序接头，片段大小选择是用AMPure XP bead，最终的cDNA文库通过PCR富集得到。cDNA文库质量用Agilent 2100检测，再通过Illumina Hi-Seq2500平台完成RNA-Seq测序分析。

1.2.3 De nove组装

获得的原始数据Raw reads经过滤得到Clean reads。使用Trinity对clean reads组装得到contigs，之后再用Tgicl对转录组样本进行两次聚类去冗余得到最终Unigene。Unigene分为两部分，一部分是clusters，是进一步去冗余后的结果；其余的是singletons，指没有聚类、单独的Unigene。每个样品测序产出不少于6 Gb Clean data，Q30碱基百分比达到85%。

1.2.4 生物信息学分析

使用BLAST对Unigene序列与Nr、Nt、KOG、KEGG、GO、Swiss-Prot及InterPro数据库进行比对（Evalue < 1e-5），获取与旱柳Unigene具有最高序列相似性的蛋白，得到Unigene的蛋白功能注释信息。

本研究从两个角度分析旱柳在Cd胁迫后的差异基因表达情况：（1）同一时间不同浓度Cd胁迫的差异基因表达情况的对比分析，命名为L1-CK1（即L1和CK1两个样品的基因表达对比，下同）、H1-CK1、L7-H7、L30-CK30、H30-CK30；（2）同一浓度不同时间的差异基因表达情况对比分析，分为2.5 mg/L组和50 mg/L组，命名为L1-CK0、H1-CK0、L7-L1、H7-H1、L30-L7、H30-H7。

根据差异表达分析结果，使用Blast2GO和WEGO软件对差异表达的Unigene进行GO功能分类^[13]。将差异表达的Unigene注释到KEGG数据库中得到注释信息，之后进行pathway分析^[14]。得到2.5 mg/L和50 mg/L浓度Cd胁迫后旱柳的生理生化响应结果。

1.2.5 差异表达基因q-PCR验证

利用primer premier 5.0对Cd胁迫响应的相关基因和内参基因肌动蛋白基因（ACT）进行q-PCR引物设计。在BIO-RAD公司的C1000^TM荧光定量PCR仪进行real-time PCR反应，每个样品重复3次，反应体系20 μL。以内参基因ACT为对照，利用公式2^-△△Ct计算其相对表达量。

2 结果 2.1 总RNA质检结果

10个样品中除CK7外9个样品的RNA质量较好，能满足建库测序要求，可用于RNA转录组测序，且总量满足1次但不足2次建库需要。CK7因RNA提取质量不佳，所以没有进行转录组测序。

2.2 转录组测序与De nove组装

9个样品的raw reads在72.22 Mb到85.59 Mb之间，过滤后总clean reads在65.17 Mb到67.15 Mb之间；Q20、Q30以及过滤前后reads的比率分别是95.04%到97.38%，86.40%到90.91%，72.75%到89.4%（表 1），这表明低质量的碱基比率低，测序质量好，可以用于de novo组装。用Trinity对其进行组装，所有样品的总reads均超过10万，总长度超过7千万，平均长度超过600，N50超过1 000，GC含量42%左右。之后使用Tgicl对转录本进行聚类去冗余得到Unigene，聚类后的Unigene质量（表 2），10个样品的Unigene数为53 705到62 933个，总长度从43 247 382 nt到50 853 969 nt，N50值都在1 200之上，这表明组装质量高，可用于后续分析。

2.3 对Unigene的功能注释

得到的Unigene有83.93%在7个数据库（Nr、Nt、KOG、KEGG、GO、Swiss-Prot及InterPro）中任意一个注释到。Nr数据库进行比对筛选，筛选条件为E值小于10^-5，结果显示：有74 919条Unigenes在Nr数据库中找到同源蛋白序列，占总Unigenes数的73.02%。Nt数据库注释到的Unigenes最多，占总Unigenes的81.12%，但被全部数据库注释到的Unigenes仅为29.04%，说明有许多旱柳Unigenes的功能还没有完全明确。

2.4 Unigene GO功能分类

根据Nr数据库注释的信息结果显示，有51 494条Unigene映射到GO不同的功节点（Term）上，根据生物过程（Biological process）、细胞组成（Cellular component）、分子功能（Molecular function）进行分类。由于经常有同一个转录本映射到不同节点现象，所以在生物过程中有106 517条，占42.69%，参与代谢过程（Metabolic process）（22 955，21.54%）、细胞过程（Cellular process）（21 034，19.75%）的Unigene最多；在细胞组成中有92 900条，占37.24%，参与膜（Membrane）（15 182，16.34%）、细胞器（Organelle）（13 010，14.01%）、细胞（Cell）（16 820，18.11%）和细胞部分（Cell part）（16 671，17.95%）的Unigene最多；在分子功能中有50 054条，占20.06%，参与催化活性（Catalytic activity）（20 318，40.59%）和结合蛋白（Binding）（20 533，41.02%）的Unigene最多（图 1）。

2.5 Unigene的KEGG注释和代谢通路分析

将旱柳差异表达基因在KEGG数据库中进行比对分析显示，53 966个Unigene分别注释到KEGG数据库的6个一级层级和21个二级层级中，共包含137条通路。其中在同一个pathway通路中10个以上差异表达基因主要涉及细胞过程、环境信息处理、遗传信息处理、人类疾病、代谢以及有机系统6个方面。这6个方面中代谢类的涉及最多，包括碳水化合物、氨基酸、多糖生物合成、能量、次生物质等；涉及人类疾病的最少，只有2个，分别是内毒素和代谢疾病（Endocine and metabolic diseases）和耐药性:抗菌素（Drug resistance : Antimicrobial）。在旱柳Cd胁迫后KEGG分析中发现涉及人类疾病的基因变化，这也表明Cd对旱柳中有关人疾病的内毒素和抗菌性的代谢有影响。同时，pathway富集结果也表明，前述的25个Cd胁迫响应候选基因主要参与代谢途径（Metabolic pathway）（10541，19.53%）、次生代谢产物的合成（Synthesis of secondary metabo-lites）（5473，10.14%）、植物MAPK信号通路（Plant MAPK signaling pathway）（1592，2.95%）、植物病原互作（Plant pathogen interaction）（1761，3.26%）、植物激素信号转导等通路（Plant hormone signal transduction）（1862，3.45%）（表 3）。由此推测旱柳可能主要是通过这些代谢通路调控Cd抗逆性过程。

2.6 差异表达基因分析

FDR（False discovery rate）值越小，表明基因表达差异越大，差异越显著。本研究在差异基因时设置的阈值为Fold Change≥2.00，FDR≤0.001，Q-value < 0.05。

由表 4可知，相同时间不同浓度的差异表达基因（Differentially expressed genes，DEG）分析显示：第1天组（L1-CK1和H1-CK1）中的差异基因数为8 550个（3 660个上调和4 890个下调）；第7天组（L7-H7）的差异基因数为3 444个（1 307个上调和2 137个下调）；第30天组（L30-CK30和H30-CK30）中的差异基因数为11 428个（5 055个上调和6 373个下调）。相同浓度不同时间的DEG分析显示：2.5 mg/L组（L1-CK0，L7-L1，L30-L7）中的差异基因数为26 623个（14 531个上调和12 092个下调；50 mg/L组（H1-CK0，H7-H1，H30-H7）中的总差异基因数为32 154个（18 318个上调和13 836个下调），这说明50 mg/L浓度的Cd胁迫比2.5 mg/L的对旱柳的影响更加明显。此外，在L1-CK1、H1-CK1、L30-CK30和H30-CK30这4组中有共同的170个DEG（其中46个共同上调和29个共同下调），它们的基因表达量差异达2倍（|Log₂ration|≥1）以上。170个DEG中又有25个其表达量差异是大于1 000倍（|Log₂ration|≥10）以上的，其中19个是基因表达上调的，6个是基因表达下调的（表 5）。这25个基因表达变化如此之大，说明它们在旱柳受Cd胁迫耐受过程中等起到了重要作用，可以作为旱柳Cd胁迫响应的重要候选基因。

2.7 差异基因的q-PCR

黄酮类化合物生物合成和油菜素内酯合成的关键基因：FLS，F3H和ROT3，这3个基因在旱柳Cd胁迫响应代谢通路中表达量显著上升，钙依赖性蛋白激酶2基因表达量显著下调，因此我们选择这4个基因进行q-PCR验证转录组测序分析的准确性和可靠性，肌动蛋白基因（ACT）因其具有稳定性而被选择为内参基因^[15]。图 2-6的q-PCR结果显示：5个差异表达基因在转录组测序中的表达变化与荧光定量结果基本一致，说明转录组测序分析的结果真实可靠。

3 讨论

旱柳是中国土生土长的高大乔木，具有易扦插、易成活、生长快、枝叶发达和抗逆性强等特点，是一种非常适合用修复土壤重金属污染的物种。为研究旱柳响应不同浓度Cd胁迫的基因表达情况，选择Cd胁迫浓度分别为2.5 mg/L和50 mg/L的旱柳叶片作为转录组测序样品，通过旱柳基因表达的差异来研究旱柳响应Cd胁迫的分子机制。转录组数据进行过滤、拼接、组装后共获得102 595个基因。

近年来对植物响应重金属胁迫后的基因表达研究已取得一系列成果，鉴定和推测出不少关于重金属响应的基因，包括ABC转运家族、锌和锰转运蛋白家族、重金属ATP酶（HMAs）家族、阳离子扩散促进者（CDF）家族、植物螯合肽合成基因（PCS）、金属硫蛋白基因（MT）和金属忍耐蛋白基因（MTP）等^[16]。本研究发现在旱柳响应Cd胁迫的过程中，金属硫蛋白基因、ABC转运家族、锌和锰转运蛋白家族的基因表达会随着Cd浓度的增高和胁迫时间的增加而有较大的变化，L7-L1、L30-L7、H7-H1、H30-H7组中金属硫蛋白共9个基因表达上调、11个下调；ABC转运家族共90个上调、88个下调，其中H30-H7中32个上调，50个下调；锌和锰转运蛋白家族有28个下调，8个上调，上面列的其它重金属响应基因的差异表达在本研究没有发现。这些结果表明旱柳在应对Cd胁迫的过程中金属硫蛋白基因、ABC转运家族及锌和锰转运蛋白起了很大的作用，也说明Cd胁迫会影响旱柳对其他金属离子的吸收，锌和锰转运蛋白表达量的变化可能是因为Cd在生命体内会置换锌，导致生命体需锌的蛋白质发生“饥饿”，为了生命正常活动的进行，旱柳加快了锌的转运过程，从而使锌、锰蛋白转运蛋白家族的表达发生变化。ABC转运家族基因表达量发生明显变化，说明旱柳在Cd胁迫响应中细胞物质跨膜转运主要是通过ABC转运家族蛋白来实现的，也说明ABC转运家族在旱柳细胞运送重金属相关物质中起重要作用。

黄酮类化合物是植物体内一类在抗菌、抗逆性上具有重要的作用的化合物，是一种良好的抗氧化剂，能够清除生物膜周围和细胞内的H₂O₂^[17-18]。有研究表明，黄酮化合物的代谢通路在水稻响应Cu²⁺和Cd²⁺的胁迫以及组培旱柳苗响应Cd²⁺中起了积极作用^[19-20]。孙爱清等^[21]利用Solexa高通量测序对花生干旱响应基因表达基因分析，发现有9个类黄酮代谢相关基因在干旱胁迫下显著表达；宋中邦^[22]通过烟草基因组数据库中搜索获得两个类黄酮合成酶基因，即NtFLS1和NtFLS2，二者表达模式一致，均在叶片中高水平表达。本研究通过KEGG代谢通路图进行分析发现，旱柳黄酮类化合物的生物合成通路在样品中均有大量的差异基因上调，L1-CK1、H1-CK1、L30-CK30、H30-CK30组中一共170个差异基因，其中有46个共同上调基因，29个共同下调基因，而这46个共同上调基因中包含了FLS、F3H基因。这说明旱柳在应对Cd胁迫时，黄酮类化合物起了重要作用。

同时，本研究发现油菜素内酯合成通路的3，6-脱氧油菜素淄酮酶（ROT3）在除L1-CK1外的其他样品中都是上调。肖瑞雪等^[23]研究发现油菜素内酯是一种被誉为第六激素的天然激素，其生理活性甚至超过现在5种激素，能提高植物的抗逆性。靳开川等^[24]结合遗传学、基因与蛋白质组学、细胞生物学等多学科方法和手段，发现油菜素内酯在植物的抗逆性（抗干旱、高盐、高温、低温、重金属）过程中的起重要作用。王喆等^[25]研究表明油菜素内酯可能会与膜蛋白结合，然后通过减轻重金属胁迫来提高代谢活性，从而降低植物对重金属的摄取能力^[26]。本研究在H1-CK1、L30-CK30、H30-CK30中油菜素内酯合成通路的3，6-脱氧油菜素淄酮酶（ROT3）上调是由于Cd²⁺在旱柳叶片中积累，旱柳为减轻Cd²⁺带来的损害，而促进了油菜素内酯的合成。

本研究对旱柳受Cd胁迫的转录组进行了差异基因的分析，找到了受胁迫的相关基因以及与之相关的代谢通路，但这些基因是如何互相作用调节旱柳Cd胁迫，还需要在分子方面进行更加深入的分析。

4 结论

旱柳在受重金属Cd胁迫后展现出一定程度的胁迫抗性，同时基因表达与代谢过程发生了重大的变化。（1）镉胁迫后不同浓度不同时间的旱柳转录组测序共获得102 595个Unigenes；（2）筛选得到与Cd胁迫响应密切相关的基因25个，其中金属硫蛋白、ABC转运蛋白、锌和锰转运蛋白基因在旱柳植物抗逆性过程中发挥重要作用；（3）Cd胁迫的GO条目主要集中在代谢过程、细胞过程、膜、细胞器、细胞、细胞部分、催化活化和结合蛋白上；（4）油菜素内酯合成通路的3，6-脱氧油菜素淄酮酶（ROT3）和黄酮类化合物合成通路的黄酮醇合成酶（FLS）、黄烷酮-3-羟化酶（F3H）均明显上调。