泡桐根系应答铅锌矿胁迫的转录组分析

2020-08-31韩良泽陈明利陈永华贾惠雁欧阳林男

中南林业科技大学学报 2020年8期

韩良泽，陈明利，陈永华，陈昊，贾惠雁，欧阳林男

（1. 中南林业科技大学 a. 环境科学与工程学院；b. 林学院，湖南长沙 410004；2. 国家林业和草原局桉树研究开发中心，广东湛江 524000）

植物对重金属的耐性往往是多基因控制的复杂性状，响应重金属胁迫的基因通常涉及到金属离子转运、螯合作用、离子区隔化、抗氧化系统等诸多方面[1]。植物应答重金属胁迫研究有助于揭示植物对重金属的应答机理，从而可以采取相应的措施去降低重金属对植物的毒害。转录组学作为基因组与蛋白质功能组的重要桥梁，其研究已经成为揭示生物生长发育规律以及寻找逆境胁迫机制等方面的最佳手段。转录组所测定的对象是某些细胞在特定状态下的mRNA，而泡桐尚未有公布的参考基因组，因此，对于无参基因的转录组，可以通过高通量测序技术对某种物种某一器官在特定状态下进行转录组测序，组装得到基因，再进行功能注释。这种方法相比较常规的分子手段，具有高效、价格低、覆盖面广等优点[2-3]。因此可以用此技术去发掘特定状态下的关键基因。

白花泡桐隶属于玄参科泡桐属，是一种喜光、耐寒性不强、对瘠薄土壤有较强适应性的速生树种。课题组前期在野外调查中发现泡桐在铅锌矿中可以存活[4]，说明泡桐有较强的适应能力。在矿上生长的泡桐，其体内的重金属含量明显高于背景区[5]，可当作金属富集植物，用于尾矿的生态修复。泡桐根系作为与污染源直接接触的部位，其对重金属的吸收、积累及耐性机制等方面已有研究[6]，但是，在重金属胁迫下，泡桐根转录组调节机制及相关基因表达尚未有明确研究，因此，以泡桐为植物材料，以铅锌矿、红壤为基质，进行盆栽实验，并对所得的泡桐根进行转录组测序，分析泡桐根部响应重金属胁迫的相关差异表达基因，以期进一步研究受重金属胁迫差异表达的基因与自然矿床或开采矿山污染土壤的关系，并为泡桐应用于矿山生态修复提供基础理论依据。

1 材料与方法

1.1 实验材料

供试材料白花泡桐购置于河南某苗圃，选取粗细一致的泡桐根置于红壤中进行土培生芽，待幼苗长至株高10 cm 左右时选取高度一致的幼苗移栽。铅锌矿取自于郴州市资兴铅锌矿尾库，红壤取自于中南林业科技大学生态园0 ～20 cm 表层土壤，钙镁磷肥购置于湖南长沙花卉市场。

1.2 实验设计

实验在中南林业科技大学苗圃内进行，将铅锌矿、红壤按照质量比分别进行搅拌成2 个实验组R（红壤土+0.3 kg 钙镁磷肥）、T（铅锌矿渣+0.3 kg 钙镁磷肥），每盆储量10 kg，每盆1 株，每个处理3 次重复，定期浇水，除虫，实验周期为1 a。土壤重金属含量见表1。

表1 盆栽土壤重金属质量分数†Table 1 Heavy metal content in potted soil

1.3 实验方法

1.3.1 RNA 提取与鉴定

将种植一年后的泡桐根部洗净，经液氮速冻后，用于总RNA 的提取。使用Qubit®2.0 Flrometer（加利福尼亚州Life Technologies）中的Qubit®RNA 分析试剂盒检测RNA 浓度。使用纳米光度计®分光光度计（Implen，CA，USA）检测RNA 纯度，使用Agilent Bioanalyzer 2100 系统（美国加利福尼亚州Agilent Technologies）的RNA nano 6000 分析试剂盒评估RNA 完整性。

1.3.2 文库构建及转录组测序

