百合(Lilium spp.)是百合科(Liliaceae)百合属(Lilium)所有种类的总称，通常为具有根状茎、块茎或鳞茎的多年生草本。西伯利亚百合(Lilium ‘Siberia’)属于东方百合，由天香百合(L. auratum Lindl.)、药百合(L. speciosum var. gloriosoides Baker)等品种杂交而成，由于其花大而美丽，开放时气味芳香，是非常受消费者欢迎的鲜切花之一^[1]。切花的生产需要优质的商品种球，鳞片繁殖成本低、操作简便、繁殖系数高，是目前种球商品化生产最常用的技术途径^[2]。百合种球从鳞片扦插到成球的生产周期较长，小鳞茎培养成合格的种球一般需要2~3年，如何提高小鳞茎的繁殖系数、促进鳞茎膨大对百合种球的生产至关重要。因此，研究促进百合小鳞茎萌发的分子机制及挖掘调控萌发的潜在基因具有重要意义。

AP2/ERF转录因子是植物中的一大类因子，该家族参与多种生物学过程，包括植物形态建成、病菌防御、对各种胁迫的响应机制、激素信号转导等^[3]。AP2/ERF家族转录因子主要分为脱水反应元件结合蛋白(dehydration-responsive element binding protein，DREB)、乙烯响应因子(ethylene response factor，ERF)、APETALA2 (AP2)和ABI3/VP1相关因子(related to ABI3/VP1，RAV)四大亚家族以及少数未分类因子(Soloists)^[4]。ERF亚家族不仅参与植物发育和生理过程的多种功能，还参与种子萌发和芽休眠。SlPti4是番茄(Solanum lycopersicum L.)果实成熟、种子萌发和非生物胁迫响应的重要调节因子^[5]；油桃(Prunus persica var. nectarina)中的ERF家族成员PpeEBB1可调控桃芽萌发^[6]；在砀山酥梨中，PpyERF060、PpyABF3和PpyMADS71互作可整合乙烯与脱落酸的信号通路，进而调控梨芽休眠进程^[7]；樱桃(Prunus pseudocerasus)中的PpcERF5表达被脱落酸(abscisic acid，ABA)诱导上调，这可能是其参与樱桃花芽休眠调控的方式之一，过表达PpcERF5基因还可提高转基因拟南芥(Arabidopsis thaliana)种子的萌发率^[8]。

褪黑素在调节植物生长、种子发芽、叶片衰老、光合效率、果实成熟等生理生长过程以及调控植物增强抗逆性、抵御细菌和真菌病害方面具有重要作用^[9-11]。褪黑素是与植物激素相关的基因表达的重要调节剂，它与植物激素之间表现出复杂的交叉作用，如参与生长素、吲哚乙酸、赤霉素、细胞分裂素、脱落酸、乙烯、水杨酸、茉莉酮、油菜素类固醇的代谢^[12-14]。拟南芥中的褪黑素可通过降低IAA17及衰老相关基因的表达抑制植物的衰老进程^[15]；褪黑素处理可以通过影响氨基酸代谢来提高冷藏果实的可溶性糖含量、香气水平和果实品质^[16]；番茄中低浓度的外源褪黑素处理可以促进不定根的形成^[17]；相反，高浓度的褪黑素处理会抑制生长素的生物合成以及转运途径中许多基因的表达^[18]；用褪黑素处理多年生黑麦草(Lolium perenne L.)，增加了热应激下褪黑素的内源含量，从而降低ABA含量^[19]；外源褪黑素可促进Cr胁迫下百日草(Zinnia elegans)的种子萌发，提高其对Cr胁迫的抗性^[20]；在黄瓜(Cucumis sativus L.)幼苗中，褪黑素通过调节ABA和赤霉素(gibberellin，GA₄)的生物合成和分解代谢以减轻NaCl的胁迫抑制作用，促进发芽^[21]。

