镉胁迫下桑树幼苗的转录组分析

王晓红，韩世玉，孙运朋，张芳，罗泽虎，陈彪

（贵州省农业科学院蚕业研究所，贵阳 550025）

0 引言

来自自然、农业和工业来源的重金属对耕地的污染已成为世界性的公共卫生问题[1]。镉(Cd)是一种重要的重金属污染物，被认为是最具毒性的重金属之一。镉污染土壤的修复是当前农业面临的一个重要问题。用于这一目的的不同的物理和化学方法都存在成本高、劳动强度大、土壤性质改变和土壤原生微生物区系干扰等严重的局限性。相比之下，植物修复技术是一种环境友好、经济高效、人们接受度较高的技术[2]。植物可以提取或稳定土壤中的Cd，在污染土壤修复中得到广泛应用，一些Cd超富集植物可以在不影响其正常生长的情况下从污染土壤中吸收和提取Cd。然而，多数Cd超富集植物是草本植物，其生物量小、根系浅、生长缓慢。相比之下，人们逐渐提出使用具有大生物量和多分枝生根系统的木本植物来修复Cd污染土壤[3]。

桑树是一个具有较高经济和生态价值的木本植物，具有适应性强、生长快、根系深、地上生物量大、耐剪伐等优点，应用于石漠化、盐碱地、戈壁沙漠等多种植物修复系统中[4]。同时，桑树还被证明能够耐受Cd、Cr、Pb、Co、Cu等重金属，具有修复重金属污染土壤的极好潜力[5-7]。此外，利用桑树进行Cd污染土壤修复也是一种经济可行的做法，因为桑叶是家蚕的饲料，可用于蚕业发展。目前，多数研究集中在桑树对Cd胁迫生理反应或部分耐Cd 基因的功能研究[8-9]。全面阐明桑树中Cd 吸收、积累、转运和解毒的机制对进一步提高桑树在Cd污染土壤中的植物修复价值具有重要意义。

Cd 通过抑制光合作用、呼吸和氮代谢，扰乱水和矿物质养分的吸收和分配，阻碍植物的生长和发育[10-11]。为减轻Cd 的毒害作用，植物已经进化出一系列的解毒策略，包括螯合作用、细胞壁固定、区室化作用和抗氧化系统解除重金属毒害[12]。这些解毒机制的阐明需要广泛的转录、蛋白质组学和代谢研究。转录组测序技术是从分子水平阐明其机制的有效方法，该技术具有效率高、通量大、运行时间快、准确性高的特点，常用于发现和表征基因、分析功能基因变异、识别和量化转录本[13]。利用这种方法已鉴定出柳树[3]、白菜[14]和萝卜[15]等的Cd胁迫相关基因。此外，比较转录组学分析已用于揭示植物中Cd积累的机制[14,16]。这些研究极大地提高了人们对Cd耐受和积累的分子调控机制的理解。本研究将桑树幼苗在Cd 处理12、24、48 h后，采用Illumina HiSeq 2000平台进行测序及分析，旨在了解桑树对Cd的吸收、转运、隔离和耐受机制，并在植物修复研究中发挥作用。

1 材料与方法

1.1 试验材料及处理

试验于2021 年3—8 月在贵州省蚕业研究所植物生理实验室进行。试验桑树品种为‘桂桑优12’，购自广西蚕业推广总站。挑选完整、饱满的种子，75%酒精灭菌60 s，无菌水冲洗5 次。将灭菌后的桑树种子接种到MS 培养基上，16 h/8 h 光周期、24℃条件下萌发培养。待胚根长至2 cm左右时取出幼苗，接种到添加50 μmol/L CdCl2的MS培养基中，每瓶接种10棵幼苗；同时，以未添加CdCl2的幼苗为对照。分别在处理12、24、48 h取样，迅速用无菌水将幼苗冲洗干净后放入液氮中速冻，保存于-80℃备用。

1.2 RNA提取、文库构建及转录组测序

分别收集Cd 处理和对照组12、24、48 h 的桑树幼苗提取RNA用于转录组测序。总RNA提取、cDNA文库构建和转录组测序由北京百迈克生物公司完成。主要流程如下：用带有Oligo(dT)的磁珠富集真核生物mRNA；加入Fragmentation Buffer 将mRNA 进行随机打断；以mRNA 为模板，用六碱基随机引物合成第一条cDNA 链，然后加入缓冲液、dNTPs、RNase H、DNA polymerase I 合成第二条cDNA 链，利用AMPure XP beads纯化cDNA；纯化的双链cDNA再进行末端修复、加A 尾并连接测序接头，然后用AMPure XP beads 进行片段大小选择；最后通过PCR 富集得到cDNA 文库。对检测合格的文库进行测序。