每个样品取总RNA1.5 μg，用带有Oligo(dT)的磁珠富集mRNA。随后加入Fragmentation buffer将mRNA 打断成短片段，以此mRNA 为模板，用六碱基随机引物（Random hexamers）合成cDNA 第一链，随后加入缓冲液、dNTPs、DNA polymerase I 和RNase H 合成cDNA 第二链，再用AMPure XP beads 纯化双链cDNA。将纯化的双链cDNA 进行末端修复、加A 尾并连接测序接头后，再用AMPure XP beads 选择用于测序的片段，并进行PCR 扩增，用AMPure XP beads 纯化PCR 产物，即获得最终的测序文库。测序由北京诺禾致源生物信息科技有限公司完成。

1.3.3 测序数据质量评估及转录本拼接

将基于测序平台所得到的原始数据进行内部Perl 脚本处理，其中包括去除有接头的reads、去除低质量的reads、去除含有超过10%的不确定碱基的reads。同时，通过Q20、Q30、GC 含量对clean reads 进行评估。使用Trinity[7]软件对质量合格的reads 进行组装，默认情况min_kmer_cov 设置为2，所有其他参数设置为默认值。

1.3.4 基因功能注释

将组装获得的基因序列与NCBI 官方非冗余蛋白质序列数据库（NCBI non-redundant protein sequence database，Nr）、NCBI 官方非冗余核苷酸序列数据库（NCBI non-redundant nucleotide sequence database，Nt）、KEGG 直系同源数据库（KEGG Orthology database，KO）、Swiss-Prot 数据库（Swiss-Prot）、蛋白质家族数据库（Protein family database，Pfam)、基因功能分类体系数据库（Gene Ontology database，GO）和真核生物直系同源序列数据库（euKaryoticOrthology Groups database，KOG)进行 Blast 比对（其中 Nr、Nt、Swiss-Prot 数据库的比对e值≤1e-5，KOG 比对e值≤1e-3，Go比对e值 =1e-6，KEGG 比对e值 =1e-10）。

1.3.5 SSR 监测及引物设计

