郭 强，王英哲，王 昆，徐 博

（1. 吉林农业大学动物科学技术学院，吉林 长春 130118；2. 吉林省农业科学院，吉林 长春 130118；3. 黑龙江省二九一农场畜牧兽医科，黑龙江 集贤 155923）

金属污染物如铅(Pb)、镉(Cd)和汞(Hg)是非必需元素，可以通过人类活动释放到环境中，对动物和植物有毒害作用并且可以通过食物链威胁人类健康[1-2]。铅是一种剧毒的重金属，能够被植物吸收，抑制根和芽的生长，导致叶片萎黄，影响其他生理过程，如有丝分裂、蒸腾和DNA合成[3-4]等。因此，恢复土壤环境状态是当务之急，土壤清洗、土壤淋溶和植物修复等技术已被用于修复重金属污染的土壤，其中植物修复被认为是一种高效、经济、环保的清洁方法[5]。

了解基因的表达和调控在植物应对重金属毒性反应以及植物修复涉及的遗传和分子机制中显得尤为重要。植物重金属胁迫相关研究已经确定了一组基因家族和相关蛋白质参与金属运输，包括铁响应转运蛋白(IRT)、阳离子扩散促进剂(CDFs)、P型ATP酶、ATP结合盒(ABC)转运蛋白、天然抗性相关巨噬细胞蛋白(NRAMP)和Fbox家族[6]。参与重金属反应的基因在不同植物中的研究较广泛[7-10]；然而，基因组网络的基因表达对重金属胁迫的反应仍然未知。最近的研究使用新一代测序技术鉴定了许多microRNA，确认为重金属胁迫敏感的相关基因[11-13]。这些发现表明，植物中重金属诱导的反应是复杂的，并且涉及多种生理过程和代谢途径。

紫花苜蓿(Medicago sativa)是豆科植物，具有高度分枝的根系、高生物量，拥有对多种环境的广泛适应性。这些特点使紫花苜蓿成为一种理想的土壤改良植物。了解紫花苜蓿的分子基础对重金属的反应将有助于其在植物修复中的应用。

1 材料和方法

1.1 试验材料和处理条件

