水稻芽期响应镉胁迫的转录组分析

王燕宁，黄 涛，吴光亮，黄诗颖，钟 奇，王 鹏，杨朦朦，李才敬，程 琴，贺浩华，2，金 伟，郭 凌，边建民，2

(1.作物生理生态与遗传育种教育部重点实验室，江西 南昌 330045；2.江西省水稻高水平工程研究中心，江西 南昌 330045；3.南昌市农业技术推广中心，江西 南昌 330200)

水稻是我国的主要粮食作物之一，具有富集镉的特性[1-2]。水稻根、茎、叶、籽粒容易富集镉，严重影响水稻的光合作用、蒸腾作用和其他生理过程，而且会产生过量的氧自由基影响水稻的正常生长发育[3]，镉还会通过稻米进入人体，危害人类健康[4-5]。近年来，水稻镉危害引起人们广泛关注，因此，了解水稻镉胁迫的遗传机制，培育耐镉水稻品种对于降低水稻镉危害、保障稻米安全具有重要作用。

水稻镉胁迫相关性状属于复杂性状，受到多个遗传位点共同调控，近年来研究人员通过外源镉处理的方法，定位了多个与根长、株高、干质量和镉含量相关的QTLs，大多数效应较小。Yan等[6]通过一个重组自交系群体(RILs)定位到与镉含量相关的5个主效QTLs(scc10、gcc3、gcc9、gcc11和sgr5)，其中sgr5是一个新的潜在QTL，对镉从地上部到籽粒的分布影响较大。Ishikawa等[7]利用Sasanishiki和Habataki构建的回交自交系(BILs)定位群体在7号染色体上鉴定到一个新型主效QTLqGCd7，qGCd7是提高水稻籽粒镉含量的特异性QTL。骆旭添[8]利用Lemont(美国)和Dular(印度)杂交构建的RILs，以叶绿素含量、根长、株高、叶长为指标共检测到16个加性QTL，涉及水稻第1，2，3，4，6，7，11号等7条染色体。Pan等[9]对生长在不同程度镉污染土壤上的338份不同稻种资源的镉积累能力进行了评估，通过全基因组关联分析(GWAS)研究发现了35个与镉积累显著相关的QTL。Li等[10]利用CH891和02428构建的BILs定位了镉胁迫下种子发芽率相关的6个加性QTL和16对上位性QTL。

镉相关基因的克隆也取得了一定的进展，在水稻12条染色体上均有分布。Ueno等[11]利用镉高积累品种Anjana Dhan和日本晴构建的定位群体，克隆了水稻低镉积累基因OsHMA3，该基因主要影响镉从地下向地上部的转运效率以及地上部的镉积累，此后发现该基因在品种间存在多个自然变异类型。Luo等[12]利用镉积累差异显著的2个品种台中1号和春江06，构建杂交群体克隆到了一个特异性调控水稻叶片镉积累的主效QTL—CAL1基因。该基因编码植物防御素类似蛋白，定向调控镉在叶片等营养器官的积累。锌铁调控转运蛋白基因OsZIP3能降低水稻根中的镉含量[13]。Chen等[14]分离出一个水稻耐镉突变体，并认为OsCADT1是造成突变体耐镉的基因。LCD是一个水稻低镉基因，定位于细胞质和细胞核，被认为是低镉水稻育种中的潜在可利用基因。镉胁迫下，其功能缺失型突变体植株的籽粒镉含量比野生型植株减少1/2[15]。OsCd1是Yan等[16]利用127份自然群体材料通过GWAS定位到的QTL中分离到的水稻籽粒镉积累基因，OsCd1属于MFS超家族转运蛋白，能转运镉离子，参与根系镉吸收。

多年来，国内外学者已经对水稻镉胁迫的遗传机理进行了许多研究，但主要集中在成熟期镉的吸收和转运，而水稻幼芽期对镉胁迫较为敏感[17-18]，因此，如何降低镉胁迫对水稻幼芽的危害也是目前抗逆育种急需解决的问题。研究还表明，不同水稻品种对镉的吸收和积累能力存在明显差异[19]，籼稻和粳稻在镉胁迫下的抗逆性也有显著差异，但是造成这种差异的机制还不清楚[20]。

