【研究意义】百合（Lilium）是国家卫生部首批公布的药食同源植物，也是著名的观赏植物［1］。食用百合是营养价值很高的特种蔬菜，其鳞茎含碳水化合物、蛋白质、果胶、脂肪、钾、磷以及钙、铁、百合甙等抗癌性植物碱和多种维生素，兼具保健功能，产品畅销中国及东南亚地区，在世界范围享有盛誉，市场前景十分广阔。兰州百合〔Lilium davidii Duch.var.unicolor(Hong) Cotton〕微甜、口感好，是食用百合中的精品，生长在 103°54´ 54´´ E、36°06´ 12´´ N、海拔2 000多米的七里河区，是食用百合中的耐寒品种代表［2］。温度是植物生长最重要的环境因子之一，通过对兰州百合抗冻转录组的研究，筛选出百合抗冻的基因及抗冻相关通路，为后续兰州百合SSR标记开发、抗寒性相关基因的发掘及分子育种提供借鉴。
【前人研究进展】目前，对百合抗寒性研究集中于亚洲百合［3］、东方百合［4］、卷丹百合［5］等观赏百合品种，对百合耐寒的研究主要集中在细胞膜透性和保护酶系统方面，相对电导率、脂膜氧化产物丙二醛、抗氧化酶类（SOD、CAT、POD）活性通常被作为植物耐寒性生理指标［6-7］。HOSHI等［8］对百合在农杆菌介导下转录组进行研究，建立了东方百合杂交品种的转基因生产体系；杜方等［9］对百合不同器官的转录数据进行测序并开发了SSR标记。
【本研究切入点】关于抗寒食用百合代表——兰州百合的抗寒转录因子研究至今仍是空白。本研究运用先进的Illumina高通量测序平台的转录组测序技术，分别对正常温度条件下和零下超低温条件下生长的兰州百合苗期叶片的转录产物mRNA进行测序，比较常温与零低温下差异表达的基因及通路，了解兰州百合超低温胁迫下的基因表达状况。【拟解决的关键问题】进一步对差异表达基因进行GO、KEGG(Kyoto Encyclopedia of Genes and Genomes)注释和富集分析，以期进一步筛选出与兰州百合抗冻相关的基因［10］。
1 材料与方法
1.1 试验材料
在兰州市七里河区采集露地自然生长成熟的兰州百合种球，贮藏于2℃冰箱110 d［11］后，种植于15 cm口径的花盆中，置于20℃智能低温人工气候箱（RXZ-0288型）中培养，相同的水肥条件管理。将生长至15 cm左右的百合幼苗放在智能低温人工气候箱中进行-4℃低温处理［12］，降至目标温度后维持24 h［13-14］，以室温20℃作为对照，每个处理3次重复。
1.2 试验方法
1.2.1 转录测序 分别取3份20℃与-4℃的兰州百合叶片样品，采用TRIzol法提取总RNA［15］，并对RNA浓度、纯度、完整性进行检测。样品检测合格后，用带有Oligo (dT)的磁珠富集真核生物mRNA［16］。随后加入fragmentation buffer将mRNA打断成短片段，以打断后的 mRNA为模板合成一链cDNA，再合成二链cDNA。然后纯化、修复，再用AMPure XP beads进行片段大小选择，最后进行PCR扩增，纯化PCR产物，得到最终文库。质检合格后使用Illumina 4000高通量测序平台(HiSeq/MiSeq)进行测序［16］。由Illumina 4000高通量测序平台测序所得的数据称为 raw reads或raw data。然后进行转录结果组装与拼接，并进行基因功能注释与表达水平分析。
1.2.2 实时荧光定量PCR验证 随机选取10个差异性表达的基因进行qRT-PCR分析验证。以20 ℃与-4 ℃的兰州百合叶片RNA为材料，采用HiScriptⅡQRT SuperMix for qPCR试剂盒反转录合成cDNA，设计引物（表1）进行qRT-PCR，以Actin基因作为内参［17］。各处理均设3次重复，计算基因的相对表达量。
表1 定量PCR引物序列
Table 1 Primer sequences for qRT-PCR
基因ID Gene ID正向引物(5′→3′)Forward primer(5′→3′)反向引物(5′→3′)Reverse primer(5′→3′)c166605_g1 TGGTCTGGATGTGTTGCG GTCATAATTGAGTTTGTGGGAG c199368_g1 GCGGGCCAGCTCCTTGTGAA GCGGGTTCAGCGTGCAGGA c171769_g1 GCAACATCATCAACTCTGG ATGGAAGCAATCACCTAAGT c164813_g2 TCCAGTCCACCACAGTCT GCCAGCAGAACAGTCAAG c153721_g1 GTTGAGCCAGCCTACATC CTAGACTCCGCACTCTGT c76799_g1 TGGAGCAATTAGGCCCGGATGA CACCAGAGCCAGGCCCAACACT c153860_g1 GGCGCGAGTACGAGCTGGTGA CGCCTCGTCGGTCTTGTCGG c149445_g1 TGAAAAGGTGAGGTCACCAG GAACAGAGATGAAAGAAAGTTGG c115377_g1 GCTTGTGGCTGTGTAGATTAAG TCCACAACCTATTATCACCTCT c155535_g1 GCAAGCAAAGCCTTGGATC ATTGCCGGAGATCAAACG Actin TGAGCACATTCCAGCAGA CCATAGACAAAGCCATCG
2 结果与分析
2.1 测序结果与Trinity拼接
兰州百合20℃转录组共得到56 882 837 raw reads，56 456 634 clean reads，Q20值为94.56%，Q30值为87.32%，序列碱基GC含量为44.84%；-4℃转录组共得到58 092 923 raw reads，57 800 442 clean reads，Q20值为94.51%，Q30值为87.32%，序列碱基GC含量为44.54%。表明兰州百合转录组测序质量较高，可为后续的数据组装与分析提供良好基础。经过Trinity 转录本拼接和CD-HIT去冗余序列处理结果见表2。Unigene序列长度＜300的最多（图1），说明转录拼接较好，利于后续分析研究。
表2 处理前后序列数量、N50、N75、N90长度和总碱基数统计
Table 2 Sequence number, length of N50, N75, N90 and total base number by different data processing
注： N50、N75、N90分别为不小于总长50%、75%、90%的拼接转录本的长度。
Note: N50, N75 and N90 represent length of the transcripts of not less than 50%, 75% and 90%, respectively.
序列Sequence转录本/unigene长度Transcript /unigene length N50 N75 N90 总碱基数Total base number最短Shortest 中间值Median 最长Longest Transcript 201 341 21073 829 376 252 168637667 Unigene 201 327 21073 757 350 246 133179110
图1 拼接后Unigene 长度分布
Fig.1 Length distribution of all unigenes after splicing
2.2 基因功能注释结果
对获得的转录组序列进行数据库比对注释可知，注释的Unigene数量共238 301个，Nr（NCBI 蛋白质非冗余数据库）注释的Unigene数量为65 314个、占27.41%，KOG（真核生物直系同源基因数据库）注释的Unigene数量为42 138个、占17.68%，GO库注释的Unigene数量为35 244个、占14.79%，Swiss-Prot注释的Unigene数量为46 337个、占19.44%，Pfam注释的Unigene 数量为58 605个、占24.59%（表3）。
2.2.1 GO 分类 对Unigene进行GO注释后，按照成功注释到生物学过程（Biological Process，BP）、细胞组分（Cellular Component，CC）和分子功能（Molecular Function，MF）的基因进行统计。其中生物学过程中较多Unigene注释有高分子代谢过程28 326个、初级代谢过程16 483个、细胞代谢过程38 883个；分子功能中较多Unigene注释有转移酶活性5 090个、水解酶活性5 153个、跨膜转运功能3 069个；细胞组分中较多Unigene注释有细胞壁8 663个、细胞膜7 832个、细胞器9 581个（图2）。
表3 Unigenes 注释结果统计
Table 3 Annotation statistics for all Unigenes
数据库Database基因数量Number of Unigenes占比Proportion (%)NR 65314 27.41 NT 37916 15.91 KO 26291 11.03 SwissProt 46337 19.44 Pfam 58605 24.59 GO 35244 14.79 COG/KOG 42138 17.68所有数据库都有注释All databases 9695 4.07至少在一个数据库里有注释At least one database 84559 35.48合计Total 238301 100
图2 GO分类
Fig.2 GO function classification
2.2.2 KEGG分类 由表3可知，兰州百合测序结果有26 291条Unigene注释到KEGG数据库，涉及5大类代谢通路。其中代谢作用中注释4 914条、占40.66%；有机系统中注释905条、占7.49%；细胞过程中注释923条、占7.64%；遗传信息处理中注释3 609条、占29.86%；环境信息处理中注释1 734条、占14.35%。
2.3 基因表达差异分析结果
2.3.1 差异基因GO富集分析 将所有差异表达基因映射到Gene Ontology 数据库的各个功能分类中，计算每个分类的基因数目，与整体基因组背景相比，找出在差异表达基因中显著富集的功能分类。
对上调的差异基因进行富集分析，结果（图4）发现，在细胞成分中，差异基因富集效果较明显的有细胞核，基因占比13.81%；细胞内细胞器，基因占比16.06%；面上细胞器，基因占比16.06%。分子功能中差异基因富集效果较明显的功能有钙离子结合，基因占比3.52%；离子束缚，基因占比8.81%；阳离子绑定，基因占比8.81%；金属离子结合，基因占比8.81%。生物过程中上调基因富集效果较明显的功能有水响应，基因占比1.27%；非生物刺激反应，基因占比1.37%；磷酸化作用，基因占比8.62%。结果表明，含量增加的差异基因在细胞成分、分子功能、生物过程三方面主要调节物质的合成代谢。
图3 KEGG 代谢通路
Fig.3 KEGG metabolic pathway
对下调的差异基因富集分析，结果（图4）发现，差异基因主要富集在生物过程中。其中富集较为明显的是光合作用、光收获，基因占比1.65%；光合作用、光反应，基因占比1.65%。结果表明，含量相对减少的差异基因主要调节生物过程中的光合作用过程。
图4 差异基因（20℃ Vs -4℃）的GO富集结果
Fig.4 GO enrichment result of differential genes (20℃ Vs -4℃)
2.3.2 差异基因KEGG富集分析 差异基因KEGG富集结果显示，常温（20℃）与超低温（-4℃）下差异表达基因共有5 848个，其中上调的差异基因有3 478个，下调的差异基因有2 370个（图5）。经分析后发现，与蛋白激酶、蛋白磷酸酶、碳代谢以及ABA等相关的基因显著变化；此外，吲哚-3-甘油磷酸合成酶、蛋白磷酸酶、已糖激酶、钙结合蛋白、叶绿素a/b结合蛋白等相关基因显著上调或下调。
图5 差异表达基因（20℃ Vs -4℃）火山图
Fig.5 Differential expressed genes（20℃ Vs -4℃）
Pathway显著性富集分析以KEGG Pathway为单位，应用超几何检验，找出差异基因相对于所有注释的基因显著富集的Pathway。处理后含量增加的差异基因显著富集的通路有植物病原菌互作、氨基酸和核苷酸糖代谢、植物昼夜节律、亚油酸代谢等（图6）。如上调的差异基因c156798_g2、c163371_g2 (calcium-dependent protein kinase，EC:2.7.11.1)、c166423_g2、c145559_g1 (respiratory burst oxidase，EC:1.6.3.- 1.11.1)、c162064_g1、c164545_g2、c166855_g1(3.1027)c145934_g1(8.2812) c101051_g1(4.9356)、c144017_g1(Inf) c144017_g2(Inf) c159000_g1(2.9702)c140553_g1(9.9651) c170503_g1(Inf)、c140013_g1(3.7728)、c154815_g1(Inf) c173122_g1(7.8456)c168962_g1(Inf) c163915_g1(Inf)、c156536_g2(2.9685) c156536_g1(1.2021)（calmodulin）、c147877_g2（mitogen-activated protein kinase kinase 4/5，EC:2.7.12.2）主要调节植物病原互作，病原互作相关基因与植物的抗性息息相关；上调的差异基因c154320_g1、c125968_g1(5.6972)c171764_g2(5.6213)、c172816_g3(5.9018)、c237479_g1(3.5347) c161991_g1(2.2688)、c172816_g4(2.2231)、c113669_g1(5.4648) c161991_g2(1.9362)、c106411_g1调节氨基酸和核苷酸糖代谢；上调的差异基因c151010_g1、c156146_g2(2.1536) c171563_g1(2.3117)、c156146_g1、c167910_g3、c152959_g1、c131089_g1(Inf)c144654_g1(6.9885)主要调节植物昼夜节律。
处理后下调的差异基因显著富集的通路有光合作用、光合作用-天线蛋白、卟啉和叶绿素代谢、光合生物体中的碳固定、脂肪酸生产等。如下调的差异基因c237057_g1、c13699_g1、c237057_g1、c76799_g1、c158296_g1、c152833_g1、c152881_g1、c158839_g1、c134148_g1、c140011_g1、c153193_g1、c131998_g1、c48634_g1、c71692_g1、c42744_g1、c153649_g1调节光合作用；下调的差异基因c149517_g1、c139230_g1、c152435_g1、c149445_g1、c153860_g1、c115377_g1、c164577_g1、c155535_g1、c174985_g1、c157047_g1调节光合作用-天线蛋白；下调的差异 基 因 c167743_g1、c167947_g1、c167743_g1、c71809_g1、c165450_g1c168519_g1、c169028_g1c133188_g1、c123480_g1、c185147_g2、c157598_g1、c157598_g1、c155686_g1、c166557_g2(-1.4097) c166557_g1(-1.5029)调节卟啉调节叶绿素代谢；下调的差异基因c104889_g2(-1.5635)c174574_g3(-3.7395)、c104889_g1、c71483_g1、c159323_g1、c172966_g1、c170442_g1、c162112_g2、c154502_g4、c43883_g1、c198353_g1、c168133_g3、c164307_g1调节光合生物体中的碳固定；下调的差异基因c134603_g2、c164323_g1、c166224_g1(-1.1149) c150017_g1(-1.3389)、c163278_g2、c148702_g1、c164323_g2、c134778_g1、c117604_g1调节脂肪酸伸长（图7）。
2.4 qRT-PCR 验证
为验证RNA-seq测序结果的可信度，随机选取10个差异表达基因，其中含上调的5条差异基因（c166605_g1、c199368_g1、c171769_g1、c164813_g2、c153721_g1）和下调的5条差异基 因（c76799_g1、c153860_g1、c149445_g1、c115377_g1、c155535_g1），以 Actin为内参进行qRT-PCR验证。由图8可知，10个基因在低温胁迫下的表达程度不同，但表达趋势与高通量测序结果基本一致，即表明测序结果真实可靠。
图6 上调差异基因（20℃ Vs -4℃）参与的KEGG代谢通路富集结果
Fig.6 Differential up-regulation genes（20℃ Vs -4℃）assignment to KEGG metabolic pathway
横坐标表示该通路下差异表达基因与注释到该通路的所有基因的比值，纵坐标为通路名称。其中点的大小表示富集的基因数目，越大表示富集的基因数目越多，颜色表示富集的可信度，其值越大表示富集效果越明显
The x-coordinate represents the ratio of differentially expressed genes in this pathway to all genes annotated into this pathway, and the y-coordinate represents the pathway name.The size of the center point indicates the number of enriched genes, the larger the point the more genes are enriched, and the color indicates the credibility of the enrichment, the larger the value the more obvious of the enrichment effect
图7 下调差异基因（20℃ Vs -4℃）参与的KEGG代谢通路富集结果
Fig.7 Differential down-regulation genes（20℃ Vs -4℃）assignment to KEGG metabolic pathway
横坐标表示该通路下差异表达基因与注释到该通路的所有基因的比值，纵坐标为通路名称。其中点的大小表示富集的基因数目，越大表示富集的基因数目越多，颜色表示富集的可信度，其值越大表示富集效果越明显
The x-coordinate represents the ratio of differentially expressed genes in this pathway to all genes annotated into this pathway, and the y-coordinate represents the pathway name.The size of the center point indicates the number of enriched genes, the larger the point the more genes are enriched, and the color indicates the credibility of the enrichment, the larger the value the more obvious of the enrichment effect
图8 差异基因的qRT-PCR验证
Fig.8 Validation of DEGs by qRT-PCR
3 讨论
前人关于百合抗冻分子机制研究较少，不良生长环境尤其是零下低温胁迫是植物生长发育中常遇到的自然灾害。本研究采用先进的测序平台Illumina 4000高通量测序平台对常温（20℃）与零下低温（-4℃）处理的兰州百合进行测序，确保测序数据的准确性和序列的高度丰富性［18］。测序结果在NCBI数据库中共注释到了238 301条Unigene，但注释到每个库中的基因最大不超过27.41%，说明还有很多基因没有被注释到NCBI库中，原因可能是数据库信息不足、测序结果中有新基因和基因片段短等。因此，在后续的研究中，我们可以进一步比对测序结果中差异极显著的基因片段，进一步深入发掘兰州百合抗冻转录新基因。
根据KEGG富集分析结果，可以看出共有25条差异基因参与病原互作途径，10条差异基因参与氨基酸和核苷酸糖代谢途径，8条差异基因参与植物昼夜节律途径，15条差异基因参与卟啉调节叶绿素代谢途径，16条差异基因参与光合作用途径，10条差异基因调节光合作用-天线蛋白途径，13条差异基因调节光合生物体中的碳固定途径，9条差异基因调节脂肪酸伸长途径。因此我们可以深入研究以上通路以进一步揭示兰州百合抗冻的分子机制。另外也可以将兰州百合相关的代谢产物作为研究对象进行深入探究。
4 结论
本研究对兰州百合转录组测序后通过KEGG通路分析，结果检测到24 465个差异表达的基因，筛选后获得8 558条差异表达基因，共占差异表达基因总数的10.27%。经富集分析后发现，与蛋白激酶、蛋白磷酸酶、碳代谢以及ABA等相关的基因可以作为研究兰州百合寒冷响应机制的主要研究对象。qRT-PCR分析结果表明，随机选出的10个差异性表达基因的表达趋势与高通量测序结果相一致。此外，还候选了吲哚-3-甘油磷酸合成酶、蛋白磷酸酶、已糖激酶、钙结合蛋白、叶绿素a/b结合蛋白等基因作为与兰州百合零下低温响应相关的应答候选基因，为揭示兰州百合抗冻分子机制奠定了基础。
［1］ 万珠珠，牛来春，谭秀梅，刘敏，吴亮，董草.三种食用百合种质资源限制生长保存及遗传稳定性研究［J］.北方园艺，2016(13)：89-92.doi:10.11937/bfyy.201613024.WAN Z Z, NIU L C, TAN X M, LIU M, WU L, DONG C.Studies on the limited growth preservation and genetic stability of three edible lily germplasm resources［J］.North Garden, 2016(13)：89-92.doi:10.11937/bfyy.201613024.
［2］ 李玉帆，明军，刘新艳，王良桂，袁素霞，刘春，王莹，徐雷锋，袁迎迎.15个百合种和品种的食用性比较研究［J］.园艺学报，2013，40(S)：2693.LI Y F, MING J, LIU Y Y, WANG L G, YUAN S X, LIU C, WANG Y,XU L F, YUAN Y Y.Comparative study on edibility of 15 lily species and varieties［J］.Journal of Horticulture, 2013, 40(S): 2693.
［3］ HUANG J, LIU X H, WANG J M, LU Y M.Transcriptomic analysis of Asiatic lily in the process of vernalization via RNA-seq［J］.Molecular Biology Reports, 2014(41): 3839-3852.doi：10.1007/s11033-014-3250-2.
［4］ BAKHSHAIE M, KHOSRAVI S, AZADI P.Biotechnological advances in Lilium［J］.Plant Cell Reports, 35 2016, 8: 1799-1826.doi：10.1007/s00299-016-2017-8.
［5］ 王晶懋.野生百合卷丹（Mlium lancifolium）响应低温信号的分子机制研究［D］.北京：北京林业大学，2015.WANG J M.The molecular mechanism of the cold response and signaling pathways in Milium lancifolium［D］.Beijing: Beijing Forestry University, 2015.
［6］ 鲁中爽，刘思言，李广隆，王蕊，刘明明，么梦凡，关淑艳，姚丹.大豆油脂的合成途径及关键酶GPAT基因的研究进展［J］.广东农业科学，2018，45(11)：7-13.LU Z S, LIU S Y, LI G L, WANG R, LIU M M, YAO M F, GUAN S Y,YAO D.The synthetic pathway of soybean oil and the research progress of its key enzyme GPAT gene［J］.Guangdong Agricultural Science，2018，45(11)：7-13.
［7］ WANG M, ZHANG X N, LIU J H.Deep sequencing-based characterization of transcriptome of trifoliate to cold stress［J］.BMC Genomics, 2015, 16: 555.doi：10.1186/s12864-015-1629-7.
［8］ HOSHI Y, KONDO M, MORI S, ADACHI Y, NAKANO M,KOBAYASHI H.Production of transgenic lily plants by Agrobacteriummediated transformation［J］.Plant Cell Reports, 2004, 22: 359-364.doi：10.1007/s00299-003-0700-z.
［9］ 杜方.百合不同器官转录组分析及SSR标记开发应用［J］.杭州：浙江大学，2014.DU F.Transcriptome analysis of differert lily organs and development and applications of SSR markers［J］.Hangzhou: Zhejiang University,2014.
［10］ 储俊，许娜，钱立生，赵岩.百合转基因体系及分子育种相关基因的研究进展［J］.安徽科技学院学报，2013, 27( 6) :26-33.CHU J, XU N, QIAN L S, ZHAO Y.Research progresses of transgenic system and molecular breeding related genes in lily［J］.Journal of Anhui University of Science And Technology, 2013, 27(6) :26-33.
［11］ 李明玺，王敏，甘玉迪，刘洋，黄莹捷，万春鹏.靖安白茶芽和叶的转录组数据组装及基因功能注释［J］.现代食品科技, 2018, 34:35.LI M X, WANG M, GAN Y D, LIU Y, HUANG Y J, WAN C P.Transcriptome analysis and gene function annotation of buds and leaves of Camellia sinensis Cultivar Jing’anbaicha［J］.Modern Food Technology, 2018, 34:35.
［12］ 冷暖，刘晓巍，张娜，许立新.草地早熟禾干旱胁迫转录组差异性分析［J］.草业学报，2017,26(12): 128-137.doi:10.11686/cyxb2017130.LENG N, LIU X W, ZHANG N, XU L X.Differential gene analysis of Poa pratensis in response to drought stress［J］.Journal of Grass Industry, 2017, 26(12): 128-137.doi:10.11686/cyxb2017130.
［13］ ZHANG X Y, ZHANG J Z, ZHANG W W.Transcriptome sequencing and de novo analysis of Rosa multi fl ora under cold stress［J］.Acta Physiol Plant , 2016, 38: 164.
［14］ 曾宪海，焦云飞，廖子荣，潘登浪，林位夫.广东不同地区引种油棕叶片解剖结构对油棕抗寒力的影响［J］.广东农业科学,2018,45(8)：50-58.ZENG X H, JIAO Y F, LIAO Z R, PAN D L, LIN W F.Effects of leaf anatomical structure features on cold resistance of oil palm (Elaeis guineensis Jacq.) introduced in different regions in Guangdong Province ［J］.Guangdong Agricultural Science, 2018,45(8)：50-58.
［15］ SHAFIEE-MASOULEH S S, HATAMZADEH A, SAMIZADEH H.Enlarging bulblet by magnetic and chelating structures of nanochitosan as supplementary fertilizer in Lilium［J］.Horticulture Environment，and Biotechnology, 2014, 55(6): 437-444.doi：10.1007/s13580-014-0175-6.
［16］ 李桂荣，连艳会，程珊珊，辛董董，牛生洋，扈惠灵.低温胁迫下不同无核葡萄品种抗寒性的分析［J］.西南农业学报, 2018, 31 (11):2399-2406.LI G R, LIAN Y H, CHENG S S, XIN D D, NIU S Y, HU H L.Analysis of cold resistance of different seedless grape varieties under low temperature stress［J］.Southwest Agricultural Journal, 2018, 31 (11):2399-2406.
［17］ 王茂辉，钟春燕，罗文龙，聂金泉，郭涛，王慧，陈志强.基于分子标记和高通量测序的基因精细定位［J］.广东农业科学，2018,45(10)：1-8.WANG M H, ZHONG C Y, LUO W L, NIE J Q, GUO T, WANG H,CHEN Z Q.Fine mapping of genes based on molecular markers and high-throughput sequencing［J］.Guangdong Agricultural Science,2018, 45(10)：1-8.
［18］ REVILLA P, RODRIGUEZV M, ORDASA.Association mapping for cold tolerance in two large maize inbred panels［J］.BMC Plant Biology, 2016, 16:127.doi:10.1186/s12870-016-0816-2.