试验采取先发芽后移栽的方式进行苗期培养。选择20～30粒种子(陇东苜蓿，由吉林省农业科学院提供)放到铺有3层滤纸的培养皿中，适量浇水后放入HPG-400HX人工气候箱中培养，培养条件为(光照12 000 lx，14 h，温度25 ℃，湿度75%；无光照，10 h，温度20 ℃，湿度75%)。待种子发芽高度为1 cm时进行移栽，移栽到水培箱中(水培箱由箱体加浮板组成，浮板规格为28 cm × 18 cm，箱体规格为32 cm × 22 cm × 17 cm。营养液为Hoagland营养液[14]，每箱用量为3 L稀释液，每两箱栽种一份材料，每份材料栽种12棵苗，即每箱为6棵。将栽种好的材料分为两组，一组为对照组，一组为处理组，培养30 d。30 d后进行胁迫处理，处理组以200 mg·L-1的硝酸铅进行胁迫，用量为每3 L Hoagland营养稀释液浇灌200 mL硝酸铅溶液，处理96 h (记为Pb96)。对照组用蒸馏水处理，用量为每3 L Hoagland营养稀释液浇灌200 mL蒸馏水，处理96 h (记为C96)。处理后取植株根部(预试验得知根部较地上部效果明显) 0.15 g进行试验指标测定。取下的根经过液氮处理后保存于-80 ℃备用，各处理进行3次重复。本研究选取水培方法参考丁振中等[15]、程贝等[16]的研究确定。

1.2 RNA提取以及转录组测序

取Pb处理和对照植物根部的总RNA，使用RNAprep Pure Plant(Bio-Teke，China)试剂盒提取Pb96和C96。使用Qubit 3.0 (Thermo Scientific，USA)和Agilent 2100生物分析仪(Agilent Technologies，USA)评估RNA的浓度和质量。

使用NEBNext1 Ultra TM Directional RNA Library Prep Kit (NEB，UK)试剂盒根据说明制备cDNA文库。分别使用Qubit2.0和Agilent 2100对文库进行检测，使用q-PCR方法对文库的有效浓度进行准确定量，以保证文库质量，准备好的cDNA文库在Hiseq 4000(Illumina)上进行测序。

1.3 测序数据的过滤、组装

在Q20条件下使用FastQC工具(http://www.bioinformatics.bbsrc.ac.uk/projects/fastqc/)去除其中的接头序列及低质量序列，获得高质量的去冗余序列 (Clean Reads)。采用 Trinity(版本2.2.0，默认参数)软件进行序列组装，获得紫花苜蓿的功能基因(Unigenes)。

1.4 差异表达基因(differentially expressed genes，DEGs)的功能注释与富集分析

使用FPKM方法计算基因表达水平并标准化。使用Student's t-test检验不同文库之间的差异，P值根据Benjamini和Yekutieli[17]的方法进行调整。最后，以 FDR ＜ 0.005 且|log2(fold-change)| ≥ 1 为条件，筛选出Pb96和C96文库中的DEGs。

利用BLAST程序(E-value = 10-5)将DEGs序列分别与COG、GO、KEGG(数据库进行比对，获得DEGs的注释信息。

1.5 q-PCR验证

PCR反应体系包括：SYBR premix Ex TaqTM(2X)10 μL，上游引物和下游引物各1 μL，cDNA 模板2 μL，RNase Free dH2O 6 μL。PCR 程序为：95 ℃ 预变性10 min；95 ℃ 变性 15 s，60 ℃ 退火 20 s，72 ℃ 延伸40 s；40个循环；80 ℃ 1 s收集荧光，读板；之后72 ℃延伸10 min，最后从65 ℃开始，以每步0.5 ℃的速度升高至95 ℃，每个温度保持1 s，绘制溶解曲线。每个样品3个生物学重复。选用和ac-tin(CpActin)作为内参基因。

2 结果与分析

2.1 转录组测序与组装

提取的RNA总浓度在324.8 ng·μL-1以上，表明RNA质量合格。Pb胁迫处理96 h与对照组的根部进行转录组测序后共获得100 007条Unigene序列，所有序列N50值为1 097，在所有的Unigene中，长度在200～300 nt的最多，占Unigene序列总数的37.21%，长度超过2 000 nt的序列达到6 297条，占Unigene序列总数的6.30%(表1)。研究结果表明随着Unigene序列长度的增加，Unigene序列呈现减少的趋势，表明测序结果良好。

2.2 差异表达基因的鉴定

以FDR ≤ 0.01和|log2Ratio| ≥ 1作为DEG筛选条件，Pb96处理组与C96对照组相比，96 h Pb处理后共得到2 242个DEGs，占总Unigene序列的2.24%。其中，1 321个DEGs上调表达，921个DEGs下调表达。

表1 功能基因长度以及数量Table 1 Functional gene length and quantity

2.3 差异表达基因的功能注释

通过BLAST(Evalue = 10-5)比对分析，2 242个DEGs序列GO数据库中注释到的DEGs数目最多，达到2 162个，占DEGs总数的96.4%；最少的为KEGG数据库，注释到的DEGs数目为448个，占DEGs总数的19.9%(表2)。

表2 差异表达基因的功能注释结果Table 2 Functional annotation results of differentially expressed genes

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

COG数据库的构建源自于细菌、藻类和真核生物的系统进化关系，利用该数据库可以对基因产物进行直系同源分类。被注释到COG数据库的696条DEGs覆盖所有 21个类别，其中“碳水化合物的运输和新陈代谢”(Carbohydrate transport and metabolism)数量最多，为70个(10.06%)；具有“翻译、核糖体结构和生物发生”(Translation，ribosomal structure，and biogenesis)、“次生代谢产物的生物合成、转运和分解代谢”(Secondary metabolites biosynthesis，transport，and catabolism)和“仅通用功能预测”(General function prediction only)功能的DEGs数量较高，分别为 63 个(9.05%)、64个(9.20%)和68个(9.77%)(图1)。可以看出，测序获得的转录本参与了多数代谢过程，这些测序数据将是研究紫花苜蓿功能基因的可利用资源。

KEGG数据库可以用来预测DEGs在细胞中参与的代谢途径和功能，了解基因参与的生物学过程。将DEGs进行通路分析，发现共有448个DEGs参与到97个通路当中(部分展示，图2)。其中DEGs参与最多的是“核糖体”(Ribosome)通路ko03010，共有48个DEGs参与；其次是“光合作用”(Photosynthesis)通路ko00195、“碳代谢”(Carbon metabolism)通路ko01200、“淀粉和蔗糖代谢”(Starch and sucrose metabolism)通路ko00500等。由此可知，植物在铅胁迫状态下，光合作用发生改变，同时体内的淀粉蔗糖代谢以及碳代谢均有明显变化。

GO数据库是一个结构化的标准生物学注释系统，其信息适用于各物种。GO分类体系分为3 个大类，包括生物学过程(Biological process)、细胞元件(Cellular component)和分子功能(Molecular function)，以及57 个亚类。本研究中共有 2 242个DEGs被注释到GO数据库，分布于3个大类的24个亚类中，其中生物学过程中有9个亚类，细胞元件中有6个亚类，分子功能中有9个亚类(图3)。被富集到物质代谢 (Metabolic process)245个 (GO:0008152)、蛋白结合(binding)214个(GO:0005488)、催化活性(Catalytic activity)202个(GO:0003824)、细胞过程 (Cellular process)162个(GO:0009987)、细胞(Cell)122个(GO:0005623)5个亚类的DEGs数量明显较其他亚类多。说明植物在受到铅胁迫时，植物体内蛋白含量和酶活性有明显的变化。

通过注释得到主要调节耐铅机制的差异表达基因如表3所列，其中下调表达有12个，上调表达有10个。选择其中差异最为显著的c53338.graph_c0、c40437.graph_c1、c60665.graph_c0、c44469.graph_c0、c56258.graph_c0、c54935.graph_c1和c46893.graph_c0为候选基因，有待进一步分析。

2.5 qRT-PCR验证结果

选取参与调节耐铅性的10个差异表达基因，进行qRT-PCR验证。结果表明，C96表达量均小于Pb96，这与高通量测序结果表达趋势一致(图4)。

图1 注释到COG数据库的基因数目Figure 1 Notes the number of genes in the COG database

图2 注释到KEGG数据库的基因数目Figure 2 Notes the number of genes in the KEGG database

3 讨论与结论

铅是世界范围内广泛存在的土壤污染物，对食物安全与人类健康构成巨大威胁[18-19]。不同的物理、化学和生物方法都有被用来解决环境中的重金属污染问题，其中，植物修复被认为是一种高效、经济、环保的方法[20-21]。但是，这些基因参与了重金属毒性和分子的保护，这种保护的潜在机制在植物中仍未定义。

本研究通过BLAST软件与GO、COG、KEGG数据库对比最终获得100 007条Unigene序列，鉴定出差异表达基因2 242个。在GO数据库中注释到2 162条，在COG数据库中注释到696条，在KEGG数据库中注释到448条。其余未被注释的序列可能是紫花苜蓿所特有的基因序列，仍有待于进一步的分析。在GO注释中的物质代谢、催化活性2个通路在黄海燕[22]的研究中同样被注释到，说明紫花苜蓿在铅胁迫条件下，酶活性起到相对较高的作用。KEGG注释中的“碳代谢”、“淀粉和蔗糖代谢”在宋丽莉[23]的紫花苜蓿抗铅研究中也同样被注释到，说明淀粉蔗糖的代谢与合成能够调节紫花苜蓿的耐铅性。而在COG数据库中注释到的“翻译、核糖体结构和生物发生”、“次生代谢产物的生物合成、转运和分解代谢”功能的DEGs数量较高，这与齐晓[24]的研究结果相似，表明紫花苜蓿在抗铅反应中主要调节作用在于物质的合成与分解。

图3 注释到GO数据库的基因数目Figure 3 Notes the number of genes in the GO database

图4 实时荧光定量PCRFigure 4 Real-time PCR

王少峡等[25]研究发现，DREB转录因子由逆环境胁迫诱导产生后，可激活其他多达 12个依赖DRE顺式作用元件的抗逆功能基因，引起脯氨酸及蔗糖含量提高，从而增强植株对多种逆境的抵抗性。而在本研究中同样注释到“淀粉和蔗糖代谢”通路，表明差异表达基因促进了淀粉蔗糖代谢过程，也充分证明了淀粉蔗糖代谢对紫花苜蓿的抗铅性至关重要。

次生代谢产物是植物长期进化中适应环境的产物，其公认的生态学功能主要是抗病、抗虫、抗环境胁迫[26]等。而在COG数据库中同样注释到“次生代谢产物的生物合成、转运和分解代谢”，且注释到的差异基因数量为64个，位居第二，表明次生代谢产物对于紫花苜蓿适应逆境环境具有较为重要的作用。

n e n t n i o p o a l r t e r，a t e r c h a n n e l a c t i v i n t e g r a l c o m 性活i t y，酶接合泡液连结o A合糖央-C结多a l s u b u n i t s中酸丁离子，分运运，转o t p e p t i d e)属成i c l a r g e r i b o s o m程e t a b o l i s m，金转体膜s t a s i s，w分过跨、整谢膜p h o s p h o r y l a t i o n，成r t e r a c t i v i t y，的跨)键i c l a r g e r i b o s o m子e o 体膜代g l y o x y l a t e c y c l e，盐，整肽，e s，c y t o p l a s m，的G e n e f u n c t i o n 酸活白s t r a t e-s p e c i f i c t r a n s m e m b r a n e t r a n s p o d i n g 膜离分r a n e t r a n s p o化是r a n e t r a n s p o r t，a m i n o a c i d t r a n s m e m b 亚动a t e r h o m不中a c t s o n c a r b o n-n i t r o g e n (b u t n亚e s，c y t o p l a s m 丁活，基p o l y s a c c h a r i d e b i n (但蛋，酸体磷成体道e t a b o l i s m，程环-氮质运整通过性糖基核循转白谢的膜n e n t，s u b碳能膜酸代p o u n d s，c a r b o h y d r a t e m d i n g水体i t y的蛋功大f r i b o s o m ，质跨分c t u r a l c o m 醛p o大p o n e n t s o f r i b o s o m，态物因性i t y，i n乙胞结A T 糖性稳性合l c o m合脒，P b i n d i n g，p r o t e i n活核基，异水，a c e t i c a c i d m化g e n e s 性T P 特胞程水成线活，y d r o l a s e a c t i v r a l c o m b r a n e，l y c o s y物构s t r u 蛋质细过l i g a s e a c t i v碳于，分p o n e n t s o谢e t a l i o n b，A白物f O 酶表r e s s e d ，胞o s p h a t e a l d o l a s e a c t i v i t y N u t r i e n t s t o r a g e a c t i v i t i e s，运b l y 缩成b r a n e i n t e g用的b l y，o A 底a c t i v i t y，t r a n s m e m b r a n e t r a n s p o r t 结膜p l a n t-t y p e t o n o p l a s t，c e l l w，r t，c e n t r a l v a c u o l e性代醋l i g a s e a c t i v i t y，作-g e m 活醛分泡膜酸合体e m酸液构化性磷解r o c e s s，h a n d m糖型p o n e n t s o f t h e m n d s i n l i n e a r r u t h e n i u m成跨列酶酸，基激基因- 二结l l y e x p活r a l c o m核体e m整p o n e n t s，物性e A b u t y r a t e-C d r o l y s i s o b r a n e基异i f f e r e n t i a 酸转酶氨的，u n i t a s s e m的植活O-糖h y o f t h e m达糖氨，体装膜，酶解表果水水H y d r o l a s e a c t i v i t y, f r u c t o s e-d i p h/苏运糖动组b r a n e ，分接o e n z y m水，，核r e o n i n e k i n a s e a c t i v i t y，酸N u c l e o s o m e a s s e m活基膜成b r a n e c o m连o l i s m，，p o s i t i o n性3 差程氨表L i s t o f d 膜，c o l y s i s p存装亚泡A e t a b性跨活过e m 丝合储成e m组大液p e t o n o p l a s t，m组b r a n e t r a n s p o酶活酸酶解结酵e g l y 质物体体型t-t y 的-辅酶r a t e m解糖T h转b o 白基蛋P r o t e i n s e r i n e/t h i n t e g r a l c o m氨A m i n o a c i d t r a n s m e m b A A b i n d i n g，s t r u c t u养小糖物膜酸t y 解e r a l l c o m b l e 3 o f t h e m R N R N 营核核R i b o s o m e l a r g e s u b s u b u n i t s植P l a n质P l a s m a m e m t r a n s m乙A c e t i c a c i d - C b u 水H y d r o l a s e a c t i v i t y，o v T a比g 2 r a t i o 1 8值4 4 2 0 7 4 6 0 1 4 5 3 4 6 9 4 1 8 6 5 8 8 4 5 6 5 8 8 4 5 3 1 8 0 2 8 3 1 8 0 2 8 0 9 5 6 2 4 9 9 2 4 5 9 6 2 6 4 6 5 6 9 5 6 6 1 l o-1.0 6 3-1.6 3 3-3.6 3 1-3.6 2 5-3.6 2 5-3.5 5 8 5 8 2 8 3-3.5-3.5 5 8-3.4 5 1-2.9 5 2-2.6 9 3-1.9 6 1 M P K P b 9 6_F 6.1 5 1 0 6.6 5 0.1 3 0.2 5 0.0 8 0.1 4 0.0 4 0.0 5 0.2 3 0.0 7 0.7 1 0.7 9 K M F P C 9 6_1 9 8.6 9 1 8.9 7 2.2 2 6 1.9 1 2.9 9 0.8 7 1.2 6 3.3 2 0.7 4 4.5 3.0 7 e n e名U n i g i g e n e I D r a p h_c 0.g r a p h_c 1.g r a p h_c 0.g r a p h_c 0.g r a p h_c 0.g r a p h_c 0.g r a p h_c 0.g r a p h_c 0.g r a p h_c 0.g r a p h_c 0.g r a p h_c 0.g r a p h_c 0.g U n 4 7 8 2 c 5 0 4 3 7 c 4 3 9 9 2 c 4 0 6 6 5 c 6 0 9 0 7 c 4 3 8 1 1 c 3 1 6 9 4 c 3 0 9 8 4 c 3 5 5 0 0 c 3 4 6 3 8 c 4 4 4 6 9 c 4 0 2 9 5 c 4

?

在上下调基因中(表3)基因的主要功能是调节质膜、酶活以及跨膜转运和离子结合，其中在酶活调节过程中，主要调节辅酶、水解酶以及氧化还原酶。在调节跨膜转运过程中主要调节水的转运和离子转运。同时在基因调解过程中还有其他部分细胞器的参与。综上所述，通过数据库的比对发现，物质的合成与分解、淀粉蔗糖代谢以及次生代谢产物3大过程均对于植物耐铅性起到调节作用。

本研究中将差异表达基因与3大数据库对比发现许多与耐铅相关的通路，在紫花苜蓿中获得了大量响应铅胁迫的DEGs，并在转录组水平鉴定了参与调控紫花苜蓿铅胁迫反应的基因，从而有利于从整体水平加深了解紫花苜蓿耐铅反应表达特性，为进一步研究这些响应铅胁迫的基因的功能提供参考，也可为紫花苜蓿耐铅功能研究提供新的基因资源。