本研究以籼稻品种CH891和粳稻品种02428为供试材料，利用RNA-Seq技术对2个品种在正常条件和镉处理3 d后幼苗的转录本进行测序，并对转录组测序结果进行GO富集分析和KEGG通路富集分析，旨在为阐明籼、粳稻响应镉胁迫的遗传差异提供理论依据。

1 材料和方法

1.1 试验材料

供试材料为镉敏感的籼稻恢复系CH891与耐镉的粳稻品种02428。

1.2 水稻种子镉处理与样品采集

挑选健康饱满的谷粒用1%的次氯酸钠溶液浸泡消毒15 min，用蒸馏水冲洗3～4次。将消毒过的谷粒浸种24 h，25 ℃催芽1～2 d。分别选取CH891和02428大小、露白基本一致的30粒种子放置于铺有单层滤纸的培养皿中，处理组用浓度为200 mg/L的CdCl2溶液处理，对照组用相同体积的蒸馏水处理，处理组和对照组均设置3次生物学重复。处理组和对照组放入同一个光照培养箱内进行培养。白天温度设置为(27.0±0.5)℃，夜间温度设置为(25 .0±0.5)℃，光期时长设为14 h，暗期时长设为10 h。第3天时将2个品种水稻镉处理和正常条件下的幼苗分别取下，用锡箔纸包裹置于液氮中速冻，放置在-80 ℃冰箱保存备用。

1.3 总RNA提取

利用TRIzol试剂提取第3天镉处理和对照组CH891和02428的RNA，利用1%的琼脂糖凝胶电泳检测RNA质量，利用2100 Bioanalyzer检测RNA纯度和浓度，RNA质量合格后送去联川生物公司。

1.4 cDNA文库构建和转录组测序

使用mRNA片段作为模板构建cDNA文库，cDNA文库构建由联川生物公司完成，最后在Illumina HiSeq 4000平台上进行双末端测序[21]。

1.5 生物信息学分析

将原始读长数据(Raw read)进行过滤，将过滤掉低质量reads的高质量序列(Clean read)比对到参考基因组，对所有比对上的基因进行表达水平分析，并用FPKM值表示各基因表达水平，再进一步进行基因表达水平差异研究和差异基因功能分析[22-23]。

1.6 差异基因表达分析

利用edgeR对StringTie组装和定量完毕的基因进行差异分析(显著差异的阈值为|Log2Fold Change|≥1，P<0.05)。Fold Change是基因的变化倍数值，一般认为数值在2以上的表示该基因的表达量在2个比较组之间存在显著差异。P-value是各个基因差异表达概率，其数值一般控制在0.05水平内，即<0.05表示在95%的置信区间内该基因是差异表达基因(可过滤假阳性差异表达基因)。此外，通过绘制火山图可以了解差异表达基因的整体分布情况，根据样品基因表达谱的相近程度，将差异基因进行聚类分析，直观地展示基因在不同样品(或是不同处理)中的表达情况，由此获取生物学相关信息[24]。

1.7 差异表达基因功能注释和分类富集分析

所有差异表达基因均在Gene Ontology(GO)和KEGG(Kyoto Encyclopedia of Genes and Genomes)中分别进行GO和KEGG注释，网址分别为http：//www.geneontology.org/page/download-go-annotatios和www.genome.jp/kegg。利用WEGO和OMICSHARE进一步对注释得到的差异表达基因进行GO和KEGG分类富集分析，网址分别为http：//wego.genomics.org.cn和http：//www.omicshare.com/tools/index.php/。

1.8 qRT-PCR验证

为了验证RNA-Seq结果，通过定量实时RT-PCR(qRT-PCR)证实了几个基因的不同表达模式。按照TaKaRa反转录试剂盒说明书对本试验样品RNA进行反转录。按照cDNA 1 μL，Premix Ex TaqTMⅡ 10 μL，ddH2O 7 μL，Rox Reference Dye 0.4 μL，上下引物10 μmol/L各0.8 μL的比例配置qRT-PCR反应体系(表1)。反应程序设定为：95 ℃ 1 min；95 ℃ 15 s，55 ℃ 30 s，72 ℃ 60 s，65 ℃ 5 s，95 ℃ 5 s，40个循环，每个基因设置3次生物学重复和3次技术性重复。利用qRT-PCR所得Ct值计算样品中基因的相对表达量(2-ΔΔCt)[25-26]。

