萱草根茎低温胁迫转录组分析

黄长兵，程培蕾，杨绍宗，张焕朝，姜正之，金立敏

(1.南京林业大学 南方现代林业协同创新中心/林学院，江苏 南京 210037； 2.苏州农业职业技术学院 园艺科技学院，江苏 苏州 215008； 3.苏州市华冠园创园艺科技有限公司，江苏 苏州 215557)

萱草(Hemerocallisfulva)又被称为蘐草、忘忧等，属于百合科萱草属多年生植物，其根茎较短但其纺锤形肉质根较厚[1-3]，在中国有着悠久的栽培历史。萱草既可以作为园林植物作观赏之用，也可以作为蔬菜食用，亦可以作为药物治病；由于其品种多样，色彩丰富，抗逆性强等因素，成为道路绿化、公园绿地、屋顶花园等生态城市建设广泛使用的园林绿化材料[4]。作为萱草属中较少几种可食用植物，萱草拥有较高的经济价值。此外，萱草中富含大量次级代谢产物，这些代谢物对于人体的焦虑、肿胀等症状具有良好的治疗效果[5-6]。目前，研究者们已经对萱草的生理生化、分类鉴定、组织培养和育种等领域进行了较多研究[7-11]，但有关于萱草逆境胁迫的报道和研究相对较少。低温胁迫作为最常见的非生物胁迫因素之一，能够导致植物的生长和发育出现延迟，同时也会影响作物的产量[12-14]。冷害是影响萱草观赏价值和经济价值的重要因素，因此，研究萱草在低温胁迫下基因的表达变化对于解决萱草的冷害具有重要意义。

转录组测序与生物信息学分析技术相结合可以大大促进功能基因组学的研究，通过测序得到一定条件下细胞的整体转录信息，为从转录水平来发掘生物体的生长发育、次级代谢与调控机制[15-16]提供了便利。近些年，研究者们用转录组分析策略来研究植物抗低温胁迫基因组学，已经获得拟南芥[17]、番茄[18]、紫花苜蓿[19]、黄瓜[20]和水稻[21]等众多植物转录组数据，这些都为阐明植物抗低温胁迫机制研究奠定了基础。本研究利用高通量测序平台Illumina HiSeq对低温胁迫条件下2种萱草的根茎进行转录组测序分析，期望揭示萱草根茎转录组在低温胁迫条件下发生的表达差异性特征，为萱草冷应激功能基因的鉴定、次生代谢途径的解析及其低温胁迫应激调控机制的研究提供参考。

1 材料与方法

1.1 材料

萱草Jinyan 7(耐冷型品种，冬季地上部常绿)和Lucretius 4(冷敏感型品种，冬季地上部枯萎)由苏州农业职业技术学院相城科技园提供。在通风温室中(20 ℃，光周期：16 h光照8 h黑暗)，植物于塑料盆中生长至6片叶子，然后移入生长室进行处理。将温度从20 ℃逐渐降低至-4 ℃开始低温处理，具体方法：让植株20 ℃生长48 h，然后分别在10 ℃和4 ℃诱导48 h，最后将它们在-4 ℃处理8 h。收集在20 ℃(作为对照组，CK)和-4 ℃(作为处理组，T)生长的萱草根茎，在液氮中冷冻用于RNA提取和生理分析。收集的4个样品分别命名为4CK(20 ℃时Lucretius 4的根茎)、4T(-4 ℃处理的Lucretius 4根茎)、7CK(20 ℃时Jinyan 7根茎)、7T(-4 ℃处理的Jinyan 7根茎)。每3株植物根茎作为1个生物学重复，共设置2组重复。

1.2 RNA提取与文库构建

采用TRNzol Universal Reagent植物总RNA提取试剂盒，分别提取冷处理前后耐冷型和冷敏感型萱草根茎总RNA，Agilent 2100 Bioanalyzer检测总RNA的质量，NanoDrop TM2000检测总RNA的浓度。利用Oligo(dT)磁珠富集每个萱草文库的mRNA，加入适量打断剂高温下片断化mRNA，并以片段化的mRNA为模板合成cDNA，选择合适的片段进行PCR扩增，构建文库。

