硒处理下冰菜的转录组响应及相关基因功能分析

姚运法，练冬梅，林碧珍，赖正锋，洪建基

(福建省农业科学院亚热带农业研究所，福建 漳州 363005)

【研究意义】冰叶日中花(MesembryanthemumcrystallinumLinn.)俗称冰菜，属于番杏科日中花属(MesembryanthemumL.)的一年生或二年生草本植物，其茎匍匐、叶肉质，具有较高食用价值[1]。原产非洲南部，世界各地均有引种栽培，近年来引种中国[2]。冰菜在环境适应性方面，具有耐干旱、盐碱、喜阳光、怕寒及水涝等特点，其适合生长的温度是5～25 ℃[3]。冰菜作为一种新型特菜资源，研究冰菜外源硒对冰菜生长和生理变化机制具有重要意义。【前人研究进展】硒是一种重要的微量元素，适量硒能清除植物体过量自由基[4]，促进根系生长发育[5]，有利于増强植物体新陈代谢，提高生命力，抵御逆境伤害，延缓组织成熟衰老等[6-8]；高浓度硒会促进植物过氧化作用，易造成植物体损伤[9]。研究表明，硒能改善多种作物的可食用部分品质[4]。蔬菜食用部位硒的含量浓度在0.0008～5.37 mg/kg范围内，均值为0.067 mg/kg[10]，其中叶菜中硒的平均浓度低于推荐值0.025 mg/kg[11]。蔬菜中的硒含量随着土壤中硒含量的增加而増加，对外源硒的响应能力因植物种类而异，薛瑞玲等[12]发现，低浓度硒盐能提高小白菜的抗氧化作用，促进叶绿素合成和叶片生长，高浓度时则相反。因此，适宜浓度的硒对植物的生长发育有益。【本研究切入点】近年来，在中国的冰菜设施栽培越来越广泛，开展冰菜富硒栽培模式具有重要经济价值。【拟解决的关键问题】本试验对冰菜以不同浓度(0、10和100 mmol/L)外源硒处理，利用Illumina Hi-seq 2500高通量测序技术，分析不同浓度硒处理下冰菜响应机制及其植物激素信号转导途径相关差异基因表达情况，筛选冰菜硒响应的相关基因及其表达机制，在转录组水平阐述冰菜对硒处理的响应机制，挖掘发现富硒的冰菜资源，为富硒冰菜育种积累理论基础。

1 材料与方法

1.1 供试材料

试验材料为冰菜，种植于福建省农业科学院亚热带农业研究所试验农场。2019年12月11日育苗，2020年1月10日定植。2020年2月15日，以亚硒酸钠为外源硒，不同浓度(M0: 0 mmol/L; M1:10 mmol/L；M2: 100 mmol/L)亚硒酸钠处理冰菜。2020年2月22日取嫩茎叶，每处理设置3个重复，共9个样品，置于-80 ℃超低温冰箱备用。

1.2 试验方法

1.2.1 转录组测序 利用Trizol(Invitroge)试剂分别提取硒处理(M1和M2)和对照(M0)冰菜嫩茎叶总RNA，采用Nanodrop(ThermoFisher)、Qubit 2.0、Aglient 2100(Bioanalyzer)技术检测RNA样品的纯度、浓度和完整性等，构建cDNA文库，再分别使用QubiM2.0、Agilent 2100和Q-PCR对文库的质量进行检测。合格后，基于SBS技术，使用Illumina Hiseq高通量测序平台对cDNA文库进行测序，产出高质量Raw Data，采用截除Reads中的测序接头以及引物序列方式，过滤低质量值数据，获得高质量Clean Data。利用Trinity软件对高质量的Clean Data进行序列组装。

1.2.2 测序数据组装与分析 将冰菜转录组测序Clean Data与组装得到的Transcript或Unigene库进行序列比对，比对到Transcript或Unigene的Reads称为Mapped Reads，Mapped Reads将用于后续的分析。

1.2.3 Unigene表达量估计与差异表达分析 采用Bowtie将测序得到的Reads与Unigene库进行比对，根据比对结果，结合RSEM进行表达量水平估计。利用FPKM值(Fragments Per Kilobase of transcript per Million mapped reads)表示对应Unigene的表达丰度。使用EBSeq进行差异表达分析。采用Benjamini-Hochberg方法对原假设试验得到的显著性P值进行校正，校正后P值，即伪发现率(False discovery rate)小于0.01且差异倍数(Fold Change,FC)大于等于2作为筛选标准。其中，FC表示两样品(组)间FPKM的比值。

1.2.4 DEGs筛选和韦恩图统计 根据对照与处理样品之间表达水平的相对高低，DEGs可分为下调基因(Down gene)和上调基因(Up gene)。然后对两组 (M0 与 M1; M0 与 M2) 数据DEGs进行维恩图分析。

1.2.5 Unigene功能注释与DEGs富集分析 使用Blast软件将Unigene序列与COG、GO、KEGG、KOG、Pfam、Swiss-Prot、eggNOG、NR等8大数据库比对，使用KOBAS 2.0得到Unigene在KEGG中的KEGG Orthology结果，预测完Unigene的氨基酸序列之后使用HMMER软件与Pfam数据库比对，获得Unigene的注释信息。系统分析基因产物代谢途径及功能，并将DEGs比对到KEGG数据库，得到DEGs的代谢途径。

1.2.6 硒处理对冰菜植物激素信号转导途径分析 对硒处理条件下冰菜植物激素信号转导途径DEGs(关键酶基因)进行转录组数据库检索，分析植物激素信号转导途径关键基因表达情况，研究不同硒处理浓度对植物激素信号转导途径中间产物影响。

1.2.7 植物激素信号转导途径关键差异基因验证 表1显示，取1 μg冰菜叶片总RNA，利用反转录试剂盒反转录成cDNA，利用Primer Premier 6软件，设计并合成引物，使用qRT-PCR分别检测M0、M1、M2处理与植物信号转导代谢途径关键差异表达基因，设置3个重复；统计9个基因(c38220为内参)在待测样品中的Ct值，计算相对表达量。

2 结果与分析

2.1 数据组装及分析

表2显示，经高通量测序和质量控制分析，共获得23.63 Gb Clean data，其中M0 Clean data为8.47 Gb，M1 Clean data为7.61 Gb，M2 Clean data为7.55 Gb，碱基百分比均达到92.0%以上。本次测序数据质量好，利于后续分析。由表3显示，转录本序列(Transcipt)组装出223 477条转录本序列，平均长度为1843.9 bp，N50长2916 bp；由单基因序列(Unigene)组装出76 806条，平长度为759.0 bp，N50长1513 bp。

表1 实时荧光定量PCR引物

表2 有效数据评估统计

表3 组装结果统计分析

表4 单基因序列注释统计表

表5 差异基因表达情况

2.2 Unigene功能注释

如表4所示，在76 806条单基因序列中，共获得28 172条Unigene的注释结果，占单基因序列总数36.68%。其中与COG数据库比对，获得9096条基因注释，占总注释基因的32.29%；与KOG数据库比对，获得16 137条基因注释，占总注释基因的57.28%；与GO数据库比对，获得10 493条基因注释，占总注释基因的37.25%；与KEGG数据库比对，获得15 636条基因注释，占总注释基因的55.50%；与Pfam数据库比对，获得18 324条基因注释，占总注释基因的65.04%；与Swiss-Prot数据库比对，获得17371条基因注释，占总注释基因的61.66%；与eggNOG数据库比对，获得25 240条基因注释，占总注释基因的89.59%；与NR数据库比对，获得26 405条基因注释，占总注释基因的93.73%。

2.3 DEGs的筛选与功能注释

2.3.1 DEGs筛选 如表5所示，通过DEGs的筛选，获得冰菜硒处理DEGs 570个(M0 与 M1)和580个(M0 与 M2)，其中M0 与 M1组表达量上调的基因493个，下调基因77个；M0 与 M2组表达量上调的基因293个，下调基因287个。通过对两组(M0 与 M1，M0 与 M2)DEGs维恩图分析，硒处理(M0 与 M1)非共表达基因468个，其中上调表达基因447个，下调表达基因21个；硒处理(M0 与 M2)非共表达基因478个，其中上调表达基因247个，下调表达基因231个；硒处理共表达基因102个，其中上调表达基因为46个，下调表达基因56个(图1)。

表6显示，将M0与M1和M0与M2两组DEGs单基因序列分别注释到COG、GO、KEGG、KOG、Pfam、Swiss-Prot、eggNOG、NR8个数据库，分别获得426和428个基因，其中COG数据库分别208和144个，GO数据库分别241和239个，KEGG数据库分别237和138个，KOG数据库分别309和168个，Pfam数据库分别366和330个，Swiss-Prot数据库分别289和325个，eggNOG数据库分别404和381个；NR数据库分别367个和422个，其中M0与M1组在eggNOG数据库注释比最高，超过94.8%；M0与M2组在NR数据库注释比最高，超过98.6%。

