Identification of the durian LEA gene family and expression analysis under low temperature stress

榴莲的基因组序列和注释文件来自figshare在线网站（https://figshare.com/ articles/dataset/Durian_genome_annotation/25237591）^[12]，拟南芥LEAs家族的基因组序列和注释文件来自tair数据库（https://www.arabidopsis.org/）^[6]，利用TBtools软件^[13]（版本：2.326）通过BLAST方法对榴莲LEA基因家族成员进行鉴定。同时，从Pfam数据库（http://pfam.xfam.org/）中下载了LEA1（PF03760）、LEA2（PF03168）、LEA3（PF03242）、LEA4（PF02987）、LEA5（PF00477）、LEA6（PF10714）、DHN（PF00257）、 和SMP（PF04927）亚家族的隐马尔可夫模型（Hidden Markov Model，HMM）。通过TBtools的simple HMM search程序得到榴莲候选LEAs基因。通过NCBI conserved Domains（http://www. ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi）数据库、Pfam数据库（http://pfam-legacy.xfam.org/），对BLASTP和HMM两种方法得到的榴莲LEAs候选基因进行结构域的确定。根据LEA家族成员在染色体上的位置，将其依次编号为DzLEA1～DzLEA46用于后续分析。

通过TBtools软件Protein Parameter Calc程序^[13]，分析DzLEAs蛋白的氨基酸长度（Amino acid length）、理论等电点（Theoretical Isoelectric Point）、分子量（Molecular Weight）、不稳定系数（Instability Index）、亲水性平均系数（Grand Average of Hydropathicity，GRAVY）等理化性质，使用WoLF PSORT网站（https://wolfpsort.hgc.jp/）预测DzLEAs蛋白的亚细胞定位^[14]。

通过MEGA12软件中最大似然法（Maximum Likelihood, ML）构建系统进化树，Bootstrap重复值设置为1 000，使用iTOL（https://itol.embl.de/）在线工具对生成的系统进化树进行美化^[15]。

在NCBI保守结构域数据库（https://www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb.cgi）分析榴莲LEA基因家族的保守结构域；使用MEME套件（https://meme-suite.org/meme/）分析榴莲LEA基因家族的保守基序，搜索motif数量设置为20个；根据榴莲的基因组注释，分析DzLEAs的基因结构。最后，通过TBtools的Gene Structure view（Advanced）程序对DzLEAs的进化树、保守基序、保守结构域和基因结构进行可视化^[13]。

通过TBtools软件获取DzLEAs基因起始密码子上游2 000 bp的序列，在PlantCARE网站（https://bioinformatics.psb.ugent.be/webtools/plantcare/html/）对获取的序列进行分析，最后通过TBtools软件对分析结果进行可视化^[13]。

基于榴莲和拟南芥的基因组文件和注释文件^[12]，通过TBtools软件^[13]获取其染色体长度、基因密度和位置信息，通过MCScanX对榴莲与拟南芥LEA基因家族进行物种内和物种间共线性分析，采用默认参数；最后通过Advanced Circos程序对共线性分析结果进行可视化。

在eggNOG-mapper网站（http://eggnog-mapper.embl.de/）获取榴莲的GO和KEGG注释文件^[16]，之后通过TBtools软件中GO Enrichment和KEGG Enrichment Analysis程序中对DzLEAs基因进行富集分析，最后通过Enrichment Bar Plot程序对分析结果进行可视化^[13]。

以‘猫山王’和‘金枕’2个榴莲品种为研究对象，榴莲苗高为1 m，苗龄为12个月，其中，处理组置于恒温培养箱中进行培养，光照设置为12 h光照和12 h黑暗循环，温度设置为10℃；对照组置于室内空调屋进行培养，空气相对湿度约为45%，温度设置为26℃。取26℃常温、10℃低温处理1 d、10℃低温处理5 d叶片样品进行转录组数据分析，每个处理重复3次。

分别取‘猫山王’和‘金枕’在26℃常温、10℃处理1 d、10℃处理3 d、10℃处理5 d、10℃处理7 d的叶片样品进行qRT-PCR实验，生物学重复3次取样，混样分3次技术重复进行后续qRT-PCR实验。采用十六烷基三甲基溴化铵法提取叶片样品RNA，通过超微紫外分光光度计（ND-100C）测定提取的RNA浓度。根据所测样品RNA浓度统一其含量，使用含gDNA清除剂的SPARKscript Ⅱ RT Plus Kit试剂盒（SparkJade）合成第1链互补DNA（cDNA）。qRT-PCR的引物（表1）通过NCBI中的Primer-BLAST（https://www.ncbi.nlm.nih.gov/tools/primer-blast/）工具设计^[17]。通过Applied Biosystems 7500实时定量PCR系统（Life Technologies，USA）进行qRT-PCR，反应体系包括10 μL 2×SYBR qPCR Mix、1 μL cDNA、0.4 μL Forward Primer、0.4 μL Reverse Primer、0.4 μL ROX Reference Dye Ⅱ、7.8 μL RNase Free H₂O，总体积为20 μL。采用相对定量法（ΔΔCt法）计算样品的相对表达量。

基于BLAST和HMM两种方法，从榴莲基因组中共鉴定出46个LEA基因家族成员，根据其在染色体上的位置分别命名为DzLEA1～DzLEA46（图1）。

Figure 1.  Chromosomal localization of the DzLEAs gene

由表2可知，46个榴莲LEA基因家族成员编码的蛋白氨基酸序列长度在65～467个之间，蛋白质的理论分子质量范围在6.965 8～52.080 58 kDa之间。理论等电点在4.33～9.79之间，其中，27个（59%）DzLEA蛋白的等电点小于7，为酸性蛋白；19个（41%）DzLEA蛋白的等电点大于7，为碱性蛋白。46个DzLEA蛋白中，稳定蛋白和不稳定蛋白均占50%（23个）。脂肪族系数在16.67～107.81之间，且大多数DzLEA蛋白为亲水性蛋白（44个），仅DzLEA20和DzLEA30为疏水性蛋白。

DzLEA蛋白的亚细胞定位预测表明：LEA1亚家族的1个成员（DzLEA1）、LEA5亚家族的3个成员（DzLEA13、DzLEA15、DzLEA29）和LEA6亚家族的2个成员（DzLEA6、DzLEA11）均位于细胞核内；LEA2亚家族的1个成员（DzLEA20）位于质膜内，其余3个成员（DzLEA4、DzLEA19、DzLEA40）位于细胞质内；LEA3亚家族的1个成员（DzLEA9）位于线粒体内，其余5个成员（DzLEA21、DzLEA26、DzLEA32、DzLEA39、DzLEA41）位于叶绿体内；LEA4亚家族的1个成员（DzLEA43）位于细胞外，5个成员（DzLEA10、DzLEA12、DzLEA35、DzLEA36、DzLEA37）位于线粒体内，其余均位于细胞核内；DHN亚家族的1个成员位于细胞质内，其余DHN成员位于细胞核内；SMP亚家族的成员在细胞核、细胞质、线粒体和叶绿体中均有存在。

将榴莲与拟南芥LEA家族成员共同构建系统进化树（图2），参考拟南芥中LEA蛋白亚家族分类结果^[6]，发现DzLEAs共分为8个亚家族：LEA1、LEA2、LEA3、LEA4、LEA5、LEA6、DHN和SMP。其中，LEA4亚家族成员最多为15个，LEA2、LEA3、LEA5、LEA6、DHN和SMP亚家族成员分别有4、6、4、2、5和9个，LEA1亚家族成员最少为1个。榴莲与拟南芥的蛋白分类几乎一致，说明榴莲LEA基因在长期进化的过程中，这些蛋白质没有遗失。

Figure 2.  Phylogenetic tree and subfamily classification of LEA proteins in durian and Arabidopsis

本研究对46个DzLEAs基因的保守基序进行分析，共鉴定出20个motif（motif1～motif20）（图3-b）。结果表明，LEA1亚家族成员不包含任何motif；LEA2亚家族成员除DzLEA20外，其余成员均包含motif10；LEA3亚家族成员主要包含motif2和motif7；LEA4亚家族中，除了4位成员（DzLEA2、DzLEA12、DzLEA33、DzLEA43），其余成员均包含motif4，且绝大多数成员包含motif1，少部分成员包含motif20；LEA5亚家族成员均包含motif11，且除了DzLEA3外，其余成员均包含motif18；LEA6亚家族成员均包含motif15；DHN亚家族成员主要包含motif5、motif9、motif13和motif17；SMP亚家族成员主要包含motif3、motif5、motif6、motif8、motif9和motif12。DzLEAs基因的保守结构域分析显示（图3-c），不同亚家族成员均含有其特定保守结构域。本研究发现，不同亚家族成员的motif类型和分布存在较大差异，而同一亚家族成员之间motif结构较为类似，表明相同亚家族成员之间功能具有相似性。

