Transcriptome Analysis of Leaves of Seven Hemiptelea davidii Clones under Salt Stress
-
摘要:目的
从转录组水平上探究刺榆叶片对盐胁迫的响应机制，筛选耐盐基因。
方法本研究以7个不同耐盐能力刺榆无性系的扦插苗为试验材料，对0、100、150、200 mmol·L−1 NaCl溶液处理15 d的叶片进行了转录组测序与分析。
结果（1）除DW1上调数量略多于下调数量外，其他无性系核心差异基因下调数量均多于上调数量，且各无性系在200 mmol·L−1盐溶液处理时差异表达基因数量最多，对盐胁迫的响应最为强烈。GO富集分析显示细胞壁等功能类别基因对盐胁迫的响应较强。KEGG富集分析显示植物激素信号转导等代谢通路在7个无性系响应盐胁迫时发挥了重要作用。（2）强耐盐型的无性系DJ1、DJ3与弱耐盐型的无性系DW1、DY1在200 mmol·L−1盐溶液处理时的4个对比组合共有的差异基因数量为1 406个，共同上调基因393个、共同下调基因996个，其中脱落酸调节因子蛋白磷酸酶PP2C基因AHG1、ABI2，乙烯受体基因ERS1、ETR2，以及环核苷酸门控离子通道基因CNGC10和CNGC15可能与DJ1、DJ3较强的耐盐能力有关。
结论盐胁迫引起刺榆转录水平上功能类别基因与代谢通路的响应，其中6个显著差异表达基因将作为候选耐盐基因。本研究为揭示刺榆响应盐胁迫的转录组特征与筛选耐盐基因提供了理论参考，为进一步开展刺榆种质资源的研究与应用奠定了基础。
Abstract:ObjectiveTo investigate the response mechanism of Hemiptelea davidii leaves under salt stress at the transcriptome level and screen for salt tolerance genes.
MethodsBased on seven different H. davidii cuttings, transcriptome sequencing and analysis were performed on leaves treated with 0, 100, 150, and 200 mmol·L−1 NaCl solution for 15 days.
Results(1) The screening results of core differentially expressed genes（DEGs） in seven clones showed that except for the DW1, which had slightly more upregulated genes than downregulated genes, all other clones had more downregulated genes than upregulated genes. Moreover, when treated with 200 mmol·L−1 salt solution, each strain had the highest number of DEGS and the strongest response to salt stress. GO enrichment analysis showed that functional genes such as cell wall, had strong responses to salt stress. KEGG enrichment analysis showed that metabolic pathways such as hormone signal transduction, played important roles in the response of 7 clones to salt stress. (2) The four comparative combinations of stronger salt tolerant clones DJ1 and DJ3 and weak salt tolerant clones DW1 and DY1 had a total of
1406 DEGs, including 393 upregulated genes and 996 downregulated genes. Among them, the abscisic acid regulatory factor protein phosphatase 2C (PP2C) genes AHG1, ABI2, ethylene receptor genes ERS1, ETR2, and cyclic nucleotide gated channel genes CNGC10 and CNGC15 might be related to the strong salt tolerance of DJ1 and DJ3.ConclusionSalt stress causes the response of functional genes and metabolic pathways in H. davidii. Six significantly DEGs will be served as candidate salt tolerant genes. This study provides a theoretical reference for revealing the transcriptome characteristics of H. davidii in response to salt stress and screening salt tolerant genes, laying a foundation for further research and application of H. davidii germplasm resources.
-
Keywords:
- Hemiptelea davidii /
- salt stress /
- differentially expressed genes /
- qRT-PCR /
- salt tolerance genes
-
盐渍土在我国分布广泛，其类型多、面积大，且在不断扩张，严重威胁生态安全与农业可持续发展[1-3]。诸多研究表明，植物具有改良土质、改善生态环境、实现有效治理与开发利用盐渍化土地的重要途径[4-5]，开展耐盐植物种质资源的挖掘、筛选与应用研究对耐盐植物栽培利用和盐碱地改良具有重要意义。植物在受到盐胁迫时首先会在基因和转录组水平上进行调控，进而控制代谢以应对逆境。转录组学研究技术可以从植物整体转录水平揭示植物体中复杂的生物学通路和调控网络分子，对研究植物响应盐胁迫的复杂机制具有重要意义[6]，目前在金丝楸（Catalpa bungei C.A.Mey）[7]、多枝柽柳（Tamarix ramosissima Ledeb.）[8]、平邑甜茶（Malus hupehensis Rehd.）[9]等树种响应盐胁迫机理研究中得到较广泛的应用。
刺榆（Hemiptelea davidii (Hance) Planch.）为落叶小乔木，榆科（Ulmaceae）刺榆属（Hemiptelea）单种，耐瘠薄、耐低温、适应性强，在我国西北、西南、东北及华东等地分布较广泛，其中在山东主要分布在济南南部山区、淄博鲁山、胶东山区丘陵等地[10]，近年来的调查发现在刘公岛等近海岛屿中也有分布。目前，刺榆的研究主要集中在种苗繁育与生态作用方面，抗逆研究报道较少。其中，耐盐性研究报道仅见杨颖[11]对刺榆种子萌发及其幼苗进行的盐胁迫研究，该研究显示在NaCl<30 mmol·L−1的盐胁迫下，刺榆种子能够正常发芽，幼苗内MDA含量和SOD、POD酶活性未发生显著变化，表明刺榆种子萌发与幼苗生长对低浓度盐环境具有一定的耐受性，为刺榆在盐碱地的种植和推广提供了理论依据。刺榆响应胁迫的分子机制研究报道亦较少，无性系响应盐胁迫的转录组学特征与耐盐基因筛选目前尚未见报道，刘华波[12]利用干旱模拟试验和转录组数据探讨了刺榆核质基因响应干旱胁迫的分子机制，研究表明核质基因中光合作用相关通路、植物激素信号转导通路和次生代谢物合成通路受干旱胁迫影响最大，筛选并初步证实TrxG-like 基因参与刺榆干旱胁迫响应、提高植株的抗旱性。本研究以7个刺榆无性系为试验材料，对不同浓度NaCl溶液处理下各无性系叶片的转录组学特征开展研究，以期为刺榆对盐胁迫的响应机制研究及耐盐基因的筛选与利用提供参考。
1. 材料与方法
1.1 试验材料与处理
试验于山东省林草种质资源中心温室中进行。7个供试刺榆无性系（DJ1、DJ2、DJ3、DJ4、DW1、DY1、DH1）为2021年春季采集野生母株1年生枝条进行扦插繁育得到的无性系，其中DJ1～DJ4采自济南市柳埠街道锦阳川、DW1采自潍坊仰天山、DY1采自烟台昆嵛山、DH1采自威海刘公岛，前期研究表明DJ1与DJ3为强耐盐型，DJ2为较强耐盐型，DH1与DJ4为中等耐盐型，DW1与DY1为弱耐盐型。
试验苗以泥炭土、蛭石、珍珠岩按10∶4∶1的体积比混匀配制的营养土为基质，栽植于直径18 cm、高度15 cm的花盆中。2022年7月，将根系发育良好、长势基本一致的扦插苗置于温室内，室内温湿度与自然条件基本一致：试验期间实测夜间温度19～22 ℃、相对湿度73%～90%，日温23～32 ℃、相对湿度50%～74%。
2022年8月，参照杨颖的工作[11]及在前期预试验的基础上，设置清水（0 mmol·L−1）为对照，50、100、150、200 mmol·L−14种浓度NaCl溶液。为避免产生盐激效应，在试验开展前先用100 mmol·L−1 NaCl溶液预处理150、200 mmol·L−1两个梯度的试验苗3 d，预处理后与其他浓度一起开始盐胁迫处理。处理时，花盆置于托盘上，每盆分2次浇施1 L相应浓度盐溶液，浇施的盐溶液从花盆自然渗出到托盘中，保留渗出的盐溶液使其浸泡花盆1夜（12 h）以使栽培基质充分吸收盐溶液，次日取出后用塑料薄膜包裹花盆减少盆口水分蒸发，每5 d处理1次，各处理重复3次，共处理105株。试验处理至15 d时，各无性系叶片出现不同程度的盐害症状，即进行盐害情况调查，并取相同叶位的叶片，去离子水冲洗后用滤纸吸干水分，采集后立即放到液氮中速冻，然后置于−80 ℃冰箱中保存，用于转录组测序。
1.2 指标测定方法
1.2.1 盐害情况调查
参照白榆耐盐性评价研究中对盐害程度分级的方法[13]，根据刺榆叶片损伤症状，将刺榆盐害分为0、Ⅰ、Ⅱ、Ⅲ4个等级（图1），盐害等级评价标准为：0级：叶片无盐害症状；Ⅰ级：叶缘轻度萎蔫，萎蔫面积不足1/5；Ⅱ级：叶缘卷曲干枯，叶片萎蔫面积1/5～1/2；Ⅲ级：叶缘卷曲干枯，叶片萎蔫面积1/2以上至干枯脱落。
![]()
图 1 刺榆不同盐害等级的叶片症状注：从左到右依次为盐害等级0级、Ⅰ级、Ⅱ级、Ⅲ级Figure 1. Symptom of different salt damage levels on Hemiptelea davidii leavesNote: From left to right, the levels of salt damage are 0,Ⅰ,Ⅱand Ⅲ1.2.2 转录组测序和分析
使用Trizol试剂（Invitrogen，CA，USA）提取样本总RNA，利用Bioanalyzer 2100和RNA 1000 Nano LabChip Kit （Agilent,CA,USA）对提取的样本总RNA质量进行检测，质检合格后，利用附着有poly-T寡核苷酸的磁珠纯化mRNA，然后使用Fragmentation Buffer对其进行片段化，利用mRNASeq sample preparation kit （Illumina,San Diego,USA）逆转录裂解的RNA片段创建最终的cDNA文库，文库质检合格后采用Illumina Novaseq™6000进行测序，测序读长为双端2×150 bp。
在去除测序接头、低质量等不合格序列后，使用FastQC验证序列质量，包括有效数据的Q20、Q30和GC含量，所有的下游分析均基于高质量的有效数据。使用Trinity 2.4.0进行转录组的从头组装，归一化后得到非重复基因序列（unigenes）[14]。使用DIAMOND将组装完的unigenes与六大数据库（NCBI_NR、GO、Swiss-Prot、KEGG、Pfam、eggNOG）进行比对、分类和相应的功能注释，E-value阈值小于0.000 01[15]。
通过Salmon计算样品基因TPM值以评估unigenes的表达水平。使用R包edgeR对样本之间进行显著差异分析，将|log2 fold change| > 1，且P值＜0.05的unigenes作为差异表达基因，并对其进行GO和KEGG富集分析[16-18]。
1.2.3 实时荧光定量PCR分析
随机选择4个差异表达基因进行实时荧光定量PCR分析，进一步验证转录组数据的可靠性。利用Trizol法提取样品的总RNA，使用NovoScript® Plus All-in-one 1st Strand cDNA Synthesis SuperMix （gDNA Purge）反转录试剂盒将RNA反转录为cDNA。以较为稳定的GAPDH作为内参基因[19]，按照TB Green® Premix Ex Taq™ （Tli RNaseH Plus）荧光定量试剂盒的操作方法，进行实时荧光定量PCR，测定样品中相关基因的表达量，按照
$2^{-\Delta \Delta C_T} $ 法[20]计算基因的相对表达量，基因扩增引物序列见表1。表 1 实时荧光定量PCR引物序列Table 1. Primer sequences for qRT-PCR基因名称与方向
Unigene name and direction引物序列
Primer sequence（5′-3′）TRINITY_DN61310_c1_g2F CGGAGACTAAGGCCAACT TRINITY_DN61310_c1_g2R GCTTTCCTGGCAACTGTACTG TRINITY_DN38595_c1_g1F CGCTCAGATGGGACTTTCAAC TRINITY_DN38595_c1_g1R GAGGAGACCAGTAGATCGATCC TRINITY_DN63193_c0_g2F TGTGAGAGTCCAGGCTGC TRINITY_DN63193_c0_g2R CCTTCCTCGATGTTCACCTC TRINITY_DN65673_c0_g2F TCTGTCCTTCGTGCTCAATTACC TRINITY_DN65673_c0_g2R CACCTACCACTCCAGAGCACTT GAPDH-F TCAACAATGCCAAACCTG GAPDH-R GTGTCAACGAGCACGAAT 2. 结果与分析
2.1 盐胁迫对各株系叶片表型的影响
如图2所示，在50 mmol·L−1浓度处理下各无性系均未表现出盐害症状；在100 mmol·L−1盐溶液处理时，各无性系均表现为Ⅰ级盐害症状；在150 mmol·L−1盐溶液处理时，各无性系进一步表现出受害症状，且受盐害程度出现较大分化，其中DJ1、DJ2、DJ3、DJ44个无性系盐害程度为Ⅰ级，DW1、DH1盐害程度为Ⅱ级，DY1无性系受害情况最重、达到Ⅲ级；在200 mmol·L−1盐溶液处理时，各无性系健康分值均为处理间最低，不同无性系的受害程度也进一步分化，其中DJ1无性系、DJ3无性系盐害程度虽然加重但仍为Ⅰ级，DJ2、DJ4、DH1 3个无性系盐害程度为Ⅱ级，DW1、DY1两个无性系受害情况最严重、达到Ⅲ级。
![]()
图 2 不同浓度盐溶液处理下刺榆各无性系中上部成熟叶片盐害情况Figure 2. Salt damage performance of the upper part leaves of Hemiptelea davidii clones under different salt stress2.2 测序质量
测序得到的原始序列（Raw Reads）数量介于39 993 532条～44 393 618条，共1 169 469 262条，过滤后得到的有效序列（Valid Reads）数量介于37 318 336条～43 091 940条，共1 123 282 434条，有效序列占比介于92.28%～97.07%，平均占比96.05%；原始碱基（Raw Bases）数量6.00 G～6.66 G，共175.42 G，有效碱基（Valid Bases）数量5.19 G～5.99 G，共156.03 G；Q20介于96.64%～97.14%，平均97.01%，Q30介于90.82%～91.88%，平均91.62%；GC含量介于43.61%～44.61%，平均44.14%。有效数据、Q20、Q30及GC含量表明测序质量可靠，可进行深层次分析。
2.3 基因组装结果
将获得的有效数据使用Trinity 2.4.0进行组装，共组装转录本174 852条、碱基192 456 811个，其长度介于201～16 721 bp，中值为715 bp，200～500 bp占比为39%，GC含量为39.97%，N50为
1799 bp；归一化后共得到基因87 712条、碱基78 334 647个，长度范围与转录本相同，中值为483 bp，200～500 bp占比达51%，GC含量为40.03%，N50为1590 bp，组装结果较好。2.4 基因功能注释结果
使用DIAMOND对组装好的Unigenes进行功能注释，6种数据库的注释结果统计见表2，各数据库注释的基因数量由多到少依次为eggNOG、GO、Pfam、Swissprot、NR、KEGG，总注释87 712个基因。
表 2 注释结果统计Table 2. Annotation result statistics数据库
Database基因数量
Number of Unigenes比率/%
RatioGO 49 522 56.46 KEGG 31 171 35.54 Pfam 45 055 51.37 Swissprot 42 904 48.91 eggNOG 55 769 63.58 NR 42 350 48.28 Total Unigenes 87 712 100.00 2.5 差异表达基因分析
2.5.1 各无性系核心差异表达基因筛选与分析
2.5.1.1 各无性系核心差异表达基因筛选
为明确获得刺榆各株系响应盐胁迫的差异表达基因，对各无性系中与对照组差异明显的100、150、200 mmol·L−13个处理浓度（依次标记为2、3、4）与对照0 mmol·L−1（标记为1）进行比较，以|log2 fold change| > 1且 P值 < 0.05为标准筛选差异表达基因，将3个对比组共有的差异基因作为该无性系核心差异表达基因。
如图3（A-G）所示，在DJ1株系共筛选出703个核心差异基因，其中在3个处理浓度对比组中均上调的基因184个、均下调的514个；DJ2株系共筛选出430个核心差异基因，其中在3个处理浓度对比组中均上调的基因154个、均下调的263个；DJ3株系共筛选出
1021 个核心差异基因，其中在3个处理浓度对比组中均上调的基因178个、均下调的830个；DJ4株系共筛选出960个核心差异基因，其中在3个处理浓度对比组中均上调的基因166个、均下调的780个；DW1株系共筛选出702个核心差异基因，其中在3个处理浓度对比组中均上调的基因351个、均下调的319个；DY1株系共筛选出1757 个核心差异基因，其中在3个处理浓度对比组中均上调的基因629个、均下调的1097 个；DH1株系共筛选出323个核心差异基因，其中在3个处理浓度对比组中均上调的基因108个、均下调的194个。除DW1株系上调基因数量略多于下调基因数量外，其他无性系下调基因数量均多于上调基因数量。此外，各无性系的3个对比组中差异表达的基因数量均为200 mmol·L−1VS对照0 mmol·L−1的对比组最多，表明各无性系在200 mmol·L−1高浓度盐溶液处理下受到的胁迫压力最重，对盐胁迫的响应也最强烈。![]()
图 3 差异表达基因韦恩图组注：（1）A-G依次为DJ1、DJ2、DJ3、DJ4、DW1、DY1、DH1无性系的差异表达基因韦恩图，H为200 mmol·L−1处理时强耐盐型与弱耐盐型无性系之间的差异表达基因韦恩图。（2）对照与100、150、200 mmol·L−13个处理浓度分别标记为1、2、3、4附在无性系名称后，如DJ11、DJ12、DJ13、DJ14分别为DJ1无性系的对照与100、150、200 mmol·L−13个处理，其他无性系编号规则相同Figure 3. Group of Venn diagrams of DEGsNotes: (1) A-G represents the venn diagram of DEGs of DJ1, DJ2, DJ3, DJ4, DW1, DY1, and DH1 clones in sequence, H represents the venn diagram of DEGs of between strong salt tolerant and weak salt tolerant clones treated with 200 mmol·L−1. (2) CK and 100, 150, and 200 mmol·L−1 treatments were labeled as 1, 2, 3, and 4, respectively, and attached to the clone names. For example, DJ11, DJ12, DJ13, and DJ14 represent CK and 100, 150, and 200 mmol·L−1 treatments for DJ1, respectively. The numbering rules for other clones are the same2.5.1.2 各无性系差异表达基因富集分析
如图4（A-G）所示，GO功能富集分析显示DJ1株系核心差异表达基因功能主要被显著富集在细胞壁、细胞外区、细胞壁组织、钙调素结合、对过氧化氢的反应、质膜等类别。DJ2株系核心差异表达基因功能主要被显著富集在叶绿体类囊体膜、叶绿体包膜、类囊体、叶绿体类囊体、细胞壁、叶绿体等类别。DJ3株系核心差异表达基因功能主要被显著富集在细胞壁、细胞外区、细胞壁组织、植物型细胞壁、细胞分裂、胞间连丝等类别。DJ4株系核心差异表达基因功能主要被显著富集在细胞壁、DNA结合、细胞外区、核、细胞分裂、DNA结合转录因子活性等类别。DW1株系核心差异表达基因功能主要被显著富集在细胞外泌体、膜、细胞壁、核糖体的结构成分、细胞外区、蛋白质结合等类别。DY1株系核心差异表达基因功能主要被显著富集在DNA结合转录因子活性、质膜、细胞外区、转录与DNA模板、蛋白质丝氨酸/苏氨酸激酶活性、转录调控与DNA模板等类别。DH1株系核心差异表达基因功能主要被显著富集在对缺水的反应、对过氧化氢的反应、细胞外区、序列特异性DNA结合、细胞壁等类别。其中，细胞壁、细胞外区、细胞壁组织、植物型细胞壁、细胞壁生物发生、膜锚定构件、质膜锚定构件、α-L-阿拉伯呋喃糖苷酶活性、DNA结合转录因子活性在7个无性系中均被显著富集，表明上述功能类别基因，尤其是与细胞壁和细胞膜相关的基因对盐胁迫的响应较强。
![]()
图 4 差异表达基因的GO富集分析结果图组注：（1）A-G依次为DJ1、DJ2、DJ3、DJ4、DW1、DY1、DH1无性系的GO富集分析结果图，H为200 mmol·L−1处理时强耐盐型与弱耐盐型无性系之间的GO富集分析结果图Figure 4. Group of GO enrichment analysis results of DEGsNotes: (1) A-G represents the GO enrichment analysis results of DEGs of DJ1, DJ2, DJ3, DJ4, DW1, DY1, and DH1 clones in sequence, H represents the GO enrichment analysis result of DEGs of between strong salt tolerant and weak salt tolerant clones treated with 200 mmol·L−1如图5（A-G）所示，KEGG富集分析显示DJ1株系核心差异表达基因主要被显著富集在淀粉和蔗糖代谢、植物激素信号转导、苯丙烷生物合成、戊糖和葡萄糖醛酸相互转化、植物MAPK信号通路、氧化磷酸化等通路。DJ2株系核心差异表达基因主要被显著富集在光合作用、植物MAPK信号通路、植物激素信号转导、苯丙烷生物合成、戊糖和葡萄糖醛酸相互转化、植物与病原体的相互作用等通路。DJ3株系核心差异表达基因主要被显著富集在戊糖和葡萄糖醛酸相互转化、淀粉和蔗糖代谢、植物激素信号转导、氨基糖与核苷酸糖代谢、半乳糖代谢、苯丙烷生物合成等通路。DJ4株系核心差异表达基因主要被显著富集在植物激素信号转导、DNA复制、植物MAPK信号通路、半乳糖代谢、淀粉和蔗糖代谢、嘧啶代谢等通路。DW1株系核心差异表达基因主要被显著富集在核糖体、氧化磷酸化、吞噬体、内质网中的蛋白质加工、植物激素信号转导、淀粉和蔗糖代谢等通路。DY1株系核心差异表达基因主要被显著富集在植物激素信号转导、苯丙烷生物合成、植物MAPK信号通路、淀粉和蔗糖代谢、戊糖和葡萄糖醛酸相互转化、半乳糖代谢等通路。DH1株系核心差异表达基因主要被显著富集在植物MAPK信号通路、淀粉和蔗糖代谢、半乳糖代谢、植物激素信号转导、氨基糖与核苷酸糖代谢、核糖体等通路。其中，植物激素信号转导、淀粉和蔗糖代谢、苯丙烷生物合成、植物与病原体互作、植物MAPK信号通路在各无性系中富集基因数量较多且在多个无性系中呈显著富集，表明盐胁迫对上述代谢通路产生了较大影响。
![]()
注：（1）A-G依次为DJ1、DJ2、DJ3、DJ4、DW1、DY1、DH1无性系的KEGG富集分析结果图，H为200 mmol·L−1处理时强耐盐型与弱耐盐型无性系之间的KEGG富集分析结果图Notes: (1) A-G represents the KEGG enrichment analysis results of DEGs of DJ1, DJ2, DJ3, DJ4, DW1, DY1, and DH1 clones in sequence, H represents the KEGG enrichment analysis result of DEGs of between strong salt tolerant and weak salt tolerant clones treated with 200 mmol·L−1图5 差异表达基因的KEGG富集分析结果图组Figure 5. Group of KEGG enrichment analysis results of DEGs2.5.2 强耐盐型与弱耐盐型的无性系之间差异表达基因筛选与分析
由2.5.1可知，各无性系在最高浓度200 mmol·L−1盐溶液处理下对盐胁迫响应最强烈、差异表达基因最多，为进一步了解耐盐性好的刺榆无性系响应盐胁迫的机制及筛选耐盐基因，对强耐盐型无性系DJ1、DJ3与弱耐盐型无性系DW1、DY1在最高浓度200 mmol·L−1盐溶液处理时的叶片转录组数据进行比较分析。
如图3（H）所示，DJ14 vs DW14有6 248个差异显著表达基因，3 209个上调、3 039个下调；DJ14 vs DY14有2 954个差异表达基因，1 057个上调、1 897个下调；DJ34 vs DW14有6 930个差异显著表达基因，3 704个上调、3 226个下调；DJ34 vs DY14有
3 913 个差异表达基因，1 687个上调、2 226个下调；4个对比组中共有的差异基因数量为1 406个，其中共同上调基因393个、共同下调基因996个。对4个对比组共有的
1406 个差异基因进行富集分析，如图4（H）GO富集分析所示，差异基因主要富集在细胞外泌体、DNA结合转录因子活性、膜、转录，DNA模板、转录调控-DNA模板、质膜、细胞壁等功能类别。图5（H）KEGG富集分析显示，差异基因主要富集在植物激素信号转导、淀粉和蔗糖代谢、内质网中的蛋白质加工、半乳糖代谢、植物MAPK信号通路、苯丙烷生物合成等代谢通路。植物激素信号转导通路在2.5.2各无性系差异基因分析和本节强耐盐型与弱耐盐型的无性系间差异表达基因分析中均为主要富集通路，且富集基因数量较多，反映其与刺榆响应盐胁迫有较大关系，本研究对其在本节富集的22个差异基因进行分析。如图6代谢通路图所示，生长素信号转导途径中富集6个基因，其中上调节点中2个基因注释为生长素响应蛋白IAA、1个基因注释为生长素响应蛋白SAUR，下调节点中1个基因注释为生长素响应因子ARF、2个基因注释到GH3生长素响应基因家族；细胞分裂素信号转导途径富集2个下调基因，分别注释到细胞分裂素受体CRE1和转录因子ARR-B；赤霉素信号转导途径富集的1个上调基因注释为赤霉素受体GID1，富集的2个下调基因均注释为DELLA 蛋白；脱落酸信号转导途径富集到5个下调基因，均注释为蛋白磷酸酶PP2C；乙烯信号转导途径富集的2个下调基因均注释为乙烯受体ETR；油菜素内酯富集的2个上调基因分别注释到油菜素内酯受体BRI1和转录因子CYCD3，富集的1个下调基因注释为木葡聚糖：木葡糖基转移酶TCH4；茉莉酸信号转导途径富集到1个下调基因，注释为转录因子MYC2。上述途径中相关基因的差异表达表明多个激素参与了DJ1与DJ3对盐胁迫的响应，且这些差异基因的表达可能与两个无性系有较好的耐盐性有关。
![]()
图 6 植物激素信号转导通路中富集的差异表达基因Figure 6. DEGs enriched in plant hormone signal transduction pathways此外，在对4个样本中基因表达量的分析中发现一个在耐盐性好的无性系中大幅表达而在耐盐性差的无性系中表达量较低的基因CNGC15，注释显示其为植物环核苷酸门控通道（cyclic nucleotide gated channel,plant），KEGG分析显示该基因在植物与病原体互作途径中通过调控Ca2 + 进而影响下游防御响应相关基因表达，其大幅表达可能与刺榆DJ1和DJ3株系有较好的耐盐性有关。
2.6 差异表达基因qRT-PCR 验证
如图7所示，TRINITY_DN61310_c1_g2、TRINITY_DN38595_c1_g1、TRINITY_DN63193_c0_g2、TRINITY_DN65673_c0_g24个基因的qRT-PCR结果与转录组测序数据趋势一致、相关性较强，表明本研究的测序结果具有可靠性。
![]()
图 7 差异表达基因的实时荧光定量PCR验证结果注：A-D为4个验证基因的荧光定量PCR与转录组测序数据的趋势对比图，E为荧光定量PCR与转录组测序数据的相关性分析图Figure 7. qRT-PCR verification results of DEGsNotes: A-D shows the trend comparison between qRT-PCR data and transcriptome sequencing data for four verification genes, while E shows the correlation analysis between qRT-PCR data and transcriptome sequencing data for those verification genes3. 讨论
通过转录组测序得到87 712条基因，对各无性系及处理浓度之间进行比较显示，除DW1株系差异基因中下调基因数量略低于上调基因数量外，其他无性系下调差异基因数量均多于上调基因数量，表明在刺榆响应盐胁迫的过程中下调基因发挥了重要作用。盐胁迫会引起植物细胞壁与细胞膜发生变化[21]，本研究中各无性系核心差异基因GO富集分析亦显示与细胞壁和细胞膜相关的功能类别在7个无性系均被富集到。KEGG富集分析显示植物激素信号转导、淀粉和蔗糖代谢、苯丙烷生物合成、植物与病原体互作、植物MAPK信号通路在各无性系中富集基因数量较多且在多个无性系中呈显著富集，与其他学者对金丝楸（C. bungei C.A.Mey）[7]、无瓣海桑（Sonneratia apetala Buchanan-Hamilton）[22]等植物响应盐胁迫的研究结果相似，表明上述通路在7个刺榆无性系响应盐胁迫时发挥了重要作用。
对耐盐性好的DJ1、DJ3株系与耐盐性差的DW1、DY1株系在高盐胁迫下的叶片转录组进行比较，共有的差异基因数量为1 406个，其中共同上调基因393个、共同下调基因996个。对4个对比组共有的1 406个差异基因进行富集分析，GO富集分析显示，差异基因主要富集在细胞外泌体、DNA结合转录因子活性、膜、转录与DNA模板、转录调控与DNA模板、质膜、细胞壁等功能类别；KEGG富集分析显示，差异基因主要富集在植物激素信号转导、淀粉和蔗糖代谢、内质网中的蛋白质加工、半乳糖代谢、植物MAPK信号通路、苯丙烷生物合成等代谢通路，GO富集分析与KEGG富集分析表明200 mmol·L−1高浓度盐溶液胁迫下DJ1与DJ3株系有更强的调节信号转导、蛋白质、碳水化合物、次生代谢物与维护细胞壁与细胞膜结构和功能的能力，上述过程对DJ1和DJ3株系响应盐胁迫、提高抵御盐害的能力具有重要作用。
本研究中多个差异基因被注释到生长素、赤霉素、细胞分裂素、油菜素内酯、脱落酸、乙烯等植物激素信号途径，研究表明这些激素不仅可以调节植物生长发育，在植物响应盐胁迫的过程中也发挥着重要的作用[23-25]。其中脱落酸（ABA）是植物响应盐胁迫非常重要的应激激素，在ABA信号通路中包含多种关键的调节因子，蛋白磷酸酶2C（protein phosphatases 2C，PP2C）作为核心调节因子之一，以负调控的方式发挥作用，其机理为PP2C与ABA正调控因子SnRK2蛋白激酶结合后抑制了其酶活性，阻止了其对下游转录因子的激活，盐胁迫下，PP2C活性会被抑制从而激活ABA信号传导，调控植株耐盐性[26-27]，李梦霜[28]对盐碱胁迫下茶树（Camellia sinensis (L.) O. Ktze.）ABA信号转导系统PP2C基因家族的研究结果显示CsPP2CA6呈现显著的下调表达。本研究中强耐盐型无性系相较于弱耐盐型无性系的差异基因中AHG1、ABI2等5个注释为PP2C的基因均显著下调，与前述研究结果一致。AHG1、ABI2是PP2C 家族中的A组成员，作为核心负调控因子在ABA信号传导通路中发挥着重要作用[29-30]，这2个基因可能与DJ1株系和DJ3株系有较好的耐盐性相关。乙烯是一种气态植物激素，在植物响应盐胁迫过程中同样起着非常重要的调控作用[31-32]。张弢等[33]阐述了在乙烯信号转导途径中ETR1、ERS1、ETR2等乙烯受体可与下游的蛋白激酶CTR1结合阻断乙烯信号传导，即以负反馈形式控制乙烯信号的输出，在对强耐盐型无性系与弱耐盐型无性系的差异基因分析中发现2个乙烯受体基因ERS1和ETR2表达量均显著下调，这2个基因表达量的下调可能会增强乙烯信号传导，进而影响了DJ1株系与DJ3株系对盐胁迫的的耐受能力。此外，注释为CNGC15的基因在强耐盐型无性系中大幅上调表达而在弱耐盐型无性系中表达量较小，其为环核苷酸门控离子通道（CNGCs）家族成员，CNGCs是植物细胞中非选择性阳离子通道和分子开关，参与植物的生长发育，抵御生物和非生物胁迫[34-37]。Saand等[38]以番茄（Solanum lycopersicum L.）沉默株为材料对CNGC15基因的进行了研究，结果显示其与抗旱性有关，与耐盐性相关性较小，但Saand也阐述了该基因在沉默株中有表达，认为需要进一步研究CNGC15基因在植物耐盐胁迫中的作用，在本研究中CNGC15基因在耐盐性好的DJ1株系、DJ3株系大幅上调表达，而在耐盐性较差的DY1株系、DW1株系中表达量较低，可能与刺榆DJ1株系、DJ3株系较好的耐盐性机制有关。CNGC10是CNGCs家族另一重要成员，Guo等[39]、Jin等[40]的研究均显示CNGC10对拟南芥（Arabidopsis thaliana (L.) Heynh.）的耐盐性具有负调控作用，在本研究中，CNGC10在强耐盐型无性系中表达量较弱耐盐型无性系显著下调，根据上述研究推测，差异基因CNGC10亦可能与DJ1株系、DJ3株系较好的耐盐能力相关。
4. 结论
盐胁迫引起刺榆细胞壁、细胞外区、DNA结合转录因子活性等功能类别基因与植物激素信号转导、淀粉和蔗糖代谢、苯丙烷生物合成等代谢通路的响应，强耐盐型无性系较弱耐盐型无性系对上述过程具有更强的调节能力，其中6个显著差异表达基因将作为候选耐盐基因。本研究为揭示刺榆响应盐胁迫的转录组特征与筛选耐盐基因提供了理论参考，为进一步开展刺榆种质资源的研究与应用奠定了基础。
-
![]()
图 1 刺榆不同盐害等级的叶片症状
注：从左到右依次为盐害等级0级、Ⅰ级、Ⅱ级、Ⅲ级
Figure 1. Symptom of different salt damage levels on Hemiptelea davidii leaves
Note: From left to right, the levels of salt damage are 0,Ⅰ,Ⅱand Ⅲ
![]()
图 2 不同浓度盐溶液处理下刺榆各无性系中上部成熟叶片盐害情况
Figure 2. Salt damage performance of the upper part leaves of Hemiptelea davidii clones under different salt stress
![]()
图 3 差异表达基因韦恩图组
注：（1）A-G依次为DJ1、DJ2、DJ3、DJ4、DW1、DY1、DH1无性系的差异表达基因韦恩图，H为200 mmol·L−1处理时强耐盐型与弱耐盐型无性系之间的差异表达基因韦恩图。（2）对照与100、150、200 mmol·L−13个处理浓度分别标记为1、2、3、4附在无性系名称后，如DJ11、DJ12、DJ13、DJ14分别为DJ1无性系的对照与100、150、200 mmol·L−13个处理，其他无性系编号规则相同
Figure 3. Group of Venn diagrams of DEGs
Notes: (1) A-G represents the venn diagram of DEGs of DJ1, DJ2, DJ3, DJ4, DW1, DY1, and DH1 clones in sequence, H represents the venn diagram of DEGs of between strong salt tolerant and weak salt tolerant clones treated with 200 mmol·L−1. (2) CK and 100, 150, and 200 mmol·L−1 treatments were labeled as 1, 2, 3, and 4, respectively, and attached to the clone names. For example, DJ11, DJ12, DJ13, and DJ14 represent CK and 100, 150, and 200 mmol·L−1 treatments for DJ1, respectively. The numbering rules for other clones are the same
![]()
图 4 差异表达基因的GO富集分析结果图组
注：（1）A-G依次为DJ1、DJ2、DJ3、DJ4、DW1、DY1、DH1无性系的GO富集分析结果图，H为200 mmol·L−1处理时强耐盐型与弱耐盐型无性系之间的GO富集分析结果图
Figure 4. Group of GO enrichment analysis results of DEGs
Notes: (1) A-G represents the GO enrichment analysis results of DEGs of DJ1, DJ2, DJ3, DJ4, DW1, DY1, and DH1 clones in sequence, H represents the GO enrichment analysis result of DEGs of between strong salt tolerant and weak salt tolerant clones treated with 200 mmol·L−1
![]()
图 6 植物激素信号转导通路中富集的差异表达基因
Figure 6. DEGs enriched in plant hormone signal transduction pathways
![]()
图 7 差异表达基因的实时荧光定量PCR验证结果
注：A-D为4个验证基因的荧光定量PCR与转录组测序数据的趋势对比图，E为荧光定量PCR与转录组测序数据的相关性分析图
Figure 7. qRT-PCR verification results of DEGs
Notes: A-D shows the trend comparison between qRT-PCR data and transcriptome sequencing data for four verification genes, while E shows the correlation analysis between qRT-PCR data and transcriptome sequencing data for those verification genes
表 1 实时荧光定量PCR引物序列
Table 1 Primer sequences for qRT-PCR
基因名称与方向
Unigene name and direction引物序列
Primer sequence（5′-3′）TRINITY_DN61310_c1_g2F CGGAGACTAAGGCCAACT TRINITY_DN61310_c1_g2R GCTTTCCTGGCAACTGTACTG TRINITY_DN38595_c1_g1F CGCTCAGATGGGACTTTCAAC TRINITY_DN38595_c1_g1R GAGGAGACCAGTAGATCGATCC TRINITY_DN63193_c0_g2F TGTGAGAGTCCAGGCTGC TRINITY_DN63193_c0_g2R CCTTCCTCGATGTTCACCTC TRINITY_DN65673_c0_g2F TCTGTCCTTCGTGCTCAATTACC TRINITY_DN65673_c0_g2R CACCTACCACTCCAGAGCACTT GAPDH-F TCAACAATGCCAAACCTG GAPDH-R GTGTCAACGAGCACGAAT 
下载: 导出CSV
表 2 注释结果统计
Table 2 Annotation result statistics
数据库
Database基因数量
Number of Unigenes比率/%
RatioGO 49 522 56.46 KEGG 31 171 35.54 Pfam 45 055 51.37 Swissprot 42 904 48.91 eggNOG 55 769 63.58 NR 42 350 48.28 Total Unigenes 87 712 100.00
下载: 导出CSV
-
[1] 杨劲松. 中国盐渍土研究的发展历程与展望[J]. 土壤学报,2008,45(5):837-845. DOI: 10.3321/j.issn:0564-3929.2008.05.010 [2] 范晓梅, 刘高焕, 唐志鹏, 等. 黄河三角洲土壤盐渍化影响因素分析[J]. 水土保持学报,2010,24(1):139-144. [3] 孙小茹. 近四十年河套平原土壤盐渍化时空变异及其影响因素分析[D]. 武汉: 华中农业大学, 2022. [4] 云雪雪, 陈雨生. 国际盐碱地开发动态及其对我国的启示[J]. 国土与自然资源研究,2020(1):84-87. [5] 刘小京, 郭 凯, 封晓辉, 等. 农业高效利用盐碱地资源探讨[J]. 中国生态农业学报(中英文),2023,31(3):345-353. DOI: 10.12357/cjea.20220967 [6] 刘 晨, 徐浩博, 斯钰阳, 等. 基于转录组学的植物响应盐胁迫调控机制研究进展[J]. 浙江农业学报,2022,34(4):870-878. DOI: 10.3969/j.issn.1004-1524.2022.04.24 [7] 郜新强, 王小艳, 焦 伟, 等. 金丝楸幼苗响应盐碱胁迫的生理和转录组分析[J]. 林业科学研究,2023,36(1):166-178. DOI: 10.12403/j.1001-1498.20220180 [8] 陈亚辉, 张 师, 杨庆山, 等. 多枝柽柳叶片响应NaCl胁迫的转录组分析[J]. 江苏农业学报,2022,38(5):1188-1202. [9] 王新亮, 彭 玲, 王 健, 等. 盐碱胁迫下平邑甜茶的转录组分析[J]. 江苏农业科学,2022,50(22):29-37. [10] 李法曾, 李文清, 樊守金. 山东木本植物志(上卷)[M]. 北京: 科学出版社, 2016: 195-196. [11] 杨 颖. 刺榆种子萌发及幼苗耐盐性研究[D]. 呼和浩特: 内蒙古农业大学, 2019. [12] 刘华波. 刺榆细胞器基因组分析及耐旱分子机制[D]. 北京: 北京林业大学, 2020. [13] 刘炳响, 王志刚, 梁海永, 等. 盐胁迫对不同生境白榆生理特性与耐盐性的影响[J]. 应用生态学报,2012,23(6):1481-1489. [14] GRABHERR M G, HAAS B J, YASSOUR M, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome[J]. Nature Biotechnology,2011,29(7): 644-652. DOI: 10.1038/nbt.1883
[15] BUCHFINK B, XIE C, HUSON D H. Fast and sensitive protein alignment using DIAMOND[J]. Nature Methods,2015,12(1): 59-60. DOI: 10.1038/nmeth.3176
[16] PATRO R, DUGGAL G, LOVE M I, et al. Salmon provides fast and bias-aware quantification of transcript expression[J]. Nature Methods,2017,14(4): 417-419. DOI: 10.1038/nmeth.4197
[17] MORTAZAVI A, WILLIAMS B A, MCCUE K, et al. Mapping and quantifying mammalian transcriptomes by RNA-Seq[J]. Nature Methods,2008,5(7): 621-628. DOI: 10.1038/nmeth.1226
[18] ROBINSON M D, MCCARTHY D J, SMYTH G K. EdgeR: a bioconductor package for differential expression analysis of digital gene expression data[J]. Bioinformatics,2010,26(1): 139-140. DOI: 10.1093/bioinformatics/btp616
[19] ZHONG H Y, CHEN J W, LI C Q, et al. Selection of reliable reference genes for expression studies by reverse transcription quantitative real-time PCR in litchi under different experimental conditions[J]. Plant Cell Reports,2011,30(4): 641-653. DOI: 10.1007/s00299-010-0992-8
[20] LIVAK K J, SCHMITTGEN T D. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) method[J]. Methods,2001,25(4): 402-408. DOI: 10.1006/meth.2001.1262
[21] ZHANG H M, ZHU J H, GONG Z Z, et al. Abiotic stress responses in plants[J]. Nature Reviews Genetics,2021,23(2): 104-119.
[22] 梁锐涛, 韩维栋, 杨少瑕, 等. 无瓣海桑根响应盐胁迫的转录组分析[J]. 林业科学研究,2023,36(1):68-78. DOI: 10.12403/j.1001-1498.20230142 [23] WAADT R, SELLER C A, HSU P K, et al. Plant hormone regulation of abiotic stress responses[J]. Nature Reviews Molecular Cell Biology,2022,23(10): 680-694. DOI: 10.1038/s41580-022-00479-6
[24] 韦建辰. NaCl胁迫下越橘组培苗转录组与代谢组分析[D]. 长春: 吉林农业大学, 2021. [25] YU Z P, DUAN X B, LUO L, et al. How Plant hormones mediate salt stress responses[J]. Trends in Plant Science,2020,25(11): 1117-1130. DOI: 10.1016/j.tplants.2020.06.008
[26] ZHU J K. Abiotic stress signaling and responses in plants[J]. Cell,2016,167(2): 313-324. DOI: 10.1016/j.cell.2016.08.029
[27] 陈素梅, 王新慧, 李 菲, 等. ABA调控植物盐胁迫应答机制研究进展[J]. 南京农业大学学报,2022,45(5):856-863. DOI: 10.7685/jnau.202103010 [28] 李梦霜. 茶树ABA信号转导家族基因鉴定及功能分析[D]. 合肥: 安徽农业大学, 2020. [29] RODRIGUEZ P L, BENNING G, GRILL E. ABI2, a second protein phosphatase 2C involved in abscisic acid signal transduction in Arabidopsis[J]. FEBS Letters,1998,421(3): 185-190. DOI: 10.1016/S0014-5793(97)01558-5
[30] NISHIMURA N, TSUCHIYA W, MORESCO J J, et al. Control of seed dormancy and germination by DOG1-AHG1 PP2C phosphatase complex via binding to heme[J]. Nature Communications,2018,9(1): 2132. DOI: 10.1038/s41467-018-04437-9
[31] 赵 赫, 陈受宜, 张劲松. 乙烯信号转导与植物非生物胁迫反应调控研究进展[J]. 生物技术通报,2016,32(10):1-10. [32] RIYAZUDDIN R, VERMA R, SINGH K, et al. Ethylene: a master regulator of salinity stress tolerance in plants[J]. Biomolecules,2020,10(6): 959. DOI: 10.3390/biom10060959
[33] 张 弢, 董春海. 乙烯信号转导及其在植物逆境响应中的作用[J]. 生物技术通报,2016,32(10):11-17. [34] LU Z Y, YIN G, CHAI M, et al. Systematic analysis of CNGCs in cotton and the positive role of GhCNGC32 and GhCNGC35 in salt tolerance[J]. BMC Genomics,2022,23(1): 560. DOI: 10.1186/s12864-022-08800-5
[35] 李可新, 赵 芊, 丁 一, 等. 植物内环核苷酸门控离子通道的研究进展[J]. 植物生理学报,2021,57(12):2235-2246. [36] KUGLER A, KÖHLER B, PALME K, et al. Salt-dependent regulation of a CNG channel subfamily in Arabidopsis[J]. BMC Plant Biology,2009,9(1): 140. DOI: 10.1186/1471-2229-9-140
[37] GOBERT A, PARK G, AMTMANN A, et al. Arabidopsis thaliana cyclic nucleotide gated channel 3 forms a non-selective ion transporter involved in germination and cation transport[J]. Journal of Experimental Botany, 2006, 57(4): 791-800.
[38] SAAND M A, XU Y P, MUNYAMPUNDU J P, et al. Phylogeny and evolution of plant cyclic nucleotide-gated ion channel (CNGC) gene family and functional analyses of tomato CNGCs[J]. DNA Research,2015,22(6): 471-483. DOI: 10.1093/dnares/dsv029
[39] GUO K M, BABOURINA O, CHRISTOPHER D A, et al. The cyclic nucleotide-gated channel, AtCNGC10, influences salt tolerance in Arabidopsis[J]. Physiol Plant,2008,134(3): 499-507. DOI: 10.1111/j.1399-3054.2008.01157.x
[40] JIN Y K, JING W, ZHANG Q, et al. Cyclic nucleotide gated channel 10 negatively regulates salt tolerance by mediating Na + transport in Arabidopsis[J]. Journal of Plant Research, 2015, 128(1): 211-220.
计量
- 文章访问数: 17
- HTML全文浏览量: 4
- PDF下载量: 14
下载:
