Transcriptome Analysis and Gene Discovery of Abscisic Acid Signaling Pathway in Kandelia obovata under Low Temperature Stress

以浙江省亚热带作物研究所选育的耐寒红树植物‘龙港’秋茄1年生容器苗为植物材料，挑选生长良好，长势一致，无病害植株作为供试材料（图1）。试验于2020年12月在浙江省农业科学院蔬菜研究所内进行，将供试材料置于人工气候室（型号：SAFE-DG-ZJNKY-6）中，光照时间为11 h·d⁻¹，15 ℃ 处理12 h为对照组（CK），−5 ℃处理12 h为低温组（LT），各处理设置4次生物学重复。各处理结束后采集完全展开的功能叶片用锡箔纸包裹并编号（CK 组编号为：CK-1、CK-2、CK-3 和 CK-4；LT组编号为：LT-1、LT-2、LT-3 和 LT-4），将包裹好的叶片迅速置于液氮中速冻，最后转移至−80 ℃超低温冰箱中保存备用。

Figure 1.  The Kandelia obovata material used in the experiment

将8份样品（CK-1、CK-2、CK-3、CK-4、LT-1、LT-2、LT-3和LT-4）委托上海美吉生物医药科技有限公司进行样品RNA的提取、质控、建库及转录组测序，用Illumina Novaseq 6000平台进行测序。

将测序所获得的clean reads与秋茄参考基因组（https://bigd.big.ac.cn/gwh/Assembly/990/show）进行比对并统计比对率，比对分析软件为HISAT。

FPKM（fragments per kilobase million）即每100万条序列中，每个基因以1 000个碱基为单位，比对上的reads。以FPKM值计算基因表达量，计算软件为RSEM。

基于表达量定量结果进行组间差异基因分析，获得2组间发生差异表达的基因，差异分析软件为DESeq2，筛选阈值为|log₂FC| ≥1和padjust <0.05。同时利用美吉生物云平台（https://cloud.majorbio.com）对筛选得到的DEGs进行GO功能注释和KEGG富集分析。

为验证测序结果准确性，本研究选择了6个KoHB-others和5个ABA信号途径DEGs进行实时荧光定量PCR（Quantitative Real-time PCR）验证。利用primer 5.0软件对所挑选基因进行特异性引物设计，内参基因为KoACT2^[19]，引物及内参序列见表1，以2^−△△CT法计算DEGs的相对表达水平。

对秋茄叶片样品进行RNA提取后，经转录组测序获得样品原始数据。通过排除低质量数据和序列拼接，获得了48 669 538 ～ 61 694 828的Clean reads（即质控后测序数据的总条目数），错误率仅为0.02％ 左右，Q20均高于98.14%，Q30均高于94.41%，GC含量在45.11% ～ 45.60%（表2）。同时，样品主成分分析结果（图2）表明：CK组与LT组各处理分别聚为一类，说明样品间的生物学重复性较好。因此，本研究的转录组测序结果较理想，后续的分析较可信。

Figure 2.  PCA analysis of samples

本研究参考基因为Kandelia obovata（https://bigd.big.ac.cn/gwh/Assembly/990/show），所参考基因组版本为ko01。通过将质控后的clean data与参考基因组进行对比发现，实验所产生的测序序列的比对率均高于96%（表3），包括总比对（能定位到基因组上的clean reads数目）、多方比对（在参考序列上有多个比对位置的clean reads数目）和唯一比对（在参考序列上有唯一比对位置的Clean reads数目）。

本研究转录组测序共鉴定到148个转录因子，分属于25个转录因子家族，其中，ERF、NAC、WRKY、bHLH、MYB、bZIP、HB-other和MYB-related等转录因子家族所包含的基因数目较多，分别占11.49%、9.46%、8.11%、8.11%、6.76%、6.08%、4.05%和4.05%（表4）。其中，HB-other家族转录因子包含6个DEGs，KoHB-other1（Maker00002036）、KoHB-other2 （Maker00003789）和KoHB-other3 （Maker00008096）下调表达，KoHB-other4（Maker00014773）、KoHB-other5（Maker00016058） 和KoHB-other6 （Maker00017922）上调表达（图3）。