Figure 3.  The phylogenetic tree (a), conserved motifs (b), protein domains (c), and gene structures (d)of DzLEA genes.

榴莲LEA基因家族成员基因结构分析显示（图3-d），LEA1亚家族成员包含1个内含子；LEA2亚家族成员含有0～2个内含子，其中DzLEA4基因不含有内含子，DzLEA19和DzLEA40均含有1个内含子，DzLEA20含有2个内含子；LEA3亚家族成员均含有1个内含子；LEA4亚家族成员除3个成员外（DzLEA2、DzLEA10、DzLEA33）均含有2个内含子；LEA5亚家族成员中，DzLEA29含有5个内含子，DzLEA3含有4个内含子，其余2个成员（DzLEA13、DzLEA15）含有1个内含子；DHN亚家族成员均含有1个内含子；SMP亚家族成员含有0～4个内含子，其中DzLEA45不含有内含子，DzLEA38含有4个内含子，其余成员均含有2个内含子。在同一亚家族成员中，具有相似的基因结构，并且进化关系密切的成员之间表现出更为相似的外显子-内含子结构、外显子数量和基因长度。

对46个榴莲LEA基因进行了顺式作用元件分析（图4），结果表明，发现了5个与光响应相关的元件，分别为3-AF1 binding site（5个成员）、ACE（9个成员）、G-box（37个成员）、GT1-motif（28个成员）和Sp1（8个成员）；3个与赤霉素响应相关的元件，分别为GARE-motif（10个成员）、P-box（13个成员）和TATC-box（7个成员）；茉莉酸甲酯响应元件CGTCA-motif（36个成员）和TGACG-motif（37个成员）；生长素响应元件AuxRR-core（4个成员）和TGA-element（11个成员）；脱落酸响应元件ABRE（38个成员）；水杨酸响应元件TCA-element（25个成员）；与低温响应的元件LTR（16个成员）；参与防御和应激反应的元件TC-rich repeats（15个成员）；参与干旱响应的元件MBS（22个成员）。除此之外，还发现了与厌氧诱导相关的元件（ARE）、与分生组织表达相关的元件（CAT-box）、MYB结合位点（MRE、CCAAT-box）、昼夜节律控制相关的元件（circadian）、胚乳表达相关的元件（GCN4_motif）、玉米醇溶蛋白代谢调控相关的元件（O2-site）、种子特异性调控元件（RY-element）和伤口响应元件（WUN-motif）。这些结果表明，DzLEAs基因可能参与多种植物生长发育过程，并且在逆境胁迫中发挥重要作用。

Figure 4.  Distribution of cis-regulatory elements in the promoters of DzLEA genes.

对DzLEA家族成员在染色体上的位置进行分析（图1），46个DzLEA家族成员分布在18条染色体上，且呈现不均匀分布。其中，17号和24号染色体均含有7个DzLEAs基因；1号和21号染色体上分别含有5个和4个DzLEAs基因；10号和18号染色体上均含有3个DzLEAs基因；2号、6号、14号、25号和27号染色体均含有2个DzLEAs基因；4号、9号、12号、16号、19号、26号和28号染色体各含有1个DzLEAs基因。并且，通过染色体定位图显示，DzLEA22和DzLEA23、DzLEA27和DzLEA28、DzLEA36和DzLEA37、DzLEA44和DzLEA45，这四组基因距离相近，推测以上成员存在遗传连锁和多种复制事件的可能性。

同时对DzLEAs基因进行了物种内共线性分析（图5），共发现28组共线性基因对，结果显示榴莲LEA基因家族至少发生了19次片段复制事件，但并未发现串联复制事件，这导致榴莲LEA家族数量上的扩增。同时，发现部分DzLEAs基因与多个基因存在共线性关系，如DzLEA5和DzLEA7、DzLEA17、DzLEA18；DzLEA21和DzLEA26、DzLEA32、DzLEA39；DzLEA22和DzLEA27、 DzLEA30、DzLEA38，这些基因分别聚类于DHN、LEA3和SMP亚家族，说明其在功能与结构上可能具有相似性。

Figure 5.  Intra-species (left)and inter-species (right)synteny analysis of the LEA gene family in durian.