百合萌发的相关作用机制很复杂，高通量转录组测序技术为探究萌发分子机制提供了途径和方法，可用于系统研究百合萌发过程的相关代谢通路。AP2/ERF转录因子在植物生长发育、代谢、生物和非生物胁迫方面扮演着重要角色，而关于促进百合鳞茎萌发的报道较少。为深入研究百合萌发阶段的分子调控机制，通过对褪黑素处理下的百合小鳞茎进行转录组测序分析，并对差异表达的基因进行功能注释，筛选受褪黑素诱导的LoERF96基因；对该基因进行克隆及生物信息学分析，利用实时荧光定量PCR (RT-qPCR)分析其在小鳞茎4个不同生长发育时期的基因表达量，为探究百合AP2/ERF转录因子在萌发过程中的生物学功能提供理论基础和参考。

1.   材料与方法

1.1   植物材料

供试百合小鳞茎由西伯利亚百合大种球(云南沃德球根生物科技有限公司)鳞片包埋繁殖获得。将西伯利亚百合大种球中外层的鳞片剥下，包埋于椰糠基质中，在28 ℃恒温培养箱中遮光培养45 d，然后将获得的小鳞茎置于4 ℃环境中低温冷藏70 d，冷藏结束后将小鳞茎洗净、晾干待用。

1.2   样本处理

选取生长状态良好、大小一致、鲜质量为(0.20±0.10) g的百合小鳞茎，随机平均分为2组，每组72个(含3次重复)，其中一组为处理组，用100 μmol/L褪黑素溶液浸泡2 h；另一组为对照组，用蒸馏水浸泡2 h。将浸泡处理后的小鳞茎放在滤纸上晾干，将小鳞茎栽培到椰糠为基质的营养钵中，每个营养钵栽培3个；将营养钵放入25 ℃的光照恒温培养箱中培养，小鳞茎破土长出芽即认为萌发，分别在第0、7、14和21天取样，液氮速冻后置于−80 ℃超低温冰箱保存，用于后续试验。

1.3   总RNA提取

按照RNAprep Pure多糖多酚植物总RNA提取试剂盒(天根生化科技有限公司，DP441)说明书，分别提取4个不同时期处理组和对照组的西伯利亚百合小鳞茎总RNA，用NanoDrop 2000分光光度计检测RNA的纯度和浓度，用1%琼脂糖凝胶电泳检测RNA的质量；检测合格后，取部分试验第21天的RNA进行转录组测序及目的基因克隆试验；4个时期的全部RNA材料均用于后续的基因表达分析。

1.4   转录组测序

RNA检测合格后进行转录组文库构建，用含有Oligo (dT)的磁珠富集真核生物mRNA将mRNA随机打断，以mRNA为模板合成cDNA第1链和第2链并进行cDNA纯化，然后进行末端修复、加A尾并连接测序接头，PCR富集得到cDNA文库。用N50数值对组装文库进行质量评估检测，Unigene N50越长，说明组装质量越好。 N50的定义为：将所有Unigene从长到短排序，并依次累加长度，当累加片段长度达到总片段长度的50%时，对应的片段长度即为N50。文库质检合格后，用Illumina NovaSeq6000测序平台进行测序(委托北京百迈客生物科技有限公司完成)。

1.5   Unigene功能注释

将Unigene序列与非冗余蛋白质序列(non-redundant protein sequence，NR)数据库、Swiss-Prot蛋白质序列 (Swiss-Prot protein sequence，Swiss-Prot)数据库、原核生物同源蛋白簇 (clusters of orthologous groups，COG)、真核生物同源蛋白簇 (eukaryotic orthologous groups，KOG)、 eggNOG4.5数据库、Pfam数据库(protein family，http://pfam.sanger.ac.uk/)、京都基因与基因组百科全书 (Kyoto encyclopedia of genes and genomes，KEGG)和(gene ontology，GO，http://www.geneontology.org/)数据库比对，获得Unigene的注释信息。

1.6   差异基因的筛选

将差异倍数|log₂FC|≥2和P<0.01作为标准筛选差异表达基因，使用DESeq2软件进行差异分析。

1.7   LoERF96序列克隆

将褪黑素诱导21 d的百合小鳞茎总RNA用All-In-One 5× RT MasterMix (abm，G592)逆转录为cDNA，基于转录组数据获得差异基因LoERF96的序列，用Premier 6.0设计扩增全长序列的特异性引物(表1)，委托擎科生物科技股份有限公司昆明分公司进行引物合成，然后进行LoERF96基因的cDNA扩增。PCR总反应体系为25.0 µL，包括2×T5 Super PCR Mix (Basic) 12.5 µL，正、反向引物各1.0 µL，cDNA模板1.0 µL，ddH₂O 9.5 µL。PCR反应程序为：98 ℃预变性3 min，98 ℃变性10 s，63 ℃退火10 s，72 ℃延伸15 s，共35个循环，循环结束后72 ℃终延伸2 min，4 ℃保存。用1%琼脂糖凝胶电泳检测PCR产物，对扩增得到的目的片段进行切胶，再使用DNA凝胶回收试剂盒(擎科生物科技股份有限公司，TSP602-200)进行胶回收和纯化；将目的片段与pClone007-T载体连接，并转化大肠杆菌DH5α感受态细胞，涂到含Amp抗生素的LB培养基平板后，倒置于37 ℃培养箱过夜培养；次日挑取单菌落并对菌液进行PCR验证，将阳性克隆菌液送至擎科生物科技股份有限公司测序。

表  1  引物序列

Table  1.  Primer sequence

基因名称 gene name	正向引物序列 (5'→3') forward primer sequence	反向引物序列 (5'→3') reverse primer sequence	用途 usage
LoERF96	ATGCACCAGCTCATGGAAGGAAGC	TTATTTGTCCCTGTCGTCGCTCTCCT	基因克隆 gene cloning
LoERF96	GGAGCCCAAACACCAACCCAA	AAGTGCCGAGCCAGACACGTA	荧光定量PCR real-time PCR
ERF105	CAGTTCTAGCACCTGCCTCCCT	GTATCGCTCTGGCGTCCGTATCT
CDC5	CCGTCTCAAAGGACACCCTACCA	TGTGGAGTGGGAGTGTGAAGAGG
MYB36	CTGCCGTTGCCGAAGATGAACT	CACCTCCACCACCACTACCAGAA
WRKY28	TGCGATGGCAACATCTGCATCT	ATGCGGGCATTGCAGGAATTGA
WOX8	TGGCTCCGAACACCTTGTGGA	GCTCCTCGAACTTCCCTTGTGC
LOB30	GCTCTCATCTCCGCCGAAACAC	TCATCGTGCCTCGCTCTGTTCT
LoFP	TCGCCTACATCGCTAACC	TTCCCAATAATCGCAAGACC	荧光定量内参 reference gene of real-time

下载: 导出CSV 

| 显示表格

1.8   RT-qPCR验证转录组及不同时期LoERF96的表达分析

为了验证转录组测序结果的准确性，随机挑选CDC5、LOB30、LoERF105、MYB36、WOX8和WRKY28共6个差异基因进行RT-qPCR验证。用Premier 6.0设计RT-qPCR基因特异性引物(表1)，按照BlasTaq™ 2× qPCR MasterMix (abm，G891)试剂盒说明书，在Applied Biosystems QuantStudio 5实时荧光定量PCR仪上进行实时荧光检测，选择百合F-box家族蛋白LoFP基因作为内参基因^[22]。RT-qPCR反应体系为20.0 µL，包括 BlasTaq^TM 2× qPCR MasterMix 10.0 µL，10 µmol/L正、反向引物各0.5 µL，cDNA 1.0 µL和ddH₂O 8.0 µL。 RT-qPCR扩增程序为：95 ℃ 3 min， 95 ℃ 15 s，60 ℃ 1 min，共40个循环。每个样品设置3次生物学重复。采用2^−ΔΔCt法计算基因相对表达量；采用Excel统计数据；采用Origin 2021绘制图表；采用SPSS 26进行显著性分析。

1.9   LoERF96的生物信息分析

使用Expasy (https://web.expasy.org/protparam/)分析LoERF96蛋白的序列长度、分子质量、等电点、不稳定指数、脂溶性指数、亲水性等理化性质；使用DeepTMHMM (https://dtu.biolib.com/DeepTMHMM)解析蛋白跨膜区域；使用NCBI中CD-Search (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi)工具预测氨基酸序列的保守结构域；使用SOPMA (https://npsa-prabi.ibcp.fr/cgi-bin/npsa_automat.pl?page=/NPSA/npsa_sopma.html)预测蛋白质的二级结构；使用SWISS-MODE预测蛋白质的三级结构 (https://swissmodel.expasy.org/)；使用Plant-mPLoc (http://www.csbio.sjtu.edu.cn/bioinf/plant-multi/)在线分析工具预测蛋白质的亚细胞定位。使用NCBI数据库中的Protein Blast工具寻找LoERF96蛋白氨基酸的同源序列；使用MEGA 11.0进行多序列比对，并采用邻接法构建系统进化树，Bootstraps设置重复1000次。

2.   结果与分析

2.1   序列的组装和质量控制

完成6个样品的转录组测序，获得的有效数据总量为38.95 Gb，各样品的有效数据量均达到6.03 Gb，GC碱基含量为50.51%，Q30碱基的百分比大于87.88%。组装共获得132750个转录本和78100条Unigene，其中长度大于1 kb的Unigene有20431条，Unigene的N50为1 355 bp，表明组装的完整性较高(表2)。

表  2  转录本和Unigene序列的长度分布

Table  2.  Length distribution of transcripts and Unigenes

序列长度/bp sequence length	转录本 transcripts			Unigenes
	数量 number	比例/% percentage		数量 number	比例/% percentage
>200~300	18515	13.95		15442	19.77
>300~500	34749	26.18		23329	29.87
>500~1000	35451	26.70		18898	24.20
>1000~2000	28747	21.65		13313	17.05
>2000	15288	11.52		7118	9.11
总数量 total number	132750			78100
总长度/bp total length	131646756			66907113
N50长度/bp N50 length	1512			1355
平均长度/bp mean length	991.69			856.69

下载: 导出CSV 

| 显示表格

2.2   Unigene功能注释统计

在NR、Swiss-Prot、COG、KOG、eggNOG4.5、Pfam、KEGG和GO共8个数据库获得42458条Unigene的注释结果，每个数据库注释到的基因数量分别为40403、21492、10357、21880、29112、26857、24294和31149。

2.3   差异基因的筛选结果

对褪黑素处理组和对照组的转录组数据进行差异表达分析，得到780个差异表达基因，其中上调表达基因392个，下调表达基因388个(图1)。

图  1  差异表达基因火山图(a)和聚类热图(b)

注：图a)中，FC为差异倍数；图b)中，CK1、CK2和CK3为蒸馏水对照组的3个重复，MT1、MT2和MT3为褪黑素处理组的3个重复。

Figure  1.  Volcano map and cluster heatmap of differentially expressed genes

Note: In Fig. a), FC indicates fold change; in Fig. b), CK1, CK2, and CK3 are three repeats of the distilled water control group; MT1, MT2, and MT3 are three repeats of the melatonin treated group.

下载: 全尺寸图片 幻灯片

2.4   差异表达基因的GO分类富集分析和代谢通路KEGG富集分析

所有差异表达基因的GO富集分析结果(图2)显示：在生物过程中富集到差异表达基因主要包括代谢过程、细胞过程、生物调节、刺激响应、信号、生殖过程等；在细胞组分中富集到差异基因功能主要包括细胞、细胞器、膜、大分子复合物等；在分子功能中富集到差异基因功能主要包括结合活性、催化活性、转运活性、核酸结合转录因子活性、分子功能调节等。

图  2  差异表达基因的GO富集图

注/Note：生物过程(biological process)：代谢过程(metabolic process)、细胞过程(cellular process)、单体过程(single-organism process)、生物调节(biological regulation)、定位(localization)、刺激响应(response to stimulus)、细胞成分组织或生物发生(cellular component organization or biogenesis)、发育过程(developmental process)、信号(signaling)、多细胞生物过程(multicellular organismal process)、生殖(reproduction)、生殖过程(reproductive process)、多有机体过程(multi-organism process)、解毒(detoxification)、免疫系统过程(immune system process)、生长(growth)、节律过程(rhythmic process)、移动(locomotion)、生物黏附(biological adhesion)、行为(behavior)、细胞聚集(cell aggregation)。细胞组分(cellular component)：细胞(cell)、细胞部分(cell part)、细胞器(organelle)、膜(membrane)、膜部件(membrane part)、大分子复合物(macromolecular complex)、细胞器部分(organelle part)、膜腔(membrane-enclosed lumen)、细胞外区域(extracellular region)、细胞连接(cell junction)、合胞体(symplast)、其他生物体(other organism)、其他生物体部分(other organism part)、超分子复合物(supramolecular complex)、胞外区部分(extracellular region part)、类核物质(nucleoid)、突触(synapse)、突触部分(synapse part)、病毒粒子(virion)、病毒粒子部分(virion part)。分子功能(molecular function)：结合活性(binding)、催化活性(catalytic activity)、转运活性(transporter activity)、结构分子活性(structural molecule activity)、核酸结合转录因子活性(nucleic acid binding transcription factor activity)、分子功能调节(molecular function regulator)、电子载体活性(electron carrier activity)、抗氧化活性(antioxidant activity)、信号转运活性(signal transducer activity)、分子转运活性(molecular transducer activity)、转录因子活性和蛋白质结合(transcription factor activity and protein binding)、营养活性(nutrient reservoir activity)、蛋白质标签(protein tag)、金属伴侣活性(metallochaperone activity)、翻译调节因子活性(translation regulator activity)。