1.3 转录组测序与组装

将构建的萱草根茎转录组文库提交Illumina HiSeq平台进行高通量测序。首先对Illumina测序的原始reads进行处理，去除接头污染、未知碱基N含量偏高(>5%)和低质量reads，然后将获得的高质量clean reads利用Trinity[22]进行denovo组装分析，最后通过Tgicl对其进行聚类和去冗余处理，得到unigenes。

1.4 转录组功能注释

利用Blastn对得到的unigene序列进行Nt数据库注释，通过Blastx将unigenes序列与以下蛋白数据库进行比对(E-value<1×10-5，E-value为比对的期望值)：Nr(NCBI nonredundant protein)、SwissPort、KOG(clusters of orthologous groups)和KEGG(kyoto encyclopedia of genes and genomes)，得到和相应unigenes序列相似性最高的蛋白序列，从而得到unigenes的对应注释信息。使用Blast2GO软件，同时结合Nr注释得到unigenes的GO(gene ontology)功能注释。

1.5 蛋白编码框和转录因子预测

通过软件Transdecoder(https://transdecoder.github.io)来识别unigenes的候选编码区域，具体过程主要包括通过TransDecoder.LongOrfs软件从unigenes中将最长的开放阅读框(open reading frame，ORF)提取出来，接着将其与SwissProt数据库进行Blast比对，其结果经Hmmscan搜索得到Pfam蛋白的同源序列，从而实现蛋白编码框(coding sequence，CDS)的预测。通过getorf软件来检测unigenes的ORF，接着利用Hmmsearch检索将ORF与PlantTFDB数据库中转录因子蛋白结构域进行比对，最后根据PlantTFDB中关于转录因子家族特征的描述对unigenes进行鉴定。

1.6 Unigenes基因表达水平分析

通过MISA(http://pgrc.ipk-gatersleben.de/misa)软件来检测萱草转录组的unigenes，搜索和验证简单重复序列(simple sequence repeats， SSR)，并进行相关统计分析。利用软件Bowtie2将clean reads与unigenes进行比对，然后通过RSEM软件计算得到各个样品的基因表达水平，并参照Saddhe等[23]的差异基因检测方法，通过软件PossionDis [fold change≥2.00和错误发现率(FDR)≤0.001]得到不同样本之间的差异基因。

1.7 差异表达基因聚类、GO功能、Pathway功能分析

差异基因的层次聚类分析采用R软件中的pheatmap函数分析法。依据GO功能的官方分类对差异基因进行功能注释，根据KEGG的官方分类对差异基因的生物通路进行注释分析，差异基因富集分析则利用R软件中的phyper函数分析法。

1.8 差异表达基因的qRT-PCR验证

通过对转录组中差异表达基因的进一步筛选，随机选择8个显著差异表达的基因，进行qRT-PCR验证。qRT-PCR检测使用Takara公司的SYBR®Premix Ex TaqTM试剂，在荧光定量仪中进行。内参基因为18S rRNA，正反向引物序列分别为：5′-CCATAAACGATGCCGACC-3′；5′-TTGCGACCATACTCCCCC-3′。PCR反应条件为：95 ℃ 5 min；95 ℃ 5 s，60 ℃ 30 s，72 ℃ 10 s，40个循环。qRT-PCR中使用的引物序列如表1所示。对每个基因和样本进行3个生物重复和3个技术重复。用2-ΔΔCt法计算差异基因相对表达水平。

2 结果与分析

2.1 萱草转录组组装与质量分析

从表2可以看出，4个样品4CK、4T、7CK、7T分别得到47 738 570、50 110 130、50 918 552、38 540 808条raw reads，过滤后分别产生46 171 660、47 791 962、47 411 290、35 880 734条高质量clean reads，各包含6 858 685 119、7 099 131 160、7 047 333 003、5 342 671 114个核苷酸信息，Q20(质量值≥20的碱基所占百分比)分别为98.61%、98.63%、98.72%、98.82%，Q30(质量值≥30的碱基所占百分比)分别为95.65%、95.7%、95.97%、96.23%，GC含量分别为48.83%、48.89%、48.57%、48.48%，表明clean reads的质量符合实验要求，可用于后续分析。样品4CK、4T、7CK、7T经Trinity组装后各获得55 348、51 400、5 6962、49 915个unigene，平均长度分别为796、824、763、753 nt，N50(按转录本长度从大到小排序后逐个累加至所有转录本总长度的50%)分别为1 294、1 337、1 225、1 197 nt。

表2 有效数据评估统计Table 2 Valid data evaluation statistics

Unigenes的长度分布图(图1)显示，Unigenes长度在300～>3 000 nt；长度为300 nt的unigenes数量最多，共25 011条；400 nt的unigenes数量次之，共11 341条；其他长度的unigenes均低于10 000条，大于3 000 nt的unigenes共2 194条。

2.2 萱草转录组unigenes的功能注释

为获得组装好的unigenes的表达和功能注释，将unigene序列与Nr、Nt、Swissprot、KEGG、KOG等数据库进行比对。结果表明，通过Nt数据库可对10 453条unigenes(12.00%)进行成功比对并进行注释，Nr数据库可对38 447条unigenes(44.15%)进行成功比对并进行注释，在SwissPort、KOG、KEGG、Pfam、GO数据库中分别有28 095(32.27%)、41 418(47.57%)、27 835(31.97%)、36 767(42.23%)、15 933(18.30%)条unigenes获得成功比对和注释。

对Nr数据库的比对注释结果进行分析，图2是unigenes注释同源基因的物种分布图，可以发现，在相似序列较高的物种中，芦笋(Asparagusofficinalis)所占比例最高，其数量可以达到24 115条(45.71%)；其次为油棕(Elaeisguineensis)、海枣(Phoenixdactylifera)和凤梨(Ananascomosus)，其数量分别为5 695条(10.79%)、3 673条(6.96%)和1 638条(3.1%)；注释到其他物种的unigenes共17 633条，占比为33.43%。

2.3 CDS和转录因子分析

通过BLAST比对分析，从4组萱草样品转录组的全部unigenes中共获得42 133个CDS 序列，CDS的长度全部大于200 nt，其中有16 337条CDS长度超过1 000 nt，3 923条CDS长度超过2 000 nt(图3)。通过转录因子预测分析可以发现，这些unigenes分属于56个家族(图4)，而其中Homeobox、zf-C2H2、bHLH、HMG、NGFIB-like、MYB、zf-GATA、ZBTB、TF_bZIP与CG-1类占主体，这说明众多转录因子参与萱草根茎的生理代谢。

图3 CDS长度分布图Fig.3 Distribution diagram of coding sequence length

图4 转录因子家族分类图Fig.4 Classification of transcription factor families

2.4 差异表达基因统计、GO和KEGG分析

表3 Unigenes注释结果统计Table 3 Annotation results of unigene sequences

2.4.1 差异表达基因统计

对冷敏感型和耐冷型萱草冷处理后的差异基因数量进行统计发现，与对照相比，冷敏感型萱草根茎在低温胁迫下，总共检测到6 153个基因表达变化，其中，上调表达基因有1 849个，下调表达基因共检测到4 304个。而耐冷萱草根茎在低温胁迫下，共有14 813个差异表达基因，其中，有6786个基因上调表达，8 027个基因表达量下调。对照组中，与冷敏感型萱草相比，耐冷萱草有10 282个基因上调表达，12 234个基因下调表达；而在低温处理条件下，与冷敏感型萱草相比，耐冷萱草上调表达基因有9 224个，下调表达基因有10 614个(图5)。

图5 上调与下调差异基因数量统计Fig.5 Numbers of up-regulated and down-regulated differentially expressed genes

2.4.2 差异表达基因GO富集分析

为了进一步了解低温胁迫条件下萱草的抗逆机制，对差异表达基因进行GO富集分析。GO分析结果表明(图6)，这些差异表达基因分为3大类：分子功能(molecular function)、细胞组分(cellular component)和生物学过程(biological process)，每个大类对应的功能组数量分别为16、19、29。差异基因大部分富集在生物学过程中，而参与分子功能的差异基因相对较少。在生物学过程中，富集较明显的包括细胞过程(7 216)、代谢过程(6 652)、生物调节(2 138)、生物过程调节(1 973)和定位(1 700)；参与细胞组分的差异基因主要富集在细胞(6 705)、细胞组分(6 614)、膜(5 386)、细胞器(4 947)和膜部分(4 934)；而参与分子功能的差异基因富集类别主要有催化活性(8 314)、结合(7 916)和转运活性(1 916)。

图6 差异基因GO功能分类图Fig.6 Gene ontology classification of differentially expressed genes

2.4.3 差异表达基因KEGG分析

根据KEGG代谢通路对差异基因进行功能注释，发现参与的代谢通路包含6个分支：细胞过程(cellular processes)、环境信息处理(environmental information processing)、遗传信息处理(genetic information processing)、人类疾病(human diseases)、代谢(metabolism)和有机系统(organismal systems)，具体的差异基因数量见图7。在代谢途径中，占据前5的类别分别为全局概述图(5 527)、碳水化合物代谢(1 694)、核苷酸代谢(1 673)、氨基酸代谢(1 311)和脂质代谢(1 183)；在遗传与信息处理类别中，翻译(2 612)、转录(2 165)和折叠分拣与降解(1 781)差异基因显著富集，而信号转导(5 416)在环境信息处理通路中差异基因富集较多。

图7 差异表达基因KEGG分类Fig.7 KEGG classification of differentially expressed genes

2.4.4 萱草低温胁迫相关代谢通路分析

对差异基因的代谢通路进行分析发现，差异基因富集程度比较明显的途径主要包括：异黄酮生物合成(isoflavonoid biosynthesis)、牛磺酸和亚牛磺酸代谢(taurine and hypotaurine metabolism)、光合天线作用蛋白(photosynthesis-antenna proteins)、二苯基庚酮和姜酚生物合成(stilbenoid， diarylheptanoid and gingerol biosynthesis)、植物昼夜节律(circadian rhythm-plant)、不饱和脂肪酸的生物合成(biosynthesis of unsaturated fatty acids)等，具体的代谢通路与差异基因数量见表4。

表4 与萱草低温胁迫作用相关的16条代谢通路Table 4 Sixteen metabolic pathways related to cold stress of Hemerocallis fulva

2.5 差异基因参与对低温的反应

通过对差异基因的进一步筛选发现，Ca2+信号通路相关基因表达发生明显差异(表5)。与对照组相比，低温处理后，冷敏感型萱草中CBL互作丝氨酸/苏氨酸蛋白激酶3基因显著上调；

表5 钙离子信号通路相关差异基因Table 5 Differentially expressed genes involved in Ca2+ signalling pathway in each comparison

在耐冷型萱草中多种Ca2+信号通路相关基因明显差异表达，钙依赖性蛋白激酶基因CAMK2A、钙调蛋白基因(Calm5、Calm3)表达显著上调，而钙调蛋白和肌联蛋白结合RhoGEF基因、CBL互作丝氨酸/苏氨酸蛋白激酶14基因则显著下调。对照组中，与耐冷型萱草相比，冷敏感型萱草钙调蛋白基因(Calm5、Calm3)、CBL互作丝氨酸/苏氨酸蛋白激酶14基因显著上调，钙依赖性蛋白激酶基因CAMK2A显著下调。在低温胁迫下，与耐冷型萱草相比，冷敏感型萱草钙调蛋白基因(Calm5、Calm3)、钙依赖性蛋白激酶基因CAMK2A与对照组表达趋势相同。CL7725.Contig5_All基因

编码了钙调蛋白和肌联蛋白结合RhoGEF(calmodulin and titin-interacting RhoGEF)。低温胁迫下，与对照组相比，耐冷型萱草中CL7725.Contig5_All基因显著上调；与耐冷型萱草相比，冷敏感型萱草的CL7725.Contig5_All基因显著下调；与对照组相比，冷敏感型萱草的CL7725.Contig5_All基因表达无显著差异。对照组中，耐冷型萱草与冷敏感型萱草的CL7725.Contig5_All基因表达也没有显著差异。Unigene6945_All基因编码钙调素结合转录激活子3(calmodulin-binding transcription activator 3)，低温胁迫抑制了其在耐冷萱草中的表达，在对照组中，该基因在2种萱草中的表达没有差异。

此外，也筛选到MAPK信号通路相关差异基因(表6)，低温处理前后比对结果发现，在冷敏感型萱草(4CK-VS-4T)中，促分裂原活化的蛋白激酶基因MAPKKKK3显著上调，MAPK1-like、MAPKKKK5显著下调。在耐冷型萱草中(7CK-VS-7T)，促分裂原活化的蛋白激酶基因MLK4-like、MAPKKKK3、MAPK1、MAPK4显著上调，MAPK1-like、MAPK7-like、MAPKKKK5显著下调。对照组中，与冷敏感型萱草相比，MAPK8互作蛋白2基因(MAPK8IP2)、MAPKKKK3基因、MAPK7-like基因上调表达，MAPK4、MAPKKKK5基因下调表达。低温处理后，耐冷型萱草与冷敏感型萱草的MAPKKKK3、MAPK8IP2基因表达上调，MAPK7-like、MLK4-like、MAPKKKK5表达下调。进一步分析发现，低温可以抑制MAPK1基因(CL4304.Contig1_All)的表达，常温下，2种萱草中没有检测到该基因的表达差异，而在低温下，冷敏感萱草中，其下调倍数为-1.36，在耐冷型萱草中，其表达被抑制的更加明显，下调倍数为-3.90。

表6 MAPK信号通路相关差异基因Table 6 Differentially expressed genes involved in MAPK cascades pathway in each comparison

对编码热激蛋白(表7)相关的基因进行筛选发现，冷敏感型萱草(4CK-VS-4T)中，热激蛋白基因HSP70-2、DNAJC6、HSP90显著上调，而在耐冷型萱草中(7CK-VS-7T)，仅热激蛋白基因HSP86上调表达，其他热激蛋白基因HSP90-5、HSP90均显著下调。不同的温度条件下，对2种萱草进行对比发现，对照条件下，与冷敏感型萱草相比，Hspa8基因上调表达，HSP70-2、HSP90-5基因下调表达，低温处理后，2种萱草的对比中筛选到的热激蛋白基因均下调表达。HSP90基因(Unigene6560_All)在2种萱草中表现出相反的表达趋势，冷敏感萱草中，低温处理后上调，而在耐冷型萱草中，与对照条件相比下调-1.61，与冷敏感萱草相比下调-2.61。

表7 编码热激蛋白的差异表达基因Table 7 Differentially expressed genes encoding HSPs protein in each comparison

对差异表达的转录因子(表8)进行分析发现，在冷敏感型萱草(4CK-VS-4T)中，转录因子WRKY11、MYB102显著下调表达，而在耐冷型萱草(7CK-VS-7T)中，转录因子MYB20、bHLH2上调表达，MYB102下调表达。与冷敏感萱草相比，对照组耐冷型萱草WRKY11上调表达，bHLH2、4CK-VS-4T，Lucretius 4低温与对照相比；7CK-VS-7T，Jinyan 7低温和对照相比；4CK-VS-7CK，对照条件下Jinyan 7与Lucretius 4相比；4T-VS-7T，低温处理下Jinyan 7与Lucretius 4相比。下同。

表8 差异表达的转录因子Table 8 Differential expression of transcription factors

2.6 荧光定量PCR验证

为了验证转录组测序结果，随机选择8个表达水平差异显著的基因进行qRT-PCR试验，分别为：类硫氧还蛋白基因(thioredoxin-like)、丝氨酸/精氨酸重复基质蛋白基因(serine/arginine repetitive matrix protein 3-like)、溶磷酰胆碱酰转移酶1基因(Lysophos phatidylcholine acyltransferase 1)、1-磷脂酰肌醇 3-磷酸 5-激酶基因(1-phosphatidylinositol 3-phosphate 5-kinase)、粘液素-2基因(Mucin-2)、肌球蛋白VIIa基因(myosin-VIIa-like)、

含有LOB结构域的蛋白38基因(LOB domain-containing protein 38-like)、多聚泛素-A基因(polyubiquitin-a)。结果显示：低温胁迫下，8个基因的qRT-PCR结果与转录组分析的趋势一致(图8)，只是具体表达倍数上存在差异，这些差异可能是由于转录组测序和qRT-PCR实验的检测灵敏度，以及实验结果分析方法不同造成的，表明转录组测序结果可靠。

3 结论与讨论

低温是许多植物所面临的常见威胁，其不仅可诱导细胞脱水和膜损伤[23-24]，还会改变植物细胞的正常生理代谢过程，主要表现在2个方面：呼吸作用障碍[25]、活性氧的积累和去除[26-28]。萱草是一种在许多地区广泛种植的营养植物，耐寒品种的选育对西北地区的营养植物生长具有重要意义。低温胁迫通过抑制生长和发育严重限制了植物物种的地理分布，植物可以通过调节基因表达来适应低温胁迫。抗低温胁迫的植物可以通过防御机制调节自身，包括生理反应、代谢和信号转导途径的适应，以及基因表达变化[29]。本研究对2种不同抗性萱草在低温胁迫条件下的转录组进行分析发现，在7个数据库同时注释的unigenes共有6 727条，至少有1种数据库注释成功的unigenes共有56 679条，其

qRT-PCR结果中的误差线表示平均值的标准差(n=3)。The error line in results represents the standard deviation of the mean value (n=3).图8 差异表达基因的qRT-PCR验证Fig.8 qRT-PCR validation for differentially expressed genes identified by RNA-Seq

比例达到65.09%；还有30 394条unigenes未获得注释，这一比例和多种植物的测定结果相类似[30-32]，说明萱草转录组中存在着大量特征序列，以及大量unigenes功能有待开发。与对照组相比，冷敏感型萱草根茎在低温胁迫下，上调表达基因有1 849条，下调表达基因有4 304条；而耐冷萱草根茎在低温胁迫下，有6 786条基因上调表达，8 027条基因下调表达，差异表达基因数量显著多于冷敏感型萱草品种，猜测其抗冻性可能与迅速调控多种基因差异表达有关。差异表达基因的GO分类结果显示，2种萱草根茎的DEG功能特性与生物过程、分子功能和细胞组分相关。2种萱草低温胁迫和非低温胁迫条件下差异表达基因均涉及KEGG代谢通路的五大分支。

次生代谢物在植物逆境胁迫下的合成和积累具有重要意义。已有研究发现，低温等逆境胁迫能够诱导类黄酮化合物的积累，提高植物抗逆性[33]；一些植物在低温胁迫下体内类黄酮的生物合成会发生显著变化[34]。本研究也发现，低温胁迫下萱草根茎的类黄酮合成基因表达水平显著高于对照组，这说明低温胁迫在一定程度上可能诱导了类黄酮的积累。花青素、吲哚生物碱等物质的合成也会受到低温胁迫的影响。

低温刺激植物体内Ca2+信号通路的变化，Ca2+信号通路被三大类蛋白质感知：CDPKs、CaMs和CBLs[35]。在水稻中，低温胁迫下，钙调蛋白OsMSR2基因在不同发育阶段不同组织均显著上调表达；拟南芥中外源表达该基因可以影响其耐盐性和耐旱性[36]。在耐冷型萱草中，2个钙调蛋白基因Calm5/Calm3均显著上调表达，而在冷敏感型萱草中，没有发现相应基因的变化。在低温处理的茶树中，5个钙调蛋白基因、2个CDPK基因和1个CBL基因参与了信号转导[37]。本研究中，CDPK3基因表达在冷敏感型萱草中显著上调，耐冷型萱草中CDPK14基因表达显著下调。Ca2+信号通路相关基因表达在2种萱草中差异较显著，其中calmodulin and titin-interacting RhoGEF(CL7725.Contig5_All)与calmodulin-binding transcription activator 3(Unigene6945_All)基因均在低温胁迫的耐冷萱草中被显著抑制，而在冷敏感萱草中无明显差异，猜测这2个基因在Ca2+信号通路对萱草耐冷调控中发挥重要作用。MAPKs信号通路也影响植物非生物胁迫响应，许多受体蛋白激酶(RLKs)负责感知外部环境的变化和相应的信号转导[38]。在冷敏感型萱草中，MAPKKKK3上调表达，MAPKKKK5下调表达；在耐冷型萱草中也发现这2个基因有相同的表达趋势，且上下调程度更显著。冷敏感型萱草中仅发现1个MAPK1基因下调表达，而耐冷型萱草中共检测到3个差异表达的MAPK基因(MAPK1、MAPK4、MAPK7-like)和1个差异表达的MLK基因。这些发现为Ca2+在萱草低温响应中的关键作用提供了证据，同时进一步证明MAPK通路影响萱草对低温胁迫的响应。MAPKs信号通路基因MAPK1(CL4304.Contig1_All)的表达受低温抑制，在耐冷型萱草中的抑制程度明显高于冷敏感萱草，表明MAPK1基因可能作为MAPKs信号通路的关键基因发挥作用。

转录因子在植物非生物胁迫应答中发挥重要作用，在拟南芥中过量表达水稻MYB4(R2R3-型MYB转录因子)和OsMYB3R-2(R1R2R3型MYB转录因子)基因，可显著提高拟南芥抗冻性[39]。冷敏感型萱草中，MYB家族转录因子MYB102显著下调表达，该基因在耐冷型萱草中下调倍数较低。此外，另一个MYB家族转录因子MYB20仅在耐冷型萱草中上调表达，在冷敏感型萱草未发现该基因表达变化。研究表明，WRKY和bHLH转录因子家族与植物低温胁迫密切相关[40-41]。bHLH转录因子家族基因bHLH2仅在耐冷型萱草低温处理后上调表达，而WRKY转录因子家族基因WRKY11在冷敏感型萱草中下调表达；与冷敏感型萱草相比，耐冷萱草中的WRKY11在常温和低温条件表达水平均较高，WRKY11在2种萱草中的表达差异及其在应对低温胁迫中的作用机理需要进一步探讨。以上结果表明，MYB、WRKY、bHLH转录因子家族在冷敏感型萱草和耐冷型萱草低温胁迫响应中发挥了重要而不同的作用。除了蛋白激酶和转录因子，一些功能蛋白也参与了低温胁迫过程。热休克蛋白(HSPs)最初是通过参与对高温的反应而被识别的，越来越多的研究发现，许多热休克蛋白也对低温有反应。在冷敏感型萱草和耐冷型萱草中均发现一些HSP蛋白基因呈现差异表达，但在2种萱草中筛选到的HSP家族成员不同。有趣的是，在2种萱草中均发生变化的HSP90，在冷敏感型萱草中上调表达，而在耐冷型萱草中则下调表达，推测HSP90可能是萱草响应低温胁迫的重要候选基因，其在萱草低温胁迫中的作用还需要进一步研究。

综上所述，本研究利用高通量测序鉴定和比较了萱草Jinyan 7(耐冷型品种)和Lucretius 4(冷敏感型品种)的低温响应基因，发现在低温胁迫下，Ca2+信号通路基因、MAPKs信号通路基因、转录因子和热激蛋白基因在2种萱草中表现出不同的变化趋势，这些关键路径基因表达的差异可能导致了2种萱草在低温胁迫中的不同反应。本研究为萱草低温响应的分子机制提供了思路，同时，筛选出的差异基因也可能成为培育耐冷型萱草品种的候选基因。