1.3 转录组从头组装和功能注释

fastq格式的原始读取首先通过内部的perl脚本进行处理，从原始数据中删除包含接头的序列，超过10%未知核苷酸的序列以及低质量的序列，以获得干净序列。使用Trinity软件对干净读数进行组装。随后，基于数据库NR(NCBI non-redundant protein sequences)、Swiss-Prot、GO(Gene Ontology)、KOG/COG/eggNOG(Clusters of Orthologous Groups of proteins)和KEGG(Kyoto Encyclopedia of Genes and genome)，使用BLASTx对转录本进行注释，期望值E-value≤10-5。使用HMMER 将预测的氨基酸序列与Pfam 数据库进行比较，期望值E-value≤10-10，获得注释信息。

1.4 Unigenes的差异表达分析

为了确定Cd 胁迫下的差异表达基因，用Bowtie软件将干净序列映射到组装的转录组上。使用RSEM估计基因表达水平，unigenes 的表达水平用FPKM 归一化。使用DEGSeq 鉴定Cd 处理样本和对照组之间的差异基因。2 组之间的错误发现率(FDR)≤0.01 和|log2 Fold Change|≥1 被作为差异表达的阈值。随后，使用topGO 和KOBAS 分别对差异表达的基因进行GO富集分析和KEGG通路富集分析。

1.5 qRT-PCR分析

采用实时荧光定量PCR(qRT-PCR)技术对转录组数据进行验证，随机选择4个DEGs进行qRT-PCR。以桑树MaACT基因作为内参基因[17]，根据各基因序列设计引物（表1），在Bio-Rad实时荧光定量PCR仪上检测差异基因表达水平。根据2×RealStar Green Fast Mixture（GenStar，北京）试剂盒的说明书配制qRT-PCR反应体系：2×RealStar Green Fast Mixture 10 μL、模板cDNA 1 μL、正/反引物各0.5 μL(10 μmol/L)，ddH2O 8 μL。PCR程序：94℃预变性2 min；94℃变性15 s，60℃退火30 s，循环40次。采用2ΔΔCt方法计算基因相对表达量。

表1 研究所用qRT-PCR引物

2 结果与分析

2.1 桑树幼苗转录组测序、组装与注释结果统计

为全面分析Cd胁迫下桑树幼苗的转录组，将对照组T1N (12 h)、T2N (24 h)、T3N (48 h)和Cd 处理组T1P(12 h)、T2P(24 h)、T3P(48 h)的6个样品进行转录组测序，6个样本分别产生了37398420～49460056条原始序列，去除低质量序列后，获得了36139366～47208418 条干净序列，占原始序列的94.99%以上，Q30在92.15%～93.44%之间（表2），说明转录组测序质量较高，可以进行下一步的数据整理和分析。使用Trinity软件对干净序列进行从头组装，结果得到39758个unigenes，长度在201～13312 bp之间。这些unigenes的平均长度和N50 分别为1134.51、1968 bp，平均GC含量为41.72%（表3）。

表2 Illumina转录组测序总结

表3 非冗余转录本统计

2.2 Unigenes功能注释与分类

组装后，通过查询7个公共数据库预测基因功能，39758 个unigenes 得到注释（表4）。其中，有26280 个unigenes(66.10%)在Nr数据库中获得了注释，18589个unigenes(46.76%)在Swiss-prot 数据库中获得了注释，18606个unigenes(46.8%)在GO数据库中获得了注释，只有7201个unigenes(18.11%)在KEGG数据库中获得了注释。

表4 桑树unigenes在公共数据库中的注释结果

在KOG分析中，有13734个单基因被分成25组用于功能预测和分类（图1）。其中，仅一般功能预测(3831)包含的基因数量最多，其次是翻译后修饰、蛋白质转换、伴侣(1465)、信号转导机制(1344)、碳水化合物运输和代谢(833)和转录(762)。此外，18606 个单基因被注释为一个或多个GO 术语，并被分为3 个类别和57 个组（图2）。对于生物过程(biological process，BP)类别，细胞过程(12445)代表最大的组，其次是代谢过程(11271)、生物调控(5578)和刺激反应(4620)。在分子功能(molecular function，MF)类别中，结合(12335)、催化活性(10218)、转运活性(1229)和转录调节活性(1190)是最常见的组。在细胞成分(cellular component，CC)类别中，前5 名分别是细胞组分(14509)、细胞器(9053)、膜(5491)、细胞器的部分(4674)和膜部分(4125)。