Figure  2.  GO enrichment map of differentially expressed genes

下载: 全尺寸图片 幻灯片

对差异表达基因进行KEGG代谢通路的富集分析，筛选出显著富集的代谢途径50条，主要包括植物激素信号转导、植物MAPK信号途径、核糖体、苯丙烷生物合成、淀粉和蔗糖代谢、光合作用、碳代谢、内质网中的蛋白质加工、植物病原菌互作等(图3)。

图  3  差异表达基因的KEGG代谢通路富集

注/Note：细胞过程(cellular processes)：自噬其他(autophagy-other)，内分泌(endocytosis)。环境信息加工(environmental information processing)：植物激素信号转导(plant hormone signal transduction)，植物MAPK信号途径(plant MAPK signaling pathway)。遗传信息加工(genetic information processing)：核糖体(ribosome)，内质网中的蛋白质加工(protein processing in endoplasmic reticulum)，RNA转运(RNA transport)，氨酰基tRNA生物合成(aminoacyl-tRNA biosynthesis)，泛素介导的蛋白水解(ubiquitin mediated proteolysis)，真核生物核糖体的生物发生(ribosome biogenesis in eukaryotes)。RNA降解(RNA degradation)，RNA聚合酶(RNA polymerase)，mRNA监测途径(mRNA surveillance pathway)，剪接体(spliceosome)。代谢(metabolism)：苯丙烷生物合成(phenylpropanoid biosynthesis)，淀粉和蔗糖代谢(starch and sucrose metabolism)，光合作用(photosynthesis)，色氨酸代谢(tryptophan metabolism)，半胱氨酸和蛋氨酸代谢(cysteine and methionine metabolism)，糖酵解/糖异生(glycolysis/gluconeogenesis)，二苯乙烯、二芳基庚烷和姜酚的生物合成(stilbenoid, diarylheptanoid and gingerol biosynthesis)，类黄酮生物合成(flavonoid biosynthesis)，碳代谢(carbon metabolism)，氧化磷酸化(oxidative phosphorylation)，酪氨酸代谢(tyrosine metabolism)，氨基糖和核苷酸糖代谢(amino sugar and nucleotide sugar metabolism)，缬氨酸、亮氨酸和异亮氨酸降解(valine, leucine and isoleucine degradation)，光合生物的固碳作用(carbon fixation in photosynthetic organisms)，氨基酸的生物合成(biosynthesis of amino acids)，角质、木栓和蜡的生物合成(cutin, suberine and wax biosynthesis)，甘油酯代谢(glycerolipid metabolism)，异喹啉生物碱生物合成(isoquinoline alkaloid biosynthesis)，氰基氨基酸代谢(cyanoamino acid metabolism)，维生素B6代谢(vitamin B6 metabolism)，脂肪酸降解(fatty acid degradation)，肌醇磷酸代谢(inositol phosphate metabolism)，玉米素生物合成(zeatin biosynthesis)，谷胱甘肽代谢(glutathione metabolism)，泛醌和其他萜类醌生物合成(ubiquinone and other terpenoid-quinone biosynthesis)，甜菜碱的生物合成(betalain biosynthesis)，果糖和甘露糖代谢(fructose and mannose metabolism)，乙醛酸和二羧酸代谢(glyoxylate and dicarboxylate metabolism)，丙氨酸、天冬氨酸和谷氨酸代谢(alanine, aspartate and glutamate metabolism)，甘油磷脂代谢(glycerophospholipid metabolism)，抗坏血酸和阿糖二酸代谢(ascorbate and aldarate metabolism)，甘氨酸、丝氨酸和苏氨酸代谢(glycine, serine and threonine metabolism)，糖基磷脂酰肌醇(GPl)-锚生物合成[glycosylphosphatidylinositol (GPl)-anchor biosynthesis]，类胡萝卜素生物合成(carotenoid biosynthesis)。生物系统(organismal systems)：植物病原菌互作(plant-pathogen interaction)，植物昼夜节律(circadian rhythm of plant)。

Figure  3.  Enrichment of KEGG metabolic pathway of differentially expressed genes

下载: 全尺寸图片 幻灯片

2.5   小鳞茎萌发相关基因LoERF96的克隆及测序结果

根据转录组数据分析结果，筛选出1个褪黑素诱导后显著上调表达的ERF转录因子，以西伯利亚百合小鳞茎cDNA为模板，通过PCR扩增得到ERF基因的cDNA全长，1%的琼脂糖凝胶电泳显示在250~500 bp之间有1条清晰的条带(图4a)，命名为LoERF96，切胶回收后测序并与转录组数据的基因序列比对一致，全长为435 bp，编码144个氨基酸(图4b)。

图  4  LoERF96的克隆(a)、氨基酸序列(b)、蛋白跨膜结构(c)、基因结构域(d)以及蛋白的二级结构(e)和三级结构(f)

Figure  4.  Cloning of LoERF96 (a), amino acid sequence (b), protein transmembrane structure (c), gene domain (d), and the secondary structure (e) and tertiary structure (f) of protein

下载: 全尺寸图片 幻灯片

2.6   LoERF96编码蛋白的生物信息学分析

LoERF96蛋白相对分子质量约为16.32 ku，理论等电点为4.94，分子式为C₆₉₂H₁₀₇₈N₂₁₈O₂₃₂S₅，共含有2225个原子。该蛋白由20种不同的氨基酸组成，以精氨酸的含量最高，占总量的12.5%；半胱氨酸、组氨酸和赖氨酸的含量最低，分别占总量的1.4%。LoERF96蛋白有27个带负电荷的残基(Asp+Glu)和20个带正电荷的残基(Arg+Lys)，不稳定系数为70.95，表明该蛋白属于不稳定蛋白；脂肪系数为53.61；亲水性系数为−1.026，表明该蛋白为亲水性蛋白。预测LoERF96蛋白不存在跨膜区，不属于跨膜蛋白(图4c)。LoERF96基因结构域分析结果(图4d)显示：其在氨基酸40和103之间有1个AP2保守结构域，是典型的AP2/ERF家族成员，属于ERF亚家族，其亚细胞定位于细胞核。

由LoERF96蛋白的二级结构(图4e)可知：该蛋白编码的144个氨基酸中有4种二级结构，其中，α-螺旋占31.94%，β-折叠占4.17%，无规则卷曲占52.08%，延伸链占11.81%，蛋白质以无规则卷曲为主，说明其结构不稳定。以含有1个AP2结构域的A0A0D6RAM1_ARACU (A0A0D6RAM1.1.A)基因为模型预测LoERF96蛋白质的三级结构，其覆盖率为56.00%，序列相似性为70.37%，以无规则卷曲和α-螺旋为主(图4f)，符合二级结构预测结果。

2.7   LoERF96蛋白的同源性与系统进化树分析

进化树分析结果(图5)表明：LoERF96氨基酸序列与番茄(Solanum lycopersicum A0A3Q7I5Y9.1)、拟南芥(Arabidopsis thaliana Q8VY90.1)、烟草(Nicotiana tabacum Q40478.1)、蒺藜苜蓿(Medicago truncatula A9Y0K9.1)、喜马拉雅凤仙花(Impatiens glandulifera XP_047312446.1)等物种的同源性较高。

图  5  基于LoERF96氨基酸序列构建的系统进化树

Figure  5.  Constructing a phylogenetic tree based on LoERF96 amino acid sequence

下载: 全尺寸图片 幻灯片

2.8   LoERF96在小鳞茎不同生长发育时期的表达

由图6可知：褪黑素处理后，LoERF96的相对表达量总体呈先上升后下降的趋势。第7天时，该基因的表达量比对照组高2.08倍；第14天时，其表达量较对照组极显著高10.06倍，第21天时，其表达量比对照组高2.94倍。可见，LoERF96基因受褪黑素诱导表达，且可以快速响应并表达。

图  6  不同生长时期小鳞茎的LoERF96相对表达量

注：“**”表示差异极显著 (P<0.01)；下同。

Figure  6.  Relative expression levels of LoERF96 in small bulbs at different growth stages

Note: “**” indicates extremely significant differences (P<0.01); the same as below.

下载: 全尺寸图片 幻灯片

2.9   RT-qPCR验证转录组数据

6个差异表达基因的RT-qPCR试验结果如图7所示。RT-qPCR结果和第21天转录组数据(表3)的基因表达趋势完全一致，表明转录组测序获得的数据比较准确。

图  7  RT-qPCR验证转录组数据

注：“*”表示差异显著 (P<0.05)。

Figure  7.  RT-qPCR validation of transcriptome data

Note: “*” indicates significant differences (P<0.05).

下载: 全尺寸图片 幻灯片

表  3  转录组数据的基因表达量

Table  3.  Gene expression levels of transcriptome data

基因名称 gene name	基因号 gene ID	FPKM						P值 P-value	基因表达量差异倍数 log2FC
		CK1	CK2	CK3	MT1	MT2	MT3
LoERF105	TRINITY_DN1258_c0_g1	24.52	44.71	20.79	93.13	77.71	43.91	0.0026	1.4920
CDC5	TRINITY_DN10098_c0_g1	113.32	92.5	110.89	43.36	49.1	40.51	0.0006	−1.0428
MYB36	TRINITY_DN11434_c0_g2	6.45	3.56	5.58	1.69	1.81	1.17	0.0029	−1.4653
WRKY28	TRINITY_DN11579_c0_g1	33.2	8.76	35.61	4.39	10.35	6.76	0.0037	−1.5668
WOX8	TRINITY_DN6545_c0_g1	33.92	16.71	29.25	7.78	6.44	14.02	0.0003	−1.3179
LOB30	TRINITY_DN10132_c0_g1	7.19	5.81	13.36	2.37	3.67	3.67	0.0064	−1.2463
注：FPKM. 每千个碱基的转录每百万映射读取的碎片数，用于表示基因表达量；FC. 差异倍数；CK1、CK2和CK3为蒸馏水对照组的3个重复，MT1、MT2和MT3为褪黑素处理组的3个重复。 Note: FPKM. fragments per kilobase of exon model per million mapped fragments, indicating gene expression levels; FC. fold change; CK1, CK2, and CK3 are three repeats of the distilled water control group; MT1, MT2, and MT3 are three repeats of the melatonin treated group.

下载: 导出CSV 

| 显示表格

3.   讨论