使用 MISA(http://pgrc.ipk-gatersleben.de/misa/misa.html) 鉴定转录组中SSR 序列，并用Primer3(https://sourceforge.net/projects/primer3/获取软件信息)设计扩增SSR 序列的引物。

1.3.6 差异表达基因分析

使用DESeq对2个样品组进行差异表达分析，使用q值调整P值[8]。设置了 qvalue ＜ 0.005 ＆|log2（foldchange）|＞1 作为显著差异表达的阈值。分别用GOseq[9]和KOBAS[10]对获得的差异表达基因进行通路富集性分析。

1.3.7 差异表达基因的qRT-PCR 验证

挑选转录组测序得到的差异基因，以Actin 基因作为内标基因，引物对为：Actin-F（5′-TTGTGCTTAGTGGTGGGT-3′）和 Actin-R（5′-CTGGAAGGTGCTGAGGGAT-3′）。通过使用Primer6.0 软件设计引物并进行荧光定量分析。以 2 µgRNA-seq 测序的根总 RNA 为模板，合成cDNA 第一链。qRT-PCR 按照以下步骤进行扩增：95 ℃预变性 10 min，95 ℃变性 5 s，72 ℃退火 40 s，40 个循环。按照 2−ΔΔCT[11]计算基因相对表达量。

2 结果与分析

2.1 测序数据质量评估

提取泡桐根部RNA 后由Nova-seq 高通量测序，共获得366 701 960 条reads 序列，总碱基数为53.90 Gb，平均GC 含量达44.58%。序列质量评估值Q20、Q30 均达90%以上（表2），说明通过测序所得到原始序列数据较好，可用于后续的序列组装。用Trinity 软件对得到的clean reads进行组装，共获得311 822 条转录本以及126 061条基因序列，最短转录本长度为301 bp，最长转录本长度为96 533 bp，平均单条转录本长度为1 420 bp，N50 为2 195 bp，最短基因和最长基因均和转录本长度一样，平均单条基因为1 068 bp，基因序列在长度为200 ～500 bp 区间数量最多，约为总数的36%（表3）。

表2 样品测序数据统计†Table 2 Summary of sample sequencing data

表3 泡桐根部基因和转录本长度分布Table 3 Distribution of all-genes and transcripts length of Paulownia fortunei root

2.2 基因功能注释与分类

2.2.1 功能注释序列比较

将获得的基因序列在常用数据库（Nr，Nt，Pfam，KOG/COG，Swiss-Prot，KEGG，GO）中进行序列对比，最终获得96 133 条有注释信息的基因序列，注释率为76.25%（表4）。经Nr 数据库比对后发现，泡桐根部基因序列与芝麻Sesamum indicum基因序列同源性最高，达到20.9%，其后依次为栓皮栎Quercus variabilis（12.5%），紫花风铃Handroanthus impetiginosus（10.6%）和蓖麻Ricinus communis（5.3%），与其他物种同源的有32 416 条（50.7%）。注释后相似度达60%～80%、80%～95%、95%～100%的基因分别有23 337 条（36.5%）、18 222 条（28.5%）、4 795 条（7.5%）。

表4 泡桐根部Unigene 序列在各数据库的功能注释情况Table 4 Functional annotations of Unigene sequences in databases

2.2.2 GO 数据库注释与分类

GO 作为转录组数据研究中一套国际标准化的基因功能描述的分类系统，可用来描述基因及其蛋白产物的特性。将126 061 条基因与GO 数据库对比，共有68 840 条基因序列获得了注释，涉及3 个功能大类、55 个功能小类，从生物过程（Biological process）功能分类来看，主要集中在细胞过程（Cellular process）、代谢过程（Metabolic process）、单一有机体过程（Single-organism process）。从细胞组分（Cellular component）功能分类来看，细胞（cell）和细胞部分（Cell part）所占比例较高。从分子功能（Molecular function）分类来看，泡桐根主要涉及结合（Binding）、催化活性（Catalytic activity）、转运活性（Transporter activity）等，其中，仅结合和催化两大类约占总数的80%以上（图1）。

泡桐作为一种耐性植物，其根部涉及到的响应重金属胁迫的基因众多，如金属离子响应（Stress response to metal ion）543 条、金属酶活性调节（Metalloenzyme regulator activity）263 条、金属群结合（Metal cluster binding）248 条、氧含量响应（Responseto oxygen levels）854 条、谷胱甘肽酶转移活性（Glutathione transferase activity）762条、细胞壁增厚防御响应（Defense response by cell wall thickening）942 条等。多样的通路为泡桐的抗重金属育种提供理论基础。

2.2.3 KOG 数据库注释与分类

KOG 作为NCBI 针对真核生物的蛋白同源数据库，目前有4 852 个分类。将某个蛋白质与所有KOGs 条目中的蛋白质进行比对，可将其归入适当的KOG 簇，用于直系同源基因功能分类。泡桐根部基因涉及的功能类别较为全面，可将获得到的29 886 条基因序列分为25 类（图2）。其中蛋白质翻译后修饰与转运、分子伴侣基因最多（Posttranslational modification, protein turnover,chaperones，3 854 条）；其次为一般功能基因（General function prediction only，3 619 条）；翻译，核糖体结构和生物合成基因（Translation,ribosomal structure and biogenesis，3 281 条），而细胞运动基因（Cell motility，29 条）和胞外结构类基因（Extracellular structures，60 条）较少。

2.2.4 KEGG 数据库注释分析与分类

KEGG 是系统分析基因产物功能及其所在细胞代谢途径的数据库。据KEGG 数据库的注释信息进一步将泡桐根部基因序列进行pathway 注释，共有28 747条基因序列富集于不同的代谢途径中，其中包括细胞过程、环境信息处理、遗传信息处理、代谢、有机系统5 个代谢通路大类、环境适应、翻译、碳水化合物代谢等19 个代谢途径及核糖体、氨基酸的生物合成、碳代谢等131 个代谢分支。在二级分类中，富集于碳水化合物（Carbohydrate metabolism）（4 793 条）、翻译（Transcription）（4 456 条）途径的基因占大多数。在各分支中，富集于核糖体（Ribosome）、碳代谢（Carbon metabolism）、氨基酸生物合成（Biosynthesis of amino acids）的基因相对较多（表5）。以上这些代谢通路可能在泡桐根部生理代谢过程中具有重要作用。

图1 GO 注释分类Fig. 1 GO annotation classification

2.2.5 基因表达水平分析

本实验采用RSEM 软件[12]将 Trinity 拼接得到的转录组作为参考序列（ref），将每个样品的clean reads 比对回ref（mapping），发现各样品均有较高的注释率（表6）。随后将RSEM 对bowtie的比对结果进行统计，进一步得到了每个样品比对到每个基因上的readcount 数目，并对其进行FPKM 转换，结果显示，大部分样品的FPKM 值位于0 ～0.1 区间和0.3 ～3.57 区间（表7）。

2.3 差异表达基因分析

2.3.1 差异基因的筛选

本实验采用DESeq2[13]对在红壤、铅锌矿中生长的泡桐根部细胞转录组进行基因差异表达分析，将在不同样本间进行两两比较时表达量变化阈值为 padj ＜ 0.05 且 |log2FoldChange|＞ 1 的基因定义为差异表达基因（图3）。在红壤与铅锌矿之间，共筛选到4 297 个差异表达基因，其中4 143 个基因下调，154 个基因上调。

2.3.2 差异表达基因的GO 富集性分析

图2 KOC 注释分类Fig. 2 KOC annotation classification diagram

为寻找响应铅锌矿胁迫基因，本研究利用GOseq 比对红壤和铅锌矿的差异基因发现，在铅锌矿中，大部分差异基因呈下调趋势，极少一部分基因呈上调趋势（图4），说明泡桐根部组织主要通过降低一些基因的表达量来应对铅锌胁迫。结合（Binding）作为分子功能的主要类别，富集其中的差异表达基因高达1 655 条，主要包括蛋白质结合698 条，金属离子结合381 条，过渡金属离子结合250 条，锌离子结合181 条等。在生物过程类别中，细胞蛋白代谢（Cellular protein metabolic process）所含差异表达基因最多（482 条），另外，本研究发现细胞成分组织或生物发生（Cellular component organization or biogenesis）、细胞蛋白修饰（Cellular protein modification process）在生物过程类别中含有相对较多的应答胁迫差异基因。同时在细胞组分类别中，细胞器（Organelle）、细胞内细胞器（Intracellular organelle）、膜结合细胞器（Membrane-bounded organelle）所含差异基因数量位居前3 位。

表5 注释Unigene 数量最多的10 个代谢通路Table 5 Ten metabolic pathways annotated by Unigene

表6 Reads 与参考序列比对Table 6 Reads and reference sequence alignment

2.3.3 差异表达基因KEGG 富集性分析

为揭示泡桐根部受到铅锌矿胁迫时差异表达基因代谢途径和信号转导途径的富集性，本研究将差异表达基因序列在KEGG 数据库中进行富集性分析发现，所有差异表达基因共富集于108 个KEGG 代谢通路，其中富集差异表达基因数最多的7 个代谢通路依次为：核糖体（Ribosome）、内质网蛋白加工（Protein processing in endoplasmic reticulum）、剪接体（Spliceosome）、胞吞作用（Endocytosis）、RNA 转运（RNA transport）、吞噬代谢（Phagosome）、mRNA 监测（mRNA surveillance pathway）（表8）。这些通路涉及到泡桐根在铅锌胁迫下蛋白质加工、遗传信息加工、信号防御等方面，表明泡桐根通过调控生命活动等各方面基因来应答铅锌胁迫。

表7 FPKM 区间分布Table 7 FPKM interval distribution

图3 样品间基因差异表达分析火山Fig. 3 Volcanic maps of gene differential expression analysis between samples

核糖体在蛋白质的翻译过程中扮演着重要角色，代谢通路富集性分析发现此途径中共包含111个差异表达基因，在铅锌矿胁迫下，108 条基因表达下调，3 条基因表达上调，所有下调基因均与核糖体蛋白合成有关，涉及到26 个小亚基核糖体蛋白，42 个大亚基核糖体蛋白，且有实验表明，单个核糖体蛋白减少会导致生物发育迟缓，活性降低[14]。另外在Ni 胁迫下，纤维细胞中与核糖体蛋白相关的基因表达呈下降趋势[15]，说明植物可能在重金属胁迫下减少蛋白质合成，减缓细胞生长。

2.3.4 铅锌矿胁迫下的转录因子分析

在泡桐根转录组注释信息中，一共鉴定得到37 个家族的3 506 个的转录因子（表9），通过对这些转录因子进行深入分析，发现HSF，C2H2，MYB 家族是注释信息较多的种类，此外，还有一些诸如bZIP[16]、ARF、WRKY[17]等与环境相关的转录因子家族也得到鉴定。众多的转录因子共同调控泡桐根部功能基因，形成完整网络来应答铅锌矿胁迫。

2.3.5 差异表达基因qRT-PCR 验证

为了验证转录组数据的准确性，本研究随机选择6 个差异基因进行qRT-PCR 验证，其对应的引物如表10 所示。结果显示，荧光定量PCR 与RNA-seq 检测得到的基因相对表达量存在一定差异，但是两者反应的上下调趋势一致（图5），因此表明转录组测序得到的结果真实可靠。

3 结论与讨论

新一代测序技术已经应用于许多无参考基因的植物包括红花[18]、甘薯[19]等，目前，关于泡桐的遗传背景研究较少，本研究应用RNA-seq 高通量测序技术对铅锌矿生长区和红壤生长区的泡桐根进行转录组测序，共得到366 701 960 条reads，总碱基数为53.90 Gb，经过滤得到259 921 130 条clean reads，且Q20、Q30 均达90%以上，说明本次得到的测序数据质量较好，再经Trinity 拼接组装得到126 061 个基因序列。这些数据揭示了泡桐根部转录组遗传信息，为寻找响应铅锌矿胁迫基因提供基础。

将数字基因表达谱与高通量测序技术相结合是目前针对无参物种寻找差异基因最高效、最便捷的方法[20]，为寻找响应铅锌矿胁迫基因，本研究以红壤中泡桐根基因表达水平为参照来检测铅锌矿胁迫下泡桐根中差异表达基因表达情况，发现在两个样品中共有4 297 个差异表达基因，其中4 143 个基因下调，154 个基因上调。说明铅锌矿胁迫已经影响了基因的正常表达水平，并且主要以下调基因表达模式来应对重金属胁迫。通过GO富集分析，在铅锌胁迫下，显著富集下调基因主要涉及蛋白质结合、金属离子结合、过度金属离子结合、锌离子结合、细胞蛋白代谢、细胞蛋白修饰、细胞器等，其上调基因主要和结合有关。推测铅锌胁迫引起了泡桐根中与离子转运以及蛋白合成相关的基因发生变化，导致泡桐根部金属离子稳态以及蛋白合成发生变化。

图4 铅锌矿和红壤差异表达基因GO 分类结果Fig. 4 Result of GO function classification on differentially expressed genes in lead-zinc mine and red soil

表8 差异基因数最多的KEGG termTable 8 KEGG term with the largest number of differentially expressed genes

根系作为重金属深入植物的主要途径，其对重金属胁迫涉及到的代谢途径是多种多样的。有研究发现，当萝卜暴露于含铅的环境中，铅主要积累于根部，且可以通过细胞壁增厚来抵制铅进入原生质体[21]。另外，根细胞壁被认为是活跃的代谢位点，面对重金属胁迫时，可以产生丝裂原活化蛋白激酶MAPKs（Mitogen activated protein kinase）和钙结合蛋白CDPK（Calcium-binding related protein kinase）[22-23]。近年来，有研究发现，谷胱甘肽（GSH）可以通过抑制铅诱导产生的ROS，来保护植物[24]。槐叶苹中的铅积累可以增加谷胱甘肽数量，提高谷胱甘肽酶活性，且能提高根部SmGS 的基因表达量[25]。谷胱甘肽也可以对东南景天铅胁迫进行解毒[26]。金属转运蛋白可以将金属离子转运到细胞外或者螯合至液泡内，这些转运过程往往涉及到精确基因的表达调控。目前发现的这些基因主要来自不同家族转录因子，例如一些ABC 转运蛋白可以直接参与运输离子或其他溶质，甚至还参与调节通道。NARMP 可以对金属二价阳离子进行转运[27]。本研究中，我们发现了大量与重金属相关的转录因子如C2H2、MYB 等，这些转录因子在响应并抵挡重金属胁迫中起到重要作用。泡桐根部组织内与重金属胁迫的相关的代谢通路多种多样，为了确定哪些代谢通路更有研究价值，我们对这些差异基因进行KEGG 富集分析发现，这些代谢途径主要涉及到核糖体、内质网蛋白加工、剪接体、胞吞作用、RNA 转运、mRNA 监测、吞噬代谢等，可能这些代谢通路与重金属积累、解毒有关，是以后研究的重点。