图2 桑树转录本GO注释分类

2.3 差异表达基因分析

分别对Cd 胁迫和对照处理12 h (T1N_vs_T1P)、24 h (T2N_vs_T2P)、48 h (T3N_vs_T3P) 样本的unigenes 进行差异表达分析（图3）。在T1N_vs_T1P组，共有6319 个差异表达基因，其中，4411 个基因上调表达，1908 个基因下调表达；7737 个基因在T2N_vs_T2P 中存在差异，其中，2129 个基因上调表达，5608 个基因下调表达；而在T3N_vs_T3P 中，共有6561 个差异表达基因，其中，3389 个基因上调表达，3172 个基因下调表达（图3A），表明桑树幼苗在Cd 胁迫不同时期转录水平的动态变化过程。值得注意的是，有534 个基因在不同时期均有差异（图3B），其中有4 个基因(TRINITY_DN10246_c0_g4、TRINITY_DN9070_c0_g1、TRINITY_DN12871_c0_g1、TRINITY_DN7347_c1_g2)注释为ABC 转运蛋白家族成员，其在Cd 胁迫过程中上调表达，表明它们可能在整个Cd胁迫过程中发挥作用。

图3 桑树各处理组差异表达基因

2.4 差异表达基因GO富集分析

为了进一步研究差异表达基因在Cd 胁迫下的生物学功能，进行了GO 富集分析。Cd 胁迫12、24、48 h后，每组DEGs可分为生物过程、细胞成分和分子功能3类，选取富集显著程度最高的10个GO条目分析（图4）。在生物过程类别，Cd 胁迫不同时间后，差异表达基因主要集中于能量代谢、DNA代谢和应激反应相关的组中，如ATP代谢过程(ATP metabolic process)、ADP代谢过程(ADP metabolic process)、DNA 重组(DNA recombination)、DNA 代谢过程(DNA metabolic process)和急性炎症反应(acute inflammatory response)等；在细胞组分类别，MHC 蛋白复合体(MHC protein complex)、NADPH 脱氢酶复合体[NAD(P)H dehydrogenase complex (plastoquinone)]、细胞(cell body)和细胞壁(cell wall)等是几个富集的组；在分子功能类别，主要富集的组与能量转换、DNA 复制与转录过程酶活性有关，例如，NADH 脱氢酶活性(NADH dehydrogenase activity)、RNA 聚合酶II 转录因子活性(RNA polymerase II transcription factor activity)、DNA结合转录因子活性(DNA binding transcription factor activity)和DNA 聚合酶活性(DNA polymerase activity)等；其他富集的组还包括酰基载体蛋白质合成酶活性(3-oxoacyl-[acyl-carrier-protein] synthase activity)和7-脱氧马钱酸葡糖基转移酶活性(7-deoxyloganetic acid glucosyltransferase activity)等。这些GO 富集结果表明，Cd胁迫引起桑树幼苗能量代谢、酶催化活性，以及DNA复制和转录相关途径基因的差异表达，这些过程可能在桑树幼苗抵抗Cd胁迫中起到重要的作用。

2.5 差异表达基因KEGG通路富集分析

KEGG数据库是系统分析基因产物在细胞中的代谢途径以及这些基因产物功能的数据库。以q＜0.05为标准，对KEGG中每个通路应用超几何检验进行富集分析，找出差异表达基因中显著性富集的通路。结果如图5 所示。在Cd 胁迫不同时间，谷胱甘肽代谢(glutathione metabolism)、倍半萜和三萜生物合成(sesquiterpenoid and triterpenoid biosynthesis)、光合作用-天线蛋白(photosynthesis-antenna proteins)、苯丙素的生物合成(phenylpropanoid biosynthesis)、光合作用(photosynthesis)、戊糖和葡萄糖醛酸转换(pentose and glucuronate interconversions)、异黄酮生物合成(isoflavonoid biosynthesis)、半乳糖代谢(galactose metabolism)、氧化磷酸化(oxidative phosphorylation)、类黄酮生物合成(flavonoid biosynthesis)和DNA 复制(DNA replication)等11 条代谢通路差异表达基因富集程度较高，这些代谢通路可能与桑树中Cd的代谢相关。

图5 差异表达基因KEGG富集分析

2.6 差异转录因子分析

许多类型的转录因子广泛参与非生物和生物胁迫反应，并在各种植物发育过程和胁迫反应中发挥重要的调节作用。与对照相比，在Cd 胁迫不同时期，桑树幼苗中共鉴定了148 个差异表达的转录因子，代表了36个不同的家族（图6）。其中，bHLH、MYB、B3、NAC、MYB-related和WRKY家族是拥有差异基因较多的转录因子家族。

图6 差异表达转录因子分类

2.7 qRT-PCR验证转录组数据

为进一步验证转录组的比较结果，随机选择4 个差异表达基因进行qRT-PCR 检测（图7），这些基因的表达模式与转录组检测到的基因表达变化是一致的，表明转录组谱数据是可靠的。

3 结论与讨论

Cd是一种常见的环境污染物，对人类健康构成潜在的威胁。如何有效地治理Cd 污染，降低Cd 对生物体的潜在风险，已成为一个重要的社会问题[18]。虽然已经开展了各种研究来探讨植物对Cd反应的生理、遗传和分子基础，但桑树中Cd 吸收、积累和运输的分子调控机制在很大程度上仍未被认识。作为一种重要的生态经济树种，桑树被认为是解析重金属调控网络的最重要的木本植物之一[9,19]。

最近，RNA-Seq方法已被广泛应用于许多重要的植物的从头转录组测序组装，如萝卜[15]、玉米[20]、柳树[3]和水稻[21]。在本研究中，通过从头转录组测序和组装，从Cd 胁迫下的桑树幼苗中获得了39758 个非冗余转录本。基于Pfam、Nr、GO 和KEGG 数据库，通过BLAST和功能生物信息学分析，成功地注释了全部转录本序列。本研究中的转录组数据将为进一步研究桑树中Cd胁迫响应基因和功能提供有价值的信息。

转录组测序分析使人们能够全面了解桑树响应Cd 胁迫有关的基因和过程。在Cd 胁迫下，桑树共检测到15882 个差异表达基因。其中，分别有6319、7737、6561 个差异表达基因在胁迫12、24、48 h 被鉴定，表明Cd胁迫导致桑树大量基因进行了重编程。在这些差异表达基因中有4 个转录本被注释为ABC 转运蛋白家族成员，它们可能对桑树Cd耐受性具有重要作用。ABC转运蛋白广泛存在于动植物中，负责各类底物在生物膜上的运输，包括重金属、激素和羧酸盐等[22]。已有的报道表明，ABC转运蛋白对Cd吸收转运具有重要的调节作用，AtABCC3缺陷型拟南芥幼苗对不同Cd 浓度的敏感性增加，过度表达AtABCC3基因的幼苗对Cd的耐受性增加[23]。

Cd 通过多种毒害作用抑制植物生长。已有研究表明，Cd会损害光合机构，改变光合过程，最终导致植物生长抑制[24]。Cd 还会改变氧化酶的活性来诱导氧化应激，从而对植物产生不利影响[25]。相反，植物会通过增强防御反应，如金属螯合、液泡隔离、转运蛋白调节以及强化抗氧化机制来抵抗重金属威胁[26]。本研究中，KEGG富集分析表明，具有差异表达的基因在谷胱甘肽代谢、次生代谢、光合作用、能量代谢和DNA复制等相关代谢通路富集程度较高。这些结果与构树类似，在Cd 胁迫下，构树差异表达基因在多个KEGG 途径富集，包括苯丙烷生物合成、氧化磷酸化和类黄酮生物合成等[16]。表明植物在抵抗Cd 胁迫时具有相似的机制。

已有的研究表明，许多转录因子家族在Cd胁迫中起着重要作用。Hu 等[27]报道水稻的一个转录因子OsMYB45参与水稻的Cd胁迫反应，OsMYB45突变导致对Cd 处理的超敏反应，而OsMYB45过度表达补充了突变体表型。Yang等[28]报道ThWRKY7可能通过激活ThVHAc1基因表达来提高转基因拟南芥对Cd的抗性。此外，萝卜中bZIP、MYB、ERF、C2H2、NAC 家族被证明在介导Cd 胁迫反应中起着重要的调节作用[15]。bHLH、GRAS、WRKY、bZIP 家族在构树Cd 胁迫中具有复杂的转录调控作用[16]。在本研究中，大量转录因子在Cd 胁迫下差异表达，包括bHLH、MYB、NAC、WRKY等转录因子家族。表明不同转录因子家族的调控可能会触发植物对Cd 胁迫的响应，也说明Cd信号调控的复杂性。

总之，本研究从桑树幼苗转录组中分离到39758个非冗余转录本，其中15882个基因在Cd胁迫下受到差异调控。KEGG 富集分析表明，这些差异表达基因参与了许多代谢、生物合成和非生物胁迫响应途径的过程。同时，还鉴定了一些ABC转运蛋白和大量转录因子，这些发现将为研究桑树Cd胁迫反应的分子机制提供有价值的信息。