本研究对外源褪黑素处理下西伯利亚百合小鳞茎生长发育时期的转录组数据进行分析，KEGG富集分析结果显示：差异表达基因显著富集在包括植物激素信号转导、植物MAPK信号途径、核糖体、苯丙烷生物合成、淀粉和蔗糖代谢、光合作用、碳代谢、内质网中的蛋白质加工、植物病原菌互作等途径。植物激素代谢和信号转导途径中的不同激素可以单独或相互调控发挥作用，调控植物生长发育、开花结果、衰老等过程^[23]。褪黑素与大多数植物激素的变化有关，它在植物生命周期的每个发育阶段起关键作用，大量研究表明：植物内源激素水平与休眠解除和萌发密切相关，但其在百合萌发过程中的调控机制尚不明确。

通过筛选褪黑素处理下的西伯利亚百合转录组数据，从小鳞茎中成功克隆出可能与萌发有关的乙烯响应因子LoERF96基因， 其cDNA全长为435 bp，编码144个氨基酸。对其蛋白理化性质进行分析，推测该蛋白属于不稳定的亲水性蛋白，且无跨膜结构；其含有1个AP2保守结构域， 是典型的AP2/ERF家族成员， 这与杧果(Mangifera indica L.)的MiERF2基因也只存在单一AP2结构域^[24]的研究结果一致。蛋白二级结构预测结果显示：LoERF96蛋白以无规则卷曲和α-螺旋为主，β-折叠和延伸链占比最少，无规则卷曲为主说明其蛋白质结构不稳定，这与苦荞[Fagopyrum tataricum (L.) Gaertn]^[25]FtERF基因以及木薯(Manihot esculenta Crantz)^[26]MeERF5基因的编码蛋白结构类似，也符合LoERF96蛋白的三级结构预测结果。亚细胞定位预测LoERF96基因位于细胞核，氨基酸多序列同源比对发现LoERF96蛋白与番茄、拟南芥、烟草、蒺藜苜蓿、喜马拉雅凤仙花等物种的同源性较高。相关研究表明：小乙烯响应因子ERF96参与拟南芥ABA反应的调控，是ABA反应的正调节因子^[27]。因此，推测LoERF96在百合中也具有类似的功能。

褪黑素主要参与植物发芽、生长、根茎生成、促进果实成熟、延缓衰老等过程，且褪黑素与大多数植物激素的变化有关^[28]。外源高浓度的褪黑素可以抑制或延迟种子发芽，表明其在维持种子休眠和防止过早发芽方面发挥作用^[29]。褪黑素可以调控植物信号转导通路中ABA和GA相关基因的表达，诱导种子萌发并解除休眠^[30]。AP2/ERF转录因子参与植物生长发育、胁迫响应、防御反应、激素信号转导、代谢物调控等，还可以反馈调节乙烯、细胞分裂素、赤霉素、脱落酸等植物激素的生物合成^[31]。杏仁(Prunus dulcis)中的PdERF24可能在休眠杏仁对超低温冷冻胁迫的响应中发挥重要作用；在蓝莓(Vaccinium spp.)花芽休眠解除过程中，多数VcERF基因均能表达^[32]。本研究中，在褪黑素处理下，4个不同生长发育时期的西伯利亚百合小鳞茎的LoERF96基因表达量均高于对照组，且都高于第1天的基因表达量，说明LoERF96在褪黑素处理后可以快速响应并表达。该基因的相对表达量总体呈先上升后下降的趋势，第14天时达到峰值，且极显著高于对照组；虽然之后其基因表达量下降，但仍高于对照组，这与前人的研究结果^[33]一致。本研究表明：乙烯响应因子LoERF96受褪黑素诱导表达，但该基因在萌发、生长发育等过程中的生物学功能还需进一步研究。

4.   结论

本研究对褪黑素处理下的西伯利亚百合小鳞茎进行转录组测序，通过转录组数据筛选得到780个差异表达基因，其中，上调表达基因392个，下调表达基因388个。这些差异表达基因的KEGG代谢途径主要包括植物激素信号转导、植物MAPK信号途径、苯丙烷生物合成、淀粉和蔗糖代谢、光合作用、碳代谢、内质网中的蛋白质加工等。LoERF96基因是典型的AP2/ERF家族成员，其编码蛋白为亲水性且无跨膜结构。本研究证明了LoERF96基因受褪黑素诱导表达，为探究百合AP2/ERF转录因子在萌发和生长发育过程中的生物学功能提供了理论基础和参考。