另外，对榴莲和拟南芥LEA基因家族进行了物种间的共线性分析（图5），共鉴定出46对同源基因，具体表现为：LEA1亚家族中，DzLEA1与AT5G06760.1；LEA2亚家族中，DzLEA4与AT2G44060.1，DzLEA19与AT1G01470.1、AT2G46140.1，DzLEA20与AT3G62580.1，DzLEA40与AT1G01470.1、AT2G46140.1；LEA3亚家族中，DzLEA21与AT4G02380.3、AT1G02820.1，DzLEA26与AT1G02820.1、AT4G02380.3，DzLEA32与AT1G02820.1、AT4G15910.1，DzLEA39与AT1G02820.1、AT4G15910.1、AT4G02380.3；LEA4亚家族中，DzLEA2与AT2G36640.1、AT3G53040.1，DzLEA12与AT2G36640.1、AT3G53040.1，DzLEA31与AT4G13560.2，DzLEA42与AT1G52680.1，DzLEA46与AT1G15415.1、AT1G52680.1；LEA5亚家族中，DzLEA13与AT2G40170.1，DzLEA15与AT2G40170.1、AT3G51810.1，DzLEA29与AT2G40170.1、AT5G05560.2；LEA6亚家族中，DzLEA11与AT2G23110.1、AT5G67350.1；DHN亚家族中，DzLEA5与AT1G76180.1、AT1G20440.1，DzLEA7与AT1G76180.1、AT1G20440.1、AT4G38410.1，DzLEA17与AT4G38410.1，DzLEA18与AT1G76180.1、AT1G20440.1；SMP亚家族中，DzLEA14与AT5G27980.1，DzLEA22与AT1G03120.1、AT3G22490.1，DzLEA27与AT1G03120.1、AT3G22490.1，DzLEA38与AT1G03120.1、AT3G22490.1。这表明这些基因可能具有某些特定功能和进化的保守性。

为了进一步研究LEA基因的具体功能，对榴莲DzLEAs基因进行基于超几何分布的GO和KEGG分析（图6）。GO富集分析结果显示，在生物进程类别中，6个LEA基因家族成员均在“对寒冷的反应”以及“对缺水的反应”显著富集，分别是DzLEA5、 DzLEA16、 DzLEA18、 DzLEA21、 DzLEA26、 DzLEA39。

Figure 6.  Bar graphs of GO (top)and KEGG (bottom)enrichment analysis for the durian LEA gene family.

在KEGG富集分析中，LEA基因家族成员主要富集在泛素介导的蛋白水解，剪接体、泛素系统，蛋白质的折叠、分选与降解，染色体及相关蛋白质功能通路上。

转录组学分析表明（图7），部分DzLEAs基因对低温胁迫具有较高敏感性，具体表现为，与常温下榴莲对照相比10℃低温胁迫1d条件下，猫山王榴莲品种有7个基因显著上调表达，3个基因显著下调表达，这其中两个基因（DzLEA21和DzLEA26）随着低温胁迫处理时间的增长，呈现持续显著上调表达，另有两个基因（DzLEA11和DzLEA29）随着低温胁迫处理时间的增长，呈现持续显著下调表达；与常温下榴莲对照相比10℃低温胁迫1d条件下，金枕榴莲品种有5个基因显著上调表达，3个基因显著下调表达，这其中3个基因（DzLEA18、DzLEA26和DzLEA44）随着低温胁迫处理时间的增长，呈现持续显著上调表达，另有2个基因（DzLEA11和DzLEA32）随着低温胁迫处理时间的增长，呈现持续显著下调表达。

Figure 7.  Transcriptomic expression profiles of LEA gene family members in durian cultivars Musang King (left)and Golden Pillow (right)under cold stress.