表6 差异表达基因注释数量统计

2.3.2 GO功能注释 如表7所示，对不同硒浓度处理M1、M2与M0(CK)DEGs进行GO分类统计，两组DEGs均归到43个功能小类。生物学过程中DEGs主要集中在代谢过程、细胞过程、单一生物过程、生物调节、刺激应答和细胞成分或生物合成等6个功能小类占比高，两组DEGs差值来看，单一生物过程和刺激应答生物学过程，随着硒处理浓度增加，DEGs分别增加29和22个；细胞过程和细胞成分或生物合成生物过程，随着硒处理浓度增加，DEGs分别减少45和24个。细胞组分过程中DEGs主要在细胞器、复杂大分子、细胞成分、细胞、细胞器成分、膜成分、膜结构7个功能小类占比高；两组DEGs差值来看：随着硒处理程度增加，细胞器、复杂大分子、细胞成分、细胞、细胞器成分、细胞组分过程DEGs减少，且DEGs均在50个以上；膜成分、膜结构和细胞器成分等过程比值较高，分别差42、35和18个DEGs。分子功能过程中DEGs在催化活性、结合活性、转运活性和分子结合活性4个功能小类占比最高。两组DEGs比值差异来看：其中催化活性过程与硒处理浓度呈正相关，DEGs个数达69个；结构分子活性过程与硒处理浓度呈负相关，DEGs数目达59个。

表7 硒处理差异表达基因GO功能注释比较

续表7Continued table 7

2.3.3 KOG功能注释 将M0与M1和M0与M2两组DEGs注释到KOG数据库，对注释结果进行直系同源分类，均获得22个功能分类，其中未知功能(S)、无机离子转运与代谢(P)、胞质运输、分泌和囊泡运输(U)、碳水化合物转运与代谢(G)、能源的产生和转化(C)、次生代谢产物的合成、转运和代谢(Q)、脂质转运与代谢(I)、翻译、核糖体结构和生物合成(J)、氨基酸转运与代谢(E)、转录(K)、无机离子转运与代谢(P)共11个功能分类DEGs差异较多。

两组KOG功能注释差值分析，一般功能预测(R)、次生代谢产物的合成、转运和代谢(Q)、碳水化合物转运与代谢(G)、氨基酸转运与代谢(F)功能分类DEGs数分别为34、14、13、12个；翻译、核糖体结构和生物合成(J)、能源的产生和转化(C)、无机离子转运与代谢(P)、细胞骨架(Z)、翻译后修饰、蛋白质转换、伴侣(O)功能分类DEGs数分别为70、27、13、10、10个(表8)。

表8 硒处理差异表达基因KOG功能注释比较

2.3.4 KEGG功能注释 表9显示，将两组DEGs与KEGG数据库比对，M0与M1和M0与M2分别有431和351条基因得到注释，分别富集在104和101条代谢通路，其中M0与M1组中，硒处理对核糖体(68个)、碳代谢(35个)、氧化磷酸化(30个)、氨基酸的生物合成(27个)、植物激素信号转导(26个)、苯丙烷类生物合成(19个)、植物—病原互作(18个)、光合作用的固碳作用(17个)、内质网蛋白加工(16个)、嘌呤代谢(15个)、丙酮酸代谢(15个)等11个代谢途径具有主要影响；M0与M2组中，硒处理对植物激素信号转导(34个)不饱和脂肪酸的生物合成(30个)、苯丙烷类生物合成(30个)、脂肪酸代谢(18个)4个代谢途径具有主要影响。M0与M1和M0与M2两组DEGs差值分析表明，不饱和脂肪酸的生物合成、脂肪酸代谢、苯丙烷类生物合成、植物激素信号转导代谢途径其DEGs增加分别为27、12、11和8个；核糖体、碳代谢、氨基酸的生物合成、氧化磷酸化、植物—病原互作、内质网蛋白加工代谢途径其DEGs分别减少66、32、25、21、9和9个。

表9 硒处理差异表达基因KEGG功能注释比较

2.4 植物激素信号转导途径分析

生长素(Auxin)信号转导途径中，生长素吲哚乙酸蛋白(Auxin/indole-3-Acetic Acid，AUX/IAA)有2个差异拷贝，M1基因表达无显著差异，M2处理均表现为下调表达；生长素上调小RNA(Small auxin up RNA,SAUR)共16个差异拷贝，M1和M2处理均具有上调和下调基因；生长素响应家族基因(Gretchen Hagen3，GH3)共有2个拷贝，M1和M2均表现下调(c32701)或上调(c32788)。细胞分裂素(Cytokinine，CK)信号转导过程中，磷酸转运蛋白(Arabidopsis histidine phosphotransfer proteins，AHP)(c31251、c36579)在M1和M2处理下均表现为上调表达，AHP(c39901)对M1处理表现下调，M2处理表现不敏感。A型拟南芥响应调节因子(Type-A Arabidopsis response regulators，A-ARR)对M1和M2处理显著上调作用。赤霉素(Gibberellin，GA)信号转导过程中，赤霉素不敏感矮秆蛋白1(Gibberellin insensitive Dwarf 1，GID1)(c39849)在M1和M2处理均表现为上调表达，转录因子(Transcription Factor，TF)(c37269)在M1和M2处理均表现为下调表达。脱落酸(Abscisic acid，ABA)信号转导途径中：蔗糖非发酵1相关蛋白激酶2(Sucrose nonfermenting 1-related protein kinase2，SnRK2)(c37312)在M1处理表现为下调表达，M2处理表现不显著表达，转导受体蛋白PYR/PYL/RCAR家族基因影响比较复杂，PYR/PYL/RCAR复合体(c26211)在M1处理表现为上调表达，M2处理表现为下调表达；PYR/PYL/RCAR复合体(c34313)在M1处理表达不显著，M2处理显著下调表达；PYR/PYL/RCAR(c28738)在M1和M2均表现上调表达；2C类蛋白磷酸酯酶(Protein phosphatase 2C,PP2C)有7个拷贝具有表达量差异表现，无规律；ABRE结合因子(ABRE binding factors，ABF)有2个差异表达拷贝，其中c34215在M1处理表现下调表达，M2处理表达差异不显著，c32284在M1处理表达差异不显著，M2处理表达显著下调。乙烯(Ethylene，ET)信号转导途径中，ETR和ERF1/2各有1个拷贝，在M1和M2处理下均表现为显著下调表达。茉莉酸(Jasmonic acid，JA)信号转导途径中，骨髓细胞瘤病蛋白2(Myeolcytomatosis proteins 2，MYC2)(c36305)和茉莉酮酸酯ZIM结构域蛋白(Jasmonate ZIM-domain，JAZ)(c28853)在M1处理下均表现显著上调表达，M2处理下表达无显著差异。水杨酸(Salicylic acid，SA)信号转导途径中，TGACG基序结合因子(TGACG motif-binding factor，TGA)和病程相关蛋白-1(Pathogenesis-related protein-1，PR-1)在M1处理条件下，表达量无显著差异表达，M2处理下，TGA表达量下调，PR-1表达量上调，上调表达量增加12倍。另外，油菜素类固醇(Brassinosteroid，BRs)信号转导途径在M1和M2处理下，均无显著差异表达基因(表10)。

表10 硒处理下冰菜植物激素信号转导途径关键基因差异表达

续表10Continued table 10

2.5 植物信号转导途径关键差异基因验证分析

对表10中部分DEGs进行荧光实时定量验证，共验证AUX/IAA、SAUR、GH3、AHP、MYC2、JAZ、PR-1等7个。以冰菜Actin (c38220) 基因为内参基因，进行qRT-PCR验证。由图2可知，7个DEGs基因qRT-PCR分析得到的相对表达量与转录组表达谱分析趋势完全一致，但表达的变化大小存在一定差异，说明基因表达谱的分析结果基本可靠，其中PR-1(c41466)基因分别M0和M2表达量差异极大，为冰菜水杨酸信号转导途径硒处理响应关键基因克隆提供研究基础。

3 讨 论