Figure 3.  Expression of HB-Other transcription factors DEGs

以P<0.05、上下调差异倍数FC>2或者FC<0.5为筛选条件，进一步对DEGs数量进行统计学分析。结果表明，LT组和CK组的DEGs共有1 330个，其中，698个DEGs上调表达（52.48%），632个DEGs下调表达（47.52%）（图4）。

Figure 4.  DEGs volcano

从基因数量上看，DEGs主要富集在植物激素信号转导、苯丙素生物合成、植物与病原体互作、淀粉与蔗糖代谢、MAPK信号通路-植物、半乳糖代谢、光合作用-天线蛋白和α-亚麻酸代谢等通路中（表5）。以P值作为参考，DEGs则显著富集在苯丙素的生物合成、植物激素信号转导、光合作用-天线蛋白、半乳糖代谢、"角质、亚伯碱和蜡的生物合成"和α-亚麻酸代谢等通路中（表6）。综合基因数量和P值看，低温胁迫下秋茄叶片DEGs显著富集在植物激素信号转导、苯丙素生物合成、半乳糖代谢、光合作用-天线蛋白和α-亚麻酸代谢这5条KEGG通路上，其中，植物激素信号转导通路中富集到的DEGs数量最多，且差异显著性最高。

ABA信号途径是植物激素信号转导通路中一条重要激素途径，对该途径中的DEGs进行深入挖掘发现，KoPYL1（Maker00013579）、KoABF1（Maker00005324）和KoABF2（Maker00011633）基因上调表达，KoPP2C1（Maker00008505）和KoABF3（Maker00004855）基因下调表达（图5）。

Figure 5.  Expression of DEGs in ABA signaling pathway

对重点关注的6个 KoHB-others和5个脱落酸信号途径DEGs进行了qRT-PCR验证，结果表明这11个基因的变化趋势与转录组测序结果一致（图6），表明转录组数据准确可靠。

Figure 6.  qRT-PCR validation of DEGs

转录组测序（RNA-seq）技术对基因挖掘具有较高的准确性和灵敏性^[20]，被广泛用于拟南芥^[10]、油菜（Brassica rapa L.）^[21]、新疆野苹果（Malus sieversii (Ledeb.) Roem.）^[22]等多种植物耐寒基因的研究中，使植物抗寒分子机制被不断地深入发掘。本研究中，各样品Q20碱基百分比均大于98%，且Q30碱基百分比均大于94%，GC含量介于45% ～ 46%，测序序列与参考基因组的比对率均高于96%。此外，qRT-PCR结果显示，所挑选DEGs的变化趋势与转录组结果基本一致。因此，本研究转录组测序数据结果良好，可满足后续分析需求。

转录因子通过与胁迫应答基因启动子中的特定顺式作用元件结合，成为植物抗非生物胁迫的主要中介因子，从而使植物能够抵御不利的环境条件^[23]。ERF、NAC、WRKY、bHLH、MYB、bZIP和MYB-related 等家族转录因子参与了多种植物的非生物胁迫响应过程^[24]，在对拟南芥^[25]、沙冬青（Ammopiptanthus Mongolicus (Maxim. ex Kom.) Cheng f.）^[26]、烟草（Nicotiana tabacum L.）^[27]、山荆子（Malus baccata (L.) Borkh.）^[28]、玉米（Zea mays L.）^[29]和水稻^[30]等的研究发现，相关家族转录因子均能一定程度上提升植株的耐寒性，有些甚至是适应冷胁迫所必须的。值得关注的是，HB-other 转录因子响应非生物胁迫的研究较少，仅在脐橙（Citrus sinensis Osbeck）^[31]和地梢瓜（Cynanchum thesioides (Freyn) K. Schum）^[32]的干旱胁迫响应中被鉴定到，而在低温胁迫的相关研究中尚未见报道。在本研究中，共鉴定到148个转录因子，分属于25个转录因子家族，其中，ERF、NAC、WRKY、bHLH、MYB、bZIP、HB-other和MYB-related等转录因子家族所包含的基因数目较多。这些转录因子中除HB-other在抗寒相关的研究中未被鉴定到外，其余均被鉴定到并被证明是与抗寒相关的重要转录因子家族，特别是bHLH、MYB-related和WRKY在前期秋茄抗寒相关的研究中被认定为低温应激反应中有价值的候选调节因子^[11-12]在本次秋茄低温胁迫中首次鉴定出HB-other转录因子，并且KoHB-others的表达情况与qRT-PCR验证结果一致，推测其可能对秋茄响应低温胁迫起重要调控作用，在后期的研究中还有待被证实。