表1 qRT-PCR引物序列Tab.1 The sequences for qRT-PCR primers

2 结果与分析

2.1 镉胁迫下CH891和02428的表型差异

正常情况下籼稻CH891幼芽显著长于粳稻02428，镉处理后，CH891和02428的生长均受到明显抑制，其中CH891受到抑制的程度大于02428，芽长显著变短，但镉处理3 d后CH891和02428的芽长没有显著差异，表明CH891幼芽生长对镉的敏感性高于02428(图1)。

2.2 样品RNA质量检测与测序序列统计和质量评估

使用高通量 Illumina HiSeq 4000(CA，美国)对CH891和02428在镉处理3 d和正常生长的幼芽进行RNA-Seq，总共产生545 010 200个原始测序读长数(Raw read)，平均每个样本产生45 417 517个原始测序读长数。过滤掉不合格序列后共得到539 524 490个有效读长(Valid read)，占原始读长的99.01%以上。所有样品测序读长与参考基因组比对率为94.81%～96.82%，其中唯一比对率为62.11%～74.21%，所有样品的GC含量均在49%以上，测序质量较好，可进行下一步试验分析(表2)。

表2 测序数据统计和参考基因组上的比对率Tab.2 Sequencing data statistics and mapped ratios to the reference genome

2.3 差异表达基因数量分析

各比较组中差异表达基因数目如图2所示。对照条件下，CH891与02428(CK3_CH891 vs.CK3_02428)之间存在4 911个差异表达基因，其中上调表达基因1 647个，下调表达基因3 264个；镉处理条件下，CH891与02428(Cd3_ CH891 vs.Cd3_02428)之间存在2 542个差异表达基因，其中有1 499个上调表达基因，1 043个下调表达基因；CH891镉处理后与对照相比(Cd3_CH891 vs.CK3_CH891)共存在1 952个差异表达基因，其中上调表达基因1 203个，下调表达749个基因；02428镉处理后与对照相比(Cd3_02428 vs.CK3_02428)共有1 792个差异表达基因，其中有516个上调表达基因，1 276个下调表达基因。

A、B、C和D.火山分布图；E.上下调差异表达基因个数。A，B，C and D.Volcano plot；E.The number of up and down-regulated DEGs.

2.4 qRT-PCR验证

为验证RNA-Seq数据的准确性，以OsActin作为内参基因，随机挑选了12个差异表达基因(LOC_Os03g08470、LOC_Os04g39880、LOC_Os09g36700、LOC_Os10g40720、LOC_Os02g02400、LOC_Os04g3830、LOC_Os06g21590、LOC_Os09g31430、LOC_Os01g18744、LOC_Os06g13280、LOC_Os03g03600、LOC_Os08g06110)，使用特异性引物进行qRT-PCR分析，将qRT-PCR结果与差异表达基因的Log2(Fold Change)值进行比较，结果表明，qRT-PCR结果与RNA-Seq趋势一致(图3)，表明转录组测序结果准确可靠。

图3 差异表达基因转录组与qRT-PCR比较Fig.3 Comparison of differentially expressed gene transcriptome and qRT-PCR results

2.5 差异表达基因分类

通过对同一水稻品种在不同环境和不同水稻品种在相同环境的差异表达基因进行比较，共获得4个比较组，分别为Cd3_CH891 vs.Cd3_02428、CK3_CH891 vs.CK3_02428、Cd3_CH891 vs.CK3_CH891和Cd3_02428 vs.CK3_02428。在4个比较组中共存在7 204个差异表达基因，这些差异表达基因可以分为15个亚组(图4)。去掉CK3_CH891 vs.CK3_02428与镉无关的基因，剩余的7个亚组可以分为3个类别：CH891特异表达的基因(SCR3)、02428特异表达的基因(RCR3)、2个品种均有表达但表达量存在差异的基因(CCR3)。其中，SCR3有849个，RCR3有676个，CCR3有770个(表3)。