硒处理对植物激素信号转导途径中除油菜素类固醇外，其它7个代谢途径均产生重要影响。生长素IAA通过细胞质膜上生长素输入载体AUX1转运至细胞内，其受体TIR1与AUX/IAA互作，促进AUX/IAA泛素化而降解，解除生长素对ARF活性的抑制作用[13-14]。随着硒处理程度加深，冰菜叶片中AUX/IAA被泛素化降解，解除AUX对下游ARF活性的抑制作用，促进生长素信号转导；另外，SAUR和GH3对硒处理的响应较为复杂，且不具有规律性。胞外细胞分裂素与质膜上受体蛋白(Cytokinin response1,CRE1)结合，使得CRE1发生自磷酸化，然后其磷酸基团被传递至细胞质中的AHP上，磷酸化的AHP进入到细胞核中，再将磷酸基团转移至B型拟南芥应答调节子(B-Arabidopsis Response Regulators,B-ARR)，作为正调控转录因子，磷酸化的B-ARR激活下游基因A-ARR等表达，促进细胞分裂，A-ARR作为负调控因子，其累积反过来抑制B-ARR活性，从而抑制细胞分裂[15-16]。冰菜硒处理过程中，AHP对硒处理敏感，但对其浓度变化不敏感；A-ARR对硒处理显著上调作用，但对浓度变化不敏感；硒处理M1与M2处理均能显著性的提高A-ARR的表达量，显著的抑制细胞分裂素的细胞分裂作用，抑制冰菜生长和发育作用。赤霉素信号转导过程中，随着硒处理浓度增加而显著提高，GID1与DELLA蛋白结合，异构化DELLA蛋白空间构象，异构化后的DELLA蛋白在GID2作用下，进行泛素化修饰后被26S蛋白酶降解，其下游TF得到激活[17-18]，促进冰菜茎叶生长，但硒处理也对赤霉素信号转导途径下游TF(c37269)表达量显著抑制作用。脱落酸信号转导途径中，低浓度硒处理(10 mmol/L)导致SnRK2表达量显著降低，从而负调控ABA信号转导下游转录因子ABF，抑制气孔关闭，导致冰菜易脱水[19-20]；高浓度硒处理(100 mmol/L)SnRK2蛋白表达量显著恢复。另外，硒处理对ABA信号转导受体蛋白PYR/PYL/RCAR家族基因影响比较复杂，无显著规律。乙烯信号转导途径中，乙烯信号分子传到内质网乙烯受体蛋白(Ethylene receptor,ETR)上，ETR与下游乙烯抑制蛋白(Constitutive triple response 1,CTR1)互作，通过激酶磷酸化下游乙烯不敏感蛋白(Ethylene insensive 2,EIN2)，从而抑制EIN2下游的乙烯响应途径；此外，EIN2还受到EIN3结合F盒蛋白(EIN3-binding F-box protein1/2,EBF1/2)介导的泛素化蛋白酶降解调控，进而抑制ERF等下游乙烯响应基因[21-22]。随着硒处理浓度的提高，ETR和ERF1/2表达量降低，抑制EIN2和ERF下游的乙烯响应途径，抑制冰菜成熟和衰老[23]。茉莉酸信号转导途径中，JA在JAR1酶催化下将异亮氨酸(Isoleucine,Ile)加到JA上，获得活化性JA-Ile，活性JA-Ile大量积累，结合JA受体蛋白—冠毒素不敏感蛋白(Coronatine-insensitive protein 1,COI-1)，COI-1结合JAZ形成复合体，促进JAZ泛素化而降解，进而下游MYC2等转录因子释放，激活JA下游响应基因表达，调控植物非生物逆境响应[24]。低浓度硒处理(10 mmol/L)促进JAZ和MYC2表达量上调，调节硒对冰菜的逆境胁迫；高浓度硒处理(100 mmol/L)对JAZ和MYC2表达量影响不显著，表明高浓度硒处理对茉莉酸信号转导途径无显著影响。水杨酸信号转导途径中，SA与细胞质中受体蛋白—非发病相关蛋白表达子(Nonexpresser of Pathogenesis-Related genes,NPR1)受体结合，以寡聚体形式的NPR1解离为单体形式，单体NPR1进入细胞核，与转录因子TGA互作，共同激活下游基因表达，激活植物抗非生物逆境能力[25-26]。本实验中，TGA和PR-1在低硒处理(10 mmol/L)条件下，表达量无上下调表达，高硒处理(100 mmol/L)下，TGA表达量下调，PR-1表达量上调，其中PR-1上调表达量增加12倍，激活冰菜硒响应能力。经实时定量PCR分析，其中7个DEGs得到的相对表达量与转录组表达谱分析趋势完全一致。

4 小 结

利用高通量转录组测序技术，前期项目组初步对冰菜盐胁迫相关代谢分析和关键基因挖掘[27]，在此基础上，继续开展冰菜硒处理试验，研究植物信号转导途径响应硒处理的分子机制，为下一步研究硒调控过程的转录因子及其关键功能基因、信号转导途径提供参考。以上盐胁迫、硒处理逆境冰菜分子机制的研究，有助于冰菜抗逆品种的改良及选育，并指导富硒冰菜生产实践上具有重要意义。