京都基因和基因组百科全书（Kyoto Encyclopedia of Genes and Genomes, KEGG）是用于分析代谢途径、生物系统和基因功能的信号通路数据库^[33]。许多 KEGG 通路与植物响应低温胁迫有关，如橄榄（Canarium album (Lour.) Rauesch.）在低温胁迫下的DEGs富集到了α-亚麻酸代谢、类胡萝卜素生物合成、光合-天线蛋白和昼夜节律相关等通路上^[34]，文冠果（Xanthoceras sorbifolium Bunge）在低温胁迫下DEGs显著富集在淀粉和蔗糖代谢、半乳糖代谢和氨基酸代谢等通路上^[35]。本研究中，秋茄叶片DEGs显著富集在植物激素信号转导、苯丙素生物合成、半乳糖代谢、光合作用-天线蛋白和α-亚麻酸代谢等多条KEGG通路上，且这些通路已经被大量研究证实与抗寒相关，其中，植物激素信号转导、α-亚麻酸代谢和苯丙酸生物合成等通路在前期秋茄抗寒相关的研究中已经得到证明^[12]。因此，植物激素信号转导、苯丙素生物合成、半乳糖代谢、光合作用-天线蛋白和α-亚麻酸代谢等是秋茄响应低温胁迫的重要KEGG通路。

植物激素是一种小的内源性信号分子，如ABA、茉莉酸和生长素等，在响应非生物胁迫过程中参与信号转导途径^[36]。ABA是植物适应低温胁迫的重要信号分子^[37]，并且其响应低温胁迫的“PYR/PYL/RCAR-|PP2Cs|-SnRK2-ABF/AREB-靶基因”途径已经研究的较为清晰^[38]。ABA信号途径基因已在许多植物中被鉴定和研究。Zhang 等^[39]研究发现，拟南芥过表达的RCAR12或RCAR13可以通过诱导低温响应基因CBFs的表达来耐受低温胁迫。Ren等^[40]在葡萄（Vitis vinifera L.）愈伤组织和拟南芥中过表达VaPYL4基因可增强转基因株系的耐寒性。大量研究表明，PP2Cs负向调控ABA信号途径^[41]；但Hu等^[42]在烟草中过表达玉米根部的ZmPP2C2基因增强了烟草对低温胁迫的耐受性，表明ZmPP2C2基因可能是玉米抗寒性的正向调控基因，这与大部分抗寒相关的研究结论相反，这可能是ZmPP2Cs的调控方式具有特殊性，具体原因还有待进一步的研究。高山离子芥（Chorispora bungeana Fisch. et Mey.）的CbABF1在烟草中表达提高了植株对冷胁迫的耐受性，且转基因植株的存活率高于野生型^[43]。在本研究中，KoPYL1（Maker00013579）、KoABF1（Maker00005324）和KoABF2（Maker00011633）上调表达，KoPP2C1（Maker00008505）和KoABF3（Maker00004855）下调表达，这和其他植物的相关研究结果一致，且这些基因的表达情况与qRT-PCR验证结果一致。因此，KoPYL1、KoPP2C1、KoABF1、KoABF2和KoABF3基因可作为后期研究秋茄响应低温胁迫的重要候选基因。

ERF、NAC、WRKY、bHLH、MYB、bZIP、HB-other和MYB-related等家族转录因子对秋茄响应低温胁迫起重要调控作用；植物激素信号转导、苯丙素生物合成、半乳糖代谢、光合作用-天线蛋白和α-亚麻酸代谢等是秋茄响应低温胁迫的重要KEGG通路；ABA信号途径中的KoPYL1、KoPP2C1、KoABF1、KoABF2和KoABF3等基因可作为后期研究秋茄响应低温胁迫的重要候选基因。