为了进一步分析DzLEAs基因在低温胁迫中的作用，基于低温处理下的榴莲叶片转录组数据，挑选出6个低温处理（10℃）下表达量升高的DzLEAs基因（图8）和4个非DzLEAs基因（图9），通过qRT-PCR测定榴莲叶片在四个低温处理时期下的表达量变化。结果显示，DzLEAs基因均在10℃低温处理下呈现不同程度上调表达，其中不同品种的DzLEA18和DzLEA39均在10℃低温处理1d表达量最高，而后表达量逐渐降低；猫山王品种DzLEA5表达量随低温处理时间增加而逐步上升，并在10℃低温处理7d表达量最高，金枕品种DzLEA5表达量呈现先升高而后降低趋势，在10℃低温处理3 d表达量最高；不同品种的DzLEA21表达量均随低温处理时间增加而逐步上升，且在10℃低温处理7 d表达量最高；猫山王品种DzLEA26表达量随低温处理时间增加而逐步上升，并在10℃低温处理7 d表达量最高，金枕品种DzLEA5表达量呈现先升高而后降低趋势，在10℃低温处理3 d表达量最高；不同品种的DzLEA45表达量均呈现先升高而后降低趋势，且在10℃低温处理5 d表达量最高。非DzLEAs基因在10℃低温处理下表达模式有所不同，evm.model.GRLG04.19和evm.model.GRLG20.232表达模式类似，在不同品种中均呈现先升高后降低而后再升高的趋势，且均在10℃低温处理7 d表达量最高；evm.model.GRLG25.635和evm.model.GRLG27.910表达模式类似，猫山王品种中呈现先升高而后降低的趋势，在10℃低温处理5 d表达量最高，金枕品种中呈现逐步升高的趋势，在10℃低温处理7 d表达量最高。综上，这些基因的表达模式在不同品种之间相近，且与转录组表达结果基本一致，呈现较为明显的上调趋势，表明这些基因可能在榴莲响应低温胁迫过程中发挥了关键作用。

Figure 8.  Relative expression levels of DzLEA genes in the two durian cultivars under cold stress.

Figure 9.  Relative expression levels of non-DzLEA genes in the two durian cultivars under cold stress.

本研究对榴莲LEA基因家族进行了鉴定和生物信息学分析，共鉴定了46个DzLEAs基因，分为8个亚族，分布在18条染色体上。榴莲LEA基因家族成员数量相较于茶树（48个）^[9]、甜橙（72个）^[18]和拟南芥（51个）^[6]经历了一定程度的收缩，相较于木薯（26个）^[19]、香蕉（23个）^[20]和石斛（17个）^[21]经历了一定程度的扩张。通过对其进行理化性质分析，结果显示，大多数DzLEA蛋白为亲水性蛋白（44个），这与前人的研究一致，高度亲水性有利植物在失水胁迫下高效捕获水分，从而保护细胞膜和蛋白质^[22]。且有研究表明，大多数LEA蛋白被预测为固有无序蛋白，且在干燥过程中形成两亲性α螺旋，这些特性能够使LEA蛋白在冷冻胁迫下阻止保护酶聚集和失活，并且提高膜脂的稳定性^[23]。对榴莲LEA基因家族的结构分析显示，2个成员缺乏内含子、20个成员只有1个内含子、20个成员有两个内含子、4个成员有3个或3个以上内含子，较少的内含子结构缩短了从转录到翻译的防御胁迫时间，这有利于植物对环境胁迫引起的刺激迅速作出反应^[24]。

在研究过程中发现，榴莲在10℃低温胁迫下，短时间（1d）就会发生大面积正常颜色叶片脱落情况。我们推测这与其长期进化过程中的生存策略有关，其生存策略是“快速生长、高效繁殖”，低温环境下，榴莲会迅速响应并引起脱落酸和乙烯水平增高，丢弃耗能器官（叶片），将资源集中于维持主干生命，等待环境回暖。榴莲叶片具有厚角质层、厚蜡质层和高糖代谢特征，且油脂成分主要为结构性脂质，不饱和脂肪酸比例较低，相变温度较高。其结构和成分特征与抵御热带地区强光、暴雨、高温和病虫害相适应，但表现出极强的冷敏感性，其厚角质层和蜡质层使其在面对低温时不易回温，易发生光合系统受损和膜脂过氧化。

对榴莲LEA基因家族进行生物信息学分析过程中，发现4个基因（DzLEA5、DzLEA21、DzLEA26和DzLEA39）启动子上均存在与低温响应相关的元件LTR，且在GO分析中，在生物进程类别中，该4个基因均与“对寒冷的反应”以及“对缺水的反应”有显著关联。将该四个基因与转录组中低温胁迫下显著上调表达的两个基因（DzLEA18和DzLEA45）共同进行表达模式分析。发现这6个基因均在低温胁迫条件下呈现不同程度的上调表达，表明其可作为后续研究榴莲低温胁迫响应的分子机制及榴莲抗寒品种选育的候选基因。针对这些基因，后续可通过在模式植物（拟南芥）异源过表达，结合表型和生理指标数据进行进一步验证；并通过多品种的重测序，发掘其单核苷酸多态性（SNP）或插入/缺失（InDel）变异，开发与低温响应性状关联的分子标记。这些标记可用于榴莲种质资源的抗寒性鉴定，加速抗寒品种的筛选与选育进程。