** 同等贡献
以玉米种质POB21为材料, 利用数字基因表达谱技术对15% PEG渗透胁迫处理后的玉米叶片cDNA文库进行差异基因表达谱分析。结果表明, 转录组基因序列与参考基因组高度一致, 基因表达呈现出高度不均一性和部分冗余性。从中筛选出1097个有效差异表达基因, 其中795个表达上调, 302个表达下调。GO功能显著性富集分析表明, 3种细胞组分和分子功能中的3种糖基转移酶活性表现为富集。KEGG代谢分析表明差异表达基因广泛涉及糖、蛋白、核酸和脂类等生物大分子代谢、次生代谢、激素代谢以及能量代谢过程。脯氨酸代谢相关基因的差异表达分析表明, 谷氨酸途径是胁迫诱导脯氨酸积累的主要方式。表达谱分析结果为玉米渗透胁迫响应分子机制的深入研究和功能基因的筛选奠定了基础。
** Contributed equally to this work
In this study, a germplasm POB21 was used to analyze leaf cDNA library of maize treated with 15% PEG by
我国玉米种植面积仅次于美国, 其中约三分之二是旱作, 从播种到收获的全生育期都可遇到干旱。干旱脱水产生渗透胁迫可致玉米减产25%~30%, 严重时甚至绝收[ 1]。盐胁迫等其他非生物逆境也可以间接产生渗透胁迫[ 2]。提高玉米渗透胁迫抗性是改善玉米对干旱为主的多种非生物胁迫抗性的重要手段, 阐明玉米对渗透胁迫的抗性机制是抗旱育种和栽培技术改良的理论基础。
在渗透胁迫下, 水势下降引起植物各种盐离子和代谢物浓度相对升高, 造成细胞结构破坏、蛋白损伤和代谢异常, 同时也能引发植物保护机制的响应[ 3, 4]。植物根系在短期胁迫时生长加速, 而长期胁迫时生长受抑制; 地上部分由于气孔关闭和光合酶活性下降导致光合作用受抑制, 生长减缓。活性氧积累造成膜脂和某些功能蛋白的过氧化, 使膜脂饱和度增加、蛋白聚合和交联等结构异常。细胞内大部分多糖和蛋白合成受到抑制。在防御反应中, 植物通过积累K+、Ca2+、脯氨酸等渗透调节物质的来调节渗透平衡; 通过提高SOD、CAT、POD等酶活性来减轻活性氧伤害; 通过改变ABA等激素水平来调节基因表达。 ZmMKK4、 ZmbZIP72、 ZmLEA3等众多基因的表达水平在渗透胁迫下发生变化[ 5, 6, 7]。玉米在生理、生化和分子水平上对渗透胁迫的响应因胁迫强度、胁迫时间和遗传基础的差异表现出广泛的多样性。
以高通量测序技术为基础的基因表达谱分析是从转录组学水平上系统研究整体基因表达差异的有效方法。与差减杂交、基因芯片等差异表达分析方法相比, 高通量测序筛选效率高, 没有对基因组测序的依赖性, 同时还提供了基因片段的真实序列信息[ 8, 9]。近年来, 玉米基因组测序的完成和基因注释信息的不断丰富, 为玉米表达谱分析中的基因比对和筛选提供了有利条件[ 10]。本研究通过高通量测序分析玉米在渗透胁迫下的数字基因表达谱变化, 以阐述玉米响应渗透胁迫的分子基础, 为抗性基因的发掘和利用提供理论依据。
试验玉米材料为POB21改良热带晚熟马齿型白粒群体。取饱满一致的种子30粒, 6次重复, 播于沙床上, 25℃, 12 h d-1光照培养12 d, 将生长整齐一致的玉米幼苗分为2组, 一组加入15% PEG-6000溶液模拟渗透胁迫处理10 h, 另一组正常供水作为对照。渗透胁迫处理结束后取完全展开叶片混合, 于液氮中保存, 用于基因表达谱分析和差异表达基因的定量RT-PCR分析。
以TRIzol法提取10 µg以上的叶片总RNA, 加热打开二级结构, 用带有Oligo(dT)的磁珠富集mRNA, 加入fragmentation buffer使其成为短片段, 再以短片段mRNA为模板, 用六碱基随机引物合成cDNA第一链, 并加入缓冲液、dNTPs、RNase H和DNA polymerase I合成cDNA第2链, 用QiaQuick PCR试剂盒纯化, 进行末端修复, 加碱基A, 加测序接头, 以琼脂糖凝胶电泳回收目的片段后进行PCR扩增, 从而完成整个文库制备工作。用Agilent 2100 Bioanalyzer和ABI StepOnePlus Real-time PCR System检测文库的质量和产量, 文库质控合格后使用Illumina HiSeq 2000测序, 测序读长(read)为49 bp。测序仪产生的图像数据经Base Calling转化为raw reads序列数据, 去除接头序列、空载序列、含氮比例大于10%和低质量序列, 获得去除杂质的reads (Clean reads)。对原始Clean reads做标准化处理(每个基因包含的原始Clean reads数/该样本中总Clean reads数×1 000 000), 获得标准化的基因表达量TPM (transcript per million clean reads)。用reads比对软件SOAPaligner/SOAP2将Clean Reads分别比对到参考基因组和参考基因序列, 获得唯一比对上参考基因的reads (unique reads)。利用唯一比对上单一基因的reads数和唯一比对上参考序列的总reads数来计算基因表达量, 使用RPKM算法对测序深度和基因长度进行归一化。
B73参考基因数据库和参考基因组数据库的数据信息来源于MaizeGDB数据库(http://ftp.maize-sequence.org/current/)。差异表达基因的GO分析数据来自Gene Ontology数据库(http://www.eneontol-ogy.org/), 基因序列首先跟NCBI的Nr库用BLAST比对(参数: -p blastx -e 1e-5 -m 7), 再用Blast2GO软件(默认参数)注释到GO的各级term下。以校正之后的差异检验值p-value≤0.05为阈值, 满足此条件的GO term定义为在差异表达基因中显著富集。基因序列通过BLAST(参数: -p blastx -e 1e-5 -m 8)比对到KEGG数据库(kyoto encyclopedia of genes and genomes)[ 11]进行注释, Q value≤0.05的代谢途径定义为在差异表达基因中显著富集。
荧光定量PCR(qRT-PCR)中应用2×SYBRGreen master mix荧光染料。使用Beacon Designer7.5软件设计引物, 扩增长度为75~110 bp, 基因上下游引物分别为GRMZM2G125813: 5°-TCTTATCCATCTT GTGTTC-3°和5°-GACCTTAACAGGAATGAC-3°; GRMZM2G178074: 5°-TAATCCTGTGCTCGTCTTA-3°和5°-GTCTCTTCCTCAGCCATA-3°; GRMZM2G339327: 5°-CACTTGTAATGTTGATGT-3°和5°-ATATTGATG GAACGGTAT-3°。扩增程序为95℃预变性2 min, 95℃ 20 s, 58℃ 30 s, 72℃ 30 s, 共35个循环。 actin1基因为内参, 引物为5°-AATGACGCAGATTATGTT- 3°和5°-GAATCCATCACAATACCA-3°。
对测得的Reads分析表明, 对照样本(97.53%)和渗透胁迫处理样本(98.14%) Raw reads去杂后获得的Clean reads比例接近, 说明cDNA文库构建和测序质量较高。由于遗传背景的差异和注释信息量的限制, 对照样本Clean reads中近79.15%可被参考基因组注释, 20.85%的Clean reads未能注释到对应的基因(表1), 胁迫处理样本的reads注释情况相近。在可注释的reads中, 几乎全部reads序列与基因组序列差异都小于2个碱基, 其中大部分reads序列(对照样本73.8%, 渗透胁迫处理样本82.5%)可以与基因组完全匹配, 说明所用玉米材料POB21与参考基因组之间遗传基础差异较小。另一方面, 多种reads对应同一基因的情况较少, 大部分注释基因只有一种read (对照样本81.9%, 渗透胁迫处理样本84.3%)。此外, 测序随机性分析表明, reads在基因中的分布均一, 说明测序样品的mRNA完整性和片段化均一性很好, 序列分析中没有mRNA降解和偏向性的影响。
正常玉米叶片基因表达谱分析表明, 在2 410 116个共计27 665种unique reads中, 拷贝数(同一种read的测得数量)超过1000的高表达reads的种类虽然很少(0.99%), 但其表达量占mRNA总表达量的44.81%; 而拷贝数小于5的unique reads虽然表达量小, 只占mRNA总量0.89%, 但在种类上非常丰富, 占mRNA种类的40.23%。在渗透胁迫处理样品unique reads拷贝数分级中, 同样表现为>1000拷贝reads种类最少(0.93%), 而表达量最高(37.44%); 而<5拷贝reads表达量最低(0.89%)而种类最多(38.78%)。这说明细胞中基因表达表现出高度的不均一性和部分冗余性, 少数基因呈现高丰度表达, 而多数基因表达水平很低(表2)。
在进行差异表达基因分析过程中, 选取了参考基因39 354个, 可被注释的reads产生的有效参考基因22 996个, 占参考基因总数的87.7%。以差异检验FDR值≤0.001且差异倍数不低于2倍为差异表达基因(DEGs, differentially expressed genes)的选择标准, 渗透胁迫处理样本与对照样本之间有1097个基因表现为差异表达, 其中302个DEGs在渗透胁迫条件下表达下调, 795个DEGs表达上调, 后者明显高于前者(图1)。
由图2可以看出, 基因表达差异高度不均一, 表达水平变化主要集中在10倍以内, 占差异表达基因总数的89.15%。其中, 表达上(下)调2~5倍的基因占差异表达上(下)调基因总数的71.82% (82.45%), 表达上(下)调5~10倍的基因占14.97% (12.91%), 表达上(下)调10~100倍的基因占10.19% (4.30%), 表达上(下)调100倍以上的基因占3.02% (0.33%)。相比而言, 上调基因的表达增强幅度更大, 表达上调10倍以上的基因所占比例是表达下调10倍以上基因所占比例的2.85倍。
在GO功能显著性富集分析中, 通过与参考基因比较, 从GO数据库中筛选出差异表达基因的GO分类条目, 可以寻找与基因表达差异相关的主要特征功能分类。在从数据库中查询到校正 P≤1.0的496个GO分类信息条目中, 细胞位置组分(cellular component)、基因的分子功能(molecular function)和参与生物过程(biological process)三类条目分别占31.82%、11.87%和56.31%, 校正 P值小于0.5的部分GO分类条目见图3。以校正 P≤0.05为显著性富集标准, 在细胞组分分类信息中, 有3个分类条目达到显著性富集, 分别为膜内在组分(GO: 0031224)、膜组分(GO: 0044425)和类囊体组分(GO: 0044436)。在分子功能分类信息中, 也有3个条目达到显著性富集, 分别为糖基转移酶活性(GO: 0016757)、葡萄糖基转移酶活性(GO:0046527)和己糖基转移酶活性(GO: 0016758), 说明渗透胁迫下玉米叶片中的糖基转移反应受影响较大。在生物过程分类信息中, 未获得显著性富集的条目。
在渗透胁迫下高水平差异表达100倍以上的25个注释基因中, 只有GRMZM2G092747一个假定蛋白(hypothetical protein)基因表达下调(表3)。在表达上调基因中, 除8个假定蛋白基因外, 其他16个基因的功能涉及以下3个方面: (1)基础代谢。参与糖代谢的差异表达基因有4个, 分别编码呋喃果糖苷酶(GRMZM2G089836)、含淀粉结合结构域的CBM家族蛋白(GRMZM2G161534)、棉籽糖合酶(GRMZM2 G340656)和PEP羧激酶(GRMZM2G178074)。参与脂代谢的差异表达基因有3个, 即脂基转移蛋白(GRMZM2G058208)、甘油-3-磷酸转运蛋白(GRMZM 2G124136)和羧酸酯酶(GRMZM2G391795)。另外, EDL3蛋白(GRMZM2G389301)参与调节蛋白降解。(2)转录因子。包括NAC转录因子(GRMZM2 G146380)、AP2/EREBP基因家族(GRMZM2G319786)和HD-ZIP因子(GRMZM2G117164)。(3)其他逆境响应因子。有蛋白磷酸化酶PP2C(GRMZM2G159811)、BAX1抑制因子(GRMZM2G014672)、ABC(ATP- binding cassette)转运体(GRMZM2G167658)、生长素响应因子(GRMZM2G460861)和功能未知的胁迫响应膜蛋白(GRMZM2G179462)。
在KEGG数据库中对差异表达基因进行功能注释和分类, 通过显著性富集能确定差异表达基因参与的主要生化代谢和信号转导途径。结果表明, 差异表达基因分布在110个分类代谢途径中。这些差异表达基因广泛涉及基础物质代谢、次生代谢、能量代谢、激素代谢及信号转导途径(图4)。受KEGG数据库注释信息量的限制, 只有68.09%的差异表达基因可以详细注释到具体的代谢路径图中。其中与糖类、蛋白(含氨基酸)、脂类和核酸等4类生物大分子代谢相关的基因共占差异表达基因总量的32.19%, 次生代谢相关的基因占16.71%, 激素代谢相关基因也较多(7.88%), 氧化磷酸化、光合作用和细胞活动相关基因所占比例相对较低。在 Q≤0.05为标准的显著性富集分析中, 达到显著水平的代谢途径可以分成以下5类: (1)次生代谢。包括类黄酮、花色素、生物碱、萜烯类次生代谢物等的合成和降解。(2)激素代谢和信号转导途径。包括传统五大类激素、油菜素内酯及茉莉酸代谢。(3)转运体代谢。多种ABC转运体代谢发生改变。(4)光合作用。如天线蛋白、叶绿素和光合电子传递体代谢。(5)渗透调节途径。如糖代谢和脯氨酸合成代谢途径。
脯氨酸是渗透胁迫下积累的重要渗透调节物质, 高等植物体内脯氨酸合成根据底物不同主要有谷氨酸途径和鸟氨酸途径两条途径。差异表达基因的KEGG代谢途径分析表明, 在精氨酸和脯氨酸代谢途径(ko00330)中, 有5个脯氨酸代谢相关基因的表达超过2倍以上(表4)。吡咯啉-5-羧酸还原酶(P5CR)和吡咯啉-5-羧酸合成酶(P5CS)是谷氨酸向脯氨酸转化的2个关键酶, 这两个基因均被胁迫诱导表达, 且限速酶P5CS上调达10倍以上; 而逆向代谢中的脯氨酸脱氢酶基因表达下调4.29倍。鸟氨酸途径的关键酶鸟氨酸转氨酶(OAT)基因表达没有检测到明显变化, 谷氨酸-N-乙酰转移酶基因表达下调, 说明谷氨酸向鸟氨酸转化过程受到抑制。精氨酸脱羧酶基因的表达上调暗示精氨酸向鸟氨酸的转化也间接受到抑制。这些脯氨酸代谢相关基因的差异表达表明, 在玉米响应渗透胁迫过程中, 谷氨酸途径可能是脯氨酸积累的主要合成途径。
随机筛选3个差异表达上调基因GRMZM 2G125-813、GRMZM2G178074、GRMZM2G339327进行qRT-PCR验证。玉米幼苗按表达谱胁迫条件进行PEG渗透胁迫处理后提取总RNA, 反转录后进行定量PCR扩增。结果表明, 3个基因在胁迫处理后发生不同程度的表达(图5)。尽管qRT-PCR分析得到的表达差异与表达谱分析的表达差异倍数不一致, 但胁迫诱导表达的变化趋势相同, 说明基因表达谱的分析结果可靠。
在干旱、盐碱等产生的渗透胁迫下, 植物通过改变基因的表达, 调控不同代谢途径和信号转导途径, 从转录和翻译等不同水平作出响应。通过高通量测序产生的数字基因表达谱分析获得了转录水平上32 208个基因的真实序列信息, 通过其拷贝数的不同直接反映相应基因的表达差异, 筛选出1097个差异表达基因。尽管由于read测序长度较短, 大多数注释基因的覆盖度较低, 对照样本和胁迫处理样本分别有近66%和75%的注释基因的覆盖度低于一半, 覆盖度超过80%的注释基因只有10%左右。但对于参考已知的基因组进行基因注释, 基本能够满足要求。尽管由于基因型的差异, 测序所得reads中有近1/5与参考基因组差异较大, 但有近60%的reads与参考基因组序列完全相同, 说明物种内遗传基础相对比较保守。
在测序所得基因中, 虽然渗透胁迫下只有3.41%的基因的表达水平发生变化, 但这些差异表达基因表现出广泛的功能多样性。在GO功能分析中, 差异表达基因涉及线粒体、叶绿体、细胞核和染色体、核糖体等几乎所有细胞组分, 包含了各种转移酶、水解酶、氧化还原酶、合成酶和异构酶等各种酶活性, 参与四类生物大分子和各种次生物质代谢、能量转换、跨膜运输、信号传递及其调节等各种生物过程。这说明植物转录系统是相对稳定的, 同时又是高度灵敏的一个开放式表达调控体系。
在渗透胁迫下高水平诱导表达(>100倍)的基因中, 参与糖代谢、脂类代谢和蛋白代谢等基础代谢过程的基因占很大比重。含有淀粉结合结构域的CBM家族蛋白、棉籽糖合酶、呋喃果糖苷酶和PEP羧激酶基因的表达上调可以改变各种可溶性糖组成和含量, 是植物响应胁迫进行渗透压调节的一种重要方式[ 12, 13, 14, 15]。脂基转移蛋白通过调节脂质合成和运输影响植物逆境适应性, 如水稻超表达 OsDIL基因后抗旱能力明显增强[ 16, 17]; 甘油-3-磷酸转运蛋白属于MSF (major facilitator superfamily)家族成员, 可以通过调节磷酸甘油向细胞质的运输影响磷脂代谢和糖酵解[ 18]; 羧酸酯酶通过调节脂质水解水平改变多种生物和非生物胁迫抗性[ 19]。作为泛素蛋白降解系统中的SCF泛素连接酶复合物组分, EDL3蛋白通过调节蛋白降解参与ABA依赖的干旱、盐等逆境信号转导过程[ 20]。其次, NAC等3种高水平表达上调的转录因子都是典型的逆境胁迫响应因子。NAC转录因子是植物特有的一类转录因子, 不仅参与植物生长发育的调控, 而且在各种植物抗逆反应中具有重要的调控作用, 如玉米 ZmSNAC1基因能显著增强转基因拟南芥的脱水耐性, 水稻 SNAC基因在小麦中超表达也可以提高其干旱和盐抗性[ 21, 22]。AP2/EREBP家族转录因子在植物中分布广泛, 在调控植物发育和低温、干旱、盐等胁迫抗性中发挥重要作用, 如过量表达 OPBP1基因可以显著提高烟草的抗病和抗盐能力, 而棉花 GhDREB基因过量表达可以提高小麦的干旱、高盐和冻害等多种胁迫抗性[ 23, 24, 25]。HD-ZIP是响应渗透胁迫的代表性转录因子, 在玉米基因组中至少有四类共55个HD-Zip基因, 其中17个玉米I类HD-Zip基因中有12个受干旱诱导上调[ 26]。在向日葵中过量表达 HaHB1和 AtHB13两种HD-Zip基因可以通过诱导表达膜稳定蛋白从而提高其对干旱和盐的抗性[ 27]。此外, 蛋白磷酸化酶PP2C是各种逆境胁迫信号激活的激酶信号通路中一种代表性的调控因子, 如拟南芥 AtPP2CG1过量表达后可以提高耐盐能力, 玉米 ZmPP2C2基因在烟草中超表达可以提高其低温抗性[ 28, 29, 30]。BAX1抑制因子是内质网胁迫响应标记基因, 可以抑制各种生物或非生物胁迫导致的细胞程序化死亡[ 31]。ABC转运体可以通过改变植物激素和次生代谢物运输、气孔调节等方式影响植物逆境胁迫抗性, 如拟南芥超表达 AtABCG36基因可以提高干旱和盐等胁迫抗性[ 32]。在玉米基因组中至少有130个ABC转运体基因[ 33]。上述基因的功能分析表明, 胁迫响应基因具有广泛的功能多样性。另一方面, 其中大多数基因可以参与多种生物或非生物逆境响应过程, 也说明不同逆境胁迫响应途径有一定的共通性。
植物可以通过多种代谢调节增加对逆境的抗性。在玉米响应渗透胁迫的代谢响应过程中, 这些防御机制主要表现为以下几个方面: (1)通过激素受体及其调控的转录因子的差异表达, 系统调节植物激素平衡[ 34]。生长素调控转录因子ARF、细胞分裂素受体CRE1、ABA调控转录因子ABF、赤霉素受体GID1和油菜素内酯受体BRI1基因表达上调, 而生长素受体AUX/IAA、ABA受体PYR/PYL、赤霉素调控DELLA蛋白、乙烯调控转录因子ERF1和茉莉酸调控转录因子MYC2基因表达抑制, 说明渗透胁迫响应是受多种激素共同调控的结果。(2)渗透调节剂的积累是植物提高渗透胁迫逆境适应能力最直接的方式。吡咯啉-5-羧酸还原酶和吡咯啉-5-羧酸合成酶基因的表达上调说明脯氨酸在胁迫诱导下快速积累[ 35]。淀粉酶、β-呋喃果糖苷酶、β-葡萄糖苷酶、β-D-木聚糖合成酶、蔗糖合成酶、乙醛还原酶、半乳糖醛酸转移酶等基因表达上调, 说明渗透胁迫诱导多糖降解, 而棉籽糖、蔗糖、果糖、甘露糖、阿拉伯糖、木糖等可溶性寡糖和单糖大量积累。脯氨酸和可溶性糖的积累都有利于调节胁迫下渗透压的稳定[ 36]。(3)通过提高抗氧化能力减轻胁迫产生的活性氧伤害[ 37]。过氧化氢酶(CAT)等抗氧化酶及抗坏血酸、维生素E、类黄酮、GSH等非酶抗氧化剂合成代谢途径中的关键酶基因均表现为诱导表达, 说明渗透胁迫下多种抗氧化机制协同发挥作用。(4)通过蛋白质量控制减轻胁迫伤害[ 38]。HSP、DnaJ和Clp等众多分子伴侣蛋白基因的诱导表达说明蛋白损伤修复系统能够对渗透胁迫作出快速响应, 而EDL3等蛋白降解系统的激活可以及时清除损伤蛋白。
上述代谢响应的多样性说明玉米对渗透胁迫的应答是多种代谢途径和信号通路共同作用的结果。要对玉米的渗透胁迫抗性进行遗传改良, 不是调控单个基因的功能可以独立完成的, 这可能也是目前通过基因工程手段进行玉米抗旱、抗盐等遗传改良不能进入生产应用的主要限制因素。通过玉米响应渗透胁迫的数字基因表达谱分析, 有助于系统理解渗透胁迫条件下各代谢途径的相互作用关系, 为阐述玉米渗透胁迫抗性形成的分子机制, 筛选主效基因提供了理论基础。
通过玉米渗透胁迫表达谱分析鉴定出1097个差异表达基因, 其中795个在渗透胁迫下表达上调, 302个在渗透胁迫下表达下调, 在表达上调基因中82.45%表达上调2~5倍。差异表达基因广泛涉及糖、蛋白、核酸和脂类等生物大分子代谢、次生代谢、激素代谢和能量代谢等过程。谷氨酸途径是渗透胁迫下脯氨酸积累的主要合成途径。
[1] |
|
[2] |
|
[3] |
|
[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] |
|