图4 4个比较组差异表达基因的韦恩图Fig.4 Venn diagrams for DEGs in the four comparison groups

以|log2(Fold change)| ≥2并且至少在一个比较组的一组数据中FPKM≥2为筛选条件，分别对3个类别的基因进行筛选，共获得554个主效的镉相关基因，其中194个来自SCR3，154个来自RCR3，206个来自CCR3。

表3 差异表达基因的分类Tab.3 Classification of three categories of DEGs

2.6 基因功能注释

2.6.1 SCR3基因功能注释 GO富集分析显示，转录调节和DNA模板、细胞核、膜的整体结构、质膜相关的基因、RNA结合、转录活性和序列特异性DNA结合、金属离子结合相关的基因在SCR3中均有富集(图5)。KEGG富集分析中，异黄酮生物合成、磷酸肌醇代谢、黄酮类生物合成和花青素生物合成等途径在SCR3中富集，其中异黄酮生物合成和花青素生物合成途径显著富集(图6)。

1.生物进程；2.转录调节、DNA模板；3.转录、DNA模板；4.蛋白磷酸化；5.代谢过程；6.伤害反应；7.氧化还原反应；8.脱落酸反应；9.油菜素类固醇反应；10.脂质转运；11.缺氧反应；12.蛋白质泛素化；13.蛋白酶体介导的泛素化依赖性蛋白质分解代谢过程；14.冷反应；15.多细胞生物发育；16.种子休眠中组织胚胎发育；17.叶发育；18.光刺激反应；19.蛋白质折叠；20.类黄酮葡萄糖醛酸化；21.胞外分泌；22.次级代谢生物合成；23.防御反应；24.赤霉素反应；25.甲基化作用；26.核；27.膜整体结构；28.质膜；29.细胞质；30.细胞外区域；31.胞液；32.叶绿体；33.膜；34.高尔基体；35.线粒体；36.细胞壁；37.植物细胞壁；38.内质网；39.质膜固定组成成分；40.胞间连丝；41.分子功能；42.RNA结合；43.转录活性、序列特异性DNA结合；44.金属离子结合；45.蛋白结合；46.ATP结合；47.羧基酯水解酶；48.蛋白丝氨酸/苏氨酸激酶活性；49.锌离子结合；50.水解酶活性。

1.α亚麻酸代谢；2.萜类骨架生物合成；3.类固醇生物合成；4.淀粉和蔗糖代谢；5.倍半萜和三萜类生物合成；6.自噬调节；7.RNA聚合酶；8.RNA降解；9.嘧啶代谢；10.卟啉和叶绿素代谢；11.植物激素信号转导；12.光合作用；13.磷脂酰肌醇信号系统；14.类苯丙酸生物合成；15.苯丙氨酸、酪氨酸和色氨酸生物合成；16.磷酸戊糖途径；17.N-聚糖生物合成；18.异黄酮生物合成；19.磷酸肌醇代谢；20.糖基磷脂酰肌醇生物合成；21.硫代葡萄糖苷生物合成；22.类黄酮生物合成；23.醚类脂代谢；24.二帖生物合成；25.角质、木栓质和蜡质的生物合成；26.类胡萝卜素；27.光合生物的碳固定；28.生物合成；29.碱基切除；30.花青素生物合成。

2.6.2 RCR3基因功能注释 GO富集分析显示，转录调控和DNA模板、细胞核、质膜、膜整体结构、和胞外区域相关的基因、转录因子活性、序列特异性DNA结合、蛋白质结合相关的基因在RCR3都有富集(图7)。KEGG富集分析显示，核糖体、自噬调节和胰岛素抗性等途径在RCR3中富集，其中核糖体和自噬调节富集程度最为显著(图8)。

2.6.3 CCR3基因功能注释 GO富集分析显示，转录调节和DNA模板、细胞核、细胞质、膜整体结构、质膜相关的基因、转录因子活性、序列特异性的DNA结合、蛋白质结合、DNA结合、锌离子结合相关的基因在CCR3中富集(图9)。KEGG富集分析显示，α-亚麻酸代谢、吞噬体、过氧化物酶体、柠檬烯和蒎烯的降解以及脂肪酸降解途径在CCR3中富集，其中α-亚麻酸代谢和脂肪酸降解途径极显著富集(图10)。

1.生物过程；2.转录调控、DNA模板；3.转录、DNA模板；4.氧化还原反应；5.防御反应；6.蛋白磷酸化；7.细胞壁组织；8.脱落酸反应；9.伤害反应；10.盐胁迫反应；11.蛋白质水解；12.多细胞生物发育；13.一维细胞生长；14.氧化应激反应；15.冷反应；16.独脚金内酯生物合成；17.茉莉酸反应；18.蛋白酶体介导的泛素依赖性蛋白质分解代谢过程；19.信号转导；20.脂质转运；21.类黄酮生物合成过程；22.卡里金反应；23.涉及多维细胞生长的植物型细胞壁修饰；24.蛋白转运体；25.胞吐作用中的囊泡对接；26.核；27.质膜；28.膜整体结构；29.胞外区域；30.细胞质；31.线粒体；32.膜；33.叶绿体；34.细胞溶质；35.胞间连丝；36.细胞壁；37.细胞组分；38.液泡；39.植物细胞壁；40.细胞内有膜的细胞器；41.分子功能；42.转录活性、序列特异性DNA结合；43.蛋白结合；44.DNA结合；45.蛋白丝氨酸/苏氨酸激酶活性；46.ATP结合；47.金属离子结合；48.激酶活性；49.核糖体的结构成分；50.锌离子结合。

1.mRNA监控途径；2.β-丙氨酸代谢；3.玉米素生物合成；4.缬氨酸、亮氨酸和异亮氨酸降解；5.囊泡运输中的SNARE相互作用；6.核糖体；7.自噬调节；8.丙酸代谢；9.植物-病原体相互作用；10.植物激素信号转导；11.苯丙酸类生物合成；12.苯丙氨酸、酪氨酸和色氨酸的生物合成；13.磷酸戊糖；14.氧化磷酸化；15.其他类型O-聚糖生物合成；16.其他聚糖降解；17.N聚糖生物合成；18.胰岛素抗性；19.糖苷磷酸肌醇生物合成；20.醣酵解/糖异生；21.果糖和甘露糖；22.二帖生物合成；23.角质、木栓质和蜡质的生物合成；24.TCA循环；25.类胡萝卜素生物合成；26.碳代谢；27.碱基切除修复；28.抗坏血酸盐和醛酸盐；29.精氨酸和脯氨酸代谢；30.花青素生物合成。

1.转录调节、DNA模板；2.生物进程；3.转录、DNA模板；4.伤害反应；5.蛋白质泛素化；6.氧化还原反应；7.冷反应；8.脱落酸反应；9.水杨酸反应；10.寡肽转运；11.防御反应；12.氧化应激反应；13.脂肪酸生物合成；14.脱水反应；15.参与细胞内蛋白质代谢过程的蛋白质水解；16.光刺激反应；17.一维细胞生长；18.乙炔激活信号转导；19.赤霉酸介导的信号转导；20.次级代谢生物合成；21.盐胁迫反应；22.几丁质反应；23.调控防御反应；24.其他有机物反应；25.木质素生物合成；26.核；27.细胞质；28.膜整体结构；29.质膜；30.膜；31.叶绿体；32.细胞内；33.胞外区域；34.线粒体；35.质膜固定成分；36.细胞溶质；37.高尔基体；38.质体；39.叶绿体基质；40.细胞壁；41.转录因子活性、序列特异性DNA结合；42.分子功能；43.蛋白质结合；44.DNA结合；45.锌离子结合；46.序列特异性DNA结合；47.转运蛋白活性；48.转移酶活性，转移糖酰基；49.蛋白泛素化转移酶活性；50.转录调控区DNA结合。

1.α亚麻酸代谢；2.缬氨酸、亮氨酸和异亮氨酸降解；3.锑烯类、二芳基庚烷类和姜酚生物合成；4.植物病原体互作；5.吞噬体；6.过氧化物酶体；7.烟酸盐和烟酰胺代谢；8.N聚糖生物合成；9.赖氨酸降解；10.亚油酸代谢；11.柠檬烯和蒎烯降解；12.异黄酮类生物合成；13.糖磷脂生物合成；14.糖磷脂生物合成-神经节系列；15.糖酵解/糖异生；16.甘氨酸、丝氨酸和苏氨酸代谢；17.甘油磷脂代谢；18.半乳糖代谢；19.脂肪酸延伸；20.脂肪酸降解；21.脂肪酸生物合成；22.二帖生物合成；23.角质、木栓质和蜡质的生物合成；24.TCA循环；25.植物昼夜节律；26.碳代谢；27.油菜素类固醇生物合成；28.精氨酸生物合成；29.氨酰基tRNA生物合成；30.丙氨酸、天冬氨酸和谷氨酸代谢。

3 结论与讨论

镉作为一种重金属污染物，严重影响水稻的正常生长发育。水稻镉抗性机制十分复杂，因品种、遗传背景不同而存在差异。本研究通过对耐镉粳稻品种02428和镉敏感籼稻品种CH891进行镉处理，并对第3天的幼苗进行差异表达基因分析。结果表明，CH891和02428的代谢途径在镉处理3 a后存在差异，表明代谢途径的差异可能是籼稻和粳稻对镉胁迫反应不同的原因。因此，深入分析这些代谢途径中差异表达的基因对于解释籼稻和粳稻镉胁迫差异是必要的。

3.1 SCR3中的差异表达基因

差异表达基因的代谢途径富集分析表明，异黄酮生物合成和花青素生物合成途径是SCR3的主要富集通路。有研究表明，异黄酮在植物防御反应中起到重要作用[27]。在镉胁迫下，植物中的细胞壁受到破坏，使其易受到病原体的侵害，从而诱导了异黄酮的表达以抵抗病害[28]。本研究结果可能为水稻在异黄酮水平上对镉胁迫的响应提供新的见解。花青素是一种抗氧化剂，能清除植物细胞内的氧自由基，减轻活性氧的毒性[29]。Roychoudhury等[30]发现，花青素生物合成能被CdCl2诱导，与本研究结果相一致，表明花青素作为一种抗氧化剂，与水稻镉胁迫反应相关，可能减轻镉对水稻的毒害作用。

3.2 RCR3中的差异表达基因

就RCR3而言，不同水稻品种的KEGG代谢通路在镉胁迫后也存在不同。Bai等[31]研究发现，与核糖体相关的DEGs在冷应激信号转导中起着关键作用。核糖体途径在本研究中富集于RCR3中，表明核糖体途径很可能在水稻镉胁迫下被诱导并发挥作用。自噬可以将应激引起的氧化蛋白和损伤的胞内物质运输到空泡中降解，从而减少毒性物质在细胞质中的积累[32]，它在植物的发育和对胁迫的响应中起着重要的作用[33]。本研究发现，RCR3中诱导了自噬途径的调节，进一步验证了自噬调节在水稻镉胁迫中发挥了一定的作用。这2条通路是本研究在水稻镉胁迫下发现的未被大量报道的途径。

3.3 CCR3中的差异表达基因

代谢途径差异同样存在于CCR3中。α-亚麻酸代谢和脂肪酸降解可概括为脂肪酸代谢，它与CCR3中的镉反应有关。已有研究表明，重金属能特异性地改变植物体内脂肪酸代谢，镉能提高脂质过氧化物水平，导致水稻体内脂肪酸含量发生强烈变化[34-35]。本研究中，α-亚麻酸代谢和脂肪酸降解途径在CCR3中显著富集，进一步证明脂肪酸代谢可能是水稻对镉胁迫的应激反应之一。

综上所述，本研究表明，水稻镉响应机制存在品种的特异性。SCR3中富集的信号通路大多与次级代谢过程相关，而RCR3中则多为蛋白质代谢。这表明CH891可能通过诱导次生代谢过程对镉产生应答，从而增强其免疫力和逆境抗性。而02428中蛋白质的产生和降解过程受到影响，从而保护自身不受镉的侵害。